The bootstrap explained: standard errors by resampling

Every estimate needs a measure of its precision. For a mean or an OLS coefficient there are formulas. For a median, a ratio or a regression whose error variance is not constant, the formulas are complicated, rest on doubtful assumptions or do not exist. The bootstrap, introduced by Efron (1979), replaces the formula with computation: it resamples the data you have to imitate the samples you could have drawn.

The plug-in idea

A standard error is the standard deviation of an estimate across repeated samples from the population. We have only one sample, and we cannot draw more.

The bootstrap plugs in the best available picture of the population: the sample itself. Its empirical distribution puts probability 1/n on each of the n observations, so drawing from it means picking rows of your data at random, with replacement. Without replacement, every resample would be the original data in a different order and the statistic would never change. With replacement, some rows appear twice or more and others not at all, so the statistic varies from resample to resample, roughly as it varies across real samples.

The recipe:

1. Draw n rows with replacement from the data.
2. Compute the statistic on that resample, θ̂*_b.
3. Repeat B times.
4. The standard deviation of the B values is the bootstrap standard error. Their 2.5th and 97.5th percentiles form a 95% percentile interval.

SE_boot = √( Σ_b (θ̂*_b − θ̄*)² / (B − 1) ),   where θ̄* = (1/B) Σ_b θ̂*_b

The statistics lab Resampling the observed baskets works through this recipe for a mean, where a formula can check the answer. This article applies it where formulas are harder to come by.

A worked example: median income

We simulate 200 people with years of schooling and annual income in thousands. Given schooling, income is lognormal, so a few high earners stretch the right tail. The data are synthetic.

Simulated annual incomes of 200 people are skewed to the right

The median is a natural summary for skewed data, but its usual standard error formula needs the density of incomes at the median, which is hard to estimate. Efron (1979) used the median as one of his first examples. The bootstrap needs only numpy:

import numpy as np


def simulate(rng, n=200):
    """Synthetic years of schooling and annual income in thousands."""
    schooling = np.clip(np.round(rng.normal(12, 3, n)), 4, 20)
    income = np.exp(1.8 + 0.15 * schooling + rng.normal(0, 0.5, n))
    return schooling, income


schooling, income = simulate(np.random.default_rng(1979))

boot = np.random.default_rng(1)
rows = boot.integers(0, 200, size=(4000, 200))  # 4,000 resamples of row numbers
medians = np.median(income[rows], axis=1)
low, high = np.percentile(medians, [2.5, 97.5])
print(f"median {np.median(income):.2f}, mean {income.mean():.2f}")
print(f"bootstrap SE {medians.std(ddof=1):.2f}, 95% percentile interval {low:.2f} to {high:.2f}")

median 35.08, mean 44.40
bootstrap SE 2.77, 95% percentile interval 31.08 to 42.67

Medians of 4,000 bootstrap resamples of the same 200 incomes

Two features stand out. The histogram is lumpy, because the median of a resample can only be an observed income or the average of two of them. And the interval is asymmetric: it reaches 4.00 below the sample median and 7.59 above it, because incomes just above the middle are more spread out than those just below. A symmetric interval, estimate ± 1.96 × SE, would hide this.

Bootstrapping a regression: the pairs bootstrap

For a regression, resample whole rows, so that each person's income travels with their schooling. This is the pairs bootstrap. Freedman (1981) contrasted it with resampling residuals around fixed regressors. The residual bootstrap needs errors with a common distribution; the pairs bootstrap treats rows as independent draws from a population and stays valid when the error variance depends on the regressors.

Here the spread of income grows with schooling, so the errors are heteroskedastic. OpenEconometrics runs the pairs bootstrap with covariance="bootstrap". The output below is from openecon 0.3.20a1:

import openecon as oe

data = {"income": income, "schooling": schooling}
for covariance, options in [("nonrobust", {}), ("HC3", {}),
                            ("bootstrap", {"reps": 999, "seed": 1})]:
    fit = oe.ols(data=data, y="income", x=["schooling"], covariance=covariance, **options)
    slope = fit.coefficients[1]
    print(f"{covariance:9} slope {slope.estimate:.2f}  SE {slope.std_error:.3f}  "
          f"95% CI {slope.ci_low:.2f} to {slope.ci_high:.2f}")
keys = ["scheme", "reps", "successful_reps", "failed_reps", "seed",
        "random_engine", "random_count_algorithm", "distribution"]
print({key: fit.inference[key] for key in keys})

nonrobust slope 6.52  SE 0.610  95% CI 5.32 to 7.72
HC3       slope 6.52  SE 0.748  95% CI 5.04 to 8.00
bootstrap slope 6.52  SE 0.729  95% CI 5.09 to 7.95
{'scheme': 'pairs', 'reps': 999, 'successful_reps': 999, 'failed_reps': 0, 'seed': 1, 'random_engine': 'torch.Generator', 'random_count_algorithm': 'torch.randint with bincount', 'distribution': 'normal'}

The slope is the same in all three fits; only its uncertainty changes. The classical standard error assumes a constant error variance and is the smallest. The HC3 and bootstrap standard errors, which both allow the variance to change with schooling, are within 0.02 of each other.

The last line is what the fit records about its draws. scheme says whole rows were resampled. successful_reps and failed_reps count resamples that did or did not give a full-rank fit; failed ones are excluded and reported. random_engine and random_count_algorithm name the generator and the method that drew how many times each row appears. distribution is a reminder that this option estimates a standard error and then uses z statistics and a normal interval, estimate ± 1.96 × SE. For a percentile interval, pass a fit with an analytical covariance to oe.bootstrap(fit, data, reps=999, seed=1, ci="percentile"), which refits the model on every resample and stores the limits in extra["bootstrap_ci"]. For the classical fit of this slope it gives 5.13 to 8.12.

Three standard errors for the same schooling slope

Checking against the truth

Because the data are simulated, we can draw new samples from the same population and measure the real sampling variation:

fresh = np.random.default_rng(99)
samples = [simulate(fresh) for _ in range(4000)]
print(f"median SD {np.std([np.median(y) for _, y in samples], ddof=1):.2f}")
print(f"slope SD {np.std([np.polyfit(x, y, 1)[0] for x, y in samples], ddof=1):.3f}")

median SD 2.17
slope SD 0.871

The classical standard error, 0.610, understates the true variability of the slope by 30%. HC3 and the bootstrap come much closer, although both are still below 0.871: a standard error computed from one sample of 200 is itself an estimate.

The median makes the same point. Its bootstrap standard error, 2.77, is above the true 2.17. When we repeated the exercise on 1,000 fresh samples, with 1,000 resamples each, the bootstrap standard error averaged 2.24, close to the truth, but its 10th and 90th percentiles were 1.72 and 2.82. Our sample sits near the top of that range. The bootstrap is right on average, not in every sample.

Grouped data: the cluster bootstrap

Resampling rows assumes the rows are independent. Households in the same village share prices, weather and labour markets, so their incomes move together. The cluster bootstrap resamples whole villages with replacement and keeps each village's households together. In oe.ols, add cluster= to covariance="bootstrap". In this simulation, both schooling and income have a village component:

def simulate_villages(rng, villages=40, size=5):
    """Five households in each of 40 villages; each village shares an income shock."""
    village = np.repeat(np.arange(villages), size)
    schooling = np.clip(np.round(np.repeat(rng.normal(12, 2, villages), size)
                                 + rng.normal(0, 2, villages * size)), 4, 20)
    shock = np.repeat(rng.normal(0, 0.4, villages), size)
    income = np.exp(1.8 + 0.15 * schooling + shock + rng.normal(0, 0.4, villages * size))
    return {"income": income, "schooling": schooling, "village": village}


grouped = simulate_villages(np.random.default_rng(40))
choices = {
    "HC3": {"covariance": "HC3"},
    "CR1 cluster": {"covariance": "cluster", "cluster": "village"},
    "pairs bootstrap": {"covariance": "bootstrap", "reps": 999, "seed": 1},
    "cluster bootstrap": {"covariance": "bootstrap", "cluster": "village", "reps": 999, "seed": 1},
}
for label, options in choices.items():
    fit = oe.ols(data=grouped, y="income", x=["schooling"], **options)
    print(f"{label:17} SE {fit.coefficients[1].std_error:.3f}")
print({key: fit.inference[key] for key in ("scheme", "cluster_count")})

HC3               SE 0.716
CR1 cluster       SE 0.950
pairs bootstrap   SE 0.683
cluster bootstrap SE 0.959
{'scheme': 'cluster', 'cluster_count': 40}

The pairs bootstrap agrees with HC3, and the cluster bootstrap with the cluster-robust CR1 standard error. The two cluster-based numbers are about a third larger, because they account for the shared village shock. A bootstrap is only as good as its resampling unit, which should be the unit that was sampled independently.

Few clusters are a harder case. A resample contains only about two-thirds of the distinct clusters, some of them repeated, and with a treatment assigned by cluster, some resamples contain no treated or no untreated clusters at all. Cameron, Gelbach and Miller (2008) show that conventional cluster-robust tests can over-reject with five to thirty clusters, and that cluster bootstrap-t procedures, which bootstrap the t statistic rather than only the standard error, bring rejection rates close to the nominal level. The wild cluster bootstrap is the standard tool for that situation.

When the bootstrap fails or needs care

The bootstrap is justified by large-sample arguments. MacKinnon (2006) notes that well-chosen bootstrap methods generally work very well in regressions with independent, identically distributed errors, but not always in other cases, such as models with dependent errors.

• Very small samples. Ten observations are a poor picture of a skewed population. In 2,000 simulated samples of 10 incomes from the population above, the 95% percentile interval for the mean covered the population mean in only 86.2% of samples, against 90.0% for the t interval. With 200 incomes, it covered it in 93.7%.
• Heavy tails. When the variance is infinite, the bootstrap of the mean is inconsistent: its distribution stays random however large the sample becomes (Athreya 1987).
• Dependence. Resampling single observations breaks serial and within-group correlation and understates uncertainty. Resample clusters, or blocks of consecutive periods for time series; oe.bootstrap offers scheme="moving_block" for OLS with a regular time column.
• Parameters at a boundary. When the true parameter lies on the edge of the parameter space, such as a mean restricted to be non-negative whose true value is zero, the bootstrap is inconsistent (Andrews 2000).
• A bootstrap standard error is not a refinement. Used with a normal critical value, as covariance="bootstrap" does, it has the same large-sample accuracy as an analytical robust standard error. Bootstrapping the t statistic itself can give more accurate tests (Cameron, Gelbach and Miller 2008).

Reproducibility: seeds and replications

Two kinds of randomness are at work. Sampling variation comes from the data, and the standard error measures it. Monte Carlo variation comes from using B random resamples instead of all possible ones. A seed fixes the second, so the same call returns the same answer.

• Report B and the seed. For the slope above, seeds 1 to 5 give bootstrap standard errors between 0.710 and 0.742 with 999 resamples, and between 0.693 and 0.762 with 199. If you omit both, oe.ols uses 199 resamples and seed 0, and records them.
• Use more draws for intervals than for standard errors. A percentile depends on the few draws in each tail.
• A seed belongs to a generator. numpy's default_rng(1) and the torch generator that openecon records as random_engine produce different streams, and oe.bootstrap draws differently from covariance="bootstrap". Name the function as well as the seed.

A short checklist

• Resample the unit that was sampled independently: rows, clusters or blocks.
• Look at the bootstrap distribution, not only its standard deviation.
• Compare with an analytical robust standard error; a large gap is information about the data.
• Take care with small samples, heavy tails, parameters near a boundary and few clusters.
• Report the scheme, B, the seed and the software version. See p-values in economics and the guide to robust and clustered standard errors for reporting the resulting inference.

References

• Andrews, D. W. K. (2000). Inconsistency of the bootstrap when a parameter is on the boundary of the parameter space. Econometrica, 68(2), 399–405. doi:10.1111/1468-0262.00114
• Athreya, K. B. (1987). Bootstrap of the mean in the infinite variance case. Annals of Statistics, 15(2), 724–731. doi:10.1214/aos/1176350371
• Cameron, A. C., Gelbach, J. B. and Miller, D. L. (2008). Bootstrap-based improvements for inference with clustered errors. Review of Economics and Statistics, 90(3), 414–427. doi:10.1162/rest.90.3.414
• Efron, B. (1979). Bootstrap methods: Another look at the jackknife. Annals of Statistics, 7(1), 1–26. doi:10.1214/aos/1176344552
• Freedman, D. A. (1981). Bootstrapping regression models. Annals of Statistics, 9(6), 1218–1228. doi:10.1214/aos/1176345638
• MacKinnon, J. G. (2006). Bootstrap methods in econometrics. Economic Record, 82(s1), S2–S18. doi:10.1111/j.1475-4932.2006.00328.x