How to estimate the acceleration term used in boot.ci()

Question

I am trying to compute a bootstrapped CI for the Kendall's concordances statistic and need to present the acceleration term used.

What is the acceleration term used in the boot.ci() from R's boot package when using the BCa (adjusted bootstrap percentile) method? Say with the following code:

library(boot)
library(DescTools)
Create function to compute my estimator
my.estimator = function(data, i){ KendallW(data[i, c("var1", "var2")], correct=TRUE) }
R = 1000 #number of bootstrap resamples
Get the bootstrap object
b = boot(data, my.estimator, R)
Get confidence intervals
boot.ci(b, conf = 0.95, type = c("bca"))

It is not entirely obvious from the package description which method is used to estimate the acceleration term, but I think it is the usual jackknife. If so, does the following code (taken from a previous SE post) present the correct method to estimate the BCa confidence interval manually? These two methods did not provide the same intervals.

theta_hat = KendallW(data, correct=TRUE)
n = nrow(data)
I = rep(NA, n)
for(i in 1:n){
   #Remove ith data point
   xnew = data[-i, ]
   #Estimate theta
   theta_jack = KendallW(xnew, correct=TRUE)
   I[i] = (n-1)*(theta_hat - theta_jack)
}
#Estimate a
a_hat = (sum(I^3)/sum(I^2)^1.5)/6
Use this acceleration constant in own bootstrap algorithm
Desired quantiles
alpha = 0.05
u = c(alpha/2, 1-alpha/2)
B = 1000 #number of bootstrap resamples
theta_boot = rep(NA, B)
for(i in 1:B){
  #Select a bootstrap sample
  xnew = sample(data, length(data), replace=TRUE)
  #Estimate index
  theta_boot[i] = KendallW(xnew, correct=TRUE)
}
#Compute constants
z0 = qnorm(mean(theta_boot <= theta_hat))
zu = qnorm(u)
#Adjusted quantiles
u_adjusted = pnorm(z0 + (z0+zu)/(1-a_hat*(z0+zu)))
#Accelerated Bootstrap CI
quantile(theta_boot, u_adjusted)

A mock data is:

data = structure(list(var1 = structure(c(3, 1, 1, 1, 3, 0, 3, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 2, 1, 0, 2, 0, 0, 1, 1, 0, 0, 2, 1, 1, 0), label = "Variable 1", class = c("labelled", "numeric")), 
                      var2 = structure(c(1, 0, 0, 0, 3, 0, 3, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 2, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 3, 0, 0, 0, 0, 1, 0, 2, 1, 2, 0, 0, 0), label = "Variable 2", class = c("labelled", "numeric"))), 
                 row.names = c(NA, -50L), class = c("tbl_df", "tbl", "data.frame"))

Answer 1

The BC$_a$ confidence interval uses the percentiles of the bootstrap distribution, but corrects for the bias in the estimate i.e $\hat{\theta}$ as well as estimating the rate of change of the standard error.

It is important to note that the BC$_a$ confidence interval adjusts the percentiles, therefore the significance level ($\alpha$) you chose will be adjusted.

We introduce the bias-correction factor $\hat{z}_0$ as well as the acceleration factor $\hat{a}$

"The acceleration parameter estimates the rate of change of the standard error of $\hat{\theta}$ with respect to the true parameter $\theta$"

The acceleration factor can be denoted as follows: \begin{align*} \hat{a} = \frac{1}{6} \frac{\sum_{i=1}^{n}\left(\hat{\theta}_{(i)}-\hat{\theta}_{(\cdot)}\right)^3}{\left[\sum_{i=1}^{n}\left(\hat{\theta}_{i}-\hat{\theta}_{(\cdot)}\right)^2\right]^{3/2}} \end{align*}

There are some excellent notes on bootstrap confidence intervals available here. The explanations are clear and easy to understand.

In the code below is an example where I uses the patch data in the bootstrap library. I calculated the BC$_a$ interval manually. The above mentioned notes also have examples with implementation in R.

library(bootstrap)
library(boot)
Functions for the bootstrap and confidence interval
theta.func <- function(x){
y.val <- patch[x, "y"]
z.val <- patch[x, "z"]
mean(y.val)/mean(z.val)
}

confint90 <- function(x){
quantile(x,probs=c(0.05,0.95))
}


Bootstrapping
patch.ci = bootstrap(1:n, nboot=B, theta=theta.func,func = confint90)


Jackknife
    xdata <- patch
theta.jack <- function(x,xdata){
    y.jack = xdata[x,"y"]
    z.jack = xdata[x,"z"]
mean.jack = mean(y.jack)/mean(z.jack)
}

patch.jack <- jackknife(1:n, theta=theta.jack, xdata)


Bias correction
theta.hat &lt;- mean(patch<span class="math-container">$y)/mean(patch$</span>z)

z0 &lt;- qnorm(sum(patch.ci$thetastar &lt;= theta.hat)/B)
alpha = 0.05
alpha.l.b &lt;- pnorm(z0+z0 + qnorm(alpha))
alpha.u.b &lt;- pnorm(z0+z0 + qnorm(1-alpha))

confint.bias &lt;- quantile(patch.ci$thetastar, probs = c(alpha.l.b,alpha.u.b))


Accelearation
acc.num &lt;- sum((mean(patch.jack<span class="math-container">$jack.values)-patch.jack$</span>jack.values)^3)

acc.denom &lt;- 6*(((patch.jack$jack.se^2)*n/(n-1))^(3/2))

accelerate &lt;- acc.num/acc.denom 

alpha.l &lt;- pnorm(z0+(z0 + qnorm(alpha))/(1-accelerate*(z0 + qnorm(alpha))))

alpha.u &lt;- pnorm(z0+(z0 + qnorm(1-alpha))/(1-accelerate*(z0 + qnorm(1-alpha))))

confint.bca &lt;- quantile(patch.ci$thetastar, probs = c(alpha.l,alpha.u))

I calculated a confidence interval for a ratio in the patch data. You can see some the adjustment in the bias-corrected and BC$_a$ compared to the original percentile interval.

percentile      -0.210967  0.113875 
bias-corrected  -0.206047  0.124526 
BCa             -0.201898  0.133085 

If you don't want to calculate BC$_a$ manually, you can use the bca function in the coxed library.