Microsimulation engine

This tutorial shows how to build an individual-level model with MicrosimModel, the model type to use when outcomes depend on a patient’s history or on heterogeneity between patients, neither of which a cohort average can carry. It runs the same run_psa, cost-effectiveness, and value-of-information calls as the Markov cohort model tutorial, so the two are worth reading as a pair. This page walks through examples/microsim.py step by step.

Specifying individual-level transition probabilities and utilities

The model has three states, Healthy, Sick, and Dead, over a 30-year horizon. Two features rule out a single cohort transition matrix: each individual carries a frailty attribute that scales their sick-state cost and mortality, and the Sick to Dead risk rises with time already spent sick. The engine maintains time_in_state and passes it, alongside the sampled attributes, to the transition_probabilities and state_rewards functions. Both also receive the intervention name as their second argument.

import numpy as np
import pandas as pd

STATES = ("Healthy", "Sick", "Dead")
BACKGROUND_MORTALITY = 0.01
WORSENING = 0.05  # extra Sick -> Dead risk per year already spent sick

def population(rng, n):
    return pd.DataFrame({"frailty": rng.lognormal(mean=0.0, sigma=0.3, size=n)})

def transition_probabilities(params, intervention, state, attrs, rng):
    n = len(state)
    probs = np.zeros((n, 3))
    on_tx = intervention == "Treatment"
    p_hs = params["p_hs"] * (params["rr_tx"] if on_tx else 1.0)
    healthy = state == 0
    probs[healthy, 1] = p_hs
    probs[healthy, 2] = BACKGROUND_MORTALITY
    probs[healthy, 0] = 1.0 - p_hs - BACKGROUND_MORTALITY
    sick = state == 1
    frailty = attrs["frailty"].to_numpy()[sick]
    time_sick = attrs["time_in_state"].to_numpy()[sick]
    p_sd = np.clip(params["p_sd"] * frailty * (1.0 + WORSENING * time_sick), 0.0, 1.0)
    probs[sick, 2] = p_sd
    probs[sick, 1] = 1.0 - p_sd
    probs[state == 2, 2] = 1.0  # Dead is absorbing
    return probs

def state_rewards(params, intervention, state, attrs):
    n = len(state)
    cost, qaly = np.zeros(n), np.zeros(n)
    frailty = attrs["frailty"].to_numpy()
    tx_cost = params["c_treat"] if intervention == "Treatment" else 0.0
    healthy = state == 0
    cost[healthy] = params["c_well"] + tx_cost
    qaly[healthy] = 1.0
    sick = state == 1
    cost[sick] = params["c_sick"] * frailty[sick] + tx_cost
    qaly[sick] = params["u_sick"]
    return cost, qaly

Configuring once, evaluating on draws

The model structure goes to the engine when you build it; the draw matrix goes to run_psa. The two interventions are bare names, and the model functions branch on the intervention name (intervention == "Treatment") to decide whether treatment is on. The two interventions share one population through common random numbers (the engine default), so the incremental result reflects the treatment effect rather than sampling noise.

from heormodel.models import MicrosimModel
from heormodel.params import Beta, Fixed, Gamma, ParameterSet
from heormodel.run import SeedManager, run_psa

seeds = SeedManager(20260704)
parameters = ParameterSet({
    "p_hs": Beta.from_mean_se(0.08, 0.02),
    "p_sd": Beta.from_mean_se(0.12, 0.03),
    "rr_tx": Beta.from_mean_se(0.60, 0.05),
    "u_sick": Beta.from_mean_se(0.65, 0.05),
    "c_well": Fixed(500.0),
    "c_sick": Gamma.from_mean_se(9_000.0, 1_500.0),
    "c_treat": Gamma.from_mean_se(2_000.0, 300.0),
})
draws = parameters.sample(256, seed=seeds.generator())

engine = MicrosimModel.discrete(
    states=STATES, transition_probabilities=transition_probabilities,
    state_rewards=state_rewards,
    population=population, n_individuals=800,
    interventions=["Standard care", "Treatment"],
    n_cycles=30,
)
outcomes = run_psa(engine, draws, seed=seeds.entropy, sequential=True).outcomes
outcomes.summary().round(3)

	cost	qaly
intervention
Standard care	36811.928	10.535
Treatment	57495.416	12.461

Runs are parallel by default and give the same numbers either way: each iteration is seeded by its position in the draw matrix, so a parallel and a sequential run agree exactly. This small run passes sequential=True.

Branching on the intervention name versus using decision levers

An intervention carries a name and, optionally, parameter decision levers. Which to use depends on how the arms differ:

• Branch on the name when the arms differ in structure or in which model function runs, as here: treatment changes a utility and a cost, so the model reads intervention == "Treatment". The comparison is not a number that lives in the ParameterSet.
• Use an Intervention with decision levers for a numeric scenario knob the model already reads as a parameter, such as a server count Intervention("Expanded", {"n_servers": 2}). The discrete-event tutorial does this for clinic capacity.

Do not encode an arm as a fake parameter, a flag like {"on_treatment": 1.0}. Such a flag is not part of the ParameterSet, so evppi_ranking and the deterministic sensitivity analyses never see it, and it would clutter the ranking with a value that is not a real source of uncertainty. Branching on the name keeps the ranking below to genuine parameters.

Analyzing cost-effectiveness and value of information

From here nothing is engine-specific: the same icer_table, evpi, and evppi_ranking calls used for a cohort model apply to this simulation’s Outcomes. At a willingness-to-pay threshold near the base-case incremental cost-effectiveness ratio, the decision is genuinely uncertain, so the expected value of perfect information is positive and the per-parameter ranking below is informative.

from heormodel.cea import icer_table
icer_table(outcomes).round(3)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
Standard care	36811.928	10.535	NaN	NaN	NaN	ND
Treatment	57495.416	12.461	20683.488	1.927	10735.334	ND

from heormodel.voi import evpi, evppi_ranking

WTP = 11_000.0
print(f"EVPI at WTP {WTP:,.0f}: {evpi(outcomes, WTP):,.1f}")
evppi_ranking(outcomes, draws, WTP).round(1)

EVPI at WTP 11,000: 2,390.5

c_treat    1510.4
rr_tx      1238.9
p_hs       1022.2
c_sick      309.5
p_sd         96.9
u_sick       36.4
c_well        0.0
Name: evppi, dtype: float64

The treatment cost and its effect rr_tx lead the ranking: they determine whether treatment’s incremental cost-effectiveness ratio stays below the threshold.

Validating against the cohort closed form

With constant transitions and no heterogeneity, a microsimulation must converge to the exact cohort calculation as the population grows; this section checks that it does, since a mismatch here would point to a bug rather than sampling noise. The engine and the cohort share the same discounting and half-cycle convention, so the only expected gap is Monte Carlo error.

P = np.array([[0.80, 0.15, 0.05], [0.0, 0.90, 0.10], [0.0, 0.0, 1.0]])
cost_vec, eff_vec = np.array([1_000.0, 3_000.0, 0.0]), np.array([1.0, 0.6, 0.0])

def cohort(horizon=40, rate=0.03):
    p, total_c, total_e = np.array([1.0, 0.0, 0.0]), 0.0, 0.0
    for c in range(horizon + 1):
        w = 0.5 if c in (0, horizon) else 1.0
        total_c += w * (1 + rate) ** -c * float(p @ cost_vec)
        total_e += w * (1 + rate) ** -c * float(p @ eff_vec)
        p = p @ P
    return total_c, total_e

check = MicrosimModel.discrete(
    states=STATES,
    transition_probabilities=lambda params, intervention, state, attrs, rng: P[state],
    state_rewards=lambda params, intervention, state, attrs: (cost_vec[state], eff_vec[state]),
    population=60_000, interventions=["cohort"], n_cycles=40,
)
sim = run_psa(check, draws.iloc[[0]], seed=1).outcomes.summary().loc["cohort"]
ref_cost, ref_eff = cohort()
pd.DataFrame({
    "microsim": [sim["cost"], sim["qaly"]],
    "cohort": [ref_cost, ref_eff],
}, index=["cost", "qaly"]).assign(rel_error=lambda d: (d["microsim"] / d["cohort"] - 1))

	microsim	cohort	rel_error
cost	19405.229563	19329.141622	0.003936
qaly	7.084982	7.048307	0.005203

The two agree to under 1% at 60,000 individuals, and the gap shrinks with the population.

Switching to a continuous-time clock

MicrosimModel.continuous is set up the same way. Instead of per-cycle probabilities, event_times returns sampled times to each competing event; the engine races them, advances to the earliest, and accrues continuously between events, valuing each segment through state_reward_rates (per-year flows rather than per-cycle amounts). On constant hazards it reproduces the exponential cohort solution. See the engines concept page for how the two clocks compare.

Next: Markov vs microsimulation models builds this microsimulation and a cohort twin from the same rates, cross-validates them, and shows what the individual model carries that the cohort averages away. ```