Introduction to `MortalityLaws`

In MortalityLaws: Parametric Mortality Models, Life Tables and HMD

library(MortalityLaws)
library(knitr)
opts_chunk$set(collapse = TRUE)

Package structure

MortalityLaws in an R package which exploits the available optimization methods to provide tools for fitting and analyzing a wide range of complex mortality models. The package can be used to construct full and abridged life tables given various input indices and to download data from Human Mortality Database as well. The main functions in the package are: MortalityLaw, LifeTable, LawTable, convertFx and ReadHMD. The package also provides support functions like availableLaws, availableLF and availableHMD that return information about the mortality laws implemented in the package, the loss functions used in optimization processes, and HMD data availability. Generic functions like predict, plot, summary, fitted, residuals are created for MortalityLaws objects. Small data set for testing purposes ahmd is saved in the package.

All functions are documented in the standard way, which means that once you load the package using library(MortalityLaws) you can just type, for example, ?MortalityLaw to see the help file.

Downloading HMD Data

Download data form Human Mortality Database [-@university_of_california_berkeley_human_2016] using the ReadHMD function:

library(MortalityLaws)

# Download HMD data - death counts
HMD_Dx <- ReadHMD(what = "Dx",
                  countries = "SWE",            # HMD country code for Sweden
                  interval  = "1x1",            # specify data format
                  username  = "[e-mail]", # here add your HMD username 
                  password  = "password",       # here add your HMD password
                  save = FALSE)                 # save data outside R

Here we download all the registered death counts Dx in Sweden from 1751 until 2014. In the same way one can download the following records: birth records births, exposure-to-risk Ex, deaths by Lexis triangles lexis, population size population, death-rates mx, life tables for females LT_f, life tables for males LT_m, life tables both sexes combined LT_t, life expectancy at birth e0, cohort death-rates mxc and cohort exposures Exc for over 40 countries and regions in different formats.

Model fitting and diagnosis

Once we have data from HMD or other sources we can start analyzing it. For example, let's fit a Heligman-Pollard [-@heligman1980] model under a Poisson setting which is already implemented as one of the standard models in the package. We have to use the MortalityLaw function in this regard.

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,   # vector with death counts
                    Ex  = exposure, # vector containing exposures
                    law = "HP",
                    opt.method = "LF2")

# inspect the output object
ls(fit) 

A summary can be obtained using the summary function:

summary(fit)

The standard plot helps us to investigate visually the goodness of fit.

plot(fit)

A model can be fitted using a subset of the data only by specifying in fit.this.x age range to be covered:

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

The gray area on the plot showing the fitted value indicates the age range used in fitting the model.

\newpage

Mortality laws in the package

\begin{table}[h] \renewcommand{\arraystretch}{1.5} \caption{Parametric functions build in the MortalityLaws package} \begin{tabular}{lll} Mortality Laws & Predictor & Code \ \hline Gompertz & $\mu(x) = Ae^{Bx}$ & gompertz \ Gompertz & $\mu(x) = \frac{1}{\sigma }exp\left{\frac{x-M}{\sigma }\right}$ & gompertz0 \ Inverse-Gompertz & $\mu(x) = {\frac{1}{\sigma }exp\left{\frac{x-M}{\sigma }\right}}/{\left(exp\left{e^{\frac{-(x-M)}{\sigma }}\right}-1\right)}$ & invgompertz\
Makeham & $\mu(x) = Ae^{Bx} + C$ & makeham \ Makeham & $\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_1e^{-B_1x} + A_2e^{-\frac{1}{2}B_2{\left(x-C\right)}^2} + A_3e^{B_3x}$ & thiele \ Wittstein \& Bumstead & q(x) = $\frac{1}{B}A^{{-(Bx)}^N} + A^{{-(M-x)}^N}$ & wittstein \ Perks & $\mu(x) = (A + BC^x) / (BC^{-x} + 1 + DC^x)$ & perks \ Weibull & $\mu(x) = \frac{1}{\sigma }{\left(\frac{x}{M}\right)}^{\frac{M}{\sigma }-1}$ & weibull \ Inverse-Weibull & $\mu(x) = {\frac{1}{\sigma }{\left(\frac{x}{M}\right)}^{-\frac{M}{\sigma }-1}}/{\left(exp\left{{\left(\frac{x}{M}\right)}^{-\frac{M}{\sigma }}\right}-1\right)}$ & invweibull \
Van der Maen & $\mu(x) = A + Bx + Cx^2 + I/(N - x)$ & vandermaen \ Van der Maen & $\mu(x) = A + Bx + I/(N - x)$ & vandermaen2 \ Quadratic & $\mu(x) = A + Bx + Cx^2$ & quadratic \ Beard & $\mu(x) = Ae^{B^x} / (1 + KAe^{B^x}) $ & beard \ Makeham-Beard & $\mu(x) = Ae^{B^x} / (1 + KAe^{B^x}) + C$ & makehambeard \ Gamma-Gompertz & $ \mu(x) = Ae^{B^x} / [1 + \frac{AG}{B} (e^{B^x} - 1)] $ & ggompertz \ Siler & $\mu(x) = A_1e^{-B_1x} + A_2 + A_3e^{B_3x}$ & siler \ Heligman - Pollard & $q(x)/p(x) = A^{{\left(x+B\right)}^C}+De^{-E{\left({\mathrm{ln} x\ }-{\mathrm{ln} F\ }\right)}^2}+GH^x$ & HP \ Heligman - Pollard & $q(x) = A^{{\left(x+B\right)}^C}+De^{-E{\left({\mathrm{ln} x\ }-{\mathrm{ln} F\ }\right)}^2} + \frac{GH^x}{1+GH^x}$ & HP2 \ Heligman - Pollard & $q(x) = A^{{\left(x+B\right)}^C}+De^{-E{\left({\mathrm{ln} x\ }-{\mathrm{ln} F\ }\right)}^2} + \frac{GH^x}{1+KGH^x}$ & HP3 \ Heligman - Pollard & $q(x) = A^{{\left(x+B\right)}^C}+De^{-E{\left({\mathrm{ln} x\ }-{\mathrm{ln} F\ }\right)}^2}+ \frac{GH^{x^K}}{1+GH^{x^K}}$ & HP3 \ Rogers-Planck & $q(x) = A_0 + A_1e^{-Ax} + A_2e^{\left{B(x - U) - e^{-C(x - U)}\right}} + A_3e^{Dx}$ & rogersplanck \ Martinelle & $\mu(x) = (Ae^{Bx} + C)/(1 + De^{Bx}) + Ke^{Bx}$ & martinelle \ Carriere & $S(x) = {\psi}1S_1(x)+\ {\psi}_2S_2\left(x\right)+\ {\psi }_3S_3\left(x\right)$ & carriere1 \ Carriere & $S(x) = {\psi}_1S_1(x)+\ {\psi}_4S_4\left(x\right)+\ {\psi }_3S_3\left(x\right)$ & carriere2 \ Kostaki & $q(x)/p(x) = A^{{(x+B)}^C}+De^{-E{i}{({\mathrm{ln} x\ }-{\mathrm{ln} F\ })}^2}+GH^x$ & kostaki \ Kannisto & $\mu(x) = Ae^{Bx} / (1+Ae^{Bx})+C$ & kannisto \ \hline \end{tabular} \end{table}

In R one can check the availability of the implemented models using availableLaws:

availableLaws()

See table 1.

Loss function in the package

A parametric model is fitted by optimizing a loss function e.g. a likelihood function or a function that minimizes errors. In MortalityLaws 8 such functions are implemented and can be used to better capture different portions of a mortality curve. Check availableLF for more details.

availableLF()

Custom mortality laws

Now let's fit a mortality law that is not defined in the package, say a reparametrize version of Gompertz in terms of modal age at death [@missov2015],

\begin{equation} \mu_x = \beta e^{\beta (x-M)}. \end{equation}

We have to define a function containing the desired hazard function and then using the custom.law argument it can be used in the MortalityLaw function.

# Here we define a function for our new model and provide start parameters
my_gompertz <- function(x, par = c(b = 0.13, M = 45)){
  hx  <- with(as.list(par), b*exp(b*(x - M)) )
  # return everything inside this function
  return(as.list(environment())) 
}

# Select data
year     <- 1950
ages     <- 45:85
deaths   <- ahmd$Dx[paste(ages), paste(year)]
exposure <- ahmd$Ex[paste(ages), paste(year)]

# Use 'custom.law' argument to instruct the MortalityLaw function how to behave
my_model <- MortalityLaw(x = ages,
                         Dx = deaths, 
                         Ex = exposure, 
                         custom.law = my_gompertz)

summary(my_model)

plot(my_model)

Full life tables

Using LifeTable function one can build full or abridged life table with various input choices like: death counts and mid-interval population estimates (Dx, Ex) or age-specific death rates (mx) or death probabilities (qx) or survivorship curve (lx) or a distribution of deaths (dx). If one of these options are specified, the other can be ignored.

# Life table for year of 1900
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)
LT2 <- LifeTable(x, mx = LT1$lt$mx)
LT3 <- LifeTable(x, qx = LT1$lt$qx)
LT4 <- LifeTable(x, lx = LT1$lt$lx)
LT5 <- LifeTable(x, dx = LT1$lt$dx)

LT1

ls(LT1)

Abridged life tables

# Example
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

Citation

citation(package = "MortalityLaws")

References

MortalityLaws documentation built on Aug. 8, 2023, 5:10 p.m.