That's very helpful John. Did you happen to run a test to make sure that boot.ci(..., type='bca') in fact gives the BCa intervals or that they at least disagree with percentile intervals?
Frank John Fox wrote > Dear Frank, > > I'm not sure that it will help, but you might look at the bootSem() > function > in the sem package, which creates objects that inherit from "boot". Here's > an artificial example: > > ---------- snip ---------- > > library(sem) > for (x in names(CNES)) CNES[[x]] <- as.numeric(CNES[[x]]) > model.cnes <- specifyModel() > F -> MBSA2, lam1 > F -> MBSA7, lam2 > F -> MBSA8, lam3 > F -> MBSA9, lam4 > F <-> F, NA, 1 > MBSA2 <-> MBSA2, the1 > MBSA7 <-> MBSA7, the2 > MBSA8 <-> MBSA8, the3 > MBSA9 <-> MBSA9, the4 > > sem.cnes <- sem(model.cnes, data=CNES) > summary(sem.cnes) > > set.seed(12345) # for reproducibility > system.time(boot.cnes <- bootSem(sem.cnes, R=5000)) > class(boot.cnes) > boot.ci(boot.cnes) > > ---------- snip ---------- > > I hope this helps, > John > >> -----Original Message----- >> From: > r-help-bounces@ > [mailto:r-help-bounces@r- >> project.org] On Behalf Of Frank Harrell >> Sent: Tuesday, March 12, 2013 11:27 AM >> To: > r-help@ >> Subject: [R] Bootstrap BCa confidence limits with your own resamples >> >> I like to bootstrap regression models, saving the entire set of >> bootstrapped >> regression coefficients for later use so that I can get confidence >> limits >> for a whole set of contrasts derived from the coefficients. I'm >> finding >> that ordinary bootstrap percentile confidence limits can provide poor >> coverage for odds ratios for binary logistic models with small N. So >> I'm >> exploring BCa confidence intervals. >> >> Because the matrix of bootstrapped regression coefficients is generated >> outside of the boot packages' boot() function, I need to see if it is >> possible to compute BCa intervals using boot's boot.ci() function after >> constructing a suitable boot object. The function below almost works, >> but >> it seems to return bootstrap percentile confidence limits for BCa >> limits. >> The failure is probably due to my being new to BCa limits and not >> understanding the inner workings. Has anyone successfully done >> something >> similar to this? >> >> Thanks very much, >> Frank >> >> bootBCa <- function(estimate, estimates, n, conf.int=0.95) { >> require(boot) >> if(!exists('.Random.seed')) runif(1) >> w <- list(sim= 'ordinary', >> stype = 'i', >> t0 = estimate, >> t = as.matrix(estimates), >> R = length(estimates), >> data = 1:n, >> strata = rep(1, n), >> weights = rep(1/n, n), >> seed = .Random.seed, >> statistic = function(...) 1e10, >> call = list('rigged', weights='junk')) >> np <- boot.ci(w, type='perc', conf=conf.int)$percent >> m <- length(np) >> np <- np[c(m-1, m)] >> bca <- tryCatch(boot.ci(w, type='bca', conf=conf.int), >> error=function(...) list(fail=TRUE)) >> if(length(bca$fail) && bca$fail) { >> bca <- NULL >> warning('could not obtain BCa bootstrap confidence interval') >> } else { >> bca <- bca$bca >> m <- length(bca) >> bca <- bca[c(m-1, m)] >> } >> list(np=np, bca=bca) >> } >> >> >> >> >> ----- >> Frank Harrell >> Department of Biostatistics, Vanderbilt University >> -- >> View this message in context: http://r.789695.n4.nabble.com/Bootstrap- >> BCa-confidence-limits-with-your-own-resamples-tp4661045.html >> Sent from the R help mailing list archive at Nabble.com. >> >> ______________________________________________ >> > R-help@ > 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. > > ______________________________________________ > R-help@ > 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. ----- Frank Harrell Department of Biostatistics, Vanderbilt University -- View this message in context: http://r.789695.n4.nabble.com/Bootstrap-BCa-confidence-limits-with-your-own-resamples-tp4661045p4661077.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ 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.