Health Economic Modeling Methods Using RWE

In plain language

A health economic model is a mathematical structure that uses real-world data to project what happens to a group of patients over many years and what those outcomes cost — information no single study can supply on its own. The most common type is a Markov model, which sorts patients into a small number of health states (for example, Stable, Hospitalized, and Dead), estimates how fast patients move between those states using data from claims or registries, and attaches a dollar cost and a quality-of-life weight to each state. Running the model forward in time produces an incremental cost-effectiveness ratio (ICER) — the extra cost per extra year of good health gained by one treatment compared with another — which is the number health technology assessment bodies use to decide whether a therapy is worth paying for.

Health economic (HE) modeling using real-world evidence

is the practice of parameterizing and running a decision-analytic model — most often a cohort-state-transition (Markov) model, a partitioned-survival model (PSM), or a discrete-event simulation (DES) — from real-world data so that observed event rates, resource use, costs, and survival are projected over a lifetime horizon into an incremental cost-effectiveness ratio (ICER) and net monetary benefit (NMB). The model is the scaffold; RWE supplies the numbers that fill it. Done well, the model inherits the strengths of RWE (real populations, real costs, long follow-up) and disciplines its weaknesses (selection, censoring, unmeasured confounding) through structural assumptions, extrapolation methods, and probabilistic sensitivity analysis (PSA). Done badly, it launders confounded or right-censored RWE into a single deterministic ICER that looks authoritative and is wrong.

Core conceptual distinction

The decision that does the work is which model structure the available RWE can actually support, and that choice flows from the evidence, not the other way around. (1) Cohort Markov assumes the population is fully described by mutually exclusive health states with memoryless transition probabilities per cycle; it is the natural fit when RWE yields incidence-density transition rates between clinically meaningful states (e.g., stable -> progressed -> dead) and the Markov assumption is tenable after enough states/tunnels are added. (2) Partitioned-survival sidesteps transition probabilities entirely: it derives state membership from independently modeled, extrapolated survival curves (e.g., overall survival and progression-free survival from a cohort), so it lives or dies on survival extrapolation beyond observed follow-up and on the internal coherence of the curves (PFS must never exceed OS). (3) DES drops the Markov memorylessness restriction and models time-to-event and individual histories explicitly — the right tool when RWE shows strong time- or history-dependence (e.g., transplant waitlists, repeated hospitalizations) that a cohort Markov cannot represent without an explosion of tunnel states. The estimand is always incremental (treatment vs comparator over a stated horizon and perspective), and three transformations are non-negotiable and routinely botched: converting RWE rates to per-cycle probabilities (p = 1 - exp(-rate cycle_length); never use rate as a probability), applying a half-cycle correction so events are not all charged to cycle ends, and discounting costs and effects to present value at the jurisdiction's rate (commonly 3.0% or 3.5%). The single-transition formula p = 1 - exp(-rate cycle_length) is exact only when a state has one exit; when a state has multiple competing exits, applying it independently per exit misallocates probability (the per-exit probabilities can sum past 1 and ignore the chance of a competing event pre-empting the modeled one). The correct multi-exit conversion uses the matrix exponential of the rate (generator) matrix, P = exp(QΔt) (Welton & Ades 2005), which is the formally consistent map from competing cause-specific rates to a per-cycle transition matrix.

Pros, cons, and trade-offs

- vs trial-data-only economic models: RWE-parameterized models capture real-world adherence, off-protocol resource use, comparator mix, and costs in the target health system, and can extend follow-up far beyond a trial's horizon. Cost: RWE transition rates and survival are confounded and informatively censored in ways a randomized trial is not, so the model inherits that bias unless the upstream RWE is itself causal (active-comparator new-user design, PS/IPTW, g-methods). Prefer RWE inputs for costs, utilities in routine care, and long-horizon survival; prefer trial inputs (or RWE explicitly de-confounded) for the comparative treatment effect. - Cohort Markov vs partitioned-survival: Markov forces an explicit, clinically interpretable transition structure and handles competing risks and multiple absorbing states cleanly, but requires defensible per-cycle transition probabilities and can need many tunnel states to defeat the memoryless assumption. PSM is simpler to fit from Kaplan-Meier-style RWE survival and communicates well, but its independent curve extrapolation can produce incoherent state occupancy and is highly sensitive to the parametric survival family chosen. Prefer Markov when transitions are the natural evidence and competing risks matter; prefer PSM when survival curves are the primary RWE deliverable and the horizon is not far beyond observed data — and always cross-check PSM against a state-transition model (NICE TSD19 expectation). - Deterministic vs probabilistic (PSA): A single deterministic ICER hides parameter uncertainty and is uninterpretable for a reimbursement decision. PSA propagates input distributions (Beta/Dirichlet for probabilities, Gamma/log-normal for costs, Beta for utilities) into a distribution of ICERs and a cost-effectiveness acceptability curve (CEAC). Always run PSA for an HTA submission; deterministic one-way analysis is a supplement (tornado), never the headline.

When to use

Building or critiquing a cost-effectiveness/cost-utility model for an HTA submission (NICE, ICER-US, CDA-AMC (formerly CADTH), G-BA) where lifetime outcomes must be projected from finite RWE; estimating budget impact or value when transitions, costs, or survival are best observed in claims/EHR/registry data rather than a trial; updating a trial-based model with real-world comparator survival, real-world costs, or real-world treatment patterns; or running value-of-information / scenario analyses that require a transparent, re-runnable model.

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

- The comparative effect is confounded and you treat it as causal. Plugging naive RWE hazard ratios (prevalent users, drug-vs-non-user, immortal-time-contaminated cohorts) into the relative-effect parameter produces a precise, fully discounted ICER built on a biased treatment effect — the most dangerous failure because PSA will not reveal it (PSA propagates parameter uncertainty, not structural confounding). Fix it upstream with an active-comparator new-user design and proper confounding control before the number ever reaches the model. - Extrapolation far beyond observed follow-up with no anchor. When 18 months of registry survival is projected to a 40-year horizon, the parametric family (exponential vs Weibull vs generalized gamma vs spline) drives the result more than the data do. Reporting one curve as if it were known is misleading; show the family, justify it with hazard plots and external/general-population mortality constraints, and treat the choice as structural uncertainty. - Forcing a Markov model onto strongly history-dependent RWE. If real-world risk depends on time-in-state or prior events, a memoryless cohort Markov misstates occupancy; either add tunnel states (and accept the parameter burden) or move to DES. Pretending the Markov assumption holds is a structural error PSA cannot detect. - A deterministic point ICER presented to a payer. Without PSA/CEAC, the decision-maker cannot see decision uncertainty; several HTA bodies will reject the submission outright.

Data-source operational depth

- Claims (FFS vs MA vs commercial): Excellent for resource use and costs (PMPM/PMPY by health state), reasonable for coded events that mark transitions, weak for clinical severity and true progression. Transition rates are estimated as events / person-time within a continuously enrolled cohort, then converted to per-cycle probabilities. Failure modes: (i) Medicare Advantage (MA) encounter records lack adjudicated paid amounts, so cost parameters built on MA person-time are unreliable — restrict costing to FFS (Parts A/B/D) person-time or use validated cost-to-charge bridges; (ii) differential competing risks by exposure in elderly claims — if one arm's patients die of other causes faster, naive cause-specific transition rates overstate the at-risk denominator for the modeled event, so use cumulative-incidence / competing-risk transitions, not 1-minus-KM; (iii) immortal time in procedure or progression studies inflates survival in the treated state if state-entry is mis-timed relative to the index date. - EHR: Strong for the clinical variables that define states and progression (labs, stage, ECOG, problem lists) and for utility-relevant severity, but visit-driven and porous — patients who leave the system are differentially censored, so transition rates and survival are biased unless loss to follow-up is modeled and linkage to a death index firms up the absorbing state. Costs are usually charges, not paid amounts; convert before using as economic inputs. - Registry: Best source for adjudicated progression/death and disease severity that drive transitions and survival extrapolation, but typically thin on full cost capture and on out-of-registry care. Link to claims for complete costs and resource use and to vital records to close out mortality. - Linked claims-EHR-vital records: The ideal substrate — EHR severity to define states, claims for paid costs and person-time, vital records for the absorbing state — but linkage selection (only the linkable subset) and date reconciliation across order/fill/service/death dates must be resolved before transitions and time zero are assigned, or the model is parameterized on a non-representative, mis-timed cohort.

Worked claims example

Question: lifetime cost-utility of a new heart-failure (HF) therapy vs standard care, modeled as a 3-state cohort Markov (Stable -> Hospitalized-HF -> Dead) with a 1-month cycle and 20-year horizon, discounted at 3%. (1) Cohort: adults with >=2 HF diagnoses and 365 days of continuous A/B enrollment before `index_date`; restrict costing to FFS person-time (drop MA-only spans where paid amounts are missing). (2) Transitions from RWE: count first hospitalized-HF events and deaths and divide by person-months at risk to get monthly cause-specific rates; because Stable has two competing exits (Hospitalized, Dead), convert the rate matrix to a per-cycle transition matrix with the matrix exponential P = exp(QΔt) rather than applying p = 1 - exp(-rate) per exit, so the Stable->Dead and Stable->Hospitalized transitions share the correct at-risk denominator and the row sums to 1. (3) Costs: aggregate paid amounts from `claim_lines` to PMPM by state (Stable PMPM from outpatient + pharmacy; Hospitalized-HF from the DRG-paid inpatient claim plus the index-month outpatient costs), inflation-adjusted to a common base year. (4) Utilities: assign a per-cycle QALY weight to each state (e.g., Stable 0.80, Hospitalized 0.60 for the event month) from an EQ-5D mapping. (5) Effect: take the treatment's hospitalization and mortality hazard ratios from an active-comparator new-user, PS-weighted RWE analysis (not a naive contrast), and apply them to the standard-care transition probabilities. (6) Run the trace with a half-cycle correction, discount costs and QALYs monthly, and compute incremental cost, incremental QALYs, the ICER, and NMB at a $100,000/QALY threshold. (7) PSA: redraw transition probabilities from Dirichlet, costs from Gamma, utilities from Beta, and the hazard ratios from their sampling distribution across 5,000 iterations; summarize the CEAC across willingness-to-pay thresholds. Report per CHEERS 2022, stating the perspective, horizon, discount rate, and structural assumptions explicitly.

Interpreting the output

The worked example produces an ICER of $50,000/QALY: the treatment arm costs $72,000 and yields 5.20 QALYs; the standard-care arm costs $52,000 and yields 4.80 QALYs. Incremental cost = $20,000; incremental QALYs = 0.40; ICER = $20,000 / 0.40 = $50,000/QALY, below the $100,000/QALY threshold.

(1) Formal interpretation. This result routes through the health-economic modeling parent concept: the right model type (cohort Markov, partitioned-survival, or DES) was selected based on clinical pathway complexity; transition probabilities were sourced from RWE and adjusted for confounding; costs were applied per state and discounted; QALYs were calculated with utility weights from mapping or direct measurement; and PSA characterized decision uncertainty. The ICER of $50,000/QALY is the ratio of incremental quantities from both arms — it is not the average cost per QALY for the treatment arm alone. A CEAC from the PSA completes the interpretation by quantifying the probability that the treatment is cost-effective across a range of willingness-to-pay values.

(2) Practical interpretation. The $50,000/QALY ICER is a routing result that should prompt the reader to understand which model structure generated it and what the key RWE inputs were. A Markov model parameterized with biased RWE transition probabilities will produce a confident-looking ICER that is systematically wrong; a DES that misrepresents event sequencing will similarly mislead. The model structure, data sources, and key assumptions must be reported transparently alongside the ICER so decision-makers can assess how much structural uncertainty — which the PSA does not capture — surrounds the headline number. Consult the sibling entries (Markov, DES, PSA, discounting, scenario analysis) for the mechanics that produced this output.

Worked example

Scenario

A researcher wants to know whether a new heart-failure pill is worth its cost compared with standard care. She builds a 3-state Markov model — Stable, Hospitalized, Dead — with monthly cycles over 10 years. She estimates transition probabilities from a Medicare claims cohort, assigns per-state costs from paid claims amounts, and assigns per-state quality-of-life weights from a published EQ-5D mapping. The new drug reduces hospitalization and death rates by 20% but costs $500 extra per month. She wants to know the ICER.

Dataset

Model structure: one row per health state with the inputs RWE supplies

Steps

• Start with 1,000 simulated patients in the Stable state at month 0.

• Each month, apply transition probabilities: about 12 Stable patients move to Hospitalized, 4 move to Dead, and the rest stay Stable.

• Patients in Hospitalized either recover to Stable, die, or remain (recovery rate ~0.55/month from claims).

• Multiply how many patients are in each state each month by that state's cost and QALY weight, then sum across all 120 months (10 years).

• Discount future costs and QALYs at 3% per year so that a dollar spent next year counts slightly less than one spent today.

• Run the model twice: once with standard-care transition probabilities and once with the new drug's probabilities (both hospitalization and death rates multiplied by 0.80), adding $500/month drug cost in the treatment arm.

• Subtract: incremental cost = total treatment cost minus total standard-care cost; incremental QALYs = total treatment QALYs minus total standard-care QALYs.

• Divide incremental cost by incremental QALYs to get the ICER.

Result

Standard-care arm: discounted total cost ~$52,000, discounted total QALYs ~4.80. Treatment arm: discounted total cost ~$72,000, discounted total QALYs ~5.20. Incremental cost = $20,000; incremental QALYs = 0.40. ICER = $20,000 / 0.40 = $50,000 per QALY — below the common $100,000/QALY threshold, so the treatment looks cost-effective. Note: the transition probabilities come from RWE; if the drug effect (the 20% rate reduction) were taken from a biased observational comparison instead of a carefully designed study, this ICER would be wrong even though all the arithmetic is correct.

Runnable example

import numpy as np

CYCLE_LEN = 1 / 12          # 1-month cycle expressed in years
HORIZON_M = 240             # 20-year horizon in monthly cycles
DISCOUNT_ANNUAL = 0.03      # 3% per year for costs and effects
WTP = 100_000.0             # willingness-to-pay per QALY
STATES = ["stable", "hosp", "dead"]

def rate_to_prob(rate, t=1.0):
    # Single-transition conversion: EXACT ONLY for a state with one exit.
    # With multiple competing exits this misallocates probability (per-exit
    # values can sum past 1 and ignore competing-event pre-emption); use the
    # matrix-exponential conversion below instead.
    return 1.0 - np.exp(-rate * t)

def _expm(M, terms=40):
    # Dependency-free matrix exponential via scaling-and-squaring + Taylor
    # series, so P = exp(Q*dt) stays runnable on numpy alone.
    norm = np.max(np.sum(np.abs(M), axis=1))
    s = max(0, int(np.ceil(np.log2(norm + 1e-12))))
    A = M / (2 ** s)
    E = np.eye(M.shape[0]); term = np.eye(M.shape[0])
    for k in range(1, terms):
        term = term @ A / k
        E = E + term
    for _ in range(s):
        E = E @ E
    return E

def transition_matrix(rates, t=1.0):
    # Correct multi-exit conversion: build the cause-specific generator
    # (rate/intensity) matrix Q with off-diagonal rates and row sums of 0,
    # then P = exp(Q*t) (Welton & Ades 2005). This is row-stochastic by
    # construction and respects competing risks; 'dead' is absorbing.
    r_sh = rates["stable_to_hosp"]
    r_sd = rates["stable_to_dead"]
    r_hd = rates["hosp_to_dead"]
    r_hs = rates["hosp_to_stable"]
    Q = np.array([
        [-(r_sh + r_sd), r_sh,           r_sd],   # from stable
        [r_hs,           -(r_hs + r_hd), r_hd],   # from hosp
        [0.0,            0.0,            0.0 ],    # from dead (absorbing)
    ])
    P = _expm(Q * t)
    return P

def apply_treatment(rates, hr):
    # Apply hazard ratios to the comparator's event rates for the treated arm.
    r = dict(rates)
    r["stable_to_hosp"] *= hr["hosp"]
    r["stable_to_dead"] *= hr["dead"]
    r["hosp_to_dead"]   *= hr["dead"]
    return r

def run_arm(rates, cost_pmpm, utility, extra_cost_pmpm=0.0):
    P = transition_matrix(rates)
    trace = np.zeros((HORIZON_M + 1, 3))
    trace[0] = [1.0, 0.0, 0.0]                       # everyone starts Stable
    for t in range(HORIZON_M):
        trace[t + 1] = trace[t] @ P
    # Half-cycle correction: average start/end occupancy within each cycle.
    occ = (trace[:-1] + trace[1:]) / 2.0
    disc_m = (1 + DISCOUNT_ANNUAL) ** (-(np.arange(HORIZON_M) * CYCLE_LEN))
    cost_vec = np.array([cost_pmpm["stable"] + extra_cost_pmpm,
                         cost_pmpm["hosp"]   + extra_cost_pmpm, 0.0])
    qaly_vec = np.array([utility["stable"], utility["hosp"], 0.0]) * CYCLE_LEN
    total_cost = float((occ @ cost_vec * disc_m).sum())
    total_qaly = float((occ @ qaly_vec * disc_m).sum())
    return total_cost, total_qaly

def evaluate(rates, hr, cost_pmpm, utility, tx_cost_pmpm):
    c0, q0 = run_arm(rates, cost_pmpm, utility)                          # standard care
    c1, q1 = run_arm(apply_treatment(rates, hr), cost_pmpm, utility,     # treatment
                     extra_cost_pmpm=tx_cost_pmpm)
    d_cost, d_qaly = c1 - c0, q1 - q0
    icer = d_cost / d_qaly if d_qaly != 0 else np.inf
    nmb = WTP * d_qaly - d_cost
    return {"d_cost": d_cost, "d_qaly": d_qaly, "icer": icer, "nmb": nmb}

def psa(rates, hr, cost_pmpm, utility, tx_cost_pmpm,
        n_iter=5000, wtp_grid=np.arange(0, 200_001, 10_000), seed=42):
    # Draw correlated transitions (Dirichlet), costs (Gamma), utilities (Beta),
    # and hazard ratios (log-normal) from RWE-implied distributions.
    rng = np.random.default_rng(seed)
    prob_ce = {w: 0 for w in wtp_grid}
    for _ in range(n_iter):
        s = transition_matrix(rates)[0]                                 # Dirichlet on Stable-row probs
        s_draw = rng.dirichlet(s * 200 + 1e-6)
        r = dict(rates)
        r["stable_to_hosp"] = -np.log(1 - s_draw[1])
        r["stable_to_dead"] = -np.log(1 - s_draw[2])
        cd = {k: rng.gamma(shape=(v / (0.2 * v)) ** 2,                  # Gamma, CV ~ 0.2
                           scale=v / ((v / (0.2 * v)) ** 2))
              for k, v in cost_pmpm.items()}
        ud = {k: rng.beta(v * 50, (1 - v) * 50) for k, v in utility.items()}
        hd = {k: float(np.exp(rng.normal(np.log(v), 0.10))) for k, v in hr.items()}
        res = evaluate(r, hd, cd, ud, tx_cost_pmpm)
        for w in wtp_grid:
            if (w * res["d_qaly"] - res["d_cost"]) > 0:
                prob_ce[w] += 1
    return {w: prob_ce[w] / n_iter for w in wtp_grid}                   # CEAC

CYCLE_LEN <- 1 / 12          # 1-month cycle in years
HORIZON_M <- 240L            # 20-year horizon in monthly cycles
DISCOUNT_ANNUAL <- 0.03
WTP <- 100000

# Single-transition conversion: EXACT ONLY for a state with one exit. With
# multiple competing exits this misallocates probability; use the matrix
# exponential of the generator (transition_matrix) instead.
rate_to_prob <- function(rate, t = 1) 1 - exp(-rate * t)

# Dependency-free matrix exponential (scaling-and-squaring + Taylor series),
# so P = exp(Q*t) stays runnable on base R alone.
expm_base <- function(M, terms = 40L) {
  norm <- max(rowSums(abs(M)))
  s <- max(0L, as.integer(ceiling(log2(norm + 1e-12))))
  A <- M / (2 ^ s)
  n <- nrow(M); E <- diag(n); term <- diag(n)
  for (k in seq_len(terms - 1L)) { term <- (term %*% A) / k; E <- E + term }
  for (i in seq_len(s)) E <- E %*% E
  E
}

transition_matrix <- function(rates, t = 1) {
  # Correct multi-exit conversion: build the cause-specific generator (rate)
  # matrix Q (off-diagonal rates, row sums 0), then P = exp(Q*t)
  # (Welton & Ades 2005). Row-stochastic by construction; respects competing
  # risks; 'dead' is absorbing. Rows: from stable / hosp / dead.
  r_sh <- rates[["stable_to_hosp"]]; r_sd <- rates[["stable_to_dead"]]
  r_hd <- rates[["hosp_to_dead"]];   r_hs <- rates[["hosp_to_stable"]]
  Q <- matrix(c(-(r_sh + r_sd), r_sh,           r_sd,
                r_hs,           -(r_hs + r_hd), r_hd,
                0,              0,              0),
              nrow = 3, byrow = TRUE)
  expm_base(Q * t)
}

apply_treatment <- function(rates, hr) {
  rates[["stable_to_hosp"]] <- rates[["stable_to_hosp"]] * hr[["hosp"]]
  rates[["stable_to_dead"]] <- rates[["stable_to_dead"]] * hr[["dead"]]
  rates[["hosp_to_dead"]]   <- rates[["hosp_to_dead"]]   * hr[["dead"]]
  rates
}

run_arm <- function(rates, cost_pmpm, utility, extra_cost_pmpm = 0) {
  P <- transition_matrix(rates)
  trace <- matrix(0, nrow = HORIZON_M + 1, ncol = 3)
  trace[1, ] <- c(1, 0, 0)                     # all start Stable
  for (t in seq_len(HORIZON_M)) trace[t + 1, ] <- trace[t, ] %*% P
  occ <- (trace[-(HORIZON_M + 1), ] + trace[-1, ]) / 2          # half-cycle correction
  disc <- (1 + DISCOUNT_ANNUAL) ^ (-(seq_len(HORIZON_M) - 1) * CYCLE_LEN)
  cost_vec <- c(cost_pmpm[["stable"]] + extra_cost_pmpm,
                cost_pmpm[["hosp"]]   + extra_cost_pmpm, 0)
  qaly_vec <- c(utility[["stable"]], utility[["hosp"]], 0) * CYCLE_LEN
  list(cost = sum((occ %*% cost_vec) * disc),
       qaly = sum((occ %*% qaly_vec) * disc))
}

evaluate <- function(rates, hr, cost_pmpm, utility, tx_cost_pmpm) {
  sc <- run_arm(rates, cost_pmpm, utility)
  tx <- run_arm(apply_treatment(rates, hr), cost_pmpm, utility, tx_cost_pmpm)
  d_cost <- tx$cost - sc$cost; d_qaly <- tx$qaly - sc$qaly
  list(d_cost = d_cost, d_qaly = d_qaly,
       icer = d_cost / d_qaly, nmb = WTP * d_qaly - d_cost)
}

psa_ceac <- function(rates, hr, cost_pmpm, utility, tx_cost_pmpm,
                     n_iter = 5000, wtp_grid = seq(0, 2e5, 1e4)) {
  ce <- numeric(length(wtp_grid))
  for (i in seq_len(n_iter)) {
    rd <- rates
    rd[["stable_to_hosp"]] <- -log(1 - rbeta(1, 12, 988))     # Beta from event/person-time counts
    rd[["stable_to_dead"]] <- -log(1 - rbeta(1, 4, 996))
    cd <- list(stable = rgamma(1, 25, scale = cost_pmpm[["stable"]] / 25),
               hosp   = rgamma(1, 25, scale = cost_pmpm[["hosp"]]   / 25))
    ud <- list(stable = rbeta(1, utility[["stable"]] * 50, (1 - utility[["stable"]]) * 50),
               hosp   = rbeta(1, utility[["hosp"]]   * 50, (1 - utility[["hosp"]])   * 50))
    hd <- list(hosp = exp(rnorm(1, log(hr[["hosp"]]), 0.10)),
               dead = exp(rnorm(1, log(hr[["dead"]]), 0.10)))
    res <- evaluate(rd, hd, cd, ud, tx_cost_pmpm)
    ce <- ce + ((wtp_grid * res$d_qaly - res$d_cost) > 0)
  }
  data.frame(wtp = wtp_grid, prob_cost_effective = ce / n_iter)         # CEAC
}