[R] Function to Sort and test AIC for mixed model lme?

[Douglas Bates]

On 5/24/07, Ken Nussear <knussear at usgs.gov> wrote:
> > Ken Nussear <knussear <at> mac.com> writes:
> >
> > > I'm running a series of mixed models using lme, and I wonder if
> > there
> > > is a way to sort them by AIC prior to testing using anova
> > > (lme1,lme2,lme3,....lme7) other than by hand.
> >
> > You can try something like the following. However, also consider
> > using dropterm or stepAIC in MASS.
> >
> > Dieter
> >
> > #---------------------
> >
> > library(nlme)
> > fmlist = vector("list",2)
> > fmlist[[1]] = lme(distance ~ age, data = Orthodont,method="ML")
> > fmlist[[2]] = lme(distance ~ age + Sex, data = Orthodont,   random
> > = ~ 1,method="ML")
> > aic0 = unlist(lapply(fmlist,AIC))
> > aic0 # larger first
> > fmlist1 = fmlist[order(aic0)]
> > unlist(lapply(fmlist1,AIC)) # smaller first
>
> I looked at stepAIC, but wanted to specify my own models.
>
> I have come pretty close with this
>
> # grab all lme objects
> tst1 <- ls(pat=".ml")
>  > tst1
> [1] "lme.T972way.ml"  "lme.T97FULL.ml"  "lme.T97NOINT.ml"
> "lme.T97NULL.ml"  "lme.T97fc.ml"    "lme.T97min.ml"
> [7] "lme.T97ns.ml"
>
> #create array of AIC for all in tst1
> tst2 <- lapply(lapply(tst1,get),AIC)
>  > tst2
> [[1]]
> [1] 507.0991
>
> [[2]]
> [1] 508.7594
>
> [[3]]
> [1] 564.8574
>
> [[4]]
> [1] 624.3053
>
> [[5]]
> [1] 502.5878
>
> [[6]]
> [1] 569.8188
>
> [[7]]
> [1] 504.8971
>
> #sort tst1 by order of tst2
> tst3 <- tst1[order(unlist(tst2))]
>
>  > tst3
> [1] "lme.T97fc.ml"    "lme.T97ns.ml"    "lme.T972way.ml"
> "lme.T97FULL.ml"  "lme.T97NOINT.ml" "lme.T97min.ml"
> [7] "lme.T97NULL.ml"
>
>
> The problem comes with getting the final anova statement to run.
>
>  >anova(tst3)
> Error in anova(tst3) : no applicable method for "anova"
>
> I also tried
>
> tst4 <- paste(toString(tst3),collapse="")
>
>  > tst4
> [1] "lme.T97fc.ml, lme.T97ns.ml, lme.T972way.ml, lme.T97FULL.ml,
> lme.T97NOINT.ml, lme.T97min.ml, lme.T97NULL.ml"
>  >
>  > anova(tst4)
> Error in anova(tst4) : no applicable method for "anova"
>
> Any ideas on getting the last part to work?


tst3 is a character vector of names so you would need to use

do.call(anova, lapply(tst3, get))

I think it is better to create a list of fitted models initially
instead of working with names.  It would look something like this
(this code is untested)

tst2 <- lapply(tst1, get)
names(tst2) <- tst1
do.call(anova, tst2[order(unlist(lapply(tst2, AIC)))])