Discrete-Event Simulation Using RWE Inputs

In plain language

Discrete-event simulation (DES) is a way to model a patient's medical journey by tracking what happens to each individual person one event at a time — things like a cancer progressing, a treatment being switched, or a patient dying — and measuring the cost and health impact of each step. Instead of pushing everyone through the same calendar-based clock, the model fast-forwards each person to their next meaningful event and asks: what did that cost, and how many healthy days did it buy? This lets researchers compare two treatments over a lifetime, accounting for the fact that a patient who had a bad side effect on line 1 therapy tends to do worse on line 2 — something simpler models cannot easily represent.

Discrete-event simulation (DES)

is an individual-level health-economic modeling method in which each virtual patient is tracked as a distinct entity carrying its own attributes (age, line of therapy, prior adverse-event history, accumulated cost and time-in-state). The patient's future is governed by a future-event list: the model samples a time to each competing next event from time-to-event distributions (or hazard functions) that may depend on current state and history, processes the earliest event in chronological order, updates the patient's attributes and accumulated cost/QALYs, then re-samples downstream events. Between events the patient "does nothing" from the model's perspective, so DES avoids both the fixed-cycle granularity and the memoryless (first-order Markov) restriction of cohort models.

In an RWE-parameterized DES, the inputs are estimated directly from longitudinal data using the same survival, multistate, and cost-regression methods catalogued elsewhere: time-to-progression or time-to-next-treatment by line and subgroup (parametric or flexible-parametric survival fits), adverse-event hazards, post-progression survival, monthly paid amounts conditional on state and time since the last event (gamma/log-link GLM, winsorized), and utility decrements per event. The simulation propagates those distributions — with their parameter uncertainty — over a lifetime horizon for a large cohort of synthetic patients sampled to match the target population's baseline covariate joint distribution.

Core conceptual distinction

DES is not an estimation method; it is a projection structure into which RWE-derived parameters are loaded. The estimand it produces is a comparative cost-effectiveness contrast (incremental cost per QALY, NMB) under a counterfactual policy or technology, conditional on the causal validity of the inputs. The structural innovation versus a cohort Markov (fixed states, memoryless per-cycle transitions) or a partitioned-survival model (PSM: area under independent PFS/OS curves) is that DES carries individual event history forward, so a patient who has already failed two lines, suffered a grade-3 toxicity, and had a dose reduction can have a genuinely different hazard for the next event than a newly treated patient — without enumerating dozens of "tunnel" states. It also handles competing events natively (the first event scheduled wins; death automatically precludes later progression) and can model resource queues (infusion chairs, transplant organs) that cohort occupancy fractions cannot represent. The price is computational cost (thousands of patients × replications × PSA draws), stochastic (Monte Carlo) output requiring many runs for stable means, and reduced transparency — the model is a program, not an auditable transition matrix or three survival curves.

Pros, cons, and trade-offs

(naming the alternatives explicitly). - vs. cohort Markov (markov-transition-probabilities-rwe): DES represents time-in-state and history-dependent hazards without the combinatorial explosion of tunnel states a Markov needs to relax the memoryless assumption; it also captures resource contention and continuous time. Cost: a Markov matrix is transparent, deterministic given its transitions, trivial to audit, and far cheaper to run. Prefer DES when prior events (number of lines, a past AE, accumulating comorbidity) materially shift downstream hazards and a faithful Markov would need an unmanageable state space. - vs. partitioned-survival model (partitioned-survival-models-rwe): DES makes post-progression pathways, subsequent lines, AE branching, and path-dependent cost/HCRU explicit, rather than leaving them as the residual area between independent PFS and OS curves (a PSM weakness that can produce implausible state occupancy and ignores the PFS–OS dependence). Cost: a PSM is simpler, communicates in two or three curves, and is often sufficient and preferred by HTA reviewers when the horizon is close to observed follow-up. Prefer DES for complex multi-line oncology, transplant, or chronic-disease pathways where independence of PFS/OS or memorylessness is clinically untenable. - vs. direct extrapolation of observed costs/survival (no structural model): "what happened in the data, blown out to 10 years" cannot answer "what if second-line uptake improved or the AE rate fell?" DES supplies the structural what-if while remaining parameterized entirely from RWE. Cost: it layers structural and distributional assumptions on top of the data, which must be justified and validated.

When to use

- Diseases with complex, highly individual pathways: advanced oncology with multiple lines and AE-driven branching; transplant or CAR-T pathways with waitlist/queue dynamics; chronic disease with accumulating comorbidity and repeated hospitalizations. - When history dependence (a prior event changing the next hazard) or competing risks are first-order drivers of cost-effectiveness and a Markov would need an unmanageable state space. - When the decision turns on capacity constraints, queues, or time-dependent resource use that cohort occupancy fractions approximate poorly. - When rich patient-level RWE (ideally linked claims + EHR + registry) supports credible individual-level time-to-event and cost distributions, and the audience (HTA body, payer) will accept a simulation.

When NOT to use — and when it is actively misleading or dangerous

- Sparse or short RWE. If follow-up is too short or events too rare to estimate the many time-to-event distributions and effect modifiers stably, DES becomes an elaborate vehicle for unverifiable extrapolation assumptions wearing the costume of patient-level realism. Prefer a parsimonious Markov/PSM whose few assumptions are visible. - The audience requires an auditable structure. Many HTA submissions still prefer Markov or PSM precisely because a reviewer can trace every number; a stochastic program is harder to interrogate. - The dominant uncertainty is structural, not parametric. If the open question is which states or transitions exist, DES does not resolve it and can hide it inside code that looks precise — the most dangerous failure mode, because Monte Carlo precision is mistaken for evidential strength. - Confounded or immortal-time-contaminated inputs for the comparative effect. DES propagates whatever treatment effect you feed it. A hazard ratio estimated from a prevalent-user or naive cohort will be reproduced with high precision across 10,000 patients — bias laundered through a sophisticated engine. The comparative input must come from an active-comparator new-user / target-trial-emulated analysis, not a raw RWE contrast.

Data-source operational depth

- Claims (FFS or commercial): Excellent for time-to-next-treatment, time-to-discontinuation, time-to-hospitalization, and paid amounts (`days_supply`, NDC, place-of-service, paid columns). Progression is not directly observed — use validated proxies (regimen change, new imaging + diagnosis cluster) and acknowledge the misclassification. Two specific traps: (1) Medicare Advantage person-time lacks fee-for-service claims, so a "no event" stretch in an MA enrollee can be unobserved care, not true event-free time — restrict survival/cost estimation to enrollees with complete FFS (Parts A/B/D) or commercial medical+pharmacy coverage and censor at MA switch; (2) differential competing risks by exposure in elderly claims — if one arm's patients are older/sicker, they die sooner, truncating observed progression and HCRU; estimate cause-specific or Fine-Gray hazards (competing-risks-cause-specific-fine-gray-rwe) and feed competing event distributions to the simulation rather than a single naive time-to-event. - EHR: Rich for progression timing (imaging, pathology, notes), ECOG/labs, and dose reductions, which sharpen both the event hazards and the utility decrements. But observation is visit-driven: sicker patients generate more encounters and therefore more recorded events, inflating apparent event rates. Use inverse-intensity-of-visit weighting or a multistate model that respects irregular observation before exporting hazards to the DES, and treat leaving the system as potentially informative censoring. - Registry: Gold standard for adjudicated progression, response, and vital status — the cleanest source for the clinical event hazards — but typically limited in duration and to a selected population, and usually missing complete cost/HCRU. Link to claims for the economic inputs. - Linked claims–EHR–vital records: The ideal substrate for a credible DES — claims for complete longitudinal cost and utilization, EHR for clinical event timing and severity, registry/death index for endpoints. Beware linkage selection (only the linkable subset is modeled) and order/fill/service date discrepancies that must be reconciled before assigning each patient's event clock. - Immortal time in procedure-defined states. When a state is entered only by surviving to a procedure (transplant, second-line start), the time from index to procedure is "immortal." Model the procedure as a scheduled event with its own time-to-event hazard from index, never as a baseline split — otherwise the procedure arm inherits guaranteed survival time and the cost-effectiveness is biased toward it.

Pre-specify, before running: the event list; the distribution family and covariates/modifiers for each time-to-event and cost process (and exactly which RWE analysis produced each); the history-dependence rules; resource modules if any; random-seed handling; the number of patients, replications, and PSA draws; and the summary measures (mean cost, QALYs, ICER, CEAC, EVPPI). Cross-validate intermediate outputs against the source data — simulated vs. observed distribution of lines of therapy, simulated vs. observed cumulative cost at 24 months, simulated vs. KM survival — before trusting any extrapolated result.

Worked claims/EHR example (advanced solid tumor, 1L → 2L → death with AE-driven history dependence). Population: Medicare FFS + commercial adults initiating first-line systemic therapy for a solid tumor, 2018–2022, with 365 days of continuous, FFS-observable enrollment before the index fill (NDC + `fill_date` + `days_supply`) and linked EHR for ECOG and labs. Index date = first qualifying fill; continuous-enrollment and washout rules establish a clean incident cohort, and MA-only person-time is excluded so absence of claims means absence of events. From these data: (1) time-to-progression on 1L (proxy = regimen switch or new metastatic-workup cluster) fitted with a Weibull AFT model with covariates age, comorbidity score, and baseline ECOG; (2) post-1L AE flag (grade-3+ from diagnosis/procedure algorithms) carried forward as an attribute; (3) time-to-progression on 2L, fitted with the same Weibull but including the prior-AE indicator — observed to shorten 2L TTP (history dependence the DES reproduces by reading each patient's `prior_ae` attribute); (4) competing time-to-death estimated as a cause-specific hazard so that, in the elderly arm, death correctly truncates further progression and cost; (5) monthly cost from a gamma GLM with log link on paid amounts conditional on current line, months since line start, and a recent-hospitalization flag, winsorized at the 99th percentile. Ten thousand synthetic patients are simulated to death or a 10-year horizon, discounting cost and QALYs at 3%. The counterfactual technology improves 1L TTP (HR 0.75) and cuts the grade-3+ AE rate 30% (with linked cost and utility effects). Output: mean life-years, QALYs, lifetime cost, and the ICER vs. standard of care; a 1,000-draw PSA over the joint parameter posterior yields the CEAC, and EVPPI flags 1L TTP and post-progression survival as the highest-value parameters for further RWE collection. Before reporting, the simulated 24-month cumulative cost and the simulated 1L/2L sojourn-time distributions are checked against the observed claims to confirm the engine reproduces the data it was built from.

Interpreting the output

The worked example simulates three patients through a sequence of line-of-therapy events. Patient A (no adverse event history) accumulates 3.0 life-years and 2.10 QALYs; Patient B (prior AE) accumulates 2.0 life-years and 1.22 QALYs; Patient C (early progression) accumulates 0.7 life-years and 0.53 QALYs. Mean across three patients: 1.90 life-years and 1.28 QALYs.

(1) Formal interpretation. DES outputs are individual-level simulated event pathways aggregated across a large number of simulated patients (typically 1,000–100,000). Mean QALYs of 1.28 and mean life-years of 1.90 are estimates of expected values under the model's parameterization. The difference in QALYs between Patient A (2.10) and Patient B (1.22) — a gap of 0.88 QALYs — illustrates history dependence: Patient B's prior adverse event triggered a utility decrement and altered subsequent event timing in a way that a Markov model, which has no memory of past events, could not represent without adding tunnel states. This history dependence is the defining advantage of DES over cohort Markov: each patient carries a full event history as attributes, so any downstream event can be conditioned on prior events without structural changes to the model.

(2) Practical interpretation. When interpreting DES results, the analyst should examine both the aggregate means (the ICER input) and the simulated event-pathway distributions (sojourn times on each line of therapy, event frequencies). If the simulated 24-month cumulative cost distribution does not match the observed claims distribution, the model requires recalibration before the ICER can be trusted. DES results are stochastic — running the same model twice with a different random seed will produce slightly different means. Report the simulation seed and confirm that the number of simulated patients is large enough that mean QALY and cost estimates have stabilized (typically checked by running at 500, 1,000, and 5,000 patients and confirming convergence).

Worked example

Scenario

Imagine three virtual cancer patients each starting first-line (1L) chemotherapy. We want to trace what happens to each one through treatment events — and tally up how many years of healthy life and how many dollars each path costs. Patient A has an uncomplicated 1L course and dies of disease progression on 2L. Patient B has a grade-3 adverse event (AE) on 1L that shortens their subsequent 2L response. Patient C dies on 1L before ever reaching 2L. This tiny three-patient trace shows how DES captures individual pathways — including the history-dependence effect where Patient B's prior AE worsens their 2L outcome.

Dataset

Simulated event log: one row per event per patient (the kind of output a DES engine produces internally, based on time-to-event distributions estimated from real-world claims and EHR data).

Steps

• For each patient, the simulation draws a time to each possible next event (progression, AE, death) from distributions estimated via real-world data, then schedules whichever event comes first.

• Patient A progresses at 1.5 years with no prior AE, moves to 2L, then dies at 3.0 years; life-years = 3.0, QALYs = (1.5 x 0.75) + (1.5 x 0.65) = 1.125 + 0.975 = 2.10.

• Patient B has a grade-3 AE at 0.8 years on 1L; the model flags prior_ae = 1 on their record, which shifts the 2L progression hazard upward so they progress sooner (1.2 years) and die earlier (2.0 years); life-years = 2.0, QALYs = (0.8 x 0.75) + (0.4 x 0.55) + (0.8 x 0.50) = 0.60 + 0.22 + 0.40 = 1.22.

• Patient C dies on 1L at 0.7 years, never reaching 2L; life-years = 0.7, QALYs = 0.7 x 0.75 = 0.53.

• A Markov cohort model would need separate tunnel states for on_1L_post_AE to capture this history effect; DES carries it naturally as a patient attribute.

Result

Mean life-years across 3 patients = (3.0 + 2.0 + 0.7) / 3 = 1.90 years. Mean QALYs = (2.10 + 1.22 + 0.53) / 3 = 1.28 QALYs. Patient B demonstrates history dependence: their prior AE reduced their mean QALYs by 0.88 compared to Patient A (2.10 - 1.22 = 0.88), a downstream cost and quality-of-life effect the simulation captured through the prior_ae attribute — no extra model states required.

Runnable example

import simpy
import numpy as np
import pandas as pd

DISCOUNT = 0.03          # annual discount rate for cost and QALYs
HORIZON_YEARS = 10.0
P_AE_2L = 0.30           # RWE-estimated grade-3+ AE probability on starting 2L (PSA parameter)

def _weibull_time(rng, params, covars):
    """Draw a Weibull AFT event time (years). scale = exp(b0 + beta.x)."""
    lin = params["scale_intercept"] + sum(
        params["beta"].get(c, 0.0) * covars.get(c, 0.0) for c in params["beta"]
    )
    scale = np.exp(lin)
    return float(scale * rng.weibull(params["shape"]))

def _monthly_cost(params, covars):
    """Expected monthly cost from a gamma GLM with log link."""
    lin = params["intercept"] + sum(
        params["beta"].get(c, 0.0) * covars.get(c, 0.0) for c in params["beta"]
    )
    return float(np.exp(lin))

def _disc(amount, t_years):
    return amount / ((1.0 + DISCOUNT) ** t_years)

def patient_process(env, attrs, tte_params, cost_params, results, rng, p_ae_2l=P_AE_2L):
    # ----- Line 1: compete progression vs death; earliest event wins -----
    t_prog1 = _weibull_time(rng, tte_params["tt_prog_1l"], attrs)
    t_death = _weibull_time(rng, tte_params["tt_death"], attrs)
    t_to_first = min(t_prog1, t_death, HORIZON_YEARS - env.now)
    yield env.timeout(t_to_first)
    cost = _disc(_monthly_cost(cost_params["1L"], attrs) * 12 * t_to_first, env.now)
    qaly = _disc(attrs["util_1l"] * t_to_first, env.now)

    if env.now >= HORIZON_YEARS or t_death <= t_prog1:
        results.append({"person_id": attrs["person_id"], "ly": env.now,
                        "qaly": qaly, "cost": cost})
        return

    # ----- Line 2: prior_ae attribute shifts the 2L hazard (history dependence) -----
    attrs = {**attrs}                                  # progressed; carry an AE realization forward
    attrs["prior_ae"] = 1 if rng.random() < p_ae_2l else attrs["prior_ae"]
    t_prog2 = _weibull_time(rng, tte_params["tt_prog_2l"], attrs)   # reads attrs['prior_ae']
    t_death2 = _weibull_time(rng, tte_params["tt_death"], attrs)
    t_to_next = min(t_prog2, t_death2, HORIZON_YEARS - env.now)
    yield env.timeout(t_to_next)
    cost += _disc(_monthly_cost(cost_params["2L"], attrs) * 12 * t_to_next, env.now)
    qaly += _disc(attrs["util_2l"] * t_to_next, env.now)

    results.append({"person_id": attrs["person_id"], "ly": env.now,
                    "qaly": qaly, "cost": cost})

def run_des(baseline: pd.DataFrame, tte_params: dict, cost_params: dict,
            seed: int = 1, p_ae_2l: float = P_AE_2L) -> pd.DataFrame:
    env = simpy.Environment()
    rng = np.random.default_rng(seed)
    results: list[dict] = []
    for _, row in baseline.iterrows():
        env.process(patient_process(env, row.to_dict(), tte_params,
                                    cost_params, results, rng, p_ae_2l))
    env.run(until=HORIZON_YEARS)
    return pd.DataFrame(results)

library(simmer)
library(data.table)

DISCOUNT <- 0.03
HORIZON_YEARS <- 10
P_AE_2L <- 0.30   # RWE-estimated grade-3+ AE probability on starting 2L (PSA parameter)

weibull_time <- function(params, covars) {
  lin <- params$scale_intercept +
    sum(vapply(names(params$beta),
               function(c) params$beta[[c]] * (covars[[c]] %||% 0), numeric(1)))
  scale <- exp(lin)
  scale * rweibull(1, shape = params$shape, scale = 1)  # unit-scale draw, AFT-scaled
}
`%||%` <- function(a, b) if (is.null(a) || is.na(a)) b else a

run_des_r <- function(baseline, tte_params, cost_params, seed = 1L, p_ae_2l = P_AE_2L) {
  set.seed(seed)
  out <- vector("list", nrow(baseline))
  for (i in seq_len(nrow(baseline))) {
    a <- as.list(baseline[i, ])
    # Line 1: progression vs death compete
    t_prog1 <- weibull_time(tte_params$tt_prog_1l, a)
    t_death <- weibull_time(tte_params$tt_death,   a)
    t1 <- min(t_prog1, t_death, HORIZON_YEARS)
    disc1 <- 1 / (1 + DISCOUNT)^0
    cost <- disc1 * exp(cost_params$`1L`$intercept) * 12 * t1
    qaly <- disc1 * a$util_1l * t1
    if (t_death <= t_prog1 || t1 >= HORIZON_YEARS) {
      out[[i]] <- data.table(person_id = a$person_id, ly = t1, qaly = qaly, cost = cost)
      next
    }
    # Line 2: prior_ae attribute shifts the 2L hazard (history dependence)
    a$prior_ae <- if (runif(1) < p_ae_2l) 1 else a$prior_ae
    t_prog2 <- weibull_time(tte_params$tt_prog_2l, a)   # uses a$prior_ae
    t2 <- min(t_prog2, weibull_time(tte_params$tt_death, a), HORIZON_YEARS - t1)
    disc2 <- 1 / (1 + DISCOUNT)^t1
    cost <- cost + disc2 * exp(cost_params$`2L`$intercept) * 12 * t2
    qaly <- qaly + disc2 * a$util_2l * t2
    out[[i]] <- data.table(person_id = a$person_id, ly = t1 + t2,
                           qaly = qaly, cost = cost)
  }
  rbindlist(out)
}