Bootstrap and Resampling Methods

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 NOT to use

Do not use the bootstrap as the primary method when:

1. n is very small (< 20–30 per cluster) — the discrete bootstrap distribution gives unreliable tail coverage;
2. 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;
3. the naive implementation resamples rows rather than patients in claims data — this is a common error that produces CIs that are too narrow;
4. the goal is a sharp null test without the need for CIs — a permutation test is better powered and more interpretable in that case;
5. 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.

Decision diagram

flowchart TD
  A[Original data<br/>n observations] --> B[Draw B resamples<br/>with replacement]
  B --> C1[Resample 1 → stat_1]
  B --> C2[Resample 2 → stat_2]
  B --> C3[... → ...]
  B --> CB[Resample B → stat_B]
  C1 --> D[Bootstrap distribution<br/>of stat_1 ... stat_B]
  C2 --> D
  C3 --> D
  CB --> D
  D --> E1[Bootstrap SE<br/>= SD of replicates]
  D --> E2[Percentile CI<br/>= 2.5th, 97.5th pct]
  D --> E3[BCa CI<br/>= adjusted percentiles]
  A --> F{Clustered data?}
  F -->|Yes - claims or EHR| G[Resample PATIENTS<br/>not rows<br/>= Cluster bootstrap]
  F -->|No - already patient-level| B
  G --> B

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

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.

Trade-offs

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]))

/* ── 1. Create the patient-level cost dataset ── */
data work.patients;
  input patient_id $ annual_cost_usd;
  datalines;
P01 1200
P02 3400
P03  800
P04 5600
P05 2100
;
run;

/* ── 2. Standard bootstrap: resample patient rows with replacement ──
   PROC SURVEYSELECT METHOD=URS draws n units with replacement in each replicate.
   REPS= sets the number of bootstrap replicates.
   SEED= fixes the random seed for reproducibility.                               */
proc surveyselect data=work.patients
      method=urs         /* unrestricted random sampling = with replacement      */
      n=5                /* resample the same number of units as the original    */
      reps=2000          /* B bootstrap replicates                               */
      seed=20240101
      out=work.boot_sample
      outall;            /* keep original sample in replicate 0                  */
run;
/* Each replicate is identified by the variable Replicate (1 to 2000).           */

/* ── 3. Compute mean cost in each replicate ── */
proc means data=work.boot_sample noprint;
  by Replicate;
  var annual_cost_usd;
  output out=work.boot_means (drop=_type_ _freq_)
         mean=boot_mean;
run;
/* Remove replicate 0 (original sample) for the bootstrap distribution           */
data work.boot_dist;
  set work.boot_means;
  where Replicate >= 1;
run;

/* ── 4. Read percentile CI from the bootstrap distribution ── */
proc univariate data=work.boot_dist noprint;
  var boot_mean;
  output out=work.boot_ci
         mean=boot_mean_of_means
         std=boot_se
         pctlpts=2.5 97.5
         pctlpre=pct_;
run;
proc print data=work.boot_ci; run;
/* pct_2_5 = lower bound; pct_97_5 = upper bound; boot_se = bootstrap SE        */

/* ── 5. Cluster bootstrap: resample patients, keep all claim rows per patient ──
   When data are at claim level (multiple rows per patient), cluster by patient_id
   using PROC SURVEYSELECT with the CLUSTER statement.                            */
data work.claims;
  input patient_id $ claim_cost;
  datalines;
P01  400
P01  300
P01  500
P02 1200
P02 2200
P03  250
P03  250
P03  150
P03  150
P04 5600
P05  900
P05 1200
;
run;

/* Sort by cluster variable before SURVEYSELECT                                  */
proc sort data=work.claims; by patient_id; run;

proc surveyselect data=work.claims
      method=urs
      n=5               /* number of clusters (patients) to resample            */
      reps=2000
      seed=20240101
      out=work.clust_boot
      outall;
  cluster patient_id;   /* resample whole patients, keeping all their rows      */
run;
/* Each drawn patient's rows are duplicated in the output dataset when that
   patient was sampled multiple times.                                            */

/* Compute per-patient totals in each replicate, then mean across patients       */
proc means data=work.clust_boot noprint nway;
  by Replicate;
  class patient_id SelectionProb;   /* patient_id + replicate-specific weight   */
  var claim_cost;
  output out=work.pt_totals sum=patient_total;
run;
proc means data=work.pt_totals noprint;
  by Replicate;
  var patient_total;
  output out=work.clust_means mean=clust_mean;
run;
data work.clust_dist;
  set work.clust_means;
  where Replicate >= 1;
run;
proc univariate data=work.clust_dist noprint;
  var clust_mean;
  output out=work.clust_ci
         std=clust_se
         pctlpts=2.5 97.5
         pctlpre=cpct_;
run;
proc print data=work.clust_ci;
  title "Cluster bootstrap 95% CI for mean per-patient cost (USD)";
run;
/* cpct_2_5 and cpct_97_5 are the percentile CI bounds.
   This correctly handles within-patient claim correlation.                       */

Citations

1. [1]Efron B. Bootstrap methods: another look at the jackknife. Annals of Statistics. 1979;7(1):1-26.
2. [2]Carpenter J, Bithell J. Bootstrap confidence intervals: when, which, what? A practical guide for medical statisticians. Statistics in Medicine. 2000;19(9):1141-1164.

1. [3]Barber JA, Thompson SG. Analysis of cost data in randomized trials: an application of the non-parametric bootstrap. Statistics in Medicine. 2000;19(23):3219-3236.