Devo impostare refit = FALSE durante il test degli effetti casuali nei modelli lmer() con anova()?

Question

Attualmente sto verificando se dovrei includere determinati effetti casuali nel mio modello di LM o no. Io uso la funzione anova per quello. La mia procedura finora è quella di adattare il modello con una chiamata di funzione a lmer() con REML=TRUE (l'opzione predefinita). Quindi chiamo anova() sui due modelli in cui uno di questi include l'effetto casuale da testare e l'altro no. Tuttavia, è noto che la funzione anova() aggiorna il modello con ML, ma nella nuova versione di anova() è possibile impedire a anova() di farlo impostando l'opzione refit=FALSE. Al fine di verificare gli effetti casuali Devo impostare refit=FALSE nella mia chiamata a anova() or not? (Se faccio a impostare refit=FALSE i valori di p tendono ad essere più bassi. Sono i valori di p anti-conservatore quando ho impostato refit=FALSE?)

metodo 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) 

Questo si tradurrà in anova() refitting modelli con ML invece di REML. (Le versioni più recenti della funzione anova() anche uscita un'informazioni su questo.)

Metodo 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) 

Ciò comporterà anova() dei propri calcoli sui modelli originali, cioè con REML=TRUE.

Quale dei due metodi è corretto per verificare se devo includere un effetto casuale o no?

Grazie per qualsiasi aiuto

Answer 1

In generale direi che sarebbe opportuno utilizzare refit=FALSE in questo caso, ma cerchiamo di andare avanti e cercare un esperimento di simulazione.

Prima montare un modello senza una pendenza casuale per il set di dati sleepstudy, quindi simulare i dati da questo modello:

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) 

Ora rimontare i dati nulli con il pieno e il modello ridotto, e salvare la distribuzione di valori p generati da anova() con e senza refit=FALSE. Questo è essenzialmente un test di boot parametrico dell'ipotesi nulla; vogliamo vedere se ha le caratteristiche appropriate (ad esempio, la distribuzione uniforme dei valori di 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) 
} 

mi piace plyr::laply per la sua comodità, anche se si potrebbe utilizzare la stessa facilità di un ciclo for o di uno degli altri *apply approcci.

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) 

tipo I tasso d'errore per alfa = 0.05:

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

Si può vedere che in questo caso le due procedure danno praticamente la stessa risposta e entrambe le procedure sono fortemente conservative, per ragioni ben note hanno a che fare con il fatto che il valore nullo del test di ipotesi è al limite del suo spazio possibile. Per il caso specifico di testare un singolo effetto casuale semplice, dimezzare il valore p dà una risposta appropriata (vedere Pinheiro e Bates 2000 e altri); questo in realtà sembra dare risposte ragionevoli qui, anche se non è in realtà giustificato perché qui stiamo calando due parametri effetti casuali (l'effetto casuale di pendenza e la correlazione tra la pendenza e intercetta effetti casuali):

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

Altri punti:

• si potrebbe essere in grado di fare un esercizio simile con la funzione PBmodcomp dal pacchetto pbkrtest.
• Il pacchetto RLRsim è progettato proprio per randomizzazione veloce (bootstrap parametrico) test di ipotesi nulle relative termini effetti casuali, ma non sembra funzionare in questo leggermente più complessa situazione
• vedi relativo GLMM faq section per informazioni simili, compresi argomenti per il motivo per cui potresti non voler testare il significato degli effetti casuali ...
• per un credito extra puoi ripetere le esecuzioni bootstrap parametriche usando le deviazioni (-2 log verosimiglianza) piuttosto che i valori p come output e verificare se i risultati sono conformi a una miscela tra uno chi^2_0 (massa punto a 0) e una distribuzione chi^2_n (dove n è probabilmente 2, ma non vorrei essere sicuro per questa geometria)