Beta Distribution for Proportions and Utilities

In plain language

The beta distribution is a bell-shaped probability curve that always stays between 0 and 1, making it the standard choice for representing uncertainty about probabilities, adherence rates, and quality-of-life scores in health-economics models. To set one up, you take a published mean and standard error for any rate or proportion, plug them into two simple formulas (the method of moments), and get two shape numbers that fully define the distribution. In probabilistic sensitivity analysis — where a model is run thousands of times with slightly different assumed values — drawing from a beta distribution guarantees the model never receives an impossible input like a probability of 1.08 or a quality-of-life weight of −0.30, errors that a normal distribution would routinely produce for parameters near the 0 or 1 boundaries.

What is the beta distribution?

The beta distribution, written Beta(α, β) with shape parameters α > 0 and β > 0, is a continuous probability distribution whose support is strictly the open interval (0, 1). Its probability density is proportional to x^(α−1) × (1−x)^(β−1), and the two parameters control shape completely: when α = β = 1 the distribution is flat (uniform); when both exceed 1 it is unimodal and bell-shaped; when both are below 1 it is U-shaped. For health-economic modelling, the bell-shaped regime (α > 1, β > 1) is the usual case, representing a parameter that is most plausibly near an interior value rather than at the extremes. The mean is α/(α+β) and the variance is αβ/((α+β)² × (α+β+1)), which implies that the standard deviation of a beta-distributed quantity is fully determined by the mean and the total count α+β — a useful property when fitting from summary statistics.

Why [0,1] support matters in RWE

Many parameters that enter health-economic models are bounded between 0 and 1: event probabilities (response rate, transition probability, mortality risk), adherence proportions such as the Proportion of Days Covered (PDC), QALY utility weights, and quality-of-life index scores from instruments such as the EQ-5D-3L. The central reason to use a beta distribution for these quantities — rather than a normal distribution — is that the normal distribution has infinite support on both sides. In probabilistic sensitivity analysis (PSA), where thousands of random draws are taken from parameter distributions, a normal distribution centred near 0 or 1 will routinely produce draws below 0 or above 1 for parameters near the boundaries. A utility of −0.12 or a transition probability of 1.08 is not a legitimate model input and will either crash the model, trigger silent errors in transition matrices, or produce inadmissible cost-effectiveness results. The beta distribution eliminates this risk by construction: every draw lies in (0, 1). Near the centre of the parameter range, normal approximations often produce similar results; the difference is most consequential for probabilities or utilities in the range [0.05, 0.15] or [0.85, 0.95], where the normal tail already overlaps the inadmissible region at realistic standard errors.

Method of moments: fitting Beta(α, β) from a mean and standard error

In practice, an analyst seldom has access to patient-level data to fit a beta distribution by maximum likelihood. The typical input is a published point estimate and its standard error, extracted from a clinical trial, registry analysis, or systematic review. The method of moments provides a closed-form solution that requires no optimisation:

Step 1 — convert SE to variance: var = SE². Step 2 — compute total shape: α + β = mean × (1 − mean) / var − 1. Step 3 — compute α: α = mean × (α + β). Step 4 — compute β: β = (α + β) − α = (1 − mean) × (α + β). Step 5 — verify: mean = α/(α + β) and var = αβ/((α+β)² × (α+β+1)).

This procedure is unambiguous whenever var < mean × (1 − mean), which must hold for any proper beta distribution. If the reported variance exceeds this bound — for example, because a small study reported a wide confidence interval — the method of moments yields a non-positive shape parameter, signalling that the beta family cannot represent these statistics and that the analyst should either investigate the source of the variance estimate or consider a different distributional form.

The PSA workhorse: beta for probabilities and utilities

The ISPOR-SMDM Modelling Good Research Practices Task Force (Briggs et al., 2012) codified the distributional assignment rule now standard in submissions to NICE, CADTH, ICER, and other HTA bodies: beta distributions for parameters bounded in [0, 1] (probabilities, utilities, adherence rates); gamma or log-normal for strictly positive parameters (costs, length of stay); log-normal for relative risks and hazard ratios; normal only for parameters with genuine support across the real line (regression coefficients on an unconstrained scale). Violating this mapping — for example, assigning a normal distribution to a utility near 0.90 — produces PSA results where a nontrivial fraction of draws exceed 1 or fall below 0, and the expected net benefit estimator will be biased. Every PSA parameter table in an HTA submission should declare the distributional family for each parameter and justify it against this schema. In practice, the two-parameter set that covers the vast majority of HTA parameter tables is beta (probabilities, utilities) and gamma (costs, resource use) — the two distributions appear together across the published ISPOR good-modelling-practice literature.

Interpreting the output

Consider a fitted Beta(12, 3) for the utility parameter of a treatment responder, derived from an EQ-5D-3L registry analysis reporting mean utility 0.80 with SE 0.10.

Formal interpretation. The mean of Beta(12, 3) is 12/(12+3) = 12/15 = 0.80, reproducing the observed mean exactly by construction of the method-of-moments fit. The distribution encodes parameter uncertainty — specifically, the second-order uncertainty about the true population mean utility — not patient-to-patient heterogeneity in utility scores. In a PSA second-order Monte Carlo loop, each of the 5,000 iterations draws a single candidate value u_i from Beta(12, 3); this u_i is treated as the true utility for that iteration and all patients in that model run are assigned that utility. The 95% credible interval is read from the distribution's quantiles: qbeta(0.025, 12, 3) ≈ 0.55 and qbeta(0.975, 12, 3) ≈ 0.96. Interpreted through a conjugate Bayesian lens, Beta(12, 3) corresponds to having observed 12 pseudo-successes and 3 pseudo-failures in a prior dataset of 15 observations — the shape parameters play the role of posterior sufficient statistics when the beta is used as a prior for a binomial proportion.

Practical interpretation. In plain language: we believe the average utility of a treatment responder is around 0.80, but given our data we cannot rule out values as low as roughly 0.55 or as high as roughly 0.96. Each PSA iteration draws one candidate truth from this range. The width of the draw distribution directly contributes to the width of the cost-effectiveness acceptability curve and to the expected value of perfect information.

The critical confusion: parameter uncertainty vs patient heterogeneity. This is the most common PSA error in submitted HTA models, and it is worth stating explicitly. The beta distribution in PSA encodes uncertainty about the mean parameter value for the population, not individual patient variation around that mean. In a cohort Markov model, the correct PSA procedure is to draw u_i once from Beta(12, 3) per PSA iteration, then apply u_i uniformly to all patients in that model run. Drawing a different beta-distributed value for each simulated patient in a cohort model is a conceptual error that inflates PSA variance and will draw criticism from NICE Technical Support Documents and CADTH reviewers. In a microsimulation or discrete-event simulation that explicitly models individual patients, patient-level heterogeneity is captured by a separate distributional component over patients; the PSA outer loop remains a loop over uncertainty about the mean parameter. These two sources of variation — parameter uncertainty and patient heterogeneity — are distinct and must not be conflated.

Bayesian conjugacy: beta prior + binomial likelihood → beta posterior

The beta distribution is the conjugate prior for the binomial likelihood. If the prior on a probability p is Beta(α₀, β₀) and n independent Bernoulli trials produce k successes, the posterior is exactly Beta(α₀ + k, β₀ + n − k) — no numerical integration required. This conjugacy is valuable in rare-event settings. Suppose a disease has an annual event probability of approximately 0.03 from a published meta-analysis, and a new registry study observes 2 events in 50 person-years. A naive maximum-likelihood estimate of 2/50 = 0.04 is unstable at this sample size. Using the published estimate to set a Beta(3, 97) prior (encoding a prior mean of 3/100 = 0.03 with 100 pseudo-observations), the posterior after observing 2 events in 50 trials is Beta(5, 145), with posterior mean 5/150 ≈ 0.033 — a reasonable shrinkage toward the prior that respects existing evidence. This framework also allows the analyst to explicitly represent the strength of prior belief and to show HTA reviewers how sensitive results are to the prior specification.

Beta regression for bounded continuous outcomes

When the outcome variable itself is bounded in (0, 1) — rather than being a model parameter — beta regression (Ferrari & Cribari-Neto, 2004) is the appropriate regression framework. Beta regression models E[Y | X] via a link function (typically logit) on a linear predictor and models the conditional distribution of Y as beta, respecting the [0, 1] boundary without transformation. Applications in RWE and HEOR include: regressing PDC on patient characteristics, drug class, and plan design while preserving the adherence proportion scale; regressing EQ-5D index scores on disease severity and comorbidities without the boundary artifacts that arise from OLS near 0 and 1; and analysing health-plan quality metrics (HEDIS rates) as continuous proportions. Beta regression requires outcomes strictly in (0, 1); for outcomes with a mass at 0 or 1 (for example, complete non-initiation or perfect adherence), a zero-one-inflated beta (ZOIB) model adds mixture components at the boundaries. An alternative that is often pragmatically adequate when the fraction at the boundary is small is to recode exact 0 values to 0.001 and exact 1 values to 0.999, with this choice documented in the analysis plan. Beta regression is available in the betareg package in R and in the Smithson-Verkuilen (2006) formulation implemented in Stata and R; a statsmodels BetaModel is available in Python.

Negative utilities: states worse than death

Standard utility weights from the EQ-5D can take negative values for health states judged worse than death; the UK EQ-5D-3L value set includes utilities as low as approximately −0.594. A raw negative utility breaks the [0, 1] support of the beta distribution. The standard rescaling approach is to transform: u_rescaled = (u − u_min) / (u_max − u_min), where u_min is the minimum plausible utility (e.g., −0.594 for the UK EQ-5D-3L value set) and u_max = 1. The beta distribution is then fitted to u_rescaled, PSA draws are taken on the rescaled scale, and each draw is back-transformed before entering the model. Analysts should document the rescaling convention explicitly in the model technical report.

Pros, cons, and trade-offs

Pros of the beta distribution: - Strictly bounded to (0, 1), eliminating inadmissible PSA draws by construction. - Two-parameter family covering a wide range of shapes (symmetric, left-skewed, right-skewed, near-uniform) depending on α and β, fitting any mean and variance satisfying the feasibility constraint. - Conjugate with the binomial likelihood, enabling Bayesian updating in closed form with no numerical integration. - Closed-form moments and quantile functions available in all statistical packages; PSA draws require a single call (rbeta in R, scipy.stats.beta.rvs in Python, RAND('BETA') in SAS). - ISPOR-SMDM endorsed distributional assignment for probabilities and utilities in HTA submissions to NICE, CADTH, and other bodies. - Beta regression (Ferrari & Cribari-Neto 2004) extends the family to regression of bounded continuous outcomes, with full R, Stata, and Python support.

Cons and limitations: - Cannot accommodate exact 0 or 1 values without inflation extensions; boundary observations require a ZOIB model or recoding, each of which adds a modelling decision to document. - Cannot directly represent negative utilities; a rescaling step is required before fitting and after drawing, adding complexity and a parameter choice (u_min) to justify. - Method-of-moments fit requires var < mean × (1 − mean); published SEs that are implausibly large relative to the mean will yield non-positive shape parameters. - In large individual-level datasets, fitting by method of moments discards distributional shape information beyond mean and variance; maximum-likelihood or Bayesian estimation from patient-level data is preferred when those data are available. - Beta regression can exhibit convergence issues with boundary-heavy data and requires diagnostic inspection (randomised quantile residuals, pseudo-R²).

When to use

• PSA parameter distributions for probabilities, utilities, and proportions: whenever a
• Bayesian prior on a proportion or event probability: particularly valuable in rare-
• Beta regression for bounded continuous outcomes such as PDC, EQ-5D index scores, and
• Overdispersed count-denominator proportions: the beta-binomial distribution (a beta

When NOT to use

• Exact 0 or 1 values present without inflation handling: standard beta regression will
• Utilities below 0 that have not been rescaled: the beta distribution cannot represent
• Modelling binary outcomes themselves: if each observation is a 0/1 event indicator,
• Patient-level heterogeneity in a cohort Markov model: drawing a different beta-
• Outcomes with genuinely unbounded or positive-only support: log-normal or gamma

Worked example

Scenario

A health economist is building a Markov cost-effectiveness model for a new treatment in relapsing multiple sclerosis. One key input is the utility weight for a patient who responds to treatment, taken from a published EQ-5D-3L registry analysis that reported a mean utility of 0.80 with a standard error of 0.10 (from 120 patients). The analyst must fit a beta distribution to this parameter for PSA so that all 5,000 PSA iterations draw a plausible value between 0 and 1. Only the published summary statistics are available — no patient-level data — so the method of moments is used.

Dataset

Summary statistics from the published registry analysis. Only the mean and SE are available; these two numbers are the complete input to the method-of-moments fitting procedure.

Steps

• Step 1 — Convert SE to variance: var = 0.10*0.10 = 0.01.

• Step 2 — Compute the total shape parameter alpha+beta. The method-of-moments formula is mean(1-mean)/var - 1 = 0.80.2/0.01 - 1 = 0.16/0.01 - 1 = 16 - 1 = 15. So alpha+beta = 15.

• Step 3 — Compute alpha from the mean: alpha = mean(alpha+beta) = 0.815 = 12.

• Step 4 — Compute beta as the remainder: beta = (alpha+beta) - alpha = 15 - 12 = 3.

• Step 5 — Verify the mean: alpha/(alpha+beta) = 12/(12+3) = 12/15 = 0.80. Correct — the fitted distribution reproduces the published mean exactly.

• Step 6 — Verify the variance: alphabeta/((alpha+beta)(alpha+beta)(alpha+beta+1)) = 123/(151516) = 36/3600 = 0.01. This matches the target var = 0.10*0.10 = 0.01 — the fit is exact by construction of the method of moments.

• Step 7 — In PSA, call rbeta(5000, 12, 3) in R or scipy.stats.beta(12, 3).rvs(5000) in Python to draw 5,000 plausible utility values. All draws fall in (0,1). The 2.5th and 97.5th percentiles of Beta(12,3) are approximately 0.55 and 0.96, spanning the plausible range for this utility parameter in the model.

Result

Fitted distribution: Beta(alpha=12, beta=3). Theoretical mean = 12/(12+3) = 12/15 = 0.80. Variance = 123/(151516) = 36/3600 = 0.01, matching target SE squared = 0.100.10 = 0.01. Every PSA draw lies in (0,1), eliminating inadmissible model inputs. The 95% interval of roughly [0.55, 0.96] conveys to HTA reviewers how much uncertainty about this single utility parameter feeds into the overall cost-effectiveness uncertainty.

Runnable example

import numpy as np
from scipy import stats

# ── Method-of-moments fit: Beta(alpha, beta) from mean and SE ──────────────────
mean_util = 0.80
se_util   = 0.10
var_util  = se_util ** 2                        # 0.10 ** 2 = 0.01

# Feasibility check: var must be < mean*(1-mean)
max_var = mean_util * (1 - mean_util)           # 0.80 * 0.20 = 0.16
if var_util >= max_var:
    raise ValueError(
        f"var={var_util} >= mean*(1-mean)={max_var}; "
        "Beta distribution cannot represent these statistics."
    )

alpha_plus_beta = mean_util * (1 - mean_util) / var_util - 1  # 0.16/0.01 - 1 = 15
alpha = mean_util * alpha_plus_beta                             # 0.8 * 15 = 12
beta  = (1 - mean_util) * alpha_plus_beta                      # 0.2 * 15 = 3

print(f"Fitted Beta({alpha:.4g}, {beta:.4g})")
print(f"  Theoretical mean:     {alpha / (alpha + beta):.4f}  (target: {mean_util})")
print(f"  Theoretical variance: "
      f"{alpha * beta / ((alpha + beta)**2 * (alpha + beta + 1)):.6f}"
      f"  (target: {var_util})")

# ── Distributional summary ──────────────────────────────────────────────────────
dist = stats.beta(alpha, beta)
print(f"  Mode:  {(alpha - 1) / (alpha + beta - 2):.4f}")
print(f"  95% CI: [{dist.ppf(0.025):.4f}, {dist.ppf(0.975):.4f}]")

# ── PSA Monte Carlo: 5,000 draws, all guaranteed in (0, 1) ───────────────────
np.random.seed(42)
psa_draws = dist.rvs(5_000)

assert psa_draws.min() > 0 and psa_draws.max() < 1, "All PSA draws must be in (0,1)"
print(f"\nPSA draws (n=5,000): mean={psa_draws.mean():.4f}, "
      f"SD={psa_draws.std():.4f}, "
      f"min={psa_draws.min():.4f}, max={psa_draws.max():.4f}")

# ── Bayesian conjugate update: Beta prior + binomial data -> Beta posterior ────
# Prior: Beta(3, 97) — prior mean = 3/100 = 0.03 (from published meta-analysis)
# Observed: 2 events in 50 trials (new registry)
prior_a, prior_b     = 3, 97
obs_k, obs_n         = 2, 50
post_a = prior_a + obs_k                         # 3 + 2 = 5
post_b = prior_b + (obs_n - obs_k)               # 97 + 48 = 145
posterior = stats.beta(post_a, post_b)
print(f"\nBayesian conjugate update:")
print(f"  Prior:     Beta({prior_a}, {prior_b}), mean = {prior_a/(prior_a+prior_b):.4f}")
print(f"  Observed:  {obs_k} events in {obs_n} trials (MLE = {obs_k/obs_n:.4f})")
print(f"  Posterior: Beta({post_a}, {post_b}), "
      f"mean = {post_a/(post_a+post_b):.4f}")

# ── Beta regression for PDC or EQ-5D outcomes (statsmodels) ──────────────────
# from statsmodels.othermod.betareg import BetaModel
# mod = BetaModel.from_formula("pdc ~ age + cci", data=df)
# res = mod.fit()
# Outcomes must be strictly in (0,1). Recode: pdc = pdc.clip(0.001, 0.999)
# Coefficients are on the logit scale; back-transform with scipy.special.expit().

# ── Method-of-moments fit ─────────────────────────────────────────────────────
mean_util <- 0.80
se_util   <- 0.10
var_util  <- se_util^2              # 0.01

# Feasibility check
stopifnot(var_util < mean_util * (1 - mean_util))

alpha_plus_beta <- mean_util * (1 - mean_util) / var_util - 1   # 15
alpha <- mean_util * alpha_plus_beta                              # 12
beta  <- (1 - mean_util) * alpha_plus_beta                       # 3

cat(sprintf("Fitted Beta(%.4g, %.4g)\n", alpha, beta))
cat(sprintf("  Theoretical mean:     %.4f  (target %.2f)\n",
            alpha / (alpha + beta), mean_util))
cat(sprintf("  Theoretical variance: %.6f  (target %.4f)\n",
            alpha * beta / ((alpha + beta)^2 * (alpha + beta + 1)),
            var_util))

# ── Quantiles and mode ────────────────────────────────────────────────────────
cat(sprintf("  Mode:  %.4f\n", (alpha - 1) / (alpha + beta - 2)))
cat(sprintf("  95%% CI: [%.4f, %.4f]\n",
            qbeta(0.025, alpha, beta), qbeta(0.975, alpha, beta)))

# ── PSA Monte Carlo ───────────────────────────────────────────────────────────
set.seed(42)
psa_draws <- rbeta(5000, alpha, beta)

stopifnot(all(psa_draws > 0) && all(psa_draws < 1))
cat(sprintf("\nPSA draws (n=5000): mean=%.4f, SD=%.4f, min=%.4f, max=%.4f\n",
            mean(psa_draws), sd(psa_draws), min(psa_draws), max(psa_draws)))

# ── Bayesian conjugate update ─────────────────────────────────────────────────
prior_a <- 3;  prior_b <- 97       # Beta(3,97): prior mean = 0.03
obs_k   <- 2;  obs_n   <- 50       # 2 events in 50 trials
post_a  <- prior_a + obs_k         # 5
post_b  <- prior_b + obs_n - obs_k # 145
cat(sprintf(
  "\nBayesian update: Beta(%d,%d) + Bin(n=%d,k=%d) -> Beta(%d,%d), mean=%.4f\n",
  prior_a, prior_b, obs_n, obs_k, post_a, post_b,
  post_a / (post_a + post_b)
))

# ── Beta regression with betareg ──────────────────────────────────────────────
# install.packages("betareg")  # run once if not installed
# library(betareg)
#
# Simulate a dataset: PDC as outcome, age and CCI as predictors
set.seed(1)
n   <- 200
df  <- data.frame(
  pdc = rbeta(n, 8, 2),                 # mean PDC ≈ 0.80
  age = rnorm(n, mean = 55, sd = 12),
  cci = rpois(n, lambda = 2)
)
# Recode exact boundary values (betareg requires strictly open (0,1))
df$pdc <- pmin(pmax(df$pdc, 0.001), 0.999)
#
# Fit beta regression (logit link for mean, log link for precision phi):
# mod <- betareg(pdc ~ age + cci, data = df)
# summary(mod)
# Marginal effect on PDC scale: plogis(coef(mod)["age"])