Wartość P dla terminu interakcji w modelach efektów mieszanych przy użyciu lme4

Ja analizując niektóre dane behawioralne używając lme4w R, głównie po doskonałych tutoriali Bodo zimowy , ale nie rozumiem, jeśli mam właściwie obsługi interakcji. Co gorsza, nikt inny, kto bierze udział w tych badaniach, nie używa mieszanych modeli, więc jestem nieco zaniepokojony, jeśli chodzi o upewnienie się, że wszystko jest w porządku.

Zamiast po prostu poprosić o pomoc, pomyślałem, że powinienem dołożyć wszelkich starań, aby zinterpretować problem, a następnie poprosić o zbiorowe poprawki. Kilka innych oprócz:

• Podczas pisania znalazłem to pytanie , pokazując, że nlmebardziej bezpośrednio podam wartości p dla terminów interakcji, ale myślę, że nadal warto zadawać pytania w odniesieniu do lme4.
• Livius'odpowiedź na to pytanie zawierała linki do wielu dodatkowych lektur, które będę starał się przejść w ciągu najbliższych kilku dni, więc będę komentować wszelkie postępy, które przyniosą.

W moich danych mam zmienną zależną dv, conditionmanipulację (0 = kontrola, 1 = warunek eksperymentalny, co powinno skutkować wyższą dv), a także warunek wstępny, oznaczony appropriate: próby zakodowane 1w tym celu powinny wykazywać efekt, ale próby zakodowane 0mogą nie, ponieważ brakuje kluczowego czynnika.

I przypadkowe również dwa przecięcia, o subject, i target, co odzwierciedla skorelowanych dvwartości w ramach każdego przedmiotu, a w każdej z 14 problemy rozwiązane (każdy uczestnik rozwiązane zarówno na kontrolną i doświadczalną wersji każdej problem).

library(lme4)
data = read.csv("data.csv")

null_model        = lmer(dv ~ (1 | subject) + (1 | target), data = data)
mainfx_model      = lmer(dv ~ condition + appropriate + (1 | subject) + (1 | target),
                         data = data)
interaction_model = lmer(dv ~ condition + appropriate + condition*appropriate +
                              (1 | subject) + (1 | target), data = data)
summary(interaction_model)

Wynik:

## Linear mixed model fit by REML ['lmerMod']
## ...excluded for brevity....
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.006594 0.0812  
##  target   (Intercept) 0.000557 0.0236  
##  Residual             0.210172 0.4584  
## Number of obs: 690, groups: subject, 38; target, 14
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    0.2518     0.0501    5.03
## conditioncontrol               0.0579     0.0588    0.98
## appropriate                   -0.0358     0.0595   -0.60
## conditioncontrol:appropriate  -0.1553     0.0740   -2.10
## 
## Correlation of Fixed Effects:
## ...excluded for brevity.

ANOVA wykazuje wtedy interaction_modelznacznie lepsze dopasowanie niż mainfx_model, z czego dochodzę do wniosku, że występuje znacząca interakcja (p = 0,035).

anova(mainfx_model, interaction_model)

Wynik:

## ...excluded for brevity....
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## mainfx_model       6 913 940   -450      901                          
## interaction_model  7 910 942   -448      896  4.44      1      0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Stamtąd izoluję podzbiór danych, dla których appropriatewymóg jest spełniony (tj. appropriate = 1), I dla tego pasuje model zerowy, a model obejmujący conditionjako efekt, porównaj dwa modele za pomocą ANOVA ponownie, i oto znajdź, że conditionjest znaczącym predyktorem.

good_data = data[data$appropriate == 1, ]
good_null_model   = lmer(dv ~ (1 | subject) + (1 | target), data = good_data)
good_mainfx_model = lmer(dv ~ condition + (1 | subject) + (1 | target), data = good_data)

anova(good_null_model, good_mainfx_model)

Wynik:

## Data: good_data
## models:
## good_null_model: dv ~ (1 | subject) + (1 | target)
## good_mainfx_model: dv ~ condition + (1 | subject) + (1 | target)
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## good_null_model    4 491 507   -241      483                          
## good_mainfx_model  5 487 507   -238      477  5.55      1      0.018 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nie widzę tu zbyt wiele do powiedzenia. Myślę, że wykonałeś dobrą robotę.

Istnieje kilka sposobów omawiania przez ludzi testowania efektów i uzyskiwania wartości p dla skomplikowanych modeli efektów mieszanych. Jest to dobry przegląd tutaj . Najlepiej jest stosować metody intensywnie obliczeniowe (ładowanie początkowe lub metody bayesowskie), ale jest to bardziej zaawansowane dla większości ludzi. Drugą najlepszą (i najwygodniejszą) metodą jest zastosowanie testu współczynnika wiarygodności. To właśnie anova()(technicznie ? Anova.merMod () ) robi. Ważne jest, aby stosować tylko test współczynnika wiarygodności modeli, które były zgodne z pełnym maksymalnym prawdopodobieństwem , zamiast ograniczonego maksymalnego prawdopodobieństwa(REML). Z drugiej strony, do ostatecznego modelu i do interpretacji, chcesz użyć REML. Jest to mylące dla wielu osób. W danych wyjściowych widzimy, że pasujesz do swoich modeli za pomocą REML (dzieje się tak, ponieważ opcja jest TRUEdomyślnie ustawiona na lmer(). To znaczy, że twój test jest nieprawidłowy, ponieważ jest to tak częsty błąd, który anova.merMod()zawiera refitargument, który wartość domyślna jest ustawiona na TRUE, a ty jej nie zmieniłeś. Więc foresight twórców pakietów cię tam ocalił.

Jeśli chodzi o twoją strategię rozpakowywania interakcji, to co zrobiłeś jest w porządku. Pamiętaj tylko, że interakcja wykorzystuje wszystkie dane do swojego testu. Możliwe jest znaczące oddziaływanie, ale żaden z testów warstwowych nie jest znaczący, co dezorientuje niektóre osoby. (Wygląda jednak na to, że ci się nie zdarzyło.)

Sam jestem nowicjuszem i postępuję zgodnie z radami Zuura i innych. Korzystam lmez nlmepakietu zamiast, lme4gdy muszę dodać hierarchiczną strukturę błędów do skądinąd liniowego modelu. Moja odpowiedź może być daleka.

Dwa komentarze:

(1) Nie jestem pewien, czy sensowne jest testowanie tylko conditionw podzbiorze appropriate==1. Jeśli chcesz uzyskać wartości p dla głównych efektów, możesz użyć Anovaz pakietu „samochód”:

require(car)
Anova(M,type="III")# type III Sum of Squares. M was fitted using method=REML

Jeśli chcesz rozwiązać interakcję, możesz uruchomić porównania w parach bezpośrednio (?) Lub zrobić to, co zrobiłeś, ale na obu podzbiorach (tj. Także z podzestawem gdzie appropriate==0).

(2) Możesz najpierw wybrać strukturę błędów zamiast zakładać, że (1 | subject) + (1 | target)jest to najlepsza struktura błędów. Z tego, co napisałeś, rozumiem, że conditionjest to czynnik wewnątrz podmiotu, podczas gdy appropriatejest to albo czynnik między podmiotem, albo czynnik między celem. Możesz dodać nachylenia dla czynników wewnątrz podmiotu i / lub wewnątrz celu, na przykład: dv ~ condition + appropriate + (1+condition | subject) + (1 | target)dodaje losowe nachylenie dla czynnika wewnątrz podmiotu condition. Nie są potrzebne nachylenia dla czynników między podmiotami / celami.

Twoje zdrowie