Worked Example With Continuous Confounding

Why Continuous Confounding Changes the Calculation

When a time-varying confounder is continuous, it is not enough to divide the data into a few confounder categories and compare outcomes within those groups. The parametric g-formula instead models how the continuous confounder evolves, then averages predicted outcomes over the trajectories that would occur under each treatment regime. This example follows two treatment decisions and a binary outcome.

Let $L_0$ be a continuous baseline covariate, $A_0$ the first treatment, $L_1$ a continuous confounder measured after $A_0$, $A_1$ the second treatment, and $Y$ a binary outcome measured after $A_1$. A fixed regime $g=(a_0,a_1)$ sets treatment to $a_0$ at the first decision and $a_1$ at the second, regardless of the observed treatment choices. Write $Y^g$ for the outcome that would be observed under this regime.

Definition: The regime-specific risk is $R(g)=P(Y^g=1)$. In this two-decision setting, the g-formula expresses it by averaging the outcome probability over the baseline distribution and the distribution of $L_1$ that would arise under $A_0=a_0$: $$ R(a_0,a_1)=\int E[Y\mid L_0=l_0,A_0=a_0,L_1=l_1,A_1=a_1]\, f(l_1\mid l_0,a_0)\,f(l_0)\,dl_1\,dl_0. $$

The expression requires that the observed-data conditional distributions represent the corresponding counterfactual distributions. Sufficient standard conditions are consistency; sequential exchangeability at each treatment decision, conditional on the measured history available then; and sequential positivity for the regimes and histories being considered. In this setup, exchangeability requires that, given $L_0$, $A_0$ is independent of the relevant future counterfactuals, and that, given $(L_0,A_0,L_1)$, $A_1$ is independent of the relevant future counterfactual outcome. Positivity requires positive probability of the regime’s treatment choices at histories occurring under that regime. Identification also presumes that the variables needed to control confounding have been measured and correctly ordered in time.

Models and the Calculation

For a parametric analysis, specify models for the baseline distribution, the transition distribution of $L_1$, and the conditional outcome probability. For example, one may use a normal linear model for $L_1$ and a logistic regression for $Y$. The outcome model must include the treatment and confounder terms needed to describe their conditional relationship; interactions or nonlinear terms may be necessary. Fit these models using observed data, then set treatment to the regime values when predicting.

Treatment models are not needed to generate trajectories under a fixed regime: the intervention supplies $A_0$ and $A_1$. They are still relevant to the identification assumptions and to diagnosing whether the observed data support the regimes. The baseline distribution can be represented by the empirical distribution of observed $L_0$, or by a suitable parametric model. For each baseline value, generate $L_1$ from its fitted conditional distribution under the chosen $A_0$, then predict the outcome probability with $A_0,A_1$ fixed at their regime values. Averaging these probabilities estimates risk; simulating a binary outcome for every trajectory is unnecessary if the target is its mean.

Fit the conditional models.
Estimate $f(L_0)$, $f(L_1\mid L_0,A_0)$, and $P(Y=1\mid L_0,A_0,L_1,A_1)$ from the observed data.

Set the intervention.
For each regime, replace the observed treatment values by its specified $a_0,a_1$ when generating and predicting.

Generate and average trajectories.
Draw baseline values and then $L_1$ values from the fitted transition model; average the predicted outcome probabilities to estimate $R(g)$.

Contrast the risks.
For two regimes $g_1,g_0$, calculate $R(g_1)-R(g_0)$, and, if useful, $R(g_1)/R(g_0)$.

Worked Example: Simulating Risks Under Two Regimes

Suppose fitted models from a hypothetical study are $$ L_1\mid L_0,A_0\sim N(0.5+0.6L_0+0.4A_0,\;0.25), \qquad \operatorname{logit}P(Y=1\mid L_0,A_0,L_1,A_1) =-2+0.5A_0+0.8A_1+0.4L_0+0.7L_1. $$ The second parameter of the normal distribution is its variance. To show the arithmetic transparently, use four baseline values and four standard-normal draws as a small illustrative simulation. These are not intended to provide precise inferential estimates.

Take $L_0=(-1,0,1,2)$ and transition errors $\epsilon=(-1,1,-1,1)$. Under $A_0=0$, the simulated confounder values are $L_1=0.5+0.6L_0+0.5\epsilon$, giving $(-0.6,1,0.6,2.2)$. Under $A_0=1$, add $0.4$ to each value, giving $(-0.2,1.4,1,2.6)$. The same baseline values and errors are used for both regimes so that the comparison is not obscured by different random draws.

For regime $(0,0)$, the four outcome linear predictors are $-2.82,-1.30,-1.18,0.34$. Applying $\operatorname{expit}(x)=1/(1+\exp(-x))$ gives probabilities approximately $0.0563,0.2142,0.2350,0.5842$. Their mean is $$ \widehat R(0,0)=\frac{0.0563+0.2142+0.2350+0.5842}{4} =\frac{1.0897}{4}\approx0.2724. $$

For regime $(1,1)$, the linear predictors are $-1.24,0.28,0.40,1.92$, giving probabilities approximately $0.2244,0.5695,0.5987,0.8721$. Thus $$ \widehat R(1,1)=\frac{0.2244+0.5695+0.5987+0.8721}{4} =\frac{2.2647}{4}\approx0.5662. $$

The estimated risk difference is $0.5662-0.2724=0.2938$, or about 29.4 percentage points. The risk ratio from the displayed rounded risks is $0.5662/0.2724\approx2.08$. These contrasts summarize the two regimes in the target population represented by the baseline distribution; they are not conditional effects for one particular value of $L_0$.

Why Averaging the Trajectories Matters

A nonlinear outcome model makes the order of prediction and averaging important. In general, the outcome probability evaluated at the mean confounder is not the average outcome probability across the confounder distribution. The g-formula averages predictions over the distribution of $L_1$, rather than simply substituting its conditional mean into the outcome model.

Worked Example: The Mean Confounder Is Not a Substitute for Integration

Suppose under a regime $L_1\sim N(0,1)$, and the fitted outcome probability is $\operatorname{expit}(-1+L_1)$. A simple two-point integration approximation uses $L_1=-1$ and $L_1=1$, each with weight one-half. The corresponding risks are $\operatorname{expit}(-2)=0.1192$ and $\operatorname{expit}(0)=0.5000$, so their weighted average is $(0.1192+0.5000)/2=0.3096$. In contrast, predicting at the mean $L_1=0$ gives $\operatorname{expit}(-1)=0.2689$. The two-point calculation is an approximation to the integral, not its exact value, but it clearly shows why substituting the mean can change the answer.

For continuous $L_1$, integration can be done by Monte Carlo simulation, numerical quadrature, or an analytic calculation when the models permit it. A Monte Carlo average approximates the integral; its simulation error can be reduced by increasing the number of generated trajectories. This does not correct misspecified models, unmeasured confounding, or unsupported extrapolation.

Result and justification: Under consistency, sequential exchangeability, and positivity, the displayed g-formula identifies $R(a_0,a_1)$. To see why, condition first on $L_0$ and then on $L_1$ under the intervention. By consistency, among individuals receiving the specified treatment at each decision, their observed outcome equals the corresponding counterfactual outcome. Sequential exchangeability allows those observed conditional outcome and confounder distributions to stand in for the counterfactual distributions at the same histories. Positivity ensures these conditional quantities are defined for histories generated under the regime. Iterated expectation therefore gives $$ E[Y^{a_0,a_1}] =\int E[Y^{a_0,a_1}\mid l_0,l_1]\, f(l_1^{a_0}\mid l_0)f(l_0)\,dl_1\,dl_0 =\int E[Y\mid l_0,a_0,l_1,a_1]\, f(l_1\mid l_0,a_0)f(l_0)\,dl_1\,dl_0. $$ Since $Y$ is binary, $E[Y^{a_0,a_1}]=P(Y^{a_0,a_1}=1)$, yielding the risk formula. $\square$

Diagnostics, Alternatives, and Limits

Continuous covariates make model checking especially important: a model can appear adequate on average while missing curvature, tail behavior, or treatment-specific relationships. Examine residuals from the continuous transition model for departures from its assumed conditional distribution, and compare simulated $L_1$ distributions with observed distributions in relevant treatment and baseline-history groups. For the outcome model, assess calibration and whether the fitted probabilities capture nonlinear associations and plausible interactions. Also compare regime-generated covariate values with the range of values observed at corresponding histories. Large regions outside that range are warnings that the result depends on extrapolation.

A practical positivity warning can be seen in the treatment probability. If $P(A_0=1\mid L_0=l)=\operatorname{expit}(-2+1.5l)$, then at $l=-2$ this probability is $\operatorname{expit}(-5)\approx0.0067$. There is very little observed treatment experience at that history. The g-formula may still produce a numerical prediction for a regime setting $A_0=1$, but the data provide little support for that prediction there. A precise-looking estimate does not remove this limitation.

Worked Example: Recognizing Sparse Treatment Support

With the treatment model above, at $L_0=0$, the estimated probability of $A_0=1$ is $\operatorname{expit}(-2)=1/(1+\exp(2))\approx0.1192$. At $L_0=-2$, it is about $0.0067$. The latter corresponds to roughly 7 treated individuals per 1,000 with that covariate value if the model is calibrated and the sample is large. A regime setting $A_0=1$ for people in this region consequently relies heavily on the fitted outcome and transition models rather than on close treated analogues. The appropriate response is to report the support limitation and assess sensitivity to plausible model choices, not to interpret the model’s extrapolated risk as directly observed evidence.

Compared with inverse-probability weighting, g-computation models the outcome and covariate transitions and then standardizes predictions, rather than constructing the regime contrast primarily through weights. It can be more stable when treatment probabilities are extreme, but it is vulnerable to misspecification of the models used for standardization. Analytic integration can avoid simulation error when available; Monte Carlo integration is flexible for complicated trajectories but adds a numerical error component. None of these computational choices repairs failures of identification assumptions or lack of empirical support.

Monte Carlo result: If the fitted models and baseline values are held fixed, and generated trajectories are independent draws from the specified transition distribution, the average predicted risk is unbiased for the corresponding finite-sample model-based integral. Indeed, if $Q_i$ is the predicted probability for generated trajectory $i$, then $E[Q_i]$ equals that integral by the definition of expectation, so $E[n^{-1}\sum_{i=1}^n Q_i]=n^{-1}\sum_{i=1}^n E[Q_i]$. This statement concerns simulation error conditional on the fitted models; it does not establish that those models or the causal assumptions are correct. $\square$

The most frequent implementation error is to predict the outcome using the observed $A_0,A_1$ values after generating trajectories, rather than setting both treatments to the intervention values. Another is to predict at a single mean confounder value and treat that as a population risk. Both depart from the target regime calculation. Check that every simulated trajectory uses the prescribed treatments and that the final risk averages over the baseline and time-varying confounder distributions.

Check Your Understanding

Answer the questions, then compare with the solutions below.

What distribution for $L_1$ is needed to calculate the risk under regime $(1,0)$?
Why does a fitted treatment model not itself generate treatment values in a fixed-regime g-computation?
For two regimes with estimated risks $0.40$ and $0.25$, calculate the risk difference and risk ratio.
Give one reason a very large number of simulated trajectories may still yield a misleading risk estimate.

Solutions. (1) Generate $L_1$ from $f(L_1\mid L_0,A_0=1)$, then predict the outcome with $A_0=1,A_1=0$. (2) The intervention fixes treatment to its specified values, so the treatment model is not used to draw treatment; it is relevant to identifying assumptions and support diagnostics. (3) The risk difference is $0.40-0.25=0.15$, and the risk ratio is $0.40/0.25=1.60$. (4) A large simulation only reduces Monte Carlo error. Misspecified transition or outcome models, unmeasured confounding, or positivity violations can still make the estimate misleading.

Worked Example With Continuous Confounding

What You'll Learn