--- title: "Introduction to `MortalityLaws`" subtitle: "A guided tour of parametric mortality models and life tables" author: "Marius D. Pascariu" date: "2026-10-06" output: html_document: toc: true toc_float: collapsed: false toc_depth: 3 number_sections: false fig_width: 7 fig_height: 4.2 fig_align: center code_folding: hide bibliography: Mlaw_References.bib biblio-style: apalike link-citations: true reference-section-title: References vignette: > %\VignetteIndexEntry{Intro} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(MortalityLaws) ``` # 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 [@tabeau2001]. 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 [@heligman1980] 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 [@hmd2026]. 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](https://mpascariu.github.io/MortalityLaws/#installation) carries the instructions for both the CRAN release and the development version. The short version is `install.packages("MortalityLaws")`. ```{r} 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](Life-tables.html) 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 [@hmd2026], `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: ```{r ReadHMD, eval=FALSE} 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: ```{r RegionalReaders, eval=FALSE} 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: ```{r availableHMD, eval=FALSE} 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`. ```{r} names(HMD_sample) ``` 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 [@heligman1980], 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: ```{r} 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](Mortality-models.html) 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. ```{r} summary(fit) ``` 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. ```{r} 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. ```{r} 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: ```{r} A <- availableLaws()$table A[as.logical(A$SCALE_X), c("NAME", "CODE")] ``` 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: ```{r} 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 ``` 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: ```{r} 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))) by_hand_unscaled <- p["A"] * exp(p["B"] * ages.makeham) + p["C"] max(abs(by_hand_unscaled - fitted(fit.makeham))) ``` 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](Mortality-models.html#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$ [@missov2015]: $$\mu(x) = \beta \exp\{\beta (x - M)\}.$$ First the hazard function: ```{r} 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: ```{r} 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 ) ``` ```{r} summary(my_model) ``` 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. ```{r, warning = FALSE, message = FALSE} lt <- LawTable(x = 0:100, par = fit$coefficients, law = "HP") head(lt$lt) ``` 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](Life-tables.html) 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](Mortality-models.html) - 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](Life-tables.html) - 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](https://mpascariu.github.io/MortalityLaws/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: ```{r} citation(package = "MortalityLaws") ``` # sessionInfo() ```{r} sessionInfo() ``` The receipt for everything loaded while you were reading. If your numbers diverge from ours, this is the first thing to compare.