Calibration workflow

This tutorial shows how to mix two parameter sources in one analysis: some parameters are calibrated to observed data, and the rest come from the literature. Both feed one probabilistic sensitivity analysis that flows through cost-effectiveness and value-of-information analysis exactly as a single-source analysis would. It builds on the full pipeline tutorial and walks through examples/calibration_workflow.py step by step. Calibration needs the calibration extra: uv pip install 'heormodel[calibration]'.

Specifying the disease model

A three-state continuous-time Markov chain runs over Healthy, Sick, and Dead from age 40. Two rates are unknown and will be calibrated: the Healthy to Sick onset hazard and the excess Sick to Dead progression hazard. Treatment scales the progression hazard by rr_progression, a literature parameter, buying time in the Sick state. Discounted state occupancy is the closed-form integral of the transition matrix, e0 @ inv(r I - Q), so no cycle loop is needed.

import numpy as np
import pandas as pd
from scipy.linalg import expm

DISCOUNT, MORTALITY, START_AGE = 0.03, 0.012, 40

def generator(onset, progression):
    return np.array([
        [-(onset + MORTALITY), onset, MORTALITY],
        [0.0, -(MORTALITY + progression), MORTALITY + progression],
        [0.0, 0.0, 0.0],
    ])

def natural_history(params):
    """Calibration simulator: rates -> Sick prevalence at ages 50 and 70."""
    q = generator(params["onset"], params["progression"])
    out = {}
    for label, age in (("prev_age50", 50), ("prev_age70", 70)):
        p = np.array([1.0, 0.0, 0.0]) @ expm(q * (age - START_AGE))
        out[label] = float(p[1] / (p[0] + p[1]))
    return out

def occupancy(onset, progression):
    """Discounted years in Healthy and Sick, vectorized over iterations."""
    n = len(onset)
    q = np.zeros((n, 3, 3))
    q[:, 0, 0] = -(onset + MORTALITY)
    q[:, 0, 1] = onset
    q[:, 0, 2] = MORTALITY
    q[:, 1, 1] = -(MORTALITY + progression)
    q[:, 1, 2] = MORTALITY + progression
    m = np.linalg.inv(DISCOUNT * np.eye(3) - q)
    return m[:, 0, 0], m[:, 0, 1]

Calibrating the rates

abc_calibrate runs approximate Bayesian computation (sequential Monte Carlo) against the registry targets and returns an iteration-indexed posterior draw matrix. The targets here were generated from onset 0.02 and progression 0.05, so calibration has a known answer to recover.

import os
os.environ["ABC_LOG_LEVEL"] = "WARNING"  # quiet pyabc's per-population logging

from heormodel.calibrate import abc_calibrate
from heormodel.params import Uniform
from heormodel.run import SeedManager

N = 2_000
seeds = SeedManager(20260704)
rng_abc, rng_lit, rng_mix = seeds.spawn(3)

calibration = abc_calibrate(
    natural_history,
    priors={"onset": Uniform(0.005, 0.04), "progression": Uniform(0.01, 0.1)},
    observed={"prev_age50": 0.147, "prev_age70": 0.283},
    population_size=200,
    max_populations=6,
    n_posterior=N,
    seed=int(rng_abc.integers(2**32)),
)
posterior = calibration.posterior
posterior.mean().round(4)

onset          0.0201
progression    0.0503
dtype: float64

The onset and progression posteriors are strongly correlated: both raise prevalence, so the data constrain their combination more tightly than either alone.

round(posterior["onset"].corr(posterior["progression"]), 3)

np.float64(0.963)

Sampling literature parameters

Utilities and costs come from the literature as mean and standard-error distributions instead, with a negative correlation between the sick-state utility and its cost.

from heormodel.params import Beta, Gamma, ParameterSet

literature = ParameterSet(
    {
        "u_sick": Beta.from_mean_se(0.60, 0.05),
        "c_sick": Gamma.from_mean_se(8_000, 1_500),
        "c_treat": Gamma.from_mean_se(4_000, 500),
        "rr_progression": Beta.from_mean_se(0.70, 0.06),
    },
    correlation={("u_sick", "c_sick"): -0.3},
)
lit_draws = literature.sample(N, seed=rng_lit)

Mixing the two sources

mix_draws joins the sources into one matrix. It resamples whole rows within each source, so the posterior’s joint correlation survives, and resamples sources independently, so no correlation is invented across them.

from heormodel.params import mix_draws

draws = mix_draws(posterior, lit_draws, n=N, seed=rng_mix)
list(draws.columns)

['onset', 'progression', 'u_sick', 'c_sick', 'c_treat', 'rr_progression']

Running the model and analyzing the outcomes

From here the workflow does not depend on where a column came from: the decision model runs over the mixed draws exactly as it would over a single-source draw matrix, and the cost-effectiveness and value-of-information analyses that follow treat every column alike.

from heormodel.models import Outcomes
from heormodel.run import run_psa

def disease_model(d):
    _, sick_std = occupancy(d["onset"].to_numpy(), d["progression"].to_numpy())
    healthy, sick_treat = occupancy(
        d["onset"].to_numpy(), d["progression"].to_numpy() * d["rr_progression"].to_numpy()
    )
    costs = pd.DataFrame({
        "Standard care": sick_std * d["c_sick"].to_numpy(),
        "Treatment": sick_treat * (d["c_sick"].to_numpy() + d["c_treat"].to_numpy()),
    }, index=d.index)
    effects = pd.DataFrame({
        "Standard care": healthy + sick_std * d["u_sick"].to_numpy(),
        "Treatment": healthy + sick_treat * d["u_sick"].to_numpy(),
    }, index=d.index)
    return Outcomes.from_wide(costs, effects)

outcomes = run_psa(disease_model, draws, sequential=True).outcomes
from heormodel.cea import icer_table
icer_table(outcomes).round(3)

	cost	effect	inc_cost	inc_effect	icer	status
intervention
Standard care	28249.718	18.229	NaN	NaN	NaN	ND
Treatment	50543.939	18.643	22294.221	0.414	53830.73	ND

from heormodel.voi import evpi, evppi_ranking

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

EVPI at WTP 50,000: 1,173.8

rr_progression    665.8
c_treat           282.8
u_sick            252.9
c_sick            121.2
onset              28.2
progression        17.8
Name: evppi, dtype: float64

The treatment effect rr_progression leads the ranking. The calibrated rates rank at the bottom: tight calibration left them little residual uncertainty, so learning them further would barely change the decision.

Recording provenance across sources

capture_run takes a draw_sources map so the run report records where each parameter came from, calibration or literature, which matters here since the two sources carry different kinds of uncertainty and a reviewer will want to tell them apart.

from heormodel.report import capture_run

record = capture_run(
    seed=seeds,
    outcomes=outcomes,
    draw_sources={
        **{name: "ABC posterior" for name in posterior.columns},
        **{name: "literature (mean/SE)" for name in lit_draws.columns},
    },
    note="Calibration workflow tutorial.",
)
print(record.to_markdown("heormodel calibration workflow run report"))

# heormodel calibration workflow run report

- **Created:** 2026-07-11T20:18:05+00:00
- **Note:** Calibration workflow tutorial.
- **Root seed entropy:** 20260704
- **Iterations:** 2000
- **Interventions:** Standard care, Treatment

## Draw sources

| Parameter | Source |
|---|---|
| onset | ABC posterior |
| progression | ABC posterior |
| u_sick | literature (mean/SE) |
| c_sick | literature (mean/SE) |
| c_treat | literature (mean/SE) |
| rr_progression | literature (mean/SE) |

## Software versions

| Package | Version |
|---|---|
| python | 3.12.3 |
| heormodel | 0.7.1 |
| numpy | 2.5.0 |
| scipy | 1.18.0 |
| pandas | 3.0.3 |
| joblib | 1.5.3 |
| scikit-learn | 1.9.0 |

Next: parameter inputs from data covers the base-case, imported, and posterior draw matrices that feed the same analysis.