library(languageR) pfad = "Verzeichnis wo diese Objekte unten gespeichert wurden" lex = read.table(paste(pfad, "lex.txt", sep="/")) anna = read.table(paste(pfad, "anna.txt", sep="/")) yuki = read.table(paste(pfad, "yuki.txt", sep="/")) alc = read.table(paste(pfad, "alcdata.txt", sep="/")) ## TEIL A: Abhaengige Variable ist kontinuierlich ######## Daten von Baayen # visual lexical decision latencies # Haben Sprache und Wortlaenge einen signifikanten Einfluss auf # die Reaktionszeiten? head(lex) lex.lmer = lmer(rt ~ Length * Lang + (1 | Subj) + (1 | Word), data=lex, REML=F) lex.ranef = ranef(lex.lmer, postVar=T) dotplot(lex.ranef) anova(lex.lmer) # Die Freiheistgrade: Anzahl der Stichproben - Anzahl der df im Nenner nrow(lex) - 3 # Wahrscheinlichkeiten fuer Lang: 1 - pf(6.2976, 1, 1553) # NB fast das gleiche Ergebnis mit 70 dfs 1 - pf(6.2976, 1, 70) # Trial koennte eine Rolle spielen. d.h. es gibt vielleicht # eine Vpn- und Word-Variation im Trial. # Trial ist jedoch einen Fixed- und keinen Random-effect. xylowess.fnc(rt ~ Trial | Subj, data=lex) xylowess.fnc(rt ~ Trial | Word, data=lex) # Trial als Fixed-Faktor und Herausrechnnen des Einflusses # von Trial auf Subj lex.lmerb = lmer(rt ~ Length * Lang * Trial + (1 + Trial | Subj) + (1 | Word), data=lex, REML=F) # (Lang, Length) anova(lex.lmerb) # Vielleicht gibt es auch eine Wort-bedingte Variation in Trial lex.lmerc = lmer(rt ~ Length * Lang * Trial + (1 + Trial | Subj) + (1 + Trial | Word), data=lex, REML=F) lex.lmerc = lmer(rt ~ Length * Lang * Trial + (1 + Trial | Subj) + (1 + Trial | Word), data=lex, REML=F) # Unterscheiden sich lex.lmerb und lex.lmerc? anova(lex.lmerb, lex.lmerc) # Nein. Wir bleiben dann bei dem einfacheren Modell lex.lmerb, der sich # von lex.lmera unterscheidet anova(lex.lmer, lex.lmerb) # Length ist signifikant, Lang ist signifikant und Length:Lang Interaktion anova(lex.lmerb) # Graphische Darstellung von Length, Lang, und Interaktion xyplot(rt ~ Length | Lang, lex, groups=Subj) xyplot(rt ~ Length | Lang, lex, groups=Word) bwplot(rt ~ Lang | Length, data=lex) densityplot(~rt | Length, lex, groups=Lang, auto.key=T) densityplot(~rt | Length, lex, layout=c(1, 8), groups=Lang) # Unterscheidet sich die Gruppe English in Length? temp = with(lex, Lang == "English") lexengl = lex[temp,] lmer.engl = lmer(rt ~ Length * Trial + (1 + Trial | Subj) + (1 + Trial | Word), data=lexengl, REML=F) # Length ist signifikant anova(lmer.engl) 1 - pf(7.5958, 1, nrow(lexengl)-3) # Unterscheidet sich die Gruppe Other in Length? lexother = lex[!temp,] lmer.other = lmer(rt ~ Length * Trial + (1 + Trial | Subj) + (1 + Trial | Word), data=lexother, REML=F) anova(lmer.other) ########### TEIL B: Abhaengige Variable ist diskret (Logistisch) ### Daten von Anna head(anna) anna.lmer1 = lmer(cbind(Incorrect, Correct) ~ Cond * Vow * Cons + (1 | Sp) + (1 | W), data=anna, family=binomial) print(anna.lmer1, corr=F) anna.lmer2 = lmer(cbind(Incorrect, Correct) ~ Cond + Vow + Cons + Cond:Vow + (1 | Sp) + (1 | W), data=anna, family=binomial) # Keine signifikanten Unterschiede. Wir waehlen daher # das einfachere Modell anna.lmer2 anova(anna.lmer1, anna.lmer2) # Post-hoc tests durchfuehren # Wir wollen wissen: unterscheiden sich "p" und "q" # im Faktor Cond im Context "i" und "a" vom Faktor Vow # Hier hat Vow die Basis "a" und es gibt einen # signifikanten Unterschied zwischen "p" und "q" vom Faktor Cond print(anna.lmer2, corr=F) Condq 1.69203 0.29771 5.683 1.32e-08 *** # Jetzt fuehren wir den Test nochmal durch, diesmal mit "i" # vom Faktor Vow als Basis Vow2 = with(anna, relevel(Vow, "i")) anna.lmer3 = lmer(cbind(Incorrect, Correct) ~ Cond + Vow2 + Cons + Cond:Vow2 + (1 | Sp) + (1 | W), data=anna, family=binomial) # Und hier gibt es wieder einen signifikanten Unterschied zwischen "p" und "q" Condq 0.77783 0.26259 2.962 0.00306 ** # Daher Schlussfolgerung: Trotz Interaktion gibt es # einen signifikanten Unterschied zwischen "p" und "q" # sowohl im Kontext "i" als auch im Kontext "a" #################################################################### ### Methode von JMH. Wir binden Cons und Vow zusammen. #################################################################### # Cond und Vow zusammenlegen Condvow = factor(paste (with(anna, Cond), with(anna, Vow), sep=".")) anna.lmerp = lmer(cbind(Incorrect, Correct) ~ Condvow + Cons + (1 | Sp) + (1 | W), data=anna, family=binomial) print(anna.lmerp, corr=F) # Wir wollen wissen: unterscheiden sich p.i von q.i und p.a von q.a # (d.h. gibt es signifikante Unterschiede in der Klassifikation # zwischen Prae- (p) und Post- (q) aspiration im [i] und im [a] Kontext? # Was ist die Basis-Stufe? levels(Condvow) [1] "p.a" "p.i" "q.a" "q.i" # Daher zeigt uns diese Zeile: Condvowq.a 1.6920 0.2977 5.683 1.32e-08 *** # dass sich p.a von q.a signifikant unterscheiden. # UND UEBRIGENS IST DIESE ZEILE GENAU DIESELBE WIE DIEJENIGE # OBEN MIT Vow = "a" ALS BASIS Condq 1.69203 0.29771 5.683 1.32e-08 *** # Wir koennten denselben Test nochmal durchfuehren, diesmal mit p.i als Basis. # Dann muessten wir den p.i vs. q.i Unterschied auf dieselbe Weise # ablesen koennen Condvow = relevel(Condvow, "p.i") anna.lmerp2 = lmer(cbind(Incorrect, Correct) ~ Condvow + Cons + (1 | Sp) + (1 | W), data=anna, family=binomial) print(anna.lmerp2, corr=F) # Condvowq.i 0.77784 0.26259 2.962 0.00305 ** # UND UEBRIGENS IST DIESE ZEILE GENAU DIESELBE WIE DIEJENIGE # OBEN MIT Vow = "i" ALS BASIS Condq 0.77783 0.26259 2.962 0.00306 ** ## DAHER: ES GIBT SIGNIFIKANTE UNTERSCHIEDE ZWISCHEN ## "p" und "q" SOWOHL IM "a" ALS AUCH IM "i" KONTEXT ######################################################################################################################################## # Abbildung. p = anna[,2]/apply(anna[,1:2], 1, sum) # Arc-Sinus Transformation p = (2/pi) * asin(sqrt(p)) anna = cbind(p, anna) bwplot(p ~ Cond | Vow, data=anna, auto.key=T) densityplot(~p | Vow, anna, groups=Cond, auto.key=T) ######### Daten von Yuki head(yuki) yuki.lmer = lmer(cbind(P, Q) ~ Stim + Mode + L1 + H1 + Stim:Mode + Stim:L1 + Stim:H1 + (1|VP), data=yuki, family=binomial) print(yuki.lmer, corr=F) Stim:ModeL -0.06417 0.03006 -2.13 0.032782 * Stim:L1L 0.04668 0.03054 1.53 0.126377 Stim:H1L 0.08551 0.03018 2.83 0.004608 ** ### 1. Inwiefern sind die Interaktionen mit Stim signifikant? E vs L vom Faktor Mode unterscheiden sich in der Verteilung der Urteile auf dem Stimulus-Kontinuum (z = 2.13, p < 0.05) H vs L vom Faktor L1 unterscheiden sich nicht in der Verteilung der Urteile auf dem Stimulus-Kontinuum H vs L vom Faktor H1 unterscheiden sich (p < 0.01) in der Verteilung der Urteile auf dem Stimulus-Kontinuum # 2. Jetzt berechnen wir die Umkipppunkte und wenden # einen gepaarten t-test an, um zu pruefen, ob sich # die Stufen der Faktoren diesbezueglich unterscheiden. # In diesem Beispiel machen wir das nur fuer E vs L vom # Faktor Mode - aber das Prinzip ist genau das gleiche fuer # die anderen Faktoren # *** Siehe bitte einige Zeilen unten fuer weitere Details temp = with(yuki, Mode)=="E" ye = yuki[temp,] yl = yuki[!temp,] ye.lmer = lmer(cbind(P, Q) ~ Stim + (1|VP), data=ye, family=binomial) yl.lmer = lmer(cbind(P, Q) ~ Stim + (1|VP), data=yl, family=binomial) # Hier sind die Koeffiziente der logistischen Regression, getrennt pro Vpn. fuer "E" coef.e = coef(ye.lmer)$VP # Hier ist der Umkipppunkt getrennt pro Vpn. fuer "E" umkip.e = -coef.e[,1]/coef.e[,2] # Wie oben fuer "L" coef.l = coef(yl.lmer)$VP umkip.l = -coef.l[,1]/coef.l[,2] # Wir vegleichen hier die Umkipppunkten mit einem gepaarten t-test t.test(umkip.e, umkip.l, paired=T) # Hier ein Boxplot derselben Sache umkip = c(umkip.e, umkip.l) umkip.labs = c(rep("E", length(umkip.e)), rep("L", length(umkip.l))) boxplot(umkip ~ umkip.labs) # Oder bwplot(umkip ~ umkip.labs) ########## ***Weitere Details nur zur Information! # coef.e variiert pro Vpn nur im Intercept # nicht in der Neigung. (coef.e[,1] variiert, coef.e[,2] # ist gleich pro Vpn). Dies ist weil wir (1|VP) in # der Formel haben = eine by-Subject-Anpassung # fuer das Intercept. Wenn wir eine by-Subject-Anpassung # auch fuer die Neigung # wollen, dann (1+Stim|VP), wie unten. # Die by-Subject-Anpassung in diesem Fall # ist in unlist(ranef(ye.lmer)). Der Intercept # selbst (ueber alle Vpn berechnet) ist # wie immer unter Intercept in print(ye.lmer, corr=F) # zu finden, also 7.66439. Daher ist die # by-Subject-Anpassung 7.66439 - unlist(ranef(ye.lmer)) # Genau diese Werte sind in der ersten Spalte von coef.e # wie round(coef.e[,1] - (7.66439 + unlist(ranef(ye.lmer))), 4) # bestaetigt # 3. Hier wollen wir jetzt versuchen, eine durch die logistische Regression # geglaettete Kurve der Urteile als Funktion vom StimulusNr. zu erzeugen # Hier einmal fuer "E" und dann getrennt fuer "L" # Response-Kurve, Vpn 1:10 fuer "E" und "L" for(j in 1:10){ curve(exp(coef.e[j,2]*x + coef.e[j,1])/(1+ exp(coef.e[j,2]*x+coef.e[j,1])), xlim=c(1, 10), ylab="", xlab="", ylim=c(0, 1)) curve(exp(coef.l[j,2]*x + coef.l[j,1])/(1+ exp(coef.l[j,2]*x+coef.l[j,1])), xlim=c(1, 10), ylab="", xlab="", add=T, col=2, ylim=c(0, 1)) abline(v=umkip.e[j], lty=1) abline(v=umkip.l[j], col=2, lty=1) locator(1) } # Response-Kurven aller Vpn auf einander ueberlagern ohne Umkipppunkt # Man wird hier den Unterschied zwischen E und L nicht unbedingt # deutlich wegen der hohen Hoerervariabilitaet sehen koennen # (trotzdem ist der Unterschied zwischen E und L hoch signifikant # wie wir oben sahen, *weil rot und schwarz pro Vpn gepaart sind* add = F for(j in 1:nrow(coef.e)){ if(j !=1) add = T curve(exp(coef.e[j,2]*x + coef.e[j,1])/(1+ exp(coef.e[j,2]*x+coef.e[j,1])), xlim=c(1, 11), ylab="", xlab="", ylim=c(0, 1), add=add) curve(exp(coef.l[j,2]*x + coef.l[j,1])/(1+ exp(coef.l[j,2]*x+coef.l[j,1])), xlim=c(1, 11), ylab="", xlab="", add=T, col=2, ylim=c(0, 1)) } # Fuer den Mittelwert nimmt man einfach den Intercept und die Neigung # aus ye.lmer und yl.lmer, daher: coef.m.e = c(7.66439, -0.98843) coef.m.l = c(7.7271, -1.0436) # Umkipp = -k/m um.m.e = -coef.m.e[1]/coef.m.e[2] um.m.l = -coef.m.l[1]/coef.m.l[2] curve(exp(coef.m.e[2]*x + coef.m.e[1])/(1+ exp(coef.m.e[2]*x+coef.m.e[1])), xlim=c(1, 11), ylab="", xlab="", ylim=c(0, 1)) curve(exp(coef.m.l[2]*x + coef.m.l[1])/(1+ exp(coef.m.l[2]*x+coef.m.l[1])), xlim=c(1, 11), ylab="", xlab="", ylim=c(0, 1), add=T, col=2) abline(v=um.m.e, lty=2) abline(v=um.m.l, lty=2, col=2) ## 3b. Wie waere es, wenn wir nicht nur ein by-Subjects-Intercept # sondern auch eine by-Subjects-Neigung berechnen wollen. Dann: ye.lmer2 = lmer(cbind(P, Q) ~ Stim + (1 + Stim | VP), data=ye, family=binomial) yl.lmer2 = lmer(cbind(P, Q) ~ Stim + (1 + Stim | VP), data=yl, family=binomial) coef.e2 = coef(ye.lmer2)$VP umkip.e2 = -coef.e2[,1]/coef.e2[,2] coef.l2 = coef(yl.lmer2)$VP umkip.l2 = -coef.l2[,1]/coef.l2[,2] t.test(umkip.e2, umkip.l2, paired=T) umkip2 = c(umkip.e2, umkip.l2) boxplot(umkip2 ~ umkip.labs) ## usw, der Rest genau wie oben. ### Daten von Sabine alc.lmer = lmer(as.matrix(alc[, 1:2]) ~ 0 + Mode + Alc * G + (1 | Vpn) + (1 | Mode:Vpn), data=alc, family=binomial) print(alc.lmer, corr=F) alc.lmer = lmer(as.matrix(alc[, 1:2]) ~ Mode + Alc * G + (1 + Mode | Vpn), data=alc, family=binomial) # Der na-a Alc Unterschied ist signifikant Alcna -0.20571 0.05334 -3.86 0.000115 *** # aber es gibt eine Interaktion mit Gender Alcna:GM 0.20225 0.07065 2.86 0.004198 ** # Eindeutig ist F die Basis. Daher trifft dieses signifikante # na-a Ergebnis nur fuer Frauen zu. # Nochmal den Test durchfuehren mit M als basis G2 = with(alc, relevel(G, "M")) alc.lmer2 = lmer(as.matrix(alc[, 1:2]) ~ 0 + Mode + Alc * G2 + (1 | Vpn) + (1 | Mode:Vpn), data=alc, family=binomial) print(alc.lmer2, corr=F) # Der na-a Unterschied ist fuer Maenner NS Alcna -0.003432 0.046365 -0.07 0.94100 ##### Bestaetigung durch Methode JMH #################################################################### ############################################ post-hoc tests both = factor(paste(with(alc, Alc), with(alc, G), sep=".")) alc.lmerp = lmer(as.matrix(alc[, 1:2]) ~ 0 + Mode + both + (1 | Vpn) + (1 | Mode:Vpn), data=alc, family=binomial) # a.F ist Basis levels(both) # wir wollen vergleichen: a.M mit na.M sowie a.F mit na.F # man sieht schon von print(alc.lmerp, corr=F), dass diese NS und S sind #################################################################### # Abbildung. p = alc[,1]/alc[,2] # Arc-Sinus Transformation p = (2/pi) * asin(sqrt(p)) alc = cbind(p, alc) temp = with(alc, G=="M") alcm = alc[temp,] alcf = alc[!temp,] bwplot(p ~ Alc | Mode, data=alcm) bwplot(p ~ Alc | Mode, data=alcf) densityplot(~p | Mode, alcm, groups=Alc, auto.key=T) densityplot(~p | Mode, alcf, groups=Alc, auto.key=T)