Analizuję wyniki eksperymentu czasu reakcji w R.

Przeprowadziłem ANOVA z powtarzanymi pomiarami (1 czynnik wewnątrz podmiotu z 2 poziomami i 1 czynnik między podmiotem z 2 poziomami). Uruchomiłem podobny liniowy model mieszany i chciałem podsumować wyniki lmera w formie tabeli ANOVA lmerTest::anova.

Nie zrozumcie mnie źle: nie spodziewałem się identycznych wyników, jednak nie jestem pewien co do stopnia swobody lmerTest::anovawyników. Wydaje mi się, że raczej odzwierciedla ANOVA bez agregacji na poziomie przedmiotu.

Zdaję sobie sprawę z tego, że obliczanie stopni swobody w modelach z efektami mieszanymi jest trudne, ale lmerTest::anovajest wspomniane jako jedno z możliwych rozwiązań w zaktualizowanym ?pvaluestemacie ( lme4pakiecie).

Czy to obliczenie jest prawidłowe? Czy wyniki lmerTest::anovapoprawnie odzwierciedlają określony model?

Aktualizacja: Zwiększiłem różnice indywidualne. Stopnie swobody lmerTest::anovaróżnią się bardziej od prostych anova, ale wciąż nie jestem pewien, dlaczego są tak duże dla czynnika / interakcji wewnątrz podmiotu.

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

Wyniki powyższego kodu [ zaktualizowane ]:

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

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

Myślę , że lmerTestto robi dobrze i ezanovaźle robi się w tym przypadku.

• wyniki są lmerTestzgodne z moją intuicją / zrozumieniem
• lmerTestzgadzają się dwa różne obliczenia w (Satterthwaite i Kenward-Roger)
• oni również się z tym zgadzają nlme::lme
• kiedy go uruchamiam, ezanovadaje ostrzeżenie, którego nie do końca rozumiem, ale którego nie należy lekceważyć ...

Przykład ponownego uruchomienia:

library(ez); library(lmerTest); library(nlme)
data(ANT)
ANT.2 <- subset(ANT, !error)
set.seed(101)  ## for reproducibility
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

Opracuj eksperymentalny projekt

with(ANT.2,table(subnum,group,direction))

Wygląda więc na to, że osoby ( subnum) są umieszczone w grupach kontrolnych lub leczonych, a każda z nich jest testowana w obu kierunkach - tj. Kierunek może być testowany w obrębie osób (mianownik df jest duży), ale grupa i grupa: kierunek można badać tylko pomiędzy osoby fizyczne

(anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
    within = .(direction), between = .(group)))
## $ANOVA
##            Effect DFn DFd         F          p p<.05          ges
## 2           group   1  18 2.4290721 0.13651174       0.1183150147
## 3       direction   1  18 0.9160571 0.35119193       0.0002852171
## 4 group:direction   1  18 4.9169156 0.03970473     * 0.0015289914

Tutaj otrzymuję Warning: collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate. mianownik DF wygląda trochę funky (wszystkie równe 18): Myślę, że powinny być większe dla kierunku i grupy: kierunek, który można testować niezależnie (ale byłby mniejszy, gdybyś dodał (direction|subnum)do modelu)?

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
##                 Df  Sum Sq Mean Sq F value Denom Pr(>F)
## group            1 12065.7 12065.7  2.4310    18 0.1364
## direction        1  1952.2  1952.2  0.3948  5169 0.5298
## group:direction  1 11552.2 11552.2  2.3299  5169 0.1270

Dfkolumna odnosi się tutaj do df liczniku Denom(drugi przedostatniej) przedstawia szacowaną mianownik df; zgadzają się z klasyczną intuicją. Co ważniejsze, otrzymujemy również różne odpowiedzi na wartości F ...

Możemy również dwukrotnie sprawdzić za pomocą Kenwarda-Rogera ( bardzo wolno, ponieważ wymaga to kilkukrotnego dopełnienia modelu)

lmerTest::anova(model,ddf="Kenward-Roger")

Wyniki są identyczne.

W tym przykładzie lme(z nlmepaczki) faktycznie wykonuje idealnie dobrą pracę zgadując odpowiedni mianownik df (wartości F i p są bardzo nieznacznie różne):

model3 <- lme(rt ~ group * direction, random=~1|subnum, data = ANT.2)
anova(model3)[-1,]
##                 numDF denDF   F-value p-value
## group               1    18 2.4334314  0.1362
## direction           1  5169 0.3937316  0.5304
## group:direction     1  5169 2.3298847  0.1270

Jeśli dopasuję interakcję między directioni subnumdf dla directioni group:directionsą one znacznie mniejsze (pomyślałbym, że będą mieć 18 lat, ale może coś źle rozumiem):

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
lmerTest::anova(model2)
##                 Df  Sum Sq Mean Sq F value   Denom Pr(>F)
## group            1 20334.7 20334.7  2.4302  17.995 0.1364
## direction        1  1804.3  1804.3  0.3649 124.784 0.5469
## group:direction  1 10616.6 10616.6  2.1418 124.784 0.1459

Generalnie zgadzam się z analizą Bena, ale dodam kilka uwag i odrobinę intuicji.

Po pierwsze, ogólne wyniki:

1. Wyniki lmerTest przy użyciu metody Satterthwaite są prawidłowe
2. 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.