Bias-corrected percentile confidence intervals

I'm trying to estimate bias-corrected percentile (BCP) confidence intervals in R on a vector from a simple for loop used for resampling. I am primarily looking for help implementing the calculation on a vector in R.

I am attempting to follow steps in Manly. 1998. Randomization, bootstrap and monte carlo methods in biology. 2nd edition. Pg. 48. PDF of the approach/steps should be available here: https://wyocoopunit.box.com/s/9vm4vgmbx5h7um809bvg6u7wr392v6i9

I aware of boot::boot but am hoping to avoid it for the current analysis I am working on (the boot function became quite challenging, for me, for a few reasons). I did post this question to R-help a little over a day ago but got no responses, presumably because it is too stats focused, so I'm reposting in what is probably a more appropriate forum.

I cannot figure out where I'm going wrong but the estimates from my attempt at the BCP CI are different enough from other methods that I assume I'm doing something wrong.

1) Bootstrap 95% CI for R-Squared via boot::boot

require(boot)
data("mtcars")

# Function for boot
rsq <- function(formula, data, indices) {
d <- data[indices,]
fit <- lm(formula, data=d)
return(summary(fit)$r.square) } # Bootstrap with 1000 replications and estimate CI's results <- boot(data=mtcars, statistic=rsq, R=1000, formula = mpg ~ wt + disp) rsq.ci <- boot.ci(results, type="all")  2) Bootstrap via for loop and estimate bias-corrected pecentile CI's iterations <- 1000 rsq.out <- vector() for(i in 1:iterations){ boot.df <- mtcars[sample(nrow(mtcars), nrow(mtcars), replace = TRUE), ] boot.lm <- lm(mpg ~ wt + disp, boot.df) boot.rsq <- summary(boot.lm)$r.squared
rsq.out <- c(rsq.out, boot.rsq) }

hist(rsq.out)


Fit LM with all the data (need for BCP steps)

lm.obs <- lm(mpg ~ wt + disp, mtcars)


Bias-corrected percentile CI
Step 1: Proportion of times that the boot estimate exceeds observed estimate (obviously are many ways to calculate this, mine is clearly not the most succint).

mtcar.boot.rsq <- data.frame(boot.rsq = rsq.out, obs.rsq = summary(lm.obs)$r.squared) # If boot mean is > observed mean, code 1, otherwise 0. mtcar.boot.rsq$boot.high <- ifelse(mtcar.boot.rsq$boot.rsq > mtcar.boot.rsq$obs.rsq, 1, 0)

# Mean is the proportion of times boot mean > obs mean
mean(mtcar.boot.rsq$boot.high) head(mtcar.boot.rsq)  Step 2: Use proportion to get Z score, then use that to calculate the new bias-correct Z score to look up the new proportion to use in quantile()  rsq.bc <- quantile(mtcar.boot.rsq$boot.rsq,
c(pnorm((2*pnorm(mean(mtcar.boot.rsq$boot.high))) - 1.96), pnorm((2*pnorm(mean(mtcar.boot.rsq$boot.high))) + 1.96)))


Estimates
Boot package

rsq.ci
Intervals :
Level      Normal              Basic
95%   ( 0.6690,  0.8659 )   ( 0.6818,  0.8824 )

Level     Percentile            BCa
95%   ( 0.6795,  0.8800 )   ( 0.6044,  0.8511 )


Simple percentile CI on for loop output

quantile(rsq.out, c(0.025, 0.975))
2.5%     97.5%
0.6871754 0.8809823


My attempt at BCP CI on for loop output (ehhhh)

rsq.bc
31.59952% -- 0.7747525
99.97103% -- 0.9034800


You almost had it. Change your Step 2 code as shown below. You want the value of Z associated with the proportion you computed in Step 1 and that is what qnorm will give you.

 rsq.bc <- quantile(mtcar.boot.rsq$boot.rsq, c(pnorm((2*pnorm(mean(mtcar.boot.rsq$boot.high))) - 1.96),
pnorm((2*pnorm(mean(mtcar.boot.rsq$boot.high))) + 1.96)))  becomes  rsq.bc <- quantile(mtcar.boot.rsq$boot.rsq,
c(pnorm((2*qnorm(mean(mtcar.boot.rsq$boot.high))) - 1.96), pnorm((2*qnorm(mean(mtcar.boot.rsq$boot.high))) + 1.96)))