[R] lmer for two models followed by anova to compare the two models

[Peter Dalgaard]

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