A short note on a bug that is easy to write, hard to see, and lives in a lot of hand-rolled survival and life-table code.
Fit a parametric mortality model (Gompertz, Gompertz–Makeham, Siler, …) to a life
table by maximum likelihood and you eventually need the probability that an
individual alive at the start of an age interval dies during it. If you have
Nx individuals at risk and observe Dx deaths, the cohort likelihood is
binomial:
dbinom(Dx, Nx, p, log = TRUE)The third argument p must be a probability in [0, 1]. The natural thing to
reach for is the model's hazard, so you compute the cumulative hazard over the
interval and pass it in:
dbinom(Dx, Nx, H, log = TRUE) # <- looks right, is wrongThat is the bug.
H is the cumulative hazard over the interval — the integral of the
instantaneous hazard μ(x):
H = ∫ μ(x) dx over the interval
It is a perfectly good quantity, but it is not a probability: it runs from 0 to ∞. The probability you actually want comes from it through the survival function:
- probability of surviving the interval:
S = exp(−H) - probability of dying in the interval:
q = 1 − exp(−H)
Because exp(−H) maps [0, ∞) onto (0, 1], q = 1 − exp(−H) is always a legal
probability. Raw H is not — not because a number above 1 is meaningless (a
cumulative hazard above 1 is fine), but because it is being used as a
probability where only [0, 1] is allowed.
For small H, 1 − exp(−H) ≈ H, so at young and middle ages the two agree to
several decimal places and every test you run looks fine. The disagreement only
grows in the tail — which, if you study ageing, is exactly the part you care
about, and exactly the part your quick sanity checks tend to miss.
With a Gompertz hazard μ(x) = a·e^(bx), a = 3×10⁻⁵, b = 0.13 per year, the raw
cumulative hazard crosses 1 at about age 80 and keeps climbing — H ≈ 3.9 at
90, ≈ 14 at 100. Fed to dbinom as a probability, each of those returns NaN,
the log-likelihood becomes NaN, and the optimiser stalls or wanders off. The
corrected curve bends over and saturates just below 1, as a death probability
must. Even below age 80, where nothing crashes, feeding raw H quietly
overstates the death probability and biases the parameter estimates.
The reason this is so easy to write is a life-table likelihood usually has two branches, and a binomial slot and a Poisson slot want genuinely different things:
# Cohort data (binomial): Dx deaths out of Nx alive at the start of the interval
# WRONG: dbinom(Dx, Nx, H, log = TRUE) # H used as a probability -> NaN once H > 1
# RIGHT: dbinom(Dx, Nx, 1 - exp(-H), log = TRUE) # q in [0, 1]
# Period / event-count data (Poisson): Dx deaths given exposure
# dpois(Dx, Lambda, log = TRUE) # Lambda = expected deaths = ∫ Y(t) μ(t) dt
# e.g. dpois(Dx, Ex * mx, log = TRUE) with person-years Ex and rate mxA binomial needs a bounded probability; a Poisson needs an expected count (a rate integrated over exposure), which is unbounded. Neither slot takes the raw cumulative hazard as-is.
One caution worth stating precisely, because it is itself a common shortcut: if
you write the Poisson mean as Nx * H with Nx the number alive at the start
of the interval, that is only a low-mortality approximation. The exact expected
death count is Nx * (1 − exp(−H)), and — no surprise — it diverges from Nx * H
in the same tail, for the same reason. Nx * H is exact only under particular
exposure assumptions (e.g. constant number at risk over the interval). Use
person-time exposure Ex * mx when you have it.
A cumulative hazard is not a probability, and it is not an expected count either.
When you convert a hazard to a probability, the bridge is the survival function:
q = 1 − exp(−H) (in code, -expm1(-H) is a touch more accurate for small H).
If a hazard is going straight into anything that expects a probability, that's
the line to check.
gompertz_death_prob.R in this gist regenerates the figure from scratch in base R.