[R-sig-ME] AIC in Cross-Classified Model
Doran, Harold
HDoran at air.org
Tue Feb 24 16:01:36 CET 2009
Kyle:
You need to use REML=FALSE in your code to compare models with different fixed effects. The ANOVA is giving the loglik for the ML. But, your using REML.
The AIC is -2*LogLik + 2*number of parameters. In your m0 model you have 1 fixed effect, the residual variance, and I think the relative standard deviations are both scalars here, giving a total of 4 parameters. So, the AIC in the model summary is computed as -2*-1161 + 2*4= 2330 (although, the AIC is 2329 in the model summary)
Now, in the ANOVA output, the loglik for the same model is -1158.9, so -2*-1158.9+2*4 = 2325.8. The ANOVA gives an AIC of 2325.9. Both seem to be using the same number of parameters, but clearly not using the same loglik. Look at the following example to see why.
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(fm2 <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy))
anova(fm1,fm2)
### Compare loglik in ANOVA to loglik in model summary for fm2
### notice they are different
(fm2 <- lmer(Reaction ~ 1 + (Days|Subject), sleepstudy, REML=FALSE))
anova(fm1,fm2)
### Compare loglik in ANOVA to loglik in model summary for fm2
### Now they are the same
HTH,
Harold
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org
> [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
> Of Roberts, Kyle
> Sent: Tuesday, February 24, 2009 9:29 AM
> To: 'Douglas Bates'; r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] AIC in Cross-Classified Model
>
> Dear All,
>
> I was working an example for my class and came upon a
> problem. This example is the cross-classified model presented
> by Hox on p 127 of "Multilevel Analysis." The problem is that
> the AIC and BIC in the summary of the model are not the same
> AIC and BIC when the anova function is used. Any idea as to
> why? I think that the anova function gives the correct AIC,
> but I don't know why. I threw the dataset up on my website if
> anyone wants to replicate.
>
> > sessionInfo()
> R version 2.7.2 (2008-08-25)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> States.1252;LC_MONETARY=English_United
> States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods
> base
>
> other attached packages:
> [1] foreign_0.8-29 mlmRev_0.99875-1 lme4_0.999375-27
> Matrix_0.999375-16
> [5] lattice_0.17-15
>
> loaded via a namespace (and not attached):
> [1] grid_2.7.2 tools_2.7.2
>
> >
> pupils<-read.table("http://www.jkyleroberts.com/rfiles/mlm/HoxData/pup
> > ils.txt", header=T)
> >
> > m0 <- lmer(ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL), pupils)
> > summary(m0)
> Linear mixed model fit by REML
> Formula: ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL)
> Data: pupils
> AIC BIC logLik deviance REMLdev
> 2329 2349 -1161 2318 2321
> Random effects:
> Groups Name Variance Std.Dev.
> PSCHOOL (Intercept) 0.171900 0.41461
> SSCHOOL (Intercept) 0.066652 0.25817
> Residual 0.513128 0.71633
> Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 6.34861 0.07891 80.45
> >
> > m1 <- lmer(ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 |
> SSCHOOL),
> > pupils)
> > summary(m1)
> Linear mixed model fit by REML
> Formula: ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | SSCHOOL)
> Data: pupils
> AIC BIC logLik deviance REMLdev
> 2270 2299 -1129 2243 2258
> Random effects:
> Groups Name Variance Std.Dev.
> PSCHOOL (Intercept) 0.171530 0.41416
> SSCHOOL (Intercept) 0.064809 0.25458
> Residual 0.475236 0.68937
> Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 5.75548 0.10576 54.42
> PUPSEXgirl 0.26132 0.04569 5.72
> PUPSES 0.11407 0.01612 7.08
>
> Correlation of Fixed Effects:
> (Intr) PUPSEX
> PUPSEXgirl -0.254
> PUPSES -0.641 0.075
> >
> > m2 <- lmer(ACHIEV~PUPSEX+PUPSES+PDENOM+SDENOM+(1 | PSCHOOL) + (1 |
> > SSCHOOL), pupils)
> > summary(m2)
> Linear mixed model fit by REML
> Formula: ACHIEV ~ PUPSEX + PUPSES + PDENOM + SDENOM + (1 |
> PSCHOOL) + (1 | SSCHOOL)
> Data: pupils
> AIC BIC logLik deviance REMLdev
> 2273 2312 -1128 2238 2257
> Random effects:
> Groups Name Variance Std.Dev.
> PSCHOOL (Intercept) 0.165721 0.40709
> SSCHOOL (Intercept) 0.058446 0.24176
> Residual 0.475207 0.68935
> Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 5.51888 0.14275 38.66
> PUPSEXgirl 0.26309 0.04568 5.76
> PUPSES 0.11352 0.01612 7.04
> PDENOMyes 0.20400 0.12623 1.62
> SDENOMyes 0.17572 0.09625 1.83
>
> Correlation of Fixed Effects:
> (Intr) PUPSEX PUPSES PDENOM
> PUPSEXgirl -0.200
> PUPSES -0.466 0.075
> PDENOMyes -0.526 0.021 0.004
> SDENOMyes -0.429 0.004 -0.025 -0.014
> >
> > anova(m0, m1, m2)
> Data: pupils
> Models:
> m0: ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL)
> m1: ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | SSCHOOL)
> m2: ACHIEV ~ PUPSEX + PUPSES + PDENOM + SDENOM + (1 | PSCHOOL) +
> m2: (1 | SSCHOOL)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> m0 4 2325.9 2345.5 -1158.9
> m1 6 2255.5 2284.9 -1121.7 74.3678 2 < 2e-16 ***
> m2 8 2253.5 2292.8 -1118.8 5.9684 2 0.05058 .
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> *********************************************************
> Dr. J. Kyle Roberts
> Department of Teaching and Learning
> Annette Caldwell Simmons School of Education
> and Human Development
> Southern Methodist University
> P.O. Box 750381
> Dallas, TX 75275
> 214-768-4494
> http://www.hlm-online.com/
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list