9
$\begingroup$

I have noticed different software use different method of calculating the acceleration parameter, $a$, for computing the bias-corrected and accelerated (BCa) bootstrap confidence interval (CI). In general, the formula for $a$ is $$ a = \frac{1}{6}\frac{\sum_i{I_i^3}}{(\sum_i{I_i^2})^{3/2}} $$ where $i$ index units in the sample and $I$ is the "empirical influence function". Calculating $I$ is where the approaches differ.

The boot package in R offers four options:

  1. The infinitesimal jackknife
  2. The usual (leave one out) jackknife
  3. The positive (include one twice) jackknife
  4. A regression estimation method, which involves use the coefficients from a regression of the estimates on the bootstrap array (which computes the effect of the number of times each unit is present in the bootstrap sample on the bootstrap estimate in that sample)

A fifth option is in the bayestestR package, which seems to use the bootstrapped statistics directly, i.e., $I_j=(R-1)(\bar{t}-t_j)$, where $R$ is the number of bootstrap replications, $t_j$ is the $j$th bootstrap estimate, and $\bar{t}$ is their mean. In this method, it appears the sums are over the bootstrap replications $j$ rather than over the units $i$. This method is much faster than the others because all the values required for calculating the acceleration parameter come straight from the bootstrap estimates, whereas the other methods require additional calculations (e.g., performing the jackknife). See bayestestR:::.bci() for the implementation.

Although it is clear that this method was developed for Bayesian inference (given its inclusion in a package for Bayesian inference), other packages use it for frequentist bootstrap inference as well. Is this valid?

$\endgroup$

1 Answer 1

7
$\begingroup$

I've been wondering the same, and while I cannot provide any theoretical justification either way the answer seems to be 'no'.

The key part is that $a$ is a function of $I$ for each observation, and in a bootstrap sample this information was essentially discarded (it may or may not contain an individual observation any number of times). The regression approach within boot gets around this by counting the original observations that occur in each sample and using that information for a design matrix.

Empirically these calculations also do not match, the bootstrap estimate is always off by at least an order of magnitude or two - see code below. Probably not exactly the answer you (or I) were looking for, hopefully someone can shed some light on the theory here.

set.seed(1)

## Statistic of interest
statfun <- mean

## Original sample
orig_n <- 42
orig_obs <- rexp(orig_n)

## Jackknife estimate of acceleration
jack_dist <- rep(NA, orig_n)

for (i in seq_len(orig_n)) {
   jack_dist[i] <- statfun(orig_obs[-i])
}

jack_mean <- mean(jack_dist)
jack_I <- jack_mean - jack_dist
jack_a <- sum(jack_I**3) / (6*sum(jack_I**2)**(3/2))

## Bootstrap estimate of acceleration
boot_n <- 1e4
boot_dist <- rep(NA, boot_n)

for (i in seq_len(boot_n)) {
   boot_dist[i] <- statfun(sample(orig_obs, orig_n, replace=TRUE))
}

boot_mean <- mean(boot_dist)
boot_I <- (boot_n - 1) * (boot_mean - boot_dist)
boot_a <- sum(boot_I**3) / (6*sum(boot_I**2)**(3/2))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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