--- title: "Die Varianzanalyse" author: "Jonathan Harrington" date: "29. Mai 2018" output: html_document --- ```{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")) ``` #$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? Mit der Varianzanalyse wird (durch einen $F$-Test) ein Verhältnis zwischen zwei Varianzen berechnet: **innerhalb von Stufen** und **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{randomisierte-Varianz-zwischen-den-Stufen}{systematische-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$ /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. #### Kein gepaarter $t$-Test möglich, aber $ANOVA$ /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. ### 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. . ### 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} 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() ggplot(v.df) + aes(x = F2, col = Sprache) + geom_density() t.test(F2 ~ Sprache, data = v.df) ``` **F2 wird signifikant von der Sprache beeinflusst (`t[11.8] = 2.7, p < 0.05`).** ```{r} 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} 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} #zunächst den Median pro FAktorstufenkombination ermitteln: dg.m = dg %>% group_by(Gen,Region) %>% summarise(F2 = median(F2)) 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: $6!/(4! × 2!) = 15$ Testpaare: ```{r} factorial(6)/(factorial(4) * factorial (2)) ``` ```{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. Der Grund: 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. ```{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, 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. zB. jede Vpn. erzeugte /i, e, a/ zu einer langsamen und schnellen Sprechgeschwindigkeit jeweils 10 Mal. Wiederholungen in derselben within-Stufe sind in einem ANOVA 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 insgesamt in der $ANOVA$ wird es geben? 36! ```{r} with(ssb, table(Vpn, interaction(Wort, Alter))) ``` Über die Wort-Wiederholungen mit `group_by()` mitteln: ```{r} ssbm = ssb %>% group_by(Wort,Alter,Vpn) %>% summarise(F2 = mean(F2)) dim(ssbm); head(ssbm) with(ssbm, table(Vpn, interaction(Wort, Alter))) 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* multipliziert, wenn es unter 0.75 liegt, - sonst mit dem *Huynh-Feldt-Epsilon*: - sollte in diesem letzten Fall das H-F-Epsilon > 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: ```{r} 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.**