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

From: Jacob Wegelin <jawegelin_at_ucdavis.edu>
Date: Fri 26 Jan 2007 - 20:26:04 GMT

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



R-help@stat.math.ethz.ch mailing list
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. Received on Sat Jan 27 07:29:18 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 31 Jan 2007 - 19:31:40 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.