I am having trouble to understand how to use bootstrapping to calculate prediction intervals for a linear regression model. Can somebody outline a step-by-step procedure? I searched via google but nothing really makes sense to me.
I do understand how to use bootstrapping for calculating confidence intervals for the model parameters.
Confidence intervals take account of the estimation uncertainty. Prediction intervals add to this the fundamental uncertainty. R's predict.lm will give you the prediction interval for a linear model. From there, all you have to do is run it repeatedly on bootstrapped samples.
n <- 100
n.bs <- 30
dat <- data.frame( x<-runif(n), y=x+runif(n) )
plot(y~x,data=dat)
regressAndPredict <- function( dat ) {
model <- lm( y~x, data=dat )
predict( model, interval="prediction" )
}
regressAndPredict(dat)
replicate( n.bs, regressAndPredict(dat[ sample(seq(n),replace=TRUE) ,]) )
The result of replicate is a 3-dimensional array (n x 3 x n.bs). The length 3 dimension consists of the fitted value for each data element, and the lower/upper bounds of the 95% prediction interval.
Gary King method
Depending on what you want, there's a cool method by King, Tomz, and Wittenberg. It's relatively easy to implement, and avoids the problems of bootstrapping for certain estimates (e.g. max(Y)).
I'll quote from his definition of fundamental uncertainty here, since it's reasonably nice:
A second form of variability, the fundamental un- certainty represented by the stochastic component (the distribution f ) in Equation 1, results from innumerable chance events such as weather or illness that may influ- ence Y but are not included in X. Even if we knew the ex- act values of the parameters (thereby eliminating esti- mation uncertainty), fundamental uncertainty would prevent us from predicting Y without error.
Bootstrapping does not assumed any knowledge of the form of the underlying parent distribution from which the sample arose. Traditional classical statistical parameter estimates are based on the normality assumption. Bootstrap deals with non-normality and is more accurate in practice than the classical methods.
Bootstrapping substitutes computers’ raw computing power for rigorous theoretical analysis. It is an estimate for the sampling distribution of a data set error term. Bootstrapping includes: re-sampling the data set a specified number of times, calculating the mean from each sample and finding the standard error of the mean.
The following “R” code demonstrates the concept:
This practical example demonstrates the usefulness of bootstrapping and estimates the standard error. The standard error is required to calculate confidence interval.
Let us assume you have a skewed data set "a":
a<-rexp(395, rate=0.1) # Create skewed data
visualization of the skewed data set
plot(a,type="l") # Scatter plot of the skewed data
boxplot(a,type="l") # Box plot of the skewed data
hist(a) # Histogram plot of the skewed data
Perform the bootstrapping procedure:
n <- length(a) # the number of bootstrap samples should equal the original data set
xbarstar <- c() # Declare the empty set “xbarstar” variable which will be holding the mean of every bootstrap iteration
for (i in 1:1000) { # Perform 1000 bootstrap iteration
boot.samp <- sample(a, n, replace=TRUE) #”Sample” generates the same number of elements as the original data set
xbarstar[i] <- mean(boot.samp)} # “xbarstar” variable collects 1000 averages of the original data set
##
plot(xbarstar) # Scatter plot of the bootstrapped data
boxplot(xbarstar) # Box plot of the bootstrapped data
hist(xbarstar) # Histogram plot of the bootstrapped data
meanOfMeans <- mean(xbarstar)
standardError <- sd(xbarstar) # the standard error is the standard deviation of the mean of means
confidenceIntervalAboveTheMean <- meanOfMeans + 1.96 * standardError # for 2 standard deviation above the mean
confidenceIntervalBelowTheMean <- meanOfMeans - 1.96 * standardError # for 2 standard deviation above the mean
confidenceInterval <- confidenceIntervalAboveTheMean + confidenceIntervalBelowTheMean
confidenceInterval
