0
$\begingroup$

I'd like to ask about mixed effects models. I'm modeling based on the significance of the likelihood ratio, within R. As a primary criterion I use the p-significance of the model comparison (anova()). Yet, with any significant comparison, I'll also check the AIC() difference separately. As I understand it, this last check helps to avoid overfitting. For a couple of times already, I have encountered the infamous discordance between the two, with a model comparison that is significant and then a general AIC difference that falls below 2 points (see image attached). 

So, I wonder how come that the model comparison doesn't pick up the overfitting... That is, do I definitely have to lighten up the model?

Thanks a lot for any advice

Output from the anova comparison, and from AIC model comparison

$\endgroup$
1
$\begingroup$

There are both statistical and computational issues involved here (I would not recommend doing both p-value and AIC comparisons - stick to one or the other!).

The first thing to do is to refit the original models with REML=FALSE; I'm pretty certain that the AIC comparison does not prevent you doing the silly thing of comparing the 'log-likelihoods' (more properly 'REML criteria', as what you get when REML=TRUE is not a log-likelihood), as anova() does ...

We really should have considered making logLik() return NA for REML-fitted models to help people avoid this trap ...

Here's an example.

Fit full and reduced models with ML and REML:

library(lme4)
fm1_REML <- lmer(Reaction~Days+(Days|Subject),sleepstudy)
fm0_REML <- update(fm1_REML,. ~ . - Days)
fm1_ML <- update(fm1_REML,REML=FALSE)
fm0_ML <- update(fm0_REML,REML=FALSE)

LRTs:

anova(fm1_ML,fm0_ML)  ## p-value: 1.226e-06
anova(fm1_REML,fm0_REML)  ## 'refitting models with ML': same p-value

AIC comparisons (I like bbmle::AICtab for readability)

library(bbmle)
AICtab(fm1_ML,fm0_ML)  ## delta-AIC=21.5
AICtab(fm1_REML,fm0_REML)  ## delta-AIC=24.2

Here's how you would override the anova() safeguards to do a bogus likelihood ratio test on the REML-fitted models ...

pchisq(2*(logLik(fm1_REML)-logLik(fm0_REML)),
       df=1,lower.tail=FALSE)  ## 3.05e-07

This p-value is an order of magnitude smaller than the correct LRT, which matches up well with the larger delta-AIC in the bogus AIC comparison based on the REML-fitted models ...

$\endgroup$
  • $\begingroup$ Thank you so much. I really appreciate the elaboration, even though I'm too new to grasp all of it. Could I then stick to the anova comparison, because it is indeed sensitive to a general overfitting of the model? $\endgroup$ – Pablo Bernabeu Oct 25 '16 at 0:33
  • $\begingroup$ Or should I also try the AIC comparisons through bbmle::AICtab? $\endgroup$ – Pablo Bernabeu Oct 25 '16 at 0:37
  • $\begingroup$ If you want to do the AIC comparisons, make sure all the models you're fitting are fitted with REML=FALSE. If that doesn't solve your immediate problem, update/edit your question accordingly and I'll have a look ... $\endgroup$ – Ben Bolker Oct 25 '16 at 3:35
  • $\begingroup$ I'll do with the anova comparison. Thanks you very much! $\endgroup$ – Pablo Bernabeu Oct 25 '16 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.