Mixed models

1. Mixed models Jonathan Harrington Die R-Befehle: mixed.txt pfad = "Das Verzeichnis, wo die Daten gespeichert sind" attach(paste(pfad, "anova1"), sep="") library(car) library(languageR) library(multcomp) alc = read.table(paste(pfad, "alcdata.txt", sep="")) soa = read.table(paste(pfad, "soa.txt", sep="")) vot = read.table(paste(pfad, "vot.txt", sep=""))

2. Mixed models Baayen, R.H. (in press) Analyzing Linguistic Data: A practical introduction to Statistics. Kapitel 7http://www.ualberta.ca/~baayen/publications/baayenCUPstats.pdf Artikel in einem Special Issue im Journal of Memory and Language, Vol. 59. insbesondere: Baayen, Davidson & Bates (2008); Quene & van den Bergh (2008); Jaeger (2008). 2 Präsentationen hier vorhanden http://hlplab.wordpress.com/2009/05/03/multilevel-model-tutorial/ Levy & Jaeger (2009) A Brief and Friendly Introduction to Mixed-Effects Models in Psycholinguistics Frank & Jaeger (April, 2009) Post hoc comparisons Additional Issues: Random effects diagnostics, multiple comparisons Erste Veröffentlichung: Pinheiro & Bates (2000). http://www.amazon.com/Mixed-Effects-Models-S-S-Plus/dp/0387989579

3. Subject und Items als Random Factors 3 Vpn ( produzierten 3 Wörter (Bart, Pfad, Start) jeweils 2 Mal: einmal phrasenmedial, einmal phrasenfinal. Unterscheiden sich phrasenmedial und –final in F1? soa: eine modifizierte Datei von Baayen et al (2008) (selbe Werte, andere Namen) F1 (F1-Werte im Vokal); Vpn (s1, s2, s3), W (Bart, Pfad, Start), Pfinal (medial, final) Faktor W: between/within: within Faktor Pfinal: between/within: within Zuerst eine Abbildung

4. Zuerst eine Abbildung attach(soa) boxplot(F1 ~ Pfinal * W) Signifikanter Unterschied in Pfinal (zwischen long und short)? Wort * Pfinal Interaktion?

5. soa.t = Anova.prepare(soa, c("s", "w", "w", "d")) soa.lm = lm(soa.t$d ~ 1) Anova(soa.lm, idata=soa.t$w, idesign = ~ W * Pfinal) Type III Repeated Measures MANOVA Tests: Pillai test statistic Df test stat approx F num Df den Df Pr(>F) (Intercept) 1 1.00 1509.78 1 2 0.0006617 *** W 1 0.98 26.67 2 1 0.1356518 Pfinal 1 0.78 7.24 1 2 0.1147994 W:Pfinal 1 0.62 0.81 2 1 0.6186755

6. Subject und Items als Random Factors Jedoch hat Pfinal eindeutig denselben Einfluss pro Wort. d.h. wir sollten hier die Variation zwischen den Wörtern ausklammern. (Dass der Mittelwert von Bart kleiner ist im Vgl. zu Pfad oder Start ist für diese Fragestellung uninteressant - genauso uninteressant wie die Variation zwischen Sprechern).

7. Subject und Items als Random Factors Wir wollen also 2 Sorten von Variation gleichzeitig ausklammern: wegen der Sprecher (Subject as a random factor) wegen der Wörter (Item as a random factor) Dadurch lösen wir gleichzeitig ein großes Problem in der Statistik (Clark, 1973): 'language as fixed effect fallacy'.

8. Subject und Items als Random Factors soa.t = Anova.prepare(soa, c("s", "w", "w", "d")) soa.lm = lm(soa.t$d ~ 1) Anova(soa.lm, idata=soa.t$w, idesign = ~ W * Pfinal) Random: Subject bedeutet: signifikante Ergebnisse sind nicht nur für diese Vpn sondern allgemein für ähnliche Sprecher gültig. Fixed: W, Pfinal (beide within) signifikante Ergebnisse in Pfinal sind nur für W gültig (Bart, Pfad, Start): d.h. damit die Ergebnisse allgemein für ähnliche Wörter gültig sind, müsste noch einen RM-Manova durchgeführt werden mit Wort als Random Faktor.

9. Subject und Items als Random Factors ein by-subject RM-Manova (Subject als Random Faktor) ein by-item RM-Manova (Wort als Random Faktor) Die Ergebnisse werden kombiniert in einer Statistik minF' Das Verfahren ist schrecklich (siehe Johnson, 2008) Keine solchen Komplikationen in einem Mixed-Model

10. Weitere Vorteile von einem Mixed-Model Es muss nicht gemittelt werden pro Stufen-Kombinationen Die Zellen müssen nicht vollständig pro Vpn sein. Sprache engl. oder span. between Vpn within Sprechtempo lang. schnell Vokal i e a i e a w1.1 w2 w3 w4 w5 w6 w1.2 w1.3 Mittelwert ... w1.10

11. Nachteile von einem Mixed-Model library(lme4) und mixed modelling (MM) überhaupt sind noch in der Entwicklungsphase. Daher bugs, häufige code Änderungen und einige Teile des Verfahrens sind in R noch nicht ganz vollständig. Siehe: http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests Mit MM können zwar Werte aus der t- und F-Verteilung berechnet werden, aber diese lassen sich nur schwierig und vielleicht sogar ungenau in Wahrscheinlichkeiten umsetzen – weil die Freiheitsgrade nicht eindeutig berechnet werden können.

12. ^ y = bx + k + eVpn + eW lmer() lmer() ist die Funktion für mixed modelling. berechnet ein lineares Modell, in dem der Abstand zwischen eingeschätzen und tatsächlichen Werten minimiert wird minimiert den Abstand durch REML (restricted maximum likelihood) muss mindestens einen Random-Faktor enthalten vereinheitlicht RM-Anova und logistische Regression 3. by-subject intercept adjustment 1. Neigung Eingeschätzte Werte 2. Intercept 4. by-item intercept adjustment x ist ein Faktor-Code* (zB 0 oder 1 für 2 Stufen) *siehe: http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm

13. lmer() abhängige Variable Vpn und W als Random Factors unabhängige Variable(n) 1. by-subject and by-item intercept adjustment F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W)) = Sprecher- und Wortvariationen werden für Pfinal (ohne zwischen den Stufen zu differenzieren) ausgeklammert. 2. by-subject and by-item intercept and slope adjustment F1.lmer = lmer(F1 ~ Pfinal + (1 + Pfinal | Vpn) + (1 +Pfinal | W)) = Sprecher- und Wortvariationen werden getrennt für die Stufen (long, short) von Pfinal berechnet und ausgeklammert.

14. lmer() ^ y = bx + k + eVpn + eW print(F1.lmer, corr=F) Fixed effects (bezieht sich daher auf Pfinal) Linear mixed model fit by REML Fixed effects: Estimate Std. Error t value (Intercept) 522.111 19.604 26.633 Pfinalshort -18.889 5.525 -3.419 Random effects (bezieht sich daher auf Vpn und W) ranef(F1.lmer) auch BLUPS genannt (best linear unbiased predictor). Ein BLUP pro Vpn und pro Wort $Vpn (Intercept) s1 -20.557646 s2 22.948070 s3 -2.390424 $W (Intercept) Bart -27.94979 Pfad 14.13553 Start 13.81427 summieren auf 0. fitted(F1.lmer) Eingeschätzte Werte [1] 473.6037 515.6890 515.3677 454.7148...

15. ^ y = bx + k + eVpn + eW lmer() Random Fixed $Vpn (Intercept) s1 -20.557646 Linear mixed model fit by REML Fixed effects: Estimate Std. Error t value (Intercept) 522.111 19.604 26.633 Pfinalshort -18.889 5.525 -3.419 $W (Intercept) Bart -27.94979 Eingeschätzer Wert für phrasenmedialer (Stufe, short) Bart, Vpn. s1 -18.889 * 1 + 522.111 -20.557646-27.949791 contrast(Pfinal) [1] 454.7146 short long 0 short 1 fitted(F1.lmer)[4] [1] 454.7148

16. lmer(), t-Werte und Basis-Kodierung Linear mixed model fit by REML Fixed effects: Estimate Std. Error t value (Intercept) 522.111 19.604 26.633 Pfinalshort-18.889 5.525 -3.419 Die Stufe short vom Faktor Pfinalist 3.419 Standard-Abweichung von der Basis-Stufe desselben Faktors entfernt. Basis-Stufe levels(Pfinal) [1] "long" "short" (Also: der absolute Abstand zwischen long und short ist 3.419 Standard-Abweichungen). Die Basis-Stufe kann man mit relevel() auswählen.

17. lmer(), t-Werte und Basis-Kodierung Linear mixed model fit by REML Fixed effects: Estimate Std. Error t value (Intercept) 522.111 19.604 26.633 Pfinalshort -18.889 5.525 -3.419 Die Basis-Stufe kann man mit relevel() auswählen. Pfinal2 = relevel(Pfinal, "short") Kein Freiheitsgrad, keine p-Werte F1b.lmer = lmer(F1 ~ Pfinal2 + (1 | Vpn) + (1 | W)) print(F1b.lmer, corr=F) Fixed effects: Estimate Std. Error t value (Intercept) 503.222 19.604 25.670 Pfinal2long 18.889 5.525 3.419

18. lmer(), t-Werte, p-Werte Linear mixed model fit by REML Fixed effects: Estimate Std. Error t value (Intercept) 522.111 19.604 26.633 Pfinalshort -18.889 5.525 -3.419 Die Wahrscheinlichkeit könnte mit der höchst möglichen Anzahl der Freiheitsgrade von der t-Verteilung berechnet werden. Freiheistgradanzahl = Anzahl der Werte - (Anzahl der Stufen) = 16 [1] 0.003516309 2 * ( 1 - pt(3.419, 16)) Aber der p-Wert ist anti-konservativ (ggf. zu niedrig) und daher nicht zuverlässig.

19. lmer() und MCMC sampling Um genauere Wahrscheinlichkeiten der im lmer() entstandenen t-Werte zu berechnen, gibt es Markov Chain Monte Carlo sampling (MCMC). kann zZt. für solche Modelle angewandt werden F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W)) nicht für diese F1.lmer = lmer(F1 ~ Pfinal + (1 + Pfinal | Vpn) + (1 +Pfinal | W))

20. lmer(), t-Werte und Basis-Kodierung Bitte beachten: es gibt einen Bug in pvals.fnc(), sodass die Werte von F1.lmer (!!) nach der Durchführung der Funktion geändert werden. Daher bitte immer gleich nach pvals.fnc() nochmal die lmer() Funktion duchführen! F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W)) F1.fnc = pvals.fnc(F1.lmer) F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W))

21. lmer() und MCMC sampling F1.fnc = pvals.fnc(F1.lmer) F1.fnc$fixed Highest Posterior Density Estimate MCMCmean HPD95lower HPD95upper pMCMC Pr(>|t|) (Intercept) 522.11 522.13 491.88 554.614 0.0001 0.0000 Pfinalshort -18.89 -18.88 -35.94 -1.013 0.0356 0.0035 2 * ( 1 - pt(3.419, 16)) Die Daten wurden mit einem Mixed Model mit Subject und Word als Random Faktoren und phrasenfinale Längung als unabhängige Variable analysiert. Alle Wahrscheinlichkeiten wurden mit Markov-Chain-Monte-Carlo sampling (MCMC) berechnet (Baayen et al, 2008). Die phrasenfinale Längung hatte einen signifikanten Einfluss auf F1 (t = 3.42, MCMC: p < 0.05).