The bootstrap has a dangerous user interface.

It looks like this:

sample with replacement

That sentence is so small that it hides the whole model. What is being sampled? Rows? Residuals? Clusters? Blocks of time? Scores? Treatment assignments? If the answer is wrong, the bootstrap does not complain. It keeps drawing samples and returns precise-looking intervals from the wrong imagined world.

The bootstrap is not a machine for avoiding assumptions. It is a machine for making assumptions runnable.

The Tiny Spell With a Hidden Clause

Efron’s 1979 paper introduced the bootstrap as a general way to estimate the sampling distribution of a statistic from the empirical distribution of the observed data.1 If the data are an i.i.d. sample

\[X_1,\ldots,X_n \sim F,\]

then the empirical distribution

\[\widehat F_n = \frac{1}{n}\sum_{i=1}^n \delta_{X_i}\]

is a natural plug-in estimate of \(F\). A bootstrap sample draws \(n\) points from \(\widehat F_n\) with replacement, computes the statistic again, and repeats that many times. The resulting cloud of statistics approximates the sampling distribution.

That is beautiful because it turns an analytic question into a computational one:

what would this estimator do if the empirical distribution were the population?

But the sentence has a hidden clause:

... and if the original sampling unit really was exchangeable

The trouble starts when observations are not exchangeable in the way your code pretends.

Regression Has Multiple Worlds

Freedman’s regression bootstrap paper is still a wonderful antidote to casual resampling.2 The important line is conceptual: the bootstrap simulation has to reflect the relevant features of the stochastic model that generated the data.

In regression, there are at least two common stories.

First, the design matrix may be fixed. The analyst chose the \(x_i\) values, and randomness enters through homoskedastic errors:

\[Y_i = \alpha + \beta X_i + \epsilon_i,\qquad \epsilon_i \stackrel{i.i.d.}{\sim} G.\]

If this is the story, resampling centered residuals and attaching them back to the fitted values can make sense. The residuals are treated as exchangeable draws from one common error distribution.

Second, the rows may be random:

\[(X_i,Y_i) \stackrel{i.i.d.}{\sim} P.\]

Now the row is the sampling unit. Resampling pairs \((X_i,Y_i)\) better matches the data collection process. It lets the distribution of \(X\) vary across bootstrap samples, and it preserves any relationship between high-leverage points and their outcomes.

Those are not cosmetic choices. They are different worlds.

Where Residuals Betray You

The simulator below creates a simple linear regression with true slope 1. The error variance can grow with \(|x|\), and a fraction of observations can live at high-leverage \(x\) values. It then runs many experiments and asks how often different nominal 95% intervals actually cover the true slope.

The five intervals are:

  • classic OLS normal interval with homoskedastic standard errors;
  • HC1 normal interval using White-style heteroskedasticity-consistent standard errors;
  • residual bootstrap, which resamples centered residuals;
  • pairs bootstrap, which resamples whole rows;
  • wild bootstrap, which keeps each row’s residual scale and multiplies the residual by random signs.
Classic HC1 Residual Pairs Wild

Deterministic Monte Carlo. Each setting generates repeated synthetic regression datasets with true slope 1, then computes nominal 95% intervals. The bootstrap intervals use the basic bootstrap form around the observed slope. Results are pedagogical, not a universal ranking of methods.

At the default setting, the residual bootstrap tends to undercover. It treats residuals from high-variance, high-leverage regions as if they could have come from anywhere. The resampling step smears the error distribution across the design and invents a more uniform world than the one that generated the data.

The pairs bootstrap usually does better here because the row is the sampling unit. High-leverage rows stay high-leverage rows. The wild bootstrap usually does better because it keeps each fitted residual attached to its row and only randomizes its sign. It says: “I believe the conditional mean structure enough to keep the fitted design, but I do not believe all residuals share one variance.”

Now set Heteroskedasticity to zero and Leverage share to zero. The methods move closer together. That is not an accident. Under a fixed-design, homoskedastic error story, the residual bootstrap’s exchangeability assumption is much less offensive.

The point is not “always use pairs” or “always use wild.” The point is:

the resampling unit is the assumption

The Sandwich Answers in Closed Form

White’s heteroskedasticity-consistent covariance estimator attacked a related problem analytically.3 Ordinary least squares coefficients can remain consistent under heteroskedasticity, but the usual covariance matrix can be wrong. The sandwich estimator replaces the constant-variance middle of the variance formula with residual-weighted score products.

For a simple slope with intercept, the robust variance has the same spirit as:

\[\widehat{\operatorname{Var}}(\widehat \beta) \propto \frac{\sum_i (x_i-\bar x)^2 \widehat \epsilon_i^2} {\left(\sum_i (x_i-\bar x)^2\right)^2}.\]

High-leverage observations with large residuals now matter more. The estimator is not trying to pretend all residuals came from the same urn.

This is why the lab includes HC1 next to the bootstraps. Sometimes a robust analytic standard error is the right simple answer. Sometimes a bootstrap is useful because the statistic is complicated, the finite-sample distribution is awkward, or a design-specific resampling scheme captures the experiment better than a first-order approximation. They are tools for the same question:

what randomness should this inference procedure replay?

Wild Bootstrap Keeps the Row’s Scale

The wild bootstrap is a useful compromise for heteroskedastic regression. In its simple Rademacher form, after fitting the model it generates

\[Y_i^\star = \widehat Y_i + v_i \widehat \epsilon_i,\]

where \(v_i\) is independently \(+1\) or \(-1\) with equal probability. The residual magnitude remains tied to the original row. The sign changes.

Wu studied resampling methods for regression and emphasized that methods justified in the i.i.d. case may fail in regression settings with special structure.4 Later work associated with Liu and Mammen developed the wild bootstrap as a way to handle non-i.i.d. errors and heteroskedastic regression problems.5 A modern score-bootstrap view says the wild bootstrap perturbs the fitted scores while keeping the Hessian fixed; that interpretation makes the method easier to generalize beyond ordinary residuals.6

The implementation in the lab is intentionally plain. It does not studentize the bootstrap statistic, does not use Mammen’s two-point weights, and does not claim second-order refinements. It just shows the core mechanism:

do not detach residual scale from the row when residual scale belongs to the row

Rows Are Not Always the Unit

Rows are not always the unit either.

In a time series, resampling individual observations usually destroys serial dependence. If returns cluster by volatility regime, if users arrive in sessions, if server latency has bursts, if panel observations share a household, school, cluster, desk, or geography, the i.i.d. row bootstrap may create datasets that no real sampling process could have produced.

Kunsch’s block bootstrap for stationary observations is one classic response: resample blocks rather than individual observations so short-range dependence is preserved.7 Cluster bootstraps and randomization inference make the same kind of move in other designs. The unit changes because the source of randomness changed.

This is the bootstrap’s deepest practical lesson. You do not ask:

can I bootstrap this?

You ask:

what would have been exchangeable before I saw the data?

Then you resample that.

Before You Resample

Before I trust a bootstrap interval, I want written answers to:

  • What is the sampling unit?
  • Are observations independent, clustered, stratified, paired, blocked, or serially dependent?
  • Are regressors fixed by design or random draws from a population?
  • Is variance plausibly constant conditional on the regressors?
  • Is the statistic smooth enough for ordinary bootstrap approximations?
  • Does the resampling scheme preserve treatment assignment, matching, stratification, censoring, or survey weights?
  • Does the interval use percentile, basic, studentized, BCa, or another form?
  • How many bootstrap draws are enough for the reported tail probability?
  • What diagnostic simulation would fail if the resampling contract is wrong?

The last question is the one I like most. If you cannot name a world where your bootstrap fails, you probably have not named the world where it succeeds.

The Line I Would Keep

The bootstrap is often taught as a recipe. It is better understood as a tiny simulator of the experiment you believe happened.

If the data really are exchangeable rows, row resampling is natural. If the design is fixed and errors are homoskedastic, residual resampling may be natural. If the errors are heteroskedastic, wild bootstrap or robust covariance estimates may be natural. If the data are dependent, blocks or clusters may be the unit.

The code can be one line. The assumption is not.

resample the unit that nature randomized

That is the part worth remembering.

Paper Trail

  1. Bradley Efron, “Bootstrap Methods: Another Look at the Jackknife”, The Annals of Statistics, 1979. DOI: 10.1214/aos/1176344552. ↩

  2. David A. Freedman, “Bootstrapping Regression Models”, The Annals of Statistics, 1981. DOI: 10.1214/aos/1176345638. ↩

  3. Halbert White, “A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity”, Econometrica, 1980. RePEc entry: EconPapers. ↩

  4. C. F. J. Wu, “Jackknife, Bootstrap and Other Resampling Methods in Regression Analysis”, The Annals of Statistics, 1986. DOI: 10.1214/aos/1176350142. ↩

  5. Enno Mammen, “Bootstrap and Wild Bootstrap for High Dimensional Linear Models”, The Annals of Statistics, 1993. DOI: 10.1214/aos/1176349025. ↩

  6. Patrick Kline and Andres Santos, “A Score Based Approach to Wild Bootstrap Inference”, NBER Working Paper 16127, 2010. The paper reviews the wild bootstrap and interprets it through score perturbations. ↩

  7. Hans R. Kunsch, “The Jackknife and the Bootstrap for General Stationary Observations”, The Annals of Statistics, 1989. DOI: 10.1214/aos/1176347265. ↩