[R] Re strict AIC comparison to succesful models?
Manuel Morales
Manuel.A.Morales at williams.edu
Thu Jun 11 22:54:25 CEST 2009
On Thu, 2009-06-11 at 16:10 -0400, Ben Bolker wrote:
> Manuel Morales wrote:
> > On Thu, 2009-06-11 at 12:08 -0700, Ben Bolker wrote:
> >>
> >> Manuel Morales wrote:
> >>> Hello list,
> >>>
> >>> I'm doing a bootstrap analysis where some models occasionally fail to
> >>> converge. I'd like to automate the process of restricting AIC to the
> >>> models that do converge. A contrived example of what I'd like to do is
> >>> below:
> >>>
> >>> resp <- c(1,1,2)
> >>> pred <- c(1,2,3)
> >>>
> >>> m1 <- lm(resp~pred)
> >>> m2 <- lm(resp~poly(pred,2))
> >>> m3 <- lm(resp~poly(pred,3)) # Fails, obviously
> >>>
> >>> ## Some test to see which models were successful
> >>> models <- test(m1,m2,m3)
> >>>
> >> How about
> >>
> >> models <- list(m1,m2,m3)
> >> is.OK <- sapply(models,test)
> >> do.call(AIC,models[is.OK])
> >>
> >> ?
> >
> > Good idea but unfortunately I get an error message with:
> >
> > models <- list(m1,m2,m3)
> > test <- function(x) length(x)>1
> > is.OK <- sapply(models,test)
> > do.call(AIC,models[is.OK])
> >
> > I think the issue is that AIC is expecting a ... of model objects (not a
> > list or a vector) and I'm not sure how to construct that.
> >
> > Manuel
> >
>
> The problem is that when your fitting function fails,
> NOTHING gets created.
>
> Try this:
>
> #####################
> resp <- c(1,1,2)
> pred <- c(1,2,3)
>
> m1 <- try(lm(resp~pred))
> m2 <- try(lm(resp~poly(pred,2)))
> m3 <- try(lm(resp~poly(pred,3))) # Fails, obviously
>
> models <- list(m1,m2,m3)
> test <- function(x) !inherits(x,"try-error")
> is.OK <- sapply(models,test)
>
> AIC(m1,m2) ## works
>
> ## all these work but lead to objects with ugly row names
> do.call(AIC,models[is.OK])
> do.call("AIC",models[is.OK])
> do.call(stats:::AIC.default,models[is.OK])
> dd <- do.call(AIC,list(m1,m2))
> rownames(dd) <- NULL
> dd
>
I see, it wasn't an error message just a very ugly output. Reassigning the row names works great!
Thanks,
Manuel
--
http://mutualism.williams.edu
More information about the R-help
mailing list