# How to: Prediction intervals for linear regression via bootstrapping

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.

• This is discussed in detail in the book by Davison and Hinkley, Bootstrap Methods and Their Application, along with an explicit algorithm (Algorithm 6.4). They explain concepts, pitfalls and details at greater length than is possible in a reasonable answer here. Dec 1, 2012 at 23:05
• @Glen_b Thanks for the reference. Unfortunately I am not in a university or company so I do not have the resources to acquire the book.
– Max
Dec 2, 2012 at 10:30
• It can be ordered from amazon; a full explanation of the algorithm and all the associated caveats and issues is not really the sort of thing you can cover in a few hundred words or even a one page answer. Dec 3, 2012 at 22:51
• @Glen_b I wrote up the algorithm from Davison and HInkely---there are several questions on CV about this, so I thought it was worth the effort. Any comments you have would be appreciated. stats.stackexchange.com/questions/226565/…
– Bill
Jan 3, 2017 at 15:17

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.

• Not sure how you build a confidence interval from this matrix of n.bs prediction intervals. Dec 3, 2012 at 21:06
• Also unsure how you then use the array of lwr and upr values to get your 'final' boostrapped prediction interval. Oct 20, 2020 at 16:57

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

• Thanks Ragy for the example. However from what I can see your answer did not cover the calculation of prediction intervals using bootstrap. As I said in my answer I already understand how to use bootstrapping to calculate confidence intervals - which your code appears to do.
– Max
Dec 2, 2012 at 9:27