[R-sig-ME] Model comparison with anova and AIC: p=0 and different AIC-values
Ben Bolker
bbolker at gmail.com
Sat Nov 16 19:10:10 CET 2013
On 13-11-16 12:52 PM, Stefan Th. Gries wrote:
> Hi
>
> BB> 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?
> The results of the two anova functions are the same (and come with warnings):
>
>> anova(m.01a.reml, m.01b.reml) # one order
> 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
>
>> anova(m.01b.reml, m.01a.reml) # other order
> 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
>
The warnings are unrelated to the anova display; they come from the
ML-refitting done by anova.merMod.
> BB> 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
> Ah ok, I did not know that. I do not get any AIC-values though from
> the summary output:
>
>> summary(m.01a.reml) # one model
> Linear mixed model fit by REML ['lmerMod']
> Formula: OVERLAParcsine ~ USED * poly(SAMPLE, 2) + (USED * SAMPLE | NAME)
>
> REML criterion at convergence: -13.038
>
> Random effects:
> Groups Name Variance Std.Dev. Corr
> NAME (Intercept) 0.011097 0.10534
> USEDyes 0.026935 0.16412 0.56
> SAMPLE 0.000392 0.01980 0.80 0.17
> USEDyes:SAMPLE 0.005344 0.07310 -0.80 -0.17 -1.00
> Residual 0.003081 0.05551
> Number of obs: 72, groups: NAME, 12
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 0.11711 1.54844 0.076
> USEDyes 0.33570 5.61965 0.060
> poly(SAMPLE, 2)1 0.41849 8.24332 0.051
> poly(SAMPLE, 2)2 -0.03079 0.07850 -0.392
> USEDyes:poly(SAMPLE, 2)1 0.55875 30.43556 0.018
> USEDyes:poly(SAMPLE, 2)2 -0.25532 0.11102 -2.300
>
> Correlation of Fixed Effects:
> (Intr) USEDys p(SAMPLE,2)1 p(SAMPLE,2)2 USED:(SAMPLE,2)1
> USEDyes -1.000
> p(SAMPLE,2)1 1.000 -1.000
> p(SAMPLE,2)2 0.000 0.000 0.000
> USED:(SAMPLE,2)1 -1.000 1.000 -1.000 0.000
> USED:(SAMPLE,2)2 0.000 0.000 0.000 -0.707 0.000
>
>> summary(m.01b.reml) # other model
> Linear mixed model fit by REML ['lmerMod']
> Formula: OVERLAParcsine ~ USED * poly(SAMPLE, 2) + (USED + SAMPLE | NAME)
>
> REML criterion at convergence: -146.4052
>
> Random effects:
> Groups Name Variance Std.Dev. Corr
> NAME (Intercept) 1.500e-28 1.225e-14
> USEDyes 8.008e-02 2.830e-01 0.42
> SAMPLE 3.362e-09 5.798e-05 0.42 1.00
> Residual 2.777e-03 5.270e-02
> Number of obs: 72, groups: NAME, 12
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 0.117109 0.009852 11.887
> USEDyes 0.335704 0.082629 4.063
> poly(SAMPLE, 2)1 0.418488 0.078339 5.342
> poly(SAMPLE, 2)2 -0.030790 0.074527 -0.413
> USEDyes:poly(SAMPLE, 2)1 0.558745 0.105396 5.301
> USEDyes:poly(SAMPLE, 2)2 -0.255317 0.105396 -2.422
>
> Correlation of Fixed Effects:
> (Intr) USEDys p(SAMPLE,2)1 p(SAMPLE,2)2 USED:(SAMPLE,2)1
> USEDyes 0.353
> p(SAMPLE,2)1 0.140 0.305
> p(SAMPLE,2)2 0.000 0.000 0.000
> USED:(SAMPLE,2)1 0.000 0.000 -0.673 0.000
> USED:(SAMPLE,2)2 0.000 0.000 0.000 -0.707 0.000
Yes, but [I think] AIC() gives you the 'AIC' based on the REML
criterion: -2*(REML criterion)+2*(# parameters) rather than (we can
argue about what it really means)
fm1 <- lmer(Reaction~Days+(Days|Subject),sleepstudy,REML=TRUE)
fm2 <- update(fm1,.~Days+(1|Subject),REML=TRUE)
fm3 <- update(fm1,REML=FALSE)
fm4 <- update(fm2,REML=FALSE)
AIC(fm1,fm2,fm3,fm4)
df AIC
fm1 6 1755.628
fm2 4 1794.465
fm3 6 1763.939
fm4 4 1802.079
I confirm that anova(fm1,fm2); anova(fm2,fm1); anova(fm3,fm4); and
anova(fm4,fm3) all give results that differ only in their labels (all
the numerical values are the same).
> BB> (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 …])
> Yes, that's what I was trying to do, following some advice in the
> literature (I think the Zuur et al. books)
>
> BB> Note that the AIC differences are almost identical (140)
> True, the anova one is 123, the separate AIC ones are 140.
>
> BB> 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 ).
> Ok, I get that now, thanks.
>
> BB> 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.
> Ok, just to clarify: m.01b.reml is the more complex one:
Yes, that's what I said ("the more complex model m.01b.reml is ... 60
log-likelihood units ... better")
>
> 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)
>
> but so you're saying I should go with m.01a.reml because ... do you
> get to the 60 LL units because of
No, I'm saying you should go with m.01b.reml, because it fits much
better. The REML criterion (analogous to twice the negative
log-likelihood [or deviance], essentially the "negative restricted
maximum likelihood") is smaller. The log-likelihood is larger.
>
> - the REML criterion at convergence: (-13.038 for m.01a.reml and
> -146.4052 for m.01b.reml), or
> - because of the logLik diff in the anova output?
>
> Thanks so much,
> STG
> --
> Stefan Th. Gries
> -----------------------------------------------
> University of California, Santa Barbara
> http://www.linguistics.ucsb.edu/faculty/stgries
> -----------------------------------------------
>
More information about the R-sig-mixed-models
mailing list