[R] Bootstrap BCa confidence limits with your own resamples

[Frank Harrell]

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