2. Bootstrap Prediction Intervals for Regression

Why Two Bootstrap Schemes Behave Differently

A regression model's prediction interval must account for two distinct sources of uncertainty: uncertainty in the fitted mean function (parameter uncertainty) and the irreducible noise around that mean. The residual bootstrap and the case bootstrap attack these sources in different ways, and that difference is exactly what we want to measure in this case.

  • Residual bootstrap holds the design matrix X fixed, refits the model on y* = ŷ + r* where r* is a resample of the fitted residuals, and rebuilds the interval from the spread of ŷ* + r* at each query point. Because the design is fixed, it isolates noise uncertainty cleanly and is stable when the design has high leverage points.
  • Case bootstrap resamples whole rows (x_i, y_i), refits, and rebuilds the interval from the spread of ŷ* plus a residual term. It captures joint uncertainty in the design and the response, but it can be unstable when the sample is small or when a few rows dominate the fit.

The practical question is not which is "correct" but which gives stable width and honest coverage on held-out data for the model you actually have. We will build both, run them on the same data, and compare.

Two Resampling Paths from the Same Fitted Model Training set (X, y), n = 200 OLS fit → ŷ, residuals r = y − ŷ Residual bootstrap (fixed X) r* = resample(r, n) with replacement y* = ŷ + r* → refit β* ŷ*(x₀) + r* → quantiles Design held fixed → isolates noise Case bootstrap (resample rows) (x*, y*) = resample rows, n refit β* on (x*, y*) ŷ*(x₀) + r* → quantiles Design varies → joint uncertainty

Building Both Intervals on the Same Data

We use a synthetic regression with a known noise level so that coverage can be checked against a large independent test set. The code below is self-contained: it generates data, fits OLS, runs both bootstrap schemes with B = 1000 replicates, and reports empirical coverage and mean interval width at a 90% nominal level.

import numpy as np

rng = np.random.default_rng(42)

def make_data(n, rng):
    x = rng.uniform(-2, 2, n)
    y = 1.5 + 0.8 * x - 0.4 * x**2 + rng.normal(0, 0.6, n)
    return x, y

def design(x):
    return np.column_stack([np.ones_like(x), x, x**2])

def fit(X, y):
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    return beta

def residual_bootstrap(X, y, x0, B=1000, alpha=0.10, rng=None):
    beta = fit(X, y)
    yhat = X @ beta
    resid = y - yhat
    X0 = design(x0)
    preds = np.empty((B, len(x0)))
    for b in range(B):
        r_star = rng.choice(resid, size=len(y), replace=True)
        y_star = yhat + r_star
        beta_star = fit(X, y_star)
        preds[b] = X0 @ beta_star + rng.choice(resid, size=len(x0), replace=True)
    lo = np.quantile(preds, alpha / 2, axis=0)
    hi = np.quantile(preds, 1 - alpha / 2, axis=0)
    return lo, hi

def case_bootstrap(X, y, x0, B=1000, alpha=0.10, rng=None):
    n = len(y)
    X0 = design(x0)
    preds = np.empty((B, len(x0)))
    for b in range(B):
        idx = rng.integers(0, n, n)
        beta_star = fit(X[idx], y[idx])
        resid_star = y[idx] - X[idx] @ beta_star
        preds[b] = X0 @ beta_star + rng.choice(resid_star, size=len(x0), replace=True)
    lo = np.quantile(preds, alpha / 2, axis=0)
    hi = np.quantile(preds, 1 - alpha / 2, axis=0)
    return lo, hi

# Train / test split
x_tr, y_tr = make_data(200, rng)
x_te, y_te = make_data(2000, rng)
X_tr, X_te = design(x_tr), design(x_te)

lo_r, hi_r = residual_bootstrap(X_tr, y_tr, x_te, rng=rng)
lo_c, hi_c = case_bootstrap(X_tr, y_tr, x_te, rng=rng)

def report(name, lo, hi, y):
    cover = np.mean((y >= lo) & (y <= hi))
    width = np.mean(hi - lo)
    print(f"{name:20s} coverage={cover:.3f}  mean_width={width:.3f}")

report("residual bootstrap", lo_r, hi_r, y_te)
report("case bootstrap",     lo_c, hi_c, y_te)

With the seed above, the residual bootstrap typically lands near 0.90 coverage with a mean width around 2.0, while the case bootstrap is close but slightly wider and more variable across seeds. Treat those numbers as expected behavior for this synthetic setup, not as a benchmark — the point is the comparison protocol, which you should rerun with your own data.

Reading the Comparison Correctly

Coverage alone is not enough. A method that returns [−∞, +∞] has perfect coverage and zero usefulness. The two diagnostics that matter together are:

DiagnosticWhat it tells youFailure symptom
Empirical coverage on held-out dataWhether the nominal level is honestCoverage far below 0.90 → intervals too narrow; far above → too conservative
Mean interval widthEfficiency at the same coverageWidth blows up at the edges of the design → leverage sensitivity
Width stability across seedsReproducibility of the resampling schemeCase bootstrap width varies by >10% across seeds on small n
Width vs. x profileWhether uncertainty tracks the designFlat width where the model is extrapolating → understated uncertainty
Run both schemes with at least three seeds and report the range of coverage and width, not a single number. A scheme whose coverage swings from 0.86 to 0.94 across seeds is not reliable at n = 200, even if its average looks fine.

When Each Scheme Is the Right Choice

Use the residual bootstrap when the design is fixed by the experiment (designed experiments, sensor sweeps at preset levels) and you want a clean read on noise-driven uncertainty. Use the case bootstrap when rows are exchangeable samples from a population and you want the interval to reflect both parameter and design variability — but check width stability first, and fall back to residual bootstrap if the case version is erratic. If neither is stable at your sample size, that is a signal to increase B, collect more data, or move to a method with finite-sample guarantees rather than to tune the bootstrap further.

Do not compare the two schemes on the training set. In-sample coverage is inflated for both and will hide exactly the instability you are trying to detect. Always evaluate on rows the model never saw.
← Previous 📖 Contents Next →