Bootstrap and resampling for statistical inference

What is resampled, what is refitted, and when the approximation works

Before reading: Means, variance, and standard errors. Large-sample reasoning for methods research supplies the approximation theory; Influence functions supports the multiplier section. Use Dependent observations and cluster-aware inference for the cluster branch.

1 A resampling distribution needs a target

A standard error describes variation under repeated sampling from a population. A bootstrap approximates that variation by generating new datasets from a fitted or empirical distribution and recomputing a stated statistic. Specify the target, estimator, independent unit, and sampling design first. More bootstrap replicates reduce numerical error in that approximation; they do not establish its validity.

For i.i.d. observations W_1,\ldots,W_n, the empirical law P_n puts mass 1/n on each observation. The ordinary nonparametric bootstrap draws n observations independently from P_n, with replacement. Resample each person’s complete vector, retaining associations among their features and covariates. Sampling each feature column separately would replace the observed joint law with another model.

“With replacement” means an observed record can appear several times in one resampled dataset, while another record can be absent. For observed values (1,3,5), one possible bootstrap sample is (3,3,5). Both datasets have size 3, but their means are 3 and 11/3. Repeating this draw-and-recompute operation produces a distribution of bootstrap means. The replicate count B is the number of resampled datasets, each still containing n observations.

There are two levels of randomness. First, the original study draws data from the unknown population. Second, with that dataset fixed, we generate resamples from its empirical law. A new original study would supply a different empirical law and a different bootstrap distribution. A star on \Pr_*, E_*, or \operatorname{Var}_* means probability, expectation, or variance over the second level, conditional on the original observations.

A parametric bootstrap instead simulates from a fitted probability model. Its justification includes that model. A permutation test calibrates a statistic through a null assignment or exchangeability mechanism; it is a different construction. See Multivariate comparisons: permutation and distance-based tests for the distinction.

2 Center the bootstrap error correctly

For one draw Y^* from the empirical distribution, E_*[f(Y^*)]=E_{P_n}(f(Y))=P_nf, using the empirical-expectation notation from Large-sample reasoning for methods research.

For an estimator \widehat\theta, let \widehat\theta^* denote its bootstrap version. A usual first-order goal is for the conditional distribution of the bootstrap error to approximate the sampling distribution of the original error, with both on the same scale:

\begin{aligned} \text{Bootstrap: }&\sqrt n(\widehat\theta^*-\widehat\theta),\\ \text{Sampling: }&\sqrt n(\widehat\theta-\theta_0). \end{aligned}

This is a claim about distributions as n increases, not equality in one finite sample. The bootstrap center is the observed estimate because the bootstrap population is fitted from that dataset. Some bootstrap methods approximate finite-sample bias using E_*[\widehat\theta^*]-\widehat\theta under appropriate conditions. That calculation does not establish that the estimator’s population limit equals the scientific target.

For the population mean \mu=E[Y], assume i.i.d. observations with 0<\sigma^2=\operatorname{Var}(Y)<\infty. The original error is \bar Y-\mu. Inside the bootstrap world, the empirical population has mean \bar Y, so the corresponding error is \bar Y^*-\bar Y. Subtracting the same unknown \mu from bootstrap estimates would not reproduce this fitted-population estimation error.

For the sample mean, let Y^* denote one bootstrap draw. Each observed value has probability 1/n, giving

\begin{aligned} E_*[Y^*]&=\bar Y,\\ \operatorname{Var}_*(Y^*)&=\frac1n\sum_i(Y_i-\bar Y)^2. \end{aligned}

The n draws in a bootstrap sample are independent conditional on the original data. Averaging them leaves the conditional mean at \bar Y and divides their conditional variance by n:

\begin{aligned} \operatorname{Var}_*(\bar Y^*)&=\frac{1}{n^2}\sum_i(Y_i-\bar Y)^2\\ &=\frac{n-1}{n}\frac{s_Y^2}{n}. \end{aligned}

Thus E_*[\bar Y^*]=\bar Y, and the conditional variance is asymptotically equivalent to the sampling variance \sigma^2/n when s_Y^2\xrightarrow{p}\sigma^2. This compares variation of \bar Y^* around \bar Y with the data fixed, and variation of \bar Y around \mu across original studies.

Variance matching is one check. Bootstrap validity also needs the distributional approximation. For the i.i.d. finite, nonzero variance mean, a conditional CLT gives, at each fixed t,

\Pr_*\{\sqrt n(\bar Y^*-\bar Y)\le t\} \xrightarrow{p}\Phi(t/\sigma).

The left side is itself a random quantity across original datasets. Convergence in probability says the fitted conditional CDF approaches the same normal CDF as the original scaled mean error. Here \Phi is the standard normal CDF. General statistics need their own smoothness, moment, and remainder arguments; the mean calculation does not establish validity for every resampling procedure.

3 Run a small example with a known reference

The following synthetic sample comes from a normal population with mean 2 and SD 1.5. The true sampling law of the scaled mean error is normal. Compare that law with the conditional bootstrap approximation before looking at the intervals.

set.seed(3101)
n <- 80; B <- 1999; y <- rnorm(n, mean = 2, sd = 1.5)
mean_hat <- mean(y)
boot_mean <- replicate(B, mean(y[sample.int(n, n, replace = TRUE)]))
root <- sqrt(n) * (boot_mean - mean_hat)
root_density <- density(root)
plot(root_density, xlim = c(-5, 5),
     ylim = c(0, 1.08*max(root_density$y, dnorm(0, sd = 1.5))),
     col = "#2F80C8", lwd = 2,
     main = "", xlab = "Scaled mean error", ylab = "Density")
curve(dnorm(x, sd = 1.5), add = TRUE, col = "#D55E00", lty = 2, lwd = 2)
legend("topright", c("Conditional bootstrap", "Known sampling law"),
       col = c("#2F80C8", "#D55E00"), lty = c(1, 2), bty = "n")
c(influence_SE = sd(y)/sqrt(n), bootstrap_SE = sd(boot_mean),
  exact_conditional_bootstrap_SE = sqrt((n-1)/n) * sd(y)/sqrt(n))

                  influence_SE                   bootstrap_SE 
                     0.1758772                      0.1693985 
exact_conditional_bootstrap_SE 
                     0.1747745 

Figure 1: One synthetic normal sample (n=80), with 1,999 ordinary bootstrap replicates. The bootstrap errors are centered at the observed mean; the reference is the known repeated-sampling law of sqrt(n)(sample mean − 2). Differences include the realized sample and bootstrap Monte Carlo error.

In the output, sd(boot_mean) estimates \sqrt{\operatorname{Var}_*(\bar Y^*)} using B simulations. The exact conditional expression follows from the algebra above; the influence-based SE estimates \sigma/\sqrt n using the original study. The two population-level quantities are close at first order, and the simulated SE also has Monte Carlo error. A new original dataset would give another approximation.

An exact normal-population t interval is also available in this particular example; the bootstrap is used to teach the construction. Its histogram is a sampling-distribution approximation, rather than a posterior distribution for \mu.

4 Interval recipes make different transformations

Let q^*_{a} be an a-quantile of the bootstrap estimates. At nominal 95% coverage:

Construction	Interval	What must be justified
Normal approximation	\widehat\theta\pm1.96\widehat{\mathrm{SE}}_*	A centered normal approximation and a consistent variance
Percentile	[q^*_{.025},q^*_{.975}]	An appropriate bootstrap approximation; this does not automatically correct bias
Basic	[2\widehat\theta-q^*_{.975},2\widehat\theta-q^*_{.025}]	Approximation of the centered error distribution
Studentized	Invert quantiles of (\widehat\theta^*-\widehat\theta)/\widehat{\mathrm{SE}}^*	A justified SE within each replicate and a valid studentized approximation

The basic interval follows by approximating the estimation-error distribution. The bootstrap error quantiles are q^*_{.025}-\widehat\theta and q^*_{.975}-\widehat\theta. If those bracket \widehat\theta-\theta_0, solving the two inequalities for \theta_0 reverses and reflects the endpoints:

\left[2\widehat\theta-q^*_{.975},\; 2\widehat\theta-q^*_{.025}\right].

For example, with observed estimate 10 and bootstrap estimate quantiles 8 and 13, the percentile interval is [8,13] and the basic interval is [7,12]. Their difference comes from the transformation, not from more data or a different replicate count. With a symmetric, centered normal approximation they can agree at first order; skewness and bias make the distinction more visible.

For studentization, calculate an SE in each resample and form T^*=(\widehat\theta^*-\widehat\theta)/\widehat{\operatorname{SE}}^*. If t^*_{.025} and t^*_{.975} are its quantiles, inversion gives

\left[\widehat\theta-t^*_{.975}\widehat{\operatorname{SE}},\; \widehat\theta-t^*_{.025}\widehat{\operatorname{SE}}\right].

The final SE here comes from the original study; each starred SE was needed to construct the resampled standardized errors. The procedure requires both the estimator and its variance calculation to behave appropriately under resampling.

The boot.ci documentation also describes BCa intervals. Their corrections and extra accuracy require additional conditions; choosing a sophisticated interval label does not remove those conditions.

q <- quantile(boot_mean, c(0.025, 0.975), names = FALSE)
rbind(normal = mean_hat + c(-1, 1)*qnorm(0.975)*sd(boot_mean),
      percentile = q, basic = 2*mean_hat-rev(q))

               [,1]     [,2]
normal     2.046090 2.710119
percentile 2.042315 2.715710
basic      2.040499 2.713894

In this particular simulated sample, all three intervals lie above the known mean 2. That is one random noncoverage event; it does not estimate a method’s coverage probability. Assess coverage across independent datasets under stated conditions.

These are pointwise intervals for one target. They do not automatically provide simultaneous or post-selection coverage for thousands of features. See Multiple testing when the tests are not independent.

5 Preserve the independent unit and the target weights

For independent people with repeated records, draw people with replacement and carry all of each sampled person’s records together. The resampled object should reproduce the sampling design. A bootstrap for a stratified design may need resampling within the original sampling strata; a complex survey or dependence across clusters needs its own scheme.

Resampling rows independently would break the pattern of shared person effects in the data. Holding people together preserves their joint records. This choice addresses dependence; the statistic recomputed on those records still determines whether the target weights people or records equally.

If the target is the mean of a person’s recorded-visit average, first form one summary per person. If it is a record-weighted ratio, retain each person’s total and record count and recompute the ratio after resampling. The following invented records demonstrate both, using the same resampled people.

set.seed(3102)
G <- 80; size <- sample(c(2, 5), G, replace = TRUE)
id <- rep(seq_len(G), size)
values <- rep(rnorm(G), size) + rnorm(sum(size), sd = 0.5)
total <- as.numeric(rowsum(values, id))
bootstrap <- replicate(1999, {
  selected <- sample.int(G, G, replace = TRUE)
  c(person = mean(total[selected]/size[selected]),
    record = sum(total[selected])/sum(size[selected]))
})
apply(bootstrap, 1, sd)

   person    record 
0.1259339 0.1324401 

Duplicating a person in a bootstrap sample is intentional. If rebuilding a raw-record analysis, give the duplicate cluster copies distinct bootstrap cluster IDs where needed. Do not replace person weighting with row weighting by accident. The Dependent observations and cluster-aware inference page derives the corresponding influence functions and few-cluster limitations.

6 Decide what must be refitted

A full bootstrap statistic repeats the estimator’s learned preprocessing, model fitting, selection, tuning, and uncertainty calculation according to its declared algorithm. Holding fitted parts fixed changes what is being approximated; it needs an argument. Repeating the full algorithm is also not sufficient for bootstrap consistency when the target or selection is nonregular.

Cross-fitting requires particular care. Copies of an original person’s outcomes can leak between training and evaluation if bootstrap duplicates are randomly sent to different folds. Specify a bootstrap weighting/identity and fold scheme whose theory applies to the estimator. Grouping duplicates can prevent that leakage, but does not itself prove conditional remainder control. There is no universal instruction to bootstrap every machine-learning pipeline in the same way.

6.1 Later extension: multiplier resampling (optional)

Return to this extension after learning influence functions. A multiplier approximation can avoid nuisance refitting when a valid asymptotically linear representation is already available. If the estimated centered contributions are \widehat D_i, independent multipliers \xi_i with mean zero and variance one, also independent of the data, generate

Z^*=\frac{1}{\sqrt n}\sum_i\xi_i\widehat D_i.

With a valid negligible estimator remainder, consistent contribution/variance estimation, and conditional CLT/non-dominance conditions, this approximates the first-order error law. For Gaussian multipliers and a scalar target it is conditionally normal with variance n^{-1}\sum_i\widehat D_i^2; that variance match alone does not prove the estimator’s expansion. With independent clusters, use cluster contributions and the cluster count. This construction approximates the first-order fluctuation; it does not measure persistent nuisance bias.

Conditional on the data, the \widehat D_i are fixed, and independent mean-zero multipliers make the covariance terms vanish. This explains the displayed variance. Z^* is on the scaled-error scale of \sqrt n(\widehat\theta-\theta_0); divide it by \sqrt n to approximate an unscaled estimation error. Inference still depends on having correctly derived the original estimator’s expansion.

7 A concrete failure: estimating a support endpoint

Let Y_i be i.i.d. uniform on [0,\theta] and estimate \theta by the sample maximum M_n. Here we study the positive gap \theta-M_n and its bootstrap counterpart M_n-M_n^*; both reverse the sign of the earlier estimate-minus-target convention. The scaling is n, rather than the usual square root. Define the scaled gap’s CDF as F_n(t)=\Pr\{n(\theta-M_n)/\theta\le t\}. For 0\le t\le n,

\begin{aligned} F_n(t)&=1-(1-t/n)^n\\ &\longrightarrow1-e^{-t}. \end{aligned}

To obtain this, note that M_n\le x requires all n independent observations to be at most x, giving probability (x/\theta)^n for 0\le x\le\theta. The scaled-error event instead asks whether M_n is at least \theta(1-t/n), so taking the complement gives the formula. Its limiting law is exponential, not the mean’s normal limit.

With a unique sample maximum, the ordinary bootstrap maximum equals M_n whenever that observed maximum is selected at least once:

\begin{aligned} \Pr_*(M_n^*=M_n)&=1-(1-1/n)^n\\ &\longrightarrow1-e^{-1}. \end{aligned}

Each draw misses that particular record with probability 1-1/n; all n draws miss it with probability (1-1/n)^n. The bootstrap cannot generate a larger value than the observed maximum. Thus it retains positive probability of exactly zero centered error even as the dataset grows.

The centered bootstrap error consequently retains an atom at zero, whereas the sampling limit is continuous. Increasing the number of bootstrap replicates estimates the wrong approximation more precisely. An ordinary percentile interval is particularly revealing: both endpoints are at most M_n<\theta almost surely, so it cannot cover the true endpoint. Specialized endpoint procedures require another argument; the failure does not say every resampling procedure fails for every endpoint problem.

set.seed(3103)
x <- runif(80); n_endpoint <- length(x); maximum <- max(x)
boot_max <- replicate(3999, max(x[sample.int(n_endpoint, n_endpoint, replace = TRUE)]))
endpoint_error <- n_endpoint*(maximum-boot_max)/maximum
plot(ecdf(endpoint_error), xlim = c(-0.2, 5), verticals = TRUE,
     do.points = FALSE, col = "#D55E00", lwd = 2, main = "",
     xlab = "Endpoint error on the n scale", ylab = "CDF")
t_grid <- seq(0, 5, length.out = 300)
lines(t_grid, 1-(1-t_grid/n_endpoint)^n_endpoint, col = "#2F80C8", lwd = 2)
legend("bottomright", c("Conditional bootstrap", "Exact sampling law"),
       col = c("#D55E00", "#2F80C8"), lty = 1, bty = "n")
c(empirical_atom = mean(boot_max == maximum),
  exact_atom = 1-(1-1/n_endpoint)^n_endpoint)

empirical_atom     exact_atom 
     0.6256564      0.6344319 

Figure 2: One uniform(0,1) teaching sample (n=80), with 3,999 ordinary bootstrap maxima. The conditional bootstrap CDF has a large jump at zero; the exact sampling CDF is continuous. This compares endpoint errors on their n scale, not sqrt(n).

Bickel and Freedman give foundational consistency results and counterexamples. Nonsmooth sample equations are not automatically failures: Quantiles and local expansions develops the regular positive-density case. Boundary behavior, unstable selection, weak identification, infinite moments, and the bootstrap’s own remainder conditions require examination.

8 Separate sampling uncertainty from resampling Monte Carlo error

The dataset size controls information about the population. The replicate count B controls simulation accuracy conditional on that dataset. At a fixed threshold, a bootstrap tail probability estimated from independent replicates has MCSE approximately \sqrt{p_*(1-p_*)/B}. An interval endpoint is a quantile; its Monte Carlo error also depends on local density and cannot be read directly from this probability-scale formula.

For example, if a conditional tail probability is 0.025, 1,999 resamples supply about 50 observations in that tail on average. Increasing B to 19,999 reduces probability-scale Monte Carlo error by about \sqrt{10} while leaving the original study unchanged. It does not divide the study’s SE by \sqrt{10}. For a cluster bootstrap, the population information is governed by the number and structure of independent clusters, rather than simply the row count.

Check endpoint stability with independent resampling seeds or a justified quantile-uncertainty calculation. Report the scheme, resampling unit, statistic, refitting/fold rules, replicate count, seed, interval construction, and failed or degenerate replicates. A failed replicate is an outcome to diagnose, not a row to silently discard.

9 Exercises

Try each checkpoint before opening its hint or solution.

Checkpoint · BOOT-01 — Derive what the bootstrap mean approximates

Task type: Derivation and interpretation.

For i.i.d. scalar observations with finite nonzero variance, derive the conditional mean and variance of the ordinary bootstrap sample mean. Explain the center and scaling of the bootstrap root, construct percentile and basic intervals from bootstrap quantiles, and distinguish n from the replicate count B. State what changes for a mean of person-level visit averages.

Ready to continue when: you can connect the empirical resampling law, the target error distribution, and the independent unit without confusing sampling error with Monte Carlo error.

Transfer: Would holding every fitted nuisance fixed approximate the same finite-sample algorithm as refitting it? When could a first-order multiplier approximation still be justified?

Hint

A bootstrap draw is uniform over the n observed values. Conditional on the dataset, the n draws are independent.

Solution outline — open after trying

The conditional mean is the observed mean and the variance is (n−1)s²/n². The bootstrap root sqrt(n)(mean*−mean) approximates sqrt(n)(mean−μ) under the mean CLT and bootstrap conditions. Percentile endpoints are the estimate quantiles; basic endpoints reverse and reflect them about the observed estimate. Increasing B reduces conditional simulation error, while increasing n adds population information. For equally person-weighted visit averages, resample whole people or their defined summaries and use the person count.

Checkpoint · BOOT-02 — Prove an ordinary bootstrap failure

Task type: Derivation and diagnosis.

For a continuous uniform sample with a unique maximum, derive the probability that an ordinary bootstrap sample contains that maximum. Compare its limiting atom with the sampling limit of n(θ−M_n)/θ. Explain why an ordinary percentile interval cannot cover θ, and whether increasing B or resampling rows instead of people could repair the respective validity problems.

Ready to continue when: you can diagnose a failure of the resampling approximation separately from limited replicate precision and an incorrect sampling unit.

Transfer: Why does a positive-density sample median require a different diagnosis from a support maximum?

Hint

Compute the probability of never drawing the maximum in n draws. Both percentile endpoints must be observed-sample-supported maxima.

Solution outline — open after trying

The probability is 1−(1−1/n)^n, tending to 1−exp(−1). The bootstrap centered error therefore has a nonvanishing atom at zero, while the sampling endpoint error converges to a continuous exponential law. Percentile endpoints are at most M_n, below θ almost surely, so coverage is zero. More replicates sharpen this invalid approximation; an independent-row bootstrap also cannot reproduce within-person dependence. Each problem requires a different justified procedure.

Additional practice (optional)

1. Use three independent bootstrap seeds with B=199 and B=1999. Compare interval endpoint stability without changing the dataset.
2. Derive the conditional variance of the Gaussian multiplier root and explain why matching it is not a proof of target consistency.
3. Write the full refitting and identity rules for bootstrapping one of your own estimators; identify which step still needs a bootstrap-validity argument.

10 Reading

• R boot documentation: statistic and resampling interfaces.
• R boot.ci documentation: interval constructions and their inputs.
• Bickel and Freedman (1981), Some asymptotic theory for the bootstrap.

11 Continue through the guide

Use Sensitivity analysis for assumptions and research decisions to investigate assumptions that resampling alone cannot assess. Use the group’s published manuscripts as case studies for identifying the target and uncertainty argument.