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 ______________________________________________ R-help@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.