Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Incorrect computation of acceleration parameter in BCA intervals #572

Open
bwiernik opened this issue Oct 21, 2022 · 0 comments
Open

Incorrect computation of acceleration parameter in BCA intervals #572

bwiernik opened this issue Oct 21, 2022 · 0 comments
Labels
Bug 🐛 Something isn't working

Comments

@bwiernik
Copy link
Contributor

Describe the bug
Currently, the BCa bootstrap intervals computed by bcai() compute the acceleration parameter from the bootstrap estimates and the mean estimate.

U <- (sims - 1) * (mean(x, na.rm = TRUE) - x)

This isn't consistent with the definition given by https://projecteuclid.org/journals/statistical-science/volume-11/issue-3/Bootstrap-confidence-intervals/10.1214/ss/1032280214.full (equations 6.5, 6.6, 6.7), which computes this parameter based on influence values/leave-one-out/jackknife values for the data points (not bootstrap replicates).

By comparison, boot::bca.ci() computes the a parameter using influence functions or jackknife procedures.

bca.ci <- function (boot.out, conf = 0.95, index = 1, t0 = NULL, t = NULL, 
  L = NULL, h = function(t) t, hdot = function(t) 1, hinv = function(t) t, 
  ...) 
{
  t.o <- t
  if (is.null(t) || is.null(t0)) {
    t <- boot.out$t[, index]
    t0 <- boot.out$t0[index]
  }
  t <- t[is.finite(t)]
  w <- qnorm(sum(t < t0)/length(t))
  if (!is.finite(w)) 
    stop("estimated adjustment 'w' is infinite")
  alpha <- (1 + c(-conf, conf))/2
  zalpha <- qnorm(alpha)
  if (is.null(L)) 
    L <- empinf(boot.out, index = index, t = t.o, ...)
  a <- sum(L^3)/(6 * sum(L^2)^1.5)
  if (!is.finite(a)) 
    stop("estimated adjustment 'a' is NA")
  adj.alpha <- pnorm(w + (w + zalpha)/(1 - a * (w + zalpha)))
  qq <- norm.inter(t, adj.alpha)
  cbind(conf, matrix(qq[, 1L], ncol = 2L), matrix(hinv(h(qq[, 
    2L])), ncol = 2L))
}

Thus, the results of bcai() are currently not correct, as they are not correctly computing the acceleration parameter.

cf. https://stats.stackexchange.com/questions/590257/using-the-bootstrap-estimates-to-compute-the-acceleration-parameter-for-bca-boot

@strengejacke strengejacke added the Bug 🐛 Something isn't working label Oct 22, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Bug 🐛 Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants