Moje pytanie opiera się na tej odpowiedzi, która wykazała, który lme4::lmer
model odpowiada dwukierunkowej ANOVA z powtarzanymi pomiarami:
require(lme4)
set.seed(1234)
d <- data.frame(
y = rnorm(96),
subject = factor(rep(1:12, 4)),
a = factor(rep(1:2, each=24)),
b = factor(rep(rep(1:2, each=12))),
c = factor(rep(rep(1:2, each=48))))
# standard two-way repeated measures ANOVA:
summary(aov(y~a*b+Error(subject/(a*b)), d[d$c == "1",]))
# corresponding lmer call:
anova(lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), d[d$c == "1",]))
Moje pytanie dotyczy teraz, jak rozszerzyć to na przypadek trójdrożnej ANOVA:
summary(aov(y~a*b*c+Error(subject/(a*b*c)), d))
## [...]
## Error: subject:a:b:c
## Df Sum Sq Mean Sq F value Pr(>F)
## a:b:c 1 0.101 0.1014 0.115 0.741
## Residuals 11 9.705 0.8822
Naturalne rozszerzenie oraz jego wersje nie pasują do wyników ANOVA:
anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject), d))
## [...]
## a:b:c 1 0.1014 0.1014 0.1500
anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) +
(1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))
## [...]
## a:b:c 1 0.1014 0.1014 0.1539
Należy zauważyć, że bardzo podobne pytanie zostało zadane wcześniej . Brakowało jednak przykładowych danych (które podano tutaj).
y ~ a*b + (1 + a*b|subject), d[d$c == "1",]
? A może coś mi brakuje?lmer
narzeka, ponieważ losowe efekty nie są już identyfikowane. Początkowo myślałem również, że to model, który chcę, ale tak nie jest. Jeśli porównasz lmer model, który proponuję dla skrzynki 2-drożnej ze standardową ANOVA, zobaczysz, że wartości F dokładnie pasują. Jak powiedziałem w odpowiedzi, połączyłem.lmer
napisany przez ciebie model (który wyklucza losowe interakcje dwukierunkowe) będzie równoważny 3-kierunkowej RM-ANOVA, ale drugi napisany przez ciebie model (który obejmuje losowy dwustronne interakcje) powinny być. Jeśli chodzi o to, dlaczego istnieje rozbieżność nawet w przypadku tego modelu, mam przeczucie, na czym polega problem, pójdę na obiad, a następnie przejrzę zestaw danych z zabawkami.Odpowiedzi:
Bezpośrednią odpowiedzią na twoje pytanie jest to, że ostatni model, który napisałeś,
Uważam, że jest „w zasadzie” poprawny, chociaż jest to dziwna parametryzacja, która nie zawsze wydaje się dobrze sprawdzać w praktyce.
Jeśli chodzi o to, dlaczego dane wyjściowe z tego modelu są niezgodne z danymi
aov()
wyjściowymi, myślę, że istnieją dwa powody.lmer()
nie pozwalają ze sobą modele mieszane (i większość innych programów modeli mieszanych).Pozwól mi najpierw zademonstrować preferowaną parametryzację na twoim początkowym dwukierunkowym przykładzie ANOVA. Załóżmy, że Twój zestaw danych
d
jest załadowany. Twój model (zwróć uwagę, że zmieniłem kody zastępcze na kontrastowe) to:które działało tutaj dobrze, ponieważ pasowało do
aov()
wyniku. Model, który preferuję, obejmuje dwie zmiany: ręczne kodowanie kontrastu czynników, abyśmy nie pracowali z obiektami współczynnika R (co zalecam robić w 100% przypadków) i inne określenie efektów losowych:Te dwa podejścia są całkowicie równoważne w prostym problemie dwukierunkowym. Teraz przejdziemy do problemu trójstronnego. Wspomniałem wcześniej, że podany przykładowy zestaw danych był patologiczny. Tak więc, co chcę zrobić przed zaadresowaniem twojego przykładowego zestawu danych, to najpierw wygenerować zestaw danych z rzeczywistego modelu komponentów wariancji (tj. Gdzie niezerowe komponenty wariancji są wbudowane w prawdziwy model). Najpierw pokażę, w jaki sposób moja preferowana parametryzacja wydaje się działać lepiej niż ta, którą zaproponowałeś. Wtedy będę wykazać w inny sposób szacowania składników wariancji, która ma nie narzucać, że muszą być nieujemne. Wtedy będziemy w stanie zobaczyć problem z oryginalnym przykładowym zestawem danych.
Nowy zestaw danych będzie identyczny pod względem struktury, z tym wyjątkiem, że będziemy mieli 50 tematów:
Współczynniki F, które chcemy dopasować, to:
Oto nasze dwa modele:
Jak widzimy, tylko druga metoda pasuje do wyniku
aov()
, chociaż pierwsza metoda jest co najmniej na boisku. Druga metoda zapewnia również wyższe prawdopodobieństwo logarytmiczne. Nie jestem pewien, dlaczego te dwie metody dają różne wyniki, ponieważ znów uważam, że są one „w zasadzie” równoważne, ale być może dzieje się tak z kilku powodów liczbowych / obliczeniowych. A może się mylę i nawet w zasadzie nie są one równoważne.Teraz pokażę inny sposób szacowania składników wariancji w oparciu o tradycyjne pomysły ANOVA. Zasadniczo weźmiemy oczekiwane równania średnich kwadratów dla twojego projektu, podstawimy obserwowane wartości średnich kwadratów i rozwiążemy składowe wariancji. Aby uzyskać oczekiwane średnie kwadraty, użyjemy funkcji R, którą napisałem kilka lat temu, zwanej
EMS()
, co jest udokumentowane TUTAJ . Poniżej zakładam, że funkcja jest już załadowana.OK, teraz wrócimy do oryginalnego przykładu. Współczynniki F, które próbujemy dopasować, to:
Oto nasze dwa modele:
W tym przypadku dwa modele dają w zasadzie te same wyniki, chociaż druga metoda ma bardzo nieznacznie większe prawdopodobieństwo logarytmiczne. Żadna z metod nie pasuje
aov()
. Ale spójrzmy na to, co otrzymujemy, gdy rozwiązujemy dla składników wariancji, jak to zrobiliśmy powyżej, stosując procedurę ANOVA, która nie ogranicza składników wariancji, aby były nieujemne (ale które można stosować tylko w zrównoważonych projektach bez ciągłych predyktorów i bez brakujące dane; klasyczne założenia ANOVA).Teraz możemy zobaczyć, co jest patologiczne w oryginalnym przykładzie. Model najlepiej dopasowany to taki, który sugeruje, że kilka losowych składników wariancji jest ujemnych. Ale
lmer()
(i większość innych programów modeli mieszanych) ogranicza szacunki składników wariancji do wartości nieujemnych. Jest to ogólnie uważane za rozsądne ograniczenie, ponieważ odchylenia oczywiście nigdy naprawdę nie mogą być ujemne. Jednak konsekwencją tego ograniczenia jest to, że modele mieszane nie są w stanie dokładnie przedstawić zestawów danych, które wykazują ujemne korelacje wewnątrzklasowe, to znaczy zestawów danych, w których obserwacje z tego samego skupienia są mniejsze(raczej niż bardziej) podobny średnio niż obserwacje losowe z zestawu danych, aw konsekwencji, gdy wariancja wewnątrz klastra znacznie przekracza wariancję między klastrami. Takie zestawy danych są całkowicie uzasadnionymi zestawami danych, które czasami można spotkać w świecie rzeczywistym (lub przypadkowo symulować!), Ale nie można ich rozsądnie opisać w modelu wariancji-składników, ponieważ implikują negatywne składniki wariancji. Można je jednak opisać „nie rozsądnie” przez takie modele, jeśli oprogramowanie na to pozwoli.aov()
pozwala na to.lmer()
nie.źródło
I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons
- czy może to zrozumieć lepiej (dwa lata później)? Próbowałem dowiedzieć się, jaka jest różnica, ale też nie rozumiem ...A
ma poziomów następnie ocenia tylko jeden parametr wariancji, natomiast szacunki parametrów wariancji i korelacji między nimi. Rozumiem, że klasyczna ANOVA szacuje tylko jeden parametr wariancji, a zatem powinna być równoważna pierwszej metodzie, a nie drugiej. Dobrze? Ale dla obie metody szacują jeden parametr i nadal nie jestem pewien, dlaczego się nie zgadzają. k - 1 k ( k - 1 ) / 2 k = 2(1|A:sub)
(0+A|sub)
lmer
wywołania dają identyczneanova()
wyjście, wariancje efektów losowych są jednak całkiem różne: patrzVarCorr(mod1)
iVarCorr(mod2)
. Nie do końca rozumiem, dlaczego tak się dzieje; czy ty? Dlamod3
imod4
widać, że cztery z siedmiu wariancji dlamod3
są w rzeczywistości równe zeru (dlamod4
wszystkich siedmiu są niezerowe); ta „osobliwość”mod3
jest prawdopodobnie powodem, dla którego tabele anova różnią się. Poza tym, jak użyłbyś swojego „preferowanego sposobu”, gdybya
ib
miał więcej niż dwa poziomy?Są
a
,b
,c
stałe lub losowe efekty? Jeśli zostaną naprawione, twoja składnia będzie po prostuźródło
subject
, dla wszystkich efektów (tjWithin
.). Zobacz Projekt eksperymentalny: Procedury dla nauk behawioralnych (2013) autorstwa Kirka, rozdział 10 (str. 458) lub mój post tutajlmer
? Mimo to dostanę moją kopię Kirka (tylko druga edycja) i zobaczę, co mówi.lmer
modeli. Najlepszym sposobem sprawdzenia dopasowania modelu jest sprawdzenie ich dfs,lmerTest
ponieważ przybliżenie KR ma dać ciexact
dfs, a zatem wartości p.