Dlaczego lme i aov zwracają różne wyniki dla ANOVA z powtarzanymi pomiarami w R?

Próbuję przejść z używania ezpakietu do lmeANOVA dla powtarzanych pomiarów (ponieważ mam nadzieję, że będę w stanie używać niestandardowych kontrastów lme).

Zgodnie z radą z tego wpisu na blogu udało mi się skonfigurować ten sam model, używając zarówno aov(jak to się dzieje ez, kiedy jest to wymagane) i lme. Jednak, podczas gdy w przykładzie podanym w tym poście , że F -values nie zgadzają się idealnie między aovi lme(sprawdziłem go i robią), to nie jest sprawa dla moich danych. Chociaż wartości F są podobne, nie są takie same.

aovzwraca wartość f 1,3399, lmezwraca 1,36264. Jestem gotów zaakceptować aovwynik jako „poprawny”, ponieważ to również zwraca SPSS (i to jest ważne dla mojego pola / przełożonego).

Pytania:

1. Byłoby wspaniale, gdyby ktoś mógł wyjaśnić, dlaczego ta różnica istnieje i jak mogę wykorzystać ją lmedo zapewnienia wiarygodnych wyników. (Byłbym również skłonny użyć lmerzamiast lmetego typu rzeczy, jeśli daje to „poprawny” wynik. Jednak do tej pory go nie używałem).

2. Po rozwiązaniu tego problemu chciałbym przeprowadzić analizę kontrastu. Byłbym szczególnie zainteresowany kontrastem łączenia dwóch pierwszych poziomów czynnika (tj. c("MP", "MT")) I porównania tego z trzecim poziomem czynnika (tj "AC".). Ponadto, testowanie trzeciego w porównaniu do czwartego poziomu współczynnika (tj. W "AC"porównaniu do "DA").

Dane:

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

I kod:

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))

Różnią się one, ponieważ model lme wymusza, aby składnik wariancji idbył większy od zera. Patrząc na surową tabelę anova dla wszystkich terminów, widzimy, że średni błąd kwadratu dla id jest mniejszy niż dla reszt.

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131   

Kiedy obliczamy składniki wariancji, oznacza to, że wariancja wynikająca z id będzie ujemna. Moja pamięć oczekiwanej pamięci średnich kwadratów jest chwiejna, ale obliczenia są podobne

(0.15052-0.16131)/3 = -0.003597.

Brzmi dziwnie, ale może się zdarzyć. Oznacza to, że średnie dla każdego identyfikatora są bliżej siebie, niż można by się spodziewać, biorąc pod uwagę resztkową zmienność w modelu.

Natomiast użycie lme zmusza tę wariancję do większej niż zero.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

Podaje odchylenia standardowe, podnosząc do kwadratu, aby uzyskać wydajność wariancji 9.553e-10dla wariancji id i 0.1586164wariancji resztkowej.

Teraz powinieneś wiedzieć, że stosowanie aovpowtarzanych miar jest właściwe tylko wtedy, gdy uważasz, że korelacja między wszystkimi parami powtarzanych miar jest identyczna; nazywa się to symetrią złożoną. (Technicznie kulistość jest wymagane, ale jest to wystarczające do teraz). Jednym z powodów do korzystania lmena aovto, że może obsługiwać różne rodzaje struktur korelacji.

W tym konkretnym zestawie danych oszacowanie tej korelacji jest ujemne; pomaga to wyjaśnić, w jaki sposób średni błąd kwadratu dla id był mniejszy niż błąd resztkowy do kwadratu. Negatywna korelacja oznacza, że ​​jeśli pierwszy pomiar u danej osoby byłby poniżej średniej, jego drugi byłby średnio powyżej średniej, dzięki czemu średnie średnie dla osób byłyby mniej zmienne, niż można by się spodziewać, gdyby istniała korelacja zerowa lub korelacja dodatnia.

Używanie lmez efektem losowym jest równoważne dopasowaniu złożonego modelu symetrii, w którym korelacja ta musi być nieujemna; możemy dopasować model, w którym korelacja może być ujemna, używając gls:

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

Ten stół ANOVA zgadza się ze stołem od aovdopasowania i od lmdopasowania.

Ok, i co? Cóż, jeśli uważasz, że wariancja idi korelacja między obserwacjami powinny być nieujemne, lmedopasowanie jest w rzeczywistości bardziej odpowiednie niż dopasowanie przy użyciu aovlub, lmponieważ jego oszacowanie wariancji rezydualnej jest nieco lepsze. Jednakże, jeśli uważasz, że korelacja pomiędzy obserwacjami może być ujemna, aovlub lmczy glsjest lepszy.

Możesz być także zainteresowany dalszym badaniem struktury korelacji; aby spojrzeć na ogólną strukturę korelacji, zrobiłbyś coś takiego

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Tutaj ograniczam tylko dane wyjściowe do struktury korelacji. Wartości od 1 do 4 reprezentują cztery poziomy factor; widzimy, że czynnik 1 i czynnik 4 mają dość silną korelację ujemną:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

Jednym ze sposobów wyboru między tymi modelami jest test współczynnika wiarygodności; pokazuje to, że model efektów losowych i model ogólnej struktury korelacji nie różnią się statystycznie znacząco; gdy tak się dzieje, zwykle preferowany jest prostszy model.

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965

aov()pasuje do modelu przy lm()użyciu najmniejszych kwadratów, lmepasuje do maksymalnego prawdopodobieństwa. Ta różnica w sposobie szacowania parametrów modelu liniowego prawdopodobnie odpowiada (bardzo małej) różnicy w wartościach F.

W praktyce (np. W przypadku testowania hipotez) szacunki te są takie same, więc nie rozumiem, w jaki sposób można uznać za „bardziej wiarygodny” niż drugi. Pochodzą z różnych modeli dopasowania modelu.

W przypadku kontrastów musisz skonfigurować matrycę kontrastu dla swoich czynników. Venebles i Ripley pokazują, jak to zrobić na s. 143, s. 144 i s. 293–294 4. edycji.