[R] Bootstrap BCa confidence limits with your own resamples

[John Fox]

Dear Frank,

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Frank Harrell
> Sent: Tuesday, March 12, 2013 2:24 PM
> To: r-help at 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 at 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 at 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.