Microsimulation replication

Reproduce the Sick-Sicker microsimulation of Krijkamp and others (2018) with MicrosimModel.

Open In Colab

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. Duration enters the death rate through a multiplier on the hazard. With dur full years already spent sick, the Sick-to-Dead probability is \(p_{SD} = 1 - e^{-\text{hr}\cdot r_{HD}\,(1 + (\text{dur}+1)\,\text{rp})}\), where rp is the per-year proportional increase. The duration count does not reset when a person progresses from Sick to Sicker. The engine therefore tracks it with duration_groups, a counter over a set of states, which time_in_state cannot do because it counts time in one state only. 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

Both arms run at the base-case parameters with common random numbers. We simulate 50,000 individuals rather than the source’s 100,000, which gives up 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 62,667 dollars and 15.28 quality-adjusted life-years (QALYs) without treatment, 117,455 dollars and 15.79 QALYs with treatment, and an incremental cost-effectiveness ratio of 107,986 dollars per QALY. The random-number sampler differs from the source, so the totals agree within Monte Carlo error rather than to the dollar, and they 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.

Back to top