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. The cost-effectiveness and value-of-information analyses that follow treat the mixed draws as they would draws from a single source. 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, so no cycle loop is needed. 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}$. 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 an onset rate of 0.02 and a progression rate of 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.0505
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 rates affect prevalence, so the data constrain their combination more tightly than either rate alone.

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

np.float64(0.965)

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, which preserves the correlation within the posterior. It resamples each source independently, which introduces no correlation between sources.

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 or from the 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 and the 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	cost_lo	cost_hi	effect	effect_lo	effect_hi	inc_cost	inc_cost_lo	inc_cost_hi	inc_effect	inc_effect_lo	inc_effect_hi	icer	icer_lo	icer_hi	status
intervention
Standard care	28180.865	18770.374	39681.601	18.237	17.294	19.193	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	ND
Treatment	50409.647	37699.904	64809.015	18.649	17.703	19.618	22228.782	16762.494	28763.469	0.413	0.231	0.661	53878.472	36377.065	89224.905	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. The expected value of perfect information (EVPI) bounds the value of resolving all parameter uncertainty. 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,154.6

rr_progression    656.9
c_treat           284.5
u_sick            258.2
c_sick            103.7
progression        37.4
onset               6.2
Name: evppi, dtype: float64

The treatment effect rr_progression leads the ranking. The calibrated rates rank at the bottom. Calibration left them little remaining uncertainty, so learning more about them would barely change the decision.

Recording provenance across sources

capture_run takes a draw_sources argument that gives each parameter’s source, calibration or literature, and the run report records it. The two sources carry different kinds of uncertainty, so a reviewer needs 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-10-26T06:33:44+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.6 |
| 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. The cost-effectiveness and value-of-information analyses that follow treat the mixed draws as they would draws from a single source. 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, so no cycle loop is needed. 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}$. 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 an onset rate of 0.02 and a progression rate of 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 rates affect prevalence, so the data constrain their combination more tightly than either rate 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, which preserves the correlation within the posterior. It resamples each source independently, which introduces no correlation between sources. ```{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 or from the 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 and the 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. The expected value of perfect information (EVPI) bounds the value of resolving all parameter uncertainty. 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. Calibration left them little remaining uncertainty, so learning more about them would barely change the decision. ## Recording provenance across sources `capture_run` takes a `draw_sources` argument that gives each parameter's source, calibration or literature, and the run report records it. The two sources carry different kinds of uncertainty, so a reviewer needs 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.