Mortality models and estimation

What a law is

A mortality law is a smooth function of age with a handful of parameters, meant to capture the age pattern of death in a population. Write \(x\) for age at the start of a one-year interval. There are two standard ways to say how dangerous that age is, and this package uses both.

• The force of mortality (the hazard) \(\mu_x\) is the instantaneous rate at which a person aged \(x\) dies. Think of it as the slope of the survival curve, scaled by how many are left to take it.
• The death probability \(q_x\) is the probability that someone aged \(x\) dies before reaching \(x + 1\). Its complement \(p_x = 1 - q_x\) is the chance of seeing the year out.

The two are readings of one curve, tied together by the cumulative hazard \(H_x\) and the survivorship \(S_x\):

\[ H_x = \int_0^x \mu_t \, dt, \qquad S_x = e^{-H_x}, \qquad q_x = 1 - \frac{S_{x+1}}{S_x}, \qquad p_x = 1 - q_x . \]

In MortalityLaws every law is a function function(x, par) returning a list whose hx element holds one value per age. The FIT column of the catalogue below says what hx contains: in mu[x] laws it is the hazard \(\mu_x\), in q[x] laws it is the death probability \(q_x\). The fitting engine and LawTable() read hx and route it to the right place, so nothing outside the law definition has to know which kind of law is in play.

One convention is worth knowing before the table: all parameters are strictly positive. The optimiser works on the log scale (see the optimisation procedure), so where a published formula wants a negative coefficient, the sign is folded into the formula itself. Oppermann’s middle term is written \(-B\), the log-quadratic families carry \(-B_2 x^2\), and so on. Two formula letters are also spelled differently in R to avoid name clashes: the coefficient \(F\) of Thiele’s formula and the hump location \(F\) of the Heligman-Pollard family are both the parameter F_, and the Rogers-Planck hump centre \(U\) is U.

The catalogue of laws

availableLaws() returns exactly the table below, together with the type legend of the next section. The columns are the code to pass as law, the model name, the formula in the parameter names the law functions use, the lifespan type, the fitted quantity, whether fitting rescales the ages, the year the catalogue assigns to the model, and where the formula comes from.

CODE	NAME	Formula	TYPE	FIT	SCALE_X	YEAR	Reference
demoivre	De Moivre	\(\mu_x = \dfrac{1}{N - x}\)	6	mu[x]	FALSE	1725	DeMoivre (1725)
gompertz	Gompertz	\(\mu_x = A e^{Bx}\)	3	mu[x]	TRUE	1825	Gompertz (1825)
gompertz0	Gompertz	\(\mu_x = \dfrac{1}{\sigma} \exp\left\{\dfrac{x - M}{\sigma}\right\}\)	3	mu[x]	TRUE	-	Gompertz (1825)
invgompertz	Inverse-Gompertz	\(\mu_x = \dfrac{\frac{1}{\sigma} \exp\left\{-\frac{x - M}{\sigma}\right\}}{\exp\left(\exp\left\{-\frac{x - M}{\sigma}\right\}\right) - 1}\)	2	mu[x]	TRUE	-	Tabeau et al. (2001)
makeham	Makeham	\(\mu_x = A e^{Bx} + C\)	3	mu[x]	TRUE	1860	Makeham (1860)
makeham0	Makeham	\(\mu_x = \dfrac{1}{\sigma} \exp\left\{\dfrac{x - M}{\sigma}\right\} + C\)	3	mu[x]	TRUE	-	Makeham (1860)
opperman	Opperman	\(\mu_x = \dfrac{A}{\sqrt{x + 1}} - B + C\sqrt{x + 1}\)	1	mu[x]	FALSE	1870	Oppermann (1870)
thiele	Thiele	\(\mu_x = A e^{-Bx} + C \exp\left\{-\tfrac{1}{2} D (x - E)^2\right\} + F e^{Gx}\)	6	mu[x]	FALSE	1871	Thiele (1871)
neggompertz	Negative-Gompertz	\(\mu_x = A e^{-Bx}\)	1	mu[x]	FALSE	1871	Thiele (1871)
wittstein	Wittstein	\(q_x = \dfrac{1}{B} A^{-(Bx)^N} + A^{-(M - x)^N}\)	6	q[x]	FALSE	1883	Wittstein and Bumsted (1883)
steffensen	Steffensen	\(\mu_x = \dfrac{A + B C^x}{B C^{-x} + 1 + D C^x}\)	6	mu[x]	TRUE	1930	Steffensen (1930)
perks	Perks	\(\mu_x = \dfrac{A + B C^x}{1 + D C^x}\)	3	mu[x]	TRUE	1932	Perks (1932)
weibull	Weibull	\(\mu_x = \dfrac{1}{\sigma} \left(\dfrac{x}{M}\right)^{\frac{M}{\sigma} - 1}\)	1	mu[x]	FALSE	1939	Weibull (1951)
pareto_2	Pareto-II	\(\mu_x = \dfrac{A}{x + C}\)	1	mu[x]	FALSE	1954	Lomax (1954)
invweibull	Inverse-Weibull	\(\mu_x = \dfrac{\frac{1}{\sigma} \left(\frac{x}{M}\right)^{-\frac{M}{\sigma} - 1}}{\exp\left(\left(\frac{x}{M}\right)^{-\frac{M}{\sigma}}\right) - 1}\)	2	mu[x]	TRUE	-	Weibull (1951)
vandermaen	Van der Maen	\(\mu_x = A + Bx + Cx^2 + \dfrac{I}{N - x}\)	4	mu[x]	TRUE	1943	Tabeau et al. (2001)
vandermaen2	Van der Maen	\(\mu_x = A + Bx + \dfrac{I}{N - x}\)	5	mu[x]	TRUE	1943	Tabeau et al. (2001)
strehler_mildvan	Strehler-Mildvan	\(\mu_x = A e^{Bx} \exp\left\{-\dfrac{V}{B}\left(1 - e^{-Bx}\right)\right\}\)	3	mu[x]	TRUE	1960	Finkelstein (2012)
quadratic	Quadratic	\(\mu_x = A + Bx + Cx^2\)	5	mu[x]	TRUE	-	Tabeau et al. (2001)
beard	Beard	\(\mu_x = \dfrac{A e^{Bx}}{1 + K A e^{Bx}}\)	4	mu[x]	TRUE	1971	Beard (1971)
beard_makeham	Beard-Makeham	\(\mu_x = \dfrac{A e^{Bx}}{1 + K A e^{Bx}} + C\)	4	mu[x]	TRUE	1971	Beard (1971)
ggompertz	Gamma-Gompertz	\(\mu_x = \dfrac{A e^{Bx}}{1 + \frac{AG}{B}\left(e^{Bx} - 1\right)}\)	4	mu[x]	TRUE	1979	Vaupel et al. (1979)
siler	Siler	\(\mu_x = A e^{-Bx} + C + D e^{E x}\)	6	mu[x]	FALSE	1979	Siler (1979)
HP	Heligman-Pollard	\(\dfrac{q_x}{p_x} = A^{(x + B)^C} + D \exp\left\{-E \log(x/F)^2\right\} + G H^x\)	6	q[x]	FALSE	1980	Heligman and Pollard (1980)
HP2	Heligman-Pollard	\(q_x = A^{(x + B)^C} + D \exp\left\{-E \log(x/F)^2\right\} + \dfrac{G H^x}{1 + G H^x}\)	6	q[x]	FALSE	1980	Heligman and Pollard (1980)
HP3	Heligman-Pollard	\(q_x = A^{(x + B)^C} + D \exp\left\{-E \log(x/F)^2\right\} + \dfrac{G H^x}{1 + K G H^x}\)	6	q[x]	FALSE	1980	Heligman and Pollard (1980)
HP4	Heligman-Pollard	\(q_x = A^{(x + B)^C} + D \exp\left\{-E \log(x/F)^2\right\} + \dfrac{G H^{x^K}}{1 + G H^{x^K}}\)	6	q[x]	FALSE	1980	Heligman and Pollard (1980)
rogersplanck	Rogers-Planck	\(q_x = A_0 + A_1 e^{-Ax} + A_2 \exp\left\{B(x - U) - e^{-C(x - U)}\right\} + A_3 e^{Dx}\)	6	q[x]	FALSE	1983	Rogers and Planck (1983)
martinelle	Martinelle	\(\mu_x = \dfrac{A e^{Bx} + C}{1 + D e^{Bx}} + K e^{Bx}\)	6	mu[x]	FALSE	1987	Martinelle (1987)
makeham_logquad	Gompertz-Makeham	\(\mu_x = A_0 + K \exp\left\{B_1 x - B_2 x^2\right\}\)	5	mu[x]	TRUE	1988	Forfar et al. (1988)
gompertz_logquad	Gompertz-Makeham	\(\mu_x = K \exp\left\{B_1 x - B_2 x^2\right\}\)	5	mu[x]	TRUE	1988	Forfar et al. (1988)
carriere1	Carriere	\(\mu_x = \dfrac{P_1 S^{w}_x \mu^{w}_x + P_2 S^{i}_x \mu^{i}_x + P_3 S^{g}_x \mu^{g}_x}{P_1 S^{w}_x + P_2 S^{i}_x + P_3 S^{g}_x}\), with the Weibull component \(\mu^{w}_x = \dfrac{(x/M_1)^{M_1/\sigma_1 - 1}}{\sigma_1}\), \(S^{w}_x = \exp\{-(x/M_1)^{M_1/\sigma_1}\}\); the inverse-Weibull component \(\mu^{i}_x = \dfrac{(x/M_2)^{-M_2/\sigma_2 - 1}}{\sigma_2\,\bigl(\exp\{(x/M_2)^{-M_2/\sigma_2}\} - 1\bigr)}\), \(S^{i}_x = 1 - \exp\{-(x/M_2)^{-M_2/\sigma_2}\}\); and the Gompertz component \(\mu^{g}_x = \dfrac{1}{\sigma_3}\exp\left\{\dfrac{x - M_3}{\sigma_3}\right\}\), \(S^{g}_x = \exp\{-\exp\{-M_3/\sigma_3\}(\exp\{x/\sigma_3\} - 1)\}\)	6	q[x]	TRUE	1992	Carriere (1992)
carriere2	Carriere	\(\mu_x = \dfrac{P_1 S^{w}_x \mu^{w}_x + P_2 S^{i}_x \mu^{i}_x + P_3 S^{g}_x \mu^{g}_x}{P_1 S^{w}_x + P_2 S^{i}_x + P_3 S^{g}_x}\), with the Weibull component \(\mu^{w}_x = \dfrac{(x/M_1)^{M_1/\sigma_1 - 1}}{\sigma_1}\), \(S^{w}_x = \exp\{-(x/M_1)^{M_1/\sigma_1}\}\); the inverse-Gompertz component \(\mu^{i}_x = \dfrac{\exp\{-(x - M_2)/\sigma_2\}}{\sigma_2\,\bigl(\exp\{\exp\{-(x - M_2)/\sigma_2\}\} - 1\bigr)}\), \(S^{i}_x = \dfrac{1 - \exp\{-\exp\{-(x - M_2)/\sigma_2\}\}}{1 - \exp\{-\exp\{M_2/\sigma_2\}\}}\); and the Gompertz component \(\mu^{g}_x = \dfrac{1}{\sigma_3}\exp\left\{\dfrac{x - M_3}{\sigma_3}\right\}\), \(S^{g}_x = \exp\{-\exp\{-M_3/\sigma_3\}(\exp\{x/\sigma_3\} - 1)\}\)	6	q[x]	TRUE	1992	Carriere (1992)
kostaki	Kostaki	\(\dfrac{q_x}{p_x} = A^{(x + B)^C} + D \exp\left\{-\left(E_i \log(x/F)\right)^2\right\} + G H^x\)	6	q[x]	FALSE	1992	Kostaki (1992)
kannisto	Kannisto	\(\mu_x = \dfrac{A e^{Bx}}{1 + A e^{Bx}}\)	5	mu[x]	TRUE	1998	Thatcher et al. (1998)
kannisto_makeham	Kannisto-Makeham	\(\mu_x = \dfrac{A e^{Bx}}{1 + A e^{Bx}} + C\)	5	mu[x]	TRUE	1998	Thatcher et al. (1998)
scholey_shifted_power	Scholey-Shifted-Power	\(\mu_x = A (x + C)^{-B}\)	1	mu[x]	FALSE	2019	Scholey (2019)
scholey	Scholey	\(\mu_x = A (x + C)^{-B} e^{-Dx}\)	1	mu[x]	FALSE	2019	Scholey (2019)

A few notes on the table.

• The formulas mirror the MODEL strings of availableLaws(), and the parameter names are the ones the law functions accept, so they can be copied straight into parS. Exceptions in spelling: \(F\) is F_ in R (in Thiele’s formula and in the Heligman-Pollard hump alike), and in siler the third exponent is the parameter \(E\), not Euler’s number.
• opperman is evaluated at ages shifted by one year, which keeps \(A/\sqrt{x}\) finite at birth. Its middle term uses the negative branch, which is the branch mortality data occupy; the published sign is free.
• weibull is undefined at \(x = 0\). Age 0 is reported as missing and takes no part in any fit.
• steffensen is attributed to Steffensen (1930), an attribution that could not be verified against the paywalled source. It is the formula this package used to ship as perks, before the published Perks form was separated out.
• carriere1 and carriere2 mix survivorship curves with weights normalised to the simplex, so \(P_3 = 1 - P_1 - P_2 > 0\). Their hx is a numerical derivative: the law function builds the mixture survivorship \(S_x\), takes \(H_x = -\log S_x\) and returns the year-to-year increments of that cumulative hazard on the age grid, which is the quantity the FIT column labels q[x]. The analytic hazard is the mixture formula in the table above, the derivative of the same \(H_x\) (Carriere 1992).
• In kostaki, \(E_i\) is \(E_1\) below the cut age \(F\) and \(E_2\) above it, so the accident hump may be asymmetric.
• The two rows dated 1988 are the Gompertz-Makeham graduation forms GM(1,3) and GM(0,3); see the log-quadratic family.
• Where the original publication is not in the bibliography (Van der Maen 1943, Strehler-Mildvan 1960, and the reparameterised Gompertz and Makeham), the Reference column names the closest work the bibliography does carry.

Lifespan types

The TYPE column says where on the lifespan a law belongs. It is not decoration: a TYPE 5 law fitted over the whole age range will happily draw nonsense at age 5. availableLaws() ships the legend as its second component.

availableLaws()$legend
#>   TYPE                       Coverage
#> 1    1               Infant mortality
#> 2    2                  Accident hump
#> 3    3                Adult mortality
#> 4    4 Adult and/or old-age mortality
#> 5    5              Old-age mortality
#> 6    6                 Full age range

Each type below comes with one representative law, fitted to observed rates to show the shape. Every figure uses the female population of England and Wales in 2010, from the bundled ahmd data, and plots the fitted curve against the observed rates with plot(fit, which = "fit"); the fitted \(R^2\) appears in the figure header. The laws can all be called by name like this, outside of any fit; each returns its hazard or death probability in hx.

TYPE 1, infant mortality. The block of neggompertz, weibull, pareto_2, the two scholey laws and opperman. The hazard starts high and falls, steeply at first. neggompertz decays at a constant relative rate, which is a decent story for the post-neonatal months and a poor one for the first day of life; a power law lets the decay slow down, and scholey_shifted_power fitted to ages 0 to 10 of the 2010 rates tracks the fall from 0.0041 at birth to 0.00007 at age 10, with \(R^2 = 0.999\). One caveat belongs to the segment rather than to the fit: the infant curve is only fully identified at day-level resolution, where the shifted-power family of Scholey (2019) is the tool, and since Scholey’s daily series is not shipped with the package, this illustration has to use the single-year rates.

x   <- 0:10
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "scholey_shifted_power")
plot(fit, which = "fit")

TYPE 2, the accident hump. invgompertz and invweibull. The inverse-Weibull rises from birth to a peak near its location parameter \(M\) and declines after it, which makes it the law to reach for when the hump is serious. It is also a hard fit on a female schedule, where the hump is a mild bump: over ages 10 to 35 of the 2010 rates it settles far from the data, with a negative \(R^2\), so the figure shows its sibling instead. The inverse-Gompertz rises steeply through the young ages and levels off at \(1/\sigma\), and over this window it climbs smoothly past the bump, giving \(R^2 = 0.927\).

x   <- 10:35
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "invgompertz")
plot(fit, which = "fit")

TYPE 3, adult mortality. gompertz, gompertz0, makeham, makeham0, perks and strehler_mildvan. Gompertz’s claim to fame is that the log hazard is a straight line in age, which is very nearly true of human adults. The hazard is unbounded, so these laws are fitted over adult ages rather than the whole lifespan. Fitted to ages 40 to 80 of the 2010 rates, the line is straight on the log scale to the eye, and \(R^2 = 0.982\).

x   <- 40:80
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "gompertz")
plot(fit, which = "fit")

TYPE 4, adult and/or old-age mortality. vandermaen, beard, beard_makeham and ggompertz. The curve rises like a Gompertz through adulthood and then bends: towards a ceiling, or towards the divergence of vandermaen’s closing term. This is the block to reach for when the old-age deceleration matters. The figure fits ggompertz to ages 40 to 110 of the 2010 rates, and the range matters: over ages 40 to 100 alone the frailty term collapses to nearly zero and the same law returns a plain Gompertz, while the longer range gives \(R^2 = 0.929\) and a visible bend away from the straight line above age 100.

x   <- 40:110
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "ggompertz")
plot(fit, which = "fit")

TYPE 5, old-age mortality. vandermaen2, quadratic, kannisto, kannisto_makeham and the two log-quadratic laws. These are fitted from the adult ages up. kannisto is the field standard for exactly this job: a Gompertz-like rise that levels off at a ceiling of one. Fitted to ages 60 to 110 of the 2010 rates it gives \(R^2 = 0.924\); the last decade of the range is thin and noisy, and the ceiling is what keeps the fitted curve steady where the data stop being informative.

x   <- 60:110
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "kannisto")
plot(fit, which = "fit")

TYPE 6, the full age range. demoivre, thiele, wittstein, steffensen, siler, the Heligman-Pollard family, rogersplanck, martinelle, the two Carriere mixtures and kostaki. These are multi-term curves with an infant part, a hump where the data have one, and a rising old-age part. Siler’s three terms are the cleanest example: immaturity, background and senescence, added together. Fitted to ages 0 to 100 of the 2010 rates, one curve covers the infant fall, the childhood trough and the adult rise, with \(R^2 = 0.975\).

x   <- 0:100
mx  <- ahmd$mx[paste(x), paste(2010)]
fit <- MortalityLaw(x = x, mx = mx, law = "siler")
plot(fit, which = "fit")

Family biographies

Thirty-eight laws are easier to face in families. What follows is who wrote what, and what each family is for. Laws whose original is not in the bibliography lean on Tabeau et al. (2001), the standard review of parametric mortality models, as their umbrella reference.

The Gompertz-Makeham lineage

Gompertz (1825) proposed the oldest law in the catalogue: the force of mortality rises exponentially with age, \(\mu_x = A e^{Bx}\). It is a one-parameter story about aging, and it is uncannily good at describing adult mortality. The code gompertz0 is the same curve written through its mode \(M\) and dispersion \(\sigma\), which are readable straight off the fitted line.

Makeham (1860) added the age-independent constant, \(\mu_x = A e^{Bx} + C\), so that accidents, infections and other age-free risks have somewhere to live; makeham0 is its mode-and-dispersion form. The full memoir of the law is Makeham (1867); see Bibliographic notes for why the catalogue dates it 1860.

Perks (1932) put a logistic denominator on the exponential, \((A + B C^x)/(1 + D C^x)\), which flattens the rise at the oldest ages. That four-parameter logistic exists in two spellings in the catalogue, perks and beard_makeham, and beard, kannisto and kannisto_makeham are its two- and three-parameter cases (Beard 1971). They fit identical curves when the parameter counts match, so pick one and not several.

The Strehler-Mildvan model of vitality decline rounds out the block (Finkelstein 2012); it predicts a negative correlation between the Gompertz intercept and slope across populations, a claim worth checking rather than assuming.

Siler and the competing risks view

Siler (1979) wrote mortality as three independent hazards added together: immaturity, \(A e^{-Bx}\); background, \(C\); and senescence, \(D e^{Ex}\). The model came out of animal mortality, where competing risks are the natural language, and it fits human curves well. Its first term is the declining hazard Thiele proposed in 1871 (Thiele 1871), which the catalogue carries separately as neggompertz. If you want one law with a story for every age and five parameters, Siler is the usual answer.

The Heligman-Pollard family

Heligman and Pollard (1980) modelled the odds of dying, \(q_x/p_x\), as three terms: the decline of infancy, a hump shaped like a lognormal density, and the Gompertz-like rise of old age. It is the standard full-range law in modern demography, and the fit is genuinely good over the whole lifespan. The catalogue carries four versions. HP is the original; HP2, HP3 and HP4 rewrite the old-age term as a logistic, with HP3 and HP4 adding a parameter \(K\) that bends the curve further (Heligman and Pollard 1980).

Kostaki (1992) extended the family to nine parameters by giving the accident hump two dispersion parameters, one below and one above the cut age \(F\), so the hump can lean. Martinelle (1987) generalised Perks the other way, adding a linear term \(K e^{Bx}\) above the logistic plateau so that very old ages can keep rising. All of these are high-parameter models and they all want opt.method = "LF2"; see Fitting.

Old age: the Kannisto family

Thatcher et al. (1998) studied the force of mortality at ages 80 to 120, and the logistic they used, \(\mu_x = A e^{Bx} / (1 + A e^{Bx})\), is kannisto. It rises like a Gompertz at the younger old ages and levels off at a ceiling of one, which is exactly the behaviour observed in the sparse, noisy data at the top of a life table. It is the field standard for closing life tables, and LifeTable() uses it for its close and omega arguments. Add the constant \(C\) and you get kannisto_makeham.

Carriere mixtures

Carriere (1992) took a different tack: instead of one hazard, mix standard survivorship curves. carriere1 combines Weibull, inverse-Weibull and Gompertz components; carriere2 swaps the middle one for an inverse-Gompertz. The weights are normalised to the simplex as they are fitted, so the third weight stays positive. The result is a flexible whole-lifespan curve where each component has a recognisable job: childhood, hump, old age.

The infant laws of Scholey

Scholey (2019) traced the age trajectory of infant mortality in US register data and found two regimes: a power law right after birth, then a constant exponential decline. The product of the two is scholey, \(A (x + C)^{-B} e^{-Dx}\); switch the truncation off and you have scholey_shifted_power, \(A (x + C)^{-B}\). The family nests the other infant laws: \(D = 0\) gives the shifted power, \(B = 1\) with \(D = 0\) gives pareto_2, and \(B = 0\) gives neggompertz. The truncation parameter is only identifiable on day- or week-level data over the first year of life, which is a limitation to know before fitting; see How to get it wrong.

Pareto, Lomax and the power hazards

pareto_2, \(\mu_x = A/(x + C)\), is the hazard of the Pareto type II, or Lomax, distribution (Lomax 1954): a shifted power law with the exponent fixed at one. It is a workhorse for the infancy and childhood years, and it is the infancy term of the delay-and-compression model of Beer and Janssen (2016). Vaupel and Yashin (1983) showed the same hazard arises from a gamma-exponential frailty model, so the smooth curve has a population story behind it.

The log-quadratic family

makeham_logquad and gompertz_logquad write the adult component as \(K \exp(B_1 x - B_2 x^2)\), with and without a Makeham constant \(A_0\). The \(-B_2 x^2\) term bends the Gompertz straight line down at the oldest ages, so the hazard decelerates instead of running away. These are the GM(1,3) and GM(0,3) forms of the Gompertz-Makeham graduation family catalogued by Forfar et al. (1988), which is what the 1988 rows of the table refer to. Whether mortality really decelerates at the oldest ages is a question with a long argument attached: Finkelstein (2012) works through it in the context of the Strehler-Mildvan model, and Missov et al. (2015) rewrites the Gompertz hazard through the modal age at death, the parameterisation used by the custom-law example in the introduction. Note that the sign of \(B_2\) is fixed to the decelerating branch; the accelerating branch cannot be fitted.

Rogers-Planck

rogersplanck is a three-term death probability covering the whole lifespan, developed as a general schedule for model life tables (Rogers and Planck 1983). Its middle term, \(A_2 \exp\{B(x - U) - e^{-C(x - U)}\}\), is a Gompertz-shaped hump on the log scale around the centre \(U\). The companion PAA paper of 1984 (Rogers and Planck 1984) is a different work and both are kept; see Bibliographic notes.

Heterogeneity: a note on what a fitted curve describes

Every law in this catalogue describes a population, not a person. Mix people with different frailties and the aggregate hazard flattens at old ages even if every individual’s hazard runs away: this is the point of Vaupel et al. (1979), whose gamma-frailty Gompertz is exactly the ggompertz law in the catalogue. The dynamics of such mixed populations are worked out in Vaupel and Yashin (1983). So a decelerating fit is evidence about the curve, and only indirectly about aging.

Choosing among the families

Pick the lifespan type first, so the law is used where it is defined. Then pick the family whose terms match the story you want to tell (immaturity, hump, senescence), and only then worry about the parameter count. A nine-parameter law will fit any single mortality schedule better than a two-parameter one; whether it describes anything is a separate question. Tabeau et al. (2001) reviews the field, and Harper (1936) remains useful background for the infant laws, though it maps to no code in the catalogue.

Bibliographic notes

Two date conflicts run through this documentation, and they are resolved here once.

• Makeham. availableLaws() prints 1860 next to the Makeham formula, so the catalogue table cites Makeham (1860), the short Assurance Magazine note. Makeham (1867), “On the law of mortality”, is the full memoir of the law, and that is the work cited in the narrative above.
• Siler. availableLaws() prints 1979 next to the Siler formula, so the catalogue table cites Siler (1979), the competing-risk paper in Ecology. The key siler1983 is kept in the bibliography but deliberately uncited, because what it contains is not verified; if it ever is, it may be cited in the Siler note above and nowhere else.

Three smaller points. The catalogue dates the Weibull law 1939; the bibliography carries the 1951 journal publication (Weibull 1951). The Rogers-Planck entry follows the 1983 IIASA working paper (Rogers and Planck 1983), with Rogers and Planck (1984) kept as the separate PAA meeting paper. And the steffensen formula is attributed to Steffensen (1930) (Steffensen 1930) on the strength of that attribution alone.

Fitting a law with MortalityLaw

Everything so far is vocabulary. MortalityLaw() is where the data comes in.

MortalityLaw(x, Dx = NULL, Ex = NULL, mx = NULL, qx = NULL,
             law = NULL, opt.method = "LF2", parS = NULL,
             fit.this.x = x, custom.law = NULL, show = FALSE, ...)

• x is the age vector, one value per age interval, at the start of the interval.
• The data is exactly one of three cases: death counts with exposures (Dx + Ex), central death rates (mx), or death probabilities (qx). Give two cases and the engine stops. A matrix or data frame of Dx (or mx, or qx) with several columns fits one model per column.
• law is a code from the catalogue, and custom.law is your own function; exactly one of the two.
• opt.method chooses what “best fit” means. There are eight options, discussed below.
• parS are starting values for the parameters, as a named, strictly positive numeric vector. Leave it alone and the built-in starting values of the law are used.
• fit.this.x restricts the optimisation to a subset of the ages. The fitted values are still returned over all of x.
• show = TRUE displays a progress bar, which is mostly useful when many columns are fitted at once.
• ... is passed to or from other methods.

The workhorse data set for the examples in these vignettes is the bundled ahmd: death counts, exposures and rates for England and Wales females, ages 0 to 110, for 1850, 1900, 1950 and 2010 (HMD 2026). Here is a Gompertz fitted to adult ages in 2010.

year     <- 2010
ages     <- 45:90
deaths   <- ahmd$Dx[paste(ages), paste(year)]
exposure <- ahmd$Ex[paste(ages), paste(year)]

fit <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "gompertz",
  opt.method = "poissonL"
)
summary(fit)
#> Gompertz model: mu[x] = A exp[Bx]
#> Fitted values: mx  |  ages 45-90  |  fitted on 45-90 (46 of 46 ages)
#> 
#> Call:
#> MortalityLaw(x = ages, Dx = deaths, Ex = exposure, law = "gompertz", 
#>     opt.method = "poissonL")
#> 
#> Coefficients:
#>   estimate
#> A   0.0009
#> B   0.1093
#> 
#> Fit:
#>   method poissonL  |  optimiser converged in 12 iterations
#>   deviance 1721 on 44 degrees of freedom  |  dispersion 40.28
#>   R-squared 0.9894  |  RMSE 0.003846
#> 
#> Goodness of fit:
#>    logLik       AIC       BIC 
#> -851178.8 1702361.6 1702365.3 
#> 
#> Residuals:
#>             Min. 1st Qu. Median   Mean 3rd Qu.    Max.
#> raw      -0.0031 -0.0012 0.0003 0.0005  0.0006  0.0217
#> deviance -8.8496 -4.9025 2.1761 0.7966  6.1234 15.7354

The eight objectives

What the optimiser minimises is chosen with opt.method. Write \(\mu\) for the fitted value and \(\nu\) for the observed one: \(\nu = D_x/E_x\) when counts and exposures are supplied, and \(\nu = m_x\) or \(\nu = q_x\) when rates or probabilities are. \(D_x\) is the death count and \(E_x\) the exposure to risk. Two of the eight are likelihoods, and six are loss functions.

\[ \begin{aligned} L_{\text{poissonL}} &= -\bigl[D_x \log \mu - \mu E_x\bigr], \\ L_{\text{binomialL}} &= -\bigl[D_x \log(1 - e^{-\mu}) - (E_x - D_x)\mu\bigr], \\ L_{\text{LF1}} &= \Bigl(1 - \frac{\mu}{\nu}\Bigr)^2, \qquad L_{\text{LF2}} = \Bigl(\log\frac{\mu}{\nu}\Bigr)^2, \\ L_{\text{LF3}} &= \frac{(\nu - \mu)^2}{\nu}, \qquad L_{\text{LF4}} = (\nu - \mu)^2, \\ L_{\text{LF5}} &= (\nu - \mu) \log\frac{\nu}{\mu}, \qquad L_{\text{LF6}} = \lvert \nu - \mu \rvert . \end{aligned} \]

The total loss is the sum over the fitted ages; a term that is not finite is replaced by a penalty of \(10^5\), and a fitted hazard larger than one is capped at one. When rates or probabilities are the input, the observed rates stand in for the death counts in the two likelihoods. availableLF() prints the same formulas at the console.

availableLF()
#> 
#> Loss functions available in the package:
#> 
#>  LOSS FUNCTION                                  CODE     
#>  L = -[Dx * log(mu) - mu*Ex]                    poissonL 
#>  L = -[Dx * log(1 - exp(-mu)) - (Ex - Dx)*mu]   binomialL
#>  L =  [1 - mu/ov]^2                             LF1      
#>  L =  log[mu/ov]^2                              LF2      
#>  L =  [(ov - mu)^2]/ov                          LF3      
#>  L =  [ov - mu]^2                               LF4      
#>  L =  [ov - mu] * log[ov/mu]                    LF5      
#>  L =  abs(ov - mu)                              LF6      
#> 
#> LEGEND:
#> Dx: Death counts
#> Ex: Population exposed to risk
#> mu: Estimated value
#> ov: Observed value
#> 
#> HINT: Most loss functions work well with 'poissonL'. However, for complex mortality laws like Heligman-Pollard (HP), a better fit can be obtained using other loss functions (e.g. 'LF2'). You are strongly encouraged to test different options before deciding on the final version. The results might be slightly different.

Choosing is mostly empirical. The two likelihoods are the principled choice for count data, and poissonL works well for most laws. The loss functions are for when the fit needs help in a particular corner of the curve, and the high-parameter laws of the Heligman-Pollard family are the reason the option exists: they fit reliably under LF2 and less reliably elsewhere, and the package says so at the console when you try. Test a couple of objectives on your data before deciding; the results will differ slightly.

The optimisation procedure

Starting values. Every catalogue law carries built-in starting values, visible by calling the law without par.

gompertz(x = 45:90)$par
#>      A      B 
#> 0.0002 0.1300

Pass parS to supply your own. The names must match the law’s parameters exactly and all values must be positive; the engine validates both before the optimiser runs. The starting point of the search is log(parS) either way.

fit_parS <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "gompertz",
  opt.method = "poissonL",
  parS       = c(A = 0.001, B = 0.05)
)
rbind(default = coef(fit), parS = coef(fit_parS))
#>                    A         B
#> default 0.0008778197 0.1093225
#> parS    0.0008778188 0.1093226

Two very different starting points, one optimum. That is the usual outcome for these laws, not a guarantee; a fit that lands somewhere else is telling you the objective has more than one basin.

The log transform. Parameters are optimised on the log scale: the objective calls the law as fn(x, exp(par)). This is why all parameters are strictly positive, why a formula needing \(-0.5\) carries the sign in its own text instead, and why an extreme probe simply underflows rather than producing a negative hazard.

The optimiser. All laws are fitted with nlminb and its PORT routines, with eval.max and iter.max at 5000. The single exception is invweibull, which uses optim with Nelder-Mead. There is no user-facing switch, and the coefficients returned are always on the original parameter scale.

Reading the convergence output. A non-zero convergence code raises a warning of the form “MortalityLaw: optimisation did not converge (code k)”, and the full optimiser object is kept in opt.diagnosis for inspection. Note that nlminb attaches a message such as “relative convergence (4)” even on success, so a message with convergence = 0 is routine. summary() reports the iteration count and the convergence status in one line.

Age scaling. Laws flagged SCALE_X = TRUE are fitted on shifted ages,

\[ x_{\text{fit}} = x - \min(x) + 1, \]

so that the youngest fitted age becomes 1 and the exponentials stay tame. The coefficients therefore refer to scaled ages and are not comparable with the unscaled parameterisation, while the fitted curves are always returned on the original age scale. The same shift is reapplied whenever the fitted law is evaluated again, including predict() and LawTable(). The one trap: for a scaled law, LawTable() is only valid from the lower bound of the fitting range upwards.

The fitting window. fit.this.x fits a subset of the ages while keeping the full fitted curve. This is how you fit an adult law to adult ages without dropping the rest of the vector.

fit_window <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "gompertz",
  opt.method = "poissonL",
  fit.this.x = 60:90
)
range(fit_window$input$fit.this.x)
#> [1] 60 90
length(fit_window$fitted.values)
#> [1] 46

Goodness of fit

A finished MortalityLaw object carries three kinds of residual (raw, deviance and Pearson), the deviance, the degrees of freedom and dispersion, and goodness.of.fit with the log-likelihood, AIC and BIC.

Information criteria. For the two likelihood objectives,

\[ \text{logLik} = -L, \qquad \text{AIC} = 2k - 2\,\text{logLik}, \qquad \text{BIC} = \log(n)\, k - 2\,\text{logLik}, \]

with \(k\) the number of parameters and \(n\) the number of fitted ages. All three are NaN for the six loss functions, because a sum of squares is not a likelihood and pretending otherwise only produces numbers that look like AICs. One caveat the engine is honest about: the log-likelihood drops data-only additive constants, so its absolute value differs from what glm reports, while comparisons between fits of the same data are unaffected.

fit$goodness.of.fit
#>    logLik       AIC       BIC 
#> -851178.8 1702361.6 1702365.3
fit$df
#>     n.param df.residual  dispersion 
#>     2.00000    44.00000    40.27795

Counts or rates. The deviance and the dispersion depend on the data case. With counts and exposures, everything follows the Poisson definitions: the Pearson residual is \((D_x - \mu E_x)/\sqrt{\mu E_x}\), the deviance is the sum of squared Poisson deviance residuals (the quantity poissonL minimises), and

\[ \text{dispersion} = \frac{\sum \text{Pearson}^2}{\text{rdf}}, \]

which is about one for a correctly specified Poisson model and larger when the data are noisier than Poisson, as single-age death counts usually are. With rates or probabilities there is no count likelihood, so the residuals are log residuals, \(\log \nu - \log \mu\); the deviance is their sum of squares, and the dispersion is their mean square.

fit_rate <- MortalityLaw(
  x          = ages,
  mx         = ahmd$mx[paste(ages), paste(year)],
  law        = "gompertz",
  opt.method = "LF2"
)
fit_rate$goodness.of.fit   # NaN: LF2 is a loss, not a likelihood
#> logLik    AIC    BIC 
#>    NaN    NaN    NaN
c(deviance = fit_rate$deviance, dispersion = fit_rate$dispersion)
#>   deviance dispersion 
#> 0.53171290 0.01208438

Residuals. Numbers tell you the fit is close; the residual panels tell you where it is not. plot(fit) draws the observed points against the fitted line plus four diagnostics (residuals against age, against the fitted values, a normal Q-Q plot and a histogram). The worked example is in the introduction to the package.

How to get it wrong

Most of the ways to get it wrong produce output that looks fine. These are the ones the package tries to warn you about, and one or two it cannot.

• Fitting weibull at birth. The hazard is undefined at \(x = 0\), so age 0 is dropped from the fit with a warning. Start the fit at age 1 and save yourself the surprise.
• Extrapolating demoivre. The hazard \(1/(N - x)\) diverges at \(N\), and the fit always places \(N\) just above the top fitted age. Predicting past the fitted range gives a negative hazard, and every demoivre fit warns about it, naming the fitted \(N\).
• Trusting scholey on year data. The truncation parameter \(D\) is only identified on day- or week-level ages over the first year of life. On single-year ages it collapses to the optimisation boundary and the model reduces to scholey_shifted_power, with a warning. The fit is still fine; it is just a different model than the one you asked for.
• Ignoring the kostaki guard. If the two hump dispersions drift more than a factor of 50 apart, the engine pulls \(E_2\) back to \(E_1/50\) before the hazard is evaluated. This is a deliberate hack to stop an artificial jump at the cut age; if the reported coefficients look oddly tidy, this is why.
• Fitting the HP family with a likelihood. HP, HP2, HP3, HP4 and kostaki are reliably estimated under opt.method = "LF2" and not under the others. The package prints a hint at the console when you try; the hint is correct.
• Plugging scaled coefficients into the unscaled formula. Laws flagged SCALE_X = TRUE are calibrated on shifted ages, \(x - \min(x) + 1\), so their coefficients are not comparable with the unscaled parameterisation and must not be read into the formula with ages as they come. The recipe is one line: shift the same way before you evaluate. The Age scaling section of the introduction shows the worked demonstration.
• Reporting AIC from a loss function. goodness.of.fit is NaN there on purpose. Do not fill it in by hand or compare such numbers across fits.
• Fitting perks and beard_makeham to different data sets and comparing them. They are the same four-parameter curve written two ways.

Custom laws

Anything not in the catalogue can be fitted by handing a function to custom.law. The contract is small.

my_law <- function(x, par) {
  # compute one mortality value per age, on the scale of your data
  list(hx = hx, par = par)      # optionally add Sx = survivorship
}

• The function takes the age vector x and a parameter vector par, and returns a list. The hx element is required and must hold one value per age; par is echoed back with the (possibly normalised) parameters, and an Sx element with the survivorship is optional. The mixture laws of the catalogue use Sx internally, and so may yours.
• hx must be on the scale of the data you fit. The engine compares hx directly with \(D_x/E_x\), \(m_x\) or \(q_x\). A hazard fitted to rate data is correct only because \(m_x\) is a rate; if you fit probabilities, return probabilities. This is what the FIT column means for your own law.
• Starting values come from my_law(1)$par. The engine evaluates the function at x = 1 and reads par out of the returned list, so the default value of par is where your search starts.
• Parameters are optimised on the log scale, so every parameter must be positive for the fit to be reachable. Fold any sign into the formula, as the built-in laws do.
• A custom law is always treated as a scaled one. Ages are shifted by \(x - \min(x) + 1\) before your function sees them, exactly as for SCALE_X = TRUE laws, and the fitted curve is reported on the original scale.
• Values of hx that are non-finite or non-positive are set aside and penalised, so a law that misbehaves part of the time drags the fit instead of crashing it. Return a clean vector and the optimiser will thank you.
• LawTable() takes catalogue codes only; to build life tables from a custom law, evaluate it yourself and pass the result to LifeTable().

A worked custom-law fit, with a Gompertz written through the modal age at death, is in the introduction to the package.

Session info and references

Everything on this page was produced with:

sessionInfo()
#> R version 4.6.0 (2026-04-24 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: Europe/Budapest
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] MortalityLaws_3.0.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60      
#>  [5] cachem_1.1.0    parallel_4.6.0  knitr_1.51      htmltools_0.5.9
#>  [9] rmarkdown_2.31  lifecycle_1.0.5 cli_3.6.6       sass_0.4.10    
#> [13] jquerylib_0.1.4 compiler_4.6.0  httr_1.4.8      tools_4.6.0    
#> [17] pbapply_1.7-4   evaluate_1.0.5  bslib_0.12.0    yaml_2.3.12    
#> [21] otel_0.2.0      rlang_1.3.0     jsonlite_2.0.0

References

Beard, Robert E. 1971. “Some Aspects of Theories of Mortality, Cause of Death Analysis, Forecasting and Stochastic Processes.” Biological Aspects of Demography 999: 57–68.

Beer, J. de, and F. Janssen. 2016. “A New Parametric Model to Assess Delay and Compression of Mortality.” Population Health Metrics 14 (1): 46. https://doi.org/10.1186/s12963-016-0113-1.

Carriere, Jacques F. 1992. “Parametric Models for Life Tables.” Transactions of the Society of Actuaries 44: 77–99.

DeMoivre, A. 1725. “Annuities on Lives: Or, the Valuation of Annuities Upon Any Number of Lives as Also of Reversions.” William Pearson, London.

Finkelstein, M. 2012. “Discussing the Strehler-Mildvan Model of Mortality.” Demographic Research 26 (9): 191–206. https://doi.org/10.4054/DemRes.2012.26.9.

Forfar, D. O., J. J. McCutcheon, and A. D. Wilkie. 1988. “On Graduation by Mathematical Formula.” Journal of the Institute of Actuaries 115 (1): 1–149.

Gompertz, Benjamin. 1825. “On the Nature of the Function Expressive of the Law of Human Mortality, and on a New Mode of Determining the Value of Life Contingencies.” Philosophical Transactions of the Royal Society of London 115: 513–83.

Harper, Floyd S. 1936. “An Actuarial Study of Infant Mortality.” Scandinavian Actuarial Journal 1936 (3-4): 234–70.

Heligman, Larry, and John H Pollard. 1980. “The Age Pattern of Mortality.” Journal of the Institute of Actuaries 107 (01): 49–80.

HMD. 2026. Human Mortality Database. Max Planck Institute for Demographic Research (Germany), University of California, Berkeley (USA), and French Institute for Demographic Studies (France). https://www.mortality.org/.

Kostaki, Anastasia. 1992. “A Nine-Parameter Version of the Heligman-Pollard Formula.” Mathematical Population Studies 3 (4): 277–88.

Lomax, K. S. 1954. “Business Failures: Another Example of the Analysis of Failure Data.” Journal of the American Statistical Association 49 (268): 847–52. https://doi.org/10.1080/01621459.1954.10501239.

Makeham, W. 1860. “On the Law of Mortality and Construction of Annuity Tables.” The Assurance Magazine and Journal of the Institute of Actuaries 8 (6): 301–10. https://doi.org/10.1017/S204616580000126X.

Makeham, William Matthew. 1867. “On the Law of Mortality.” Journal of the Institute of Actuaries (1866-1867) 13 (6): 325–58.

Martinelle, Sten. 1987. A Generalized Perks Formula for Old-Age Mortality.

Missov, Trifon I, Adam Lenart, Laszlo Nemeth, Vladimir Canudas-Romo, and James W Vaupel. 2015. “The Gompertz Force of Mortality in Terms of the Modal Age at Death.” Demographic Research 32: 1031–48.

Oppermann, Ludvig Henrik Ferdinand. 1870. “On the Graduation of Life Tables, with Special Application to the Rate of Mortality in Infancy and Childhood.” The Insurance Record Minutes from a meeting in the Institute of Actuaries.: 42.

Perks, Wilfred. 1932. “On Some Experiments in the Graduation of Mortality Statistics.” Journal of the Institute of Actuaries (1886-1994) 63 (1): 12–57.

Rogers, A., and F. Planck. 1983. MODEL: A General Program for Estimating Parametrized Model Schedules of Fertility, Mortality, Migration, and Marital and Labor Force Status Transitions. Working Paper WP-83-102. IIASA.

Rogers, A, and F Planck. 1984. “Parameterized Multistage Population Projections.” Working Paper for Presentation at the Annual Meeting of the Population Association of America, Minnesota, May 3-5.

Scholey, J. 2019. The Age-Trajectory of Infant Mortality in the United States: Parametric Models and Generative Mechanisms. PAA Annual Conference, Austin.

Siler, W. 1979. “A Competing-Risk Model for Animal Mortality.” Ecology 60 (4): 750–57. https://doi.org/10.2307/1936612.

Steffensen, JF. 1930. “Infantile Mortality from an Actuarial Point of View.” Skandinavisk Aktuarietidskrift 13: 272–86.

Tabeau, Ewa, Anneke van den Berg Jeths, and Christopher Heathcote. 2001. Forecasting Mortality in Developed Countries: Insights from a Statistical, Demographic and Epidemiological Perspective. Vol. 9. Springer Science & Business Media.

Thatcher, A Roger, Väinö Kannisto, and James W Vaupel. 1998. The Force of Mortality at Ages 80 to 120.

Thiele, Thorvald Nicolai. 1871. “On a Mathematical Formula to Express the Rate of Mortality Throughout the Whole of Life.” Journal of the Institute of Actuaries and Assurance Magazine. Translated by T.B. Sprague in 1872. 16 (5): 313–29.

Vaupel, J. W., and A. I. Yashin. 1983. The Deviant Dynamics of Death in Heterogeneous Populations. Research Report RR-83-1. IIASA.

Vaupel, J., K. G. Manton, and E. Stallard. 1979. “The Impact of Heterogeneity in Individual Frailty on the Dynamics of Mortality.” Demography 16 (3): 439–54. https://doi.org/10.2307/2061224.

Weibull, Waloddi. 1951. “A Statistical Distribution Function of Wide Applicability.” Journal of Applied Mechanics 18: 293–97.

Wittstein, Theodor, and DA Bumsted. 1883. “The Mathematical Law of Mortality.” Journal of the Institute of Actuaries and Assurance Magazine 24 (3): 153–73.