Generalnie zgadzam się z analizą Bena, ale dodam kilka uwag i odrobinę intuicji.
Po pierwsze, ogólne wyniki:
- Wyniki lmerTest przy użyciu metody Satterthwaite są prawidłowe
- Metoda Kenwarda-Rogera jest również poprawna i zgadza się z Satterthwaite
Ben przedstawia projekt, w którym subnumjest zagnieżdżony w groupjednocześnie direction
i group:directionsą skrzyżowane z subnum. Oznacza to, że termin błędu naturalnego (tj. Tak zwana „obejmująca warstwa błędu”) dla groupjest, subnumpodczas gdy obejmującą warstwą błędu dla innych terminów (w tym subnum) są reszty.
Strukturę tę można przedstawić na tak zwanym schemacie struktury czynnikowej:
names <- c(expression("[I]"[5169]^{5191}),
expression("[subnum]"[18]^{20}), expression(grp:dir[1]^{4}),
expression(dir[1]^{2}), expression(grp[1]^{2}), expression(0[1]^{1}))
x <- c(2, 4, 4, 6, 6, 8)
y <- c(5, 7, 5, 3, 7, 5)
plot(NA, NA, xlim=c(2, 8), ylim=c(2, 8), type="n", axes=F, xlab="", ylab="")
text(x, y, names) # Add text according to ’names’ vector
# Define coordinates for start (x0, y0) and end (x1, y1) of arrows:
x0 <- c(1.8, 1.8, 4.2, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 7, 5, 5, 3, 7)
x1 <- c(2.7, 2.7, 5, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 7, 7, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length=0.1)
Tutaj losowe terminy są ujęte w nawiasy kwadratowe, 0reprezentują ogólną średnią (lub punkt przecięcia), [I]reprezentują termin błędu, liczby superskryptów to liczba poziomów, a liczby subskryptów to liczba stopni swobody przy założeniu zrównoważonego projektu. Wykres wskazuje, że składnikiem błędu naturalnego (obejmującego warstwę błędu) dla groupjest subnumi że licznik df dla subnum, który jest równy mianownikowi df dla group, wynosi 18: 20 minus 1 df dla groupi 1 df dla średniej ogólnej. Bardziej kompleksowe wprowadzenie do diagramów struktury czynników jest dostępne w rozdziale 2 tutaj: https://02429.compute.dtu.dk/eBook .
Gdyby dane były dokładnie zbalansowane, bylibyśmy w stanie skonstruować testy F z rozkładu SSQ, jak zapewnia anova.lm. Ponieważ zestaw danych jest bardzo ściśle zrównoważony, możemy uzyskać przybliżone testy F w następujący sposób:
ANT.2 <- subset(ANT, !error)
set.seed(101)
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
fm <- lm(rt ~ group * direction + subnum, data=ANT.2)
(an <- anova(fm))
Analysis of Variance Table
Response: rt
Df Sum Sq Mean Sq F value Pr(>F)
group 1 994365 994365 200.5461 <2e-16 ***
direction 1 1568 1568 0.3163 0.5739
subnum 18 7576606 420923 84.8927 <2e-16 ***
group:direction 1 11561 11561 2.3316 0.1268
Residuals 5169 25629383 4958
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tutaj wszystkie wartości F i p są obliczane przy założeniu, że wszystkie terminy mają reszty jako otaczającą warstwę błędów, i jest to prawdą dla wszystkich oprócz „grupy”. Zamiast tego „zrównoważony-poprawny” test F dla grupy to:
F_group <- an["group", "Mean Sq"] / an["subnum", "Mean Sq"]
c(Fvalue=F_group, pvalue=pf(F_group, 1, 18, lower.tail = FALSE))
Fvalue pvalue
2.3623466 0.1416875
gdzie używamy subnumMS zamiast ResidualsMS w mianowniku wartości F.
Zauważ, że te wartości dość dobrze pasują do wyników Satterthwaite:
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
anova(model, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 12065.3 12065.3 1 18 2.4334 0.1362
direction 1951.8 1951.8 1 5169 0.3936 0.5304
group:direction 11552.2 11552.2 1 5169 2.3299 0.1270
Pozostałe różnice wynikają z tego, że dane nie są dokładnie zrównoważone.
OP porównuje anova.lmz anova.lmerModLmerTest, co jest w porządku, ale aby porównać jak z podobnym, musimy użyć tych samych kontrastów. W tym przypadku istnieje różnica między anova.lmi anova.lmerModLmerTestponieważ domyślnie produkują one odpowiednio testy typu I i III, a dla tego zestawu danych istnieje (mała) różnica między kontrastami typu I i III:
show_tests(anova(model, type=1))$group
(Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment 0 1 0.005202759 0.5013477
show_tests(anova(model, type=3))$group # type=3 is default
(Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment 0 1 0 0.5
Gdyby zestaw danych był całkowicie zrównoważony, kontrasty typu I byłyby takie same jak kontrasty typu III (na które nie ma wpływu obserwowana liczba próbek).
Ostatnia uwaga jest taka, że „powolność” metody Kenwarda-Rogera nie jest spowodowana ponownym dopasowaniem modelu, ale ponieważ wiąże się z obliczeniami z macierzą marginalnych wariancji-kowariancji obserwacji / reszt (w tym przypadku 5191 x 5191) przypadek metody Satterthwaite.
Dotyczy modelu 2
Jak dla MODEL2 sytuacja staje się bardziej złożona i myślę, że łatwiej jest rozpocząć dyskusję z innym modelem gdzie mam obejmowała „klasyczną” interakcji między subnumi direction:
model3 <- lmer(rt ~ group * direction + (1 | subnum) +
(1 | subnum:direction), data = ANT.2)
VarCorr(model3)
Groups Name Std.Dev.
subnum:direction (Intercept) 1.7008e-06
subnum (Intercept) 4.0100e+01
Residual 7.0415e+01
Ponieważ wariancja związana z interakcją jest zasadniczo zerowa (w obecności subnumlosowego efektu głównego), interakcja nie ma wpływu na obliczenie mianownika stopni swobody, wartości F i wartości p :
anova(model3, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 12065.3 12065.3 1 18 2.4334 0.1362
direction 1951.8 1951.8 1 5169 0.3936 0.5304
group:direction 11552.2 11552.2 1 5169 2.3299 0.1270
Jednak subnum:directionjest załączając stratum błędu subnumwięc jeśli usuniemy subnumwszystkie związane SSQ spada z powrotem dosubnum:direction
model4 <- lmer(rt ~ group * direction +
(1 | subnum:direction), data = ANT.2)
Teraz naturalny określenie błędu group, directiona group:directionto
subnum:directioni nlevels(with(ANT.2, subnum:direction))= 40 i czterech parametrów stopni swobody mianownika dla tych składników powinna wynosić około 36:
anova(model4, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 24004.5 24004.5 1 35.994 4.8325 0.03444 *
direction 50.6 50.6 1 35.994 0.0102 0.92020
group:direction 273.4 273.4 1 35.994 0.0551 0.81583
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Te testy F można także aproksymować za pomocą testów F „zrównoważonych-poprawnych” :
an4 <- anova(lm(rt ~ group*direction + subnum:direction, data=ANT.2))
an4[1:3, "F value"] <- an4[1:3, "Mean Sq"] / an4[4, "Mean Sq"]
an4[1:3, "Pr(>F)"] <- pf(an4[1:3, "F value"], 1, 36, lower.tail = FALSE)
an4
Analysis of Variance Table
Response: rt
Df Sum Sq Mean Sq F value Pr(>F)
group 1 994365 994365 4.6976 0.0369 *
direction 1 1568 1568 0.0074 0.9319
group:direction 1 10795 10795 0.0510 0.8226
direction:subnum 36 7620271 211674 42.6137 <2e-16 ***
Residuals 5151 25586484 4967
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
teraz przechodzę na model2:
model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
Ten model opisuje dość skomplikowaną strukturę kowariancji losowego efektu z macierzą wariancji-kowariancji 2x2. Domyślna parametryzacja nie jest łatwa do opanowania, a lepszym rozwiązaniem jest ponowna parametryzacja modelu:
model2 <- lmer(rt ~ group * direction + (0 + direction | subnum), data = ANT.2)
Jeśli porównamy model2to model4, mają one równie wiele losowych efektów; 2 na każdy subnum, tj. Łącznie 2 * 20 = 40. Chociaż model4określa pojedynczy parametr wariancji dla wszystkich 40 efektów losowych, model2zastrzega, że każda subnumpara efektów losowych ma dwu-zmienny rozkład normalny z macierzą wariancji-kowariancji 2x2, której parametry są podane przez
VarCorr(model2)
Groups Name Std.Dev. Corr
subnum directionleft 38.880
directionright 41.324 1.000
Residual 70.405
Wskazuje to na nadmierne dopasowanie, ale zachowajmy to na inny dzień. Ważną rzeczą jest to, że model4jest specjalnym przypadkiem model2 i że modeljest również szczególnym przypadkiem model2. Luźno (i intuicyjnie) mówienie (direction | subnum)zawiera lub rejestruje wariacje związane z głównym efektem, subnum a także interakcją direction:subnum. W kategoriach efektów losowych możemy uznać te dwa efekty lub struktury za przechwytujące odpowiednio różnice między wierszami i wierszami po kolumnach:
head(ranef(model2)$subnum)
directionleft directionright
1 -25.453576 -27.053697
2 16.446105 17.479977
3 -47.828568 -50.835277
4 -1.980433 -2.104932
5 5.647213 6.002221
6 41.493591 44.102056
W tym przypadku zarówno oszacowania efektu losowego, jak i oszacowania parametru wariancji wskazują, że tak naprawdę mamy tutaj tylko losowy efekt główny subnum(zmienność między wierszami). Wszystko to prowadzi do tego, że mianownik Satterthwaite posiada stopnie swobody
anova(model2, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 12059.8 12059.8 1 17.998 2.4329 0.1362
direction 1803.6 1803.6 1 125.135 0.3638 0.5475
group:direction 10616.6 10616.6 1 125.136 2.1418 0.1458
jest kompromisem między tymi głównymi strukturami efektu i interakcji: Grupa DenDF pozostaje na 18 (zagnieżdżona w subnumprojekcie), ale directioni
group:directionDenDF są kompromisami między 36 ( model4) a 5169 ( model).
Nie sądzę, aby cokolwiek tutaj wskazuje, że przybliżenie Satterthwaite (lub jego implementacja w lmerTest ) jest błędna.
Tabela równoważna z metodą Kenwarda-Rogera daje
anova(model2, type=1, ddf="Ken")
Type I Analysis of Variance Table with Kenward-Roger's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group 12059.8 12059.8 1 18.000 2.4329 0.1362
direction 1803.2 1803.2 1 17.987 0.3638 0.5539
group:direction 10614.7 10614.7 1 17.987 2.1414 0.1606
Nic dziwnego, że KR i Satterthwaite mogą się różnić, ale dla wszystkich praktycznych celów różnica w wartościach p jest niewielka. Moja analiza powyżej wskazuje, że DenDFfor directioni group:directionnie powinien być mniejszy niż ~ 36 i prawdopodobnie większy niż ten, biorąc pod uwagę, że w zasadzie mamy tylko losowy główny efekt directionteraźniejszości, więc jeśli cokolwiek myślę, to jest wskazanie, że metoda KR staje się DenDFzbyt niska w tym przypadku. Pamiętaj jednak, że dane tak naprawdę nie obsługują (group | direction)struktury, więc porównanie jest trochę sztuczne - byłoby bardziej interesujące, gdyby model był obsługiwany.
ezAnovaostrzeżenie, ponieważ nie powinieneś uruchamiać anova 2x2, jeśli w rzeczywistości twoje dane pochodzą z projektu 2x2x2.