Как сравнить модель без случайных эффектов для модели со случайным эффектом с использованием lme4?

Я могу использовать gls() из пакета nlme для создания mod1 без случайных эффектов. Затем я могу сравнить mod1 с помощью AIC для mod2, созданного с помощью lme(), который включает случайный эффект.

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

Есть ли что-то похожее на gls() для пакета lme4, которое позволило бы мне построить mod3 без случайных эффектов и сравнить его с mod4, построенным с использованием lmer(), который включает случайный эффект?

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

Ответы

Ответ 1

С современными ( > 1.0) версиями lme4 вы можете сделать прямое сравнение между lmer fits и соответствующей моделью lm, , но вам нужно использовать ML --- it трудно найти разумный аналог "критерия REML" для модели без случайных эффектов (потому что это будет связано с линейным преобразованием данных, которые устанавливают все фиксированные эффекты в ноль...)

Вы должны знать, что существуют теоретические проблемы с теоретико-информационными сопоставлениями моделей с компонентами дисперсии и без них: см. http://glmm.wikidot.com/faq для получения дополнительной информации.

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

Я предпочитаю вывод в этом формате (delta-AIC, а не сырые значения AIC):

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

Чтобы протестировать, пусть имитирует данные без какого-либо случайного эффекта (мне пришлось попробовать пару семян случайного числа, чтобы получить пример, где стандартный объект std был оценен как нуль):

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