---
title: Markov vs microsimulation models
description: "Build one model as both a cohort trace and an individual simulation from the same rates, to cross-validate the engines and show what heterogeneity adds."
---
<!-- colab-badge:start -->
<a href="https://colab.research.google.com/github/pedroliman/heormodel/blob/main/docs/_notebooks/markov-vs-microsim.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"></a>
<!-- colab-badge:end -->
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, 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](mdm-cohort.qmd) and the [microsimulation engine tutorial](microsim.qmd). The full script is [`examples/markov_vs_microsim.py`](https://github.com/pedroliman/heormodel/blob/main/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.
Competing hazards convert to a transition-probability row the same way in both engines, which is what `rows_from` computes. From a state with destination hazards $h_1, \dots, h_k$ and total $H = \sum_j h_j$, the probability of leaving within a cycle is $1 - e^{-H}$, split across destinations in proportion to their hazards,
$$P_j = \left(1 - e^{-H}\right)\frac{h_j}{H}.$$
```{python}
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)
```
## 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.
```{python}
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)
```
The cohort trace gives 78,614 dollars and 11.740 quality-adjusted life-years. At 40,000 individuals the microsimulation is 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.
```{python}
het = run_psa(microsim(40_000, 0.5), draws, seed=1, sequential=True).outcomes.summary()
het.round(3)
```
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 an error, and it is why a microsimulation is the right model type to use when risk is heterogeneous.
## Representing history with duration groups
The microsimulation 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.
```{python}
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)
```
The trade-off is clear. The cohort model is faster and exact under its assumptions. The microsimulation costs iterations but represents the heterogeneity and history the cohort averages away. Neither is more correct in general; each is valid under its own assumptions.
Next: the [discrete-event simulation engine](des.qmd) adds a scarce, shared resource that couples entities together, a dependence neither the cohort nor the independent microsimulation carries.
