# R/MortalityLaw_models.R In MortalityLaws: Parametric Mortality Models, Life Tables and HMD

```# -------------------------------------------------------------- #
# Author: Marius D. PASCARIU
# Last Update: Thu Jul 20 21:37:06 2023
# -------------------------------------------------------------- #

# ---- LAWS ---------------------------------------

#' Gompertz Mortality Law - 1825
#'
#' @param x vector of age at the beginning of the age classes
#' @param par parameters of the selected model. If NULL the
#' default values will be assigned automatically.
#' @examples gompertz(x = 45:90)
#' @return A list of rates and model parameters
#' @keywords internal
#' @export
gompertz <- function(x, par = NULL){
par <- bring_parameters('gompertz', par)
hx  <- with(as.list(par), A*exp(B*x) )
Hx  <- with(as.list(par), A/B * (exp(B*x) - 1) )
Sx  <- exp(-Hx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Gompertz Mortality Law - informative parameterization
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples gompertz0(x = 45:90)
#' @keywords internal
#' @export
gompertz0 <- function(x, par = NULL){
par <- bring_parameters('gompertz0', par)
hx  <- with(as.list(par), (1/sigma) * exp((x - M)/sigma) )
Hx  <- with(as.list(par), exp(-M/sigma) * (exp(x/sigma) - 1) )
Sx  <- exp(-Hx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Inverse-Gompertz Mortality Law - informative parameterization
#'
#' m - is a measure of location because it is the mode of the density, m > 0
#' sigma - represents the dispersion of the density about the mode, sigma > 0
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples invgompertz(x = 15:25)
#' @keywords internal
#' @export
invgompertz <- function(x, par = NULL){
par <- bring_parameters('invgompertz', par)
hx  <- with(as.list(par), 1/sigma * exp(-(x - M)/sigma) / (exp(exp(-(x - M)/sigma)) - 1))
Sx  <- with(as.list(par), (1 - exp(-exp(-(x - M)/sigma))) / (1 - exp(-exp(M/sigma))))
Hx  <- -log(Sx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Makeham Mortality Law - 1860
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples makeham(x = 45:90)
#' @keywords internal
#' @export
makeham <- function(x, par = NULL){
par <- bring_parameters('makeham', par)
hx  <- with(as.list(par), A*exp(B*x) + C)
Hx  <- with(as.list(par), A/B * (exp(B*x) - 1) + x*C )
Sx  <- exp(-Hx)
return(list(hx = hx, par = par))
}

#' Makeham Mortality Law - informative parameterization
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples makeham0(x = 45:90)
#' @keywords internal
#' @export
makeham0 <- function(x, par = NULL){
par <- bring_parameters('makeham0', par)
hx <- with(as.list(par), (1/sigma) * exp((x - M)/sigma) + C)
Hx <- with(as.list(par), exp(-M/sigma) * (exp(x/sigma) - 1) + x*C)
Sx <- exp(-Hx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Opperman Mortality Law - 1870
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples opperman(x = 1:25)
#' @keywords internal
#' @export
opperman <- function(x, par = NULL){
par <- bring_parameters('opperman', par)
x  <- x + 1
hx <- with(as.list(par), A/sqrt(x) - B + C*sqrt(x))
hx <- pmax(0, hx)
return(list(hx = hx, par = par))
}

#' Thiele Mortality Law - 1871
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples thiele(x = 0:100)
#' @keywords internal
#' @export
thiele <- function(x, par = NULL){
par <- bring_parameters('thiele', par)
mu1 <- with(as.list(par), A*exp(-B*x) )
mu2 <- with(as.list(par), C*exp(-.5*D*(x - E)^2) )
mu3 <- with(as.list(par), F_*exp(G*x) )
hx <- ifelse(x == 0, mu1 + mu3, mu1 + mu2 + mu3)
return(list(hx = hx, par = par))
}

#' Wittstein Mortality Law - 1883
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples wittstein(x = 0:100)
#' @keywords internal
#' @export
wittstein <- function(x, par = NULL){
par <- bring_parameters('wittstein', par)
hx  <- with(as.list(par), (1/B)*A^-((B*x)^N) + A^-((M - x)^N) )
return(list(hx = hx, par = par))
}

#' Weibull Mortality Law - 1939
#'
#' Note that if sigma > m, then the mode of the density is 0 and hx is a
#' non-increasing function of x, while if sigma < m, then the mode is
#' greater than 0 and hx is an increasing function.
#' m > 0 is a measure of location
#' sigma > 0 is measure of dispersion
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples weibull(x = 1:20)
#' @keywords internal
#' @export
weibull <- function(x, par = NULL){
par <- bring_parameters('weibull', par)
hx <- with(as.list(par), 1/sigma * (x/M)^(M/sigma - 1) )
hx[x == 0] <- 1
Hx <- with(as.list(par), (x/M)^(M/sigma) )
Sx <- exp(-Hx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Inverse-Weibull Mortality Law
#'
#' The Inverse-Weibull proves useful for modelling the childhood and teenage years,
#' because the logarithm of h(x) is a concave function.
#' m > 0 is a measure of location
#' sigma > 0 is measure of dispersion
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples invweibull(x = 1:20)
#' @keywords internal
#' @export
invweibull <- function(x, par = NULL){
par <- bring_parameters('invweibull', par)
hx <- with(as.list(par),
(1/sigma) * (x/M)^(-M/sigma - 1) / (exp((x/M)^(-M/sigma)) - 1) )
Hx <- with(as.list(par), -log(1 - exp(-(x/M)^(-M/sigma))) )
Sx <- exp(-Hx)
return(list(hx = hx, par = par, Sx = Sx))
}

#' Perks Model - 1932
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples perks(x = 50:100)
#' @keywords internal
#' @export
perks <- function(x, par = NULL){
par <- bring_parameters('perks', par)
hx  <- with(as.list(par), (A + B*C^x) / (B*(C^-x) + 1 + D*C^x) )
return(list(hx = hx, par = par))
}

#' Van der Maen Model - 1943
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples vandermaen(x = 0:100)
#' @keywords internal
#' @export
vandermaen <- function(x, par = NULL){
par <- bring_parameters('vandermaen', par)
hx  <- with(as.list(par), A + B*x + C*(x^2) + I/(N - x))
return(list(hx = hx, par = par))
}

#' Van der Maen 2 Model - 1943
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples vandermaen(x = 0:100)
#' @keywords internal
#' @export
vandermaen2 <- function(x, par = NULL){
par <- bring_parameters('vandermaen2', par)
hx  <- with(as.list(par), A + B*x + I/(N - x))
return(list(hx = hx, par = par))
}

#' Strehler-Mildvan Model - 1960
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples strehler_mildvan(x = 30:85)
#' @keywords internal
#' @export
strehler_mildvan <- function(x, par = NULL){
par <- bring_parameters('strehler_mildvan', par)
hx  <- with(as.list(par), K * exp(-V0 * (1 - B * x)/D) )
return(list(hx = hx, par = par))
}

#' Beard Model - 1971
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples beard(x = 50:100)
#' @keywords internal
#' @export
beard <- function(x, par = NULL){
par <- bring_parameters('beard', par)
hx  <- with(as.list(par), (A*exp(B*x)) / (1 + K*A*exp(B*x)) )
return(list(hx = hx, par = par))
}

#' Makeham-Beard Model - 1971
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples beard_makeham(x = 0:100)
#' @keywords internal
#' @export
beard_makeham <- function(x, par = NULL){
par <- bring_parameters('beard_makeham', par)
hx  <- with(as.list(par), A*exp(B*x) / (1 + K*A*exp(B*x)) + C)
return(list(hx = hx, par = par))
}

#' Gamma-Gompertz Model as in Vaupel et al. (1979)
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples ggompertz(x = 50:120)
#' @keywords internal
#' @export
ggompertz <- function(x, par = NULL){
par <- bring_parameters('ggompertz', par)
hx  <- with(as.list(par), (A*exp(B*x)) / (1 + (A*G/B)*(exp(B*x) - 1)) )
return(list(hx = hx, par = par))
}

#' @inheritParams gompertz
#' @inherit gompertz return
#' @keywords internal
#' @export
quadratic <- function(x, par = NULL){
hx  <- with(as.list(par), A + B*x + C*(x^2))
return(list(hx = hx, par = par))
}

#' Siler Mortality Law - 1979
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples siler(x = 0:100)
#' @keywords internal
#' @export
siler <- function(x, par = NULL){
par <- bring_parameters('siler', par)
hx <- with(as.list(par), A*exp(-B*x) + C + D*exp(E*x))
return(list(hx = hx, par = par))
}

#' Heligman-Pollard Mortality Law - 8 parameters - 1980
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples HP(x = 0:100)
#' @keywords internal
#' @export
HP <- function(x, par = NULL){
par <- bring_parameters('HP', par)
mu1 <- with(as.list(par), A^((x + B)^C) + G*H^x )
mu2 <- with(as.list(par), D*exp(-E*(log(x/F_))^2) )
eta <- ifelse(x == 0, mu1, mu1 + mu2)
hx <- eta/(1 + eta)
return(list(hx = hx, par = par))
}

#' Heligman-Pollard 2 Mortality Law - 8 parameters
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples HP2(x = 0:100)
#' @keywords internal
#' @export
HP2 <- function(x, par = NULL){
par <- bring_parameters('HP2', par)
mu1 <- with(as.list(par), A^((x + B)^C) + (G*H^x)/(1 + G*H^x) )
mu2 <- with(as.list(par), D*exp(-E*(log(x/F_))^2) )
eta <- ifelse(x == 0, mu1, mu1 + mu2)
hx <- eta
return(list(hx = hx, par = par))
}

#' Heligman-Pollard 3 Mortality Law - 9 parameters
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples HP3(x = 0:100)
#' @keywords internal
#' @export
HP3 <- function(x, par = NULL){
par <- bring_parameters('HP3', par)
mu1 <- with(as.list(par), A^((x + B)^C) + (G*H^x)/(1 + K*G*H^x) )
mu2 <- with(as.list(par), D*exp(-E*(log(x/F_))^2) )
eta <- ifelse(x == 0, mu1, mu1 + mu2)
hx <- eta
return(list(hx = hx, par = par))
}

#' Heligman-Pollard 4 Mortality Law - 9 parameters
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples HP4(x = 0:100)
#' @keywords internal
#' @export
HP4 <- function(x, par = NULL){
par <- bring_parameters('HP4', par)
mu1 <- with(as.list(par), A^((x + B)^C) + (G*H^(x^K)) / (1 + G*H^(x^K)) )
mu2 <- with(as.list(par), D*exp(-E*(log(x/F_))^2) )
eta <- ifelse(x == 0, mu1, mu1 + mu2)
hx <- eta
return(list(hx = hx, par = par))
}

#' Martinelle Model - 1987
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples martinelle(x = 0:100)
#' @keywords internal
#' @export
martinelle <- function(x, par = NULL){
par <- bring_parameters('martinelle', par)
hx  <- with(as.list(par), (A*exp(B*x) + C) / (1 + D*exp(B*x)) + K*exp(B*x))
return(list(hx = hx, par = par))
}

#' Rogers-Planck Model - 1983
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples rogersplanck(x = 0:100)
#' @keywords internal
#' @export
rogersplanck <- function(x, par = NULL){
par <- bring_parameters('rogersplanck', par)
hx  <- with(as.list(par),
A0 + A1*exp(-A*x) + A2*exp(B*(x - U) - exp(-C*(x - U))) + A3*exp(D*x))
return(list(hx = hx, par = par))
}

#' Carriere Mortality Law - 1992
#'
#' Carriere1 = weibull + invweibull + gompertz
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples carriere1(x = 0:100)
#' @keywords internal
#' @export
carriere1 <- function(x, par = NULL){
par <- bring_parameters('carriere1', par)
# Compute distribution functions
S_wei  <- weibull(x, par[c('sigma1', 'M1')])\$Sx
S_iwei <- invweibull(x, par[c('sigma1', 'M2')])\$Sx
S_gom  <- gompertz0(x, par[c('sigma3', 'M3')])\$Sx

f1 <- par['P1'] <- max(0.0001, min(par['P1'], 1))
f2 <- par['P2'] <- max(0.0001, min(par['P2'], 1))
f3 <- 1 - f1 - f2

Sx <- f1*S_wei + f2*S_iwei + f3*S_gom
Sx <- pmax(0, pmin(1, Sx))
Hx <- -log(Sx)
hx <- c(Hx[1], diff(Hx)) # here we will need a numerical solution!
return(list(hx = hx, par = par))
}

#' Carriere Mortality Law - 1992
#'
#' Carriere2 = weibull + invgompertz + gompertz
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples carriere2(x = 0:100)
#' @keywords internal
#' @export
carriere2 <- function(x, par = NULL){
par <- bring_parameters('carriere2', par)
# Compute distribution functions
S_wei  <- weibull(x, par[c('sigma1', 'M1')])\$Sx
S_igom <- invgompertz(x, par[c('sigma2', 'M2')])\$Sx
S_gom  <- gompertz0(x, par[c('sigma3', 'M3')])\$Sx

f1 <- par['P1'] <- max(0.0001, min(par['P1'], 1))
f2 <- par['P2'] <- max(0.0001, min(par['P2'], 1))
f3 <- 1 - f1 - f2

Sx <- f1*S_wei + f2*S_igom + f3*S_gom
Sx <- pmax(0, pmin(1, Sx))
Hx <- -log(Sx)
hx <- c(Hx[1], diff(Hx)) # here we will need a numerical solution!
return(list(hx = hx, par = par))
}

#' Kostaki Model - 1992
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples kostaki(x = 0:100)
#' @keywords internal
#' @export
kostaki <- function(x, par = NULL){
par <- bring_parameters('kostaki', par)
with(as.list(par), {
# Sometimes the difference between estimated parameters E1 and E2 is
# very large, in which case the resulted mortality curve will exhibit
# a significant artificial jump in one age grup. I am imposing a
# restriction below to limit this behaviour.
if (E1 >= 50*E2) E2 <- E1/50 # This hack seems to work.

L   <- x <= F_    # Logical
mu1 <- A^((x + B)^C) + G*H^x
e1  <- -(E1*log(x/F_))^2
e2  <- -(E2*log(x/F_))^2
mu2 <- D*exp(L*e1 + (!L)*e2)
eta <- ifelse(x == 0, mu1, mu1 + mu2)
hx  <- eta/(1 + eta)
return(list(hx = hx, par = par))
})
}

#' Kannisto Mortality Law - 1998
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples kannisto(x = 85:120)
#' @keywords internal
#' @export
kannisto <- function(x, par = NULL){
par <- bring_parameters('kannisto', par)
with(as.list(par), {
hx  <- A*exp(B*x) / (1 + A*exp(B*x))
Hx  <- 1/A * log((1 + B*exp(B*x)) / (1 + A))
Sx  <- exp(-Hx)
return(list(hx = hx, par = par))
})
}

#' Kannisto-Makeham Mortality Law - 1998
#' @inheritParams gompertz
#' @inherit gompertz return
#' @examples kannisto_makeham(x = 85:120)
#' @keywords internal
#' @export
kannisto_makeham <- function(x, par = NULL){
par <- bring_parameters('kannisto_makeham', par)
with(as.list(par), {
hx  <- A*exp(B*x) / (1 + A*exp(B*x)) + C
return(list(hx = hx, par = par))
})
}

#' Bring or Rename Starting Parameters in the Law Functions
#' @inheritParams MortalityLaw
#' @inheritParams gompertz
#' @return Vector or initial model parameters
#' @keywords internal
bring_parameters <- function(law, par = NULL) {
Spar <- switch(law,
demoivre    = c(A = 105),
gompertz    = c(A = 0.0002, B = 0.13),
gompertz0   = c(sigma = 7.7, M = 49),
invgompertz = c(sigma = 7.7, M = 49),
makeham     = c(A = .0002, B = .13, C = .001),
makeham0    = c(sigma = 7.692308, M = 49, C = .001),
opperman    = c(A = .04, B = .0004, C = .001),
thiele      = c(A = .02474, B = .3, C = .004, D = .5,
E = 25, F_ = .0001, G = .13),
wittstein   = c(A = 1.5, B = 1, N = .5, M = 100),
perks       = c(A = .002, B = .13, C = .01, D = .01),
weibull     = c(sigma = 2, M = 1),
invweibull  = c(sigma = 10, M = 5),
vandermaen  = c(A = .01, B = 1, C = .01, I = 100, N = 200),
vandermaen2 = c(A = .01, B = 1, I = 100, N = 200),
strehler_mildvan = c(K = .01, V0 = 2.5, B = 0.2, D = 6),
quadratic   = c(A = .01, B = 1, C = .01),
beard       = c(A = .002, B = .13, K = 1),
beard_makeham = c(A = .002, B = .13, C = .01, K = 1),
ggompertz = c(A = .002, B = .13, G = 1),
siler     = c(A = .0002, B = .13, C = .001, D = .001, E = .013),
HP        = c(A = .0005, B = .004, C = .08, D = .001,
E = 10, F_ = 17, G = .00005, H = 1.1),
HP2       = c(A = .0005, B = .004, C = .08, D = .001,
E = 10, F_ = 17, G = .00005, H = 1.1),
HP3       = c(A = .0005, B = .004, C = .08, D = .001,
E = 10, F_ = 17, G = .00005, H = 1.1, K = 1),
HP4       = c(A = .0005, B = .004, C = .08, D = .001,
E = 10, F_ = 17, G = .00005, H = 1.1, K = 1),
rogersplanck = c(A0 = .0001, A1 = .02, A2 = .001, A3 = .0001,
A = 2, B = .001, C = 100, D = .1, U = .33),
martinelle = c(A = .001, B = .13, C = .001, D = 0.1, K = .001),
kostaki    = c(A = .0005, B = .01, C = .10, D = .001,
E1 = 3, E2 = .1, F_ = 25, G = .00005, H = 1.1),
carriere1  = c(P1 = .003, sigma1 = 15, M1 = 2.7,
P2 = .007, sigma2 = 6, M2 = 3,
sigma3 = 9.5, M3 = 88),
carriere2  = c(P1 = .01, sigma1 = 2, M1 = 1,
P2 = .01, sigma2 = 7, M2 = 49,
sigma3 = 7, M3 = 49),
kannisto   = c(A = 0.5, B = 0.13),
kannisto_makeham = c(A = 0.5, B = 0.13, C = 0.001)
)
if (is.null(par)) par <- Spar
# If 'par' is provided, just give them a name anyway.
names(par) <- names(Spar)
return(par)
}
```

## Try the MortalityLaws package in your browser

Any scripts or data that you put into this service are public.

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