Time-dependent cohort model

Extend the cohort replication with age-varying mortality and one-time transition rewards.

Open In Colab

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 is time-dependent because the transition matrix changes from cycle to cycle. The death rate follows a US life table. At cycle \(t\), when the cohort has aged to \(a_0 + t\), the Healthy-state death probability is \(p_{HD}(t) = 1 - e^{-\text{mort}(a_0 + t)}\), and hazard ratios scale it in the Sick and Sicker states. The model function returns the whole sequence at once, as a per-cycle transition array of shape (n_cycles, n_states, n_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, and Intervention B and Intervention AB remain 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 116400.16 18.88 NaN NaN NaN ND
Intervention B 202562.23 20.20 86162.07 1.32 65287.64 ND
Intervention A 218815.07 19.64 16252.84 -0.56 NaN D
Intervention AB 296326.10 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 code 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.

Back to top