[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