# Aehnlich wie in der Klausur: gemischte Fragen: # pfadu = "http://www.phonetik.uni-muenchen.de/~jmh/lehre/Rdf" # Anmerkung: es werden immer alle packages und selbstgeschriebene # Funktionen bereitgestellt werden. # Das heisst aber nicht notwendigerweise, dass Sie die alle brauchen! library(ggplot2) library(dplyr) library(ez) library(lmerTest) library(emmeans) source(file.path(pfadu, "phoc.txt")) source(file.path(pfadu, "sig.fn.R")) proben <- function(unten = 1, oben = 6, k = 10, N = 50) { # default: wir werfen 10 Wuerfel 50 Mal alle <- NULL for (j in 1:N) { ergebnis = mean(sample(unten:oben, k, replace = T)) alle = c(alle, ergebnis) } alle } #Dann die benoetigten Data.frames: tv = read.table(file.path(pfadu, "tv.df.txt")) fremd = read.table(file.path(pfadu, "fremd.df.txt")) glottal = read.table(file.path(pfadu, "glottal.txt")) vcv = read.table(file.path(pfadu, "vcvC.txt")) form10 = read.table(file.path(pfadu, "form10.df.txt")) svot = read.table(file.path(pfadu, "svot.txt")) stefan = read.table(file.path(pfadu, "stefan.txt")) # 1. Diese Daten: dim(tv) # zeigen eine Messung der Sprachkompetenz (cdi) fuer # 80 Schueler (Vpn) sowie die Anzahl # der Stunden am Tag, die sie fernsehen (tv.hours). # Pruefen Sie durch eine Abbildung und statistischen Test, # ob eine Verbindung zwischen Sprachkompetenz und # der Anzahl der Fernsehstunden vorliegt. plot(cdi ~ tv.hours, data = tv) # Trotz der Nennung von "80 Schuelern" ist hier kein Mixed Model noetig, # da fuer jede Versuchsperson nur eine einzige Beobachtung (also insgesamt 80) ... dim(tv) # ...vorhanden ist # (d.h., die Daten sind unabhaengig voneinander und koennen deshalb durch # Regression bearbeitet werden) reg = lm(cdi ~ tv.hours, data = tv) abline(reg) #oder ggplot(tv) + aes(y = cdi, x = tv.hours) + geom_point() + geom_abline(intercept = coef(reg)[1], slope = coef(reg)[2]) summary(reg) # Es gibt eine lineare Beziehung zwischen cdi und tv.hours (R^2 = 0.22, F[1,78] = 23.1, p < 0.001). shapiro.test(resid(reg)) # OK plot(resid(reg)) # OK acf(resid(reg)) # OK # 2. Die Daten in diesem Data-Frame: dim(glottal) # zeigen ob ein aspirierter (t) oder glottalisierter (tQ) # Verschluss von einem Erkennungssystem # richtig (ja) erkannt wurde oder nicht (nein). # Erstellen Sie eine Abbildung, und fuehren Sie # einen statistischen Test durch, um zu pruefen, # ob die Erkennung von dem Verschluss beeinflusst wurde. head(glottal) p1 = ggplot(glottal) p2 = aes(fill = erkannt, x = verschluss) p3 = geom_bar(position = "fill") p1 + p2 + p3 #beide Variablen (die abhaengige und die unabhaengige) sind kategorial, daher glm() g = glm(erkannt ~ verschluss, binomial, data = glottal) anova(g, test = "Chisq") # Die Erkennung wird vom Verschluss nicht signifikant beeinflusst, # obwohl es eine ganz leichte Tendenz dazu gibt. # 3. Diese Daten: dim(fremd) # zeigen fuer 30 Sprecher (Vpn) eine Messung der Sprechgeschwindigkeit (tempo), # wenn sie in ihrer Muttersprache oder ihrer Zweitsprache (Sprache) reden. # Pruefen Sie durch eine Abbildung und statistischen Test, # inwiefern die Sprechgeschwindigkeit von der Sprache beeinflusst wird. ggplot(fremd) + aes(y = tempo, x = Sprache) + geom_boxplot() #Fuer den Test theoretisch drei Moeglichkeiten: #t.test(), ezANOVA(), oder lmer(); fuer letzteres ist die Anzahl der Beobachtungen #sehr gering (dim(fremd) zeigt 60 ) with(fremd, table(Vpn, Sprache)) # Sprache ist within; es gibt keine fehlenden Werte und keine Wiederholungen; # eigentlich ein Fall fuer einen gepaarten t-test (den man auf zwei Weisen realisieren kann); # der Vollstaendigkeit halber: man kann das auch als ezANOVA() mit einem within-Faktor berechnen: ezANOVA(fremd, .(tempo), .(Vpn), .(Sprache)) # tempo wurde signifikant von der Sprache beeinflusst (F[1,29] = 28.2, p < 0.001). # oder eben auf zweierlei Weise als gepaarter t-test: fremddiff = fremd %>% group_by(Vpn) %>% summarise(tempo = diff(tempo)) #dann auch die bessere Abbildung: boxplot(fremddiff$tempo) abline(h = 0) # oder komplizierter und trotzdem unschoener ggplot(fremddiff) + aes(y = tempo, x = factor(0)) + geom_boxplot() + geom_hline(yintercept = 0) shapiro.test(fremddiff$tempo) # fremddiff$tempo ist normalverteilt, und daher duerfen wir einen t-Test anwenden: t.test(fremddiff$tempo) # tempo wurde signifikant von der Sprache beeinflusst (t[29] = 5.3, p < 0.001). # oder 3. auf diese Weise (mit dem urspruenglichen data.frame): t.test(tempo ~ Sprache, paired = TRUE, data = fremd) # Es reicht natuerlich vollkommen EINER dieser drei moeglichen Loesungswege aus... #(und die Hilfslinie mit abline() duerfen Sie auch getrost weglassen) # ABER: auch wenn Sie die ezANOVA() benutzen, sollten Sie bitte die Differenzen berechnen # und davon die Abbildung erstellen sowie die Normalverteilung der Differenzen testen # (mit shapiro.test(fremddiff$tempo)) # 4. Diese Daten: dim(vcv) # zeigen Reaktionszeiten (RT) von verschiedenen Versuchspersonen (Subject) # auf 3 Vokale (Faktor Vowel) in verschiedenen davor und danach # kommenden Konsonanten-Kontexten (Context). Pruefen Sie durch eine Abbildung und # statistischen Test, ob die Reaktionszeiten vom Vokal (Vowel) und/oder # Sprachgruppe (Faktor Lang) beeinflusst wurden. # RT ist numerisch; Vowel kategorial; mehrere Random Factors (Subject und Context); # relativ viele Beobachtungen (dim(vcv) = 810;7), daher MM mit lmer() moeglich: ggplot(vcv) + aes(y = RT, x = Vowel) + geom_boxplot() + facet_wrap( ~ Lang) #oder: ggplot(vcv) + aes(y = RT, x = Vowel, col = Lang) + geom_boxplot() #oder: ggplot(vcv) + aes(y = RT, x = Vowel, linetype = Lang) + geom_boxplot() # oder beliebige Kombinationen davon with(vcv, table(Subject, interaction(Vowel, Lang))) with(vcv, table(Context, interaction(Vowel, Lang)))#oder einfacher aufgeteilt: with(vcv, table(Subject, Vowel))#-->Vowel within in Bezug zu Subject with(vcv, table(Subject, Lang))#--> Lang between in Bezug zu Subject with(vcv, table(Context, Lang))#-->Lang within in Bezug zu Context with(vcv, table(Context, Vowel))#-->Vowel within in Bezug zu Context #Formel: # RT~Vowel+Lang+Vowel:Lang-->RT~Vowel*Lang + (Vowel|Subject) + (Lang+Vowel|Context) vcv.lmer = lmer(RT ~ Vowel * Lang + (Vowel | Subject) + (Vowel + Lang | Context), data = vcv) vcv.step = step(vcv.lmer) anova(lmer(RT ~ Vowel + (Vowel | Subject) + (Vowel | Context), data = vcv)) get_model(vcv.step) anova(get_model(vcv.step)) # Die Daten zeigen einen signifikanten Einfluss vom Vokal (F[2,18.6] = 5.75, p < 0.05) # aber keine Einfluesse von der Sprachgruppe und es gab auch keine Interaktion zwischen diesen Faktoren. # da Vowel dreistufig ist, kann man post-hoc testen, # inwiefern sich a von i, i von u, und a von u unterscheiden pairs(emmeans(get_model(vcv.step), ~ Vowel)) #auch phsel(pairs(emmeans(get_model(vcv.step), ~ Vowel))) #Ein tukey-korrigierter post-hoc Test zeigte, dass fuer die Unterschiede nur die Paarungen # a vs. i (p<0.05) und a vs. u (p<0.05), aber nicht die Paarung i vs. u verantwortlich sind. # 5. Fuer diese Daten: dim(form10) # pruefen Sie durch eine Abbildung und statistischen Test, # ob F1 aus der Dauer vorhergesagt werden kann. # Was waere der vorhergesagte F1-Wert fuer eine Dauer von 50 ms? # Ist diese Vorhersage sehr verlaesslich? plot(F1 ~ Dauer, data = form10) form10.lm = lm(F1 ~ Dauer, data = form10) abline(form10.lm) #oder sehr umstaendlich ggplot(form10) + aes(y = F1, x = Dauer) + geom_point() + geom_abline(intercept = coef(form10.lm)[1], slope = coef(form10.lm)[2]) summary(form10.lm) predict(form10.lm, data.frame(Dauer = 50)) # Es gibt eine signifikante lineare Beziehung zwischen # F1 und Dauer (R^2 = 0.05, F[1,198] = 12.0, p < 0.001) #Das Modell sagt einen F1-Wert von 791.7138 fuer einen #Dauerwert von 50 vorher; wegen des geringen R^2-Wertes macht diese #Vorhersage aber so gut wie keinen Sinn. shapiro.test(resid(form10.lm)) # OK plot(resid(form10.lm)) # OK acf(resid(form10.lm)) # OK (auch wenn bei Lag=2 der blaue Bereich minimal ueberschritten wird). # 6. Fuer die Daten in diesem Data-Frame: dim(svot) # pruefen Sie mit einer Abbildung und einem # statistischen Text, inwiefern VOT von # Gender und Region beeinflusst wird. # (Vpn enthaelt das Sprecherkuerzel). # Immer, wenn wir einen Faktor mit mehr als 2 Stufen haben # oder/und mehr als einen Faktor, # koennen wir keinen t-Test anwenden. Wir benoetigen # eine Varianzanalyse oder ein Mixed Model. # Auch ggplot ist hier # dem Standard "boxplot" ueberlegen #Die geforderte Abbildung ggplot(svot)+ aes(y=VOT,x=Region,col=G)+ geom_boxplot() dim(svot) #[1] 24 4 # Nur 24 Beobachtungen sind sicher zu wenig für ein MM.! Daher ANOVA: #Der Test with(svot,table(Vpn,G)) with(svot,table(Vpn,Region)) # sowohl "G" als auch "Region" sind also (selbstverständlich) "between" # in Relation zum einzigen Random Factor "Vpn" ezANOVA(svot,.(VOT),.(Vpn),between=.(G,Region)) #$ANOVA #Effect DFn DFd F p p<.05 ges #1 G 1 18 9.810811 5.758440e-03 * 0.3527697 #2 Region 2 18 48.729730 5.439849e-08 * 0.8441011 #3 G:Region 2 18 3.972973 3.722434e-02 * 0.3062500 #wegen Interaktion: post-hoc-t-Tests (bonferroni-korrigiert) p=phoc(svot,.(VOT),.(Vpn),.(G,Region)) phsel(p$res) phsel(p$res,2) # EineVarianzanalyse zeigte jeweils einen Haupteffekt # fuer Geschlecht (F[1,18]=9.8,p<0.001) und Region (F[2,18]=48.7,p<0.001), # sowie eine Interaktion # zwischen beiden (F[2,18]=4.0,p<0.05). # Post-hoc-t-Tests # zeigten nur Unterschiede zwischen Maennern # aus den Regionen A vs. C (p<0.05),und zwischen Frauen aus A vs. B (p<0.01) # und aus A vs. C (p < 0.01) . #alternativ haette man schreiben koennen: # VOT wurde signifikant von der Region (F[2, 18] = 48.7, p < 0.001) und von # Gender (F[1, 18] = 9.8, p < 0.001) beeinflusst, und es gab eine signifikante # (F[2, 18] = 3.9, p < 0.05) Interaktion zwischen diesen Faktoren. # Post-hoc Tests zeigten keinen Einfluss von Gender auf VOT. Sie zeigten auch # signifikante Unterschiede zwischen Regionen A und C fuer Maenner (p < 0.05) und # fuer Frauen (p < 0.01) sowie zwischen Regionen A und B fuer Frauen (p < 0.01). # 7. Fuer diese Daten: dim(stefan) # erzeugen Sie eine Abbildung, # mit der (F1) als Funktion der Dauer (Dur), getrennt fuer die zwei Sprechtempi (Tempo) # gezeigt wird (also eine Abbildung, die zwei Bilder enthaelt, # eins fuer 'slow' und eins fuer 'fast'). # Schaetzen Sie anhand der Abbildung ein, ob F1 von der Dauer vorhergesagt werden kann. # Schaetzen Sie ebenfalls ein, ob eine Interaktion zwischen der Dauer und Sprechtempo vorliegt. a = ggplot(stefan) + aes(y = F1, x = Dur) + geom_point() + facet_wrap(~Tempo) a # "Offenbar steigt F1 mit zunehmender Dauer; die Dauer-Werte scheinen # unabhängig von der Geschwindigkeit zu sein; aber es scheint eine Interaktion # zu geben, d.h. bei "fast" scheint der Anstieg steiler zu sein." # Berechnen Sie mit einem linearen gemischten Modell, ob der Zusammenhang zwischen # der Dauer und F1 auch von der Sprechgeschwindigkeit abhängt. Bedenken Sie hierbei, # dass die Daten von verschiedenen Sprechern stammen (Vpn), und dass diese Sprecher # unterschiedliche Vokale (V: /e/ vs. /o/) in verschiedenen Wörtern gesprochen hatten. # Sie müssen nur die obenstehende Frage beantworten: "Haengt der Zusammenhang zwischen # der Dauer und F1 auch von der Sprechgeschwindigkeit ab?" Begruenden Sie ihre Antwort # in Worten! Welches Ergebnis muessen Sie anschauen? with(stefan,table(Vpn,Tempo)) with(stefan,table(Word,Tempo)) with(stefan,table(V,Tempo)) # jeweils "within" in Bezug zu Vpn, Word, V # Vorsicht: Dur ist numerisch, Tables machen hier keinen Sinn, da jeder Wert # einzigartig ist with(stefan,table(Vpn,Dur)) # d.h. numerische unabhängige Variablen sind auf jeden Fall wie "between" zu behandeln. stefan.lmer = lmer(F1~Dur*Tempo + (Tempo|Vpn) + (Tempo|Word) + (Tempo|V),data=stefan) stefan.step = step(stefan.lmer) anova(get_model(stefan.step)) # Das Modell ergibt nur einen Effekt für Dauer (F[1,105]=4.5, p < 0.05), aber nicht für # Sprechgeschwindigkeit; v.a. aber gibt es keine Interaktion zwischen # der Dauer und der Sprechgeschwindigkeit. Diese Interaktion wäre der entscheidende # Faktor, den man betrachten muss, um die Frage zu beantworten, # ob der Zusammenhang zwischen Dauer und F1 # bei unterschiedlichen Sprechgeschwindigkeiten unterschiedlich ist oder nicht. # Ein Extra (unnötig für die Beantwortung der Frage; d.h. das wird NICHT verlangt!!!!) # also nur zur Info und zum Vergleich zu Aufgabe 1 (in diesem Dokument): fixef(get_model(stefan.step)) ggplot(stefan) + aes(y = F1, x = Dur) + geom_point() +geom_abline(intercept=462.1840,slope=486.6206) # oder plot(F1~Dur, data=stefan);abline(fixef(get_model(stefan.step))) # != plot(F1~Dur, data=stefan) stefan.lm = lm(F1~Dur, data=stefan) abline(stefan.lm) coef(stefan.lm)!=fixef(get_model(stefan.step)) # weil bei lm() die unterschiedlichen Sprecher (von denen # jeweils mehrere Punkte stammen, die sich nach verschiedenen Vokalen und/oder Wörter aufteilen lassen) # NICHT berücksichtigt werden (daher brechen wir hier ab und ersparen uns das Testen der Voraussetzungen # für ein lm() ); # Vgl. hierzu Aufgabe 1! Dort wurden Vpn erwähnt, aber es gab genausoviele Vpn wie Beobachtungen: # es gab also nur eine Beobachtung pro Vpn; somit gibt es keine within-Vpn-Variation; # außerdem gab es keinen Grund, die Sprecher in Gruppen einzuteilen (also einen between-Vpn-Faktor anzunehmen) # (zumindest war das nicht erwähnt; es gab jedoch weitere Spalten, die eine Gruppeneinteilung erlaubt hätten, # ABER WIR BETRACHTEN BITTE NUR EINFLUSSFAKTOREN, DIE IN DER AUFGABENSTELLUNG ERWÄHNT WERDEN) # Ein Mixed Model mit fixed factor "tv.hours" und random factor "Vpn" würde bei Aufgabe 1 # auch mit folgender Fehlermeldung scheitern: tv.lmer = lmer(cdi~tv.hours + (1|Vpn),data = tv) # Fehler: number of levels of each grouping factor must be < number of observations # number of observations: nrow(tv) #80 # number of levels of each grouping factor length(levels(tv$Vpn)) # auch 80 # 8. Und noch ein kurzer Rausschmeisser! # Die voice-onset-time-Werte für /t/, # gesprochen von 16 Sprechern mit L1-Franzoesisch, # waren wie folgt: vot = c(-4,28,23,1,-13,20,21,-15,4,29,16,22,35,3,3,-4) # (a) Wie wahrscheinlich es ist, dass /t/ # voellig unaspiriert ist # ("voellig unaspiriert" = mit einem VOT-Wert von 0 ms)? # Berichten Sie hierbei auch das 95%-Konfidenzintervall # für die Daten in vot! # Folgend den Angaben in http://www.phonetik.uni-muenchen.de/%7Ejmh/lehre/sem/ss18/statistikR_files/t-Test.html # könnte man diese Frage nach dem 95%-Konfidenzintervall auf die folgende, # äußerst komplizierte Weise lösen: verwende qt() anstelle von qnorm(), # mit den Parametern # mu (= mean(vot) ), # SE (= sd(vot)/sqrt(length(vot)) ), # df (= Anzahl der Beobachtungen minus 1, # also (length(vot)-1) ), # und den Wahrscheinlichkeitswerten 2.5% (0.025) und 97.5% (0.975): a = mean(vot) + sd(vot)/sqrt(length(vot)) * qt(0.025, df=length(vot)-1) b = mean(vot) + sd(vot)/sqrt(length(vot)) * qt(0.975, df=length(vot)-1) a b # Puh, kompliziert. Viel, viel einfacher ist: #Wie wahrscheinlich es ist, dass /t/ # voellig unaspiriert ist? lässt sich übersetzen in # "die Wahrscheinlichkeit, dass die Werte in vot gleich 0 sind" # --> klassischer one-sample t.test mit mu=0, also t.test(vot) # One Sample t-test # # data: vot # t = 2.7137, df = 15, p-value = 0.01601 # alternative hypothesis: true mean is not equal to 0 # 95 percent confidence interval: # 2.26633 18.85867 # sample estimates: # mean of x # 10.5625 # t[15]=2.7, p < 0.05 (diesmal ausnahmsweise ganz genau: p = 0.01601); # Die Wahrscheinlichkeit, dass /t/ völlig unaspiriert ist, liegt bei 1.6 Prozent. # Das 95 %-Konfidenzintervall liegt zwischen 2.26633 und 18.85867 # (diese Werte werden beim t-Test-Ergebnis angegeben!!!); vgl. die händisch errechneten a b # # (b) Die VOT-Werte derselben 16 franzoesischen Sprecher # nach einem 6 monatigen Aufenthalt in Deutschland waren wie folgt: vot2 = c(2,28,26,4,-9,23,19,-6,11,33,24,31,41,10,7,8) # # Pruefen Sie anhand einer Abbildung und einem statistischen Test, # ob und falls ja, wie die voice-onset-time von dem Aufenthalt # in Deutschland beeinflusst wurde! boxplot(vot2-vot) t.test(vot2-vot) # Ein gepaarter t-Test ergab, dass nach einem sechsmonatigem Aufenthalt in Deutschland # französische Muttersprachler eine höhere VOT hatten (t[15] = 5.9, p < 0.001).