It was very statsy...I shoulda known better :) stats.stackexchange answer, for posterity: https://stats.stackexchange.com/questions/298506/bias-corrected-percentile-confidence-intervals
Joe On Wed, Aug 16, 2017 at 1:30 PM, Joe Ceradini <joecerad...@gmail.com> wrote: > Hi folks, > > I'm trying to estimate bias-corrected percentile (BCP) confidence > intervals on a vector from a simple for loop used for resampling. I am > attempting to follow steps in Manly, B. 1998. Randomization, bootstrap > and monte carlo methods in biology. 2nd edition., p. 48. PDF of the > approach/steps should be available here: > https://wyocoopunit.box.com/s/9vm4vgmbx5h7um809bvg6u7wr392v6i9 > If this is too statsy for the R list, I can try stackexchange or > crossval, but the R list provides all kinds of great help so thought > I'd try it first. 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 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. > > require(boot) > data("mtcars") > > # 1) Bootstrap 95% CI for R-Squared via boot::boot > # statmethods.net/advstats/bootstrapping.html > > # Function for boot > rsq <- function(formula, data, indices) { > d <- data[indices,] > fit <- lm(formula, data=d) > return(summary(fit)$r.square) > } > > # bootstrapping with 1000 replications > results <- boot(data=mtcars, statistic=rsq, > R=1000, formula = mpg ~ wt + disp) > > # get 95% confidence interval > rsq.ci <- boot.ci(results, type="all") > > > > # 2) Bootstrap via for loop and estimate bias-corrected pecentile CI's > > # Fit LM with all the data > lm.obs <- lm(mpg ~ wt + disp, mtcars) > > # Bootstrap with for loop > iterations <- 1000 # # of iterations > rsq.out <- vector() # to hold loop output > > 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) > > # Bias-corrected percentile CI > # Manly, B. 1998. Randomization, bootstrap and monte carlo methods in > biology. 2nd edition. > # - p. 48 > > # Proportion of times that the boot estimate exceeds observed estimate > mtcar.boot.rsq <- data.frame(boot.rsq = rsq.out, > obs.rsq = summary(lm.obs)$r.squared) > head(mtcar.boot.rsq) > > # 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) > > # 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 > # My attempt at BCP CI (ehhhh) > rsq.bc > 31.59952% -- 0.7747525 > 99.97103% -- 0.9034800 > > # 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 > quantile(rsq.out, c(0.025, 0.975)) > 2.5% 97.5% > 0.6871754 0.8809823 > > boot.ci(results, type="perc") > Level Percentile > 95% ( 0.6795, 0.8800 ) > > I appreciate any advice! > Thanks. > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
