[R-sig-ME] Model comparison with anova and AIC: p=0 and different AIC-values

[Ben Bolker]

On 13-11-15 09:30 PM, Stefan Th. Gries wrote:
> Hi
> 
> I am trying to do a mixed-effects model analysis on data that were
> published with a repeated-measures ANOVA as follows:
> 
> summary(aov(OVERLAParcsine ~ USED*SAMPLE + Error(NAME/(USED*SAMPLE))))
> 
> In order to first determine the required random-effects structure, I
> created the following two models:
> 
> m.01a.reml <- lmer(OVERLAParcsine ~ USED*poly(SAMPLE, 2) +
> (USED*SAMPLE|NAME), REML=TRUE)
> m.01b.reml <- lmer(OVERLAParcsine ~ USED*poly(SAMPLE, 2) +
> (USED+SAMPLE|NAME), REML=TRUE)
> 
> The problems begin when I try to find out which model accounts for the
> data better:
>> anova(m.01a.reml, m.01b.reml, test="F")
> Data:
> Models:
> m.01b.reml: OVERLAParcsine ~ USED * poly(SAMPLE, 2) + (USED + SAMPLE | NAME)
> m.01a.reml: OVERLAParcsine ~ USED * poly(SAMPLE, 2) + (USED * SAMPLE | NAME)
>            Df     AIC      BIC logLik deviance Chisq Chi Df Pr(>Chisq)
> m.01b.reml 13 -144.74 -115.139 85.368  -170.74
> m.01a.reml 17  -21.20   17.503 27.600   -55.20     0      4          1
> Warning message:
> In optwrap(optimizer, devfun, x at theta, lower = x at lower) :
>   convergence code 1 from bobyqa: bobyqa -- maximum number of function
> evaluations exceeded
> 
> So, there are negative AIC values. Ok, according to
> <http://r.789695.n4.nabble.com/Negative-AIC-td899943.html> and
> <http://emdbolker.wikidot.com/faq> this may not be real problematic so
> I would go with m.01b.reml because its AIC value is smaller. But the
> remaining values are of course also strange, with Chisq=0 because of
> the negative difference of the deviance values.

  I'm not sure about that, it looks like it *might* be a glitch in the
display.  Do the results of anova(m.01b.reml, m.01a.reml) look more
sensible?

> 
> Now, I also used AIC on the models and get results that are different
> from the anova comparison above:
> 
>> AIC(m.01b.reml)
> [1] -120.4052
>> AIC(m.01a.reml)
> [1] 20.96197
> 
> So, two questions:
> 
> (i) which AIC-values are correct - anova or AIC?

see http://glmm.wikidot.com/faq#error_anova_lmer_AIC:

As pointed out by several users (here, here, and here, for example), the
AICs computed for lmer models in summary and anova are different;
summary uses the REML specification as specified when fitting the model,
while anova always uses REML=FALSE (to safeguard users against
incorrectly using restricted MLs to compare models with different fixed
effect components). (This behavior is slightly overzealous since users
might conceivably be using anova to compare models with different random
effects [although this also subject to boundary effects as described
elsewhere in this document …])

  Note that the AIC differences are almost identical (140)



> (ii) so I can't do a p-value based test on which model to use?

 Well, anova() gave you a likelihood ratio test (based on a ML
refitting, which is not necessarily what you want: see
https://github.com/lme4/lme4/issues/141 ).

  I would say that the bottom line here is that since the more complex
model m.01b.reml is about 60 log-likelihood units (or "REML criterion
units") better, it doesn't really matter what test you do -- the more
complex model is overwhelmingly better.

> 
> Thanks,
> STG
> --
> Stefan Th. Gries
> -----------------------------------------------
> University of California, Santa Barbara
> http://www.linguistics.ucsb.edu/faculty/stgries
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>