20

I can use gls() from the nlme package to build mod1 with no random effects. I can then compare mod1 using AIC to mod2 built using lme() which does include a random effect.

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

Is there something similar to gls() for the lme4 package which would allow me to build mod3 with no random effects and compare it to mod4 built using lmer() which does include a random effect?

mod3 = ???(response ~ fixed1 + fixed2, REML=T, data)
mod4 = lmer(response ~ fixed1 + fixed2 + (1|random1), REML=T, data)
AIC(mod3,mod4)
Ben Bolker
  • 173,430
  • 21
  • 312
  • 389
It Figures
  • 363
  • 1
  • 2
  • 11

1 Answers1

30

With modern (>1.0) versions of lme4 you can make a direct comparison between lmer fits and the corresponding lm model, but you have to use ML --- it's hard to come up with a sensible analogue of the "REML criterion" for a model without random effects (because it would involve a linear transformation of the data that set all of the fixed effects to zero ...)

You should be aware that there are theoretical issues with information-theoretic comparisons between models with and without variance components: see the GLMM FAQ for more information.

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

I prefer output in this format (delta-AIC rather than raw AIC values):

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

To test, let's simulate data with no random effect (I had to try a couple of random-number seeds to get an example where the among-subject std dev was actually estimated as zero):

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
Ben Bolker
  • 173,430
  • 21
  • 312
  • 389
  • 1
    with the example I gave, or with your own data? I can run the example fine with the current (devel) version of lme4. If it's with your own data, then more information is required; either ask a new question on StackOverflow, or send an e-mail to `r-sig-mixed-models@r-project.org` [subscribe to the list first; you can find the info/subscription page by searching], with a reproducible example in either case (and the results of `packageVersion("lme4")`) – Ben Bolker Jul 19 '18 at 15:39
  • And what can be the approach when the random model refitted with ML results to be singular? But it is not singular when fitted with REML...any suggestion? – Asier Mar 21 '20 at 12:30
  • Not really surprising, especially with a highly parameterized/unstable model. Can you post a *reproducible* example to [CrossValidated](https://stats.stackexchange.com) or `r-sig-mixed-models@r-project.org` ? – Ben Bolker Mar 21 '20 at 17:21
  • Thanks a lot, Ben. I would try to do that... I agree with you, I think the the model is too highly parameterized. Odd enough, it still has much lower AICc than other models were one of the factors have been modelled as linear and quadratic covariates. That means that I will need to fit the whole dataset in the reproducible example. – Asier Mar 23 '20 at 21:42
  • Done, I provided a reproducible example at Stack Overflow: stackoverflow.com/q/60892398/13099627?sem=2 – Asier – Asier Mar 30 '20 at 12:30