Introduction to MortalityLaws

Marius D. Pascariu

2026-05-05

Package structure

MortalityLaws is an R package that provides tools for fitting and analyzing a wide range of parametric mortality models. It exploits various optimisation methods to handle complex models efficiently. The package can construct full and abridged life tables from different inputs and download data directly from the Human Mortality Database (HMD).

The core functions include:

Support functions such as availableLaws, availableLF, and availableHMD return information about the implemented mortality laws, loss functions, and HMD data availability, respectively. Standard generics (predict, plot, summary, fitted, residuals) work on MortalityLaw and LifeTable objects.

A small built-in dataset (ahmd) is provided for testing. All functions are documented in the usual R way; after loading the package with library(MortalityLaws), you can type ?MortalityLaw to see its help file.

Downloading HMD data

The ReadHMD function downloads data from the Human Mortality Database (University of California Berkeley, USA and Max Planck Institute for Demographic Research, Germany 2016). To use it, you need a valid HMD account.

library(MortalityLaws)
# Download Swedish death counts (ages 0–110, years 1751–2014)
HMD_Dx <- ReadHMD(
  what      = "Dx",
  countries = "SWE",
  interval  = "1x1",
  username  = "user@email.com",
  password  = "password",
  save      = FALSE
)

You can similarly download:

for over 40 countries and regions in various formats (1x1, 1x5, 5x1, 5x5).

Model fitting and diagnosis

Once you have data, fitting a model is straightforward. Below we fit a Heligman–Pollard (Heligman and Pollard 1980) model to Swedish mortality in 1950 using a Poisson likelihood (loss function "LF2").

year     <- 1950
ages     <- 0:100
deaths   <- ahmd$Dx[paste(ages), paste(year)]
exposure <- ahmd$Ex[paste(ages), paste(year)]

fit <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "HP",
  opt.method = "LF2"
)

Examining the fitted model

ls(fit)   # components of the fitted object
## [1] "coefficients"    "deviance"        "df"              "fitted.values"  
## [5] "goodness.of.fit" "info"            "input"           "opt.diagnosis"  
## [9] "residuals"

A summary of the fit gives parameter estimates and goodness‑of‑fit statistics:

summary(fit)
## Heligman-Pollard model: q[x]/p[x] = A^[(x + B)^C] + D exp[-E log(x/F)^2] + G H^x
## Fitted values: mx
## 
## Call: MortalityLaw(x = ages, Dx = deaths, Ex = exposure, law = "HP", 
##     opt.method = "LF2")
## 
## Deviance Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.0312 -0.0002  0.0000  0.0009  0.0003  0.0462 
## 
## Parameters:
##       A       B       C       D       E      F_       G       H 
##  0.0022  0.0146  0.1229  0.0009  2.7565 29.0080  0.0000  1.1141 
## 
## Residual standard error: 0.6915 on 93 degrees of freedom

Visual inspection is also easy:

plot(fit)

The default plot shows observed and fitted age‑specific death rates (or probabilities), the residuals, and the fitted parameters.

Fitting on a subset of ages

Sometimes you may want to fit the model only to a certain age range while still plotting over the full range. Use fit.this.x to specify the ages used during optimisation:

fit.subset <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "HP",
  opt.method = "LF2",
  fit.this.x = 0:65
)
plot(fit.subset)

The grey area indicates the ages that were actually used in the fit.

Age scaling for numerical stability

When fitting models with exponential or power terms (e.g., Gompertz, Makeham), large age values can cause numerical overflow or slow convergence. Centering and/or scaling the age variable helps stabilise the optimisation.

The MortalityLaw function provides the argument scale.age. When set to TRUE, ages are centered by subtracting the mean age in the fitting range and scaled by dividing by the standard deviation. This transformation often improves convergence speed and parameter identifiability without affecting the shape of the mortality curve.

# Fit with automatic age scaling
fit_scaled <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  law        = "gompertz",
  scale.age  = TRUE
)

Explanation:
- Centering removes a large constant offset, reducing correlations between the intercept and slope parameters.
- Scaling brings all ages to a comparable magnitude, which is especially beneficial for algorithms like nlm or optim that assume unit‑scale variables.
- Internally, the model is fitted on the transformed ages, and the resulting parameters are back‑transformed so that the reported values correspond to the original age scale.

You can also manually scale ages before passing them to MortalityLaw (e.g., (x - mean(x))/sd(x)), but using the built‑in scale.age argument is simpler and ensures correct back‑transformation.

Mortality laws in the package

The following table lists all parametric mortality laws currently implemented. The Code column shows the string used in the law argument of MortalityLaw.

Mortality Law Predictor Code
Gompertz \(\mu(x) = A e^{Bx}\) gompertz
Gompertz (reparam.) \(\mu(x) = \frac{1}{\sigma} \exp\left\{\frac{x-M}{\sigma}\right\}\) gompertz0
Inverse‑Gompertz \(\mu(x) = \frac{\frac{1}{\sigma}\exp\left\{\frac{x-M}{\sigma}\right\}}{\exp\left\{e^{-(x-M)/\sigma}\right\}-1}\) invgompertz
Makeham \(\mu(x) = A e^{Bx} + C\) makeham
Makeham (reparam.) \(\mu(x) = \frac{1}{\sigma} \exp\left\{\frac{x-M}{\sigma}\right\} + C\) makeham0
Opperman \(\mu(x) = \frac{A}{\sqrt{x}} + B + C \sqrt[3]{x}\) opperman
Thiele \(\mu(x) = A_1 e^{-B_1 x} + A_2 e^{-\frac12 B_2 (x-C)^2} + A_3 e^{B_3 x}\) thiele
Wittstein \(q(x) = \frac{1}{B} A^{-(B x)^N} + A^{-(M - x)^N}\) wittstein
Perks \(\mu(x) = \frac{A + B C^x}{B C^{-x} + 1 + D C^x}\) perks
Weibull \(\mu(x) = \frac{1}{\sigma} \left( \frac{x}{M} \right)^{\frac{M}{\sigma}-1}\) weibull
Inverse‑Weibull \(\mu(x) = \frac{\frac{1}{\sigma} \left( \frac{x}{M} \right)^{-\frac{M}{\sigma}-1}}{\exp\left\{\left( \frac{x}{M} \right)^{-\frac{M}{\sigma}}\right\}-1}\) invweibull
Van der Maen \(\mu(x) = A + Bx + Cx^2 + \frac{I}{N-x}\) vandermaen
Van der Maen (simpl.) \(\mu(x) = A + Bx + \frac{I}{N-x}\) vandermaen2
Quadratic \(\mu(x) = A + Bx + Cx^2\) quadratic
Beard \(\mu(x) = \frac{A e^{Bx}}{1 + K A e^{Bx}}\) beard
Makeham‑Beard \(\mu(x) = \frac{A e^{Bx}}{1 + K A e^{Bx}} + C\) makehambeard
Gamma‑Gompertz \(\mu(x) = \frac{A e^{Bx}}{1 + \frac{AG}{B}(e^{Bx}-1)}\) ggompertz
Siler \(\mu(x) = A_1 e^{-B_1 x} + A_2 + A_3 e^{B_3 x}\) siler
Heligman–Pollard \(\frac{q(x)}{p(x)} = A^{(x+B)^C} + D e^{-E (\ln x - \ln F)^2} + G H^x\) HP
Heligman–Pollard (HP2) \(q(x) = A^{(x+B)^C} + D e^{-E (\ln x - \ln F)^2} + \frac{G H^x}{1+G H^x}\) HP2
Heligman–Pollard (HP3) \(q(x) = A^{(x+B)^C} + D e^{-E (\ln x - \ln F)^2} + \frac{G H^x}{1+K G H^x}\) HP3
Rogers–Planck \(q(x) = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-U) - e^{-C(x-U)}} + A_3 e^{D x}\) rogersplanck
Martinelle \(\mu(x) = \frac{A e^{Bx} + C}{1 + D e^{Bx}} + K e^{Bx}\) martinelle
Carriere 1 \(S(x) = \psi_1 S_1(x) + \psi_2 S_2(x) + \psi_3 S_3(x)\) carriere1
Carriere 2 \(S(x) = \psi_1 S_1(x) + \psi_4 S_4(x) + \psi_3 S_3(x)\) carriere2
Kostaki \(\frac{q(x)}{p(x)} = A^{(x+B)^C} + D e^{-E_i (\ln x - \ln F)^2} + G H^x\) kostaki
Kannisto \(\mu(x) = \frac{A e^{Bx}}{1 + A e^{Bx}} + C\) kannisto

Use availableLaws() in R to see this list dynamically:

availableLaws()

Loss functions in the package

A parametric mortality model is fitted by optimising a loss function (e.g., a negative log‑likelihood or a sum of squared errors). MortalityLaws provides eight loss functions that can emphasise different portions of the mortality curve. Use availableLF() to view their descriptions:

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

Custom mortality laws

You are not limited to the built‑in laws – you can define your own hazard function and pass it via custom.law. As an example, we fit a reparametrised Gompertz model in terms of the modal age at death \(M\) (Missov et al. 2015):

\[ \mu_x = \beta e^{\beta (x - M)}. \]

First, define the hazard function. The function must return a list containing at least the hazard vector hx.

my_gompertz <- function(x, par = c(b = 0.13, M = 45)){
  hx <- with(as.list(par), b * exp(b * (x - M)))
  return(as.list(environment()))   # must return a list
}

Select data and fit:

year     <- 1950
ages     <- 45:85
deaths   <- ahmd$Dx[paste(ages), paste(year)]
exposure <- ahmd$Ex[paste(ages), paste(year)]
my_model <- MortalityLaw(
  x          = ages,
  Dx         = deaths,
  Ex         = exposure,
  custom.law = my_gompertz
)
summary(my_model)
## Custom Mortality Law
## Fitted values: mx
## 
## Call: MortalityLaw(x = ages, Dx = deaths, Ex = exposure, custom.law = my_gompertz)
## 
## Deviance Residuals:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.0029 -0.0011  0.0002  0.0009  0.0015  0.0122 
## 
## Parameters:
##       b       M 
##  0.0993 35.8014 
## 
## Residual standard error: 0.4486 on 39 degrees of freedom
plot(my_model)

Full life tables

The LifeTable function can build complete or abridged life tables from several types of input:

Only one of these inputs needs to be specified – the rest are ignored if present.

y  <- 1900
x  <- as.numeric(rownames(ahmd$mx))
Dx <- ahmd$Dx[, paste(y)]
Ex <- ahmd$Ex[, paste(y)]

LT1 <- LifeTable(x, Dx = Dx, Ex = Ex)   # primary input
LT2 <- LifeTable(x, mx = LT1$lt$mx)     # from mx
LT3 <- LifeTable(x, qx = LT1$lt$qx)     # from qx
LT4 <- LifeTable(x, lx = LT1$lt$lx)     # from lx
LT5 <- LifeTable(x, dx = LT1$lt$dx)     # from dx

LT1
## 
## Full Life Table
## 
## Number of life tables: 1 
## Dimension: 111 x 10 
## Age intervals: [0,1) [1,2) [2,3) ...   ... [108,109) [109,110) [110,+) 
## 
##      x.int   x     mx     qx   ax    lx    dx    Lx      Tx    ex
##      [0,1)   0 0.1513 0.1404 0.49 1e+05 14042 92802 4808384 48.08
##      [1,2)   1 0.0479 0.0468  0.5 85958  4023 83931 4715582 54.86
##      [2,3)   2 0.0191 0.0189  0.5 81935  1552 81157 4631651 56.53
##      [3,4)   3  0.013 0.0129  0.5 80383  1041 79861 4550494 56.61
##      [4,5)   4 0.0094 0.0093  0.5 79342   740 78971 4470633 56.35
##      [5,6)   5 0.0067 0.0067  0.5 78602   523 78340 4391661 55.87
##       <NA> ...    ...    ...  ...   ...   ...   ...     ...   ...
##  [108,109) 108   6.75 0.9988 0.15     0     0     0       0  0.15
##  [109,110) 109   6.75 0.9988 0.15     0     0     0       0  0.15
##    [110,+) 110   6.75      1 0.15     0     0     0       0  0.15
ls(LT1)   # components of the life table object
## [1] "call"         "lt"           "process_date"

Abridged life tables

Abridged life tables use wider age intervals (e.g., 0, 1–4, 5–9, …). The same LifeTable function works for any abridged age structure.

x  <- c(0, 1, seq(5, 110, by = 5))
mx <- c(.053, .005, .001, .0012, .0018, .002, .003, .004, 
       .004, .005, .006, .0093, .0129, .019, .031, .049, 
       .084, .129, .180, .2354, .3085, .390, .478, .551)

lt <- LifeTable(x, mx = mx, sex = "female")
lt
## 
## Abridged Life Table
## 
## Number of life tables: 1 
## Dimension: 24 x 10 
## Age intervals: [0,1) [1,5) [5,10) ...   ... [100,105) [105,110) [110,+) 
## 
##      x.int   x     mx     qx   ax    lx   dx     Lx      Tx    ex
##      [0,1)   0  0.053 0.0516  0.2 1e+05 5162  95878 6485467 64.85
##      [1,5)   1  0.005 0.0198 1.44 94838 1878 374547 6389590 67.37
##     [5,10)   5  0.001  0.005  2.5 92960  464 463640 6015042 64.71
##    [10,15)  10 0.0012  0.006  2.5 92496  553 461098 5551402 60.02
##    [15,20)  15 0.0018  0.009  2.5 91943  824 457653 5090304 55.36
##    [20,25)  20  0.002   0.01  2.5 91119  907 453326 4632651 50.84
##       <NA> ...    ...    ...  ...   ...  ...    ...     ...   ...
##  [100,105) 100   0.39 0.8577 1.73   407  349    896    1014  2.49
##  [105,110) 105  0.478 0.9084 1.59    58   53    110     119  2.05
##    [110,+) 110  0.551      1 1.59     5    5      8       8  1.59

Citation

citation(package = "MortalityLaws")

References

Heligman, Larry, and John H Pollard. 1980. “The Age Pattern of Mortality.” Journal of the Institute of Actuaries 107 (01): 49–80.
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.
University of California Berkeley, USA, and Max Planck Institute for Demographic Research, Germany. 2016. Human Mortality Database. https://www.mortality.org/.