How to: Bootstrap Prediction Intervals for Regression Models with non IID noise?

Question

Question:

Consider a regression model $Y=m(X)+\epsilon$ for which $\epsilon$ is neither independent of $x$ nor identically distributed.

How would we go about generating prediction intervals in such a scenario? I’m specifically asking in a non-time series context, though if the method is equally applicable there – then even better.

(I.e. A $1-\alpha$ prediction interval $[l,u]$ for $Y|X = x$ such that $\mathbb{P}(l\leq Y \leq u | X=x) = 1-\alpha$.)

Personal Attempt at Finding Answer:

My usual approach would be (for IID noise): generate $N$ bootstrapped samples. For each bootstrapped sample, fit the model and compute the residuals. Make a prediction with the fitted model, and randomly sample from the residuals and add this to the prediction. (This is the approach also used in the video in bullet 3).

However, in the non-IID case, I’ve done a lot of reading and can’t seem to find a concrete answer:

• The method offered here, based on Section 6.3.3 of Davidson and Hinckley (1997), I believe requires an OLS model and i.i.d,

• Other proposals require transforming the data (such as with heteroskedasticity) – however this would only work in certain circumstances.

• The approach provided in this video makes the most intuitive sense and works for i.i.d noise, but on testing on my own synthetic data that does not have i.i.d noise it’s not able to capture the true population distribution.

The most potentially promising bootstrap method I can find is in Cosma Shalizi’s Bootstrap Lecture Notes, see page 13, section 5.4. He describes the “Case Bootstrap” (aka. “Pairs bootstrap”, “x-y bootstrap” or “rows bootstrap”).

Resampling cases makes only very weak assumptions about the data-generating distribution, that all data points ($(x, y)$ pairs) are independent and identically distributed. It does not assume that any $m(x)$ is correct, or that the noise is independent of $x$, or has constant variance.

Now, as Shalizi mentions this x-y pairs method is extensively described by Buja et al (2014).

However, in neither Shalizi or Buja’s can I find a description of the actual method (in Shalizi’s case, I struggle to understand the R code), and more specifically in the context of prediction intervals.

Summary:

Does anyone know how to do the xy bootstrap for prediction intervals and/or know how to actually do prediction intervals for non iid noise?

Answer 1

You are correct that the case bootstrap is appropriate here. The algorithm itself is very simple. Assuming that your data set consists of n rows (observations), you want to repeat the following steps:

1. Sample with replacement n rows from your original data set, to obtain a new data set of the same size. Shalizi does this with the resample.data.frame function, which outputs df[resample(1:nrow(df)),]. He previously defines the resample function as sample(x, size=length(x), replace=TRUE), i.e. a sample of size x from the numbers {1,...,x} with replacement.

2. Fit your regression model to this new data set. In Shalizi's code, this is re.lm(cats.lm, resample.data.frame(cats)). The re.ml function just takes the equation from cats.lm (which is Hwt ~ Sex*Bwt) and fits it to the resampled data.

Every repetition will produce a different set of estimates for the model coefficients, and you can use these to approximate their sampling distribution.