import numpy as np
import pandas as pd
from heormodel.models import CohortSpec, MarkovModel
from heormodel.params import Beta, Gamma, LogNormal, ParameterSet
from heormodel.run import SeedManager, run_psa
STATES = ("H", "S1", "S2", "D")
INTERVENTIONS = ("Standard of care", "Intervention A", "Intervention B", "Intervention AB")
WTP = 100_000.0
def r2p(rate): # rate to per-cycle probability
return 1.0 - np.exp(-np.asarray(rate))
def model(p, intervention):
p_HS1, p_S1H, p_S1S2 = r2p(p["r_HS1"]), r2p(p["r_S1H"]), r2p(p["r_S1S2"])
p_HD, p_S1D, p_S2D = r2p(p["r_HD"]), r2p(p["r_HD"] * p["hr_S1"]), r2p(p["r_HD"] * p["hr_S2"])
prog = r2p(p["r_S1S2"] * p["hr_S1S2_trtB"]) if "B" in intervention else p_S1S2
P = np.zeros((4, 4))
P[0, 0], P[0, 1], P[0, 3] = (1 - p_HD) * (1 - p_HS1), (1 - p_HD) * p_HS1, p_HD
P[1, 0], P[1, 2], P[1, 3] = (1 - p_S1D) * p_S1H, (1 - p_S1D) * prog, p_S1D
P[1, 1] = (1 - p_S1D) * (1 - p_S1H - prog)
P[2, 2], P[2, 3], P[3, 3] = 1 - p_S2D, p_S2D, 1.0
add = {"Standard of care": 0.0, "Intervention A": p["c_trtA"],
"Intervention B": p["c_trtB"], "Intervention AB": p["c_trtA"] + p["c_trtB"]}[intervention]
cost = np.array([p["c_H"], p["c_S1"] + add, p["c_S2"] + add, 0.0])
u_s1 = p["u_trtA"] if intervention in ("Intervention A", "Intervention AB") else p["u_S1"]
return CohortSpec(P, cost, np.array([p["u_H"], u_s1, p["u_S2"], 0.0]))
engine = MarkovModel(states=STATES, interventions=INTERVENTIONS, transitions_and_rewards=model,
n_cycles=75, initial_state="H", cycle_correction="simpson")
params = ParameterSet({
"r_HD": Gamma(20, 1/10000), "r_HS1": Gamma(30, 1/200), "r_S1H": Gamma(60, 1/120),
"r_S1S2": Gamma(84, 1/800), "hr_S1": LogNormal(np.log(3), 0.01),
"hr_S2": LogNormal(np.log(10), 0.02), "hr_S1S2_trtB": LogNormal(np.log(0.6), 0.02),
"c_H": Gamma(100, 20.0), "c_S1": Gamma(177.8, 22.5), "c_S2": Gamma(225, 66.7),
"c_trtA": Gamma(73.5, 163.3), "c_trtB": Gamma(86.2, 150.8),
"u_H": Beta(200, 3), "u_S1": Beta(130, 45), "u_S2": Beta(230, 230), "u_trtA": Beta(300, 15),
})
draws = params.sample(10_000, seed=SeedManager(20260705).generator())
outcomes = run_psa(engine, draws).outcomesValue of information
Run the value-of-information workflow: expected value of perfect, partial perfect, and sample information.
Keywords
value of information, VOI, EVPI, EVPPI, EVSI, expected value of perfect information, expected value of partial perfect information, expected value of sample information, health economics, Python
This tutorial shows how to run the full value-of-information workflow on the Sick-Sicker cohort state-transition model from the cohort state-transition tutorial. It computes the expected value of perfect information (EVPI), of partial perfect information (EVPPI), and of sample information (EVSI). The full script is examples/voi_tutorial.py.
Specifying the decision model
The model is the four-state Sick-Sicker model from the cohort state-transition tutorial, evaluated over 75 annual cycles. It compares standard of care, Intervention A (raises the Sick-state utility), Intervention B (slows progression from Sick to Sicker), and the combination AB. That tutorial explains the model and the parameter distributions. The code below rebuilds both so this page runs on its own.
Value-of-information estimates need more draws than a mean cost-effectiveness estimate does, so this run uses 10,000 iterations against the cohort tutorial’s 1,000.
Calculating the expected value of perfect information (EVPI)
EVPI is the expected gain from making the decision with perfect information rather than current information. It is an upper bound on the value of any research.
from heormodel.voi import evpi
evpi(outcomes, WTP)2802.51210525009Ranking parameters with EVPPI
The expected value of partial perfect information (EVPPI) is the value of resolving one parameter, or a group, while the rest stay uncertain. It replaces the inner maximization of EVPI with an expectation taken over only the parameter of interest \(\phi\),
\[\text{EVPPI}_\phi = \mathbb{E}_\phi\!\left[\max_d \mathbb{E}_{\theta \mid \phi}\, \text{NMB}_d\right] - \max_d \mathbb{E}_\theta\, \text{NMB}_d,\]
so it identifies which parameters drive the decision uncertainty. evppi_ranking sweeps each of the model’s 16 parameters and sorts them into a research priority list.
from heormodel.voi import evppi, evppi_ranking
ranking = evppi_ranking(outcomes, draws, WTP)
ranking.round(1)u_S1 823.4
c_trtA 366.8
c_trtB 66.0
u_H 2.2
hr_S2 1.5
r_S1S2 1.4
hr_S1S2_trtB 0.0
hr_S1 0.0
c_H 0.0
c_S2 0.0
u_trtA 0.0
r_S1H 0.0
r_HS1 0.0
r_HD 0.0
c_S1 0.0
u_S2 0.0
Name: evppi, dtype: float64u_S1, the Sick-state utility under standard of care, dominates the ranking. Intervention A’s entire benefit is raising that utility to u_trtA, so whether the true u_S1 is already close to u_trtA decides whether Intervention A is worth its cost. The two treatment costs c_trtA and c_trtB follow, and everything else contributes almost nothing at this threshold. The default spline metamodel handles one parameter at a time well. For groups of many parameters evaluated jointly, use method="gp" (a Gaussian process).
Expected value of sample information (EVSI)
EVSI is the value of a study of finite sample size. It is the practical counterpart to the perfect-information measures above, which assume a parameter is resolved exactly. Since u_S1 drives the ranking, we consider a preference-based utility survey of Sick patients. Each patient reports a health-utility score with the sampling standard deviation typical of an instrument like the EQ-5D-5L in a chronic-disease population. simulate_summaries draws one survey summary per parameter draw from its true u_S1, and evsi_regression computes the value from the fitted conditional means.
from heormodel.voi import evsi_regression, simulate_summaries
PATIENT_SD = 0.20 # per-patient sd of a preference-based utility instrument
n_survey = 400
tau = PATIENT_SD / np.sqrt(n_survey)
summaries = simulate_summaries(
draws, lambda row, r: {"ubar": row["u_S1"] + r.normal(0.0, tau)}, seed=n_survey)
evsi_regression(outcomes, summaries, WTP)643.7317931556609Reported against the willingness-to-pay threshold, EVSI is a family of curves, one per sample size. Each EVSI curve stays below the u_S1 EVPPI, which in turn tracks EVPI closely. The survey targets the parameter EVPPI ranked highest, so resolving it captures most of the value of resolving everything. All the curves are small below about 90,000 per quality-adjusted life-year. In that range standard of care is optimal by a wide enough margin that the decision does not turn on u_S1. The curves peak sharply near 120,000, where the choice between interventions is closest, and fall off above that. Larger surveys shrink tau, which lifts each curve toward the EVPPI ceiling.
import matplotlib.pyplot as plt
from matplotlib import cm, colors
grid = np.linspace(0, 200_000, 21)
sizes = [25, 50, 100, 200, 400, 800, 1_600, 3_200, 6_400]
rng = np.random.default_rng(0)
fan = {}
for n in sizes: # the same survey, evaluated across the threshold grid
ubar = draws["u_S1"] + rng.normal(0.0, PATIENT_SD / np.sqrt(n), len(draws))
summ = pd.DataFrame({"ubar": ubar}, index=draws.index)
fan[n] = [evsi_regression(outcomes, summ, k) for k in grid]
ramp, scale = plt.get_cmap("Purples"), colors.LogNorm(min(sizes), max(sizes))
fig, ax = plt.subplots()
for n in sizes:
ax.plot(grid, fan[n], color=ramp(0.3 + 0.6 * scale(n)), lw=1.4)
ax.plot(grid, evpi(outcomes, grid).to_numpy(), color="#333333", lw=2, label="EVPI")
ax.plot(grid, [evppi(outcomes, draws, "u_S1", k) for k in grid],
color="#eb6834", lw=2, label="EVPPI (u_S1)")
ax.set(xlabel="willingness to pay", ylabel="value per person")
ax.legend(loc="upper left")
sm = cm.ScalarMappable(norm=scale, cmap=ramp); sm.set_array([])
fig.colorbar(sm, ax=ax, label="survey size (patients)");
Expected net benefit of sampling (ENBS)
How many patients should the survey enroll? The expected net benefit of sampling (ENBS) answers that. It multiplies the per-person EVSI by the discounted population the decision would inform, then subtracts the survey’s cost. Suppose 2,000 patients enter the Sick state each year for 10 years, discounted at 3% a year, and each would benefit from a decision informed by the resolved utility. The survey costs a fixed 100,000 to set up plus 300 per participant. A survey is worth running only if its ENBS is at least zero, and the size that maximizes ENBS is the one to fund.
i_wtp = int(np.where(grid == WTP)[0][0])
evsi = pd.Series({n: fan[n][i_wtp] for n in sizes}) # EVSI at the reference threshold
years = np.arange(10)
beneficiaries = 2_000 * (1.03 ** -years).sum() # discounted future patients
cost = 100_000 + 300 * np.array(sizes)
enbs = pd.Series(beneficiaries * evsi.to_numpy() - cost, index=pd.Index(sizes, name="patients"))
best = enbs.idxmax()
print(f"optimal size: {best} patients, ENBS {enbs[best]:,.0f}")
enbs.round(0)optimal size: 1600 patients, ENBS 13,399,952patients
25 2062662.0
50 5072873.0
100 7661816.0
200 10393951.0
400 12101950.0
800 13340125.0
1600 13399952.0
3200 13046084.0
6400 12383773.0
dtype: float64ENBS rises with survey size while the added precision is worth more than its cost. It falls once EVSI approaches the EVPPI ceiling, because each extra participant still costs 300 but resolves little further uncertainty.
ax = enbs.plot(marker="o", xlabel="survey size (patients)", ylabel="ENBS")
ax.axhline(0, color="0.6", lw=0.8);
ax.plot(best, enbs[best], "r*", ms=14);
Keywords: value of information, VOI, EVPI, EVPPI, EVSI, expected value of perfect information, expected value of partial perfect information, expected value of sample information, health economics, Python