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@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.