Why Simulate Longitudinal Histories?
The longitudinal g-formula expresses a regime-specific mean by averaging outcomes over the covariate histories that would arise under that regime. When histories have many components, or include continuous covariates, computing this average analytically may be difficult. Monte Carlo g-computation approximates it by repeatedly drawing histories from fitted conditional distributions and averaging the resulting outcome predictions.
The simulation must follow the time order of the study. At each visit, assign the treatment specified by the regime, then draw the next covariate value conditional on the history simulated so far. After the final treatment and covariate history, predict the outcome. This is not a simulation of observed treatment choices: treatments are set by the regime being evaluated.
The Simulation Target and Its Assumptions
Write \(\bar L_t=(L_0,\ldots,L_t)\) and \(\bar A_t=(A_0,\ldots,A_t)\). For a given regime, the covariate history is generated from the conditional distributions of covariates given the past and the treatments assigned so far. For example, the distribution of \(L_{t+1}\) is conditional on \(\bar L_t\) and \(\bar A_t\), with each \(A_j\) set by \(d_j\). The outcome model gives the conditional mean of \(Y\) given the full history.
A simulation-based estimate relies on the same causal identification conditions used for the g-formula: consistency under the regime, sequential exchangeability given the measured history, and positivity for the treatment choices required by the regime. It also relies on adequate fitted models for the covariate transition distributions and the conditional outcome mean. These conditions have distinct roles: the causal conditions link observed data to regime-specific counterfactual quantities, while the models approximate the distributions and means needed to calculate those quantities.
Let \(\widehat q_t\) denote a fitted conditional distribution for \(L_{t+1}\), and let \(\widehat m\) denote a fitted conditional mean for \(Y\). With an empirical target population of \(n\) baseline records, the fitted g-formula mean for regime \(d\) averages \(\widehat m\) over each record’s regime-generated covariate histories. Monte Carlo replaces that history average with a finite number of simulated histories.
Here \(M\) is the number of simulated histories per baseline record, and \(\bar A_{K-1,ir}\) is assigned by regime \(d\) along simulated history \(r\) for record \(i\). If the target is a population rather than a fixed data set, baseline covariates can instead be sampled from an appropriate estimate of that population’s distribution. The average is then over both sampled baselines and simulated histories.
A Practical Simulation Procedure
For each target baseline record, start with its observed \(L_0\). In each replicate, apply the regime’s treatment rule, draw the next covariate from its fitted conditional distribution, and continue in chronological order. Once the full history is available, evaluate the fitted outcome mean. The final estimate averages these predictions over all baseline records and replicates.
Specify the population of baseline covariate records over which the regime-specific mean is wanted.
At each time, evaluate the static assignment or dynamic rule using only the history available at that time.
Use the fitted conditional distribution given the simulated past and the treatment assigned by the regime.
Evaluate the fitted conditional outcome mean for each completed history, then average across records and replicates.
When \(Y\) is binary, \(\widehat m\) is usually a predicted risk, not a simulated zero-or-one outcome. Averaging predicted risks integrates over outcome randomness analytically while simulating the covariate histories. One may also draw outcomes, but that adds Monte Carlo variation without changing the target mean.
Worked Example: A Sustained Two-Visit Regime
Suppose hypothetical fitted models describe two binary intermediate covariates, \(L_1\) and \(L_2\), and a binary final outcome. Consider the sustained regime \(A_0=0,A_1=1\), with one fixed baseline stratum. The fitted transition model gives \(P(L_1=1\mid A_0=0)=0.40\). For the second covariate, \(P(L_2=1\mid L_1=0,A_0=0,A_1=1)=0.20\) and \(P(L_2=1\mid L_1=1,A_0=0,A_1=1)=0.70\). The fitted outcome risks are:
| \(L_1\) | \(L_2\) | Predicted risk of \(Y=1\) |
|---|---|---|
| 0 | 0 | 0.10 |
| 0 | 1 | 0.25 |
| 1 | 0 | 0.30 |
| 1 | 1 | 0.55 |
For illustration, a Monte Carlo run draws six histories in time order. Its simulated \(L_1\) values are \(0,1,0,0,1,0\). Conditional on those values, the simulated \(L_2\) values are \(0,1,1,0,1,0\). These draws follow the specified transition probabilities; a finite run need not reproduce them exactly.
The resulting histories have predicted outcome risks \(0.10,0.55,0.25,0.10,0.55,0.10\). Thus the Monte Carlo estimate for this run is
The estimate is a predicted risk of 0.275, or 27.5%, in this hypothetical target stratum under the specified regime. It is a simulation approximation, not an observed proportion of outcomes in six people.
For comparison, the exact mean under these fitted models can be calculated by enumerating all histories. For \(L_1=0\), the conditional mean risk is \(0.80(0.10)+0.20(0.25)=0.13\). For \(L_1=1\), it is \(0.30(0.30)+0.70(0.55)=0.475\). Weighting by the fitted distribution of \(L_1\) gives \(0.60(0.13)+0.40(0.475)=0.268\). The simulation result \(0.275\) differs because it used only six draws; increasing \(M\) would generally reduce this finite-simulation discrepancy.
Worked Example: Applying a Dynamic Rule
Consider a hypothetical target population with four baseline records, two with \(L_0=0\) and two with \(L_0=1\). At the first treatment time the regime sets \(A_0=0\). It then sets \(A_1=1\) if \(L_1=1\), and \(A_1=0\) if \(L_1=0\). The fitted model gives \(P(L_1=1\mid L_0=0,A_0=0)=0.25\) and \(P(L_1=1\mid L_0=1,A_0=0)=0.65\). Suppose the fitted outcome risks for histories allowed by this rule are 0.10 for \((L_0,L_1)=(0,0)\), 0.20 for \((0,1)\), 0.25 for \((1,0)\), and 0.35 for \((1,1)\).
Simulate two histories for each baseline record. For the two \(L_0=0\) records, suppose the draws are \(L_1=0,1\) for each record. For the two \(L_0=1\) records, suppose the draws are \(L_1=1,0\) for each. Treatment \(A_1\) is assigned after each \(L_1\) draw, using the rule, so the corresponding risks are \(0.10,0.20,0.10,0.20,0.35,0.25,0.35,0.25\).
Averaging over the eight simulated histories gives
The estimated risk is 22.5% for this hypothetical empirical target population. To check the direction and scale of the result, compute the exact fitted mean: for baseline value 0 it is \(0.75(0.10)+0.25(0.20)=0.125\); for baseline value 1 it is \(0.35(0.25)+0.65(0.35)=0.315\). The four baseline records give equal weight to each baseline value, so the exact mean is \((0.125+0.315)/2=0.220\). The Monte Carlo estimate is close but not identical. Crucially, each treatment assignment followed the dynamic rule after the relevant covariate was generated.
Worked Example: Continuous Covariate Histories
Suppose a hypothetical analysis targets two baseline records with \(L_0=0\) and \(L_0=2\), equally weighted. It evaluates the fixed regime \(A_0=1,A_1=0\). The fitted models are \(L_1\mid L_0,A_0=1\sim N(0.5+0.4L_0,0.5^2)\) and \(L_2\mid L_1,A_1=0\sim N(0.2+0.5L_1,0.4^2)\). The fitted conditional outcome mean, in outcome points, is \(\widehat m=10+2L_0-L_1+3L_2-A_0+0.5A_1\).
Draw two histories for each baseline record. Suppose the simulated pairs \((L_1,L_2)\) are \((0.2,0.4)\) and \((0.8,0.5)\) for \(L_0=0\), and \((1.0,0.7)\) and \((1.8,1.0)\) for \(L_0=2\). For this regime, the outcome prediction simplifies to \(9+2L_0-L_1+3L_2\). The four predictions are therefore \(9-0.2+1.2=10.0\), \(9-0.8+1.5=9.7\), \(9+4-1.0+2.1=14.1\), and \(9+4-1.8+3.0=14.2\) points.
The Monte Carlo mean is \((10.0+9.7+14.1+14.2)/4=12.0\) outcome points. It averages the conditional outcome predictions over the simulated continuous histories; it does not require a discrete table of every possible history. The particular four-draw estimate is not the exact fitted g-formula mean. For example, the model implies \(E(L_1)=0.9\) across the two equally weighted baseline records and \(E(L_2)=0.2+0.5(0.9)=0.65\). Linearity of the stated outcome model gives an exact mean of \(10+2(1)-0.9+3(0.65)-1=12.05\) points. The small difference reflects finite simulation.
Why the Average Works
Conditional on the fitted models and a fixed baseline record, each simulated history is drawn from the fitted sequence of covariate distributions under the regime. The conditional expectation of its predicted outcome is therefore the fitted g-formula average for that record. Averaging independent replicates estimates that expectation.
More formally, let \(Z_{ir}=\widehat m(\bar L_{K,ir},\bar A_{K-1,ir})\) be the predicted outcome for replicate \(r\) and baseline record \(i\). Given the fitted models and baseline data, the replicates for record \(i\) have the same conditional mean, say \(\mu_i\). Thus
If the replicates are independent and have finite variance \(\sigma_i^2\), the conditional variance of their average is \(\sigma_i^2/M\). Increasing \(M\) reduces Monte Carlo error for the fitted-model target. It does not, by itself, correct misspecified models, remove bias from unmeasured confounding, or quantify uncertainty from estimating the models and target population.
Using conditional outcome predictions rather than simulated outcomes also has a variance advantage. If one additionally draws \(Y\) given each history, the law of total variance gives
Thus averaging conditional means avoids the extra variation from drawing outcomes while preserving the same mean target.
Common Mistakes and Practical Checks
- Drawing treatment from its observed-data model. That simulates naturally occurring treatment, not outcomes under the regime of interest. Set each treatment using the specified static assignment or dynamic rule.
- Drawing covariates from the wrong conditional distribution. A later covariate must be drawn conditional on the simulated history and the regime-assigned treatments so far, not on the original observed treatment history.
- Applying a dynamic rule too early. A rule that depends on \(L_t\) cannot assign \(A_t\) before that covariate is available. Simulate and assign in chronological order.
- Treating a larger \(M\) as a cure for model problems. More draws reduce simulation error around the fitted-model estimate. They do not repair incorrect transition or outcome models, failed causal assumptions, or lack of support for the regime.
- Reporting a simulated risk as an observed event proportion. When the final values are fitted conditional means, their average estimates a model-based regime-specific risk; it is not the fraction of simulated people who experienced the outcome.
A useful implementation check is to inspect simulated histories at every time point: verify that treatment follows the regime, that each covariate draw uses the correct preceding history, and that predictions are averaged over the intended baseline population. The number of replicates should be large enough that repeated runs give a stable Monte Carlo estimate, but stability alone is not evidence that the fitted causal models are adequate.
Check Your Understanding
Answer these questions using the definitions and calculations above.
- In a dynamic regime that assigns treatment from the current covariate, when should the treatment rule be evaluated relative to drawing that covariate?
- Why does a Monte Carlo estimate from six simulated histories generally differ from the exact fitted g-formula mean?
- For a binary outcome, why might averaging predicted risks be preferable to simulating a binary outcome for every history?
- What does increasing the number of simulated histories reduce, and what important sources of error does it not remove?
- Explain why drawing treatments from an observed-data treatment model does not implement a specified treatment regime.