Model comparison with anova and AIC: p=0 and different AIC-values
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