Skip to content
Snippets Groups Projects
Commit d0af6322 authored by Mirko Birbaumer's avatar Mirko Birbaumer
Browse files

added cohen's d Wert

parent eb52fd23
No related branches found
No related tags found
No related merge requests found
Pipeline #664227 passed with stage
in 6 hours, 52 minutes, and 8 seconds
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import scipy.stats as st import scipy.stats as st
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import warnings import warnings
warnings.filterwarnings("ignore") warnings.filterwarnings("ignore")
import arviz as az import arviz as az
import pymc as pm import pymc as pm
import scipy.stats as st import scipy.stats as st
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Wir lesen den Datensatz `students.csv` ein und betrachten die ersten Zeilen des Dataframes. Wir lesen den Datensatz `students.csv` ein und betrachten die ersten Zeilen des Dataframes.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
data = pd.read_csv("./Daten/students.csv") data = pd.read_csv("./Daten/students.csv")
data.head() data.head()
``` ```
%% Output %% Output
stud.id name gender age height weight religion \ stud.id name gender age height weight religion \
1 833917 Gonzales, Christina Female 19 160 64.8 Muslim 1 833917 Gonzales, Christina Female 19 160 64.8 Muslim
2 898539 Lozano, T'Hani Female 19 172 73.0 Other 2 898539 Lozano, T'Hani Female 19 172 73.0 Other
3 379678 Williams, Hanh Female 22 168 70.6 Protestant 3 379678 Williams, Hanh Female 22 168 70.6 Protestant
4 807564 Nem, Denzel Male 19 183 79.7 Other 4 807564 Nem, Denzel Male 19 183 79.7 Other
5 383291 Powell, Heather Female 21 175 71.4 Catholic 5 383291 Powell, Heather Female 21 175 71.4 Catholic
nc.score semester major minor \ nc.score semester major minor \
1 1.91 1st Political Science Social Sciences 1 1.91 1st Political Science Social Sciences
2 1.56 2nd Social Sciences Mathematics and Statistics 2 1.56 2nd Social Sciences Mathematics and Statistics
3 1.24 3rd Social Sciences Mathematics and Statistics 3 1.24 3rd Social Sciences Mathematics and Statistics
4 1.37 2nd Environmental Sciences Mathematics and Statistics 4 1.37 2nd Environmental Sciences Mathematics and Statistics
5 1.46 1st Environmental Sciences Mathematics and Statistics 5 1.46 1st Environmental Sciences Mathematics and Statistics
score1 score2 online.tutorial graduated salary score1 score2 online.tutorial graduated salary
1 NaN NaN 0 0 NaN 1 NaN NaN 0 0 NaN
2 NaN NaN 0 0 NaN 2 NaN NaN 0 0 NaN
3 45.0 46.0 0 0 NaN 3 45.0 46.0 0 0 NaN
4 NaN NaN 0 0 NaN 4 NaN NaN 0 0 NaN
5 NaN NaN 0 0 NaN 5 NaN NaN 0 0 NaN
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
data.columns data.columns
``` ```
%% Output %% Output
Index(['stud.id', 'name', 'gender', 'age', 'height', 'weight', 'religion', Index(['stud.id', 'name', 'gender', 'age', 'height', 'weight', 'religion',
'nc.score', 'semester', 'major', 'minor', 'score1', 'score2', 'nc.score', 'semester', 'major', 'minor', 'score1', 'score2',
'online.tutorial', 'graduated', 'salary'], 'online.tutorial', 'graduated', 'salary'],
dtype='object') dtype='object')
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
female_height = data[data["gender"]== "Female"]["height"][0:50] female_height = data[data["gender"]== "Female"]["height"][0:50]
female_height.head() female_height.head()
``` ```
%% Output %% Output
1 160 1 160
2 172 2 172
3 168 3 168
5 175 5 175
7 156 7 156
Name: height, dtype: int64 Name: height, dtype: int64
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Wir betrachten also die Körpergrösse von Frauen, von welcher wir annehmen, dass diese normalverteilt ist. Zur Ueberprüfung dieser Annahme wollen wir den Q-Q-Plot der Körpergrösse der Studentinnen im Datensatz erstellen. Wir betrachten also die Körpergrösse von Frauen, von welcher wir annehmen, dass diese normalverteilt ist. Zur Ueberprüfung dieser Annahme wollen wir den Q-Q-Plot der Körpergrösse der Studentinnen im Datensatz erstellen.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
st.probplot(female_height,plot=plt) st.probplot(female_height,plot=plt)
plt.title("Grösse Frauen") plt.title("Grösse Frauen")
#plt.savefig('qq_students.png', dpi=300) #plt.savefig('qq_students.png', dpi=300)
``` ```
%% Output %% Output
Text(0.5, 1.0, 'Grösse Frauen') Text(0.5, 1.0, 'Grösse Frauen')
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Aufgrund des QQ-Plots schliessen wir, dass die Körpergrösse von Frauen normalverteilt ist. Wir bezeichnen mit $ X $ die Körpergrösse einer zufällig ausgewählten Studentin, wobei Aufgrund des QQ-Plots schliessen wir, dass die Körpergrösse von Frauen normalverteilt ist. Wir bezeichnen mit $ X $ die Körpergrösse einer zufällig ausgewählten Studentin, wobei
$$ $$
X X
\sim \mathcal{N}(\mu, \sigma^{2}) \sim \mathcal{N}(\mu, \sigma^{2})
$$ $$
Wir wollen nun eine Aussage über $ \mu $ treffen. Zunächst gehen wir davon aus, dass $ \sigma $ bekannt ist, sagen wir $ \sigma=10 $. Wir wollen nun eine Aussage über $ \mu $ treffen. Zunächst gehen wir davon aus, dass $ \sigma $ bekannt ist, sagen wir $ \sigma=10 $.
Es stellt sich nun die Frage, wie wir die Prior-Verteilung für $ \mu $ wählen. Eine Es stellt sich nun die Frage, wie wir die Prior-Verteilung für $ \mu $ wählen. Eine
Beta-Verteilung ist für diese Situation nicht passend, da deren Parameterwerte $ \theta $ nur Werte Beta-Verteilung ist für diese Situation nicht passend, da deren Parameterwerte $ \theta $ nur Werte
von $ 0 $ bis $ 1 $ annimmt. Viele Verteilungen sind in diesem Fall denkbar. Wir gehen der Einfachheit halber vorläufig von einer gleichförmigen Verteilung aus. Wir gehen hier von grossem Unwissen aus, nämlich dass $ \mu $ zwischen 100 und 250cm liegen kann , d.h. alle Werte in diesem Bereich können mit der gleichen 'Wahrscheinlichkeit' vorkommen. Natürlich könnten wir hier mehr Vorwissen einfliessen lassen, was wir später auch machen werden. von $ 0 $ bis $ 1 $ annimmt. Viele Verteilungen sind in diesem Fall denkbar. Wir gehen der Einfachheit halber vorläufig von einer gleichförmigen Verteilung aus. Wir gehen hier von grossem Unwissen aus, nämlich dass $ \mu $ zwischen 100 und 250cm liegen kann , d.h. alle Werte in diesem Bereich können mit der gleichen 'Wahrscheinlichkeit' vorkommen. Natürlich könnten wir hier mehr Vorwissen einfliessen lassen, was wir später auch machen werden.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Als Likelihood-Funktion wählen wir eine Normalverteilungsfunktion mit dem Datensatz `female_height`, der aus 50 Messungen der Körpergrösse von Frauen besteht. Falls $ \sigma=10 $ und $ x=\{x_{1},x_{2},\dots,x_{50} \} =\{160,\ldots, 155\}$, so lautet die Likelihood-Funktion Als Likelihood-Funktion wählen wir eine Normalverteilungsfunktion mit dem Datensatz `female_height`, der aus 50 Messungen der Körpergrösse von Frauen besteht. Falls $ \sigma=10 $ und $ x=\{x_{1},x_{2},\dots,x_{50} \} =\{160,\ldots, 155\}$, so lautet die Likelihood-Funktion
\begin{align*} \begin{align*}
p(x | \mu,\sigma) p(x | \mu,\sigma)
&=\dfrac{1}{\sigma\sqrt{2\pi}}\cdot\exp{\left(-\frac{(x_{1}-\mu)^2}{2\sigma^2}\right)}\cdot\ldots\cdot \dfrac{1}{\sigma\sqrt{2\pi}}\cdot\exp{\left(-\frac{(x_{50}-\mu)^2}{2\sigma^2}\right)}\\ &=\dfrac{1}{\sigma\sqrt{2\pi}}\cdot\exp{\left(-\frac{(x_{1}-\mu)^2}{2\sigma^2}\right)}\cdot\ldots\cdot \dfrac{1}{\sigma\sqrt{2\pi}}\cdot\exp{\left(-\frac{(x_{50}-\mu)^2}{2\sigma^2}\right)}\\
&=\dfrac{1}{10\sqrt{2\pi}}\cdot\exp{\left(-\frac{(160-\mu)^2}{2\cdot 100}\right)}\cdot\ldots\cdot \dfrac{1}{10\sqrt{2\pi}}\cdot\exp{\left(-\frac{(155-\mu)^2}{2\cdot 100}\right)}\\ &=\dfrac{1}{10\sqrt{2\pi}}\cdot\exp{\left(-\frac{(160-\mu)^2}{2\cdot 100}\right)}\cdot\ldots\cdot \dfrac{1}{10\sqrt{2\pi}}\cdot\exp{\left(-\frac{(155-\mu)^2}{2\cdot 100}\right)}\\
\end{align*} \end{align*}
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Diese Likelihood-Funktion betrachten wir nun als abhängig von $ \mu $ und können nun zusammen mit der Prior-Verteilung die Posterior-Verteilung bestimmen. Dies machen wir nun natürlich mit `pymc3`. Diese Likelihood-Funktion betrachten wir nun als abhängig von $ \mu $ und können nun zusammen mit der Prior-Verteilung die Posterior-Verteilung bestimmen. Dies machen wir nun natürlich mit `pymc3`.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Der Input ist sehr ähnlich wie beim Beta-Prior, aber wir müssen hier 2 Parameter spezifizieren, nämlich $ \mu $ und $ \sigma $, da die Likelihood-Funktion von zwei Parametern abhängig ist. Wir wollen zunächst nur einen Parameter, nämlich $ \mu $, mit MCMC bestimmen. Der Input ist sehr ähnlich wie beim Beta-Prior, aber wir müssen hier 2 Parameter spezifizieren, nämlich $ \mu $ und $ \sigma $, da die Likelihood-Funktion von zwei Parametern abhängig ist. Wir wollen zunächst nur einen Parameter, nämlich $ \mu $, mit MCMC bestimmen.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
with pm.Model() as model_g: with pm.Model() as model_g:
μ = pm.Uniform('μ', lower=100, upper=250) μ = pm.Uniform('μ', lower=100, upper=250)
σ = 10 σ = 10
y = pm.Normal('y', mu=μ, sigma=σ, observed=female_height) y = pm.Normal('y', mu=μ, sigma=σ, observed=female_height)
trace_g = pm.sample() trace_g = pm.sample()
az.plot_posterior(trace_g) az.plot_posterior(trace_g)
``` ```
%% Output %% Output
Auto-assigning NUTS sampler... Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag... Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs) Multiprocess sampling (4 chains in 4 jobs)
NUTS: [μ] NUTS: [μ]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds. Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.
<Axes: title={'center': 'μ'}> <Axes: title={'center': 'μ'}>
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Wir stellen also fest, dass der Mittelwert der Posterior-Verteilung bei 163cm liegt und dass 94\% der wahrscheinlichsten $ \mu $'s im Bereich von 160 bis 165cm liegen. Wir stellen also fest, dass der Mittelwert der Posterior-Verteilung bei 163cm liegt und dass 94\% der wahrscheinlichsten $ \mu $'s im Bereich von 160 bis 165cm liegen.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Nehmen wir an, eine Zeitung schreibt, dass die durchschnittliche Körpergrösse von Frauen in der Schweiz bei 175cm liegt. Passt diese Angabe zu unseren Daten und zu unserer Prior-Verteilung? Nehmen wir an, eine Zeitung schreibt, dass die durchschnittliche Körpergrösse von Frauen in der Schweiz bei 175cm liegt. Passt diese Angabe zu unseren Daten und zu unserer Prior-Verteilung?
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Männer sind durchschnittlich eher grösser als Frauen. Aber ist dieser Unterschied auch statistisch relevant? Dazu wählen wir 65 Studenten aus und führen dieselbe Untersuchung durch. Wir stellen fest, dass die Körpergrösse dieser 65 Studenten annähernd einer Normalverteilung folgt. Männer sind durchschnittlich eher grösser als Frauen. Aber ist dieser Unterschied auch statistisch relevant? Dazu wählen wir 65 Studenten aus und führen dieselbe Untersuchung durch. Wir stellen fest, dass die Körpergrösse dieser 65 Studenten annähernd einer Normalverteilung folgt.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
male_height = data[data["gender"]== "Male"]["height"][0:65] male_height = data[data["gender"]== "Male"]["height"][0:65]
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Die Berechnung der Posterior-Verteilung der Körpergrösse von Studenten einerseits und die Berechnung der Posterior-Verteilung Körpergrösse von Studentinnen andererseits kann mittels `pymc3` Die Berechnung der Posterior-Verteilung der Körpergrösse von Studenten einerseits und die Berechnung der Posterior-Verteilung Körpergrösse von Studentinnen andererseits kann mittels `pymc3`
beides in einem Schritt durchgeführt werden. beides in einem Schritt durchgeführt werden.
Dazu müssen wir aber zwei $ \mu$'s spezifizieren, $ \mu_{1} $ für die Frauen und $ \mu_{2} $ für die Männer. In Abbildung unten sehen wir die beiden Posterior-Plots. Dazu müssen wir aber zwei $ \mu$'s spezifizieren, $ \mu_{1} $ für die Frauen und $ \mu_{2} $ für die Männer. In Abbildung unten sehen wir die beiden Posterior-Plots.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
with pm.Model() as model_mul: with pm.Model() as model_mul:
μ_1 = pm.Uniform('μ_1', lower=100, upper=250) μ_1 = pm.Uniform('μ_1', lower=100, upper=250)
μ_2 = pm.Uniform('μ_2', lower=100, upper=250) μ_2 = pm.Uniform('μ_2', lower=100, upper=250)
σ = 10 σ = 10
y_f = pm.Normal('y_f', mu=μ_1, sigma=σ, observed=female_height) y_f = pm.Normal('y_f', mu=μ_1, sigma=σ, observed=female_height)
y_m = pm.Normal('y_m', mu=μ_2, sigma=σ, observed=male_height) y_m = pm.Normal('y_m', mu=μ_2, sigma=σ, observed=male_height)
trace_mul = pm.sample(1000) trace_mul = pm.sample(1000)
az.plot_posterior(trace_mul) az.plot_posterior(trace_mul)
``` ```
%% Output %% Output
Auto-assigning NUTS sampler... Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag... Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs) Multiprocess sampling (4 chains in 4 jobs)
NUTS: [μ_1, μ_2] NUTS: [μ_1, μ_2]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds. Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 2 seconds.
array([<Axes: title={'center': 'μ_1'}>, <Axes: title={'center': 'μ_2'}>], array([<Axes: title={'center': 'μ_1'}>, <Axes: title={'center': 'μ_2'}>],
dtype=object) dtype=object)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Der 94\%-HDI von $ \mu_{1} $ ist $ [160,165] $ und für $ \mu_{2} $ ist es $ [175,180] $. Das heisst, die wahrscheinlichsten Werte für $ \mu_{1} $ und $ \mu_{2} $ überschneiden sich nicht. Wir ziehen daraus den Schluss, dass Männer statistisch relevant grösser sind als Frauen. Der 94\%-HDI von $ \mu_{1} $ ist $ [160,165] $ und für $ \mu_{2} $ ist es $ [175,180] $. Das heisst, die wahrscheinlichsten Werte für $ \mu_{1} $ und $ \mu_{2} $ überschneiden sich nicht. Wir ziehen daraus den Schluss, dass Männer statistisch relevant grösser sind als Frauen.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Wir können die Effektgrösse mit Hilfe von Cohens'd Wert berechnen: ### Cohen's $d$ Wert
%% Cell type:markdown id: tags:
Eine gängige Art, die Effektgrösse des Gruppenunterschieds zu quantifizieren, ist mit Hilfe von Cohen's $d$ Wert, welcher wie folgt definiert ist:
\begin{equation*}
\delta=\frac{\mu_2 - \mu_1}{\sqrt{\frac{\sigma_1^2+\sigma_2^2}{2}}}
\end{equation*}
Dieser Ausdruck besagt, dass die Effektgrösse die Differenz zwischen den Gruppenmittelwerten dividiert durch die gepoolte Standardabweichung beider Gruppen ist. Indem wir die gepoolte Standardabweichung nehmen, standardisieren wir die Differenzen der Gruppenmittelwerte.
Nehmen wir an, wir haben eine Differenz von 1 zwischen den Gruppenmittelwerten und eine gepoolte Standardabweichung von 0.1, dann ist die Effektgrösse grösser als bei der gleichen Differenz der Gruppenmittelwerte und einer gepoolten Standardabweichung von 10.
Cohen's $d$ Wert kann demzufolge als $z$-Score interpretiert werden. Ein $z$-Score ist die Anzahl Standardabweichungen, die der Gruppenmittelwert der ersten Gruppe vom Gruppenmittelwert der zweiten Gruppe abweicht.
Eine weitere Möglichkeit, eine Kennzahl für die Effektgrösse anzugeben, ist die _\"Uberlegenheits-Wahrscheinlichkeit_(englisch _probability of superiority_). Diese ist definiert als die Wahhrscheinlichkeit, dass ein zufällig gewählter Datenpunkt der ersten Gruppe einen grösseren Wert als ein zufällig gewählter Wert der zweiten Gruppe hat. Wenn wir annehmen können, dass die Daten normalverteilt sind, können wir die \"Uberlegenheits-Wahrscheinlichkeit mit Hilfe von Cohen's $d$ Wert berechnen:
\begin{equation*}
\text{ps}=\Phi\left(\frac{\delta}{\sqrt{2}}\right)
\end{equation*}
wobei $\text{ps}$ für _probability of superiority_ steht, $\Phi$ die kumulative Normalverteilung und $\delta$ Cohen's $d$ Wert bezeichnen.
Im Folgenden Beispiel berechnen wir Cohen's $d$ Wert und die Ueberlegenswahrscheinlichkeit mit Hilfe von `PyMC`.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Extract posterior samples # Extract posterior samples
cg_posterior = az.extract(trace_mul) cg_posterior = az.extract(trace_mul)
# Compute mean difference # Compute mean difference
means_diff = cg_posterior["μ_1"] - cg_posterior["μ_2"] means_diff = cg_posterior["μ_1"] - cg_posterior["μ_2"]
# Use observed standard deviations (since they were fixed) # Use observed standard deviations (since they were fixed)
σ_m = male_height.std() σ_m = male_height.std()
σ_f = female_height.std() σ_f = female_height.std()
# Compute pooled standard deviation # Compute pooled standard deviation
pooled_std = np.sqrt((σ_m**2 + σ_f**2) / 2) pooled_std = np.sqrt((σ_m**2 + σ_f**2) / 2)
# Compute Cohen's d # Compute Cohen's d
d_cohen = (means_diff / pooled_std).mean().item() d_cohen = (means_diff / pooled_std).mean().item()
# Compute probability of superiority (PS) # Compute probability of superiority (PS)
ps = st.norm.cdf(d_cohen / np.sqrt(2)) ps = st.norm.cdf(d_cohen / np.sqrt(2))
# Plot posterior of mean difference # Plot posterior of mean difference
fig, ax = plt.subplots(figsize=(7, 4)) fig, ax = plt.subplots(figsize=(7, 4))
az.plot_posterior(means_diff.values, ref_val=0, ax=ax) az.plot_posterior(means_diff.values, ref_val=0, ax=ax)
# Annotate plot with Cohen's d and probability of superiority # Annotate plot with Cohen's d and probability of superiority
ax.set_title("Posterior of Mean Difference") ax.set_title("Posterior of Mean Difference")
ax.plot([], [], ' ', label=f"Cohen's d = {d_cohen:.2f}\nProb sup = {ps:.2f}") ax.plot([], [], ' ', label=f"Cohen's d = {d_cohen:.2f}\nProb sup = {ps:.2f}")
ax.legend(loc=1) ax.legend(loc=1)
``` ```
%% Output %% Output
<matplotlib.legend.Legend at 0x7f2ecaa8ad40> <matplotlib.legend.Legend at 0x7f2ecaa8ad40>
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment