# AIC discordance in mixed effects modeling: which to prioritize?

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

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 ...

• 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? – Pablo Bernabeu Oct 25 '16 at 0:33
• Or should I also try the AIC comparisons through bbmle::AICtab? – Pablo Bernabeu Oct 25 '16 at 0:37
• 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 ... – Ben Bolker Oct 25 '16 at 3:35
• I'll do with the anova comparison. Thanks you very much! – Pablo Bernabeu Oct 25 '16 at 13:10