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

[Stefan Th. Gries]

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

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

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:

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

- 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