R help - Oct 2008 - lmer for two models followed by anova to compare the two models

Dear Colleagues,

I run this model:

mod1 <- lmer(x~category+subcomp+category*subcomp+(1|id),data=impchiefsrm)

obtain this summary result:

Linear mixed-effects model fit by REML
Formula: x ~ category + subcomp + category * subcomp + (1 | id)
  Data: impchiefsrm
 AIC  BIC logLik MLdeviance REMLdeviance
4102 4670  -1954       3665         3908
Random effects:
Groups   Name        Variance Std.Dev.
id       (Intercept) 0.11562  0.34003
Residual             0.22765  0.47713
number of obs: 2568, groups: id, 107


run this model (only difference is I've removed the interaction term):

mod2 <- lmer(x~category+subcomp+(1|id),data=impchiefsrm)

obtain this summary result:

Linear mixed-effects model fit by REML
Formula: x ~ category + subcomp + (1 | id)
   Data: impchiefsrm
  AIC  BIC logLik MLdeviance REMLdeviance
 3987 4151  -1966       3823         3931
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.11528  0.33953
 Residual             0.23584  0.48564
number of obs: 2568, groups: id, 107

Note that the loglik from the first model is -1954 and from the second
model loglik is -1966.

Next, to test the difference between the two models I run:

anova(mod1, mod2) and obtain this result:

Data: impchiefsrm
Models:
mod2: x ~ category + subcomp + (1 | id)
mod1: x ~ category + subcomp + category * subcomp + (1 | id)
     Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
mod2 28  3879.1  4042.9 -1911.5
mod1 97  3859.3  4426.9 -1832.7 157.72     69   6.71e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1

Note that in this result the logLiks are reported as -1832.7 and
-1911.5 respectively for models 1 and 2.

Now my question:

Why are the logLiks from the anova command that compares the two
models different from what was reported in the separate model results?

Thanks,

Larry

Lawrence Hanser wrote:
> Dear Colleagues,
> 
> I run this model:
> 
> mod1 <-
lmer(x~category+subcomp+category*subcomp+(1|id),data=impchiefsrm)
> 
> obtain this summary result:
> 
> Linear mixed-effects model fit by REML
> Formula: x ~ category + subcomp + category * subcomp + (1 | id)
>   Data: impchiefsrm
>  AIC  BIC logLik MLdeviance REMLdeviance
> 4102 4670  -1954       3665         3908
> Random effects:
> Groups   Name        Variance Std.Dev.
> id       (Intercept) 0.11562  0.34003
> Residual             0.22765  0.47713
> number of obs: 2568, groups: id, 107
> 
> 
> run this model (only difference is I've removed the interaction term):
> 
> mod2 <- lmer(x~category+subcomp+(1|id),data=impchiefsrm)
> 
> obtain this summary result:
> 
> Linear mixed-effects model fit by REML
> Formula: x ~ category + subcomp + (1 | id)
>    Data: impchiefsrm
>   AIC  BIC logLik MLdeviance REMLdeviance
>  3987 4151  -1966       3823         3931
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  id       (Intercept) 0.11528  0.33953
>  Residual             0.23584  0.48564
> number of obs: 2568, groups: id, 107
> 
> Note that the loglik from the first model is -1954 and from the second
> model loglik is -1966.
> 
> Next, to test the difference between the two models I run:
> 
> anova(mod1, mod2) and obtain this result:
> 
> Data: impchiefsrm
> Models:
> mod2: x ~ category + subcomp + (1 | id)
> mod1: x ~ category + subcomp + category * subcomp + (1 | id)
>      Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
> mod2 28  3879.1  4042.9 -1911.5
> mod1 97  3859.3  4426.9 -1832.7 157.72     69   6.71e-09 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
> 
> Note that in this result the logLiks are reported as -1832.7 and
> -1911.5 respectively for models 1 and 2.
> 
> Now my question:
> 
> Why are the logLiks from the anova command that compares the two
> models different from what was reported in the separate model results?
> 
REML likelihoods cannot be compared between models with different mean 
value structure, so it is switching to ordinary likelihood. This is 
arguably not all that smart, but at least it makes some sense. Compare 
the anova table with half the MLdeviances and you'll see the light.

-- 
    O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
   c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
  (*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907