[R] bootstrap bca confidence intervals for large number of statistics in one model; library("boot")

[Jacob Wegelin]

Sometimes one might like to obtain pointwise bootstrap bias-corrected,
accelerated (BCA) confidence intervals for a large number of statistics
computed from a single dataset.  For instance, one might like to get
(so as to plot graphically) bootstrap confidence bands for the fitted
values in a regression model.

(Example: Chiu S et al., Early Acceleration of Head Circumference in
Children with Fragile X Syndrome and Autism. Journal of Developmental &
Behavioral Pediatrics 2007;In press. In this paper we plot the
estimated trajectories of head circumference for two different
groups, computed by linear mixed-effects models, with confidence bands
obtained by bootstrap.)

Below is a toy example of how to do this using the "boot" library.
We obtain BCA CIs for all three regression parameters and for the fitted
values at 50 levels of the predictor.

set.seed(1234567)
x<-runif(150)
y<-2/3 + pi * x^2 + runif(length(x))/2
plot(x,y)
DAT<-data.frame(x,y)
NEWDATA<-data.frame(x=seq(min(x), max(x), length=50))
library('boot')
myfn<-function(data, whichrows) {
	TheseData<-data[whichrows,]
	thisLM<-lm( y~poly(x,2), data=TheseData)
	thisFit<-predict(thisLM, newdata=NEWDATA)
	c(
		coef(summary(thisLM))[,"Estimate"]
		, thisFit)
}
bootObj<-boot( data=DAT, statistic=myfn, R=1000 )
names(bootObj)
dim(bootObj$t)
sofar<-t(sapply( 1:ncol(bootObj$t), function(thiscolumn) boot.ci(bootObj, type='bca', index=thiscolumn)$bca[4:5] ))
colnames(sofar)<-c("lo", "hi")
NEWDATA<-cbind(NEWDATA, sofar[4:nrow(sofar),])
thecoefs<-sofar[1:3,]
lines( NEWDATA$x, NEWDATA$lo, col='red')
lines( NEWDATA$x, NEWDATA$hi, col='red')

Note: To get boot.ci to run with type='bca' it seems necessary to have a
large value of R.  Otherwise boot.ci returns an error, in my (limited)
experience.

Thanks in advance for any critiques.  (For instance, is there an easier or more elegant way?)

Jacob Wegelin