load("RTL_beispieldaten_clean.RData")8 Statistische Analysen
8.1 Laden von Daten und Paketen
Alle statistischen Analysen verwenden den Datensatz, den wir in den letzten Kapiteln schon gescreent und bereinigt haben.
So wird der Datensatz nun wieder geladen:
Wir benötigen auch nochmal Daten zu den Versuchspersonen und laden deshalb auch nochmal den Datensatz mit den demographischen Daten.
load("RTL_beispieldaten.RData")Zusätzlich verwenden wir ein paar Packages, die für die statistischen Analysen verwendet werden, nämlich:
- afex benötigen wir zur Berechnung von (Generalized) Linear Mixed Models.
- emmeans verwenden wir zur Berechung von durch ein Modell geschätzten Werten (z.B. für Abbildungen) und Berechnung von post-hoch Tests
library(tidyverse)
library(afex)
library(emmeans)8.2 Vorwort zu statistische Analysen
Bevor wir an die Analysen der Antwort- und RT-Daten gehen, hier ein paar allgemeine Informationen vorweg.
8.2.1 Statistische Auswertung von Prozentdaten
Anteile oder Prozentwerte sind in der kognitiven Psychologie sehr üblich. Sie entstehen z.B. dann, wenn Versuchspersonen Antoworten geben, die korrekt oder inkorrekt sein können.
Prozentwerte darf man nicht mit einer ANOVA analysieren. Details dazu finden sich in: Jaeger, T. F. (2008). Categorical data analysis: Away from ANOVAs (transformation or not) and towards logit mixed models. Journal of Memory and Language, 59(4), 434–446. https://doi.org/10.1016/j.jml.2007.11.007
Stattdessen analysiert man sie mit generalized linear mixed models (GLMM).
Diese modellieren die Wahrscheinlichkeit, dass ein Trial korrekt oder inkorrekt beantwortet war. Dies passiert auf Ebene eines einzelnen Trials: Statt Prozentwerte zu berechnen, wie viele Trials richtig waren, beantworten GLMM die Frage, wie wahrscheinlich ein gegebener Trial korrekt oder inkorrekt beantwortet wird, gegeben die Faktorausprägungen des jeweiligen Trials.
Dabei ändert also ein Faktor die Wahrscheinlichkeit, dass man richtig oder falsch antwortet; dadurch ändert sich natürlich indirekt auch die Prozentzahl richtiger/falscher Trials.
GLMM vs. ANOVA: Macht’s einen Unterschied in den Ergebnissen?
Dass man für %-Werte GLMM nutzen muss, ist vielen Psycholog:innen offenbar immer noch nicht bekannt, oder es wird ignoriert. In der Praxis ergeben beide Verfahren, ANOVA und GLMM, oft die gleichen qualitativen Antworten, obwohl Jaeger genau dies in seinem Paper bestreitet.
In unserer Gruppe verwenden wir grundsätzlich die angemessene Analyseform, da man nicht von vorneherein wissen kann, in welchen Fällen die beiden Analysen dann doch mal was Unterschiedliches ergeben. Aber keine Sorge: der erforderliche Code ist nicht kompliziert; die Ergebnisse interpretiert man analog zu ANOVAs; und es wird hier im Dokument alles erklärt.
Im Rahmen von BSc/MSc-Arbeiten sind v.a. folgende Aspekte wichtig:
Prozentwerte nicht mit ANOVA auswerten, sondern korrekt mit GLMM; wir nutzen dazu das package afex.
Dies zieht andere posthoc-Tests nach sich; wir nutzen dazu das package emmeans.
Ein wesentlicher Punkt bei GLMM für korr/inkorr Daten ist, dass Veränderungen im Bereich nahe 50% (also Zufallsniveau) sehr viel wahrscheinlicher sind als Veränderungen im Beriech 100% (also perfekte Performance). Das heißt, Veränderungen bei GLMM sind nicht linear im Prozentbereich. Stattdessen sind sie linear auf der Skala, die das GLMM benutzt - dies ist eine logit Skala. Wer sich damit nicht genauer auseinandersetzen will, der muss das auch nicht - aber es hat Konsequenzen füt die Darstellung, s.u. Die Nichtlinearität auf der Prozent-Skala wirkt sich darin aus, dass error bars nicht mehr symmetrisch sind.
Wie kann das sein? Das GLMM schätzt lineare Effekte auf der logit-Skala. Dort gibt es dann auch ganz normale, lineare s.d. und s.e. Aber wenn man diese dann in Prozentwerte transformiert, ist diese Transformation nicht linear. Wie man transformiert und plottet, findet sich hier im Skript.
8.2.2 Statistische Auswertung von Metrischen Daten mit ANOVAs und linear mixed models
Kontinuierliche oder metrische Daten werden mit einer ANOVA oder mit einem linear mixed model (LMM) ausgewertet. Dabei ist zu beachten, durch was für eine Art von Design die Daten zustande gekommen sind.
Wenn es sich um eine abhängiges, within-subject Design handelt, darf keine “normale” ANOVA gerechnet werden, sondern es muss eine repeated-measures ANOVA (rmANOVA) gerechnet werden, die berücksichtigt, dass mehre Daten (z.B. aus unterschiedlichen Bedingungen) von der gleichen Person kommen. In einem solchen Fall besteht auch die Möglichkeit, anstatt einer rmANOVA ein LMM zu rechnen.
Der Unterschied zwischen einer rmANOVA und einem LMM bestehen in den sogenannten random effects bzw. random intercepts und random slopes.
Infos: ANOVA vs. LMM
In aller Kürze: ANOVA und LMM sind beides Verfahren, die für uns berechnen, wie stark die uns interessierenden Faktoren (unsere UVs) berücksichtigt werden müssen, um die erhobenen Daten möglichst gut zu schätzen.
Eine ANOVA schätzt dabei die Mittelwerte für jede Faktorstufenkombination der UV-Faktoren. Dabei wird für jede Vp genau derselbe Wert geschätzt.
LMMs ermitteln ebenfalls einen Effekt, der für alle Vpn gleichgilt, schätzen darüber hinaus auch noch individuelle Effekte pro Vp.
Genaueres findet man z.B. hier:
Bei vielen Erklärungen zu LMM werden nicht-psychologische Beispiele verwendet. Um den Zusammenhang zur Psychologie und insbesondere zu within-subject-Designs herzustellen, ist wichtig zu verstehen, dass LMM hierarchische Modelle abbilden. Häufig wird hierzu dann das Beispiel von Schülern in verschiedenen Klassen, die wiederum zu verschiedenen Schulen gehören. Das LMM schätzt, wie die Schule alle ihre Schüler beeinflusst. In unserem Fall der within-subject-Designs ist die hierarchische Ebene die Versuchsperson. Das LMM schätzt, wie die einzelne Vp die Leistung in den Experimentalbedingungen beeinflusst.
LMMs haben gegenüber der rmANOVA einige Vorzüge, dazu zählen
- ein besserer Umgang mit fehlenden Datenpunkten
- eine hohe Flexibilität in der Modellformulierung, bei der von sehr einfachen bis komplexen Modellen alles innerhalb des gleichen Frameworks berechnet werden kann
Für die meisten BSc/MSc-Arbeiten sind diese Dinge nicht so wahnsinnig relevant. Wenn man aber für %-Werte eh schon ein GLMM berechnet, kann man ggf. konsistent sein und auch für die RTs ein LMM berechnen… hierauf bestehen wir aber nicht (nur auf die GLMMs). Bitte besprechen Sie Ihre Statistik immer mit Ihre:r Anleiter:in.
8.2.3 RTs als metrische Daten
In dem hier behandelten Beispieldatensatz sind Reaktionszeiten (RTs) enthalten. Wir werten sie als metrische Daten mit einer ANOVA aus.
Es gibt viel Diskussion darüber, ob RTs so ausgewertet werden sollten. Das liegt v.a. daran, dass RTs nicht normalverteilt sind, sondern rechtsschief - d.h. die Verteilungen haben einen langen Schwanz zu hohen Werten hin. Nach unten kann es so einen Schwanz nicht geben, da RTs nicht <0 sein können. Viele Wissenschaftler:innen befürworten Datentransformationen, die RT-Daten annähernd normalverteilt machen, bspw. eine Log-Transformation. Andere haben gezeigt, dass diese Transformationen nichts bringen oder sogar schädlich für die Analyse sind.
Infos hier: Schramm, P., & Rouder, J. (2019). Are Reaction Time Transformations Really Beneficial? [Preprint]. https://doi.org/10.31234/osf.io/9ksa6
Wir werten RTs in ihrer Rohform aus, d.h. ohne Transformation.
8.2.4 RTs bei korrekten und inkorrekten Trials
Viele berücksichtigen bei der Analyse von RTs nur Trials, in denen die Vp korrekt geantwortet hat. Dies ist aber uneinheitlich und man findet auch Artikel, bei denen einfach alle Trials in RT-Analysen eingehen. Wir gehen hier nicht auf die Argumente für diese beiden Varianten ein. Wir nutzen in unserer Gruppe immer nur korrekte Trials für RT-Analysen.
Bei der Erstellung des Datensatzes für die Analysen der RTs werden also immer die inkorrekten Antworten ausgeschlossen.
8.3 Die “normale” = between-subjects ANOVA
Bei unserem Datensatz handelt es sich um ein within-subjects Design. Daher können wir eine between-subjects ANOVA nicht an diesen Daten demonstrieren.
Ein typisches Anwendungsbeispiel für die between-subjects ANOVA ist die Frage, ob sich drei Gruppen (z.B. 3 Patientengruppen) bezüglich einer metrischen Variable (z.B. ein Testwert) unterscheiden. Jede Vp muss genau in 1 Gruppe enthalten sein; bei Patientengruppen dürfte also niemand 2 Diagnosen haben und dadurch in 2 Gruppen enthalten sein. In diesem Beispiel wird daher analysiert, ob die im Experiment verwendete Reihenfolge der Arme (welcher Arm oben lag) und das Geschlecht (weiblich/männlich) einen Effekt auf die mittlere Reaktionszeit (gemittelt über alle Bedingungen und SOAs) hat.
Unter “Geschlecht” wird dabei verstanden, welche Antwort die Vp auf die Frage gegeben hat; im vorliegenden Fall hat niemand “divers” angegeben. Geschlecht ist eine typische between-subjects-Variable, da jede Vp nur in einer der 2 oder 3 Gruppen einsortiert wird.
8.3.1 Daten vorbereiten
Zunächst generieren wir einen passenden Datensatz. Hierzu muss die Information aus dem Datensatz “bsp_demo” mitverwendet und an den andern Datensatz angehängt werden. Diese Daten wurden bei unserem Beipielexperiment miterhoben, von uns dann aber nicht weiter beachtet. Wir verwenden sie jetzt nur zur Demonstration der ANOVA.
# mittlere Reaktionszeiten der korrekten Antworten berechnen
dAnova <- dc %>% filter(respClean==1) %>% group_by(participant) %>% summarise(n=n(), meanRT=mean(rt))
# Anfangsbuchstaben der Variablennamen klein machen
bsp_demo <- bsp_demo %>% rename_all(tolower)
# Info aus den demographischen Daten hinzufügen
dAnova <- right_join(bsp_demo, dAnova)Joining with `by = join_by(participant)`# Variablen korrekt als Faktoren definieren
dAnova$sex <- factor(dAnova$sex)
dAnova$arm_top <- factor(dAnova$arm_top)8.3.2 Der Anova-Befehl und Output
Als nächstes rechnen wir die ANOVA mit dem Befehl “aov_ez” aus dem afex Paket. und betrachten den Output.
meanRTs.anova <- aov_ez(id = "participant", dv = "meanRT", data = dAnova, between = c("sex", "arm_top"))Contrasts set to contr.sum for the following variables: sex, arm_topNun betrachten wir den Output der ANOVA, den sogenannten ANOVA-Table:
meanRTs.anovaAnova Table (Type 3 tests)
Response: meanRT
Effect df MSE F ges p.value
1 sex 1, 14 17089.40 2.02 .126 .177
2 arm_top 1, 14 17089.40 1.12 .074 .308
3 sex:arm_top 1, 14 17089.40 0.00 <.001 .996
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Wie man der ANOVA-Tabelle entnehmen kann, hat keiner der beiden Faktoren Geschlecht und Armreihenfolge einen signifikanten Effekt auf die Reaktionszeit. Auch die Interaktion ist nicht signifikant.
8.3.3 Post-hoc Vergleiche
Hat einer der Faktoren mehr als 2 Level - bspw. wenn wir Geschlecht in 3 Kategorien codieren - könnten ggf. signifikante Effekte mit Posthoc-Tests weiter untersucht werden. Zu Demonstrationszwecken berechnen wir hier post-hoc Tests, obwohl unsere Faktoren nur 2 Level hatten und nicht signifikant waren. Wir verwenden “emmeans” aus dem gleichnamigen Paket.
Das Paket emmeans ist sehr mächtig, es ist aber auch nicht ganz einfach in der Handhabung. Es lohnt sich für einen vertieften Einstig entweder die Vignette des Pakets unter https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html oder diesen Blogpost https://aosmith.rbind.io/2019/03/25/getting-started-with-emmeans/ zu lesen.
Zunächst schauen wir uns den Haupteffekt von “sex” an:
emmeans(meanRTs.anova, pairwise ~ sex, adjust = "fdr")NOTE: Results may be misleading due to involvement in interactions$emmeans
sex emmean SE df lower.CL upper.CL
female 765 39.6 14 680 850
male 856 49.9 14 749 963
Results are averaged over the levels of: arm_top
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
female - male -90.6 63.7 14 -1.422 0.1768
Results are averaged over the levels of: arm_top Nun schauen wir uns die Interaktion von sex und arm_top an. Wenn ein Faktor mehr als 2 Level hat, müssen p-Werte von post-hoc Tests für multiples Testen korrigiert werden. Als Methoden bieten sich Bonferroni oder False Discovery Rate (fdr) an. Diese Methoden lassen sich bei emmeans einfach anwenden indem sie bei dem “adjust”-Parameter definiert werden.
emmeans(meanRTs.anova, pairwise ~ sex:arm_top, adjust = "fdr")$emmeans
sex arm_top emmean SE df lower.CL upper.CL
female left arm on top 799 53.4 14 684 913
male left arm on top 890 75.5 14 728 1051
female right arm on top 732 58.5 14 606 857
male right arm on top 822 65.4 14 682 962
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio
female left arm on top - male left arm on top -90.9 92.4 14 -0.984
female left arm on top - female right arm on top 67.0 79.2 14 0.847
female left arm on top - male right arm on top -23.3 84.4 14 -0.276
male left arm on top - female right arm on top 158.0 95.5 14 1.655
male left arm on top - male right arm on top 67.7 99.8 14 0.678
female right arm on top - male right arm on top -90.3 87.7 14 -1.030
p.value
0.6108
0.6108
0.7868
0.6108
0.6108
0.6108
P value adjustment: fdr method for 6 tests #alternativer code:
#emm.sex_armtop <- emmeans(meanRTs.anova, ~ sex:arm_top)
#contrast(emm.sex_armtop, method = "pairwise", adjust = "fdr")8.4 Die repeated measures ANOVA
Im Gegensatz zur “normalen” ANOVA, ist eine repeated measures anova (rmANOVA) für Daten geeignet für within-subject-Designs, d.h. Designs mit mehren Messpunkten pro Person.
Im Folgenden wird ein Datensatz erstellt, bei dem zunächst die RT über die SOAs innerhalb jeder Person gemittelt wird. Dann wird der Einfluss von factor_A (gleich/ungleich) und factor_B (gleich/ungleich) und posture (crossed/uncrossed) auf die RT analysiert.
8.4.1 Daten vorbereiten
Da auch die rmANOVA pro Person und Faktorkombination nur einen Wert erwartet, berechnen wir die mittlere Reaktionszeiten der entsprechenden Faktorkombinationen für die Vpn.
drmAnova <- dc %>% filter(respClean==1) %>%
group_by(participant, factor_A, factor_B, posture) %>%
summarise(meanRT=mean(rtClean)) # filter sucht nur die korrekten Trials raus`summarise()` has grouped output by 'participant', 'factor_A', 'factor_B'. You
can override using the `.groups` argument.Bei der rmANOVA sollte man eigentlich sicherstellen, dass es keine fehlenden Werte (missing values) gibt, weil dies zu Fehlern in der Seignifkanzberechnung führen kann. Ein missing value würde bedeuten, dass eine VP für eine Faktorenkombination keinen Mittelwert hat - in einer Bedingung also keine Daten vorliegen. Dies könnte passieren, wenn das Experiment vorzeitig abgebrochen wurde. Ein LMM könnte diese Daten verwenden, die ANOVA hingegen nicht. Eine solche Vp müsste hier also ausgeschlossen werden.
Man kann auf fehlende Werte mit dem Befehl “table” prüfen.
with(drmAnova, table(participant, factor_A, factor_B, posture)) , , factor_B = Level_1, posture = uncr
factor_A
participant Level_1 Level_2
B006b-08 1 1
B006b-09 1 1
B006b-10 1 1
B006b-11 1 1
B006b-12 1 1
B006b-13 1 1
B006b-14 1 1
B006b-16 1 1
B006b-17 1 1
B006b-18 1 1
B006b-19 1 1
B006b-20 1 1
B006b-21 1 1
B006b-22 1 1
B006b-24 1 1
B006b-25 1 1
B006b-26 1 1
B006b-27 1 1
, , factor_B = Level_2, posture = uncr
factor_A
participant Level_1 Level_2
B006b-08 1 1
B006b-09 1 1
B006b-10 1 1
B006b-11 1 1
B006b-12 1 1
B006b-13 1 1
B006b-14 1 1
B006b-16 1 1
B006b-17 1 1
B006b-18 1 1
B006b-19 1 1
B006b-20 1 1
B006b-21 1 1
B006b-22 1 1
B006b-24 1 1
B006b-25 1 1
B006b-26 1 1
B006b-27 1 1
, , factor_B = Level_1, posture = crossed
factor_A
participant Level_1 Level_2
B006b-08 1 1
B006b-09 1 1
B006b-10 1 1
B006b-11 1 1
B006b-12 1 1
B006b-13 1 1
B006b-14 1 1
B006b-16 1 1
B006b-17 1 1
B006b-18 1 1
B006b-19 1 1
B006b-20 1 1
B006b-21 1 1
B006b-22 1 1
B006b-24 1 1
B006b-25 1 1
B006b-26 1 1
B006b-27 1 1
, , factor_B = Level_2, posture = crossed
factor_A
participant Level_1 Level_2
B006b-08 1 1
B006b-09 1 1
B006b-10 1 1
B006b-11 1 1
B006b-12 1 1
B006b-13 1 1
B006b-14 1 1
B006b-16 1 1
B006b-17 1 1
B006b-18 1 1
B006b-19 1 1
B006b-20 1 1
B006b-21 1 1
B006b-22 1 1
B006b-24 1 1
B006b-25 1 1
B006b-26 1 1
B006b-27 1 1Das sieht alles richtig aus.
8.4.2 ANOVA-Befehl und Output
Der wichtige Unterschied zur between-subjects ANOVA sind die “within” Faktoren, also alle Faktoren, die jede Person durchläuft und sich nicht auf einzelne Vp-Gruppen aufteilen. In typischen Experimenten unserer Arbeitsrgruppe sind alle Faktoren within-subject. (Bei uns seltene) Ausnahmen wären z.B.:
- Untersuchung von Kindern verschiedener Altersstufen
- Untersuchung von Geschlechterunterschieden
- Vergleich von Kontroll- und Patient:innengruppen.
meanRTs.rmanova <- aov_ez(id = "participant", dv = "meanRT", data = drmAnova, within = c("factor_A", "factor_B", "posture"))
meanRTs.rmanova #Ergebnisse anzeigenAnova Table (Type 3 tests)
Response: meanRT
Effect df MSE F ges p.value
1 factor_A 1, 17 1209.01 0.05 <.001 .826
2 factor_B 1, 17 10669.09 5.74 * .019 .028
3 posture 1, 17 31739.73 30.85 *** .234 <.001
4 factor_A:factor_B 1, 17 1093.64 1.32 <.001 .266
5 factor_A:posture 1, 17 561.62 0.00 <.001 .973
6 factor_B:posture 1, 17 3091.67 33.92 *** .032 <.001
7 factor_A:factor_B:posture 1, 17 563.88 0.59 <.001 .454
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Wie man in der ANOVA-Tabelle sieht, sind “factor_B”, “posture” und auch die Interaktion aus “faktor_B” und “Posture” signifikant. Wie interpretiert man dies nun?
8.4.3 Plots von estimated marginal means zum Nachvollziehen der Interaktion
Im Folgenden wird eine Abbildung erstellt, die hilft die Ergebnisse visuell zu erfassen. Die Abbildung bricht die Interaktion auf, indem sie die Leistung der Vpn für alle Faktorkombinationen darstellt. Das Paket emmeans kann uns hierbei wieder weiterhelfen. Es ist nämlich nicht nur dazu da, post-hoc Tests zu rechnen, sondern auch, um estimated marginal means (EMM) zu berechnen. EMM sind vom Modell geschätzte Werte der Randhäufigkeiten. Unter letzteren versteht man den Wert der AV für eine Faktorstufe eines Faktors, gemittelt über alle Faktorstufen der anderen Faktoren. Konkret: Habe ich factor_A mit 2 Stufen und factor_B mit 3 Stufen, so ergeben sich eigentlich Werte für 6 Bedingungen: A1_B1, A2_B1, A1_B2 usw. Die Randhäufigkeit von A1 ist der Mittelwert von A1_B1, A1_B2 und A1_B3 - also der vom Modell geschätzte Wert von A1, wenn wir den factor_B wegrechnen.
Beim Herunterbrechen einer Interaktion mittels EMM fängt man normalerweise mit der höchsten signifikanten Interaktion an. In diesem Fall ist das die aus factor_B und posture. Unter “hoch” versteht man dabei, wie viele Faktoren involviert sind - je mehr Faktoren an einer Interaktion beteiligt sind, desto “höher” ist sie.
emm.factorBXposture <- emmeans(meanRTs.rmanova, ~factor_B:posture)
print(emm.factorBXposture) factor_B posture emmean SE df lower.CL upper.CL
Level_1 uncr 729 34.2 17 657 801
Level_2 uncr 716 23.4 17 667 766
Level_1 crossed 840 40.2 17 755 925
Level_2 crossed 935 42.6 17 846 1025
Results are averaged over the levels of: factor_A
Confidence level used: 0.95 Die EMM kann man nun plotten, um sich einen besseren Überblick zu verschaffen:
ggplot(as.data.frame(emm.factorBXposture),
aes(x=factor_B, y=emmean, color=posture)) +
geom_point(position = position_dodge(width=0.5), size=3) +
geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), position = position_dodge(width=0.5), width=0.2) +
ylab(label = "reaction time [ms]")
Man sieht, dass die Reaktionszeiten insgesamt bei der “crossed” posture höher sind, was den im ANOVA-Table sichtbaren Haupteffekt von “posture” reflektiert. Zudem kann man sehen, dass die Reaktionszeit bei “Level_2” für den factor_B bei “crossed” erhöht ist und bei uncrossed nicht (hier ist die RT sogar etwas niedriger). Dieses Verhältnis spiegelt die signifikante Interaktion wieder.
Denken Sie beim Interpretieren einer Interaktion immer an diesen Merksatz: Eine Interaktion bedeutet, dass man den Effekt eines Faktors A nicht benennen kann, ohne sich auch auf den factor_B zu beziehen. Fragt man also, “was ist der Effekt von A”, ist die Antwort “das kommt drauf an, welche Faktorstufe B hat”. Der Ausdruck “es kommt drauf an” kann Ihnen als Signal dienen - wenn er bei der Interpretation eines Faktors vorkommt, handelt es sich um eine Interaktion.
Im vorliegenden Beispiel ist das “kommt drauf an” durch Level_2 bei factor_B begründet. Wenn wir fragen “wie groß ist der Effekt der Handkreuzung”, dann ist die Antwort “da kommt drauf an, ob factor_B Level 1 oder Level 2 hat”.
8.4.4 Post-hoc paarweiser Vergleich für die factor_B x posture Interaktion
Die Interaktion wird nun noch mit post-hoc Tests heruntergebrochen, um zu prüfen, welche der Unterschiede signifikant sind. Unter “Herunterbrechen” versteht man, dass man alle Einzelvergleiche rechnet, die den übergeordneten Effekt - hier die Interaktion - beschreiben.
Man kann auch den Haupteffekt eines einzelnen Faktors herunterbrechen, wenn der Faktor 3 oder mehr Faktorstufen hat. Man rechnet dann paarweise Vergleiche aller Faktorstufen des Faktors.
contrast(emm.factorBXposture, method = "pairwise", adjust = "fdr") contrast estimate SE df t.ratio p.value
Level_1 uncr - Level_2 uncr 12.7 18.7 17 0.682 0.5042
Level_1 uncr - Level_1 crossed -111.0 24.6 17 -4.504 0.0005
Level_1 uncr - Level_2 crossed -206.2 38.6 17 -5.341 0.0002
Level_2 uncr - Level_1 crossed -123.7 29.4 17 -4.204 0.0007
Level_2 uncr - Level_2 crossed -218.9 36.4 17 -6.006 0.0001
Level_1 crossed - Level_2 crossed -95.2 20.4 17 -4.670 0.0004
Results are averaged over the levels of: factor_A
P value adjustment: fdr method for 6 tests Wie man in der Tabelle sehen kann, sind alle Vergleiche zwischen den Stufen signifikant mit Ausnahme des Vergleiches zwischen “Level_1”” und “Level_2”” bei “uncrossed”. Dies bedeutet: obwohl sich die Größe des Kreuzungseffekts unterscheidet, wenn factor_B Level 1 und 2 hat, ist in beiden Fällen der Kreuzungseffekt signifikant. Die Interpretation der Interaktion ist also: es gibt immer einen Kreuzungseffekt, seine Größe kommt aber auf factor_B an.
8.5 Das Linear Mixed Model (LMM)
Grundidee und Vorteile des LMM
Linear Mixed Models (LMM) können alternativ zu ANOVAs verwendet werden. Warum würde man das tun? Die Idee ist, dass LMM die Daten adäquater abbilden (modellieren) können als eine ANOVA.
| ANOVA | LMM |
|---|---|
|
|
Falls Sie es noch nicht getan haben, könnten Sie jetzt die o.g. Quellen zu LMM lesen. Sie können hier allerdings auch einfach weiterlesen, wenn Sie LMM und GLMM einfach nur anwenden wollen, ohne sich weiter darin zu vertiefen - das empfehlen wir zwar nicht, ist aber zeitsparend.
8.5.1 afex-Paket und “mixed”-Befehl
Warum afex statt lme4?
LMM und GLMM können über verschiedene Befehle/Packages/Programme berechnet werden, so wie auch ANOVAs. Das bekannteste R-Package ist lme4 (steht für linear mixed effects v4), das die Befehle lmer und glmer (steht für [generalized] linear mixed effects regression) enthält. Allerdings gibt es bei der Berechnung von (G)LMM mit (g)lmer einige Schwierigkeiten, die nicht direkt sichtbar sind, und die zu Fehlern führen, wenn man sie nicht kennt. Dies hängt damit zusammen, dass man bei Regressionsanalysen die Varianz auf verschiedene Weisen den unterschiedlichen Regressoren - in der ANOVA dann als Faktoren bezeichnet - zuweisen kann. Die Grundeinstellungen von R, und damit auch der meisten LMM Packages in R nutzen andere Methoden als für ANOVAs in der Psychologie typisch. Beachtet man dies nicht, macht man Fehler bei der Interpretation der Ergebnisse. Das Paket afex (steht für: Analysis For EXperiments) hat all diese Schwierigkeiten so behoben, dass man mit ihm (G)LMM ohne weitere Einstellungen so berechnet, wie man es erwarten würde, wenn man aus der Psychologie kommt und ANOVAs berechnet hat.
Wir führen dazu (G)LMMs mit dem Befehl “mixed” aus dem afex-Package aus. mixed stellt alles richtig für uns ein und ruft dann im Hintergrund lmer aus lme4 auf. Der Aufruf erfolgt etwas anders als der von ANOVAs mit ez. Man gibt hier die Faktorstruktur über eine Formel an: ein “+” steht dabei für einen Faktor, ein “:” für eine Interaktion, und ein “*” für Faktoren mit Interaktion.
Hier ein Überblick
| Modellformel | Bedeutung |
|---|---|
| factor_1 + factor_2 | berücksichtigt nur die beiden Haupteffekte |
| factor_1 + factor_2 + factor_1:factor_2 | berücksichtigt Haupteffekte und Interatkion, so wie man es von ANOVAs her gewohnt ist |
| factor_1 * factor_2 | bedeutet dasselbe wie die Zeile drüber, ist aber kürzer hinzuschreiben |
8.5.2 Random effects und ihre Notierung in R bzw. (g)lmer/mixed
Im mixed-Befehl kommen einige Dinge vor, die Sie noch nicht kennen.
Eines ist die Struktur der random effects. Random effects sind in (G)LMM das, was die Variabilität über das hinaus beschreibt, was eine ANOVA kann. Inhaltlich beziehen sich die random effects auf die hierarchisch über den Faktoren liegende Einheit. Oben hatten wir das Beispiel von Schülern in Schulklassen. testen wir z.B. den Effekt von zwei Lehrstilen, zeigen sich vielleicht in jeder Schulklasse spezifische Effekte (z.B. aufgrund der etwas anderen sozialen Situation jeder einzelnen Schulklasse); d.h. der Effekt des Lehrstils unterscheidet sich von Klasse zu Klasse. “Klasse” ist daher die hierarchische Einheit über den Schülern. In unseren within-subjects-Designs haben wir nun experimentelle Faktoren, die sich in jeder Person anders auswirken können. Entsprechend ist “Versuchsperson” die hierarchische Einheit über den Experimentalfaktoren.
Im mixed/lmer-Befehl schreibt man dies als “(1|participant)”. Der Teil hinter dem | bezieht sich also auf die hierarchische Einheit, über die die abhängige Variable variieren kann.
Der Teil vor dem | definiert, über welche Faktoren des Experiments sich die Gruppe hinter dem | laut Modell unterscheiden darf. Die 1 hier bedeutet, dass wir lediglich annehmen, dass sich die Versuchspersonen im Gesamtniveau der abhängigen Variable unterscheiden - hier also, dass jede Vp im Mittel mehr oder weniger schnell ist. Diese Art von random effect nennt man ein random intercept für den Faktor Versuchsperson.
Welche Random Effects bei LMM?
Man könnte auch andere Terme angeben, bspw. (factor_A | participant). Dies würde bedeuten, dass jede Versuchsperson nicht nur ein individuelles Gesamtniveau in der RT hat, sondern dass sich auch der factor_A von Person zu Person unterschiedlich auswirkt. Dies ist meistens eine plausible Annahme: wann möchte man schon davon ausghehen, dass eine experimentelle Manipulation sich bei allen Vpn genau gleich auswirkt? Allerdings macht diese Annahme das Modell komplexer, was zwei Nachteile hat:
- Das Modell dauert länger in der Berechnung. Je nach Komplexität der random factor Struktur kann ein LMM schon mal mehrere Stunden oder sogar Tage rechnen.
- Das Modell kann möglicherweise nicht berechnet (“gefittet”) werden. Dies liegt daran, dass jeder zusätzliche random factor etwas ist, was geschätzt werden muss. Je mehr Parameter man schätzen will, desto mehr Daten benötigt man. Daher kommt es ziemlich oft vor, dass komplexe random structures mit einer Fehlermeldung enden - wenn man Pech hat, nach ziemlich langer Wartezeit. s Es gibt verschiedene Ideen dazu, wie man mit diesem Problem umgehen soll. Die korrekte Formulierung der random effects ist eine große Kontroverse zu (G)LMMs. Im Rahmen von BSc/MSc-Modellen empfehlen wir, die Struktur (also die Formel) der Random Effects mit Ihrer Anleiter*In zu besprechen bzw. sie sich sagen zu lassen. Das heißt, wir verlangen nicht, dass Sie sich mit diesem Thema im Detail auseinandersetzen. Wenn Sie es trotzdem tun wollen, können Sie sich diesen Artikel ansehen:
Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. https://doi.org/10.1016/j.jml.2012.11.001
8.5.3 Methoden zur Berechnung der Freiheitsgrade
Ein weiterer neuer Parameter des Aufrufs vom LMMs ist “method”. Hier wird eine von mehreren zur Verfügung stehenden Methoden für die Schätzung der Freiheitsgrade definiert; dies beeinflusst entsprechend die p-Werte. Es gibt:
- Kenward-Rogers “KR”: langsam
- Satterthwaite “S”: etwas schneller
- Asymptotic: das schnellste, schätzt df auf Inf
Meistens ergeben die verschiedenen Methoden keine großen Unterschiede für die p-Werte. Für den Zweck einer BSc/MSc-Arbeit kann man daher die schnelleren Methoden verwenden. Wir nehmen hier “S”.
8.5.4 LMM in R berechnen
8.5.4.1 Daten vorbereiten
Wie oben beschreiben, wird nicht aggregiert; wir sortieren aber wieder die falschen Antworten raus, indem wir nur Trials mit richtigen Antworten behalten.
dLMM <- dc %>% filter(respClean==1)8.5.4.2 Mixed-Befehl und Output
meanRTs.lmm <- mixed(rtClean~factor_A*factor_B*posture + (1|participant), data = dLMM, method = "S") Contrasts set to contr.sum for the following variables: factor_A, factor_B, posture, participantDen Output schaut man wie zuvor an:
meanRTs.lmmMixed Model Anova Table (Type 3 tests, S-method)
Model: rtClean ~ factor_A * factor_B * posture + (1 | participant)
Data: dLMM
Effect df F p.value
1 factor_A 1, 24290.02 0.58 .444
2 factor_B 1, 24290.20 191.66 *** <.001
3 posture 1, 24290.37 2827.75 *** <.001
4 factor_A:factor_B 1, 24290.01 5.52 * .019
5 factor_A:posture 1, 24290.01 0.03 .858
6 factor_B:posture 1, 24290.07 330.74 *** <.001
7 factor_A:factor_B:posture 1, 24290.01 1.93 .165
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Afex liefert für (G)LMMs eine ANOVA-Tabelle, so dass man die Ergebnisse direkt interpretieren kann wie bei einer ANOVA. Auch hier bekommen wir gezeigt, dass es einen signifikanten Haupteffekt von “factor_B” und “posture” gibt. Und eine Interaktion aus “factor_B” und “posture”, wie bei der rmANOVA. Einziger Unterschied ist hier, dass das LMM auch noch eine (gerade so) signifikante Interaktion zwischen “factor_A” und “factor_B” anzeigt. Der folgende Plot soll aufzeigen, was es mit den Interaktionen auf sich hat.
8.5.4.3 Plots von EMM für die factor_B x posture Interaktion
Wir beginnen zunächst wieder mit der Interaktion aus “factor_B” und “posture”, die wir schon aus der rmANOVA kennen.
Bevor wir die EMMs berechnen, sagen wir R noch, welche Methoden verwendet werden sollen, um die Freiheitsgrade zu schätzen. Das Paket emmeans kann mehrere verwenden, die sich mathematisch unterscheiden und teilweise sehr unterschiedlich lange Berechnungsdauern haben (ähnlich wie oben). Wir wählen die Methode, die “asymptotic” heißt.
#options für df schätzung bei post hoc tests
emm_options(lmer.df = "asymptotic") # options: 'satterthwaite', 'kenward-roger', 'asymptotic'
#Wenn Sie 'satterthwaite' verwenden wollen, müssen Sie das Limit an Tests hochsetzen
#emm_options(lmerTest.limit = 30000)emm.lmm.factorBXPosture <- emmeans(meanRTs.lmm, ~factor_B:posture)NOTE: Results may be misleading due to involvement in interactionsggplot(as.data.frame(emm.lmm.factorBXPosture),
aes(x=factor_B, y=emmean, color=posture)) +
geom_point(position = position_dodge(width=0.5), size=3) +
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL),
position = position_dodge(width=0.5), width=0.2) +
ylab(label = "reaction time [ms]")
Das sieht im Wesentlichen aus wie oben bei der rmANOVA.
Sowohl rmANOVA als auch LMM sind adäquate Methoden. Sie schätzen beide Parameter aus den Daten und berechnen so statistische Effekte. Es ist nicht selten, dass eine rmANOVA und ein LMM äquivalente Ergebnisse bzgl. der statistischen Signifikanz ergeben.
Welches ist nun die richtige Methode?
Es ist ein bisschen so wie wenn Sie im Baumarkt vor einem Regal von Zangen stehen und sich fragen, welche die richtige ist. Alle haben die gleiche Grundfunktion, aber welche sich am besten eignet, hängt vom Einsatzzweck ab. Bei statistischen Methoden unterscheiden sich die Grundannahmen, die den jeweiligen Verfahren zugrunde liegen. Wie oben beschrieben, geht ein LMM von Variabilität auf einer hierarchisch höheren Ebene aus (bei uns: Versuchsperson). Die ANOVA sieht dies einfach nicht vor. Ist man nun der Meinung, dass die Annahmen der LMMs die Wirklichkeit besser beschreiben als die Annahmen der ANOVA, dann setzt man vielleicht besser das LMM ein. Dass bei der rm-Aonva letztlich dasselbe rauskommt, zeigt dann an, dass diese spezifische Annahme offenbar für die statistische Aussage, die wir hier gesucht haben, keine wesentliche Rolle gespielt hat.
Sie kennen diese Thematik auch aus anderen Bereichen der Statistik: bspw. bei der Frage, ob man Verfahren mit oder ohne Verteilungsannahmen verwenden soll - T-test vs. Wilcoxon-Test. Wenn Ihre Daten nicht normalverteilt sind, sollten Sie eigentlich keinen t-Test einsetzen. Man weiß aber, dass der t-Test recht “robust” gegen Verletzungen dieser Annahme ist. D.h., der Unterschied in den Annahmen der beiden Verfahren hat dann in der Praxis keine so große Bedeutung.
8.5.5 Post-hoc paarweiser Vergleich für die factor_B x posture Interaktion
Das Vorgehen für Post-hoc-Tests ist identisch wie bei der rmANOVA:
contrast(emm.lmm.factorBXPosture, method = "pairwise", adjust = "fdr") contrast estimate SE df z.ratio p.value
Level_1 uncr - Level_2 uncr 13.6 4.33 Inf 3.138 0.0017
Level_1 uncr - Level_1 crossed -109.9 4.37 Inf -25.126 <.0001
Level_1 uncr - Level_2 crossed -210.3 4.49 Inf -46.884 <.0001
Level_2 uncr - Level_1 crossed -123.5 4.39 Inf -28.142 <.0001
Level_2 uncr - Level_2 crossed -223.9 4.50 Inf -49.779 <.0001
Level_1 crossed - Level_2 crossed -100.4 4.53 Inf -22.152 <.0001
Results are averaged over the levels of: factor_A
Degrees-of-freedom method: asymptotic
P value adjustment: fdr method for 6 tests Auch die Vergleiche fallen von der Interpretation her ähnlich aus wie oben bei der rmANOVA. Es wird aber auch der erste Vergleich zwischen Level_1 und Level_2 bei “uncrossed” signifikant. Insgesamt sind die p-Werte überall niedriger. Hier drückt sich die etwas höhere Power des LMM aus, das alle Trials mit einbezieht (die rmANOVA war ja auf der Basis von Mittelwerten berechnet worden).
8.5.6 Plots von EMM für die factor_B x factor_A Interaktion
Im Folgenden soll noch einmal die Abhängigkeit von “factor_A”von “factor_B” verdeutlicht werden (das war die “gerade so” signifikante Interaktion, die die rmANOVA nicht zeigte).
emm.lmm.factorAXfactorB <- emmeans(meanRTs.lmm, ~factor_A:factor_B)NOTE: Results may be misleading due to involvement in interactionsggplot(as.data.frame(emm.lmm.factorAXfactorB),
aes(x=factor_B, y=emmean, shape=factor_A)) +
geom_point(position = position_dodge(width=0.5), size=3) +
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL),
position = position_dodge(width=0.5), width=0.2) +
ylab(label = "reaction time [ms]")
Die Stufen von “factor_A” über “factor_B” verlaufen in leicht unterschiedliche Richtungen; dies verursacht die Interaktion. Auch hier kann man wieder einen post-hoc-Test rechnen.
8.5.6.1 Post-hoc paarweiser Vergleich für die factor_B x factor_A Interaktion
Das Vorgehen ist wiederum analog zu oben.
contrast(emm.lmm.factorAXfactorB, method = "pairwise", adjust = "fdr") contrast estimate SE df z.ratio p.value
Level_1 Level_1 - Level_2 Level_1 4.97 4.37 Inf 1.136 0.2559
Level_1 Level_1 - Level_1 Level_2 -36.05 4.43 Inf -8.144 <.0001
Level_1 Level_1 - Level_2 Level_2 -45.81 4.43 Inf -10.332 <.0001
Level_2 Level_1 - Level_1 Level_2 -41.02 4.43 Inf -9.252 <.0001
Level_2 Level_1 - Level_2 Level_2 -50.78 4.44 Inf -11.435 <.0001
Level_1 Level_2 - Level_2 Level_2 -9.76 4.49 Inf -2.172 0.0358
Results are averaged over the levels of: posture
Degrees-of-freedom method: asymptotic
P value adjustment: fdr method for 6 tests 8.6 Das Generalized Linear Mixed Model (GLMM)
Grundidee und Vorteile des GLMM
Beim GLMM wird deutlich, warum (G)LMM auf einzelne Trials zurückgreifen statt auf Mittelwerte. Die Grundidee einer Regression ist, dass die abhängige Variable eines spezfischen Falles (hier: Trials) durch die Prädiktoren (in der Psychologie oft kategoriale Faktoren, aber auch kontinuierliche Variablen) vorhergesagt wird. Im GLMM wird nun die Wahrscheinlichkeitsverteilung, die dieser Zuweisung zugrunde liegt, flexibilisiert. Im normalen LMM sagen wir immer eine kontinuierliche abhängige Variable vorher, die normalverteilt ist (dies ist zumindest die Annahme, die das Modell zwangsweise macht). IM GLMM kann man andere Verteilungen annehmen.
Hier wird es für Fehler- bzw. Korrektheitsmaße interessant. Diese sind ja meistens binär: man hat einen Trial richtig oder falsch gemacht; links oder rechts gedrückt; etc. Das GLMM sagt daher nun keine kontinuierliche Variable vorher, sondern - für jeden Trial - ob man aufgrund der Prädkiktoren eher “richtig” oder eher “falsch” erwarten sollte. Dabei übersetzt sich “eher” in eine Wahrscheinlichkeit. Daher verwendet man im Falle von Fehlermaßen eine Binomialverteilung für die Regression.
Das bedeutet: Ob jemand den Trial richtig oder falsch beantwortet, ist wie ein Münzwurf. Aber je nach Experimentalbedingungen (also Faktorausprägungen) tendiert die Münze stärker oder schwächer in die eine oder andere Richtung.
Machen Sie sich klar, dass diese Konzeption völlig anders ist als wenn man über alle Trials einer Bedingung aggregiert (z.B. via Mittelwert) und dann so tut, als sei der resultierende %-korrekt-Wert normalverteilt und daher per ANOVA analysierbar.
8.6.1 Berechnung eines GLMM in R
8.6.1.1 Methoden zur Schätzung der p-Werte
In der Anwendung ist der wesentliche Unterschied im Code zu dem LMM, dass wir hier einen Parameter setzen: family=“binomial”. Das sagt dem “mixed”-Befehl, dass es sich um ein GLMM handelt mit einer binären AV. Zusätzlich müssen wir auch eine Methode wählen, um die p-Werte zu schäzen.
Hier hat man die Wahl zwischen:
- Parametric Bootstrap (“PB”): das ist eine sichere Wahl, die eine gute Schätzung liefert für unterschiedlichste Modelle bzw. Faktorstrukturen. Der Nachteil ist, dass die Schätzung je nach Modell sehr lange dauern kann (tw. mehrere Stunden).
- Likelihood-Ratio-Test (“LRT”): Dieses Verfahren geht deutlich schneller, aber die Schätzungen können ungenau werden, wenn eine einfache random-effect-Struktur gewählt wurde (siehe ?mixed für Informationen) dazu.
Da GLMMs besonders bei “PB” eine deutlich längere Berechnungsdauer als andere Modelle haben, lohnt es sich, mehrere Prozessoren gleichzeitig zu nutzen, und es macht unter diesen Umständen Sinn, die Ergebnisse abzuspeichern, damit man das Modell nicht erneut berechnen musss, wenn man die R-Session neu startet.
8.6.1.2 Der GLMM-Befehl und Output
Wir zeigen hier ein Beispiel, wie das Fitten mit mehreren Prozessorkernen und das Abspeichern der Daten geht:
library( parallel ) #für Multicore-Prozesse
nc <- detectCores() - 1 #Anzahl an Prozessoren minus 1 - der eine verbleibende erlaubt uns dann, am Rechner weiter normal zu arbeiten
#Um zu sehen, was der Prozess macht, wird Output in eine Textdatei geschrieben.
cl <- makeCluster(rep("localhost", nc), outfile = "cl1.log.txt")
startTime <- Sys.time() #Anfangszeit, um zu sehen, wie lange das Modell läuft
resp.glmm <- mixed( respClean ~ posture * factor_A * factor_B + ( 1 | participant ),
data = dc, family = "binomial", method = "PB",
args_test=list(nsim = 500), check_contrasts = TRUE, cl = cl ) # family =
endTime <- Sys.time() #Endzeit für den Vergleich mit der Anfangszeit
( endTime- startTime )
#Modellergebnisse abspeichern
save(resp.glmm, file = "resp_glmm.RData")Hier wählen wir noch einen Ansatz mit der “LRT”-Methode:
resp.glmm <- mixed( respClean ~ posture * factor_A * factor_B + ( 1 | participant ),
data = dc, family = "binomial", method = "LRT",
check_contrasts = TRUE ) Contrasts set to contr.sum for the following variables: posture, factor_A, factor_B, participantAls nächstes schauen wir den Output (ANOVA Table) an:
resp.glmmMixed Model Anova Table (Type 3 tests, LRT-method)
Model: respClean ~ posture * factor_A * factor_B + (1 | participant)
Data: dc
Df full model: 9
Effect df Chisq p.value
1 posture 1 649.80 *** <.001
2 factor_A 1 1.63 .201
3 factor_B 1 150.84 *** <.001
4 posture:factor_A 1 0.02 .876
5 posture:factor_B 1 14.56 *** <.001
6 factor_A:factor_B 1 0.31 .576
7 posture:factor_A:factor_B 1 1.94 .164
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1So ähnlich wie bei den RTs sehen wir auch hier einen Haupteffekt von “posture” und “factor_B”, sowie eine signifikante Interaktion “posture:factor_B”.
8.6.1.3 EMM für das GLMM
Auch hier kann man diese Interaktion am besten mit einem Plot und post-hoc Tests nachvollziehen.
Grundsätzlich geht das ähnlich wie bei dem LMM, aber es gibt einen Unterschied: das GLMM hat die Daten auf eine logit-Skala transformiert, die etwas kompliziert nachzuvollziehen ist. Der logit bildet die sog. Log-Odds ab, d.h. den Logarithmus des Verhältnisses der beiden Events. Ist in einer Bedingung die Wahrscheinlichkeit, richtig zu antworten, 2x so hoch wie falsch zu antworten, ergeben sich die Odds als 2:1. Davon zieht man dann noch den Logarithmus. Entsprechend hat man eher kein Gefühl dafür, was die Zahlen, die bei dieser Methode rauskommen, wirklich bedeuten.
Wir können aber in “emmeans” definieren, dass wir die Daten wieder auf die Einheit der Antwort (also zwischen 0 und 1) zurücktransformieren wollen. Als Resultat bekommen wird Werte, die Schätzungen für den Anteil korrekter Antworten (den Wert 1 in unserer Variable) in den Faktorkombinationen sind. Dazu defnieren wir in “emmeans” einfach das Argument: type=“response”. In anderen Worten, hier wird die GLMM-Berechnung für uns am Schluss in %-Werte umgerechnet.
Die sich ergebenden %-Werte sind völlig andere, als hätten wir eine ANOVA mit kumulierten %-Werten pro Bedingung durchgeführt!
emm.resp.glmm <- emmeans(resp.glmm, ~posture:factor_B, type = "response")NOTE: Results may be misleading due to involvement in interactionsprint(emm.resp.glmm) posture factor_B prob SE df asymp.LCL asymp.UCL
uncr Level_1 0.985 0.00390 Inf 0.975 0.991
crossed Level_1 0.952 0.01170 Inf 0.923 0.971
uncr Level_2 0.976 0.00619 Inf 0.960 0.985
crossed Level_2 0.884 0.02615 Inf 0.822 0.926
Results are averaged over the levels of: factor_A
Confidence level used: 0.95
Intervals are back-transformed from the logit scale 8.6.1.4 Plot der EMM aus dem GLMM
ggplot(as.data.frame(emm.resp.glmm),
aes(x=factor_B, y=prob, color=posture)) +
geom_point(position = position_dodge(width=0.5), size=3) +
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL),
position = position_dodge(width=0.5), width=0.2) +
ylab(label = "proportion of correct response")
8.6.2 Post-hoc paarweise Vergleiche für die Interaktion im GLMM
Für die paarweisen Vergleiche kann man identisch vorgehen wie zuvor:
contrast(emm.resp.glmm, method = "pairwise", adjust = "fdr") contrast odds.ratio SE df null z.ratio p.value
uncr Level_1 / crossed Level_1 3.310 0.3135 Inf 1 12.636 <.0001
uncr Level_1 / uncr Level_2 1.646 0.1706 Inf 1 4.808 <.0001
uncr Level_1 / crossed Level_2 8.670 0.7727 Inf 1 24.236 <.0001
crossed Level_1 / uncr Level_2 0.497 0.0407 Inf 1 -8.533 <.0001
crossed Level_1 / crossed Level_2 2.619 0.1618 Inf 1 15.591 <.0001
uncr Level_2 / crossed Level_2 5.268 0.3961 Inf 1 22.096 <.0001
Results are averaged over the levels of: factor_A
P value adjustment: fdr method for 6 tests
Tests are performed on the log odds ratio scale Hier wird dann allerdings die Odds Ratio dargestellt (ohne log), also wie oben beschrieben das Verhältnis von zwei Wahrscheinlichkeiten. Eine korrekte Antwort ist also bei uncr Level_1 3.3 mal wahrscheinlicher als bei crossed Level_1, und das ist statistisch signifikant.