Cohort state-transition model

Build a cohort state-transition model with MarkovModel and validate it against a published Sick-Sicker result.

Keywords

cohort state-transition model, Markov model, Markov cohort model, state-transition model, Sick-Sicker, cost-effectiveness modeling, health economics, Python

This tutorial shows how to build a cohort state-transition model with MarkovModel and validate it against a published result, using the introductory Sick-Sicker analysis of Alarid-Escudero and others (2023). Reproducing a known answer first, before trusting the model on a new question, is the standard check for any cost-effectiveness model. Full script: examples/mdm_cohort.py; companion replications: replication gallery.

Specifying the four-state model

The model tracks four states, Healthy, Sick, Sicker, and Dead, over 75 annual cycles starting at age 25. Four interventions compete: standard of care, Intervention A (improves the Sick-state utility), Intervention B (slows progression from Sick to Sicker), and the combination Intervention AB. The transitions_and_rewards function below returns each intervention’s transition matrix and per-state payoffs from a single parameter row.

Transition rates, not probabilities, are the quantities sampled and adjusted by hazard ratios. A hazard ratio scales a rate multiplicatively, and a constant rate $r$ becomes a per-cycle probability over a cycle of length $\Delta t$ through

\[p = 1 - e^{-r\,\Delta t}.\]

That conversion is the r2p function; hr_S1 and hr_S2 multiply the death rate before it is converted, matching how the source article specifies uncertainty. The engine then propagates the cohort: with the state occupancy written as a row vector $m_t$ holding the proportion in each state at cycle $t$, one cycle advances it by the transition matrix $P$,

\[m_{t+1} = m_t P,\]

starting from the whole cohort in Healthy. Each cycle’s expected cost and quality-adjusted life-years (QALYs) are the occupancy weighted by the per-state payoff vectors, discounted and summed over the horizon with a half-cycle correction.

import numpy as np
import pandas as pd
from heormodel.models import CohortSpec, MarkovModel

STATES = ("H", "S1", "S2", "D")
INTERVENTIONS = ("Standard of care", "Intervention A", "Intervention B", "Intervention AB")

def r2p(rate):  # rate to per-cycle probability
    return 1.0 - np.exp(-np.asarray(rate))

def model(p, intervention):
    p_HS1, p_S1H, p_S1S2 = r2p(p["r_HS1"]), r2p(p["r_S1H"]), r2p(p["r_S1S2"])
    p_HD, p_S1D, p_S2D = r2p(p["r_HD"]), r2p(p["r_HD"] * p["hr_S1"]), r2p(p["r_HD"] * p["hr_S2"])
    prog = r2p(p["r_S1S2"] * p["hr_S1S2_trtB"]) if "B" in intervention else p_S1S2
    P = np.zeros((4, 4))
    P[0, 0], P[0, 1], P[0, 3] = (1 - p_HD) * (1 - p_HS1), (1 - p_HD) * p_HS1, p_HD
    P[1, 0], P[1, 2], P[1, 3] = (1 - p_S1D) * p_S1H, (1 - p_S1D) * prog, p_S1D
    P[1, 1] = (1 - p_S1D) * (1 - p_S1H - prog)
    P[2, 2], P[2, 3], P[3, 3] = 1 - p_S2D, p_S2D, 1.0
    add = {"Standard of care": 0.0, "Intervention A": p["c_trtA"],
           "Intervention B": p["c_trtB"], "Intervention AB": p["c_trtA"] + p["c_trtB"]}[intervention]
    cost = np.array([p["c_H"], p["c_S1"] + add, p["c_S2"] + add, 0.0])
    u_s1 = p["u_trtA"] if intervention in ("Intervention A", "Intervention AB") else p["u_S1"]
    return CohortSpec(P, cost, np.array([p["u_H"], u_s1, p["u_S2"], 0.0]))

engine = MarkovModel(states=STATES, interventions=INTERVENTIONS, transitions_and_rewards=model,
                            n_cycles=75, initial_state="H", cycle_correction="simpson")

Reproducing the published base case

Before sampling any uncertainty, running the model once at the article’s point estimates checks that the transition matrix and rewards are specified correctly: the result should match the published table exactly. It does. Intervention A is dominated; the frontier runs standard of care, then Intervention B (incremental cost-effectiveness ratio, or ICER, about 73,000 per quality-adjusted life-year, or QALY), then Intervention AB (about 126,000).

from heormodel.cea import icer_table
base = dict(r_HD=0.002, r_HS1=0.15, r_S1H=0.5, r_S1S2=0.105, hr_S1=3.0, hr_S2=10.0,
            hr_S1S2_trtB=0.6, c_H=2000.0, c_S1=4000.0, c_S2=15000.0, c_trtA=12000.0,
            c_trtB=13000.0, u_H=1.0, u_S1=0.75, u_S2=0.5, u_trtA=0.95)
draws0 = pd.DataFrame([base], index=pd.RangeIndex(1, name="iteration"))
icer_table(engine.evaluate(draws0)).round(2)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
Standard of care	151579.87	20.71	NaN	NaN	NaN	ND
Intervention B	259100.41	22.18	107520.54	1.47	72987.64	ND
Intervention A	284804.51	21.50	25704.10	-0.69	NaN	D
Intervention AB	378875.20	23.14	119774.79	0.95	125763.79	ND

Running the probabilistic sensitivity analysis

The article’s uncertainty distributions for each parameter become a ParameterSet; gamma distributions bound rates and costs at zero, beta distributions bound utilities to the unit interval, and log-normal distributions keep hazard ratios positive. From there the analysis is the same as for any other model. The run below passes sequential=True because 1,000 iterations finish quickly; larger runs are parallel by default.

from heormodel.params import Beta, Gamma, LogNormal, ParameterSet
from heormodel.report import format_icer_table
from heormodel.run import SeedManager, run_psa
from heormodel.voi import evpi

params = ParameterSet({
    "r_HD": Gamma(20, 1/10000), "r_HS1": Gamma(30, 1/200), "r_S1H": Gamma(60, 1/120),
    "r_S1S2": Gamma(84, 1/800), "hr_S1": LogNormal(np.log(3), 0.01),
    "hr_S2": LogNormal(np.log(10), 0.02), "hr_S1S2_trtB": LogNormal(np.log(0.6), 0.02),
    "c_H": Gamma(100, 20.0), "c_S1": Gamma(177.8, 22.5), "c_S2": Gamma(225, 66.7),
    "c_trtA": Gamma(73.5, 163.3), "c_trtB": Gamma(86.2, 150.8),
    "u_H": Beta(200, 3), "u_S1": Beta(130, 45), "u_S2": Beta(230, 230), "u_trtA": Beta(300, 15),
})
draws = params.sample(1000, seed=SeedManager(20260705).generator())
outcomes = run_psa(engine, draws, sequential=True).outcomes
format_icer_table(outcomes)

	Cost	Effect	Incremental cost	Incremental effect	ICER	Status
Intervention
Standard of care	152,191 (122,343, 186,644)	20.55 (18.79, 22.53)				ND
Intervention B	259,880 (199,867, 321,559)	21.97 (20.30, 23.78)	107,689 (71,709, 152,397)	1.42 (1.20, 1.64)	75,816 (52,209, 103,586)	ND
Intervention A	285,012 (224,052, 353,900)	21.36 (19.68, 23.17)	25,132 (-18,161, 66,377)	-0.61 (-0.98, -0.20)		D
Intervention AB	379,375 (289,409, 477,240)	22.96 (21.43, 24.55)	119,494 (81,515, 161,215)	0.98 (0.62, 1.40)	121,595 (79,587, 187,942)	ND

Because outcomes holds 1,000 parameter draws, the table reports a 95% uncertainty interval for every estimate, written as the point estimate followed by its 2.5th and 97.5th percentiles in parentheses. format_icer_table rounds costs and the ratio to whole units and effects to two decimals for reading; icer_table returns the same numbers unrounded, in separate _lo and _hi columns, for further computation. Incremental cost and effect are differences against the comparator, the cheapest frontier intervention still above the row in cost, formed within every draw and then summarized by percentiles rather than by differencing the two interventions’ separate intervals.

The ICER is a frontier quantity, so it is blank for a dominated intervention, and the status column carries the dominance label instead. Intervention A is dominated: it costs about 25,000 more than Intervention B for about 0.6 fewer QALYs, the negative incremental effect that marks the dominance. The mean frontier matches the deterministic ranking above: standard of care, then Intervention B at about 76,000 per QALY (95% interval about 52,000 to 104,000), then Intervention AB at about 122,000 (about 80,000 to 188,000). The interval on the ratio is wider for Intervention AB because its incremental effect over Intervention B, about 0.98 QALYs against 1.42, is smaller, so a given spread in incremental cost produces a wider range of ratios.

Whether that ranking holds at a given threshold in every draw, rather than only on average, is what the expected value of perfect information (EVPI) measures next:

print(f"EVPI at WTP 100,000: {evpi(outcomes, 100_000.0):,.0f}")

EVPI at WTP 100,000: 2,757

EVPI is positive at a threshold of 100,000 because Intervention AB, the most expensive option, is the best choice in some draws but not others: the uncertainty in whether it is worth adopting has a real cost, and resolving it would be worth exactly this much per person. The value of information tutorial takes this same model further: which parameter that uncertainty attaches to, and whether a proposed study to resolve it is worth its cost.

Reference: Alarid-Escudero F, Krijkamp EM, Enns EA, Yang A, Hunink MGM, Pechlivanoglou P, Jalal H. An introductory tutorial on cohort state-transition models in R using a cost-effectiveness analysis example. Medical Decision Making. 2023;43(1):3-20.

Keywords: cohort state-transition model, Markov model, Markov cohort model, state-transition model, Sick-Sicker, cost-effectiveness modeling, health economics, Python

--- title: Cohort state-transition model description: "Build a cohort state-transition model with `MarkovModel` and validate it against a published Sick-Sicker result." keywords: [cohort state-transition model, Markov model, Markov cohort model, state-transition model, Sick-Sicker, cost-effectiveness modeling, health economics, Python] ---  <a href="https://colab.research.google.com/github/pedroliman/heormodel/blob/main/docs/_notebooks/mdm-cohort.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"></a>  This tutorial shows how to build a cohort state-transition model with `MarkovModel` and validate it against a published result, using the introductory Sick-Sicker analysis of Alarid-Escudero and others (2023). Reproducing a known answer first, before trusting the model on a new question, is the standard check for any cost-effectiveness model. Full script: [`examples/mdm_cohort.py`](https://github.com/pedroliman/heormodel/blob/main/examples/mdm_cohort.py); companion replications: [replication gallery](replication-gallery.qmd). ## Specifying the four-state model The model tracks four states, Healthy, Sick, Sicker, and Dead, over 75 annual cycles starting at age 25. Four interventions compete: standard of care, Intervention A (improves the Sick-state utility), Intervention B (slows progression from Sick to Sicker), and the combination Intervention AB. The `transitions_and_rewards` function below returns each intervention's transition matrix and per-state payoffs from a single parameter row. Transition rates, not probabilities, are the quantities sampled and adjusted by hazard ratios. A hazard ratio scales a rate multiplicatively, and a constant rate $r$ becomes a per-cycle probability over a cycle of length $\Delta t$ through $$p = 1 - e^{-r\,\Delta t}.$$ That conversion is the `r2p` function; `hr_S1` and `hr_S2` multiply the death rate before it is converted, matching how the source article specifies uncertainty. The engine then propagates the cohort: with the state occupancy written as a row vector $m_t$ holding the proportion in each state at cycle $t$, one cycle advances it by the transition matrix $P$, $$m_{t+1} = m_t P,$$ starting from the whole cohort in Healthy. Each cycle's expected cost and quality-adjusted life-years (QALYs) are the occupancy weighted by the per-state payoff vectors, discounted and summed over the horizon with a half-cycle correction. ```{python} import numpy as np import pandas as pd from heormodel.models import CohortSpec, MarkovModel STATES = ("H", "S1", "S2", "D") INTERVENTIONS = ("Standard of care", "Intervention A", "Intervention B", "Intervention AB") def r2p(rate): # rate to per-cycle probability return 1.0 - np.exp(-np.asarray(rate)) def model(p, intervention): p_HS1, p_S1H, p_S1S2 = r2p(p["r_HS1"]), r2p(p["r_S1H"]), r2p(p["r_S1S2"]) p_HD, p_S1D, p_S2D = r2p(p["r_HD"]), r2p(p["r_HD"] * p["hr_S1"]), r2p(p["r_HD"] * p["hr_S2"]) prog = r2p(p["r_S1S2"] * p["hr_S1S2_trtB"]) if "B" in intervention else p_S1S2 P = np.zeros((4, 4)) P[0, 0], P[0, 1], P[0, 3] = (1 - p_HD) * (1 - p_HS1), (1 - p_HD) * p_HS1, p_HD P[1, 0], P[1, 2], P[1, 3] = (1 - p_S1D) * p_S1H, (1 - p_S1D) * prog, p_S1D P[1, 1] = (1 - p_S1D) * (1 - p_S1H - prog) P[2, 2], P[2, 3], P[3, 3] = 1 - p_S2D, p_S2D, 1.0 add = {"Standard of care": 0.0, "Intervention A": p["c_trtA"], "Intervention B": p["c_trtB"], "Intervention AB": p["c_trtA"] + p["c_trtB"]}[intervention] cost = np.array([p["c_H"], p["c_S1"] + add, p["c_S2"] + add, 0.0]) u_s1 = p["u_trtA"] if intervention in ("Intervention A", "Intervention AB") else p["u_S1"] return CohortSpec(P, cost, np.array([p["u_H"], u_s1, p["u_S2"], 0.0])) engine = MarkovModel(states=STATES, interventions=INTERVENTIONS, transitions_and_rewards=model, n_cycles=75, initial_state="H", cycle_correction="simpson") ``` ## Reproducing the published base case Before sampling any uncertainty, running the model once at the article's point estimates checks that the transition matrix and rewards are specified correctly: the result should match the published table exactly. It does. Intervention A is dominated; the frontier runs standard of care, then Intervention B (incremental cost-effectiveness ratio, or ICER, about 73,000 per quality-adjusted life-year, or QALY), then Intervention AB (about 126,000). ```{python} from heormodel.cea import icer_table base = dict(r_HD=0.002, r_HS1=0.15, r_S1H=0.5, r_S1S2=0.105, hr_S1=3.0, hr_S2=10.0, hr_S1S2_trtB=0.6, c_H=2000.0, c_S1=4000.0, c_S2=15000.0, c_trtA=12000.0, c_trtB=13000.0, u_H=1.0, u_S1=0.75, u_S2=0.5, u_trtA=0.95) draws0 = pd.DataFrame([base], index=pd.RangeIndex(1, name="iteration")) icer_table(engine.evaluate(draws0)).round(2) ``` ## Running the probabilistic sensitivity analysis The article's uncertainty distributions for each parameter become a `ParameterSet`; gamma distributions bound rates and costs at zero, beta distributions bound utilities to the unit interval, and log-normal distributions keep hazard ratios positive. From there the analysis is the same as for any other model. The run below passes `sequential=True` because 1,000 iterations finish quickly; larger runs are parallel by default. ```{python} from heormodel.params import Beta, Gamma, LogNormal, ParameterSet from heormodel.report import format_icer_table from heormodel.run import SeedManager, run_psa from heormodel.voi import evpi params = ParameterSet({ "r_HD": Gamma(20, 1/10000), "r_HS1": Gamma(30, 1/200), "r_S1H": Gamma(60, 1/120), "r_S1S2": Gamma(84, 1/800), "hr_S1": LogNormal(np.log(3), 0.01), "hr_S2": LogNormal(np.log(10), 0.02), "hr_S1S2_trtB": LogNormal(np.log(0.6), 0.02), "c_H": Gamma(100, 20.0), "c_S1": Gamma(177.8, 22.5), "c_S2": Gamma(225, 66.7), "c_trtA": Gamma(73.5, 163.3), "c_trtB": Gamma(86.2, 150.8), "u_H": Beta(200, 3), "u_S1": Beta(130, 45), "u_S2": Beta(230, 230), "u_trtA": Beta(300, 15), }) draws = params.sample(1000, seed=SeedManager(20260705).generator()) outcomes = run_psa(engine, draws, sequential=True).outcomes format_icer_table(outcomes) ``` Because `outcomes` holds 1,000 parameter draws, the table reports a 95% uncertainty interval for every estimate, written as the point estimate followed by its 2.5th and 97.5th percentiles in parentheses. `format_icer_table` rounds costs and the ratio to whole units and effects to two decimals for reading; `icer_table` returns the same numbers unrounded, in separate `_lo` and `_hi` columns, for further computation. Incremental cost and effect are differences against the comparator, the cheapest frontier intervention still above the row in cost, formed within every draw and then summarized by percentiles rather than by differencing the two interventions' separate intervals. The ICER is a frontier quantity, so it is blank for a dominated intervention, and the status column carries the dominance label instead. Intervention A is dominated: it costs about 25,000 more than Intervention B for about 0.6 fewer QALYs, the negative incremental effect that marks the dominance. The mean frontier matches the deterministic ranking above: standard of care, then Intervention B at about 76,000 per QALY (95% interval about 52,000 to 104,000), then Intervention AB at about 122,000 (about 80,000 to 188,000). The interval on the ratio is wider for Intervention AB because its incremental effect over Intervention B, about 0.98 QALYs against 1.42, is smaller, so a given spread in incremental cost produces a wider range of ratios. Whether that ranking holds at a given threshold in every draw, rather than only on average, is what the expected value of perfect information (EVPI) measures next: ```{python} print(f"EVPI at WTP 100,000: {evpi(outcomes, 100_000.0):,.0f}") ``` EVPI is positive at a threshold of 100,000 because Intervention AB, the most expensive option, is the best choice in some draws but not others: the uncertainty in whether it is worth adopting has a real cost, and resolving it would be worth exactly this much per person. The [value of information](voi.qmd) tutorial takes this same model further: which parameter that uncertainty attaches to, and whether a proposed study to resolve it is worth its cost. Reference: Alarid-Escudero F, Krijkamp EM, Enns EA, Yang A, Hunink MGM, Pechlivanoglou P, Jalal H. An introductory tutorial on cohort state-transition models in R using a cost-effectiveness analysis example. Medical Decision Making. 2023;43(1):3-20.