Bootstrap Inference for G-Computation

Why Bootstrap the Entire g-Computation Procedure?

A g-computation estimate is not just a prediction from one fitted outcome model. It depends on the collection of conditional models used to generate counterfactual histories, the baseline population over which predictions are averaged, and the calculation of the final regime-specific contrast. Uncertainty intervals should reflect the sampling variability of that whole procedure. The nonparametric bootstrap does so by repeatedly resampling individuals and repeating every estimation step.

The resampling unit is the independent person, not an individual person-time record. If a person contributes baseline covariates and several follow-up records, all of that person’s records travel together into or out of a bootstrap sample. This preserves the within-person dependence present in the observed longitudinal data.

Definition: For a sample of $n$ independent individuals with complete observed records $O_1,\ldots,O_n$, a nonparametric bootstrap sample is a sample of size $n$ drawn with replacement from those $n$ records. A bootstrap g-computation replicate refits all conditional models using that sample, repeats g-computation for each regime, and recalculates the estimand or contrast.

For example, let $g$ denote a specified treatment regime and let $\psi(g)=E[Y^g]$ denote its mean counterfactual outcome at a fixed follow-up time. A contrast might be $\theta=\psi(g_1)-\psi(g_0)$. In each bootstrap replicate $b$, calculate both regime-specific estimates from the same resampled individuals and then form $\hat\theta^{*(b)}=\hat\psi^{*(b)}(g_1)-\hat\psi^{*(b)}(g_0)$. The paired difference is the bootstrap replicate of the contrast.

Implementation: Repeat Every Data-Dependent Step

A valid implementation follows the estimator as a complete pipeline. If the target population is the empirical baseline population, the bootstrap sample supplies that population as well as the data for fitting the models. If the analysis instead standardizes to a separately defined fixed population, that target must be handled consistently across replicates. For longitudinal g-computation, the conditional models include the models required to generate the covariate and outcome histories under intervention. Any model-selection or data-adaptive fitting step that is part of the reported procedure must also be repeated.

Resample individuals.
Draw $n$ person-level records with replacement. Retain every observed time point for each selected person; duplicate draws count as duplicate individuals.

Refit the conditional models.
Fit each model used by the g-formula in the bootstrap sample, using the same prespecified modeling procedure as in the original analysis.

Repeat g-computation.
For every regime, generate or integrate the counterfactual histories and average over the specified target population to obtain $\hat\psi^{*(b)}(g)$.

Form the contrast and store it.
Calculate $\hat\theta^{*(b)}=\hat\psi^{*(b)}(g_1)-\hat\psi^{*(b)}(g_0)$, or the corresponding ratio or other prespecified contrast.

Summarize the replicates.
Use the empirical distribution of the $B$ bootstrap contrasts to calculate a standard error or confidence interval, and inspect failed fits and extreme values.

Use enough replicates for the interval endpoints to be stable. Increasing $B$ reduces numerical error in approximating the bootstrap distribution; it does not increase the original sample size or repair a poorly supported model. When g-computation itself uses simulated trajectories, simulation error can add noise to each replicate. Analytic integration, sufficiently large simulation runs, or carefully coordinated random numbers can reduce this noise. The simulation settings should be checked so that the resulting interval reflects sampling variability rather than avoidable simulation fluctuation.

Worked Example: One Bootstrap Replicate for a Risk Difference

Consider a one-time treatment and binary outcome, with no covariates. For illustration, use a saturated outcome model: within each treatment level, the fitted risk is the observed proportion of outcomes equal to 1. The g-computation estimates under treatment and control are those two fitted risks. The six original records are:

IDs 1, 2, and 3 have treatment $A=1$ and outcomes $Y=1,0,1$. IDs 4, 5, and 6 have $A=0$ and outcomes $Y=0,0,1$. Thus the original estimated risks are $2/3$ under treatment and $1/3$ under control, and the estimated risk difference is $2/3-1/3=1/3$.

Suppose one bootstrap draw selects IDs $1,1,2,4,5,6$. This is a sample of six individuals drawn with replacement; ID 1 appears twice, so its record contributes twice. The treated outcomes in this sample are $1,1,0$. The refitted treated risk is $2/3$. The control outcomes are $0,0,1$, so the refitted control risk is $1/3$. The bootstrap risk difference is therefore $2/3-1/3=1/3$.

This small example makes the resampling unit and refitting explicit. In a longitudinal analysis, each selected ID would bring along the entire observed history, and every conditional model would be refitted before computing the regime risks. A bootstrap sample can occasionally contain no observations in a treatment-history group needed by a model; that is a fitting or support failure to record and diagnose, not a reason to silently substitute the original fit.

Confidence Intervals From Bootstrap Replicates

Let $\hat\theta$ be the original-sample estimate and let $\hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}$ be the bootstrap contrasts. The bootstrap standard error is the sample standard deviation of these replicates:

\widehat{\operatorname{SE}}_{\mathrm{boot}}(\hat\theta) = \left\{\frac{1}{B-1}\sum_{b=1}^B \left(\hat\theta^{*(b)}-\bar{\theta}^{*}\right)^2\right\}^{1/2}, \qquad \bar{\theta}^{*}=\frac{1}{B}\sum_{b=1}^B\hat\theta^{*(b)}.

A percentile interval uses the empirical quantiles of the bootstrap estimates. A basic interval reflects the bootstrap deviations around $\hat\theta$ back across the original estimate. These are distinct constructions; neither should be treated as automatically reliable when the estimator is irregular or the fitted models are unstable.

\text{Percentile interval: } \left[ \hat Q^{*}_{\alpha/2},\, \hat Q^{*}_{1-\alpha/2} \right], \qquad \text{Basic interval: } \left[ 2\hat\theta-\hat Q^{*}_{1-\alpha/2},\, 2\hat\theta-\hat Q^{*}_{\alpha/2} \right].

Worked Example: A Percentile Interval and Its Quantiles

Suppose, only to illustrate the arithmetic, that $B=10$ bootstrap risk differences sorted from smallest to largest are $-0.10,-0.04,0.00,0.02,0.05,0.07,0.09,0.12,0.16,0.20$. Use the common linear-interpolation quantile convention, in which the location for probability $p$ is $1+(B-1)p$.

For a nominal 95% percentile interval, the lower location is $1+9(0.025)=1.225$. Interpolating 0.225 of the way from the first to the second value gives $-0.10+0.225[-0.04-(-0.10)]=-0.10+0.0135=-0.0865$.

The upper location is $1+9(0.975)=9.775$. Interpolating 0.775 of the way from the ninth to the tenth value gives $0.16+0.775(0.20-0.16)=0.16+0.031=0.191$. The resulting illustrative interval is $[-0.0865,0.191]$. Ten replicates are far too few for a practical 95% interval; the calculation demonstrates the quantile rule, not a recommended number of replicates.

When Does Bootstrap Inference Work?

The usual justification requires independent, identically distributed individuals and a sufficiently regular estimator. For a parametric g-computation estimator, useful sufficient conditions include a fixed follow-up horizon; correctly specified, identifiable conditional models; parameters in the interior of their parameter spaces; nonsingular information matrices; and a smooth g-computation functional. The relevant model fits must have finite, well-behaved asymptotic variances. In particular, severe nonpositivity, separation, boundary estimates, or non-identifiability can violate the regularity needed for ordinary bootstrap approximations.

Theorem: Suppose the estimator of the regime contrast is asymptotically linear, so that for an influence function $\phi$ with $E[\phi(O)]=0$ and $E[\phi(O)^2]=V$, where $0<V<\infty$, $$ \sqrt n(\hat\theta-\theta)=\frac{1}{\sqrt n}\sum_{i=1}^n\phi(O_i)+o_p(1). $$ Suppose also that the same expansion holds conditionally for the ordinary individual-level bootstrap, with the empirical-centered influence values $\phi(O_i^*)-\bar\phi_n$, where $\bar\phi_n=n^{-1}\sum_i\phi(O_i)$, and with a conditional remainder that is $o_p(1)$ in probability. Then the conditional distribution of $\sqrt n(\hat\theta^*-\hat\theta)$ converges in probability to $N(0,V)$. If the limiting distribution is continuous, the percentile interval has asymptotic coverage $1-\alpha$.

The expansion condition is important: the result does not assert that every g-computation estimator has a valid ordinary bootstrap. It is a regularity result for estimators with the stated first-order behavior. A smooth finite-dimensional fitting and standardization procedure often has this form when the model and support conditions hold.

Proof. By the assumed conditional expansion, $\sqrt n(\hat\theta^*-\hat\theta)=n^{-1/2}\sum_{i=1}^n\{\phi(O_i^*)-\bar\phi_n\}+o_p^*(1)$, where the remainder vanishes in conditional probability, in probability over the original samples. Given the observed sample, the $O_i^*$ are independent draws from its empirical distribution. The conditional variance of one centered influence value is $n^{-1}\sum_i\{\phi(O_i)-\bar\phi_n\}^2$, which converges in probability to $V$ by the law of large numbers. Finite second moment also gives the conditional Lindeberg condition in probability: for each $\epsilon>0$, the empirical average of squared centered influence values exceeding $\epsilon\sqrt n$ tends to zero in probability. The conditional central limit theorem therefore gives convergence of the leading sum to $N(0,V)$; the conditional remainder does not change that limit. The ordinary asymptotic expansion and the central limit theorem give $\sqrt n(\hat\theta-\theta)\xrightarrow{d}N(0,V)$. Consequently the bootstrap quantiles of the centered, scaled estimate converge to the corresponding normal quantiles, yielding asymptotic coverage $1-\alpha$ for the percentile interval. $\square$

Worked Example: Why the Two Regimes Must Stay Paired

Suppose four bootstrap replicates give regime-specific risks $\hat\psi^{*(b)}(g_1)=(0.40,0.50,0.60,0.50)$ and $\hat\psi^{*(b)}(g_0)=(0.20,0.30,0.40,0.30)$. The replicate contrasts, calculated within each resample, are $0.20,0.20,0.20,0.20$; their sample standard deviation is zero.

Each set of regime-specific risks separately has sample variance $0.00667$: its mean is 0.50 for $g_1$, and the squared deviations are $0.01,0,0.01,0$, whose sum divided by $4-1$ is $0.00667$. The control risks have the same deviations and variance. Their sample covariance, using paired replicates, is also $0.00667$, because the two centered sequences are identical. Thus the variance of the paired difference is $0.00667+0.00667-2(0.00667)=0$, consistent with the four identical contrasts.

This deliberately simple set of replicates shows why estimating uncertainty for each regime separately and then combining those uncertainties can be misleading: the covariance matters. In practice, compute each contrast within each bootstrap sample and use the distribution of those contrasts directly.

Diagnostics, Alternatives, and Limitations

Examine the full set of bootstrap replicates rather than relying only on a standard error or interval. Check how many fits failed, whether estimates are extreme or highly skewed, whether particular treatment histories repeatedly produce sparse data, and whether the interval changes materially with a larger $B$. Inspect the fitted-model diagnostics already used in the original analysis. A bootstrap cannot make an unsupported extrapolation credible or correct misspecified conditional models: it measures sampling variability of the procedure being repeated.

Ordinary individual-level resampling is appropriate when individuals are independent sampling units. If sampling instead occurs in clusters, the resampling unit must reflect that design. A parametric bootstrap simulates new data from fitted parametric distributions rather than resampling observed records; it can be useful when the parametric sampling model is the intended uncertainty model, but depends on that model. Analytic variance methods can be more computationally efficient when their derivations and regularity conditions apply. Bootstrap inference is flexible, but may be computationally expensive and can be unreliable for irregular estimators, very small samples, or near-violations of positivity.

The resampling scheme also explains why whole-person records must stay intact. Each bootstrap draw is from the empirical distribution of complete records $O_i$, not from separate empirical distributions for each time point. Drawing records as units retains the observed dependence across times; independently resampling person-time rows would create histories that were never observed and would not reproduce the sampling structure of the data.

Common mistake: Resampling each person-time row independently, then treating the resulting rows as a bootstrap sample of people. This breaks within-person dependence and can create artificial histories. Instead, draw $n$ individuals with replacement and carry every selected individual’s full record into the replicate. Also refit all conditional models; holding the original fits fixed omits a major source of sampling uncertainty.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In a study with five records per person, what is the appropriate ordinary bootstrap resampling unit, and what happens to the five records when a person is selected twice?
Why should the two regime-specific estimates be calculated in the same bootstrap replicate before their difference is formed?
Give two reasons why many g-computation bootstrap replicates might fail or be unstable, and explain why increasing $B$ does not fix them.
Under what key condition does the theorem above justify ordinary bootstrap intervals for a regime contrast?

Solutions. (1) Resample individuals. If a person is selected twice, that person’s complete five-record history appears twice and contributes twice to fitting and standardization. (2) The two estimates from a shared sample are dependent; forming their difference within each replicate preserves their covariance and directly approximates the sampling distribution of the contrast. (3) Sparse treatment histories can cause failed or unstable fits, and separation or near-positivity violations can produce extreme estimates. More replicates reduce Monte Carlo error in the bootstrap distribution but do not add information to the original sample or repair these problems. (4) The estimator must have the stated regular, asymptotically linear expansion, and the bootstrap must reproduce that expansion conditionally; under the theorem’s finite-variance and continuity conditions, the centered bootstrap distribution then estimates the sampling distribution.