The original code for Bishara and Hittner (2012) is preserved on my website (http://bisharaa.people.cofc.edu/modelcode/Code14o_public1.R). After publication, I noticed that there were two bugs in this code for bootstrapping with BCa. Fortunately, the bugs affected the results so little that the overall pattern remained unchanged: bootstrapping with BCa still led to inflation of Type I error in various scenarios. Still, if you would like to adapt code from this article for BCa bootstrapping, I have written a description of the bugs and corrections below. If the BCa does not matter to you, then you can disregard the rest of this message.
The bcafun in the article code was from pieces of the "bcanon" function in the "bootstrap" package. There are actually 2 problems with it, one from the original package, and one that I unfortunately introduced when I adapted it for our research. The problem in my adaptation is in the line:
n <- length(origdata)
Origdata is an nX2 matrix, so n stored was two times the value that it really should have been. This overly large n affected the acceleration term in later parts of the code.
The problem in the original packaged code (from the "bootstrap" package) is that it occasionally rounds to the 0th percentile for the lower bound, which is interpreted in R as the 0th element of a vector, which doesn’t exist. In 2013, I contacted the package maintainer about this issue, but I never received a reply.
In a more recent project, I started with the "bcanon" function, fixed the above-mentioned issue, and then made it more general so that it produced several different types of confidence interval bounds in addition to the BCa. I’ve pasted the code below, as well as advice from the package help about how to use it. Follow the example at the bottom, and you should be in good shape. If you have any questions, feel free to contact me. I realize that it’s always difficult to interpret someone else’s code.
#This is a generalized function expanded from bcanon() of Efron & Tib., #but it also includes CIs for the percentile BS in output.
#I have also added code to replace any BS sample that is all identical.
#I have also added code to fix the round-to-zero ooo problem #that leads to equal lo and hi CIs erratically for BCa. This is labeled:
#fixedBCa in the output.
bcanongeneral=function(x, nboot, theta, ..., alpha = c(0.025,0.975)) {
call <- match.call()
n <- length(x)
thetahat<-theta(x, ...)
bootsam <- matrix(sample(x, size = n * nboot, replace = TRUE),
nrow = nboot)
if(n<20) while(sum(apply(bootsam,1,sd)==0)>0)
{
rowstoreplace=which(apply(bootsam,1,sd)==0)
for (currow in rowstoreplace)
bootsam[currow,]=sample(x, size = n , replace = TRUE)
}
thetastar <- apply(bootsam, 1, theta, ...)
z0 <- qnorm(sum(thetastar < thetahat)/nboot)
accloops=n
u <- rep(0, accloops)
for (i in 1:accloops) {
u[i] <- theta(x[-i], ...)
}
uu <- mean(u) - u
acc <- sum(uu * uu * uu)/(6 * (sum(uu * uu))^1.5)
zalpha <- qnorm(alpha)
tt <- pnorm(z0 + (z0 + zalpha)/(1 - acc * (z0 + zalpha)))
ooo <- trunc(tt * nboot)
confpoints <- sort(thetastar)[ooo]
confpoints <- cbind(alpha, confpoints)
fixedooo=ooo
fixedooo[ooo==0]=1
fixedconfpoints <- sort(thetastar)[fixedooo]
#if(length(dimnames(confpoints)[[2]]) != 2) browser()
dimnames(confpoints)[[2]] <- c("alpha", "bca point")
percCI=quantile(thetastar,alpha,names=0)
AA=sqrt((n+2)/(n-1))*diff(percCI)
AACI=mean(percCI)+c(-AA/2,+AA/2)
return(list(confpoints = confpoints, fixedBCa=fixedconfpoints,
percCI=percCI, AACI=AACI,
z0 = z0, acc = acc, u = u, call = call)) }
# To obtain bca limits for functions of more # complex data structures, write theta # so that its argument x is the set of observation # numbers and simply pass as data to bcanongeneral # the vector 1,2,..n.
# For example, find bca limits for
# the correlation coefficient from a set of 15 data pairs:
xdata <- matrix(rnorm(30),ncol=2)
n <- 15
theta <- function(x,xdata){ cor(xdata[x,1],xdata[x,2]) }
results <- bcanongeneral (1:n,100,theta,xdata)