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)

soa = read.table(paste(pfad, "soa.txt", sep="")) vot = read.table(paste(pfad, "vot.txt", sep="")) library(multcomp)

alc = read.table(paste(pfad, "alcdata.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).

Frank & Jaeger (April, 2009) Post hoc comparisons Additional Issues: Random eﬀects diagnostics, multiple comparisons

http://hlplab.wordpress.com/2009/05/03/multilevel-model-tutorial/

Levy & Jaeger (2009) A Brief and Friendly Introduction to Mixed-Eﬀects Models in Psycholinguistics

2 Präsentationen hier vorhanden

Erste Veröﬀentlichung: Pinheiro & Bates (2000).

http://www.amazon.com/Mixed-Eﬀects-Models-S-S-Plus/dp/0387989579

(3)

soa: eine modifizierte Datei von Baayen et al (2008) (selbe Werte, andere Namen)

Subject und Items als Random Factors Subject und Items als Random Factors

F1 (F1-Werte im Vokal); Vpn (s1, s2, s3), W (Bart, Pfad, Start), Pfinal (medial, final)

3 Vpn ( produzierten 3 Wörter (Bart, Pfad, Start) jeweils 2 Mal:

einmal phrasenmedial, einmal phrasenfinal. Unterscheiden sich phrasenmedial und –final in F1?

Faktor W: between/within:

Faktor Pfinal: between/within:

Zuerst eine Abbildung

within within

(4)

Signifikanter Unterschied in Pfinal (zwischen long und short)?

Wort * Pfinal Interaktion?

Zuerst eine Abbildung attach(soa)

boxplot(F1 ~ Pfinal * W)

(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)

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)

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 eﬀect fallacy'.

(8)

Subject und Items als Random Factors 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)

Fixed: W, Pfinal (beide within) Random: Subject

bedeutet: signifikante Ergebnisse sind nicht nur für diese Vpn sondern allgemein für ähnliche Sprecher gültig.

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)

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'

Keine solchen Komplikationen in einem Mixed-Model Das Verfahren ist schrecklich (siehe Johnson, 2008)

(10)

Weitere Vorteile von einem Mixed-Model Weitere Vorteile von einem Mixed-Model

Vpn

i e a

lang. schnell Sprechtempo

Vokal

Sprache engl. oder span.

i e a w_1.1w₂ w₃ w₄ w₅ w₆ between

within

w_1.2 w_1.3

w_1.10 ...

Mittelwert

Es muss nicht gemittelt werden pro Stufen-Kombinationen Die Zellen müssen nicht vollständig pro Vpn sein.

(11)

Nachteile von einem Mixed-Model 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:

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.

http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests

(12)

lmer() lmer()

lmer() ist die Funktion für mixed modelling.

3. by-subject intercept adjustment

4. by-item intercept adjustment 1. Neigung

2. Intercept

y = bx + k + e

^{^} _Vpn

+ e

_W

Eingeschätzte Werte

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

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)

F1.lmer = lmer(F1 ~ Pfinal + (1 | Vpn) + (1 | W)) abhängige Variable

unabhängige Variable(n)

Vpn und W als Random Factors lmer()

lmer()

1. by-subject and by-item intercept adjustment

= Sprecher- und Wortvariationen werden für Pfinal (ohne zwischen den Stufen zu diﬀerenzieren) ausgeklammert.

F1.lmer = lmer(F1 ~ Pfinal + (1 + Pfinal | Vpn) + (1 +Pfinal | W)) 2. by-subject and by-item intercept and slope adjustment

= Sprecher- und Wortvariationen werden getrennt für die

Stufen (long, short) von Pfinal berechnet und ausgeklammert.

(14)

lmer() lmer()

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

print(F1.lmer, corr=F)

print(F1.lmer, corr=F) Fixed eﬀects (bezieht sich daher auf Pfinal)

ranef(F1.lmer)

ranef(F1.lmer) Random eﬀects (bezieht sich

daher auf Vpn und W)

$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.

auch BLUPS genannt (best

linear unbiased predictor). Ein BLUP pro Vpn und pro Wort

y = bx + k + e

^{^} _Vpn

+ e

_W

fitted(F1.lmer) fitted(F1.lmer)

[1] 473.6037 515.6890 515.3677 454.7148...

Eingeschätzte Werte

(15)

lmer() lmer()

$Vpn

(Intercept) s1 -20.557646

$W

(Intercept) Bart -27.94979

fitted(F1.lmer)[4]

[1] 454.7148

Fixed Random

y = bx + k + e

^{^} _Vpn

+ e

_W

Eingeschätzer Wert für phrasenmedialer (Stufe, short) Bart, Vpn. s1

contrast(Pfinal) short

long 0 short 1

-18.889 * 1 + 522.111 -20.557646 -27.94979 1 [1] 454.7146

(16)

lmer(), t-Werte und Basis-Kodierung lmer(), t-Werte und Basis-Kodierung

Die Stufe short vom Faktor Pfinal ist 3.419 Standard- Abweichung von der Basis-Stufe desselben Faktors entfernt.

Basis-Stufe levels(Pfinal)

[1] "long" "short"

Die Basis-Stufe kann man mit relevel() auswählen.

(Also: der absolute Abstand zwischen long und short ist 3.419 Standard-Abweichungen).

(17)

Die Basis-Stufe kann man mit relevel() auswählen.

Pfinal2 = relevel(Pfinal, "short")

F1b.lmer = lmer(F1 ~ Pfinal2 + (1 | Vpn) + (1 | W))

Fixed effects:

Estimate Std. Error t value (Intercept) 503.222 19.604 25.670 Pfinal2long 18.889 5.525 3.419

print(F1b.lmer, corr=F)

Kein

Freiheitsgrad, keine p-Werte

(18)

lmer(), t-Werte, p-Werte lmer(), t-Werte, p-Werte

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

2 * ( 1 - pt(3.419, 16)) [1] 0.003516309

Aber der p-Wert ist anti-konservativ (ggf. zu niedrig) und daher nicht zuverlässig.

(19)

lmer() und MCMC sampling 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)

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 lmer() und MCMC sampling

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

F1.fnc = pvals.fnc(F1.lmer) F1.fnc$fixed

Highest Posterior Density

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).