Microsimulation replication

This tutorial shows how to reproduce the Sick-Sicker microsimulation of Krijkamp and others (2018) with MicrosimModel, building on the microsimulation engine tutorial. Full script: examples/mdm_microsim.py.

Modeling duration-dependent risk and a treatment effect modifier

Two features rule out a cohort average: mortality rises with the number of consecutive years spent in the Sick and Sicker states, and a per-individual effect modifier scales the treatment utility. The duration clock keeps running across the Sick to Sicker progression, so the engine tracks it with duration_groups, a counter over a set of states that time_in_state (one exact state) cannot express. The counter is 0 on the first sick cycle; mortality reads it as dur + 1 and utility reads it directly.

import numpy as np
import pandas as pd
from heormodel.models import MicrosimModel
from heormodel.run import run_psa

STATES = ("H", "S1", "S2", "D")

def population(rng, n):
    return pd.DataFrame({"x": rng.uniform(0.95, 1.05, n)})

def transition_probabilities(p, intervention, state, attrs, rng):
    n = len(state); probs = np.zeros((n, 4))
    mult = 1.0 + (attrs["dur"].to_numpy() + 1.0) * p["rp_S1S2"]
    r_HD = -np.log(1 - p["p_HD"])
    p_S1D = 1 - np.exp(-p["rr_S1"] * r_HD * mult)
    p_S2D = 1 - np.exp(-p["rr_S2"] * r_HD * mult)
    h = state == 0
    probs[h] = np.array([1 - p["p_HS1"] - p["p_HD"], p["p_HS1"], 0.0, p["p_HD"]])
    s1 = state == 1
    probs[s1, 0], probs[s1, 2], probs[s1, 3] = p["p_S1H"], p["p_S1S2"], p_S1D[s1]
    probs[s1, 1] = 1 - p["p_S1H"] - p["p_S1S2"] - p_S1D[s1]
    s2 = state == 2
    probs[s2, 2], probs[s2, 3] = 1 - p_S2D[s2], p_S2D[s2]
    probs[state == 3, 3] = 1.0
    return probs

def state_rewards(p, intervention, state, attrs):
    n = len(state); cost = np.zeros(n); util = np.zeros(n)
    dur, x, tx = attrs["dur"].to_numpy(), attrs["x"].to_numpy(), intervention == "Treatment"
    tx_cost = p["c_Trt"] if tx else 0.0
    util[state == 0] = p["u_H"]; cost[state == 0] = p["c_H"]
    s1 = state == 1
    cost[s1] = p["c_S1"] + tx_cost
    util[s1] = x[s1] * (p["u_Trt"] - dur[s1] * p["ru_S1S2"]) if tx else p["u_S1"]
    s2 = state == 2
    cost[s2] = p["c_S2"] + tx_cost; util[s2] = p["u_S2"]
    return cost, util

Running the base case and comparing interventions

The base-case parameters drive both arms through common random numbers, with 50,000 individuals rather than the source’s 100,000, which trades some Monte Carlo precision for a faster run.

from heormodel.cea import icer_table
base = dict(p_HD=0.005, p_HS1=0.15, p_S1H=0.5, p_S1S2=0.105, rr_S1=3.0, rr_S2=10.0,
            rp_S1S2=0.2, c_H=2000.0, c_S1=4000.0, c_S2=15000.0, c_Trt=12000.0,
            u_H=1.0, u_S1=0.75, u_S2=0.5, u_Trt=0.95, ru_S1S2=0.03)
engine = MicrosimModel.discrete(
    states=STATES, transition_probabilities=transition_probabilities,
    state_rewards=state_rewards, population=population,
    n_individuals=50_000,
    interventions=["No Treatment", "Treatment"],
    n_cycles=30, discount_rate=0.03,
    cycle_correction="none", duration_groups={"dur": ("S1", "S2")})
outcomes = run_psa(
    engine, pd.DataFrame([base], index=pd.RangeIndex(1, name="iteration")),
    seed=1, sequential=True,
).outcomes
icer_table(outcomes).round(2)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
No Treatment	62691.63	15.34	NaN	NaN	NaN	ND
Treatment	117468.28	15.85	54776.65	0.51	107824.62	ND

The published Table 3 reports no-treatment 62,667 dollars and 15.28 quality-adjusted life-years (QALYs), treatment 117,455 dollars and 15.79 QALYs, and an incremental cost-effectiveness ratio of 107,986 dollars per QALY. The sampler differs from the source, so the totals agree within Monte Carlo error rather than to the dollar and converge as the population grows.

Reference: Krijkamp EM, Alarid-Escudero F, Enns EA, Jalal HJ, Hunink MGM, Pechlivanoglou P. Microsimulation modeling for health decision sciences using R: a tutorial. Medical Decision Making. 2018;38(3):400-422.