MortalityLaw: Fit Mortality Laws

Description

Fit parametric mortality models given a set of input data. The data can be supplied as death counts and mid-interval population estimates (Dx, Ex), age-specific death rates (mx), or death probabilities (qx). Use the law argument to specify the model to be fitted. Over 30 parametric models are currently implemented; run availableLaws to see the full list. Models can be fitted using maximum likelihood or by optimising a loss function. See the availableLF function for the implemented options.

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, ...)

Value

An object of class "MortalityLaw", which is a list with the following components:

input: List of input arguments, stored for reproducibility.
info: Model information (name, formula, date of fitting).
coefficients: Estimated parameters of the mortality law. A named vector for a single fit, or a matrix for multiple fits.
fitted.values: Fitted hazard rates (or death probabilities) evaluated at the input ages x.
residuals: Deviance residuals, computed as observed minus fitted values.
goodness.of.fit: List or matrix of goodness-of-fit measures: AIC, BIC and log-likelihood (available only for likelihood-based methods).
opt.diagnosis: Object returned by the optimisation routine, useful for checking convergence.
df: Number of parameters and residual degrees of freedom.
deviance: Sum of squared log-residuals, used as a deviance measure.

Arguments

x

Numeric vector of ages at the beginning of each age interval. For a full life table, use single-year ages (e.g., 0:110). For an abridged life table, use the lower bound of each interval (e.g., c(0, 1, 5, 10, ..., 110)).

Dx

Death counts. Each element represents the total number of deaths during the calendar year to persons aged x to x + n (where n is the length of the age interval). Must be provided together with Ex.

Ex

Exposure-to-risk in the period. This is usually approximated by the mid-year population aged x to x + n. Must be provided together with Dx.

mx

Age-specific death rate in the age interval [x, x+n). Defined as Dx / Ex.

qx

Probability of dying within the age interval [x, x+n).

law

The name of the mortality law to be used (e.g., "gompertz", "makeham"). Run availableLaws to see all options.

opt.method

The function to optimise. Available options:

"poissonL": Poisson log-likelihood.
"binomialL": Binomial log-likelihood.
"LF1": Squared relative error (1 - mu/nu)^2.
"LF2": Squared log-ratio log(mu/nu)^2.
"LF3": Chi-squared-type ((nu - mu)^2)/nu.
"LF4": Squared error (nu - mu)^2.
"LF5": Deviance-type (nu - mu) * log(nu/mu).
"LF6": Absolute error abs(nu - mu).

See availableLF for details.

parS

Optional starting parameter values for the optimisation. If NULL, sensible defaults are automatically chosen via bring_parameters.

fit.this.x

A subset of x over which to fit the model. The default is the entire x vector. Use this to exclude, for example, advanced ages where data are sparse.

custom.law

A user-defined function for fitting a model not included in the package. The function must accept arguments x (age vector) and par (named parameter vector) and return a list containing at least an element named hx (the hazard or force of mortality). See the examples below.

show

Logical. If TRUE, a progress bar is displayed during fitting. Default: FALSE.

...

Additional arguments passed to or from other methods.

Author

Marius D. Pascariu

Details

Optimisation: The PORT routines (via nlminb) are used for unconstrained and box-constrained optimisation. Parameters are estimated on the log scale to ensure positivity, and the routine is set to allow up to 5000 iterations. When the optimisation method is "poissonL" or "binomialL", the AIC, BIC and log-likelihood are computed from the likelihood. Otherwise these are set to NaN.

Scaling of the age vector: For models that cover only a portion of the lifespan (e.g., adult or old-age mortality), the age vector x is automatically re-scaled as x = x - min(x) + 1 before fitting. This transformation improves numerical stability and helps the optimisation algorithm converge, especially when the starting age is far from zero. Models that apply this scaling are flagged with SCALE_X = TRUE in the table returned by availableLaws. When using predict.MortalityLaw or LawTable with such models, the same scaling is applied internally, so predictions remain consistent with the fitted coefficients.

Handling matrix input: If Dx, Ex, mx or qx are provided as matrices (with one column per population or time period), the function iterates over the columns and fits a separate model to each, returning a collection of results.

Examples

Run this code

# Example 1: Fitting the Makeham model --------------------------
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))

Run the code above in your browser using DataLab