--- title: "Die t-Verteilung und der t-Test" author: "Jonathan Harrington" date: "21 Mai 2019" output: word_document: default html_document: default --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ```{r,echo = FALSE} pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf" ``` ```{r,message=FALSE} library(ggplot2) source(file.path(pfadu, "proben.R")) form = read.table(file.path(pfadu, "bet.txt")) e.df = read.table(file.path(pfadu, "e.txt")) ``` ## 1. SE (Standard Error) und Konfidenzintervall Zur Erinnerung: Der Standard Error (SE) ist die Populationsstandardabweichung von Mittelwerten. Wir werfen beispielsweise 12 Würfel gleichzeitig, und berechnen den Mittelwert der Zahlen. Der SE lässt sich berechnen mit $σ/√k$: - $σ$ ist die Populationsstandardabweichung = `sd(1:6) * sqrt(5/6)` - $k$ ist die Anzahl der Würfel = 12 ```{r} SE = sd(1:6) * sqrt(5/6) / sqrt(12) SE ``` Die Bedeutung hiervon ist: wenn wir 12 Würfel gleichzeitg werfen, davon den Mittelwert berechnen, diesen Vorgang unendlich viel Mal wiederholen, sodass wir unendlich viele Mittelwerte hätten, dann davon die Stichprobenstandardabweichung berechnen, dann wäre diese Standardabweichung genau `SE = sd(1:6) * sqrt(5/6) / sqrt(12)`. Wir müssten ziemlich nah an diesen SE mit z.B. 50000 Mittelwerten kommen: ```{r} o = proben(k=12, N = 50000) sd(o) ``` ### Konfidenzintervall Wir wollen zwei Werte a und b auf eine solche Weise berechnen, sodass der Mittelwert zwischen a und b mit einer Wahrscheinlichkeit von 0.95 liegt. Hier benötigen wir den Populations-Mittelwert: ```{r} mu = mean(1:6) # und SE SE = sd(1:6) * sqrt(5/6) / sqrt(12) # a und b a = qnorm(0.025, mu, SE) b = qnorm(0.975, mu, SE) a b ``` Das heißt also, dass, wenn wir 12 Würfel zusammen werfen und davon den Mittelwert berechnen, der Mittelwert dann mit einer Wahrscheinlichkeit von 0.95 zwischen (a) 2.533725 und (b) 4.466275 (= ein Mittelwert von weniger als 2.533725 oder mehr als 4.466275 wird meistens nur in 5/100 Fällen vorkommen) fällt. Prüfen wir das: 12 Würfel werfen, davon den Mittelwert berechnen, diesen Vorgang 100 Mal wiederholen (daher 100 Mittelwerte): ```{r} m = proben(k = 12, N = 100) ``` Wieviele dieser Mittelwerte sind kleiner als 2.533725 oder größer als 4.466275? ```{r} sum(m < 2.533725 | m > 4.466275) ``` ## 2. Konfidenzintervall einer Stichprobe Hier sind 12 Dauerwerte von einem /a:/ Vokal: ```{r} d = c(119, 111, 105, 130, 133, 122, 124, 129, 95, 100, 109, 111) ``` Wir wollen aufgrund der Stichprobe ein Konfidenzintervall für den Dauer-Mittelwert von /a:/ erstellen (= dieser Mittelwert fällt zwischen a und b mit einer Wahrscheinlichkeit von 95%). Zu diesem Zweck müssen ... - i. (mu, SE) aufgrund der Stichprobe wie folgt eingeschätzt werden, angenommen dass es sich um eine randomisiert ausgewählte Stichprobe handelt: beste Einschätzung von mu: ```{r} mu = mean(d) ``` Die beste Einschätzung von SE: ```{r} SE = sd(d)/sqrt(12) ``` - ii. Wenn (mu, SE) nicht aus theoretischen Überlegungen (wie beim Würfel-Spiel) heraus berechnet werden können, dann wird das Konfidenzintervall *nicht durch die Normalverteilung*, sondern durch die **t-Verteilung** mit einer gewissen Anzahl von sogenannten *Freiheitsgraden* erstellt. Die t-Verteilung ist der Normalverteilung sehr ähnlich, und nähert sich der Normalverteilung zunehmend an, je höher die Anzahl der Freiheitsgrade ist, z.B. ```{r} # t-Verteilung # mit 3, 6, 20 Freiheitsgrade auf die Normalverteilung überlagern: ylab = "Wahrscheinlichkeitsdichte" curve(dnorm(x, 0, 1), xlim = c(-4, 4), ylab = ylab) # t-Verteilung mit 3 df curve(dt(x, 3), add = T, col = "red") # t-Verteilung mit 6 df curve(dt(x, 6), add = T, col = "blue") # t-Verteilung mit 20 df curve(dt(x, 20), add = T, col = "green") ``` Die Anzahl der **Freiheitsgrade** entspricht der **Anzahl der Stichproben minus 1**: ```{r} df = 11 ``` Das Konfidenzintervall für eine Stichprobe wie die obige wird mit `qt()` statt mit `qnorm()` berechnet: ```{r} a = mu + SE * qt(0.025, df) b1 = mu - SE * qt(0.025, df) # oder b2 = mu + SE * qt(0.975, df) a b1 b2 ``` Bedeutung: aufgrund dieser Stichprobe fällt der Dauermittelwert für /a/ zwischen 107.8 ms und 123.5 ms mit einer Wahrscheinlichkeit von 0.95 (also 95%). ## 3. Konfidenzintervall für den Unterschied zwischen 2 Stichproben Die Dauerwerte, wenn 12 Sprecher (Vpn1, Vpn2... Vpn12) ein betontes /a/ produzierten, waren wie folgt: ```{r} bet = c(119, 111, 105, 130, 133, 122, 124, 129, 95, 100, 109, 111) ``` Die Dauerwerte derselben 12 Sprecher (Vpn1, Vpn2... Vpn12), die ein unbetontes /a/ produzierten, waren: ```{r} un = c(110, 95, 108, 80, 120, 110, 120, 95, 72, 83, 90, 95) ``` Wir wollen ein 95%-Konfidenzintervall erstellen für den **Unterschied** zwischen den Dauern für *betont* und *unbetont*: ```{r} mu = mean(bet - un) SE = sd(bet - un)/sqrt(12) # 95% Konfidenzintervall mu + SE * qt(0.025, 11) mu + SE * qt(0.975, 11) ``` Das 95% Konfidenzintervall für den Unterschied zwischen den Dauerwerten ist 8.7 ms ≤ mu ≤ 26.3 ms. Bedeutung: der Unterschied zwischen dem Dauer-Mittelwert von einem betontem /a/ und dem Dauer-Mittelwert von einem unbetonten /a/ liegt zwischen 8.7 ms und 26.3 ms mit einer Wahrscheinlichkeit von 0.95. ## 4. Prüfen einer Hypothese, oder: unterscheiden sich Mittelwerte *signifikant*? $H0$ (die Null-Hypothese): - "Betonung hat keinen Einfluss auf die Dauer." + Bedeutung: der Unterschied zwischen den Mittelwerten ist Null (0). $H1$ (die Alternativ-Hypothese): - "Betonung beeinflusst die Dauer."" + Bedeutung: der Unterschied zwischen den Mittelwerten weicht von Null ab. $α$-Wert (alpha-Wert): - 0.05 ist der hier gewählte, sogenannte $α$-Wert (alpha-Wert), bei dem wir $H0$ *verwerfen* + Ein Konfidenzintervall von 0.95 bedeutet einen $α$-Wert von 1 - 0.95 = 0.05. **Prüfen:** wenn 0 außerhalb des Konfidenzintervalls 8.7 ms ≤ mu ≤ 26.3 ms fällt, verwerfen wir $H0$ und akzepierten $H1$. Schlussfolgerung für das obige Beispiel: Wir verwerfen $H0$ und akzeptieren $H1$ ("Betonung beeinflusst die Dauer."). **Berichten:** "Die Dauer wird signifikant von der Betonung beeinflusst (p < 0.05)." *Bedeutung (i)*: - die Wahrscheinlichkeit, dass die Dauer *nicht* von der Betonung beeinflusst wird, ist weniger als 0.05 (weniger als 5%). *Bedeutung (ii)*: - die Wahrscheinlichkeit, dass der Unterschied zwischen dem Dauer-Mittelwert vom betonten /a/ und dem Dauer-Mittelwert vom unbetonten /a/ 0 (Null) sein könnte, ist kleiner als 0.05. ## 5. Der t-Test Wir bekommen dieselbe Auswertung einfacher durch den sogenannten **t-Test**. Manuell noch einmal: ```{r} mu = mean(bet - un) SE = sd(bet - un)/sqrt(12) # 95% Konfidenzintervall mu + SE * qt(0.025, 11) mu + SE * qt(0.975, 11) ``` Mit dem t-Test: ```{r} t.test(bet - un) ``` Identifizieren: ```{r} mu ``` - μ (`mean of x: 17.5`) - Freiheitgrade (`df = 11`) - 95%-Konfidenzintervall (`8.739938` --- `26.260062`) - `t = 4.3969` (die $t$-Statistik). + Der Abstand zwischen μ und 0 (Null) in `SE`-Einheiten. ```{r} (mu - 0)/SE ``` - `p-value = 0.001069` + Bedeutung (i): Die Wahrscheinlichkeit, dass der Unterschied zwischen den Mittelwerten Null sein könnte. + Bedeutung (ii): Die Wahrscheinlichkeit, dass $H0$ zutrifft. + Bedeutung (iii): die Wahrscheinlichkeit, dass der Unterschied zwischen den Mittelwerten außerhalb des Konfidenzintervalls fällt, auch: ```{r} (1- pt(4.396914, 11)) * 2 ``` ### Ergebnis berichten: Wir wählen immer drei $α$-Werte aus: $p < 0.05$, $p < 0.01$, $p < 0.001$ und wählen den $α$-Wert, der am nächsten über dem $p$-Wert liegt Hier wählen wir $p < 0.01$, da der $p$-Wert ($0.001069$) über 0.001 aber unter 0.01 liegt. Wir berichten entweder: - "Betonung hatte einen signifikanten Einfluss auf die Dauer ($t[11] = 4.4, p < 0.01$)" Oder: - "Die Betonung wurde signifikant von der Dauer beeinflusst ($t[11] = 4.4, p < 0.01$)" Sollte $p > 0.05$ sein, dann ist das Ergebnis nicht signifikant ("n.s.") (Wir verwerfen nicht H0), und schreiben: - "Die Betonung hatte keinen signifikanten Einfluss auf die Dauer" oder - "Die Dauer wurde nicht signifikant von der Dauer beinflusst." ## 6. Gepaarter t-Test Das obige Beispiel ist ein sogenannter *gepaarter t-Test*, weil eine Differenz --- ob betont oder unbetont --- pro Paar berechnet wird (daher 12 Paare in dem obigen Beispiel). Ein gepaarter t-test kommt in der Phonetik meistens dann vor, wenn Stichproben-Paare pro Versuchsperson verglichen werden. Ein zweites Beispiel, bei dem (wie üblich) die Werte in einem Data-Frame stecken: - *12 Versuchspersonen produzierten jeweils ein betontes und ein unbetontes /i/. Unterscheiden sich das betonte und das unbetonte /i/ voneinander in F2?* ```{r} dim(form) head(form) summary(form) ``` Die Frage bitte immer umstellen ("Wird y von x beeinflusst?"): - *"Wurde F2 (abhängige Variable) von der Betonung (unabhängiger Faktor mit 2 Stufen: betont/unbetont) beeinflusst?"* Der Test ist gepaart: es gibt jeweils ein Wertepaar pro Versuchsperson. ### 6.1. Boxplot #### 6.1.1. Unterschiede pro Paar (hier Versuchsperson) berechnen Entweder mit `aggregate()`; dort in der Formel den unterscheidenden Faktor weglassen, und als Funktion `diff()` angeben: ```{r} unterschied = aggregate(F2 ~ Vpn, FUN = diff, data = form) unterschied ``` Oder mit dem Package `dplyr`: ```{r,message=FALSE} library(dplyr) unterschied = form%>% group_by(Vpn)%>% summarise(F2=diff(F2)) unterschied ``` #### 6.1.2. Boxplot der Unterschiede ```{r} # am einfachsten ohne ggplot2: boxplot(unterschied$F2) # mit ggplot2 (x einfach weglassen): ggplot(unterschied) + aes(y = F2) + geom_boxplot() #oder: ggplot(unterschied) + aes(x = "betont-unbetont", y = F2) + geom_boxplot() + xlab("Unterschiede in F2") ``` ### 6.2. `t.test()` Der t-Test prüft, ob der Mittelwert der Unterschiede signifikant von 0 abweicht oder nicht: ```{r} t.test(unterschied$F2) ``` #### Antwort: - "F2 wurde signifikant von der Betonung beeinflusst ($t[11] = 4.4, p < 0.01$)." ##### Alternative Berechnung - **y ~ x** unter der Berüchsichtigung, dass $x$ (hier `Bet`) im ursprünglichen data.frame `form` gepaarte Werte für $y$ (hier `F2`) kodiert hat: ```{r} t.test(F2 ~ Bet, paired = T, data = form) ``` ##### Alternative Abbildung Manche Leute bevorzugen auch, beide Datensätze als Punkte abzubilden, und die Paarungen als Linienverbindungen darzustellen: ```{r} ggplot(form) + aes(x = Bet, y = F2,col = Vpn, group = Vpn) + geom_point() + geom_line() ``` In unserem Kurs wollen wir aber lieber bei der Abbildung mit einer Box, die die Paar-Differenz-Werte statistisch-deskriptiv beschreibt, bleiben! ## 7. Ungepaarter t-test Ein *ungepaarter t-Test* liegt vor, wenn nicht Paare von Stichproben, sondern zwei Gruppen miteinander verglichen werden. In der Phonetik werden z.B. oft zwei verschiedene Sprecher-Gruppen (männlich/weiblich; Bayern/Hessen; englisch/deutsch) verglichen. - *Unterscheiden sich deutsch und englisch in F2 von /e/?* + = Wird F2 (abhängige Variable) von der Sprache (unabhängige Variable mit 2 Stufen: englisch/deutsch) beeinflusst? ```{r} head(e.df) dim(e.df) ``` Vorgang: wie oben, aber die Stichproben werden nicht paarweise voneinander subtrahiert (da sie nicht gepaart sind). Daher ist eine genaue Übereinstimmung der Gruppengrößen auch nicht erforderlich (dennoch sollten die Gruppengrößen auch nicht zu sehr voneinander abweichen), hier z.B. ```{r} table(e.df$Sprache) ``` ```{r} # 1. Boxplot oder Densitplot der Unterschiede ggplot(e.df) + aes(y = F2, x = Sprache) + geom_boxplot() # density plot ggplot(e.df) + aes(x = F2, col = Sprache) + geom_density() ``` Hier prüfen wir, ob signifikante Unterschiede zwischen den Mittelwerten der beiden Gruppen vorliegen (NB: nicht gepaart!): ```{r} t.test(F2 ~ Sprache, data = e.df) ``` ##### Antwort: - "F2 wurde signifikant ($t[21.1] = 2.3, p < 0.05$) von der Sprache beeinflusst." (In einem ungepaarten Welch Two Sample t-Test bekommt man immer Gleitkommazahlen als Freiheitsgrade, da die Freiheitsgrade mittels einer Formel (der sogenannten Welch-Satterthwaite-Formel) verändert ("approximiert") werden. Durch diese Änderung müssen wir uns weniger Gedanken darüber machen, wann ein solcher Test erlaubt ist oder nicht.) # 8. Voraussetzung für den t-Test prüfen: `shapiro.test()` zur Überprüfung von Verteilungen Der `shapiro.test()` überprüft, ob die Verteilungen der Daten signifikant von einer Normalverteilung abweichen. Sollte dies der Fall sein, dürfen wir den t-Test nicht anwenden, sondern müssen statt dessen den Wilcoxon-Vorzeichen-Rang-Test anwenden (`wilcox.test()`). Ergibt der `shapiro.test()` keine signifikante Abweichung der Verteilungen unserer Stichproben von einer Normalverteilung, bleiben wir beim t-Test. Beispiel für den gepaarten Fall: ```{r, message = FALSE} shapiro.test(unterschied$F2) # da hier p < 0.05, ist der oben angewandte t-Test eigentlich nicht zulässig ... t.test(unterschied$F2) # ... und muss daher mit dem Wilcoxon signed rank test ersetzt werden: wilcox.test(unterschied$F2) ``` Wir müssen also unsere Antwort ändern zu: "F2 wurde signifikant von der Betonung beeinflusst, wie ein Wilcoxon-Test zeigte ($V = 3, p < 0.01$)." Im Beispiel für den ungepaarten Fall war der t-Test allerdings durchaus zulässig: ```{r, message = FALSE} with(e.df, tapply(F2, Sprache, shapiro.test)) #da beide nicht signifikant --> t.test(F2 ~ Sprache, data = e.df) #wäre zumindest ein Datensatz signifikant von einer Normalverteilung abweichend, wäre die Syntax für die Alternative gewesen: wilcox.test(F2 ~ Sprache, data = e.df) ``` Da der Shapiro-Test hier für keine der Datensätze eine Abweichung von der Normalverteil gezeigt hatte, können wir beim t-Test bleiben und schreiben auch weiterhin: "F2 wurde signifikant ($t[21.1] = 2.3, p < 0.05$) von der Sprache beeinflusst."