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

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)

u_adjusted = pnorm(z0 + (z0+zu)/(1-a_hat*(z0+zu)))

#Accelerated Bootstrap CI


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"))


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 <- mean(patch$$y)/mean(patch$$z)

z0 <- qnorm(sum(patch.ci$thetastar <= theta.hat)/B) alpha = 0.05 alpha.l.b <- pnorm(z0+z0 + qnorm(alpha)) alpha.u.b <- pnorm(z0+z0 + qnorm(1-alpha)) confint.bias <- quantile(patch.ci$thetastar, probs = c(alpha.l.b,alpha.u.b))

# Accelearation
acc.num <- sum((mean(patch.jack$$jack.values)-patch.jack$$jack.values)^3)

acc.denom <- 6*(((patch.jack$jack.se^2)*n/(n-1))^(3/2)) accelerate <- acc.num/acc.denom alpha.l <- pnorm(z0+(z0 + qnorm(alpha))/(1-accelerate*(z0 + qnorm(alpha)))) alpha.u <- pnorm(z0+(z0 + qnorm(1-alpha))/(1-accelerate*(z0 + qnorm(1-alpha)))) confint.bca <- 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.

• Thank you for breaking down the manual calculation of BCa! It answers my question perfectly. I will go through the presentation also.
– Dani
Commented Jun 22, 2020 at 14:07