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

24

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))
Henrik
źródło
Wygląda na to, że właśnie odpowiedziałeś na część dotyczącą kontrastów w swojej odpowiedzi tutaj ; jeśli nie, edytuj to pytanie, abyśmy wiedzieli, jakie trudności pozostaną.
Aaron - Przywróć Monikę
2
@Aaron, o ile istnieją różnice w lmewynikach standardowego podręcznika ANOVA (podanego przez aovi właśnie tego potrzebuję), nie jest to dla mnie opcja. W mojej pracy chcę zgłosić ANOVA, a nie coś takiego jak ANOVA. Co ciekawe Venables i Ripley (2002, s. 285) pokazują, że oba podejścia prowadzą do identycznych szacunków. Ale różnice w wartościach F powodują, że czuję się źle. Ponadto Anova()(od car) zwraca tylko wartości Chi² dla lmeobiektów. Dlatego dla mnie na moje pierwsze pytanie jeszcze nie ma odpowiedzi.
Henrik,
Rozumiem (ale nie podzielam się) twoją ostrożnością lme; ale w przypadku kontrastów glhtdziała również na lmpasowania, a nie tylko lmepasowania. (Również lmewyniki są również standardowymi wynikami podręczników.)
Aaron - Przywróć Monikę
Niestety nie można określić lmanalizy powtarzanego pomiaru. Tylko aovradzi sobie z powtarzanymi pomiarami, ale zwróci obiekt klasy aovlist, który nie jest jeszcze obsługiwany przez glht.
Henrik,
3
lmwykorzystuje błąd resztkowy jako termin błędu dla wszystkich efektów; kiedy występują efekty, które powinny używać innego terminu błędu, aovjest konieczne (lub zamiast tego, używając wyników z lmdo ręcznego obliczenia statystyk F). W twoim przykładzie termin błędu dla interakcji factorjest id:factorinterakcją, która jest wartością błędu resztkowego w modelu addytywnym. Porównaj swoje wyniki z anova(lm(value~factor+id)).
Aaron - Przywróć Monikę

Odpowiedzi:

28

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
Aaron - Przywróć Monikę
źródło
2
W rzeczywistości możliwe jest użycie symetrii złożonej w lmecelu uzyskania takich samych wyników jak w przypadku aov(a tym samym włączenia lmedla wszystkich ANOVA), a mianowicie użycie argumentu korelacji w wezwaniu do lme:anova(lme(value ~ factor, data = tau.base, random = ~1|id, correlation = corCompSymm(form = ~1|id)))
Henrik
1
Niezłe znalezisko. Ale czy nie ma w tym dopasowaniu dodatkowego parametru? Ma trzy parametry wariancji; wariancja dla id, rezydualna wariancja i korelacja, podczas gdy gls ma tylko resztkową wariancję i korelację.
Aaron - Przywróć Monikę
1
Twój argument wydaje się wiarygodny, jednak wyniki się nie zgadzają. Wszystkie tabele ANOVA ( aov, lmebez związku symetrii, i lmeze związkiem symetrii) mają dokładnie taką samą liczbę DFS.
Henrik
1
Musisz mnie przekonać, że te trzy parametry są naprawdę nadparametryzacją pierwszych dwóch. Czy wiesz, jak są ze sobą powiązane?
Aaron - Przywróć Monikę
1
Nie. Ufam wynikom anova.lme(). Z twojej odpowiedzi wynika, że ​​związek między ANOVA a modelami mieszanymi leży w ich strukturze korelacji. Przeczytałem wtedy, że narzucenie złożonej symetrycznej struktury korelacji prowadzi do równości między dwoma podejściami. Dlatego narzuciłem to. Nie mam pojęcia, czy to zjada kolejny df. Wynik nie zgadza się jednak z tą interpretacją.
Henrik
2

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.

Chris
źródło
Hmm, ale dlaczego więc czasami występują różnice, a czasem wyniki są dokładnie takie same? Ponadto wydaje się, że niemożliwe jest użycie lmelub lmerobliczenie ANOVA (ściśle mówiąc), ponieważ wykorzystuje ona metodę podobną, ale nie identyczną. Czy zatem nie ma sposobu obliczania kontrastów dla ANOVA z powtarzanym pomiarem w R?
Henrik,
Jeśli układ, w którym modelujesz, jest naprawdę liniowy niż najmniejszych kwadratów, a ML powinien dawać tę samą statystykę F. Dopiero gdy w danych jest inna struktura, obie metody dadzą różne wyniki. Pinheiro i Bates opisują to w książce o modelach z efektami mieszanymi. Poza tym prawdopodobnie nie są „dokładnie” równe, gdybyś poszedł wystarczająco daleko cyframi sig, jestem pewien, że znajdziesz różnicę. Ale dla wszystkich praktycznych celów są one takie same.
Chris,