Expected survival and event probability under a gamma-distributed constant hazard

Signature

S_bar = (1 + s * t)^(-k); p_bar = 1 - S_bar

Function

Expected value of an uncertain cost, health outcome or model output

Computational function

• Computational function: expected event probability and deterministic bias from the mean and standard deviation of an uncertain hazard

  Takes the summary a source usually reports for an uncertain constant hazard, its mean and standard deviation, with a horizon, a cost per event and a cohort size, and returns the expected event probability and the bias of a run at the mean hazard. It fits the gamma shape and scale by the method of moments, applies HE-FM-EV-003 for the expected probability and HE-FM-TP-001 for the deterministic one, and adds the second-order approximation HE-FM-EV-002 as a check. The inputs differ from the formula's variables: the formula needs the gamma shape and scale, and the function derives them from the mean and standard deviation.

  Inputs and outputs: h_mean: Mean of the uncertain annual hazard; required, above zero. Unit: events per person-year.; h_sd: Standard deviation of the hazard; required, above zero. Unit: events per person-year.; t: Horizon; required, zero or above. Unit: years.; c_event: Cost per event; required. Unit: pounds.; n_cohort: Cohort size for counting events; default 1,000. Unit: people.; k: Fitted gamma shape. Unit: none.; s: Fitted gamma scale. Unit: events per person-year.; p_det: Event probability at the mean hazard. Unit: probability from 0 to 1.; p_bar: Expected event probability. Unit: probability from 0 to 1.; p_approx: Second-order approximation to the expected probability. Unit: probability from 0 to 1.; gap_events: Events overstated per cohort by the run at the mean hazard. Unit: events.; gap_cost: Cost overstated per person by the run at the mean hazard. Unit: pounds per person.

  Assumption: The hazard is constant over the horizon in each draw, and its uncertainty is represented by a gamma distribution matched to the reported mean and standard deviation. Another distribution with the same mean and standard deviation would give a somewhat different expected probability.

  Worked example (Fracture hazard with mean 0.10 and standard deviation 0.0707 over five years): The fitted gamma has shape 2 and scale 0.05. The probability at the mean hazard is about 0.3935, the expected probability 0.36 and the approximation about 0.356, so a run at the mean overstates fractures by about 33 per 1,000 people and cost by about £402 per person at £12,000 per fracture, as in the article. h_mean = 0.10; h_sd = 0.0707107; t = 5; c_event = 12000; n_cohort = 1000; k = 2; s = 0.05; p_det = 0.393469; p_bar = 0.36; p_approx = 0.355561; gap_events = 33.47; gap_cost = 401.63

  Worked example (Same hazard over ten years): Over 10 years the expected probability is about 0.5556 against 0.6321 at the mean hazard, a gap of about 77 per 1,000 people, larger than at 5 years, as the article describes for horizons up to about 20 to 25 years. The second-order approximation, about 0.540, now overshoots the correction. The figures are computed here for illustration. h_mean = 0.10; h_sd = 0.0707107; t = 10; c_event = 12000; n_cohort = 1000; k = 2; s = 0.05; p_det = 0.632121; p_bar = 0.555556; p_approx = 0.540151; gap_events = 76.57; gap_cost = 918.78

  Excel: =1-(1+(HazSD^2/HazMean)*Horizon)^(-(HazMean^2/HazSD^2)) returns the expected probability from named cells HazMean, HazSD and Horizon, held in ExpProbCF; =(1-EXP(-HazMean*Horizon))-ExpProbCF gives the bias in probability, held in GapProb, and GapProb multiplied by CostPerEvent and by CohortSize gives the cost and event gaps.

  R: ev_gamma_hazard <- function(h_mean, h_sd, t, c_event, n_cohort = 1000) { k <- h_mean^2 / h_sd^2; s <- h_sd^2 / h_mean; p_det <- 1-exp(-h_mean * t); p_bar <- 1-(1 + s * t)^(-k); p_approx <- p_det-0.5 * t^2 * exp(-h_mean * t) * h_sd^2; list(k = k, s = s, p_det = p_det, p_bar = p_bar, p_approx = p_approx, gap_events = n_cohort * (p_det-p_bar), gap_cost = c_event * (p_det-p_bar)) } Base R only; vectorised over t, so a vector of horizons traces the gap over time.

  Python: def ev_gamma_hazard(h_mean, h_sd, t, c_event, n_cohort=1000): k = h_mean**2 / h_sd**2; s = h_sd**2 / h_mean; p_det = 1-math.exp(-h_mean * t); p_bar = 1-(1 + s * t)**(-k); p_approx = p_det-0.5 * t**2 * math.exp(-h_mean * t) * h_sd**2; return {"k": k, "s": s, "p_det": p_det, "p_bar": p_bar, "p_approx": p_approx, "gap_events": n_cohort * (p_det-p_bar), "gap_cost": c_event * (p_det-p_bar)} Needs import math; returns the same values as the R function.

  Test (Simulated mean matches the closed form): With 10,000 gamma draws of the hazard in HazDraws, made with =GAMMA.INV(RAND(),HazMean^2/HazSD^2,HazSD^2/HazMean), the average of the converted draws lies within four standard errors of the closed form; the probability has a standard deviation of about 0.19 across draws in the article's example. Expected result: TRUE on almost every recalculation. Excel check: =ABS(SUMPRODUCT(1-EXP(-HazDraws*Horizon))/ROWS(HazDraws)-ExpProbCF)<4*0.19/SQRT(ROWS(HazDraws))

  Test (Fitted gamma reproduces the reported moments): The fitted shape (named FitShape, =HazMean^2/HazSD^2) times the fitted scale (named FitScale, =HazSD^2/HazMean) returns the mean hazard, and shape times scale squared returns the variance. Expected result: TRUE. FALSE shows the scale and rate confused. Excel check: =AND(ABS(FitShape*FitScale-HazMean)<1E-12,ABS(FitShape*FitScale^2-HazSD^2)<1E-12)

  Common error (Fitting the gamma with the standard deviation treated as the variance): Taking 0.0707 as the variance gives a shape of about 0.14 and a scale of about 0.71, a hazard far more dispersed than reported, and an expected five-year probability of about 0.19 instead of 0.36.

  Source: Balan TA, Putter H. A tutorial on frailty models. Statistical Methods in Medical Research. 2020;29(11):3424-3454. doi:10.1177/0962280220921889 (full text read). Section 2.3.1, Laplace transform of the gamma distribution and marginal survival as the Laplace transform of the cumulative hazard; Briggs A, Claxton K, Sculpher M. Decision Modelling for Health Economic Evaluation. Oxford: Oxford University Press; 2006. Chapter 4, Making decision models probabilistic (pp. 77-120) (chapter abstract read). Choice of distributions for parameters; fitting the gamma by the method of moments is a standard textbook step.

  k = h_mean^2 / h_sd^2; s = h_sd^2 / h_mean; p_det = 1 - exp(-h_mean * t); p_bar = 1 - (1 + s * t)^(-k); p_approx = p_det - 0.5 * t^2 * exp(-h_mean * t) * h_sd^2; gap_events = n_cohort * (p_det - p_bar); gap_cost = c_event * (p_det - p_bar)

Try this function

Implementations

• Excel

  Expected event probability under a gamma hazard from named cells

  With HazShape, HazScale and Horizon named, the formula returns the expected event probability, held in ExpProb. Excel GAMMA.DIST and GAMMA.INV take the same shape and scale arguments.

  =1-(1+HazScale*Horizon)^(-HazShape)

Assumptions

• Hazard constant over the horizon in every draw

  Each value of h applies unchanged from time 0 to t; the uncertainty is in the level of the hazard, not in its shape over time.

• Gamma distribution of the hazard in shape and scale form

  The hazard follows a gamma distribution with shape k and scale s. Excel GAMMA.DIST and Python random.gammavariate take the scale; R rgamma takes the rate as its third positional argument, so a rate needs converting to 1 over the rate before use here.

• Uncertain hazard shared by the whole cohort

  The average is over uncertainty in a hazard shared by every patient, the parameter or second-order uncertainty of the article, so the result is the expected cohort probability, not a probability known to apply to each patient. Read as variation in hazard between patients, the same formula gives the marginal survival of a gamma-frailty population.

Worked examples

• Gamma hazard with shape 2 and scale 0.05 over five years

  With a mean hazard of 0.10 a year and a variance of 0.005, the expected fracture-free survival over 5 years is 1.25 raised to the power minus 2, or 0.64, and the expected fracture probability 0.36, against 0.3935 at the mean hazard. At £12,000 per fracture the expected cost is £4,320 rather than about £4,722 per person, as in the article.

  k = 2; s = 0.05; t = 5; S_bar = 0.64; p_bar = 0.36

• Same gamma hazard over a horizon of 23 years

  Over 23 years the expected probability is about 0.7837 against about 0.8997 at the mean hazard, a gap of about 0.116 and the largest for this hazard; the article places the peak at a horizon of about 20 to 25 years.

  k = 2; s = 0.05; t = 23; S_bar = 0.21633; p_bar = 0.78367

• More precise hazard with the same mean of 0.10

  With shape 200 and scale 0.0005, the same mean of 0.10 but a standard deviation of about 0.0071, the expected probability over 5 years is about 0.3931, within 0.0004 of the deterministic 0.3935: the gap shrinks with a more precise hazard, as the article states. The figures are computed here for illustration.

  k = 200; s = 0.0005; t = 5; S_bar = 0.60691; p_bar = 0.39309

Common errors

• Converting the mean hazard to a probability and calling it the expected probability

  Applying 1 minus exp(minus h t) at the mean hazard gives 0.3935 in place of 0.36, about 393 rather than 360 fractures per 1,000 people and about £402 more cost per person, an overstatement of about 9 per cent of the expected value.

• Entering the gamma rate where the scale is expected

  With a rate of 20 entered as the scale, the formula returns an expected five-year probability of about 0.9999 instead of 0.36. The parameterisation of each software function has to be checked before the hazard distribution is set up.

• Averaging hazard draws before converting them to probabilities

  Averaging sampled hazards and converting the average reproduces the deterministic 0.3935. The conversion has to be applied to each draw before averaging, or the closed form used.

Sources

• Gamma Laplace transform behind the expected event probability under an uncertain hazard

  Balan TA, Putter H. A tutorial on frailty models. Statistical Methods in Medical Research. 2020;29(11):3424-3454. doi:10.1177/0962280220921889 (full text read). Section 2.3.1: the Laplace transform of a non-negative random variable Z is E[exp(minus cZ)], and the marginal survival function is the Laplace transform evaluated at the cumulative hazard. The gamma distribution with shape eta and rate theta has Laplace transform (theta over theta plus c) raised to the power eta, which with scale s equal to 1 over theta and c equal to t gives (1 + s t) raised to the power minus k.

  View source →

• Gamma moment generating function behind expected survival at an uncertain hazard

  Casella G, Berger RL. Statistical Inference. 2nd ed. Pacific Grove, CA: Duxbury; 2002 (reprinted Boca Raton: Chapman and Hall/CRC; 2024). Chapter 3, the gamma family: mean equal to shape times scale, variance equal to shape times scale squared, and moment generating function (1 minus scale times u) raised to the power minus shape, which at u equal to minus t gives the expected survival. Textbook result.

  View source →

Canonical Identity