Bootstrap and Resampling Methods

In plain language

The bootstrap is a way to measure how uncertain an estimate is — like a sample mean or a cost difference — without needing a complex formula. You take your observed data and repeatedly draw a new sample of the same size from it (allowing the same observation to be picked more than once), compute your statistic in each fake sample, and the spread of those results tells you how much the estimate could vary if you ran the study again. In health economics, the bootstrap is the go-to tool for costs (which are rarely bell-shaped) and for situations where the statistical pipeline is too complex for a formula — as long as you resample whole patients, not individual claim rows, to keep correlated data intact.

What the bootstrap actually does

The bootstrap, introduced by Efron (1979), exploits a simple observation: the observed data are the best available approximation to the population. If you treat the n observations as a surrogate population and draw B re-samples of size n with replacement, each re-sample mimics "what would happen if you ran the study again." The statistic computed on each re-sample is one draw from the bootstrap sampling distribution. Aggregate B such draws and you have an empirical approximation to the sampling distribution of your estimator — obtained purely from the data, without assuming normality, without a closed-form SE formula, and without the Delta method.

Three quantities come out of the bootstrap distribution: the bootstrap standard error (the standard deviation of the B replicate statistics), the bootstrap bias estimate (the mean of the replicates minus the original-sample statistic), and bootstrap confidence intervals in one of several flavors described below. The bias is usually small relative to the SE, but for ICERs and ratio estimators it can matter.

Types of bootstrap confidence intervals

Three variants are in common use, in increasing order of accuracy:

The normal-approximation CI uses the bootstrap SE as a plug-in: estimate ± 1.96 × boot_SE. This is the simplest form but assumes the bootstrap distribution is symmetric and bell-shaped — an assumption violated by many HEOR estimands (costs, ICERs).

The percentile CI reads directly from the tails of the bootstrap distribution: the 2.5th and 97.5th percentiles of the B replicate statistics. It is transformation-equivariant (working on the log scale gives the same result as exponentiating percentile bounds computed on the log scale), does not assume symmetry, and is the dominant method in health economics practice. Its main limitation is poor coverage when the bootstrap distribution is skewed or biased.

The bias-corrected and accelerated (BCa) CI adjusts the percentiles used to extract the CI, correcting for both bootstrap bias and the rate at which the SE changes with the parameter value (acceleration). BCa generally achieves better coverage than the percentile method, particularly for small n or skewed estimands, but requires B ≥ 2000 replicates and jackknife acceleration estimates. For ICER CIs reported to HTA bodies (NICE, ICER), BCa is the defensible choice.

When bootstrap beats formulas in HEOR

Four situations in which the bootstrap is the right tool, not just a convenience:

Mean cost differences from skewed claims data — the canonical HEOR application. Total annual costs per patient follow a right-skewed distribution with a spike at low values and a heavy upper tail driven by high-cost admissions. The analytic SE for the mean difference assumes an approximately normal sampling distribution of the mean, which the CLT protects only at large n. At moderate n (< 500 per arm), or when comparing subgroups, the CLT may not have kicked in and the normal approximation for the SE is unreliable. The bootstrap SE is exact for any sample size and any distribution shape. The bootstrap is so standard here that Barber and Thompson (2000) demonstrated it specifically for randomized cost analyses using claims-style data.

ICER uncertainty — a cluster of cost-effectiveness pairs on the CE plane. The ICER is a ratio of expected costs to expected effects: ICER = ΔCost / ΔEffect. The Delta-method SE for a ratio is unreliable when ΔEffect is close to zero or when costs and effects are both skewed. The bootstrap generates B pairs (ΔCost_b, ΔEffect_b) that trace out the joint sampling distribution on the cost-effectiveness plane. Reading the 95% CI from the percentile ellipse is the standard way to characterize ICER uncertainty for HTA submissions and net monetary benefit (NMB) analyses.

Medians and quantiles — there is no simple analytic SE formula for the median that works at small-to-moderate n with skewed data; the bootstrap is the default. Similarly for the 75th percentile of cost distributions or the 90th percentile of length-of-stay, neither of which has a standard formula under non-normal data.

Complex pipelines where analytic SEs are wrong — propensity-score matching, inverse probability weighting, two-part cost models, and propensity-trimmed cohorts all involve multiple estimation steps. Analytic SEs derived from the final-stage regression ignore the variance introduced by the earlier stages (selecting a caliper, estimating weights, choosing a bandwidth). The full-pipeline cluster bootstrap — re-running the entire pipeline in each replicate — correctly propagates all sources of uncertainty.

The cluster bootstrap: the RWE-critical variant

In claims and EHR data, the unit of analysis is the patient, but the unit of observation is often the claim row: a patient may contribute dozens of pharmacy, outpatient, or inpatient rows. A naive bootstrap that resamples rows treats rows as independent, underestimating the SE because it ignores within-patient correlation. The cluster bootstrap resamples whole patients: draw B samples of n patients with replacement from the patient roster; for each sampled patient, keep all of their rows intact. This preserves the longitudinal correlation structure and gives a valid SE even when patients have very different numbers of claims.

The cluster bootstrap is also the correct approach for matched cohorts, where pairs or clusters of matched patients must be resampled together to maintain balance. It extends naturally to multi-level data (patients nested in hospitals) by resampling the highest-level cluster (hospitals).

Number of replicates: conventions

B = 1000 replicates is the minimum for a reliable percentile CI at the 95% level. B = 2000 is the recommended minimum for BCa CIs (the jackknife acceleration estimate adds variance that additional replicates dampen). B = 5000–10000 is used for small-p tail probabilities, bootstrap hypothesis tests, or when the number of patients is small (< 100) and the bootstrap distribution is discrete. Runtime scales linearly with B; for large claims datasets, the cluster bootstrap is parallelizable by splitting replicate indices across cores.

Permutation tests: the testing cousin

The permutation test addresses the same question — is this difference consistent with chance — using a complementary idea. Rather than resampling with replacement, it randomly shuffles the group labels among the observed units and computes the test statistic in each shuffled dataset. The resulting null distribution is the exact distribution under the null of exchangeability (no treatment effect). The p-value is the fraction of permuted statistics more extreme than the observed. Permutation tests are valid under a weaker assumption than the bootstrap (exchangeability rather than IID sampling) and are exact in finite samples; they do not produce CIs or effect estimates. The two methods are complementary: use the bootstrap when the goal is estimation and uncertainty quantification; use permutation when the primary goal is testing a sharp null with no parametric assumptions.

Failure modes

Tiny n (< 20): the bootstrap distribution is discrete — only a finite number of distinct re-samples are possible — and tail coverage is poor. At n = 5, the exact percentile CI is dominated by the two extreme values in the sample; any additional precision claimed from the bootstrap is illusory. A minimum of n ≈ 30 is a practical threshold for reliable percentile CIs; BCa CIs require still larger n.

Extreme tails and heavy-tailed distributions: when the underlying distribution is Pareto-like (insurance catastrophic claims, rare high-cost events), the bootstrap SE is highly variable across experiments because each re-sample may or may not capture the extreme observation. Coverage of the bootstrap CI for the mean can be poor. In these cases, winsorization of costs above the 99th percentile, or a model-based approach (gamma or lognormal GLM), may be more stable.

Matching with replacement subtleties: in propensity-score matching where controls are matched with replacement, the effective sample size for controls is smaller than the nominal sample size (some controls are used many times). The cluster bootstrap should resample treated patients, then re-do the matching in each replicate — not just resample the already- matched dataset, which would double-count the matching uncertainty.

Time-series and panel data: standard bootstrap is invalid for longitudinally correlated outcomes where the correlation is across time within a patient. The cluster bootstrap (cluster = patient) handles this correctly. If the research question involves time trends, a block bootstrap (sampling blocks of consecutive time-periods) may be needed in addition.

Pros, cons, and trade-offs

Pros: valid for virtually any estimator without distributional assumptions; automatically correct for complex pipelines where analytic SEs are unavailable or wrong; produces full empirical sampling distributions that can be plotted and inspected; the cluster variant naturally handles correlated data; BCa CIs have asymptotically correct coverage under mild regularity conditions; computationally straightforward to parallelize.

Cons: computationally intensive (B × pipeline runtime); fails at very small n or under very heavy tails; does not substitute for a valid identification strategy — the bootstrap quantifies sampling uncertainty but cannot correct for unmeasured confounding or selection bias; naive implementation (resampling rows instead of patients) in claims data gives anti-conservative SEs; results can appear to differ across runs if the random seed is not fixed.

When to use

Use the bootstrap when: (a) the estimator has no tractable analytic SE (ICERs, medians, complex pipeline outputs); (b) the outcome is right-skewed and n is insufficient for the CLT to protect the normal approximation for the mean (cost analyses at n < 500 per arm, subgroup analyses); (c) the data have a clustered structure (claims, EHR rows nested within patients) and the cluster bootstrap is the correct SE estimator; (d) the analysis involves PS-matching, IPTW, or other multi-stage procedures where the SE must propagate uncertainty through all stages; (e) a full sampling distribution on the cost-effectiveness plane is needed for probabilistic sensitivity analysis.

When NOT to use

Do not use the bootstrap as the primary method when: (a) n is very small (< 20–30 per cluster) — the discrete bootstrap distribution gives unreliable tail coverage; (b) the estimand is a mean from a large balanced RCT where the CLT is fully operational and an analytic t-test CI is both valid and more efficient; (c) the naive implementation resamples rows rather than patients in claims data — this is a common error that produces CIs that are too narrow; (d) the goal is a sharp null test without the need for CIs — a permutation test is better powered and more interpretable in that case; (e) the pipeline does not re-run all estimation steps (PS model, matching, outcome model) within each replicate — partial bootstrapping understates the true uncertainty.

Interpreting the output

The worked example below uses five patients with annual costs [1200, 3400, 800, 5600, 2100] USD (n = 5, illustrative only). The original sample mean is 2620 USD. Three explicit re-samples yield bootstrap replicate means of 2180, 3060, and 2880 USD. With B = 2000 replicates run in software, the bootstrap SE ≈ 860 USD and the percentile 95% CI ≈ [1060, 4740] USD.

(1) Formal interpretation. The bootstrap percentile 95% CI [1060, 4740] is a data-derived interval: in 95% of hypothetical repeated studies of size n = 5 drawn from the same population, an interval constructed by the same percentile procedure would contain the true population mean cost. The interval is not symmetric around the point estimate (2620 USD) — it extends further to the right than to the left — because the bootstrap distribution inherits the right-skew of the underlying cost distribution, which pulls replicate means upward when the high-cost patient (5600 USD) is over-represented in a re-sample.

(2) Practical interpretation. The mean annual cost is 2620 USD with a 95% confidence interval of roughly 1060 to 4740 USD. This wide range reflects genuine uncertainty from a small sample; a decision-maker should treat the point estimate as uncertain by a factor of roughly two in either direction. In a real analysis with n in the hundreds per arm, this interval would be far narrower, but the bootstrap procedure — including the cluster variant that resamples patients rather than individual claims rows — remains the correct method for skewed cost data regardless of sample size.

Worked example

Scenario

A health economist is estimating the mean annual total cost for five patients with a rare condition in a pilot claims analysis. The five observed costs are right-skewed (one patient has very high costs). She applies the bootstrap to quantify uncertainty around the mean without assuming normality. She draws three explicit re-samples by hand to illustrate the method, then would run B = 2000 replicates in software for the final CI.

Dataset

Annual total cost per patient (USD) from five pilot patients. Costs are right-skewed; patient P04 is a high-cost outlier. These are the five observations that will be resampled with replacement.

Steps

• Compute the original sample mean. Sum of the five costs: 1200+3400+800+5600+2100 = 13100 dollars. Original mean = 13100/5 = 2620 dollars.

• Draw Resample A with replacement (five draws, order determined by random indices). One possible draw: [800, 3400, 1200, 3400, 2100]. P02 and P03 appear once each, P01 once, P02 again, P05 once; P04 is absent. Sum = 800+3400+1200+3400+2100 = 10900 dollars. Resample A mean = 10900/5 = 2180 dollars. This resample missed the high-cost patient, so its mean is below the original.

• Draw Resample B with replacement. One possible draw: [5600, 1200, 5600, 800, 2100]. P04 appears twice (sampled on two draws), P01 once, P03 once, P05 once. Sum = 5600+1200+5600+800+2100 = 15300 dollars. Resample B mean = 15300/5 = 3060 dollars. This resample overrepresents the high-cost patient, pulling the mean above the original.

• Draw Resample C with replacement. One possible draw: [2100, 2100, 3400, 1200, 5600]. P05 appears twice, P02 once, P01 once, P04 once. Sum = 2100+2100+3400+1200+5600 = 14400 dollars. Resample C mean = 14400/5 = 2880 dollars.

• The three replicate means are 2180, 3060, and 2880 USD. Their spread shows that the original estimate of 2620 USD is uncertain. With B = 2000 replicates run in software, the standard deviation of all replicate means — the bootstrap SE — is approximately 860 USD, and the percentile 95% CI spans approximately [1060, 4740] USD.

• In a real claims analysis, the cluster bootstrap would resample patient IDs rather than individual claim rows, so all claims for a drawn patient travel together and within-patient correlation is preserved.

Result

Original mean = 13100/5 = 2620 USD. Three explicit bootstrap replicate means: 10900/5 = 2180 USD (missed high-cost patient), 15300/5 = 3060 USD (over-represented high-cost patient), 14400/5 = 2880 USD. The variability across resamples — from 2180 to 3060 USD — illustrates that the true mean is uncertain by roughly plus or minus 440 USD even across just three replicates; with B = 2000 replicates the bootstrap SE is approximately 860 USD and the percentile 95% CI is approximately [1060, 4740] USD.

Runnable example

import numpy as np
from scipy import stats

rng = np.random.default_rng(seed=42)

# ── Motivating dataset: five annual costs from worked example ──
costs = np.array([1200, 3400, 800, 5600, 2100], dtype=float)
n = len(costs)
original_mean = costs.mean()
print(f"Original sample mean: {original_mean:.2f} USD")

# ── 1. Manual bootstrap loop for SE and percentile CI ──
B = 2000
boot_means = np.empty(B)
for b in range(B):
    sample = rng.choice(costs, size=n, replace=True)
    boot_means[b] = sample.mean()

boot_se = boot_means.std(ddof=1)
pct_ci = np.percentile(boot_means, [2.5, 97.5])
print(f"Bootstrap SE: {boot_se:.2f}")
print(f"Percentile 95% CI: [{pct_ci[0]:.2f}, {pct_ci[1]:.2f}]")

# ── 2. BCa CI via scipy.stats.bootstrap ──
result = stats.bootstrap(
    (costs,),
    statistic=np.mean,
    n_resamples=2000,
    confidence_level=0.95,
    method="BCa",
    random_state=42,
)
print(f"BCa 95% CI: [{result.confidence_interval.low:.2f}, {result.confidence_interval.high:.2f}]")
print(f"Bootstrap SE (scipy): {result.standard_error:.2f}")

# ── 3. Cluster bootstrap for correlated claims data ──
# Simulated patient-level claims dataset: multiple rows per patient
import pandas as pd

claims = pd.DataFrame({
    "patient_id": [1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 5, 5],
    "claim_cost":  [400, 300, 500, 1200, 2200, 250, 250, 150, 150, 5600, 900, 1200],
})

# Build per-patient totals first (these are the analysis unit)
patient_totals = claims.groupby("patient_id")["claim_cost"].sum().values
patient_ids = claims["patient_id"].unique()
n_patients = len(patient_ids)

# Cluster bootstrap: resample patient IDs, keep all rows for each patient
B_clust = 2000
clust_means = np.empty(B_clust)
for b in range(B_clust):
    sampled_ids = rng.choice(patient_ids, size=n_patients, replace=True)
    # Collect all rows for sampled patients (with repetition)
    boot_claims = pd.concat(
        [claims[claims["patient_id"] == pid] for pid in sampled_ids],
        ignore_index=True,
    )
    # Re-aggregate: for each draw of a patient, their cost appears once per draw
    # Use a running index to keep duplicate patients separate
    boot_costs = [
        claims.loc[claims["patient_id"] == pid, "claim_cost"].sum()
        for pid in sampled_ids
    ]
    clust_means[b] = np.mean(boot_costs)

clust_pct_ci = np.percentile(clust_means, [2.5, 97.5])
print(f"\nCluster bootstrap mean: {np.mean(clust_means):.2f}")
print(f"Cluster bootstrap SE: {clust_means.std(ddof=1):.2f}")
print(f"Cluster percentile 95% CI: [{clust_pct_ci[0]:.2f}, {clust_pct_ci[1]:.2f}]")
print("Note: cluster bootstrap resamples patients, not rows — correct for correlated claims.")

# ── 4. Permutation test for two-group mean difference ──
group_a = np.array([1200, 3400, 800, 5600, 2100], dtype=float)
group_b = np.array([900, 1100, 2200, 3000, 1500], dtype=float)
observed_diff = group_a.mean() - group_b.mean()
combined = np.concatenate([group_a, group_b])
n_a = len(group_a)

perm_diffs = np.empty(B)
for b in range(B):
    perm = rng.permutation(combined)
    perm_diffs[b] = perm[:n_a].mean() - perm[n_a:].mean()

p_perm = np.mean(np.abs(perm_diffs) >= np.abs(observed_diff))
print(f"\nPermutation test: observed diff = {observed_diff:.2f}, p = {p_perm:.4f}")
print("Permutation test produces a p-value only (no CI); use bootstrap for estimation.")

library(boot)

# ── Motivating dataset ──
costs <- c(1200, 3400, 800, 5600, 2100)
n <- length(costs)
cat("Original mean:", mean(costs), "USD\n")

# ── 1. Manual bootstrap loop ──
set.seed(42)
B <- 2000
boot_means <- replicate(B, mean(sample(costs, size = n, replace = TRUE)))
cat(sprintf("Bootstrap SE: %.2f\n", sd(boot_means)))
cat(sprintf("Percentile 95%% CI: [%.2f, %.2f]\n",
            quantile(boot_means, 0.025), quantile(boot_means, 0.975)))

# ── 2. boot package: percentile and BCa ──
mean_fn <- function(data, indices) mean(data[indices])
boot_obj <- boot(data = costs, statistic = mean_fn, R = B)

pct_ci <- boot.ci(boot_obj, type = "perc")
bca_ci <- boot.ci(boot_obj, type = "bca")
cat(sprintf("Percentile CI: [%.2f, %.2f]\n",
            pct_ci$percent[4], pct_ci$percent[5]))
cat(sprintf("BCa CI:        [%.2f, %.2f]\n",
            bca_ci$bca[4], bca_ci$bca[5]))
cat("BCa adjusts for bootstrap bias and acceleration (prefer for skewed estimands).\n")

# ── 3. Cluster bootstrap by patient_id ──
claims <- data.frame(
  patient_id = c(1, 1, 1, 2, 2, 3, 3, 3, 3, 4, 5, 5),
  claim_cost = c(400, 300, 500, 1200, 2200, 250, 250, 150, 150, 5600, 900, 1200)
)

cluster_bootstrap <- function(data, B = 2000, seed = 42) {
  set.seed(seed)
  patient_ids  <- unique(data$patient_id)
  n_patients   <- length(patient_ids)
  # Per-patient totals (these are the analysis units)
  totals <- tapply(data$claim_cost, data$patient_id, sum)

  replicate(B, {
    drawn_ids  <- sample(patient_ids, size = n_patients, replace = TRUE)
    mean(totals[as.character(drawn_ids)])
  })
}

clust_means <- cluster_bootstrap(claims)
cat(sprintf("\nCluster bootstrap SE: %.2f\n", sd(clust_means)))
cat(sprintf("Cluster percentile 95%% CI: [%.2f, %.2f]\n",
            quantile(clust_means, 0.025), quantile(clust_means, 0.975)))
cat("Cluster bootstrap: resamples patients not rows — correct for correlated claims data.\n")

# ── 4. Permutation test ──
group_a <- c(1200, 3400, 800, 5600, 2100)
group_b <- c(900, 1100, 2200, 3000, 1500)
obs_diff <- mean(group_a) - mean(group_b)
combined <- c(group_a, group_b)
n_a <- length(group_a)

set.seed(42)
perm_diffs <- replicate(B, {
  perm <- sample(combined)
  mean(perm[seq_len(n_a)]) - mean(perm[(n_a + 1):length(perm)])
})
p_perm <- mean(abs(perm_diffs) >= abs(obs_diff))
cat(sprintf("\nPermutation p-value: %.4f (observed diff = %.2f)\n", p_perm, obs_diff))
cat("Permutation test gives a p-value only; use boot() for a CI on the mean difference.\n")

# ── 5. Bootstrap for ICER uncertainty ──
# Patient-level incremental cost and incremental effect (illustrative)
set.seed(42)
n_pts <- 50
delta_cost   <- rnorm(n_pts, mean = 5000, sd = 3000)
delta_effect <- rnorm(n_pts, mean = 0.10,  sd = 0.08)
icer_fn <- function(data, idx) {
  mean(data[idx, 1]) / mean(data[idx, 2])
}
icer_data <- cbind(delta_cost, delta_effect)
icer_boot <- boot(data = icer_data, statistic = icer_fn, R = B)
icer_pct  <- boot.ci(icer_boot, type = "perc")
cat(sprintf("\nICER bootstrap percentile 95%% CI: [%.0f, %.0f] USD/QALY\n",
            icer_pct$percent[4], icer_pct$percent[5]))