--- title: "Die Varianzanalyse" author: "Jonathan Harrington" date: "04. Juni 2019" output: html_document: default word_document: default pdf_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) library(dplyr) library(ez) source(file.path(pfadu, "phoc.txt")) source(file.path(pfadu, "ph.step.R")) blang = read.table(file.path(pfadu, "blang.txt")) v.df = read.table(file.path(pfadu,"vokal.txt")) dg = read.table(file.path(pfadu, "dg.txt")) ssb = read.table(file.path(pfadu, "ssb.txt")) ``` # Vorbemerkung zur sogenannten "wissenschaftlichen Notation" von Zahlen Wenn ein Effekt eines Faktors sehr groß ist, wird der p-Wert sehr klein, z.B. etwas wie `0.0000000147`. Ganz generell ist eine so kleine Zahl sehr schwer zu lesen (und dasselbe gilt auch für sehr große Zahlen). Daher verwendet R die sogenannte "wissenschaftliche Notation" für sehr große und sehr kleine Zahlen. Die wissenschaftliche Notation beschreibt jede Zahl als Kommazahl $m$, die größer als 0, aber kleiner als 10 ist, multipliziert mit $10 ^e$, also $$m * 10^e$$ So ist z.B. die Zahl eine Milliarde, wenn Sie sie in der `R`-Konsole eintippen (also als 1 mit 9 Nullen), wie folgt dargestellt, ```{r} 1000000000 ``` also als $1 * 10^{9}$. Ganz genauso ist es mit einer sehr kleinen Zahl, wie z.B. einem Milliardstel, ```{r} 0.000000001 ``` das also als $1 * 10^{-9}$ dargestellt wird (denn es ist ja auch eine 1 an der neunten Nachkommastelle). Normalerweise stellt `R` nur Zahlen in dieser Schreibweise dar, die wirklich groß bzw. klein (und damit unleserlich) sind: ```{r} 0.01 0.001 0.0001 0.00001 #usw. ... ``` Einige Packages (wie das heute zu verwendende `ez`) stellen aber in den Ergebnistabellen alle p-Werte in dieser Notation dar, so dass aus `0.18` ein `1.8e-01` wird. Damit man trotzdem auf einen Blick erkennt, was signifikant ist und was nicht, also ob dieser Wert von $p$ kleiner oder größer als 0.05 ist, wird in `ez` (und anderen Packages) zusätzlich ein `*` angegeben, wenn $p < 0.05$. # $ANOVA$ und $t$-Test ## $t$-Test oder $ANOVA$ ("Analysis of Variance" = Varianzanalyse) Ein $t$-Test *oder* eine $ANOVA$ sind dann äquivalent, wenn man nur **einen** Faktor hat, der nur **zwei Stufen** hat, z.B. - Hat Geschlecht einen Einfluss auf die Dauer? ## $ANOVA$ Hat man hingegen einen **Faktor mit mehr als 2 Stufen**, oder **mehr als einen Faktor**, so kann man keinen $t$-Test mehr anwenden, und braucht stattdessen eine Varianzanalyse. Beispiele wären: - Es gibt 3 *Alter*sgruppen, *jung*, *mittel*, *alt*. Hat die Altersgruppe einen Einfluss auf die Dauer? (= ein Faktor mit 3 Stufen) - Haben Geschlecht und Dialekt einen Einfluss auf die Dauer? (= 2 Faktoren) ## Was ist die Varianzanalyse? Die *Varianz* ist ein Streuungsmaß, das sich berechnet als Summe der quadrierten Abweichungen aller Einzelwerte einer Verteilung vom arithmetischen Mittel derselben Verteilung geteilt durch die Gesamtzahl der Werte. Wir kennen dies Maß indirekt bereits: die uns schon bekannte Stichprobenstandardabweichung (`sd()`) ist nämlich die Quadratwurzel der Varianz. Mit der Varianzanalyse wird (durch einen sogenannten $F$-Test) ein Verhältnis zwischen zwei Varianzen berechnet: derjenigen **innerhalb von Stufen** und derjenigen **zwischen Stufen**. ### Beispiel: F1 von drei Vokalkategorien, /ɪ,ɛ,a/ #### Innerhalb der Stufen Es gibt eine **randomisierte Variation von F1** *innerhalb* jeder Stufe (F1 von /ɪ/ variiert, F1 von /ɛ/ variiert, F1 von /a/ variiert). #### Zwischen den Stufen F1 variiert, weil es eine **systematische Variation** zwischen den Verteilungen der Vokalkategorien gibt: die Werte von /ɪ/, /ɛ/, und /a/ liegen in ganz unterschiedlichen F1-Bereichen, und je unterschiedlicher sie sind, umso **größer wird diese Varianz im Verhältnis zu der willkürlichen, randomisierten Varianz innerhalb der Stufen sein**: $$ F = \frac{systematische-Varianz-ZWISCHEN-den-Stufen}{randomisierte-Varianz-INNERHALB-der-Stufen} $$ ### Innersubjektfaktoren (within subject factors) Man spricht von Innersubkjektfaktoren, wenn, wie beim gepaarten $t$-Test auch, Analysen innerhalb derselben Person vorgenommen wurden, so dass zumindest eine Messung pro Versuchsperson pro Stufe vorliegt: #### Gepaarter $t$-Test *oder* $ANOVA$ mit Innersubjektfaktor /ba, pa/ wurden von denselben Sprechern produziert. *Unterscheiden sich /ba, pa/ in VOT?* - Abhängige Variable: **VOT** - Within-Faktor: **Stimmhaftigkeit** (= */ba/* vs. */pa/*) Wir haben hier genau einen Wert für /ba/ und einen Wert für /pa/ pro Versuchsperson. Wir können also einen gepaarten $t$-Test durchführen, oder auch eine $ANOVA$ mit dem Innersubjektfaktor **Stimmhaftigkeit**. #### Kein gepaarter $t$-Test möglich, aber $ANOVA$ mit Innersubjektfaktor /i, y, u/ wurden von denselben Sprechern produziert. *Unterscheiden sich /i, y, u/ in F2?* - Abhängige Variable: **F2** - Within-Faktor: **Vokal** (*/i/* vs. */y/* vs. */u/*) Also ein Wert für /i/ ein Wert für /y/, ein Wert für /u/ pro Versuchsperson. Wegen der drei Stufen ist kein gepaarter $t$-Test möglich, aber eine $ANOVA$ mit dem Innersubjektfaktor **Vokal**. ### Zwischensubjektfaktoren (Between subjects factors) Ein Zwischensubjektfaktor beschreibt meistens eine Eigenschaft der Versuchsperson, z.B. Muttersprache (englisch oder deutsch oder französisch), Geschlecht (m oder w), Alter (jung oder alt) usw. . In Analogie zum ungepaarten $t$-Test haben wir es also hier mit entweder-oder-Fällen zu tun; eine Versuchsperson ist z.B. *entweder* männlich *oder* weiblich. Natürlich können Zwischensubjektfaktoren bei einer $ANOVA$ aber auch mehr als zwei Stufen haben. ### Beispieldaten zur Unterscheidung zwischen Inner- und Zwischensubjektfaktoren Die Kieferposition wurde in 3 **Vokalen** */i, e, a/* und jeweils zu 2 **Sprechtempi** (*langsam*, *schnell*) gemessen. Die Messungen (3 x 2 = 6 pro Vpn) sind von 16 Vpn erhoben worden, 8 mit **Muttersprache** *Spanisch*, 8 mit Muttersprache *Englisch*. *Inwiefern haben die Faktoren Sprache, Sprechtempo, und Vokal einen Einfluss auf die Kieferposition?* #### Between - Sprache #### Within - Sprechtempo - Vokal # Within-Factor: gepaarter $t$-Test und $ANOVA$ ```{r} head(blang); dim(blang) ``` *12 Versuchspersonen produzierten /i/ in betonter und unbetonter Position. Hat Betonung einen Einfluss auf F2?* Jede Stufe des unabhängigen within-Faktor wird von jeder Versuchsperson einmal belegt: ```{r} with(blang, table(Vpn, Betonung)) ``` ## Drei Schritte ### 1. Differenz-Berechnung ```{r} d = blang%>% group_by(Vpn)%>% summarise(F2 = diff(F2)) d ``` ### 2. Boxplot ```{r} boxplot(d$F2) ``` ### 3. Statistischer Test #### Gepaarter $t$-Test ```{r} t.test(d$F2) ``` **F2 wird signifikant von Betonung beeinflusst (`t[11] = 4.4, p < 0.01`). ** #### Oder: Lösung mit $ANOVA$ (`ezANOVA()` aus dem package `ez`, angewendet auf den ursprünglichen data.frame) ```{r} ezANOVA(blang, .(F2), .(Vpn), .(Betonung)) ``` **F2 wird signifikant von Betonung beeinflusst (`F[1,11] = 19.0, p < 0.01`).** ## Generelle Syntax des Befehls `ezANOVA()` ```{r, eval=FALSE} ezANOVA(DATAFRAME, .(Abhängige numerische Variable), .(Versuchspersonen), .( WITHIN-FAKTOR(EN) ), between = .( BETWEEN-FAKTOR(EN) )) ``` ## Vergleich $ANOVA$ und $t$-Test-Ergebnisse ### Derselbe Wahrscheinlichkeitswert: ```{r} t.test(d$F2)$p.value ezANOVA(blang, .(F2), .(Vpn), .(Betonung))$ANOVA$p ``` ### Zusammenhang zwischen den Prüfstatistikwerten: $F = t^2$ ```{r, message=FALSE} t.test(d$F2)$statistic t.test(d$F2)$statistic^2 ezANOVA(blang, .(F2), .(Vpn), .(Betonung))$ANOVA$F ``` # Between-factor: Two-sample $t$-Test und $ANOVA$ ```{r} head(v.df); dim(v.df) ``` **Unterscheiden sich Deutsche und Engländer in F2 von /e/? ** Jede Stufe des unabhängigen between-Faktors wird einmal pro Versuchsperson belegt (between-Faktor: Die Versuchspersonen sind entweder Deutsch oder Englisch): ```{r} with(v.df, table(Vpn, Sprache)) ggplot(v.df) + aes(y = F2, x = Sprache) + geom_boxplot() ``` Oder: ```{r} ggplot(v.df) + aes(x = F2, col = Sprache) + geom_density() ``` $t$-Test: ```{r} t.test(F2 ~ Sprache, data = v.df) ``` **F2 wird signifikant von der Sprache beeinflusst (`t[11.8] = 2.7, p < 0.05`).** Oder: $ANOVA$ mit between-subjects-factor: ```{r, message=FALSE} ezANOVA(v.df, .(F2), .(Vpn), between = .(Sprache)) ``` **F2 wird signifikant von der Sprache beeinflusst (`F[1,18] = 7.2, p < 0.05`).** Eine Anova mit between-Faktor wird unter der Annahme durchgeführt, dass sich die Varianzen der Stufen nicht unterscheiden. Daher wird der Levene-Test (wenn p > 0.05, dann ist die Anova berechtigt) automatisch durchgeführt. Insofern bekommt man das gleiche Ergebnis mit einem $t$-test nur unter der Annahme, dass sich die Varianzen in den Stufen nicht unterscheiden: ```{r,message=FALSE} t.test(F2 ~ Sprache, data = v.df, var.equal=T) t.test(F2 ~ Sprache, data = v.df, var.equal=T)$statistic t.test(F2 ~ Sprache, data = v.df, var.equal=T)$statistic^2 ezANOVA(v.df, .(F2), .(Vpn), between = .(Sprache))$ANOVA$F ``` # Zwei Faktoren ```{r} head(dg); names(dg) ``` **Inwiefern wird F2 vom Dialekt *und* vom Geschlecht beeinflusst?** ```{r} with(dg, table(Vpn, interaction(Region, Gen))) ``` ### Between - Gender ### Within - Region ## Bei 2 Faktoren gibt es immer 3 Fragen: - **Frage zu *Faktor 1*:** Hat *Gender* einen Einfluss auf F2? - **Frage zu *Faktor 2*:** Hat *Region* einen Einfluss auf F2? - **Frage zur *Interaktion*:** Gibt es eine *Interaktion zwischen Region und Gender*? = Ist der Unterschied zwischen männlich und weiblich derselbe in allen 3 Regionen? ```{r} ggplot(dg) + aes(y = F2, x = Gen, colour = Region) + geom_boxplot() ``` **Interaktion**: *ist der Unterschied zwischen männlich und weiblich ähnlich in den 3 Regionen?* Wenn ja, dann müsste der Abstand zwischen den *männlich*-*weiblich*-Medianen ähnlich sein (d.h. diese Linien müssten mehr oder weniger parallel zu einander sein): ```{r, message=FALSE} # zunächst den Median pro Faktorstufenkombination ermitteln: dg.m = dg %>% group_by(Gen,Region) %>% summarise(F2 = median(F2)) # diese Medianwerte dann plotten ggplot(dg.m) + aes(y = F2, x = Region, group=Gen, colour = Gen) + geom_line() ezANOVA(dg, .(F2), .(Vpn), between =.(Region, Gen)) ``` **F2 wurde signifikant von der Region (`F[2,54] = 119.6, p < 0.001`) und von Geschlecht (`F[1,54] = 106.1, p < 0.001`) beeinflusst, und es gab eine signifikante Interaktion zwischen diesen Faktoren (`F[2,54] = 12.1, p < 0.001).** # Post-hoc $t$-Tests Wenn eine Interaktion vorliegt, sollte durch $t$-Tests geprüft werden, ob sich alle Paare von Stufen-Kombinationen in der abhängigen Variable (hier F2) unterscheiden. Die Anzahl dieser Tests berechnet sich nach: - Region: 3 Stufen - Geschlecht: 2 Stufen $= 3 *2 = 6$ Stufen. Alle Paarungen davon: ```{r} factorial(6)/(factorial(4) * factorial (2)) ``` $6!/(4! × 2!) = 15$ Testpaare! Post-hoc-t-Tests lassen sich durchführen mittels der Funktion `phoc()`, die wir ganz am Anfang des Skriptes geladen hatten (mit `source(file.path(pfadu, "phoc.txt"))`): ```{r} phoc(dg, .(F2), .(Vpn), .(Region, Gen)) ``` - `$res`: die Ergebnisse der t-tests - `$name`: die Testpaare - `$paired`: ob ein gepaarter oder ungepaarter t-test durchgeführt wurde - `$bonf`: Anzahl der durchgeführten Tests (= Bonferroni-Faktor) ### Bonferroni-Korrektur `prob-adj` in `$res` ist der durch die sogenannte Bonferroni-Korrektur angepasste $p$-Wert. Hierbei wird der Wahrscheinlichkeitswert der inviduellen $t$-Tests mit der Anzahl der theoretisch möglichen Testkombinationen (= der Bonferroni-Faktor: 15 in diesem Fall) multipliziert. Sollte dies einen Wert > 1 ergeben, wird der Wert einfach auf 1 (= 100%) gesetzt. Der Grund liegt im Problem des multiplen Testens in derselben Stichprobe: Je mehr post-hoc Tests durchgeführt werden, um so wahrscheinlicher ist es, dass einer von den vielen Tests durch reinen Zufall signifikant sein wird. Die Bonferroni-Korrektur ist eine Maßnahme dagegen. ### Auswahl relevanter post-hoc $t$-Tests Nicht alle $t$-Tests werden benötigt, sondern eher nur *Vergleiche zwischen Stufen* von **einem** Faktor, wenn die Stufen aller anderen Faktoren **konstant** sind. So interessiert es uns beispielsweise *nicht*, ob sich die Werte von Männern aus der Region A von den Werten der Frauen aus Region B unterschieden oder nicht. ```{r} vok.ph = phoc(dg, .(F2), .(Vpn), .(Region, Gen)) ``` - **Unterscheiden sich die Regionen innerhalb desselben Geschlechts? (Region variiert, Geschlecht ist konstant).** ```{r} phsel(vok.ph$res, 1) ``` - **Unterscheiden sich Männer und Frauen in derselben Region? (Geschlecht variiert, Region ist konstant).** ```{r} phsel(vok.ph$res, 2) ``` - **Aber *nicht*, wenn beide Faktoren variieren.** #### Runden der adjustierten $p$-Werte Noch ersichtlicher werden die Signifikanzen der relevanten Tests, wenn auf z.B. 3 Zahlen aufgerundet wird: ```{r} p1 = phsel(vok.ph$res, 1) round(p1, 3) p2 = phsel(vok.ph$res, 2) round(p2, 3) ``` **Post-hoc Bonferroni-adjusted $t$-Tests zeigten signifikante F2-Unterschiede zwischen A vs C (`p < 0.001`) und zwischen B vs C (`p < 0.001`), jedoch nicht zwischen A vs. B. F2 von Männern und Frauen unterschieden sich signifikant für die Regionen A (`p < 0.001`) und B (`p < 0.001`), jedoch nicht für C.** # Zwei Bedingungen für die Durchführung der Varianzanalyse ## Wiederholungen in within-Stufen Wenn es $n$ within-Stufen gibt, dann müssen es $n$ Werte pro Vpn sein, also pro Vpn ein Wert pro within-Stufe z.B: *Englische und spanische Vpn produzierten /i, e, a/ zu 2 Sprechgeschwindigkeiten* **Within**: *Vokal* (3 Stufen) und *Sprechgeschwindigkeit* (2 Stufen) Daher: $3 × 2 = 6$ within-Werte pro Versuchsperson (ein Wert pro within-Stufe pro Vpn). Jedoch haben die meisten phonetischen Untersuchungen mehrere Werte pro within-Stufe, da in diesen Untersuchungen häufig Wiederholungen eingesetzt werden. Als Beispiel: jede Vpn. erzeugte /i, e, a/ zu einer langsamen und schnellen Sprechgeschwindigkeit, und das jeweils 10 Mal. Wiederholungen in derselben within-Stufe sind in einer $ANOVA$ aber nicht zulässig und müssen gemittelt werden – damit wir pro Vpn. einen Wert pro within-Stufe haben (6 Mittelwerte pro Vpn. in diesem Beispiel). ```{r} summary(ssb) ``` In einer Untersuchung zur /u/-Frontierung im Standardenglischen wurde von 12 Sprecherinnen (6 *alt*, 6 *jung*) F2 zum zeitlichen Mittelpunkt in drei verschiedenen /u/-Wörtern erhoben (*used*, *swoop*, *who'd*). Jedes Wort ist von jeder Vpn. 10 Mal erzeugt worden. **Inwiefern wird F2 vom Alter und Wort beeinflusst?** Wieviele Werte pro Vpn. dürfen in der $ANOVA$ vorkommen? 3! Wieviele Werte darf es insgesamt in der $ANOVA$ geben? 36! ```{r} dim(ssb) with(ssb, table(Vpn, interaction(Wort, Alter))) ``` Über die Wort-Wiederholungen pro Versuchsperson mit `group_by()` mitteln ("Alter" steht mit in `group_by()`, um die Altersinformation nicht zu verlieren): ```{r} ssbm = ssb %>% group_by(Wort,Alter,Vpn) %>% summarise(F2 = mean(F2)) dim(ssbm); head(ssbm) with(ssbm, table(Vpn, interaction(Wort, Alter))) ``` Erst jetzt, mit den gemittelten Daten, sollte man eine Abbildung erstellen: ```{r} ggplot(ssbm) + aes(y = F2, x = Alter, colour = Wort) + geom_boxplot() ezANOVA(ssbm, .(F2), .(Vpn), .(Wort), between = .(Alter)) ``` ## Sphericity-Korrektur *Sphericity* ist die Annahme, dass die Unterschiede zwischen den Stufen eines within-Faktors **dieselbe Varianz haben**. Wenn Sphericity nicht gegeben ist, werden **die Wahrscheinlichkeiten durch Änderungen in den Freiheitsgraden nach oben** gesetzt. Dieses Problem kommt *nur dann* vor, wenn **ein within-Faktor mehr als 2 Stufen** hat. **Man soll grundsätzlich immer für Sphericity korrigieren**, wenn `Sphericity Corrections` in der Ausgabe von `ezANOVA() erscheint. Die betroffenen Freiheitsgrade werden mit - dem *Greenhouse-Geisser-Epsilon* (GGe) multipliziert, wenn es unter 0.75 liegt, - sonst (also wenn das GGe > 0.75) mit dem *Huynh-Feldt-Epsilon* (HFe): - sollte in diesem letzten Fall das HFe > 1 sein, dann einfach die ursprünglichen Freiheitsgrade nehmen d.h. keine Korrektur einsetzen. Für ```{r} ezANOVA(ssbm, .(F2), .(Vpn), .(Wort), between = .(Alter)) ``` Wort: F[2,20] ➞ F[2 * 0.6860511, 20 * 0.6860511] = F[1.4, 13.7] Alter × Wort Interaktion: F[2,20] ➞ F[1.4, 13.7] Die neuen damit verbunden Wahrscheinlichkeiten sind `p[GG]` (wenn mit `GGe` multipliziert wurde), sonst `p[HF]`. Das sind die Wahrscheinlichkeiten mit den korrigierten Freiheitsgraden (vergleiche mit p[GG] für Wort und die Alter:Wort-Interaktion): ```{r} 1 - pf(78.505534, 2 * 0.6860511, 20 * 0.6860511) 1 - pf(9.8908882, 2 * 0.6860511, 20 * 0.6860511) ``` **Alter (`F[1,10] = 14.9, p < 0.001`), Wort (`F[1.4 , 13.7] = 78.5, p < 0.001`) sowie die Interaktion von Wort und Alter (`F[1.4, 13.7] = 9.9, p < 0.01`) hatten einen signifikanten Einfluss auf F2.**