3
$\begingroup$

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

I appreciate any advice!

$\endgroup$
3
$\begingroup$

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

You might find this page helpful: http://influentialpoints.com/Training/bootstrap_confidence_intervals.htm#bias

$\endgroup$
  • $\begingroup$ Thank you @kbiolsi! I was close :). Also, not sure how I missed that website you added which shows examples of all kinds in R. I swear I searched high and low but kept getting directed to boot package examples. Thanks again for the prompt reply. $\endgroup$ – jCeradini Aug 17 '17 at 23:34
  • $\begingroup$ Seems odd to use a Wald type interval in this case. Rsq is Not normal or symmetric. Seems like you should take HPD rather than this symmetric +- interval. $\endgroup$ – pauljohn32 Aug 18 '17 at 3:30
  • $\begingroup$ Thanks for the feedback @pauljohn32. I'm not familiar with HPD, is that Highest Posterior Density in a Bayesian context? Just to be clear, I'm not actually working with r-squared but used it as a reproducible example...guess it wasn't the best choice of statistic. $\endgroup$ – jCeradini Aug 18 '17 at 16:07
  • $\begingroup$ @pauljohn32, this is out of my depth and beyond the scope of the post, but since it was brought up, don't the bias-corrected CI's account for lack of normality in the sampling distribution? I thought that was part of the reason it was preferred to the simple percentile method. I only have a conceptual understanding of most of these concepts so, when possible, conceptual explanations are most useful :). $\endgroup$ – jCeradini Aug 18 '17 at 18:18
  • $\begingroup$ Just wanted to add that the R code from influentialpoints (see answer kbiolsi) contains errors. E.g. a=sum(je^3)/(6*sum(je^2))^(3/2) should be a=sum(je^3)/(6*(sum(je^2))^(3/2)) or alternatively a=(1/6)*sum(je^3) / (sum(je^2))^(3/2). This is not the only error, so would recommend against using that code. $\endgroup$ – Koen Pouwels Oct 2 at 21:46

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.