Welcome

MortalityLaws is a small package for a classic demography problem: you have a mortality schedule, noisy and full of warts, and you would like a smooth formula that describes it. Parametric mortality laws are compact summaries of the age pattern of mortality. A good one lets you compare populations on equal footing, graduate a series with gaps in it, describe a life in a handful of numbers, and project the curve forward, which is where most mortality forecasting begins (Tabeau et al. 2001).

This vignette is the tour. By the end of it you will have taken a real mortality schedule, the England & Wales females of 1950 that ship inside the package, and done the full round with it:

  • fitted the Heligman-Pollard law (Heligman and Pollard 1980) and read its diagnostics;
  • refitted the same law on a chosen window of ages;
  • written a mortality law of your own and fitted that too;
  • turned the fitted curve into a life table;
  • and fetched fresh data from the Human Mortality Database (HMD 2026).

Everything you need is bundled, so nothing here requires an account or a download. Let us start.

Install and load

The Installation section of the README carries the instructions for both the CRAN release and the development version. The short version is install.packages("MortalityLaws").

library(MortalityLaws)
data(ahmd)

ahmd is one of the small datasets bundled with the package: deaths, exposures and death rates for England & Wales females, ages 0 to 110, in the four census years 1850, 1900, 1950 and 2010. It is the workhorse of the examples below and of the help pages.

Mortality in four shapes

Mortality has four classic shapes: the death rate \(m(x)\), the death probability \(q(x)\), the survivorship \(l(x)\) and the death distribution \(d(x)\). Real data arrive in two further wrappings: raw death counts with their exposures, and life expectancy \(e(x)\), the summary everyone quotes. The package accepts all six as input, exactly one at a time.

Case What it holds
Dx + Ex Death counts and exposure-to-risk; their ratio is the observed death rate at age \(x\).
mx Death rates: the force of mortality, the per-year hazard of dying at age \(x\).
qx Death probabilities: the chance of dying between age \(x\) and the next birthday.
lx Survivorship: how many of a synthetic cohort of \(l_0\) newborns reach age \(x\).
dx Cohort deaths: how many of those \(l_0\) die between age \(x\) and the next birthday.
ex Life expectancy at age \(x\): the average number of years still to live.

MortalityLaw() fits on the first three cases: deaths and exposures, rates, or probabilities. LifeTable() builds a life table from any of the six, and convertFx() translates between them, rates into probabilities, survivors into rates, and so on. Think of the six cases as interchangeable currencies for the same thing; the Life tables article develops the notation, the recursions and the exchange rules at leisure.

Getting data

There are two ways to get a mortality schedule into R: download it, or use what ships with the package. Both are short.

Four national databases have a reader: ReadHMD() for the Human Mortality Database (HMD 2026), ReadJMD() for the Japanese Mortality Database, ReadCHMD() for the Canadian Human Mortality Database, and ReadAHMD() for the Australian Human Mortality Database. The HMD wants a free account; the other three are open.

A first example, Swedish death counts at single years of age:

HMD_Dx <- ReadHMD(
  what      = "Dx",
  countries = "SWE",
  interval  = "1x1",
  username  = "user@email.com",
  password  = "password",
  save      = FALSE
)

The what argument names the product: death counts (Dx), exposures (Ex) and death rates (mx); births and population; deaths and exposures split by Lexis triangle (Dx_lexis, Ex_lexis); period life tables by sex (LT_f, LT_m, LT_t) and their cohort versions; cohort rates (mxc, Exc); and life expectancy at birth (e0, e0c). The interval argument sets the format: 1x1, 1x5, 1x10, 5x1, 5x5 or 5x10. Together these products cover 50 countries and regions, so the HMD really is the reference collection of human mortality records.

The regional readers work the same way with regions in place of countries, and without a login:

JMD_LT <- ReadJMD(       # Japanese prefectures: female life tables
  what     = "LT_f",
  regions  = c("Aichi", "Tokyo"),
  interval = "1x1",
  save     = FALSE
)

CHMD_mx <- ReadCHMD(     # Canadian regions: death rates
  what     = "mx",
  regions  = "CAN",
  interval = "1x1",
  save     = FALSE
)

AHMD_Ex <- ReadAHMD(     # Australian states: exposures
  what     = "Ex",
  regions  = c("NSW", "VIC"),
  interval = "1x1",
  save     = FALSE
)

If a database does not publish what you asked for in the format you asked for, the reader tells you so and skips that piece; it never takes the whole call down with it. And if you are not sure what is on offer to begin with, availableHMD() scrapes the HMD data-availability table and lays it out:

availableHMD()

Not ready to spend a login? Each reader ships with a sample of its output, so you can write and test your code first and download later. One real download is bundled as HMD_sample, JMD_sample, CHMD_sample and AHMD_sample.

names(HMD_sample)
#> [1] "input"         "data"          "download.date" "years"        
#> [5] "ages"

Five pieces: the input arguments, the tidy data with one row per age and year, the download.date, and the years and ages covered. The Read* help pages carry the same detail for the other three readers.

Fitting a law

A parametric mortality law is a formula \(f(x; \theta)\) with a handful of parameters \(\theta\), shaped so that it can only produce plausible mortality curves. Fitting is the search for the \(\theta\) whose curve sits closest to the data.

We fit the Heligman-Pollard law (Heligman and Pollard 1980), a compact formula whose three terms trace the three famous movements of the mortality curve: the fall of infant and child mortality, the accident hump of young adulthood, and the exponential rise of adult mortality. It was built for exactly the kind of data in ahmd, and we take the England & Wales females of 1950:

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"
)

law = "HP" selects the model from the catalogue. opt.method = "LF2" picks the loss the optimiser minimises, here the squared log-ratio between fitted and observed values, which weighs the young ages as heavily as the old. The package offers eight such objectives (availableLF()), and the Mortality models article walks through the catalogue of laws and the estimation engine behind them.

What comes back

The fit is an object of class "MortalityLaw", and every piece of the estimation is kept in it.

Element What it holds
input The arguments as supplied: data, ages, law code, loss, fit window.
info The law’s metadata: name, formula, type, and the date of the fit.
coefficients The estimated parameters \(\hat\theta\), named.
fitted.values The fitted curve, one value per age, on the original age scale.
residuals Observed minus fitted, on the scale of the data.
deviance.residuals Deviance residuals, one per age.
pearson.residuals Pearson residuals: \((D_x - \hat\mu_x E_x) / \sqrt{\hat\mu_x E_x}\) for counts.
goodness.of.fit Log-likelihood, AIC and BIC, filled in for the likelihood losses only.
opt.diagnosis The optimiser’s report: convergence code, iterations, objective value.
df Degrees of freedom: parameters estimated, and residuals left over.
dispersion The dispersion statistic of the fit.
deviance The deviance at the optimum.

coef(), fitted() and residuals() pull out the pieces you use most, and summary() prints the report a modeller wants first: what was fitted on what, the parameter estimates, and the goodness-of-fit numbers.

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  |  ages 0-100  |  fitted on 0-100 (101 of 101 ages)
#> 
#> Call:
#> MortalityLaw(x = ages, Dx = deaths, Ex = exposure, law = "HP", 
#>     opt.method = "LF2")
#> 
#> Coefficients:
#>    estimate
#> A    0.0022
#> B    0.0146
#> C    0.1229
#> D    0.0009
#> E    2.7565
#> F_  29.0080
#> G    0.0000
#> H    1.1141
#> 
#> Fit:
#>   method LF2  |  optimiser converged in 23 iterations
#>   deviance 721.3 on 93 degrees of freedom  |  dispersion 7.719
#>   R-squared 0.9967  |  RMSE 0.007083
#> 
#> Residuals:
#>             Min. 1st Qu. Median    Mean 3rd Qu.   Max.
#> raw      -0.0312 -0.0002 0.0000  0.0009  0.0003 0.0462
#> deviance -8.7216 -1.5237 0.0198 -0.1061  1.2903 6.2744

The summary is quiet about likelihoods here, because there are none to report: LF2 is not a likelihood. Ask for opt.method = "poissonL" or "binomialL" and log-likelihood, AIC and BIC appear in their place.

Reading the figure

plot(fit) draws the whole diagnosis in one figure: the fit chart on top, and four residual panels underneath.

plot(fit)

The top panel is the one you can show anyone: observed points against the fitted curve on a logarithmic scale, the ages that carried the fit shaded, and \(R^2\) and RMSE in the subtitle. The four panels below are for the modeller. Deviance residuals against age, and against fitted values, show where and at what level the law misses. A normal Q-Q plot and the residual distribution show whether the misses look like noise or like structure. The first panel is the honest one: a Heligman-Pollard curve cannot track every shoulder in the data, and the residuals against age tell you exactly which shoulders those were.

which selects the figure: "both" (the default), "fit" for the chart alone, or "diagnostics" for the four panels alone. split rearranges those panels, c(1, 4) for one row of four.

Fitting on a subset of ages

Sometimes the law is only wanted for part of the lifespan, or only part of the data can be trusted. fit.this.x names the ages that take part in the optimisation; the fit is still evaluated and drawn over every age you passed in.

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

The shaded band marks the ages that carried the fit. Beyond it the curve is pure extrapolation: the model’s opinion of what mortality does next, informed by whatever the parameters learned inside the window. Used as such, it is one of the most useful moves in the package. Read as data, it is a trap.

Age scaling

Some laws are built for one stretch of the lifespan, and their exponential or power terms fall apart numerically when evaluated far outside it. A term like \(A e^{Bx}\) is at home near its own origin: the further an age strays from that origin, the more wildly a small change in \(B\) swings the value, and the optimiser ends up fighting arithmetic instead of shape. Such laws carry SCALE_X = TRUE in availableLaws(), and MortalityLaw() rescales the age vector before fitting:

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

so that the youngest age in the fitting range becomes 1. No user argument is involved: the flag travels with the law definition.

Two consequences are worth keeping in mind. The parameters refer to the scaled ages and are not transformed back, so they are not directly comparable with the parameterisation of the unscaled law: the same curve written on the original axis carries different coefficients. And the fitted and predicted curves are always returned on the original age scale, because that is the scale your data live on: predict() and LawTable() apply the same shift consistently. When you reuse scaled coefficients in LawTable(), start the table at the same lower age you fitted on.

Which laws are fitted on a rescaled age vector? This lists them:

A <- availableLaws()$table
A[as.logical(A$SCALE_X), c("NAME", "CODE")]
#>                NAME             CODE
#> 2          Gompertz         gompertz
#> 3          Gompertz        gompertz0
#> 4  Inverse-Gompertz      invgompertz
#> 5           Makeham          makeham
#> 6           Makeham         makeham0
#> 11       Steffensen       steffensen
#> 12            Perks            perks
#> 15  Inverse-Weibull       invweibull
#> 16     Van der Maen       vandermaen
#> 17     Van der Maen      vandermaen2
#> 18 Strehler-Mildvan strehler_mildvan
#> 19        Quadratic        quadratic
#> 20            Beard            beard
#> 21    Beard-Makeham    beard_makeham
#> 22   Gamma-Gompertz        ggompertz
#> 30 Gompertz-Makeham  makeham_logquad
#> 31 Gompertz-Makeham gompertz_logquad
#> 32         Carriere        carriere1
#> 33         Carriere        carriere2
#> 35         Kannisto         kannisto
#> 36 Kannisto-Makeham kannisto_makeham

Here is what that means on one of the laws from the list. Fit makeham, \(\mu(x) = A e^{Bx} + C\), to the England & Wales females of 2010 over ages 40 to 90:

ages.makeham <- 40:90
fit.makeham  <- MortalityLaw(
  x          = ages.makeham,
  Dx         = ahmd$Dx[paste(ages.makeham), "2010"],
  Ex         = ahmd$Ex[paste(ages.makeham), "2010"],
  law        = "makeham",
  opt.method = "LF2"
)
p <- coef(fit.makeham)
p
#>            A            B            C 
#> 0.0004571225 0.1109898211 0.0004955652

Three coefficients come back, and they belong to the shifted axis. The ages 40 to 90 reached the law as 1 to 51, so the fitted curve is the printed \(A e^{B x_{\text{fit}}} + C\); written on the original ages the same curve carries the leading coefficient \(A e^{-39B} \approx 6.03 \times 10^{-6}\), smaller than the printed \(A\) by a factor of \(e^{39B} \approx 75.8\). Same curve, different bookkeeping.

Now evaluate the law by hand from p, twice. Once with the shift the engine used, once on the raw ages, which is what the formula appears to ask for:

x_scaled <- ages.makeham - min(ages.makeham) + 1
by_hand  <- p["A"] * exp(p["B"] * x_scaled) + p["C"]
max(abs(by_hand - fitted(fit.makeham)))
#> [1] 0

by_hand_unscaled <- p["A"] * exp(p["B"] * ages.makeham) + p["C"]
max(abs(by_hand_unscaled - fitted(fit.makeham)))
#> [1] 9.828153

The first comparison prints zero: on the scaled axis the hand computation reproduces fitted() exactly. The second is off by up to 9.83 in absolute hazard, and at age 90 it returns 9.96 where the fit says 0.1318, about 76 times too high. Nothing is wrong with the fit; the arithmetic was simply run 39 years further along the exponential than the coefficients were estimated on.

So the recipe for a stored fit is to keep the offset in mind. Shift by the same amount before plugging any age into the formula, or let predict() do it for you: predict(fit.makeham, x = 95) returns 0.2292, the same value the hand computation gives at the shifted age 56. When you rebuild a table with LawTable() from scaled coefficients, start it at the lower age you fitted on, so that the engine’s shift lands where the coefficients expect it. A makeham table run over 0:100 would hand the engine a different origin and quietly come back as the fitted curve shifted 40 years, with age 0 carrying the mortality of age 40.

The same trap is listed in the How to get it wrong section of the Mortality models article.

Custom laws

The catalogue is not the limit. custom.law takes your own mortality law and makes it first class in the fitting engine. The function needs the signature function(x, par), must return a list containing at least the hazard vector hx, and offers its starting parameter values in its own defaults; the engine reads them by calling custom.law(1)$par.

Our example is a Gompertz law written in terms of the modal age at death \(M\) (Missov et al. 2015):

\[\mu(x) = \beta \exp\{\beta (x - M)\}.\]

First the hazard function:

missov <- 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
}

Then the data and the fit, on ages 45 to 85 of the same England & Wales 1950 schedule:

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 = missov
)
summary(my_model)
#> Custom Mortality Law
#> Fitted values: mx  |  ages 45-85  |  fitted on 45-85 (41 of 41 ages)
#> 
#> Call:
#> MortalityLaw(x = ages, Dx = deaths, Ex = exposure, custom.law = missov)
#> 
#> Coefficients:
#>   estimate
#> b   0.0993
#> M  35.8014
#> 
#> Fit:
#>   method LF2  |  optimiser converged in 11 iterations
#>   deviance 755 on 39 degrees of freedom  |  dispersion 19.36
#>   R-squared 0.9951  |  RMSE 0.003272
#> 
#> Residuals:
#>             Min. 1st Qu. Median    Mean 3rd Qu.   Max.
#> raw      -0.0029 -0.0011 0.0002  0.0009  0.0015 0.0122
#> deviance -8.5379 -3.8161 0.1859 -0.0687  2.9870 8.0085

A custom law is treated as SCALE_X = TRUE (see the previous section), so the fit ran on ages 1 to 41 and the reported \(M = 35.8\) counts from there. Back on the original scale the modal age at death is \(35.8 + 45 - 1 \approx 79.8\) years, and that is the number a demographer would recognise. Everything else about the object, the diagnostics and the plots, works exactly as it does for a catalogued law.

From fit to life table

A fitted law is a complete description of mortality, so a life table follows mechanically from it. LawTable() does the round in one call: it evaluates the law at the ages you give it and hands the result to the life table engine.

lt <- LawTable(x = 0:100, par = fit$coefficients, law = "HP")
head(lt$lt)
#>   x.int x           mx           qx        ax        lx         dx       Lx
#> 1 [0,1) 0 0.0263700569 0.0257797415 0.1316506 100000.00 2577.97415 97761.42
#> 2 [1,2) 1 0.0022229673 0.0022204992 0.5000000  97422.03  216.32553 97313.86
#> 3 [2,3) 2 0.0013104855 0.0013096274 0.5000000  97205.70  127.30325 97142.05
#> 4 [3,4) 3 0.0009449708 0.0009445245 0.5000000  97078.40   91.69292 97032.55
#> 5 [4,5) 4 0.0007450915 0.0007448141 0.5000000  96986.70   72.23706 96950.59
#> 6 [5,6) 5 0.0006191301 0.0006189385 0.5000000  96914.47   59.98410 96884.48
#>        Tx       ex
#> 1 7093709 70.93709
#> 2 6995948 71.81074
#> 3 6898634 70.96944
#> 4 6801492 70.06184
#> 5 6704459 69.12761
#> 6 6607508 68.17876

The table has one row per age and the usual columns: the age, the mortality schedule in its interchangeable forms, the survivors lx, the deaths dx, the person-years Lx and Tx, and life expectancy ex. A model curve never reaches a death probability of exactly 1 at the last age in the range, so the call notes that the table is not closed at the top; the engine closes the final interval for you and says so. Closing the input yourself is optional. (The note is switched off above to keep the output readable.)

There are no life table figures here on purpose. The Life tables article takes this table apart: the recursions that build every column, the average-years-lived terms, what to do about the open tail, indicator conversions with convertFx(), and the inverse problem of building a life table from given life expectancies.

Where to go next

This tour covered the round trip: data, law, fit, diagnosis, life table. Two companion articles go deep where this one stayed shallow.

  • Mortality models and estimation - the catalogue of every law with its formula and citation, the fitting engine in detail (the eight loss functions, the optimisers, starting values, age rescaling), and how to read every diagnostic.
  • Life tables - the life table engine: notation and recursions, the average-years-lived methods, closing and extending the tail, indicator conversions with convertFx(), LawTable(), and the inverse life table from ex.
  • The function reference - every function and dataset, with the arguments spelled out.

If the package earns a line in your paper, the citation is one call away:

citation(package = "MortalityLaws")
#> Warning in citation(package = "MortalityLaws"): could not determine year for
#> 'MortalityLaws' from package DESCRIPTION file
#> To cite package 'MortalityLaws' in publications use:
#> 
#>   Pascariu M (????). _MortalityLaws: Parametric Mortality Models, Life
#>   Tables and HMD_. R package version 3.0.0,
#>   <https://mpascariu.github.io/MortalityLaws/>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Manual{,
#>     title = {MortalityLaws: Parametric Mortality Models, Life Tables and HMD},
#>     author = {Marius D. Pascariu},
#>     note = {R package version 3.0.0},
#>     url = {https://mpascariu.github.io/MortalityLaws/},
#>   }

sessionInfo()

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

The receipt for everything loaded while you were reading. If your numbers diverge from ours, this is the first thing to compare.

References

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/.
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.
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.