R help - Apr 2009 - AICs from lmer different with summary and anova

Dear R Helpers,

I have noticed that when I use lmer to analyse data, the summary function
gives different values for the AIC, BIC and log-likelihood compared with the
anova function.

Here is a sample program

#make some data
set.seed(1);
datx=data.frame(array(runif(720),c(240,3),dimnames=list(NULL,c('x1','x2','y'
))))
id=rep(1:120,2); datx=cbind(id,datx)

#give x1 a slight relation with y (only necessary to make the random effects
non-zero in this artificial example)
datx$x1=(datx$y*0.1)+datx$x1

library(lme4)

#fit the data
fit0=lmer(y~x1+x2+(1|id), data=datx); print(summary(fit0),corr=F)
fit1=lmer(y~x1+x2+(1+x1|id), data=datx); print(summary(fit1),corr=F)

#compare the models
anova(fit0,fit1)


Now, look at the output, below. You can see that the AIC from
"print(summary(fit0))" is 87.34, but the AIC for fit0 in
"anova(fit0,fit1)"
is 73.965. There are similar changes for the values of BIC and logLik.

Am I doing something wrong, here? If not, which are the real AIC and logLik
values for the different models?

Thanks for your help,

Jonathan Williams


Output:-
> fit0=lmer(y~x1+x2+(1|id), data=datx); print(summary(fit0),corr=F)
Linear mixed model fit by REML 
Formula: y ~ x1 + x2 + (1 | id) 
   Data: datx 
   AIC   BIC logLik deviance REMLdev
 87.34 104.7 -38.67    63.96   77.34
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.016314 0.12773 
 Residual             0.062786 0.25057 
Number of obs: 240, groups: id, 120

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.50376    0.05219   9.652
x1           0.08979    0.06614   1.358
x2          -0.06650    0.06056  -1.098
> fit1=lmer(y~x1+x2+(1+x1|id), data=datx); print(summary(fit1),corr=F)
Linear mixed model fit by REML 
Formula: y ~ x1 + x2 + (1 + x1 | id) 
   Data: datx 
   AIC   BIC logLik deviance REMLdev
 90.56 114.9 -38.28    63.18   76.56
Random effects:
 Groups   Name        Variance  Std.Dev. Corr  
 id       (Intercept) 0.0076708 0.087583       
          x1          0.0056777 0.075351 1.000 
 Residual             0.0618464 0.248689       
Number of obs: 240, groups: id, 120

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.50078    0.05092   9.835
x1           0.09236    0.06612   1.397
x2          -0.06515    0.06044  -1.078
> anova(fit0,fit1)
Data: datx
Models:
fit0: y ~ x1 + x2 + (1 | id)
fit1: y ~ x1 + x2 + (1 + x1 | id)
     Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fit0  5  73.965  91.368 -31.982                         
fit1  7  77.181 101.545 -31.590 0.7839      2     0.6757

At 16:22 15/04/2009, Jonathan Williams wrote:
>Dear R Helpers,
>
>I have noticed that when I use lmer to analyse data, the summary function
>gives different values for the AIC, BIC and log-likelihood compared with the
>anova function.
I do not think I have seen a reply to this.

What happens if you fit the models originally using ML rather than REML?

>Here is a sample program
>
>#make some data
>set.seed(1);
>datx=data.frame(array(runif(720),c(240,3),dimnames=list(NULL,c('x1','x2','y'
>))))
>id=rep(1:120,2); datx=cbind(id,datx)
>
>#give x1 a slight relation with y (only necessary to make the random effects
>non-zero in this artificial example)
>datx$x1=(datx$y*0.1)+datx$x1
>
>library(lme4)
>
>#fit the data
>fit0=lmer(y~x1+x2+(1|id), data=datx); print(summary(fit0),corr=F)
>fit1=lmer(y~x1+x2+(1+x1|id), data=datx); print(summary(fit1),corr=F)
>
>#compare the models
>anova(fit0,fit1)
>
>
>Now, look at the output, below. You can see that the AIC from
>"print(summary(fit0))" is 87.34, but the AIC for fit0 in
"anova(fit0,fit1)"
>is 73.965. There are similar changes for the values of BIC and logLik.
>
>Am I doing something wrong, here? If not, which are the real AIC and logLik
>values for the different models?
>
>Thanks for your help,
>
>Jonathan Williams
>
>
>Output:-
>
> > fit0=lmer(y~x1+x2+(1|id), data=datx); print(summary(fit0),corr=F)
>Linear mixed model fit by REML
>Formula: y ~ x1 + x2 + (1 | id)
>    Data: datx
>    AIC   BIC logLik deviance REMLdev
>  87.34 104.7 -38.67    63.96   77.34
>Random effects:
>  Groups   Name        Variance Std.Dev.
>  id       (Intercept) 0.016314 0.12773
>  Residual             0.062786 0.25057
>Number of obs: 240, groups: id, 120
>
>Fixed effects:
>             Estimate Std. Error t value
>(Intercept)  0.50376    0.05219   9.652
>x1           0.08979    0.06614   1.358
>x2          -0.06650    0.06056  -1.098
> > fit1=lmer(y~x1+x2+(1+x1|id), data=datx); print(summary(fit1),corr=F)
>Linear mixed model fit by REML
>Formula: y ~ x1 + x2 + (1 + x1 | id)
>    Data: datx
>    AIC   BIC logLik deviance REMLdev
>  90.56 114.9 -38.28    63.18   76.56
>Random effects:
>  Groups   Name        Variance  Std.Dev. Corr
>  id       (Intercept) 0.0076708 0.087583
>           x1          0.0056777 0.075351 1.000
>  Residual             0.0618464 0.248689
>Number of obs: 240, groups: id, 120
>
>Fixed effects:
>             Estimate Std. Error t value
>(Intercept)  0.50078    0.05092   9.835
>x1           0.09236    0.06612   1.397
>x2          -0.06515    0.06044  -1.078
> > anova(fit0,fit1)
>Data: datx
>Models:
>fit0: y ~ x1 + x2 + (1 | id)
>fit1: y ~ x1 + x2 + (1 + x1 | id)
>      Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
>fit0  5  73.965  91.368 -31.982
>fit1  7  77.181 101.545 -31.590 0.7839      2     0.6757
Michael Dewey
http://www.aghmed.fsnet.co.uk

> Am I doing something wrong, here? If not, which are the real AIC and logLik
> values for the different models?
I don't think it's reasonable to expect that the log-likelihood
computed by different functions be should comparable.  Are the
constant terms included or dropped?

Hadley

-- 
http://had.co.nz/