Weibull Distribution for Time-to-Event Data

In plain language

The Weibull distribution is a mathematical formula for describing how quickly patients in a study reach some outcome — like death or disease progression — over time. Unlike simpler models that assume risk is constant, the Weibull lets risk rise or fall depending on a "shape" parameter: a shape above 1 means risk accumulates over time (like progressive disease), exactly 1 means constant risk (like background mortality in healthy adults), and below 1 means risk declines (like early post-surgical complications). It is especially useful in health economics because it is the only standard survival model that can simultaneously report a hazard ratio for clinicians and a time ratio for economic models from the same fit, and because regulators and health technology assessment bodies require it as one of the standard candidate models for projecting survival beyond trial follow-up.

What the Weibull distribution is and why it matters in RWE

The Weibull distribution is a two-parameter family for strictly positive, continuous outcomes — most commonly elapsed time to a clinical event such as death, hospitalization, disease progression, or treatment discontinuation. It is parameterized by a shape parameter k (also written α or γ depending on software convention) and a scale parameter λ (also written σ, θ, or λ). Together they control the full shape of the hazard function over time. The survival function is S(t) = exp(−(t/λ)^k) and the hazard is h(t) = (k/λ)(t/λ)^(k−1), giving a power-law relationship between elapsed time and instantaneous event risk.

In RWE and health economics, the Weibull occupies a central position for three reasons. First, it nests the exponential distribution as a special case (k = 1, constant hazard), so it generalizes the simplest survival model without requiring a different software pipeline. Second, it is the only member of the standard six parametric distributions used in HTA extrapolation (exponential, Weibull, Gompertz, log-normal, log-logistic, generalized gamma) that simultaneously satisfies both the proportional hazards (PH) and accelerated failure time (AFT) model structures — a property that makes it the canonical bridge between the semi-parametric Cox world and the fully parametric extrapolation world. Third, its monotone hazard is mechanistically plausible for a wide range of clinical trajectories: early-failure burn-in (k < 1), steady-state (k = 1), and progressive wear-out (k > 1).

The survival function, hazard function, and parameterizations

Given shape k > 0 and scale λ > 0, the Weibull specifies: - Survival function: S(t) = exp(−(t/λ)^k) - Hazard function: h(t) = (k/λ)(t/λ)^(k−1) - Cumulative hazard: H(t) = (t/λ)^k - Mean (expectation): E[T] = λ · Γ(1 + 1/k), where Γ is the gamma function

An alternative parameterization common in AFT software writes the model as log T = μ + σ ε, where ε follows the standard extreme-value (Gumbel) distribution, and the shape is recovered as k = 1/σ. SAS PROC LIFEREG uses this log-linear AFT convention: the intercept is μ = log(λ) and the scale output is σ = 1/k. The R survreg function also uses this convention. The flexsurvreg function in R and scipy.stats.weibull_min in Python use the direct (k, λ) parameterization. Parameterization mismatch is one of the most common implementation errors — always verify which k you are receiving before reporting it.

Hazard shapes by shape parameter k

The shape parameter k is the most clinically consequential quantity the Weibull delivers:

k < 1 (decreasing hazard): Risk starts high and declines monotonically. This matches an early-failure or burn-in pattern — post-surgical complications, acute toxicity after chemotherapy, or frailty-driven short-term mortality in elderly cohorts. The hazard approaches infinity at t = 0 and declines toward zero as t increases.

k = 1 (constant hazard): Reduces to the exponential distribution. Appropriate for memoryless processes where the risk of the event in the next short interval does not depend on how long the patient has already survived — background mortality in healthy adults over short follow-up, or Poisson-process equipment failure.

k > 1 (increasing hazard): Risk accumulates over time, as in progressive diseases such as advancing heart failure, Parkinson's disease, or accumulating cancer burden. k = 2 gives a linear hazard (the Rayleigh distribution); larger k gives a faster acceleration. Matching the shape parameter to the clinical mechanism is the primary substantive check on a Weibull fit — an increasing-hazard Weibull applied to post-MI survival data (which typically shows an early hazard spike then decline) signals model misspecification regardless of AIC.

Dual citizenship: proportional hazards and accelerated failure time simultaneously

The Weibull's unique mathematical property is that it is the only continuous distribution satisfying both the proportional hazards assumption and the accelerated failure time assumption simultaneously. This dual citizenship makes it the canonical bridge between two major modeling frameworks.

In the PH frame, a Weibull model with treatment covariate x specifies h(t|x) = h_0(t) · exp(β x), where h_0(t) is a Weibull baseline hazard. The quantity exp(β) is a hazard ratio (HR), interpreted the same way as a Cox HR — the ratio of instantaneous event rates between groups, assumed constant over time.

In the AFT frame, the same model specifies log T = μ + γ x + σ ε. The quantity exp(γ) is a time ratio (TR) — the factor by which the entire time scale is stretched or compressed for the treated group. A time ratio of 1.50 means every quantile of the survival distribution (median, 75th percentile, 90th percentile) is 1.50 times as large under treatment. The conversion between the two representations for the Weibull is: HR = TR^(−k). This means the hazard ratio can be derived from an AFT fit, and vice versa, once the shape parameter k is known. This conversion is not available for log-normal or log-logistic distributions, which satisfy only the AFT structure.

Interpreting the output

Consider a Weibull AFT model comparing treated versus control for time to disease progression. The software reports an AFT treatment coefficient of 0.405 with shape k = 1.4.

The time ratio is exp(0.405) ≈ 1.50. The formal interpretation is: event times in the treated arm are scaled by a factor of 1.50 relative to control, conditional on covariates. Every quantile of the progression-free survival distribution — the median, the 75th percentile, the 90th percentile — is 1.50 times larger under treatment. Because Weibull also satisfies PH, the equivalent hazard ratio is HR = TR^(−k). For these example values, TR = 1.50 and k = 1.4, the HR is approximately 0.59, meaning the instantaneous rate of progression in the treated arm is roughly 59% of the rate in the control arm at any given time point. The shape k = 1.4 > 1 implies an increasing hazard — the risk of progression rises over time in both arms — which is biologically plausible for a progressive condition where accumulated disease burden increases susceptibility.

The practical interpretation for a clinical or payer audience: patients on treatment typically go about 50% longer before their disease progresses, with the hazard of progression rising over time in both arms but at a substantially lower rate in the treated arm at any instant.

HTA extrapolation role and the NICE TSD candidate set

In health technology assessment, the standard workflow (Latimer 2013; NICE DSU TSD 14) requires fitting six candidate parametric distributions — exponential, Weibull, Gompertz, log-normal, log-logistic, and generalized gamma — to observed trial or registry data and selecting among them based on AIC/BIC, visual fit to the Kaplan-Meier curve, and clinical plausibility of the projected hazard shape. The Weibull is almost always in this candidate set and is often selected when the smoothed hazard is monotone.

The critical caution is that tail behavior dominates the lifetime QALY estimate because most of the area under the survival curve lies beyond the observed data. A Weibull with k slightly above 1 (modestly increasing hazard) and a log-normal (decreasing hazard after a peak) can fit 24 months of trial data with nearly identical AIC yet project mean survival differing by years. The analyst must justify the hazard shape in the unobserved tail on clinical grounds — not just on in-sample fit — and present the alternative distributions as scenarios in probabilistic sensitivity analysis. See the survival-extrapolation-hta-rwe entry for the full workflow.

Checking Weibull fit: the log(−log S(t)) vs log(t) diagnostic

The canonical graphical check is based on the cumulative hazard. If S(t) = exp(−(t/λ)^k), then log(−log S(t)) = k · log(t) − k · log(λ). Plotting the estimated log(−log S(t)) — computed from the Kaplan-Meier — against log(t) should yield an approximately straight line with slope k and intercept −k · log(λ) if the Weibull assumption holds. Systematic curvature (a convex or concave deviation) indicates the hazard is not a power law of time, and a more flexible model (log-normal, log-logistic, generalized gamma, or Royston-Parmar spline) may be needed. In R, PROC LIFETEST with the LOGLOGS option produces this plot directly; in R, survfit objects can be post-processed; in Python, lifelines provides diagnostic plots via WeibullFitter.

Pros, cons, and trade-offs

Pros of the Weibull distribution: - Parsimonious: only two parameters, straightforward to estimate by maximum likelihood. - Flexible monotone hazard shapes: k < 1, k = 1 (exponential special case), k > 1. - Dual PH/AFT structure: the only standard distribution fitting both frameworks, enabling hazard ratio reporting and AFT time-ratio reporting from one fit, with exact conversion HR = TR^(−k). - Closed-form survival function, hazard, and quantile function: computationally convenient for HTA cost-effectiveness models, probabilistic sensitivity analysis, and simulation. - Nests the exponential (k = 1): a likelihood ratio test of H_0: k = 1 is a formal test of the exponential special case. - Well-implemented in all major software: lifelines WeibullFitter and WeibullAFTFitter (Python), survreg and flexsurvreg (R), PROC LIFEREG with DIST=WEIBULL (SAS).

Cons and limitations: - Monotone hazard constraint: the Weibull hazard is strictly decreasing, constant, or strictly increasing throughout follow-up. It cannot model a non-monotone hazard (a hazard that rises then falls). When the clinical process produces a non-monotone pattern, the Weibull will misfit and produce a biased tail extrapolation. - No cure fraction: the Weibull hazard is always positive as t approaches infinity, so it cannot accommodate a long-term survivor fraction. If data show a genuine survival plateau (immunotherapy response, curative surgery), a mixture cure model is required. - Two-parameter family: less flexible than the three-parameter generalized gamma (which nests the Weibull as a special case) or Royston-Parmar splines. When in-sample fit is poor, upgrading is the natural next step. - Parameterization ambiguity: different software uses different conventions for shape and scale; always verify the parameterization in use before reporting k.

When to use

Use the Weibull distribution as the primary parametric survival model when: - The clinical or biological mechanism suggests a monotone hazard: increasing for progressive disease, decreasing for early-failure patterns, constant for memoryless processes. - You need a single model that reports both a hazard ratio (for clinical audiences) and a time ratio or mean life-years gain (for health economic models) without fitting two separate models. - You are conducting HTA survival extrapolation and the Weibull is the candidate whose projected hazard is most clinically plausible in the unobserved tail beyond trial data. - The log(−log S(t)) vs log(t) diagnostic is approximately linear, confirming the power- law hazard assumption. - Parametric efficiency matters: when a parametric form is correct, the Weibull achieves lower variance estimates than the semi-parametric Cox model.

When NOT to use — and when it is actively misleading

Do not use the Weibull distribution when: - The smoothed hazard is non-monotone: if the hazard rises then falls (e.g., post- chemotherapy toxicity then a declining relapse rate), the Weibull cannot fit this shape and will extrapolate with a systematically wrong tail. Use log-logistic (non-monotone hazard), log-normal, or generalized gamma / Royston-Parmar spline instead. - A cure fraction is biologically plausible and supported by mature follow-up data: a genuine survival plateau invalidates the Weibull's assumption that the hazard is always positive. Use a mixture or non-mixture cure model. - The log(−log S(t)) vs log(t) diagnostic shows systematic curvature: this directly falsifies the Weibull power-law assumption on the observed data. Upgrade to a more flexible model before extrapolating. - Competing risks are present and the Weibull is applied naively to the event of interest after censoring the competing event: this produces a cause-specific hazard that can overstate cumulative incidence when competing mortality differs by arm. Pair the Weibull fit with a formal competing-risks analysis. - Flexible spline models fit materially better and the decision horizon is long: if the Royston-Parmar spline AIC is substantially lower and the hazard shape is clinically complex, the Weibull's rigidity introduces extrapolation error that dominates the cost-effectiveness estimate.

Worked example

Scenario

A pharmacoepidemiology team has fitted a Weibull model to time-to-disease-progression data from a small illustrative cohort of five patients with a progressive condition. The maximum-likelihood estimates yield shape parameter k = 2 and scale parameter λ = 2 months — values chosen because the resulting arithmetic is exact and verifiable by hand. The team wants to evaluate the model-implied survival probabilities at months 2 and 4, verify the rising hazard pattern that k = 2 > 1 implies, and confirm that the hazard doubles from month 2 to month 4 under these parameter values.

Dataset

Observed event times (time_months) and event indicators (1 = progression observed, 0 = censored before progression) for five illustrative patients. A Weibull model fit to these data yields shape k = 2 and scale λ = 2 months, chosen for arithmetic transparency; a real maximum-likelihood fit would optimize these numerically.

Steps

• The Weibull survival function is S(t) = exp(-(t/λ)^k). With k = 2 and λ = 2 months, this simplifies to S(t) = exp(-(t/2)**2). The shape k = 2 > 1 tells us immediately that the hazard will be increasing over time, consistent with a progressive disease where cumulative burden raises risk.

• Compute the survival probability at month 2. The exponent is (2/2)**2 = 1, so S(2) = exp(-1) ≈ 0.368. About 37% of patients remain progression-free at 2 months under this Weibull model.

• Compute the survival probability at month 4. The exponent is (4/2)**2 = 4, so S(4) = exp(-4) ≈ 0.018. Nearly all patients have experienced progression by month 4 — only about 2% remain event-free.

• The Weibull hazard function with k = 2 and λ = 2 simplifies to h(t) = t/2. At month 2: h(2) = 2/2 = 1 event per month. At month 4: h(4) = 4/2 = 2 events per month. The hazard doubles between month 2 and month 4, confirming the monotone rising pattern that k = 2 > 1 implies.

• In an AFT two-arm comparison, the software would report a time ratio (TR) for the treatment covariate. A TR of 1.50 would mean the treated arm's progression times are stretched by 50% relative to control. The equivalent Weibull hazard ratio is HR = TR^(-k); with TR = 1.50 and k = 2, HR would be 1.50^(-2) = 1/(1.50**2) ≈ 0.44, meaning the treated arm progresses at roughly 44% of the control arm's instantaneous rate at any given moment.

Result

With k = 2 and λ = 2 months, (2/2)2 = 1 and (4/2)2 = 4, giving S(2) = exp(-1) ≈ 0.368 and S(4) = exp(-4) ≈ 0.018. Hazard at month 2: h(2) = 2/2 = 1 event per month. Hazard at month 4: h(4) = 4/2 = 2 events per month. The hazard doubles, confirming monotone increasing risk consistent with k = 2 > 1. A Weibull AFT coefficient of 0.405 implies a time ratio exp(0.405) ≈ 1.50 (treated arm goes 50% longer before progression); the equivalent PH hazard ratio is HR = TR^(-k) with the shape k used as the exponent.

Runnable example

import numpy as np
from scipy import stats
from lifelines import WeibullFitter, WeibullAFTFitter

# ── 1. Marginal Weibull fit (no covariates) ──────────────────────────────────────────
# df must have columns: time (float > 0), event (int: 1=occurred, 0=censored)
wf = WeibullFitter()
wf.fit(df["time"], event_observed=df["event"])

# lifelines parameterizes as S(t) = exp(-(t/lambda_)**rho_)
k      = wf.rho_       # shape parameter (rho in lifelines)
lam    = wf.lambda_    # scale parameter (lambda in lifelines)
print(f"Shape k: {k:.4f}   Scale λ: {lam:.4f}   AIC: {wf.AIC_:.2f}")

# Model check: is the hazard increasing (k>1), constant (k=1), or decreasing (k<1)?
if k > 1:
    print("Increasing hazard — consistent with progressive disease or wear-out.")
elif k < 1:
    print("Decreasing hazard — consistent with early-failure / burn-in pattern.")
else:
    print("Constant hazard — Weibull reduces to exponential (memoryless process).")

# ── 2. Compute S(t) and h(t) at specified time points ───────────────────────────────
times = np.array([1.0, 2.0, 4.0])
# scipy.stats.weibull_min uses (c=k, scale=lambda) for S(t) = exp(-(t/lambda)**k)
dist = stats.weibull_min(c=k, scale=lam)
surv_probs = dist.sf(times)          # survival function = 1 - CDF
hazards    = dist.pdf(times) / dist.sf(times)   # h(t) = f(t)/S(t)
for t, s, h in zip(times, surv_probs, hazards):
    print(f"  t={t:.1f}: S(t)={s:.4f}   h(t)={h:.4f}")

# ── 3. Log(−log S(t)) diagnostic: should be linear in log(t) if Weibull holds ───────
import pandas as pd
from lifelines import KaplanMeierFitter
kmf = KaplanMeierFitter()
kmf.fit(df["time"], event_observed=df["event"])
km_sf    = kmf.survival_function_
log_t    = np.log(km_sf.index.values + 1e-9)
log_log_s = np.log(-np.log(km_sf["KM_estimate"].values + 1e-9))
# Plot log_log_s vs log_t: a straight line supports Weibull; curvature suggests log-normal
# or log-logistic. Slope of the linear fit estimates k.

# ── 4. AFT model with treatment covariate ────────────────────────────────────────────
# df must additionally have column arm (0=control, 1=treated)
aft = WeibullAFTFitter()
aft.fit(df[["time", "event", "arm"]], duration_col="time", event_col="event")
aft.print_summary(decimals=4)

# Extract the time ratio (TR) for the arm coefficient
# In lifelines AFT output, lambda_ ancillary row = log(scale), covariate rows = log(TR)
# exp(arm coefficient in lambda_ rows) = time ratio
# lifelines WeibullAFTFitter stores rho_ as log(shape), so shape k = exp(rho_Intercept)
# NOT 1/exp(...); that inversion would produce 1/k and give a wrong HR.
arm_coef = aft.params_["lambda_"]["arm"]        # log time ratio for arm covariate
tr       = np.exp(arm_coef)
k_aft    = np.exp(aft.params_["rho_"]["Intercept"])   # shape: rho_ is log(shape)
hr_from_aft = tr ** (-k_aft)
print(f"\nAFT time ratio (TR): {tr:.4f}  ->  Weibull HR = TR^(-k) = {hr_from_aft:.4f}")

library(survival)
library(flexsurv)

# ── 1. AFT parameterization via survreg ──────────────────────────────────────────────
# survreg(dist="weibull") fits log T = mu + gamma*x + sigma*eps (Gumbel error)
# scale output = sigma = 1/k; intercept = log(lambda)
fit_aft <- survreg(Surv(time, event) ~ arm, data = df, dist = "weibull")
summary(fit_aft)

sigma <- fit_aft$scale           # sigma = 1/k
k_aft <- 1 / sigma               # shape parameter
lam   <- exp(coef(fit_aft)["(Intercept)"])   # scale parameter lambda
cat(sprintf("Shape k = %.4f   Scale λ = %.4f\n", k_aft, lam))

# Time ratio for arm covariate (note survreg sign convention)
# In survreg, the arm coefficient is log(TR) with positive = longer survival
arm_coef <- coef(fit_aft)["armtreated"]
TR <- exp(arm_coef)              # time ratio
HR <- TR ^ (-k_aft)             # convert to hazard ratio via HR = TR^(-k)
cat(sprintf("Time ratio TR = %.4f -> HR = TR^(-k) = %.4f\n", TR, HR))

# ── 2. PH parameterization via flexsurvreg ───────────────────────────────────────────
# dist="weibullPH" directly parameterizes by the hazard ratio for covariates
fit_ph <- flexsurvreg(Surv(time, event) ~ arm, data = df, dist = "weibullPH")
print(fit_ph)    # exp(arm coef) = HR directly

# ── 3. Compute S(t) and h(t) at specified time points ───────────────────────────────
times <- c(1, 2, 4)
S_t <- exp(-(times / lam)^k_aft)                    # survival function
h_t <- (k_aft / lam) * (times / lam)^(k_aft - 1)  # hazard function
cat("\n  t   S(t)   h(t)\n")
for (i in seq_along(times))
  cat(sprintf("  %g  %.4f  %.4f\n", times[i], S_t[i], h_t[i]))

# ── 4. Log(−log S(t)) Weibull diagnostic ─────────────────────────────────────────────
km <- survfit(Surv(time, event) ~ 1, data = df)
plot(log(km$time), log(-log(km$surv)),
     xlab = "log(t)", ylab = "log(-log S(t))",
     main = "Weibull diagnostic: should be linear if Weibull holds")
abline(lm(log(-log(km$surv)) ~ log(km$time)), col = "blue")
# Slope of the fitted line estimates k; a straight line supports Weibull assumption.

# ── 5. Mean survival to lifetime horizon for HTA ─────────────────────────────────────
HORIZON <- 30      # years (adapt to time units in the data)
DISC    <- 0.035   # annual discount rate (NICE reference case)
grid    <- seq(0, HORIZON, by = 1/12)
S_grid  <- summary(fit_ph, t = grid, type = "survival", ci = FALSE)[[1]]$est
disc    <- exp(-DISC * grid)
disc_ly <- sum(diff(grid) * (head(S_grid * disc, -1) + tail(S_grid * disc, -1)) / 2)
cat(sprintf("\nDiscounted mean life-years to %d-year horizon: %.3f\n", HORIZON, disc_ly))