[R] Bootstrap BCa confidence limits with your own resamples

[John Fox]

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 at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Frank Harrell
> Sent: Tuesday, March 12, 2013 11:27 AM
> To: r-help at r-project.org
> 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 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.