Survival Curves From Simulated Event Histories

From Simulated Histories to a Survival Curve

A survival curve is not a single fixed-time risk. It describes how the probability of remaining event-free changes over follow-up. In the parametric g-formula, the curve under a treatment regime is obtained by simulating longitudinal histories under that regime, recording each simulated event time, and averaging the resulting event-free indicators. The main practical issue is that the event process and the covariate process must be simulated in a coherent time order.

Let \(g\) denote a treatment regime and \(T^g\) the counterfactual event time under that regime. We use the convention that a person is event-free at time \(t\) when \(T^g>t\); an event occurring exactly at \(t\) therefore counts as having occurred by \(t\). The target population is the baseline population being averaged over, as in baseline averaging.

The simulation is from the fitted longitudinal data law with treatment set according to \(g\), and baseline covariates drawn from, or averaged over, the target population. Thus \(\widehat S_g(t)\) is a Monte Carlo approximation to the survival curve under that fitted law. It is not automatically a causal estimate: causal interpretation also depends on the identification assumptions and the adequacy of the specified models.

Ordering Events and Covariate Updates

For discrete-time data, specify what is known at the start of each interval and when the event indicator is assessed. One coherent convention is: at the start of interval \(k\), the person is event-free and has history \(H_{k-1}\); simulate the covariate value \(L_k\) for that interval; then simulate whether an event occurs in \((t_{k-1},t_k]\), conditional on the history and \(L_k\). If an event occurs, the person leaves the event-free process and no later covariates or event indicators need to be generated. This convention is appropriate when \(L_k\) is measured before the interval’s event assessment.

Other data structures may require another ordering. For example, if a covariate is recorded after an event assessment at a visit, it cannot predict an event that has already been assessed. The simulation order must match the measurement and risk-set definitions in the data. In continuous time, a time-varying covariate used to generate an event time must likewise be available immediately before the event time it helps predict.

The event of interest is distinct from loss to follow-up. When the target is survival from the event, administrative censoring in observed data should not be recoded as an event in simulated histories. The event-time model and the observation or censoring process have different meanings.

Why the Discrete-Time Calculation Works

Write \(D_k=1\) if the event occurs in interval \((t_{k-1},t_k]\), and \(D_k=0\) otherwise. For histories still event-free at \(t_{k-1}\), let \(p_k(H_{k-1},L_k)\) be the conditional probability of an event in interval \(k\), given the history available before the event assessment. The conditional probability of no event in that interval is \(1-p_k(H_{k-1},L_k)\).

To see the sequential calculation, condition first on the history at the start of interval \(K\), including the covariate value generated for that interval. Among histories that have remained event-free through \(t_{K-1}\), the probability of no event in interval \(K\) is \(1-p_K\). Multiply this factor by the probability of having been event-free through \(t_{K-1}\), then average over the histories that can arise under the regime. Repeating this conditioning for \(K-1,\ldots,1\) gives the probability of no event in every interval. In a simulation, generating the binary indicators performs these conditional steps directly; averaging \(I(D_1=\cdots=D_K=0)\) estimates their joint probability.

This result does not say that one can generally multiply interval survival probabilities calculated from marginal, population-level hazards. The histories of people remaining event-free change over time: people with higher event risk tend to experience events earlier and leave later risk sets. The conditional probabilities in the sequential calculation must retain the relevant history and risk-set conditioning.

Simulation Properties and a Marginal-Hazard Trap

The next result describes what the empirical curve guarantees when the histories are generated independently from the same fitted regime-specific law. It concerns Monte Carlo variation conditional on the fitted models; it does not account for uncertainty from estimating those models.

For the expectation, each summand \(I\{T^{g,(m)}>t\}\) has expectation \(P(T^g>t)=S_g(t)\). Linearity of expectation gives $$ E[\widehat S_g(t)] =\frac{1}{M}\sum_{m=1}^{M}E[I\{T^{g,(m)}>t\}] =\frac{1}{M}\sum_{m=1}^{M}S_g(t) =S_g(t). $$ For monotonicity, \(s\leq t\) implies \(\{T^{g,(m)}>t\}\subseteq\{T^{g,(m)}>s\}\) for every history \(m\). Each indicator at \(s\) is therefore at least its corresponding indicator at \(t\); summing and dividing by \(M\) preserves the inequality. \(\square\)

At a single time, the indicators are Bernoulli variables, so the conditional Monte Carlo variance is \(S_g(t)\{1-S_g(t)\}/M\leq 1/(4M)\). This quantifies simulation noise at that time, not model misspecification or sampling uncertainty in the original study. Using common simulated histories at all times also ensures the curve’s monotonicity by construction.

Checks, Comparisons, and Limitations

Several simple checks can expose implementation errors. A curve should start at one when everyone is event-free at baseline, remain between zero and one, and never increase. In discrete time, the value at \(t_k\) must exclude events in \((t_{k-1},t_k]\). Check that no post-event covariate update is used to generate an earlier event, and that a generated event is not accidentally treated as censoring. If the curve is unexpectedly jagged across independent simulation batches, increase \(M\) and compare the estimates; small \(M\) can create substantial Monte Carlo noise.

A fixed-follow-up risk reports only \(1-S_g(t)\) at one chosen time, whereas a survival curve displays that risk across times. An observed Kaplan–Meier curve summarizes survival in an observed group; it is not, by itself, a curve under a specified intervention regime. The simulated g-formula curve instead depends on the fitted covariate and outcome models, the treatment regime, and the target baseline distribution. It can be smooth or detailed in time, but simulation does not protect it from poorly specified models, inadequate support for the regime, or an incorrect event/covariate ordering.

A particularly tempting shortcut is to calculate a population-average hazard for each interval and multiply the corresponding interval survival probabilities. As the risk-group example shows, marginal risk-set composition changes over time, so that shortcut generally fails. Retain histories in the sequential simulation, or directly simulate event times and average event-free indicators.