Probabilistic Sensitivity Analysis (PSA) for Health-Economic Models

In plain language

Probabilistic sensitivity analysis (PSA) answers the question: given everything we are uncertain about in a cost-effectiveness model, how confident are we that a treatment is worth its price? Instead of plugging a single best-guess number into each model input, PSA treats each input as a range with a most-likely value and plausible spread, then re-runs the model thousands of times, each time drawing a fresh random set of inputs. The result is not one cost-effectiveness ratio but a cloud of thousands, from which you can read off the probability that the treatment is cost-effective at any given price threshold. One honest caveat: PSA captures uncertainty about the numbers we feed the model, but it cannot capture uncertainty about whether the model structure itself is correct.

Probabilistic sensitivity analysis (PSA)

treats every estimated input of a decision-analytic model (transition probabilities, hazard ratios, utilities, resource use, unit costs) as a random variable with a distribution that reflects its sampling uncertainty, then draws a large Monte Carlo sample (commonly 5,000-10,000 iterations), runs the full model once per draw, and summarizes the resulting joint distribution of incremental costs and incremental effects (QALYs). The outputs are the cost-effectiveness (CE) plane scatter, the cost-effectiveness acceptability curve (CEAC), the expected net benefit, and — increasingly required by HTA bodies — the expected value of perfect information (EVPI). In real-world evidence (RWE), the inputs are not assumed; they are estimated from claims, EHR, registry, or linked data, so PSA is the bridge between the uncertainty in your RWE estimates and the uncertainty in the reimbursement decision.

Core conceptual distinction

Three distinctions are doing the work and they are routinely conflated. (1) PSA vs deterministic (one-way/multi-way) sensitivity analysis: deterministic analysis varies one or a few parameters across plausible bounds to find drivers; PSA varies all uncertain parameters simultaneously according to their distributions and so quantifies the probability that a decision is correct given everything we do not know. Deterministic SA answers "what matters?"; PSA answers "how sure are we?" — they are complements, not substitutes. (2) Parameter uncertainty vs structural uncertainty: PSA addresses only second-order (parameter) uncertainty. Choices of model structure — time horizon, cycle length, half-cycle correction, the competing-risk specification, the survival-extrapolation family — are not captured by PSA and must be handled by scenario analysis. Presenting a tight CEAC while ignoring structural uncertainty is a known way to overstate confidence and is exactly what reviewers (and NICE/CADTH critique) flag. (3) Sampling uncertainty vs variability: distributions in PSA represent uncertainty about a mean parameter (the standard error of an estimate), not patient-to-patient heterogeneity; heterogeneity belongs in subgroup analysis or individual-level (microsimulation/DES) models, not in the PSA distributions of cohort-model means.

Pros, cons, and trade-offs

- vs deterministic sensitivity analysis (tornado diagrams): PSA yields a coherent probabilistic statement of decision uncertainty (P[cost-effective at threshold k]) and supports value-of-information analysis; deterministic SA cannot. Cost: PSA requires specifying a distribution and, critically, a correlation structure for every parameter, is computationally heavier, and can give a false sense of precision if distributions are guessed. Prefer PSA as the headline uncertainty analysis for any reimbursement submission, with deterministic SA retained to identify drivers and to communicate. - vs reporting a point ICER with a confidence interval: A single ICER ratio is statistically ill-behaved (the ratio of two uncertain quantities; CIs can straddle quadrants of the CE plane and become uninterpretable). PSA on the net-benefit scale sidesteps the ratio problem and produces a well-defined CEAC. Prefer PSA / net-benefit whenever uncertainty must be communicated. - vs bootstrapping a within-trial / within-cohort cost-effectiveness analysis: When costs and effects are observed at the patient level in one RWE dataset, nonparametric bootstrap of the patient-level data propagates the real (correlated, skewed) joint distribution directly and is often preferable. PSA on a decision model is required instead when the model synthesizes parameters from multiple sources, extrapolates beyond observed follow-up, or links exposure effects to long-term costs/QALYs. Prefer bootstrap for a self-contained patient-level CEA; prefer model-based PSA for synthesis and extrapolation. - vs the E-value and other causal sensitivity analyses: these quantify robustness of a causal estimate to unmeasured confounding; PSA quantifies robustness of a decision to parameter sampling uncertainty. They are not interchangeable — a confounded hazard ratio fed into PSA produces a precisely characterized but biased decision.

When to use

Any cost-utility or cost-effectiveness model intended for HTA submission (NICE, CADTH, ICER, IQWiG, PBAC) where decision uncertainty and value-of-information must be reported; whenever the model synthesizes RWE parameters whose standard errors are known; whenever the deterministic ICER sits near the decision threshold and the reviewer needs the probability of cost-effectiveness; and as the substrate for EVPI/EVPPI to prioritize future evidence generation.

When NOT to use — and when it is actively misleading or dangerous

- As a substitute for handling structural uncertainty. A narrow CEAC built on one fixed structure (e.g., a particular extrapolation curve) is misleadingly confident. If extrapolation or model structure drives the result, PSA understates total uncertainty — actively dangerous for a reimbursement decision. - When parameters are sampled independently but came from one regression. Treating jointly estimated coefficients (e.g., a multivariable cost model, a parametric survival fit, a multinomial transition model) as independent miscalibrates the propagated uncertainty — typically inflating it and distorting the CEAC. Use the full variance-covariance matrix (multivariate normal on the natural/log scale; Dirichlet for transition rows). - When the point estimates are biased. PSA propagates uncertainty, not bias. Garbage-in confounded RWE effects yield a confidently wrong decision. Resolve confounding/immortal-time/selection issues in the estimation step first; PSA cannot rescue them. - When distributions are invented without an evidence basis. Assuming convenient (often too-narrow or symmetric) distributions for skewed costs or bounded utilities will mis-state tail behavior and the decision probability.

Data-source operational depth

- Claims (FFS or commercial): Costs are right-skewed with heavy tails (ICU stays, dialysis, biologics, cost outliers); a Gamma or log-normal fitted to the raw arm mean understates the upper tail and the cost variance that feeds PSA. Fit on appropriately handled data (Gamma/log-normal with explicit outlier handling) and run a separate deterministic scenario on outlier trimming. Event rates that become transition probabilities should be derived from FFS-observable person-time only: Medicare Advantage enrollees lack adjudicated FFS claims, so including MA-only person-time biases both rates and per-member costs — restrict to Parts A/B/D FFS and document whether a 5% sample or 100% denominator was used. In elderly claims cohorts, differential competing risk of death by exposure distorts cause-specific transition probabilities; derive them with a competing-risk model and sample the sub-distribution, not naive Kaplan-Meier complements. - EHR: Utilities/PROs and disease-severity states are richer than in claims, but visit-driven capture makes state occupancy and resource use differentially missing for patients who leave the system; loss to follow-up is potentially informative and should be reflected as additional (scenario) uncertainty rather than ignored. Unit costs are usually absent from EHR and must be imported (and their import uncertainty represented). - Registry: Strong for adjudicated clinical states and disease-severity transitions (good source for transition probabilities and their SEs) but weak for complete cost capture; link to claims for costs and to a death index to pin the absorbing state, and propagate the linkage-completeness uncertainty. - Linked claims–EHR–vital records: The ideal substrate — EHR severity for utilities, claims for costs, vital records for the death state — but only the linkable subset is analyzable, and order/fill/service-date discrepancies must be reconciled before deriving cycle-level rates. The selection induced by linkage is itself a structural-uncertainty scenario.

Worked example (RWE → PSA → decision)

Decision: a novel oral anticoagulant (NOAC) vs warfarin for stroke prevention in non-valvular atrial fibrillation, costed to a payer over a lifetime horizon in a 3-month-cycle Markov cost-utility model with states {AF-no-event, post-ischemic-stroke, post-major-bleed, dead}. RWE inputs and their PSA distributions: (1) Treatment effect — the adjusted hazard ratio for ischemic stroke (NOAC vs warfarin) comes from an active-comparator new-user analysis in linked claims (continuous Parts A/B/D enrollment, 365-day washout, index = first NOAC/warfarin `fill_date`, first-event coding, censoring at disenrollment/death/end-of-data); sample log(HR) ~ Normal(log Ĥ R, SE) and exponentiate, drawing it once per iteration and applying it to the baseline stroke transition so effect uncertainty propagates coherently. (2) Baseline transitions — from the warfarin arm's FFS person-time, computed with a competing-risk model for stroke vs bleed vs death; sample each transition row as Dirichlet(α = event counts) so the row sums to 1 and counts drive precision. (3) Utilities — state utilities from EHR-linked PRO data, sampled Beta(a,b) (bounded 0-1). (4) Costs — acute stroke, acute bleed, and per-cycle maintenance costs from claims `paid_amount` aggregated to the cycle, fitted Gamma to respect right skew, with an outlier-trimming scenario held aside. (5) Run 10,000 iterations, each drawing the full correlated parameter vector (use the survival model's variance-covariance matrix for any jointly estimated terms), evaluate the model, and store incremental cost ΔC and incremental QALYs ΔE. Outputs: CE-plane scatter; CEAC giving P(NOAC cost-effective) at WTP = $50k, $100k, $150k/QALY; expected incremental net monetary benefit; and population EVPI to judge whether the residual uncertainty justifies a confirmatory study. Report PSA alongside deterministic one-way SA (drivers) and structural scenarios (time horizon, extrapolation family, competing-risk handling) — the PSA alone does not cover the last.

Interpreting the output

The worked example reports that 780 of 1,000 PSA iterations produce an ICER below $100,000/QALY, giving a cost-effectiveness acceptability curve (CEAC) value of 78% at that threshold.

(1) Formal interpretation. The CEAC value of 78% at λ = $100,000/QALY means that in 78% of the joint parameter draws — each representing one plausible combination of all uncertain model inputs — the intervention is estimated to be cost-effective at that willingness-to-pay level. This quantifies decision uncertainty under the model: it is not the probability that the intervention works for a given patient, nor a statement about the fraction of patients who benefit. The 22% of draws above the threshold reflect the combined effect of uncertainty across all distributional parameters simultaneously — utility weights, transition probabilities, unit costs, and relative effects — propagated jointly through the model structure.

(2) Practical interpretation. A CEAC of 78% at $100,000/QALY means the decision to adopt the intervention carries meaningful but not extreme residual uncertainty. HTA bodies typically want to see both the CEAC and the expected value of perfect information (EVPI), which converts this uncertainty into the monetary value of resolving it — a key input to the value-of-research argument. One critical error to avoid: reporting "there is a 78% probability the drug works" from a PSA CEAC. The 78% is a property of the decision model's parameter uncertainty, not a patient-level efficacy probability. Report the CEAC alongside the base-case ICER and the tornado diagram from one-way deterministic sensitivity analysis; the three together provide a complete picture of uncertainty.

Worked example

Scenario

A health economist is building a simple cost-effectiveness model comparing a new drug (Drug A) versus standard care for a chronic condition. The model has three uncertain inputs: the annual cost of Drug A, the reduction in hospitalizations Drug A produces, and the quality-of-life benefit per hospitalization avoided. Instead of running the model once with single best-guess numbers, she runs 1,000 Monte Carlo iterations, drawing each input randomly from its distribution each time. She then tallies what fraction of those 1,000 runs show Drug A to be cost-effective at a willingness-to-pay threshold of $100,000 per QALY gained.

Dataset

A sample of five simulation draws from the 1,000-iteration PSA run. Each row is one complete re-run of the model with a different randomly sampled set of inputs. The final two columns show the resulting incremental cost and incremental QALYs for that draw.

Steps

• For each of the 1,000 draws, sample drug_a_annual_cost_usd from a Gamma distribution (right-skewed because costs cannot be negative and tend to have a long upper tail), sample hosp_reduction_rate from a Beta distribution (bounded between 0 and 1 because it is a probability), and sample qaly_per_hosp_avoided from a Beta distribution for the same reason.

• Run the cost-effectiveness model once using that draw's three sampled values to produce incremental_cost_usd (the extra cost of Drug A vs standard care) and incremental_qalys (the extra health benefit).

• Compute the implied ICER for that draw: incremental_cost_usd divided by incremental_qalys. For draw 1: $14,800 / 0.31 QALYs = $47,742 per QALY.

• Classify each draw as cost-effective if its ICER is below the $100,000 per QALY threshold (equivalently, if 100,000 x incremental_qalys minus incremental_cost_usd is positive).

• After all 1,000 draws, count the fraction classified as cost-effective to get the CEAC value at $100,000 per QALY.

Result

In this illustrative 5-draw sample, all five draws have an ICER below $100,000 per QALY (ICERs range from roughly $48k to $86k), so 5/5 = 100% are cost-effective in this small sample. Across all 1,000 full-simulation draws using the file's model parameters, 780 of 1,000 draws fall below the $100,000 threshold, giving a CEAC value of 78% at that threshold. This means the model estimates a 78% probability that Drug A is cost-effective at a willingness-to-pay of $100,000 per QALY.

Runnable example

import numpy as np

RNG = np.random.default_rng(20240101)
N_ITER, N_CYCLES = 10_000, 200          # 3-month cycles over a lifetime horizon
WTP = np.array([50_000, 100_000, 150_000])
STATES = ["AF", "post_stroke", "post_bleed", "dead"]
DISC_C, DISC_E = 0.03, 0.03             # annual discount rates (apply per cycle)

def run_psa(transition_counts, util_beta, cost_gamma, loghr_stroke):
    # --- one draw of every uncertain parameter, sampled jointly across N_ITER ---
    # Transition rows: Dirichlet(counts) keeps each row a valid probability vector; counts drive precision.
    P = {arm: np.stack([RNG.dirichlet(c, size=N_ITER) for c in rows], axis=1)   # (N_ITER, n_states, n_states)
         for arm, rows in transition_counts.items()}
    # Treatment effect drawn ONCE per iteration on the log scale, then exponentiated.
    hr = np.exp(RNG.normal(loghr_stroke[0], loghr_stroke[1], N_ITER))
    u = np.stack([RNG.beta(*util_beta[s], size=N_ITER) for s in STATES], axis=1)  # utilities, bounded 0..1
    cost = np.stack([RNG.gamma(cost_gamma[s][0], cost_gamma[s][1], N_ITER)        # right-skewed per-cycle cost
                     for s in STATES], axis=1)

    results = {}
    for arm in ("warfarin", "noac"):
        Parm = P[arm].copy()
        if arm == "noac":                      # apply RWE hazard ratio to the AF->stroke transition
            Parm[:, 0, 1] *= hr
            Parm[:, 0, 0] = 1.0 - Parm[:, 0, 1:].sum(axis=1)   # renormalize the originating row
        trace = np.zeros((N_ITER, len(STATES)))
        trace[:, 0] = 1.0                       # everyone starts in AF
        tot_c = np.zeros(N_ITER); tot_e = np.zeros(N_ITER)
        for t in range(N_CYCLES):
            trace = np.einsum("is,isj->ij", trace, Parm)        # advance one cycle
            dfac_c = 1 / (1 + DISC_C) ** (t / 4); dfac_e = 1 / (1 + DISC_E) ** (t / 4)
            tot_c += (trace * cost).sum(axis=1) * dfac_c
            tot_e += (trace * u).sum(axis=1) * 0.25 * dfac_e    # quarter-cycle QALYs
        results[arm] = (tot_c, tot_e)

    dC = results["noac"][0] - results["warfarin"][0]
    dE = results["noac"][1] - results["warfarin"][1]
    ceac = {k: float(np.mean(k * dE - dC > 0)) for k in WTP}    # P(NOAC cost-effective) by WTP
    inb = WTP[1] * dE - dC                                      # incremental NMB at $100k
    evpi = float(np.maximum(inb, 0).mean() - max(inb.mean(), 0))
    return {"dC": dC, "dE": dE, "ceac": ceac, "evpi_per_patient": evpi}

library(MCMCpack)   # rdirichlet
set.seed(20240101)
N_ITER <- 10000L; N_CYCLES <- 200L
WTP <- c(50000, 100000, 150000)
STATES <- c("AF", "post_stroke", "post_bleed", "dead")
DISC_C <- 0.03; DISC_E <- 0.03

run_psa <- function(transition_counts, util_beta, cost_gamma, loghr_stroke) {
  draw_P <- function(rows)                      # list of Dirichlet-sampled transition matrices, one per iteration
    lapply(seq_len(N_ITER), function(i)
      t(vapply(rows, function(cnt) rdirichlet(1, cnt), numeric(length(STATES)))))
  P <- lapply(transition_counts, draw_P)
  hr <- exp(rnorm(N_ITER, loghr_stroke[1], loghr_stroke[2]))         # effect drawn once per iteration
  u  <- sapply(STATES, function(s) rbeta(N_ITER, util_beta[[s]][1], util_beta[[s]][2]))
  cost <- sapply(STATES, function(s) rgamma(N_ITER, cost_gamma[[s]][1], scale = cost_gamma[[s]][2]))

  arm_totals <- function(arm) {
    tot_c <- numeric(N_ITER); tot_e <- numeric(N_ITER)
    for (i in seq_len(N_ITER)) {
      Pi <- P[[arm]][[i]]
      if (arm == "noac") {                       # apply RWE hazard ratio to AF->stroke, renormalize the row
        Pi[1, 2] <- Pi[1, 2] * hr[i]
        Pi[1, 1] <- 1 - sum(Pi[1, -1])
      }
      tr <- c(1, 0, 0, 0)
      for (t in seq_len(N_CYCLES)) {
        tr <- as.numeric(tr %*% Pi)
        dC <- 1 / (1 + DISC_C)^((t - 1) / 4); dE <- 1 / (1 + DISC_E)^((t - 1) / 4)
        tot_c[i] <- tot_c[i] + sum(tr * cost[i, ]) * dC
        tot_e[i] <- tot_e[i] + sum(tr * u[i, ]) * 0.25 * dE
      }
    }
    list(c = tot_c, e = tot_e)
  }
  w <- arm_totals("warfarin"); n <- arm_totals("noac")
  dC <- n$c - w$c; dE <- n$e - w$e
  ceac <- sapply(WTP, function(k) mean(k * dE - dC > 0))
  inb  <- WTP[2] * dE - dC
  evpi <- mean(pmax(inb, 0)) - max(mean(inb), 0)
  list(dC = dC, dE = dE, ceac = setNames(ceac, WTP), evpi_per_patient = evpi)
}