Bootstrap BCa confidence limits with your own resamples

-----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?

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

R-help@

 mailing list

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