Time-dependent cohort model

This tutorial shows how to extend the cohort replication with age-varying mortality and one-time transition rewards, reproducing the time-dependent Sick-Sicker tutorial of Alarid-Escudero and others (2023). Full script: examples/mdm_cohort_timedep.py.

Specifying age-varying transitions and rewards

The model function returns a per-cycle transition array of shape (n_cycles, n_states, n_states): the Healthy-to-Dead rate follows a US life table, scaled by hazard ratios in the Sick and Sicker states. transition_cost and transition_effect attach a one-time cost of dying, a cost of onset, and a disutility of onset to the flows between states; the engine accrues them on the flow and discounts them with the state rewards.

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")
N = 75

# US all-cause mortality rate, ages 25 to 99 (period life table, total population)
MORT = np.array([
    0.001014, 0.000999, 0.00107, 0.001087, 0.001162, 0.001167, 0.001213, 0.001289, 0.001331,
    0.001375, 0.00142, 0.00149, 0.00155, 0.001616, 0.001657, 0.001747, 0.001902, 0.002052,
    0.002173, 0.002395, 0.002559, 0.002807, 0.003023, 0.003349, 0.003712, 0.004085, 0.00449,
    0.004905, 0.005364, 0.005806, 0.006253, 0.006775, 0.007395, 0.007895, 0.008418, 0.008974,
    0.009666, 0.010456, 0.011384, 0.011838, 0.012667, 0.013593, 0.0147, 0.015732, 0.01734,
    0.018758, 0.020967, 0.022917, 0.024913, 0.026767, 0.029707, 0.032412, 0.035982, 0.039238,
    0.043595, 0.048727, 0.053735, 0.059911, 0.066618, 0.074051, 0.08219, 0.090754, 0.103968,
    0.115093, 0.124341, 0.137872, 0.154177, 0.172393, 0.1941, 0.212654, 0.243752, 0.259087,
    0.287781, 0.316429, 0.339149])

def r2p(rate):
    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"])
    vHD, vS1D, vS2D = r2p(MORT), r2p(MORT * p["hr_S1"]), r2p(MORT * p["hr_S2"])
    prog = r2p(p["r_S1S2"] * p["hr_S1S2_trtB"]) if "B" in intervention else p_S1S2
    P = np.zeros((N, 4, 4))
    P[:, 0, 0], P[:, 0, 1], P[:, 0, 3] = (1 - vHD) * (1 - p_HS1), (1 - vHD) * p_HS1, vHD
    P[:, 1, 0], P[:, 1, 2], P[:, 1, 3] = (1 - vS1D) * p_S1H, (1 - vS1D) * prog, vS1D
    P[:, 1, 1] = (1 - vS1D) * (1 - p_S1H - prog)
    P[:, 2, 2], P[:, 2, 3], P[:, 3, 3] = 1 - vS2D, vS2D, 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"]
    eff = np.array([p["u_H"], u_s1, p["u_S2"], 0.0])
    tc = np.zeros((4, 4)); tc[0, 1] = p["ic_HS1"]; tc[0, 3] = tc[1, 3] = tc[2, 3] = p["ic_D"]
    te = np.zeros((4, 4)); te[0, 1] = -p["du_HS1"]
    return CohortSpec(P, cost, eff, transition_cost=tc, transition_effect=te)

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

Reproducing the published base case

Running the model once at the article’s point estimates should reproduce its published table exactly, the same check used in the time-independent replication. Age-varying mortality lowers life expectancy, so every intervention accrues fewer quality-adjusted life-years than in the time-independent model. Intervention A is dominated; Intervention B and Intervention AB stay on the frontier.

from heormodel.cea import icer_table
base = dict(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,
            du_HS1=0.01, ic_HS1=1000.0, ic_D=2000.0)
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	116374.45	18.88	NaN	NaN	NaN	ND
Intervention B	202536.46	20.20	86162.01	1.32	65287.60	ND
Intervention A	218789.36	19.64	NaN	NaN	NaN	D
Intervention AB	296300.33	21.10	93763.87	0.90	104460.51	ND

Intervention B has an incremental cost-effectiveness ratio near 65,000 per quality-adjusted life-year and Intervention AB near 104,000, matching the published results. The cost-effectiveness analysis is identical to the time-independent case.

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