Cohort state-transition model

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 compare standard of care with treatment A (improves the Sick-state utility), treatment B (slows progression from Sick to Sicker), and their combination 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, so hr_S1 and hr_S2 multiply the death rate before r2p converts it to a per-cycle probability, matching how the source article specifies uncertainty.

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 math and rewards are wired 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 about 73,000 per quality-adjusted life-year), 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	NaN	NaN	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.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
icer_table(outcomes).round(2)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
Standard of care	152191.46	20.55	NaN	NaN	NaN	ND
Intervention B	259880.16	21.97	107688.69	1.42	75816.26	ND
Intervention A	285011.81	21.36	NaN	NaN	NaN	D
Intervention AB	379374.62	22.96	119494.47	0.98	121595.31	ND

The mean frontier matches the deterministic ranking above. 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.

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.