Calibration workflow

Mix calibrated parameters with literature parameters in one probabilistic analysis using mix_draws.

Keywords

model calibration, Bayesian calibration, calibrating a health economic model, posterior, prior, cost-effectiveness modeling, health economics, Python

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, and the cost-effectiveness and value-of-information analyses that follow treat the mixed draws exactly as they would a single-source analysis. 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, extending time in the Sick state. Discounted state occupancy has a closed form: for a chain with generator matrix $Q$ and discount rate $r$, the expected discounted years spent in each state, starting from the initial distribution $e_0$, are $e_0 (rI - Q)^{-1}$, so no cycle loop is needed. The occupancy function evaluates this for every draw at once.

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 observed calibration 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.0200
progression    0.0501
dtype: float64

Both posterior means are close to the onset and progression rates that generated the targets, 0.02 and 0.05, so calibration recovers the known answer.

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.969)

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 is preserved, and resamples sources independently, so no correlation is introduced 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']

One matrix now holds both sources, so the decision model reads every parameter, calibrated and literature, from a single set of columns.

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	28252.612	18.262	NaN	NaN	NaN	ND
Treatment	50484.827	18.674	22232.215	0.411	54067.921	ND

The value-of-information analysis below asks which parameters the decision is most sensitive to at a willingness-to-pay threshold of 50,000. Expected value of perfect information (EVPI) bounds the value of resolving all parameter uncertainty, and the ranking attributes that value to individual parameters by their expected value of partial perfect information.

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,138.9

rr_progression    634.3
c_treat           269.3
u_sick            241.9
c_sick            107.1
progression        21.9
onset               4.4
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 argument giving each parameter’s source, 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-19T11:14:31+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.4 |
| 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.

Keywords: model calibration, Bayesian calibration, calibrating a health economic model, posterior, prior, cost-effectiveness modeling, health economics, Python

--- title: Calibration workflow description: "Mix calibrated parameters with literature parameters in one probabilistic analysis using `mix_draws`." keywords: [model calibration, Bayesian calibration, calibrating a health economic model, posterior, prior, cost-effectiveness modeling, health economics, Python] ---  <a href="https://colab.research.google.com/github/pedroliman/heormodel/blob/main/docs/_notebooks/calibration-workflow.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"></a>  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, and the cost-effectiveness and value-of-information analyses that follow treat the mixed draws exactly as they would a single-source analysis. It builds on the [full pipeline tutorial](full-pipeline.qmd) and walks through [`examples/calibration_workflow.py`](https://github.com/pedroliman/heormodel/blob/main/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, extending time in the Sick state. Discounted state occupancy has a closed form: for a chain with generator matrix $Q$ and discount rate $r$, the expected discounted years spent in each state, starting from the initial distribution $e_0$, are $e_0 (rI - Q)^{-1}$, so no cycle loop is needed. The `occupancy` function evaluates this for every draw at once. ```{python} 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 observed calibration 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. ```{python} 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) ``` Both posterior means are close to the onset and progression rates that generated the targets, 0.02 and 0.05, so calibration recovers the known answer. The onset and progression posteriors are strongly correlated: both raise prevalence, so the data constrain their combination more tightly than either alone. ```{python} round(posterior["onset"].corr(posterior["progression"]), 3) ``` ## 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. ```{python} 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 is preserved, and resamples sources independently, so no correlation is introduced across them. ```{python} from heormodel.params import mix_draws draws = mix_draws(posterior, lit_draws, n=N, seed=rng_mix) list(draws.columns) ``` One matrix now holds both sources, so the decision model reads every parameter, calibrated and literature, from a single set of columns. ## 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. ```{python} 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) ``` The value-of-information analysis below asks which parameters the decision is most sensitive to at a willingness-to-pay threshold of 50,000. Expected value of perfect information (EVPI) bounds the value of resolving all parameter uncertainty, and the ranking attributes that value to individual parameters by their expected value of partial perfect information. ```{python} 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) ``` 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` argument giving each parameter's source, 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. ```{python} 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")) ``` Next: [parameter inputs from data](parameter-inputs.qmd) covers the base-case, imported, and posterior draw matrices that feed the same analysis.