Styjun

I am looking for a way to calculate bias-corrected accelerated confidence intervals in R using a vector of bootstrapped results (which are bootstrap estimates of population growth rate - lambda). However, the packages I find are either made to use specific object types (as in the "boot" package) or do not calculate BCa type confidence intervals. The reason I have bootstrapped the results using a for loop and then stored the results in a vector is that for each bootstrap resample I first get a 80 x 33 matrix of results which defines parameters for each population in each year of sampling, which in turn define lambda for each population. As far as I can tell this would be cumbersome in the boot package and was easy to program as a for loop. The actual set of functions is rather complex and cannot be included here.

I did try using this question as a guide to fake a "boot" object, but it did not work: How can I use pre bootstrapped data to obtain a BCa confidence interval?.

Lets say I have my observed estimate of lambda

lambda = 1.18

and that we simulate a vector of bootstrapped estimates

library(fGarch)
lambdaBS = rsnorm(999,mean=lambda-0.04,sd=0.11,xi=2.5)
plot(density(lambdaBS))

which is right skewed and biased.

I am hoping that, using this information, there is a currently existing function which calculates BCa confidence intervals or else that it is easy to program a function to do so. As of yet, I have not found this to be the case.

As would be the case with R, where certain utilities are spread far and wide between packages, there was an easy solution, but it took me hours of searching to find, so I will answer my own question for anyone who may be searching for something similar.

Using the example data in the question, the bca function in the "coxed" R package give bias-corrected and accelerated confidence intervals for a vector of bootstrapped results. And we can compare them to other confidence intervals.

library(fGarch)
library(coxed)

set.seed(15438)

#simulate bootstrap statistics
lambdaBS = rsnorm(9999,mean=lambda-0.04,sd=0.11,xi=2.5)

#bias-corrected and accelerated
bca(lambdaBS)

1.002437 1.452525

#confidence intervals using standard error (inappropriate)
c(lambda-(sd(lambdaBS)*2),lambda+(sd(lambdaBS)*2))

0.9599789 1.4000211

#percentile confidence intervals
quantile(lambdaBS, c(0.025,0.975))

 2.5% 97.5% 
0.9895892 1.4016528 

This seems to work well. I am unsure of how it corrects for bias without requiring the initial estimate of the statistic in question, but have not read the paper this method is based on, yet.

Another simulation shows how this compares to results using boot and boot.ci.

library(boot)

#generate data
set.seed(12345)
dat = rsnorm(500,mean=1.6,sd=0.5,xi=3.0)

#bootstrap the median
meanfun = function(x,id) mean(x[id])
test = boot(data=dat,R=999,statistic=meanfun)

#BCa using boot.ci
boot.ci(test,type="bca")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = test, type = "bca")

Intervals : 
Level BCa 
95% ( 1.537, 1.626 ) 
Calculations and Intervals on Original Scale


#BCa using bca function from coxed package
bca(test$t)

 1.536888 1.625524

And in this case both functions give identical results.