Sprawdzanie założeń mieszanych modeli lmer / lme w R.

Przeprowadziłem powtarzający się projekt, w którym przetestowałem 30 mężczyzn i 30 kobiet w trzech różnych zadaniach. Chcę zrozumieć, jak różni się zachowanie mężczyzn i kobiet i jak to zależy od zadania. Użyłem zarówno pakietu lmer, jak i lme4, aby to zbadać, jednak utknąłem przy próbie sprawdzenia założeń dla każdej z metod. Kod, który uruchamiam, to

lm.full <- lmer(behaviour ~ task*sex + (1|ID/task), REML=FALSE, data=dat)
lm.full2 <-lme(behaviour ~ task*sex, random = ~ 1|ID/task, method="ML", data=dat)

Sprawdziłem, czy interakcja była najlepszym modelem, porównując ją z prostszym modelem bez interakcji i uruchamiając anova:

lm.base1 <- lmer(behaviour ~ task+sex+(1|ID/task), REML=FALSE, data=dat)
lm.base2 <- lme(behaviour ~ task+sex, random= ~1|ID/task), method="ML", data=dat)
anova(lm.base1, lm.full)
anova(lm.base2, lm.full2)

P1: Czy można używać tych predyktorów jakościowych w liniowym modelu mieszanym?
P2: Czy dobrze rozumiem, że dobrze jest, że zmienna wynikowa („zachowanie”) nie musi być normalnie rozkładana sama (na płeć / zadania)?
P3: Jak mogę sprawdzić jednorodność wariancji? Używam prostego modelu liniowego plot(LM$fitted.values,rstandard(LM)). Czy używanie jest plot(reside(lm.base1))wystarczające?
P4: Aby sprawdzić normalność, czy używasz następującego kodu ok?

hist((resid(lm.base1) - mean(resid(lm.base1))) / sd(resid(lm.base1)), freq = FALSE); curve(dnorm, add = TRUE)

P1: Tak - jak każdy model regresji.

Q2: Podobnie jak ogólne modele liniowe, twoja zmienna wyniku nie musi być normalnie dystrybuowana jako zmienna jednowymiarowa. Jednak modele LME zakładają, że reszty modelu są zwykle rozłożone. Tak więc transformacja lub dodanie wag do modelu byłoby sposobem na załatwienie tego (i sprawdzenie za pomocą wykresów diagnostycznych, oczywiście).

P3: plot(myModel.lme)

Q4: qqnorm(myModel.lme, ~ranef(., level=2)). Ten kod pozwoli ci tworzyć wykresy QQ dla każdego poziomu efektów losowych. Modele LME zakładają, że nie tylko reszty wewnątrz klastra są normalnie rozmieszczone, ale także każdy poziom efektów losowych jest również. Zmieniaj levelod 0, 1, do 2, abyś mógł sprawdzić szczury, zadania i resztki wewnątrz badanego.

EDYCJA: Powinienem również dodać, że chociaż zakłada się normalność i że transformacja prawdopodobnie pomaga zmniejszyć problemy z niestandardowymi błędami / efektami losowymi, nie jest jasne, czy wszystkie problemy zostały faktycznie rozwiązane lub czy uprzedzenie nie zostało wprowadzone. Jeśli Twoje dane wymagają transformacji, zachowaj ostrożność przy szacowaniu losowych efektów. Oto artykuł na ten temat .

Wydajesz się dość wprowadzać w błąd co do założeń dotyczących modeli wielopoziomowych. Nie zakłada się jednorodności wariancji w danych, tylko to, że reszty powinny być w przybliżeniu normalnie rozłożone. I predyktory jakościowe są stosowane przez cały czas w regresji (podstawową funkcją w R, która uruchamia ANOVA, jest polecenie regresji liniowej).

Szczegółowe informacje na temat sprawdzania założeń znajdują się w książce Pinheiro i Bates (s. 174, sekcja 4.3.1). Ponadto, jeśli planujesz używać lme4 (którego nie ma w książce), możesz replikować ich wykresy za pomocą wykresu z lmermodelem ( ?plot.merMod).

Aby szybko sprawdzić normalność, wystarczy qqnorm(resid(myModel)).

W odniesieniu do drugiego kwartału:

Według książki Pinheiro i Batesa możesz zastosować następujące podejście:

„ lmeFunkcja umożliwia modelowanie heteroscesdastyczności grupy błędów wewnętrznych za pomocą weightsargumentu. Temat ten zostanie szczegółowo omówiony w § 5.2, ale na razie wystarczy wiedzieć, że varIdentstruktura funkcji wariancji pozwala na różne wariancje dla każdego poziomu współczynnik i może być stosowany do dopasowania modelu heteroscedastycznego [...] ”

Pinheiro i Bates, s. 1 177

Jeśli chcesz sprawdzić równe wariancje sex, możesz użyć tego podejścia:

plot( lm.base2, resid(., type = "p") ~ fitted(.) | sex,
  id = 0.05, adj = -0.3 )

Jeśli wariancje są różne, możesz zaktualizować swój model w następujący sposób:

lm.base2u <- update( lm.base2, weights = varIdent(form = ~ 1 | sex) )
summary(lm.base2u)

Co więcej, możesz spojrzeć na robustlmmopakowanie, które również stosuje metodę ważenia. Rozprawa doktorska Kollera na temat tej koncepcji jest dostępna jako otwarty dostęp („Solidna ocena liniowych modeli mieszanych”). Abstrakt stwierdza:

„Nowa ocena skali, projektowa ocena skali adaptacyjnej, została opracowana w celu zapewnienia solidnej podstawy dla kolejnych solidnych testów. Dokonuje tego poprzez wyrównanie naturalnej heteroskedastyczności reszt i dostosowanie się do solidnego równania szacunkowego dla samej skali Te korekty adaptacyjne projektu są kluczowe w małych próbkach, gdzie liczba obserwacji może być zaledwie pięć razy większa od liczby parametrów do oszacowania lub mniej. ”

Nie mam wystarczającej liczby punktów na komentarze. Widzę jednak konieczność wyjaśnienia niektórych aspektów powyższej odpowiedzi @Johna. Stan Pinheiro i Bates na s. 174:

Założenie 1 - błędy wewnątrz grupy są niezależne i identycznie normalnie rozłożone, ze średnim zerem i wariancją σ2 i są niezależne od efektów losowych.

To stwierdzenie rzeczywiście nie jest jasne na temat jednorodnych wariancji i nie jestem wystarczająco głęboko w statystykach, aby poznać wszystkie matematykę leżącą u podstaw koncepcji LME. Jednak na str. 175, §4.3.1, w części dotyczącej Założenia 1 piszą:

W tej sekcji koncentrujemy się na metodach oceny założenia, że ​​błędy wewnątrz grupy są normalnie rozłożone, są wyśrodkowane na zero i mają stałą wariancję .

Ponadto w poniższych przykładach „ stałe wariancje ” są rzeczywiście ważne. Można więc spekulować, czy implikują jednorodne wariancje, pisząc „ identycznie normalnie rozłożone” na p. 174, nie zajmując się tym bardziej bezpośrednio.

P1: Tak, dlaczego nie?

Q2: Myślę, że wymaganie polega na tym, że błędy są zwykle dystrybuowane.

P3: Można przetestować na przykład testem Levena.