Search This Blog

Monday, April 02, 2012

R: Bootstrap-t confidence intervals

The theory behind bootstrapping a pivotal quantity always eluded me.  The simplest way to approach bootstrapping a statistic is to generate B bootstrap samples of the original data and calculate B statistics from these samples.  Then take the alpha/2 and 1-alpha/2 quantiles of the statistic to form a confidence interval of level alpha. But in some circumstances, more than alpha of the intervals might not contain the true value.  This can be shown through a simple simulation.
Take 20 random exponential variables with mean 3.  In R this looks like:
x = rexp(20,rate=1/3)
Then generate B=1000 bootstrap samples of x, and calculate the mean for each bootstrap sample.
s = numeric(B)
for (j in 1:B) {
boot = sample(n,replace=TRUE)
s[j] = mean(x[boot])
}

Then, for an alpha = .05 / 95% confidence interval, look at the .025 and .975 quantiles of the bootstrap statistics in the vector s:
simple.ci = quantile(s,c(.025,.975))
If I repeat this process from the start (including drawing a new x of 20 random exponential variables of mean 3) I can see how often the intervals actually contain the true mean.  Here are 100 replicates of the whole interval creation process:

11 of the intervals, highlighted in red, do not contain the true mean 3, the blue vertical line.  On average we would expect 5 if these are 95% confidence intervals.  If I repeat this 1000 times, I get 88.4% of the intervals containing the true mean.
The bootstrap-t interval is a way of dealing with this problem.  In An Introduction to the Bootstrap by Efron and Tibshirani, the authors write,
The quantity (theta-hat – theta)/se-hat is called an approximate pivot: this means that its distribution is approximately the same for each value of theta….
Some elaborate theory shows that in large samples the coverage of the bootstrap-t interval tends to be closer to the desired level than the coverage of the standard interval or the interval based on the t table….
Notice also that the normal and t percentage points are symmetric about zero, and as a consequence the resulting intervals are symmetric about the point estimate theta-hat.  In contrast, the bootstrap-t percentiles can be asymmetric about 0, leading to intervals which are longer on the left or right.  This asymmetry represents an important part of the improvement in coverage it enjoys.
Note the authors here are not comparing the bootstrap-t interval to the bootstrap percentile method. Also, they add a caveat: “The bootstrap-t method, at least in its simple form, cannot be trusted for more general problems, like setting a confidence interval for a correlation coefficient.”
But it works well for calculating the mean, as in this case. This time we calculate a pivotal quantity as the bootstrapped statistic. For the vector x of random exponential variables of mean 3, we first calculate the mean and standard deviation of the original dataset:
x = rexp(n,rate=1/true.mean)
mean.x = mean(x)
sd.x = sd(x)

Then for each bootstrap sample, calculate the difference between the mean and the bootstrap mean, divided by the standard deviation of the bootstrap.
z = numeric(B)
for (j in 1:B) {
boot = sample(n,replace=TRUE)
z[j] = (mean.x - mean(x[boot]))/sd(x[boot])
}

Then to form a confidence interval, take quantiles of our bootstrapped statistic, multiply by the standard deviation of the original data and add this to the mean of the original data:
pivot.ci = mean.x + sd.x*quantile(z,c(.025,.975))
Now there are only 4 out of 100 intervals that do not contain the true mean:

If I do this 1000 times, I get 95.7% of the intervals containing the true mean. The mean interval size has increased though, from 2.4 for the simple intervals to 3.2 for the bootstrap-t intervals.
Here is an R script for this simulation.

R: Bootstrap-t Confidence Limits


boott {bootstrap}R Documentation

Bootstrap-t Confidence Limits

Description

See Efron and Tibshirani (1993) for details on this function.

Usage

boott(x,theta, ..., sdfun=MISSING, nbootsd=25, nboott=200,
      VS=FALSE, v.nbootg=100, v.nbootsd=25, v.nboott=200,
      perc=c(.001,.01,.025,.05,.10,.50,.90,.95,.975,.99,.999))

Arguments

x a vector containing the data. Nonparametric bootstrap sampling is used. To bootstrap from more complex data structures (e.g. bivariate data) see the last example below.
theta function to be bootstrapped. Takes x as an argument, and may take additional arguments (see below and last example).
... any additional arguments to be passed to theta
sdfun optional name of function for computing standard deviation of theta based on data x. Should be of the form: sdmean <- function(x,nbootsd,theta,...) where nbootsd is a dummy argument that is not used. If theta is the mean, for example, sdmean <- function(x,nbootsd,theta,...) {sqrt(var(x)/length(x))} . If sdfun is missing, then boott uses an inner bootstrap loop to estimate the standard deviation of theta(x)
nbootsd The number of bootstrap samples used to estimate the standard deviation of theta(x)
nboott The number of bootstrap samples used to estimate the distribution of the bootstrap T statistic. 200 is a bare minimum and 1000 or more is needed for reliable α % confidence points, α > .95 say. Total number of bootstrap samples is nboott*nbootsd.
VS If TRUE, a variance stabilizing transformation is estimated, and the interval is constructed on the transformed scale, and then is mapped back to the original theta scale. This can improve both the statistical properties of the intervals and speed up the computation. See the reference Tibshirani (1988) given below. If FALSE, variance stabilization is not performed.
v.nbootg The number of bootstrap samples used to estimate the variance stabilizing transformation g. Only used if VS=TRUE.
v.nbootsd The number of bootstrap samples used to estimate the standard deviation of theta(x). Only used if VS=TRUE.
v.nboott The number of bootstrap samples used to estimate the distribution of the bootstrap T statistic. Only used if VS=TRUE. Total number of bootstrap samples is v.nbootg*v.nbootsd + v.nboott.
perc Confidence points desired.

Value

list with the following components:
confpoints Estimated confidence points
theta, g theta and g are only returned if VS=TRUE was specified. (theta[i],g[i]), i=1,length(theta) represents the estimate of the variance stabilizing transformation g at the points theta[i].

References

Tibshirani, R. (1988) "Variance stabilization and the bootstrap". Biometrika (1988) vol 75 no 3 pages 433-44.
Hall, P. (1988) Theoretical comparison of bootstrap confidence intervals. Ann. Statisi. 16, 1-50.
Efron, B. and Tibshirani, R. (1993) An Introduction to the Bootstrap. Chapman and Hall, New York, London.

Examples

#  estimated confidence points for the mean
x <- rchisq(20,1)
theta <- function(x){mean(x)}
results <- boott(x,theta)
# estimated confidence points for the mean, 
#  using variance-stabilization bootstrap-T method
results <-  boott(x,theta,VS=TRUE)
results$confpoints          # gives confidence points
# plot the estimated var stabilizing transformation
plot(results$theta,results$g) 
# use standard formula for stand dev of mean
# rather than an inner bootstrap loop
sdmean <- function(x, ...) 
    {sqrt(var(x)/length(x))}
results <-  boott(x,theta,sdfun=sdmean) 
                                     
# To bootstrap 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 boot the vector 1,2,..n. 
# For example, to bootstrap
# 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 <- boott(1:n,theta, xdata)