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

[Ken Nussear]

On May 24, 2007, at 2:12 PM, Douglas Bates wrote:

> 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)))])



I get errors with each method that I'm not sure how to solve. Any Ideas?

Method 1

 > do.call(anova, lapply(tst3, get))
Error in `row.names<-.data.frame`(`*tmp*`, value = c("structure(list 
(modelStruct = structure(list(reStruct = structure(list(",  :
	invalid 'row.names' length


Method 2
 > names(tst2) <- tst1
 > tst2
$lme.T972way.ml
[1] 507.0991

$lme.T97FULL.ml
[1] 508.7594

$lme.T97NOINT.ml
[1] 564.8574

$lme.T97NULL.ml
[1] 624.3053

$lme.T97fc.ml
[1] 502.5878

$lme.T97min.ml
[1] 569.8188

$lme.T97ns.ml
[1] 504.8971

 > do.call(anova, tst2[order(unlist(lapply(tst2, AIC)))])
Error in logLik(object) : no applicable method for "logLik"
 >