4
$\begingroup$

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"))
$\endgroup$

1 Answer 1

5
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.