Próbuję przejść z używania ez
pakietu do lme
ANOVA 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 aov
i lme
(sprawdziłem go i robią), to nie jest sprawa dla moich danych. Chociaż wartości F są podobne, nie są takie same.
aov
zwraca wartość f 1,3399, lme
zwraca 1,36264. Jestem gotów zaakceptować aov
wynik jako „poprawny”, ponieważ to również zwraca SPSS (i to jest ważne dla mojego pola / przełożonego).
Pytania:
Byłoby wspaniale, gdyby ktoś mógł wyjaśnić, dlaczego ta różnica istnieje i jak mogę wykorzystać ją
lme
do zapewnienia wiarygodnych wyników. (Byłbym również skłonny użyćlmer
zamiastlme
tego typu rzeczy, jeśli daje to „poprawny” wynik. Jednak do tej pory go nie używałem).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))
lme
wynikach standardowego podręcznika ANOVA (podanego przezaov
i 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. PonadtoAnova()
(odcar
) zwraca tylko wartości Chi² dlalme
obiektów. Dlatego dla mnie na moje pierwsze pytanie jeszcze nie ma odpowiedzi.lme
; ale w przypadku kontrastówglht
działa również nalm
pasowania, a nie tylkolme
pasowania. (Równieżlme
wyniki są również standardowymi wynikami podręczników.)lm
analizy powtarzanego pomiaru. Tylkoaov
radzi sobie z powtarzanymi pomiarami, ale zwróci obiekt klasyaovlist
, który nie jest jeszcze obsługiwany przezglht
.lm
wykorzystuje błąd resztkowy jako termin błędu dla wszystkich efektów; kiedy występują efekty, które powinny używać innego terminu błędu,aov
jest konieczne (lub zamiast tego, używając wyników zlm
do ręcznego obliczenia statystyk F). W twoim przykładzie termin błędu dla interakcjifactor
jestid:factor
interakcją, która jest wartością błędu resztkowego w modelu addytywnym. Porównaj swoje wyniki zanova(lm(value~factor+id))
.Odpowiedzi:
Różnią się one, ponieważ model lme wymusza, aby składnik wariancji
id
był 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.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
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.
Podaje odchylenia standardowe, podnosząc do kwadratu, aby uzyskać wydajność wariancji
9.553e-10
dla wariancji id i0.1586164
wariancji resztkowej.Teraz powinieneś wiedzieć, że stosowanie
aov
powtarzanych 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 korzystanialme
naaov
to, ż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
lme
z 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ącgls
:Ten stół ANOVA zgadza się ze stołem od
aov
dopasowania i odlm
dopasowania.Ok, i co? Cóż, jeśli uważasz, że wariancja
id
i korelacja między obserwacjami powinny być nieujemne,lme
dopasowanie jest w rzeczywistości bardziej odpowiednie niż dopasowanie przy użyciuaov
lub,lm
ponieważ jego oszacowanie wariancji rezydualnej jest nieco lepsze. Jednakże, jeśli uważasz, że korelacja pomiędzy obserwacjami może być ujemna,aov
lublm
czygls
jest lepszy.Możesz być także zainteresowany dalszym badaniem struktury korelacji; aby spojrzeć na ogólną strukturę korelacji, zrobiłbyś coś takiego
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ą: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.
źródło
lme
celu uzyskania takich samych wyników jak w przypadkuaov
(a tym samym włączenialme
dla wszystkich ANOVA), a mianowicie użycie argumentu korelacji w wezwaniu dolme
:anova(lme(value ~ factor, data = tau.base, random = ~1|id, correlation = corCompSymm(form = ~1|id)))
aov
,lme
bez związku symetrii, ilme
ze związkiem symetrii) mają dokładnie taką samą liczbę DFS.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ą.aov()
pasuje do modelu przylm()
użyciu najmniejszych kwadratów,lme
pasuje 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.
źródło
lme
lublmer
obliczenie 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?