Jak wybrać strukturę efektów losowych i stałych w liniowych modelach mieszanych?

Rozważ następujące dane z dwustronnego projektowania przedmiotów:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

Chciałbym to przeanalizować za pomocą modeli mieszanych liniowych. Biorąc pod uwagę wszystkie możliwe efekty stałe i losowe, istnieje wiele możliwych modeli:

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

1. Jaki jest zalecany sposób wyboru najlepiej dopasowanego modelu w tym kontekście? Jaka jest zalecana procedura podczas korzystania z testów współczynnika wiarygodności? Generowanie modeli w górę (od modelu zerowego do najbardziej złożonego modelu) czy w dół (od najbardziej złożonego modelu do modelu zerowego)? Stopniowe włączenie lub wyłączenie? Czy też zaleca się poddanie wszystkich modeli jednemu testowi współczynnika wiarygodności i wybranie modelu o najniższej wartości p? Jak porównać modele, które nie są zagnieżdżone?

2. Czy zaleca się najpierw znaleźć odpowiednią strukturę efektów stałych, a następnie odpowiednią strukturę efektów losowych lub odwrotnie (znalazłem odniesienia do obu opcji ...)?

3. Jaki jest zalecany sposób raportowania wyników? Zgłaszanie wartości p z testu logarytmu ilorazu wiarygodności porównującego pełny model mieszany (z danym efektem) z modelem zredukowanym (bez danego efektu). Czy może lepiej jest użyć testu współczynnika logarytmu, aby znaleźć najlepszy model dopasowania, a następnie użyć lmerTest, aby zgłosić wartości p z efektów w modelu najlepiej dopasowanym?

Nie jestem pewien, czy naprawdę jest na to kanoniczna odpowiedź, ale dam temu szansę.

Jaki jest zalecany sposób wyboru najlepiej dopasowanego modelu w tym kontekście? Jaka jest zalecana procedura podczas korzystania z testów współczynnika wiarygodności? Generowanie modeli w górę (od modelu zerowego do najbardziej złożonego modelu) czy w dół (od najbardziej złożonego modelu do modelu zerowego)? Stopniowe włączenie lub wyłączenie? Czy też zaleca się poddanie wszystkich modeli jednemu testowi współczynnika wiarygodności i wybranie modelu o najniższej wartości p? Jak porównać modele, które nie są zagnieżdżone?

To zależy od twoich celów.

• W ogóle należy być bardzo , bardzo ostrożny przy wyborze modelu (patrz np. Ta odpowiedź lub ten post , lub po prostu Google „Harrell krok po kroku” ...).
• Jeśli jesteś zainteresowany, aby twoje wartości p były znaczące (tj. Wykonujesz test potwierdzający hipotezę), powinieneś nie robić wybór modelu. Jednak : nie jest dla mnie tak jasne, czy procedury wyboru modelu są równie złe, jeśli dokonuje się wyboru modelu na nieogniskowych częściach modelu , np. Dokonując wyboru modelu na efektach losowych, jeśli głównym zainteresowaniem jest wnioskowanie o efektach stałych.
• Nie ma czegoś takiego jak „umieszczenie wszystkich modeli w jednym teście współczynnika wiarygodności” - testowanie współczynnika wiarygodności jest procedurą parą. Jeśli chciałbyś dokonać wyboru modelu (np.) Na efektach losowych, prawdopodobnie zaleciłbym podejście „od razu” przy użyciu kryteriów informacyjnych, jak w tym przykładzie - co najmniej pozwala uniknąć niektórych problemów z podejściem etapowym (ale nie wybór modelu bardziej ogólnie).
• Barr i in. 2013 „Keep it maximal” Journal of Memory and Language (doi: 10.1016 / j.jml.2012.11.001) zaleca stosowanie modelu maksymalnego (tylko).
• Shravan Vasishth nie zgadza się , argumentując, że takie modele będą słabe, a przez to problematyczne, chyba że zestaw danych jest bardzo duży (a stosunek sygnału do szumu jest wysoki)
• Innym rozsądnie uzasadnionym podejściem jest dopasowanie dużego ale rozsądnego modelu, a następnie, jeśli dopasowanie jest pojedyncze, usuń warunki, dopóki nie będzie już więcej
• Z pewnymi zastrzeżeniami (wymienionymi w GLMM FAQ ), możesz użyć kryteriów informacyjnych, aby porównać nie zagnieżdżone modele z różnymi losowymi efektami (chociaż Brian Ripley się nie zgadza: patrz dół str. 6 tutaj )

Czy zaleca się najpierw znaleźć odpowiednią strukturę efektów stałych, a następnie odpowiednią strukturę efektów losowych lub odwrotnie (znalazłem odniesienia do obu opcji ...)?

Nie sądzę, żeby ktokolwiek wiedział. Zobacz poprzednią odpowiedź na temat wyboru modelu bardziej ogólnie. Jeśli potrafisz wystarczająco jasno zdefiniować swoje cele (co robi niewiele osób), pytanie może być możliwe. Jeśli masz odniesienia do obu opcji, przydałaby się edycja pytania, aby je uwzględnić ... ten przykład (cytowany powyżej) wykorzystuje kryteria informacyjne do wyboru części efektów losowych, a następnie unika wyboru na część modelu z efektem stałym.

Jaki jest zalecany sposób raportowania wyników? Zgłaszanie wartości p z testu logarytmu ilorazu wiarygodności porównującego pełny model mieszany (z danym efektem) z modelem zredukowanym (bez danego efektu). Czy może lepiej jest użyć testu współczynnika logarytmu, aby znaleźć najlepszy model dopasowania, a następnie użyć lmerTest, aby zgłosić wartości p z efektów w modelu najlepiej dopasowanym?

To (niestety) kolejne trudne pytanie. Jeśli zgłaszasz efekty krańcowe zgłoszone przez lmerTest, musisz martwić się o marginesowość (np. Czy szacunki głównych efektów Ai Bsą znaczące, gdy w modelu występuje interakcja A-przez- B); jest to ogromna puszka robaków, ale jest nieco złagodzona, jeśli użyjesz jej contrasts="sum"zgodnie z zaleceniami afex::mixed(). Trochę też pomagają zrównoważone projekty. Jeśli naprawdę chcesz zapisać wszystkie te pęknięcia, myślę, że poleciłbym afex::mixed, co daje wyniki podobne do lmerTest, ale próbuje rozwiązać te problemy.

Aktualizacja maja 2017 r . : Jak się okazuje, jest wiele tego, co tu napisałem trochę niesłuszne . Niektóre aktualizacje są wprowadzane w całym poście.

Zgadzam się bardzo z tym, co już powiedział Ben Bolker (dziękuję za wiadomość afex::mixed()), ale pozwólcie, że dodam jeszcze kilka ogólnych i szczegółowych przemyśleń na ten temat.

Skoncentruj się na efektach stałych i losowych oraz na sposobie raportowania wyników

W przypadku rodzaju badań eksperymentalnych, które są reprezentowane w przykładowym zbiorze danych od Jonathana Barona, ważne jest pytanie, czy zmanipulowane czy nie czynnik ma ogólny efekt. Na przykład, czy znajdujemy ogólny główny efekt lub interakcję Task? Ważną kwestią jest to, że w tych zbiorach danych zwykle wszystkie czynniki są pod pełną kontrolą eksperymentalną i losowo przypisywane. W związku z tym zainteresowanie skupia się zwykle na ustalonych efektach.
Natomiast składowe efektów losowych można postrzegać jako parametry „uciążliwe”, które wychwytują wariancję systemową (tj. Różnice międzyosobnicze w wielkości efektu), które niekoniecznie są istotne dla pytania głównego. Z tego punktu widzenia sugerowana przez Barr i in. następuje nieco naturalnie. Łatwo sobie wyobrazić, że eksperymentalna manipulacja nie dotyczy wszystkich osób w dokładnie taki sam sposób i chcemy to kontrolować. Z drugiej strony liczba czynników lub poziomów zwykle nie jest zbyt duża, więc niebezpieczeństwo nadmiernego dopasowania wydaje się stosunkowo niewielkie.

W związku z tym podążałbym za sugestią Barra i in. i podaj maksymalną strukturę efektów losowych i zgłoś testy stałych efektów jako moje główne wyniki. Aby przetestować ustalone efekty, sugerowałbym również użycie, afex::mixed()ponieważ raportuje testy efektów lub czynników (zamiast testu parametrów) i oblicza te testy w dość rozsądny sposób (np. Używa tej samej struktury efektów losowych dla wszystkich modeli, w których pojedynczy efekt jest usuwany, używa kontrastów suma do zera, oferuje różne metody obliczania wartości p , ...).

Co z przykładowymi danymi

Problem z podanymi danymi przykładowymi polega na tym, że dla tego zestawu danych maksymalna struktura efektów losowych prowadzi do przesyconego modelu, ponieważ istnieje tylko jeden punkt danych na komórkę projektu:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

W rezultacie lmerdławiki w strukturze maksymalnych efektów losowych:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

Niestety, według mojej wiedzy nie ma uzgodnionego sposobu rozwiązania tego problemu. Ale pozwól mi naszkicować i przedyskutować niektóre:

1. Pierwszym rozwiązaniem może być usunięcie najwyższego losowego nachylenia i przetestowanie efektów dla tego modelu:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67

   Jednak to rozwiązanie jest trochę ad hoc i nie jest nadmiernie motywowane.

   Aktualizacja maja 2017 r .: To podejście, które obecnie popieram. Zobacz ten post na blogu i szkic rozdziału, którego jestem współautorem , sekcja „Struktury efektów losowych dla tradycyjnych projektów ANOVA”.

2. Alternatywnym rozwiązaniem (i takim, które może być popierane przez dyskusję Barr i wsp.) Może być zawsze usunięcie losowych nachyleń w celu uzyskania najmniejszego efektu. Wiąże się to jednak z dwoma problemami: (1) Jakiej struktury efektów losowych używamy, aby dowiedzieć się, jaki jest najmniejszy efekt, i (2) R niechętnie usuwa efekt niższego rzędu, taki jak efekt główny, jeśli efekty wyższego rzędu, takie jak interakcja tego efektu jest obecna (patrz tutaj ). W rezultacie należałoby ręcznie ustawić tę strukturę efektów losowych i przekazać tak skonstruowaną matrycę modelu do wywołania Lmer.

3. Trzecim rozwiązaniem może być zastosowanie alternatywnej parametryzacji części efektów losowych, a mianowicie takiej, która odpowiada modelowi RM-ANOVA dla tych danych. Niestety (?) lmerNie zezwala na „ujemne wariancje”, więc ta parametryzacja nie odpowiada dokładnie RM-ANOVA dla wszystkich zestawów danych , patrz dyskusja tutaj i gdzie indziej (np. Tutaj i tutaj ). „Lmer-ANOVA” dla tych danych to:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76

Biorąc pod uwagę wszystkie te problemy, po prostu nie użyłbym lmerdo dopasowania zestawów danych, dla których istnieje tylko jeden punkt danych na komórkę projektu, chyba że dostępne jest bardziej uzgodnione rozwiązanie problemu maksymalnej struktury losowych efektów.

4. Zamiast tego chciałbym, aby nadal można było użyć klasycznej ANOVA. Użycie jednego z opakowań car::Anova()w afexwynikach byłoby:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

   Zauważ, że afexteraz pozwala również zwrócić model wyposażony w aovktóry można przekazać do lsmeanstestów post hoc (ale do testu efektów zgłoszone przez car::Anovasą jeszcze bardziej rozsądne):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95