Dear Frank, > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-bounces@r- > project.org] On Behalf Of Frank Harrell > Sent: Tuesday, March 12, 2013 2:24 PM > To: r-help@r-project.org > Subject: Re: [R] Bootstrap BCa confidence limits with your own > resamples > > 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?
If you run the example I included, you'll see that the BCa intervals differ from the simple percentile intervals. OTOH, I don't believe that I ever checked that the BCa intervals are correct for objects produced by bootSem() -- I did many years ago check against manual computations for a regression model fit by robust regression, but that's another matter. Best, John > > 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. ______________________________________________ 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.
