Czy muszę ustawić refit = FALSE podczas testowania losowych efektów w modelach lmer() z anova()?

Question

Obecnie testuję, czy powinienem uwzględnić pewne efekty losowe w moim lmingu, czy nie. Używam do tego funkcji anova. Moja dotychczasowa procedura polega na dopasowaniu modelu do funkcji wywołania funkcji lmer() z REML=TRUE (opcja domyślna). Następnie zadzwonię pod numer anova() w dwóch modelach, w których jeden z nich zawiera losowy efekt do przetestowania, a drugi nie. Jednak dobrze wiadomo, że funkcja anova() zmienia model za pomocą ML, ale w nowej wersji anova() można uniemożliwić anova(), ustawiając opcję refit=FALSE. W celu zbadania efektów losowych należy ustawić refit=FALSE na moje wezwanie do anova() or not? (Jeśli ustawisz refit=FALSE p-wartości są na ogół niższe. Czy wartości p anty-konserwatywny kiedy ustawić refit=FALSE?)

metoda 1:

mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat) 
    mod1_reml <- lmer(x ~ y + z + (y | w), data=dat) 
    anova(mod0_reml, mod1_reml) 

To spowoduje anova() montażu modeli z ML zamiast REML. (Nowsze wersje funkcji anova() będzie również wydać informacji na ten temat.)

Metoda 2:

mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat) 
    mod1_reml <- lmer(x ~ y + z + (y | w), data=dat) 
    anova(mod0_reml, mod1_reml, refit=FALSE) 

To spowoduje anova() Wykonując swoje obliczenia na oryginalnych modeli, to znaczy z REML=TRUE.

Która z dwóch metod jest poprawna, aby sprawdzić, czy powinienem uwzględnić efekt losowy, czy nie?

Dzięki za wszelką pomoc

Answer 1

Ogólnie powiedziałbym, że byłoby wskazane, aby używać refit=FALSE w tym przypadku, ale chodźmy naprzód i spróbować eksperymentu symulacyjnego.

Najpierw dopasować model bez losowej zboczu do zestawu sleepstudy danych, a następnie symulować dane z tego modelu:

library(lme4) 
mod0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy) 
## also fit the full model for later use 
mod1 <- lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy) 
set.seed(101) 
simdat <- simulate(mod0,1000) 

Teraz Założyć zerowe danych z pełnego i zredukowanego modelu i zapisać dystrybucji wartości p wygenerowane przez anova() zi bez refit=FALSE. Jest to w istocie parametryczny test bootstrap hipotezy zerowej; chcemy sprawdzić, czy ma on odpowiednią charakterystykę (tj. równomierną dystrybucję wartości p).

sumfun <- function(x) { 
    m0 <- refit(mod0,x) 
    m1 <- refit(mod1,x) 
    a_refit <- suppressMessages(anova(m0,m1)["m1","Pr(>Chisq)"]) 
    a_no_refit <- anova(m0,m1,refit=FALSE)["m1","Pr(>Chisq)"] 
    c(refit=a_refit,no_refit=a_no_refit) 
} 

Lubię plyr::laply dla jego wygody, choć można równie dobrze użyć for pętli lub jeden z pozostałych *apply podejść.

library(plyr) 
pdist <- laply(simdat,sumfun,.progress="text") 

library(ggplot2); theme_set(theme_bw()) 
library(reshape2) 
ggplot(melt(pdist),aes(x=value,fill=Var2))+ 
    geom_histogram(aes(y=..density..), 
     alpha=0.5,position="identity",binwidth=0.02)+ 
    geom_hline(yintercept=1,lty=2) 
ggsave("nullhist.png",height=4,width=5) 

Type I Error Rate dla alfa = 0.05:

colMeans(pdist<0.05) 
## refit no_refit 
## 0.021 0.026 

Widać, że w tym przypadku te dwie procedury dają praktycznie tę samą odpowiedź i oba procedury są silnie konserwatywne, dla znanych powodów mających do czynienia z faktem, że wartość null z test hipotezy znajduje się na granicy jego realnej przestrzeni. Dla konkretnego przypadku testowania pojedynczego prostego efektu losowego, zmniejszenie o połowę wartości p daje odpowiednią odpowiedź (patrz Pinheiro i Bates 2000 i inne); to rzeczywiście wydaje się dać tutaj rozsądne odpowiedzi, chociaż nie jest to naprawdę uzasadnione, ponieważ tu spadają dwa efektów losowych parametrów (losowy efekt nachylenia i korelacja pomiędzy nachylenia i przecięcia losowe efekty):

colMeans(pdist/2<0.05) 
## refit no_refit 
## 0.051 0.055 

Inne punkty:

• Możecie być w stanie wykonać ćwiczenia z podobną funkcją PBmodcomp z pakietu pbkrtest.
• Pakiet RLRsim jest przeznaczony właśnie dla szybkiego randomizacji (parameteric bootstrap) prób zerowych hipotez dotyczących efektów losowych względem, ale nie wydaje się, aby pracować w tym nieco bardziej skomplikowanej sytuacji
• patrz odpowiednie GLMM faq section podobnych informacji, w tym argumenty za to, dlaczego nie chcesz testować znaczenia efektów losowych w ogóle ...
• dla dodatkowego kredytu można powtórzyć parametryczne uruchomienia początkowe przy użyciu odchyleń (-2 log prawdopodobieństwo), a nie wartości p jako wynik i sprawdź, czy wyniki są zgodne z mieszaniną pomiędzy chi^2_0 (masa punktowa przy 0) i rozkładem chi^2_n (gdzie n wynosi prawdopodobnie 2, ale ja na pewno nie będzie dla tej geometrii)