Parameter inputs from data

Turn point values, an exported draw matrix, or a weighted posterior sample into the parameter draw matrix run_psa expects.

Open In Colab

Not every analysis starts from ParameterSet distributions. Some run at one base-case set of point values, some use a draw matrix exported from another tool, and some use a weighted posterior sample. This tutorial shows how to turn each of those three sources into a parameter draw matrix with single_draw (via ParameterSet.at_means), read_draws, and resample_posterior. It walks through examples/parameter_inputs.py step by step. Bring your own outputs covers the same problem on the outcome side. The rule for all three sources is the same: once the result is a draw matrix with an iteration index and numeric columns, run_psa treats it like any other sample.

Specifying a model to run

Two interventions compete for a chronic disease: standard care and a new drug. Cost and quality-adjusted life-years (QALYs) per iteration are simple functions of three parameters.

import numpy as np
import pandas as pd
from heormodel.models import Outcomes
from heormodel.params import (
    Beta, Gamma, ParameterSet, read_draws, resample_posterior, single_draw,
)
from heormodel.run import run_psa

def model(draws: pd.DataFrame) -> Outcomes:
    base_qaly = 8.0
    effect_drug = base_qaly + draws["u_gain"] * draws["p_response"] * 10
    costs = pd.DataFrame(
        {"Standard care": 40_000.0, "New drug": 40_000.0 + draws["c_drug"]},
        index=draws.index,
    )
    effects = pd.DataFrame(
        {"Standard care": base_qaly, "New drug": effect_drug}, index=draws.index,
    )
    return Outcomes.from_wide(costs, effects)

params = ParameterSet(
    {
        "p_response": Beta.from_mean_se(0.35, 0.05),
        "c_drug": Gamma.from_mean_se(12_000, 1_500),
        "u_gain": Beta.from_mean_se(0.12, 0.03),
    },
    correlation={("p_response", "u_gain"): 0.4},
)

Running a base case at one set of point values

ParameterSet.at_means returns the analytic means as a one-row draw matrix at iteration 0, giving the deterministic run that a probabilistic analysis is usually checked against. It is equivalent to single_draw(params.means().to_dict()). Call single_draw directly when the point values come from somewhere other than a ParameterSet, such as a published base case.

base = params.at_means()
run_psa(model, base).outcomes.summary().round(2)
cost qaly
intervention
Standard care 40000.0 8.00
New drug 52000.0 8.42

The summary holds the deterministic base case, evaluated once at the parameter means rather than over the sampled distribution the later sections use.

Reading a draw matrix from another tool

read_draws validates a CSV path or a DataFrame as a draw matrix. It uses an explicit iteration column when present and otherwise assigns a fresh index. A non-numeric column raises an error before it reaches the model, so a malformed export fails immediately rather than deep inside the analysis. The code below writes a sampled matrix to CSV to stand in for an external tool’s export, then reads it back as a real export would be.

external = params.sample(2_000, seed=1)
external.to_csv("external_draws.csv")
csv_draws = read_draws("external_draws.csv", iteration="iteration")
run_psa(model, csv_draws).outcomes.summary().round(2)
cost qaly
intervention
Standard care 40000.00 8.00
New drug 52002.43 8.42

The summary is the same probabilistic analysis a native ParameterSet sample would give. read_draws changes only where the draws come from, not how they are analyzed.

Resampling a weighted posterior

A calibration or a bootstrap may return parameter rows with a weight column instead of an already-uniform sample. resample_posterior draws whole rows with replacement, each row \(i\) with probability proportional to its weight,

\[\Pr(\text{row } i) = \frac{w_i}{\sum_j w_j},\]

so any correlation between parameters in the posterior is preserved. It then drops the weight column and returns a plain draw matrix for run_psa. Resampling to an n larger than the input adds no information. It only smooths Monte Carlo noise in downstream expectations.

grid = params.sample(500, seed=2)
grid["weight"] = np.exp(4.0 * grid["p_response"])  # favor higher response
posterior = resample_posterior(grid, n=2_000, seed=7)
print(f"grid mean p_response {grid['p_response'].mean():.3f}, "
      f"resample mean {posterior['p_response'].mean():.3f}")
run_psa(model, posterior).outcomes.summary().round(2)
grid mean p_response 0.348, resample mean 0.360
cost qaly
intervention
Standard care 40000.00 8.00
New drug 51898.82 8.44

The resampled mean of p_response is higher than the unweighted grid mean because the weights favor higher values. The reweighting in turn raises the new drug’s expected QALYs.

The calibration workflow produces such a posterior with abc_calibrate and mixes it with literature draws through mix_draws. Next: deterministic sensitivity analysis sweeps parameters one at a time through the same analysis.

Back to top