Ich analysiere die Ergebnisse eines Reaktionszeitversuchs in R.
Ich führte eine ANOVA mit wiederholten Messungen durch (1 Faktor innerhalb des Subjekts mit 2 Stufen und 1 Faktor zwischen den Subjekten mit 2 Stufen). Ich habe ein ähnliches lineares gemischtes Modell ausgeführt und wollte die früheren Ergebnisse in Form einer ANOVA-Tabelle mit zusammenfassen lmerTest::anova
.
Verstehen Sie mich nicht falsch: Ich habe nicht die identischen Ergebnisse erwartet, bin mir jedoch nicht sicher, welche Freiheitsgrade die lmerTest::anova
Ergebnisse haben. Es scheint mir eher eine ANOVA ohne Aggregation auf Subjektebene zu reflektieren.
Mir ist bewusst, dass die Berechnung von Freiheitsgraden in Modellen mit gemischten Effekten schwierig ist, aber lmerTest::anova
im aktualisierten ?pvalues
Thema ( lme4
Paket) als eine mögliche Lösung erwähnt wird .
Ist diese Berechnung korrekt? Entsprechen die Ergebnisse von lmerTest::anova
korrekt dem angegebenen Modell?
Update: Ich habe die individuellen Unterschiede größer gemacht. Die Freiheitsgrade lmerTest::anova
unterscheiden sich stärker von einfachen Anova, aber ich bin mir immer noch nicht sicher, warum sie für den Faktor / die Interaktion innerhalb des Subjekts so groß sind.
# mini example with ANT dataset from ez package
library(ez); library(lme4); library(lmerTest)
# repeated measures ANOVA with ez package
data(ANT)
ANT.2 <- subset(ANT, !error)
# update: make individual differences larger
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum),
within = .(direction), between = .(group))
anova.ez
# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
# simple ANOVA on all available data
m <- lm(rt ~ group * direction, data = ANT.2)
anova(m)
Ergebnisse des obigen Codes [ aktualisiert ]:
anova.ez
$ ANOVA
Effect DFn DFd F p p<.05 ges
2 group 1 18 2.6854464 0.11862957 0.1294475137
3 direction 1 18 0.9160571 0.35119193 0.0001690471
4 group:direction 1 18 4.9169156 0.03970473 * 0.0009066868
lmerTest :: anova (Modell)
Analysis of Variance Table of type 3 with Satterthwaite
approximation for degrees of freedom
Df Sum Sq Mean Sq F value Denom Pr(>F)
group 1 13293 13293 2.6830 18 0.1188
direction 1 1946 1946 0.3935 5169 0.5305
group:direction 1 11563 11563 2.3321 5169 0.1268
Anova (m)
Analysis of Variance Table
Response: rt
Df Sum Sq Mean Sq F value Pr(>F)
group 1 1791568 1791568 242.3094 <2e-16 ***
direction 1 728 728 0.0985 0.7537
group:direction 1 12024 12024 1.6262 0.2023
Residuals 5187 38351225 7394
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
quelle
ezAnova
Warnung, da Sie 2x2 anova nicht ausführen sollten, wenn Ihre Daten tatsächlich aus 2x2x2-Design stammen.ez
umformuliert werden. Es besteht eigentlich aus zwei Teilen, die wichtig sind: (1) dass Daten aggregiert werden und (2) Informationen zu Teilentwürfen. # 1 ist für die Diskrepanz am relevantesten, da es erklärt, dass man, um eine traditionelle Anova mit nicht gemischten Effekten durchzuführen, die Daten zu einer einzigen Beobachtung pro Zelle des Entwurfs aggregieren muss. In diesem Fall möchten wir eine Beobachtung pro Subjekt pro Ebene der Variablen "Richtung" (unter Beibehaltung der Gruppenbezeichnungen für Subjekte). ezANOVA berechnet dies automatisch.summary(aov(rt ~ group*direction + Error(subnum/direction), data=ANT.2))
und es gibt 16 (?) Dfs fürgroup
und 18 fürdirection
undgroup:direction
. Die Tatsache, dass es ~ 125 Beobachtungen pro Gruppen- / Richtungskombination gibt, ist für RM-ANOVA ziemlich irrelevant, siehe z. B. meine eigene Frage stats.stackexchange.com/questions/286280 : Die Richtung wird sozusagen gegen das Subjekt getestet. Richtungswechselwirkung.lmerTest
dieser Schätzung von ~ 125 dfs nicht zu.lmerTest::anova(model2, ddf="Kenward-Roger")
18.000 df fürgroup
und17.987
df für die beiden anderen Faktoren zurück, was in hervorragender Übereinstimmung mit RM-ANOVA (gemäß ezAnova) steht. Mein Fazit ist, dass Satterthwaites Annäherung ausmodel2
irgendeinem Grund fehlschlägt .Ich stimme Bens Analyse im Allgemeinen zu, aber lassen Sie mich ein paar Bemerkungen und ein wenig Intuition hinzufügen.
Zunächst die Gesamtergebnisse:
Ben skizziert das Design, in das währenddessen
subnum
verschachtelt ist und mit dem gekreuzt wird . Dies bedeutet, dass der natürliche Fehlerterm (dh die sogenannte "einschließende Fehlerschicht") für ist, während die einschließende Fehlerschicht für die anderen Begriffe (einschließlich ) die Residuen sind.group
direction
group:direction
subnum
group
subnum
subnum
Diese Struktur kann in einem sogenannten Faktorstrukturdiagramm dargestellt werden:
Hier sind zufällige Terme in Klammern eingeschlossen,
0
stellen den Gesamtmittelwert (oder Achsenabschnitt) dar,[I]
stellen den Fehlerterm dar, die Superskriptnummern sind die Anzahl der Ebenen und die Unterskriptnummern sind die Anzahl der Freiheitsgrade, die ein ausgeglichenes Design voraussetzen. Das Diagramm zeigt, dass der natürliche Fehlerterm (einschließlich der Fehlerschicht) fürgroup
istsubnum
und dass der Zähler df fürsubnum
, der dem Nenner df für entsprichtgroup
, 18: 20 minus 1 df fürgroup
und 1 df für den Gesamtmittelwert beträgt . Eine umfassendere Einführung in Faktorstrukturdiagramme finden Sie in Kapitel 2 hier: https://02429.compute.dtu.dk/eBook .Wenn die Daten genau ausgeglichen wären, könnten wir die F-Tests aus einer SSQ-Zerlegung wie von konstruieren
anova.lm
. Da der Datensatz sehr genau ausgewogen ist, können wir ungefähre F-Tests wie folgt erhalten:Hier werden alle F- und p- Werte unter der Annahme berechnet, dass alle Terme die Residuen als ihre einschließende Fehlerschicht haben, und das gilt für alle außer 'Gruppe'. Der 'ausgeglichen-korrekte' F- Test für die Gruppe lautet stattdessen:
wo wir die
subnum
MS anstelle derResiduals
MS im F- Wert-Nenner verwenden.Beachten Sie, dass diese Werte recht gut mit den Satterthwaite-Ergebnissen übereinstimmen:
Die verbleibenden Unterschiede sind darauf zurückzuführen, dass die Daten nicht genau ausgewogen sind.
Das OP vergleicht
anova.lm
mitanova.lmerModLmerTest
, was in Ordnung ist, aber um Gleiches mit Gleichem zu vergleichen, müssen wir die gleichen Kontraste verwenden. In diesem Fall gibt es einen Unterschied zwischenanova.lm
undanova.lmerModLmerTest
da sie standardmäßig Tests vom Typ I bzw. III erzeugen, und für diesen Datensatz gibt es einen (kleinen) Unterschied zwischen den Kontrasten vom Typ I und III:Wenn der Datensatz vollständig ausgeglichen gewesen wäre, wären die Kontraste vom Typ I dieselben gewesen wie die Kontraste vom Typ III (die von der beobachteten Anzahl von Proben nicht beeinflusst werden).
Eine letzte Bemerkung ist, dass die "Langsamkeit" der Kenward-Roger-Methode nicht auf eine Modellanpassung zurückzuführen ist, sondern auf Berechnungen mit der marginalen Varianz-Kovarianz-Matrix der Beobachtungen / Residuen (in diesem Fall 5191x5191), was nicht der Fall ist der Fall für Satterthwaites Methode.
Zum Modell2
Was Modell 2 betrifft, wird die Situation komplexer und ich denke, es ist einfacher, die Diskussion mit einem anderen Modell zu beginnen, in das ich die 'klassische' Interaktion zwischen
subnum
und aufgenommen habedirection
:Da die mit der Wechselwirkung verbundene Varianz im Wesentlichen Null ist (bei Vorhandensein des
subnum
zufälligen Haupteffekts), hat der Wechselwirkungsterm keinen Einfluss auf die Berechnung der Nennerfreiheitsgrade , F- Werte und p- Werte:Jedoch
subnum:direction
ist die umschließende Fehler Stratum fürsubnum
so , wenn wir entfernensubnum
alle zugehörigen SSQ fällt zurück insubnum:direction
Nun ist die natürliche Fehlerterm für
group
,direction
undgroup:direction
istsubnum:direction
und mitnlevels(with(ANT.2, subnum:direction))
= 40 und vier Parametern die Nenner - Freiheitsgrade für diese Begriffe etwa 36 sein sollten:Diese F- Tests können auch mit den 'ausgeglichen-korrekten' F- Tests angenähert werden:
Wenden wir uns nun Modell2 zu:
Dieses Modell beschreibt eine ziemlich komplizierte Kovarianzstruktur mit zufälligen Effekten mit einer 2x2-Varianz-Kovarianz-Matrix. Die Standardparametrierung ist nicht einfach zu handhaben und wir sind besser mit einer Neuparametrisierung des Modells:
Wenn wir vergleichen
model2
zumodel4
, haben sie ebenso viele Zufallseffekte; 2 für jedensubnum
, dh 2 * 20 = 40 insgesamt. Währendmodel4
für alle 40 zufälligen Effekte einmodel2
einzelnersubnum
Varianzparameter festgelegt ist , ist festgelegt, dass jedes Paar zufälliger Effekte eine bi-variable Normalverteilung mit einer 2x2-Varianz-Kovarianz-Matrix aufweist, deren Parameter durch gegeben sindDies deutet auf eine Überanpassung hin, aber bewahren wir das für einen anderen Tag auf. Der wichtige Punkt hier ist , dass
model4
ist ein Sonderfallmodel2
und dasmodel
ist auch ein Sonderfallmodel2
. Das lose (und intuitive) Sprechen(direction | subnum)
enthält oder erfasst die Variation, die mit dem Haupteffektsubnum
sowie der Interaktion verbunden istdirection:subnum
. In Bezug auf die zufälligen Effekte können wir uns diese beiden Effekte oder Strukturen so vorstellen, dass sie Variationen zwischen Zeilen bzw. zeilenweise erfassen:In diesem Fall zeigen diese zufälligen Effektschätzungen sowie die Varianzparameterschätzungen beide, dass wir hier wirklich nur einen zufälligen Haupteffekt von
subnum
(Variation zwischen Zeilen) haben. Dies alles führt dazu, dass Satterthwaite Nenner Freiheitsgrade inist ein Kompromiss zwischen diesen Haupteffekt- und Interaktionsstrukturen: Die Gruppe DenDF bleibt bei 18 (
subnum
vom Design her verschachtelt ), aber diedirection
undgroup:direction
DenDF sind Kompromisse zwischen 36 (model4
) und 5169 (model
).Ich glaube, hier deutet nichts darauf hin, dass die Satterthwaite-Näherung (oder ihre Implementierung in lmerTest ) fehlerhaft ist.
Die äquivalente Tabelle mit der Kenward-Roger-Methode gibt
Es ist nicht überraschend, dass KR und Satterthwaite unterschiedlich sein können, aber für alle praktischen Zwecke ist der Unterschied in den p- Werten winzig. Meine obige Analyse zeigt, dass das
DenDF
fürdirection
undgroup:direction
nicht kleiner als ~ 36 und wahrscheinlich größer als das sein sollte, da wir im Grunde nur den zufälligen Haupteffekt derdirection
Gegenwart haben. Wenn ich also etwas denke, ist dies ein Hinweis darauf, dass die KR-MethodeDenDF
zu niedrig wird in diesem Fall. Beachten Sie jedoch, dass die Daten die(group | direction)
Struktur nicht wirklich unterstützen , sodass der Vergleich etwas künstlich ist - es wäre interessanter, wenn das Modell tatsächlich unterstützt würde.quelle
model3
. Es wird jedochsubnum:direction
als Fehlerbegriff zum Testen verwendetdirection
. Während Sie dies hier nur durch Ausschließen(1|subnum)
wie in erzwingen könnenmodel4
. Warum? (2b) Außerdem ergibt RM-ANOVA df = 18 fürdirection
und nicht 36, wenn Sie einsteigenmodel4
. Warum?summary(aov(rt ~ group*direction + Error(subnum/direction), data=ANT.2))
.subnum
undsubnum:direction
nicht Null sind, dannanova(lm(rt2 ~ group * direction + subnum + subnum:direction, data = ANT.2))
ergibt 18 df für alle drei Faktoren, und das ist es, was die KR-Methode aufgreift. Dies ist bereits zu sehen,model3
wenn KR die designbasierten 18 df für alle Terme angibt, auch wenn die Interaktionsvarianz Null ist, während Satterthwaite den verschwindenden Varianzterm erkennt und den df entsprechend anpasst ....