Markov vs microsimulation models

This tutorial shows how to build one Sick-Sicker-style model twice, as a MarkovModel cohort trace and as a MicrosimModel individual simulation, from the same transition rates and rewards, in order to cross-validate the two engines against each other and show what a microsimulation represents that a cohort model averages away. It builds on the cohort state-transition tutorial and the microsimulation engine tutorial. The full script is examples/markov_vs_microsim.py.

Specifying one model for two engines

The model is progressive, with four states, Healthy, Sick, Sicker, and Dead, over 40 annual cycles. Onset, progression, and death are competing annual hazards; each cycle’s transition probabilities come from the same hazards function for both engines, so any difference between the two engines’ results is Monte Carlo noise, not a coding difference between two separately written models. Frailty z multiplies the progression and mortality hazards, and the cohort runs at z = 1, the population average.

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

STATES = ("H", "S1", "S2", "D")
N_CYCLES = 40
BASE = dict(r_HS1=0.12, r_S1S2=0.10, r_HD=0.010, hr_S1=3.0, hr_S2=10.0,
            c_H=1_000.0, c_S1=4_000.0, c_S2=15_000.0, u_H=1.0, u_S1=0.75, u_S2=0.5)
draws = pd.DataFrame([BASE], index=pd.RangeIndex(1, name="iteration"))
COST = np.array([BASE["c_H"], BASE["c_S1"], BASE["c_S2"], 0.0])
EFF = np.array([BASE["u_H"], BASE["u_S1"], BASE["u_S2"], 0.0])

def hazards(p, state, z):
    haz = np.zeros((len(state), 4))
    h, s1, s2 = state == 0, state == 1, state == 2
    haz[h, 1] = p["r_HS1"]; haz[h, 3] = p["r_HD"] * z[h]
    haz[s1, 2] = p["r_S1S2"] * z[s1]; haz[s1, 3] = p["r_HD"] * p["hr_S1"] * z[s1]
    haz[s2, 3] = p["r_HD"] * p["hr_S2"] * z[s2]
    return haz

def rows_from(haz):
    total = haz.sum(axis=1)
    with np.errstate(invalid="ignore", divide="ignore"):
        share = np.where(total[:, None] > 0, haz / total[:, None], 0.0)
    return (1.0 - np.exp(-total))[:, None] * share

def cohort_model(p, intervention):
    P = rows_from(hazards(p, np.arange(4), np.ones(4)))
    P[np.arange(4), np.arange(4)] += 1.0 - P.sum(axis=1)
    P[3] = [0.0, 0.0, 0.0, 1.0]
    return CohortSpec(P, COST, EFF)

def make_pop(var):
    def population(rng, n):
        return pd.DataFrame({"z": np.ones(n) if var == 0 else rng.gamma(1 / var, var, n)})
    return population

def transition_probabilities(p, intervention, state, attrs, rng):
    probs = rows_from(hazards(p, state, attrs["z"].to_numpy()))
    probs[np.arange(len(state)), state] += 1.0 - probs.sum(axis=1)
    probs[state == 3] = [0.0, 0.0, 0.0, 1.0]
    return probs

def state_rewards(p, intervention, state, attrs):
    return COST[state], EFF[state]

def microsim(n, var):
    return MicrosimModel.discrete(
        states=STATES, transition_probabilities=transition_probabilities,
        state_rewards=state_rewards, population=make_pop(var),
        n_individuals=n, interventions=["Standard of care"], n_cycles=N_CYCLES,
        discount_rate=0.03, cycle_correction="half_cycle")

cohort = MarkovModel(states=STATES, interventions=("Standard of care",),
    transitions_and_rewards=cohort_model,
    n_cycles=N_CYCLES, initial_state="H", discount_rate=0.03, cycle_correction="half_cycle")
cohort.evaluate(draws).summary().round(3)

	cost	qaly
intervention
Standard of care	78613.52	11.74

Cross-validating the two engines

With a memoryless, homogeneous population, the microsimulation mean should approach the cohort trace as the population grows; if it did not, one of the two models would have a coding error. The check below runs the microsimulation at increasing population sizes and compares each to the cohort trace computed above.

rows = [dict(n=n, **run_psa(microsim(n, 0.0), draws, seed=1, sequential=True)
             .outcomes.summary().loc["Standard of care"]) for n in (2_000, 10_000, 40_000)]
pd.DataFrame(rows).set_index("n").round(3)

	cost	qaly
n
2000	80088.470	11.780
10000	78186.563	11.763
40000	78949.204	11.769

The cohort trace gives 78,614 dollars and 11.740 quality-adjusted life-years. At 40,000 individuals the microsimulation lands within half a percent of both, and the gap shrinks further as the population grows. The two engines agree because they are, mathematically, the same model.

Adding heterogeneous frailty

This section assigns each individual a frailty z, drawn from a Gamma distribution with mean 1 and variance 0.5, that multiplies the progression and mortality hazards. The mean hazard across the population is unchanged, so the cohort model, which tracks only the average person, produces the same result as before. The microsimulation does not, because it tracks each individual rather than the average.

het = run_psa(microsim(40_000, 0.5), draws, seed=1, sequential=True).outcomes.summary()
het.round(3)

	cost	qaly
intervention
Standard of care	80445.793	12.665

Quality-adjusted life-years rise by about 8% and cost by about 2% against the cohort trace, even though the mean rates are identical. The reason is frailty selection: individuals with higher frailty die sooner, so the survivors are increasingly drawn from the lower-frailty part of the population, and this survivor group lives longer on average than a cohort model, which tracks a single average-frailty person throughout, would predict. This gap is a property of the model, not a bug, and it is why a microsimulation is the right model type to use when risk is heterogeneous.

Representing history with duration groups

The same machinery also carries history, which the frailty example above did not need. This section makes mortality in the sick states rise 8% for every year already spent sick. A cohort model can only represent that by adding a tunnel state for each year already spent sick, since a Markov transition cannot depend on how long an individual has been in a state; the microsimulation instead tracks a duration_groups counter per individual and reads it directly in the transition function, without restructuring the state space.

def hist_transition(p, intervention, state, attrs, rng):
    haz = hazards(p, state, attrs["z"].to_numpy())
    haz[:, 3] *= np.where(np.isin(state, (1, 2)), 1.0 + 0.08 * attrs["tis"].to_numpy(), 1.0)
    probs = rows_from(haz)
    probs[np.arange(len(state)), state] += 1.0 - probs.sum(axis=1)
    probs[state == 3] = [0.0, 0.0, 0.0, 1.0]
    return probs

hist = MicrosimModel.discrete(
    states=STATES, transition_probabilities=hist_transition,
    state_rewards=state_rewards, population=make_pop(0.5),
    n_individuals=40_000, interventions=["Standard of care"], n_cycles=N_CYCLES,
    discount_rate=0.03, cycle_correction="half_cycle",
    duration_groups={"tis": ("S1", "S2")})
run_psa(hist, draws, seed=1, sequential=True).outcomes.summary().round(3)

	cost	qaly
intervention
Standard of care	63616.929	11.893

The trade is plain. The cohort model is faster and exact for its assumptions. The microsimulation costs iterations but represents heterogeneity and history the cohort averages away. Neither is more correct in general; they answer under different assumptions.

Next: the discrete-event simulation engine adds a scarce, shared resource that couples entities together, a dependence neither the cohort nor the independent microsimulation carries.