Introduction to
MortalityLaws
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 (2016)
using the ReadHMD function:
# Download HMD data - death counts
HMD_Dx <- ReadHMD(what = "Dx",
countries = "SWE", # HMD country code for Sweden
interval = "1x1", # specify data format
username = "[email protected]", # here add your HMD username
password = "password", # here add your HMD password
save = FALSE) # save data outside RHere 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 (1980) 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)
## [1] "coefficients" "deviance" "df" "fitted.values"
## [5] "goodness.of.fit" "info" "input" "opt.diagnosis"
## [9] "residuals"A summary can be obtained using the summary
function:
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.07118 on 93 degrees of freedomThe standard plot helps us to investigate visually the goodness of 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.
Mortality laws in the package
In R one can check the availability of the
implemented models using 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()
##
## 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 valueCustom 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 (Missov et al. 2015),
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)
## 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.8021
##
## Residual standard error: 0.02097 on 39 degrees of freedom
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
##
## 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.15Abridged 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
##
## 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.59Citation
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 2.1.0,
## <https://github.com/mpascariu/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 2.1.0},
## url = {https://github.com/mpascariu/MortalityLaws},
## }