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 before using 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, which matches how the source article specifies uncertainty. The engine then propagates the cohort. Writing the state occupancy as a row vector $m_t$ that holds 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, we run the model once at the article’s point estimates to check that the transition matrix and rewards are specified correctly. The result should match the published table exactly. It does. Intervention A is dominated, and 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 keep rates and costs above zero, beta distributions keep utilities between zero and one, 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, which is the next cheaper intervention on the frontier. They are computed within every draw and then summarized by percentiles, not 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 and yields about 0.6 fewer QALYs. 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 is smaller (about 0.98 QALYs over Intervention B, against 1.42 for Intervention B over standard of care). A given spread in incremental cost then produces a wider range of ratios.

The ranking above holds on average, but not necessarily in every draw. The expected value of perfect information (EVPI) measures what that disagreement across draws costs at a given threshold:

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 in others. Resolving that uncertainty would be worth the printed amount per person. The value of information tutorial takes this model further. It asks which parameter the uncertainty comes from 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 before using 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, which matches how the source article specifies uncertainty. The engine then propagates the cohort. Writing the state occupancy as a row vector $m_t$ that holds 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, we run the model once at the article's point estimates to check that the transition matrix and rewards are specified correctly. The result should match the published table exactly. It does. Intervention A is dominated, and 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 keep rates and costs above zero, beta distributions keep utilities between zero and one, 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, which is the next cheaper intervention on the frontier. They are computed within every draw and then summarized by percentiles, not 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 and yields about 0.6 fewer QALYs. 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 is smaller (about 0.98 QALYs over Intervention B, against 1.42 for Intervention B over standard of care). A given spread in incremental cost then produces a wider range of ratios. The ranking above holds on average, but not necessarily in every draw. The expected value of perfect information (EVPI) measures what that disagreement across draws costs at a given threshold: ```{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 in others. Resolving that uncertainty would be worth the printed amount per person. The [value of information](voi.qmd) tutorial takes this model further. It asks which parameter the uncertainty comes from 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.