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

#### Documented in addDetailschoose_optimMortalityLawobjective_funscale_x

# --------------------------------------------------- #
# Author: Marius D. PASCARIU
# Last update: Wed Nov 24 11:19:28 2021
# --------------------------------------------------- #

#' Fit Mortality Laws
#'
#' Fit parametric mortality models given a set of input data which can be
#' represented by death counts and mid-interval population estimates
#' \code{(Dx, Ex)} or age-specific death rates \code{(mx)} or death
#' probabilities \code{(qx)}. Using the argument \code{law} one can specify
#' the model to be fitted. So far more than 27 parametric models have been
#' implemented; check the \code{\link{availableLaws}} function to learn
#' about the available options. The models can be fitted under
#' the maximum likelihood methodology or by selecting a loss function to be
#' optimised. See the implemented loss function by running the
#' @usage
#' MortalityLaw(x, Dx = NULL, Ex = NULL, mx = NULL, qx = NULL,
#'                 law = NULL,
#'                 opt.method = "LF2",
#'                 parS = NULL,
#'                 fit.this.x = x,
#'                 custom.law = NULL,
#'                 show = FALSE, ...)
#' @details Depending on the complexity of the model, one of following
#' optimization strategies is employed: \enumerate{
#'  \item{Nelder-Mead method:}{ approximates a local optimum of a problem with n
#'   variables when the objective function varies smoothly and is unimodal.
#'   For details see \code{\link{optim}}}
#'  \item{PORT routines:}{ provides unconstrained optimization and optimization
#'  subject to box constraints for complicated functions. For details check
#'  \item{Levenberg-Marquardt algorithm:}{ damped least-squares method.
#'  For details check \code{\link{nls.lm}}}
#' }
#' @inheritParams LifeTable
#' @param law The name of the mortality law/model to be used. e.g.
#' \code{gompertz}, \code{makeham}, ... To investigate all the possible options,
#' see \code{\link{availableLaws}} function.
#' @param opt.method How would you like to find the parameters? Specify the
#' function to be optimize. Available options: the Poisson likelihood function
#' \code{poissonL}; the Binomial likelihood function -\code{binomialL}; and
#' 6 other loss functions. For more details, check the \code{\link{availableLF}}
#' function.
#' @param parS Starting parameters used in the optimization process (optional).
#' @param fit.this.x Select the ages to be considered in model fitting.
#' By default \code{fit.this.x = x}. One may want to exclude from the fitting
#' procedure, say, the advanced ages where the data is sparse.
#' @param custom.law Allows you to fit a model that is not defined
#' in the package. Accepts as input a function.
#' @param show Choose whether to display a progress bar during the fitting
#' process. Logical. Default: \code{FALSE}.
#' @param ... Arguments to be passed to or from other methods.
#' @return The output is of the \code{"MortalityLaw"} class with the components:
#'  \item{input}{List with arguments provided in input. Saved for convenience.}
#'  \item{info}{Brief information about the model.}
#'  \item{coefficients}{Estimated coefficients.}
#'  \item{fitted.values}{Fitted values of the selected model.}
#'  \item{residuals}{Deviance residuals.}
#'  \item{goodness.of.fit}{List containing goodness of fit measures like
#' AIC, BIC and log-Likelihood.}
#'  \item{opt.diagnosis}{Resultant optimization object useful for
#' checking the convergence etc.}
#' @seealso
#' @author Marius D. Pascariu
#' @examples
#' # Example 1: --------------------------
#' # Fit Makeham Model for Year of 1950.
#'
#' x  <- 45:75
#' Dx <- ahmd$Dx[paste(x), "1950"] #' Ex <- ahmd$Ex[paste(x), "1950"]
#'
#' M1 <- MortalityLaw(x = x, Dx = Dx, Ex = Ex, law = 'makeham')
#'
#' M1
#' ls(M1)
#' coef(M1)
#' summary(M1)
#' fitted(M1)
#' predict(M1, x = 45:95)
#' plot(M1)
#'
#'
#' # Example 2: --------------------------
#' # We can fit the same model using a different data format
#' # and a different optimization method.
#' x  <- 45:75
#' mx <- ahmd$mx[paste(x), ] #' M2 <- MortalityLaw(x = x, mx = mx, law = 'makeham', opt.method = 'LF1') #' M2 #' fitted(M2) #' predict(M2, x = 55:90) #' #' # Example 3: -------------------------- #' # Now let's fit a mortality law that is not defined #' # in the package, say a reparameterized Gompertz in #' # terms of modal age at death #' # hx = b*exp(b*(x-m)) (here b and m are the parameters to be estimated) #' #' # A function with 'x' and 'par' as input has to be defined, which returns #' # at least an object called 'hx' (hazard rate). #' my_gompertz <- function(x, par = c(b = 0.13, M = 45)){ #' hx <- with(as.list(par), b*exp(b*(x - M)) ) #' return(as.list(environment())) #' } #' #' M3 <- MortalityLaw(x = x, Dx = Dx, Ex = Ex, custom.law = my_gompertz) #' summary(M3) #' # predict M3 for different ages #' predict(M3, x = 85:130) #' #' #' # Example 4: -------------------------- #' # Fit Heligman-Pollard model for a single #' # year in the dataset between age 0 and 100 and build a life table. #' #' x <- 0:100 #' mx <- ahmd$mx[paste(x), "1950"] # select data
#' M4 <- MortalityLaw(x = x, mx = mx, law = 'HP', opt.method = 'LF2')
#' M4
#' plot(M4)
#'
#' LifeTable(x = x, qx = fitted(M4))
#' @export
MortalityLaw <- function(x,
Dx = NULL,
Ex = NULL,
mx = NULL,
qx = NULL,
law = NULL,
opt.method = 'LF2',
parS = NULL,
fit.this.x = x,
custom.law = NULL,
show = FALSE,
...){

info    <- addDetails(law, custom.law, parS)
law     <- info$law scale.x <- info$scale.x
parS    <- info$parS input <- c(as.list(environment())) K <- find.my.case(Dx, Ex, mx, qx) # TR: if inputs are matrix, then we have class matrix, array, and this # throws a warning. If we have a dim attribute then this won't work. Even # if it's a 1-column matrix (class =array) or a single-dim length-attribute vector # (also array!). Both of those cases are probably things we want to treat as # vectors! if (any(K$iclass == "numeric")) {

check.MortalityLaw(input) # Check input

# Set-up progress bar
if (show) {pb <- startpb(0, 4); on.exit(closepb(pb)); setpb(pb, 1)}

# Find optimal coefficients
optim.model <- choose_optim(input)
if (show) setpb(pb, 2)

fit   <- optim.model$hx dgn <- optim.model$opt #diagnosis
cf    <- optim.model$C p <- length(cf) resid <- switch( K$case,
C1_DxEx = Dx/Ex - fit,
C2_mx = mx - fit,
C3_qx = qx - fit
)
dev  <- sum(resid^2)
rdf  <- length(x) - p
df   <- c(n.param = p, df.residual = rdf)
gof  <- with(
optim.model,
c(logLik = logLik, AIC = AIC, BIC = BIC)
)
info <- list(model.info = info$model, process.date = date()) names(fit) = names(resid) <- x if (show) setpb(pb, 4) } else {# if input is a matrix then iterate here N <- K$nLT
# Set-up progress bar
if (show) {pb <- startpb(0, N + 1); on.exit(closepb(pb))}
cf = fit = gof = resid = dgn = df = dev <- NULL

for (i in 1:N) {
if (show) setpb(pb, i)
M <- suppressMessages(
MortalityLaw(
x = x,
Dx = Dx[, i],
Ex = Ex[, i],
mx = mx[, i],
qx = qx[, i],
law = law,
opt.method = opt.method,
parS = parS,
fit.this.x = fit.this.x,
custom.law = custom.law,
show = FALSE
)
)
fit      <- cbind(fit, fitted(M))
gof      <- rbind(gof, M$goodness.of.fit) dgn[[i]] <- M$dgn
cf       <- rbind(cf, coef(M))
resid    <- cbind(resid, resid(M))
df       <- rbind(df, M$df) dev <- c(dev, M$dev)
}

info <- M$info rownames(cf) = rownames(gof) = rownames(df) = names(dev) <- K$LTnames
dimnames(fit) = dimnames(resid) <- list(x, K$LTnames) if (show) setpb(pb, N + 1) } # Exit output <- list( input = input, info = info, coefficients = cf, fitted.values = fit, residuals = resid, goodness.of.fit = gof, opt.diagnosis = dgn, df = df, deviance = dev ) output$info$call <- match.call() out <- structure(class = "MortalityLaw", output) return(out) } #' Depending on the chosen mortality law, additional details need to be #' specified in order to be able to fit the models taking into account it's #' particularities. #' @inheritParams MortalityLaw #' @keywords internal addDetails <- function(law, custom.law = NULL, parS = NULL) { if (is.null(law) & is.null(custom.law)) { stop("Which mortality law do you intend to fit?", call. = FALSE) } if (!is.null(custom.law)) { law <- "custom.law" parS <- custom.law(1)$par
MI   <- "Custom Mortality Law"
sx   <- TRUE

} else {
law  <- law
parS <- parS
A    <- availableLaws(law)[["table"]]
MI   <- data.frame(A[A$CODE == law, ], row.names = "") sx <- as.logical(MI$SCALE_X)
}

out <- list(
law = law,
parS = parS,
model = MI,
scale.x = sx
)
return(out)
}

#' Function to be Optimize
#' @inheritParams MortalityLaw
#' @keywords internal
objective_fun <- function(par, x, Dx, Ex, mx, qx,
law, opt.method, custom.law){

C  <- find.my.case(Dx, Ex, mx, qx)$case mu <- eval(call(law, x, par = exp(par)))$hx
mu[is.infinite(mu)] <- 1

if (is.null(Ex))    Ex = 1
if (C == "C1_DxEx") nu = Dx/Ex
if (C == "C2_mx")   nu = Dx <- mx
if (C == "C3_qx")   nu = Dx <- qx

# compute likelihoods or loss functions
loss <- switch(
EXPR = opt.method,
poissonL  = -(Dx * log(mu) - mu*Ex),
binomialL = -(Dx * log(1 - exp(-mu)) - (Ex - Dx)*mu),
LF1       =  (1 - mu/nu)^2,
LF2       =  log(mu/nu)^2,
LF3       =  ((nu - mu)^2)/nu,
LF4       =  (nu - mu)^2,
LF5       =  (nu - mu) * log(nu/mu),
LF6       =  abs(nu - mu)
)

# Here I want to make sure that the optimisation algorithm is not returning
# NaN values when it converges (because that is possible).
loss[is.infinite(loss)] <- 10^5
if (sum(is.na(mu)) != 0) loss = loss + 10^5
out <- sum(loss, na.rm = TRUE)
# because nls.lm function requires a vector we have to do the following:
if (any(law %in% c('thiele', 'wittstein'))) out = loss

return(out)
}

#' Scaling method for x vector
#' @inheritParams MortalityLaw
#' @keywords internal
scale_x <- function(x) {
x - min(x) + 1
}

#' Select an optimizing method
#' @param input list of all inputs collected from MortalityLaw function
#' @keywords internal
choose_optim <- function(input){
with(as.list(input), {
# Subset the data
select.x <- x %in% fit.this.x

if (scale.x) {
new.fit.this.x = scale_x(fit.this.x)
d = fit.this.x - new.fit.this.x
new.x = x - d

} else {
new.fit.this.x <- fit.this.x
new.x <- x
}

if (is.null(parS)) parS <- bring_parameters(law, parS)
# Optimize
foo <- function(k) {
objective_fun(
par = k,
x = new.fit.this.x,
Dx = Dx[select.x],
Ex = Ex[select.x],
mx = mx[select.x],
qx = qx[select.x],
law,
opt.method,
custom.law)
}

if (any(law %in% c('HP', 'HP2', 'HP3', 'HP4', 'kostaki'))) {
opt <- nlminb(
start = log(parS),
objective = foo,
control = list(eval.max = 5000, iter.max = 5000)
)
opt$fnvalue <- opt$objective

} else if (any(law %in% c('thiele', 'wittstein'))) {
opt <- nls.lm(
par = log(parS),
fn = foo,
control = nls.lm.control(maxfev = 10000, maxiter = 1024)
)
opt$fnvalue <- sum(opt$fvec)

} else {
opt <- optim(
par = log(parS),
fn = foo,
)
opt$fnvalue <- opt$value
}

C <- exp(opt$par) if (law == 'kostaki') { #kostaki hack if (C >= 50*C) C <- C/50 } hx <- do.call(law, list(x = new.x, par = C))$hx
logLik <- log(opt\$fnvalue)
AIC    <- 2 * length(parS) - 2 * logLik
BIC    <- log(length(fit.this.x)) * length(parS) - 2 * logLik

if (!any(opt.method %in% c('poissonL', 'binomialL'))) {
logLik = AIC  = BIC  <- NaN
}
out <- as.list(environment())
return(out)
})
}


## Try the MortalityLaws package in your browser

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

MortalityLaws documentation built on July 1, 2022, 5:06 p.m.