mixtos - ¿Cómo comparar un modelo sin efectos aleatorios a un modelo con un efecto aleatorio utilizando lme4?

modelos lineales mixtos en r

Question

Puedo usar gls () del paquete nlme para compilar mod1 sin efectos aleatorios. Entonces puedo comparar mod1 usando AIC a mod2 construido usando lme () que incluye un efecto aleatorio.

mod1 = gls(response ~ fixed1 + fixed2, method="REML", data)
mod2 = lme(response ~ fixed1 + fixed2, random = ~1 | random1, method="REML",data)
AIC(mod1,mod2)

¿Hay algo similar a gls () para el paquete lme4 que me permitiría crear mod3 sin efectos aleatorios y compararlo con mod4 creado con lmer () que incluye un efecto aleatorio?

mod3 = ???(response ~ fixed1 + fixed2, REML=T, data)
mod4 = lmer(response ~ fixed1 + fixed2 + (1|random1), REML=T, data)
AIC(mod3,mod4)

Answer 1

Con las versiones modernas (> 1.0) de lme4 puede hacer una comparación directa entre los ajustes de lmer y el modelo de lm correspondiente, pero tiene que usar ML --- es difícil encontrar un análogo sensato del "criterio REML" para una Modelo sin efectos aleatorios (porque implicaría una transformación lineal de los datos que ponen a cero todos los efectos fijos ...)

Debe tener en cuenta que existen problemas teóricos con las comparaciones de la teoría de la información entre modelos con y sin componentes de varianza: consulte las Preguntas frecuentes de GLMM para obtener más información.

library(lme4)
fm1 <- lmer(Reaction~Days+(1|Subject),sleepstudy, REML=FALSE)
fm0 <- lm(Reaction~Days,sleepstudy)
AIC(fm1,fm0)
##     df      AIC
## fm1  4 1802.079
## fm0  3 1906.293

Prefiero la salida en este formato (delta-AIC en lugar de los valores sin procesar de AIC):

bbmle::AICtab(fm1,fm0)
##     dAIC  df
## fm1   0.0 4 
## fm0 104.2 3 

Para realizar la prueba, simulemos datos sin efectos aleatorios (tuve que probar un par de semillas de números aleatorios para obtener un ejemplo en el que la desviación estándar entre los sujetos se estimó en cero):

rr <- simulate(~Days+(1|Subject),
               newparams=list(theta=0,beta=fixef(fm1),
                         sigma=sigma(fm1)),
               newdata=sleepstudy,
               family="gaussian",
               seed=103)[[1]]
ss <- transform(sleepstudy,Reaction=rr)
fm1Z <- update(fm1,data=ss)
VarCorr(fm1Z)
##  Groups   Name        Std.Dev.
##  Subject  (Intercept)  0.000  
##  Residual             29.241
fm0Z <- update(fm0,data=ss)
all.equal(c(logLik(fm0Z)),c(logLik(fm1Z)))  ## TRUE