[R] How can I extract the AIC score from a mixed model object produced using lmer?

[Douglas Bates]

On Dec 19, 2007 9:42 AM, David Hewitt <dhewitt at vims.edu> wrote:
>
>
> David Barron-3 wrote:
> >
> > You can calculate the AIC as follows:
> >
> > (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
> > aic1 <- AIC(logLik(fm1))
> >
> >

> Is AIC() [extractAIC()] "valid" for models with random effects? I noticed
> that the help page for extractAIC() does not list models with random
> effects. I think this boils down to the difference between the likelihoods
> for models with and without random effects, and I don't know. Just
> curious...

The log-likelihood for a linear mixed model is well-defined.  Whether
this makes AIC valid or not depends on how comfortable you are with
the idea of AIC in the first place.  My impression is that the
justification for AIC is not entirely rigorous but I must admit that I
haven't gone back to look at the original literature on it.

To the best of my knowledge and ability the log-likelihood from a
model fit by lmer with method = "ML" is properly defined and
accurately evaluated.  (The default estimation criterion in lmer is
REML and models fit by REML provide a close approximation to the
log-likelihood but not an exact result.  If you really want a
log-likelihood and AIC value you should refit with method = "ML".)
What is later done to the log-likelihood to obtain the AIC value is
more problematic.  In particular, one needs to provide a value for the
number of parameters in the model and that can be tricky.  Recently I
was working with models for data on test scores by students over time.
 There were over 200,000 students.  Under one way of counting
parameters, if I incorporate a random effect for the student this
costs me only 1 parameter, corresponding to the variance component for
that random effect.  However, I am incorporating over 200,000 random
effects to help model the observed responses.  So is the number of
parameters 1 or over 200,000?  I don't know.

Regarding the fact the the extractAIC help page doesn't mention models
with random effects, it can't list all the possible methods because
any package can add methods to a generic function.

>
>
> > On 12/18/07, Peter H Singleton <psingleton at fs.fed.us> wrote:
> >>
> >> I am running a series of candidate mixed models using lmer (package lme4)
> >> and I'd like to be able to compile a list of the AIC scores for those
> >> models so that I can quickly summarize and rank the models by AIC. When I
> >> do logistic regression, I can easily generate this kind of list by
> >> creating
> >> the model objects using glm, and doing:
> >>
> >> > md <- c("md1.lr", "md2.lr", "md3.lr")
> >> > aic <- c(md1.lr$aic, md2.lr$aic, md3.lr$aic)
> >> > aic2 <- cbind(md, aic)
> >>
> >> but when I try to extract the AIC score from the model object produced by
> >> lmer I get:
> >>
> >> > md1.lme$aic
> >> NULL
> >> Warning message:
> >> In md1.lme$aic : $ operator not defined for this S4 class, returning NULL
> >>
> >> So... How do I query the AIC value out of a mixed model object created by
> >> lmer?
> >>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
> -----
> David Hewitt
> Virginia Institute of Marine Science
> http://www.vims.edu/fish/students/dhewitt/
> --
> View this message in context: http://www.nabble.com/How-can-I-extract-the-AIC-score-from-a-mixed-model-object-produced-using-lmer--tp14406832p14419438.html
> Sent from the R help mailing list archive at Nabble.com.
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>