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

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$$.)

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?

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.

• Hi Doctor Milt, thanks for taking the time to respond and the explanation. Isn’t the case bootstrap method as you describe it though only going to provide confidence intervals around the model coefficients? I’m interested in a prediction interval around a new observation Feb 7, 2023 at 12:42
• I think you'll have to move beyond standard bootstrap techniques, because you're trying to predict what will happen when you plug in a new (unseen) value of x, but ϵ depends on x. Have you looked at (split sample) conformal prediction intervals? Feb 7, 2023 at 13:47
• No I can’t say I have, I’ll definitely look into that. Do you have any helpful resources or know where I might find a step by step of the process? Thanks again Feb 7, 2023 at 14:29
• Feb 7, 2023 at 14:55
• Just coming back to this to say that this last month I’ve taken a full dive into conformal prediction thanks to your comment — it’s definitely the solution to the question I was asking. If you’re able to submit that as an answer, happy to give you the tick :) Thanks Mar 7, 2023 at 23:43