library(MortalityLaws) library(knitr) opts_chunk$set(collapse = TRUE)
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.
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 = "user@email.com", # 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.
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
\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.
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()
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)
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)
# 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(package = "MortalityLaws")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.