MortalityLaw function

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 (Dx, Ex) or age-specific death rates (mx) or death probabilities (qx). Using the argument law one can specify the model to be fitted. So far more than 27 parametric models have been implemented; check the 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 availableLF function.

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

Arguments

• x: Vector of ages at the beginning of the age interval.

• Dx: Object containing death counts. An element of the Dx object represents the number of deaths during the year to persons aged x to x+n.

• Ex: Exposure in the period. Ex can be approximated by the mid-year population aged x to x+n.

• mx: Life table death rate in age interval [x, x+n).

• qx: Probability of dying in age interval [x, x+n).

• law: The name of the mortality law/model to be used. e.g. gompertz, makeham, ... To investigate all the possible options, see availableLaws function.

• opt.method: How would you like to find the parameters? Specify the function to be optimize. Available options: the Poisson likelihood function poissonL; the Binomial likelihood function -binomialL; and 6 other loss functions. For more details, check the availableLF

  function.

• parS: Starting parameters used in the optimization process (optional).

• fit.this.x: Select the ages to be considered in model fitting. By default fit.this.x = x. One may want to exclude from the fitting procedure, say, the advanced ages where the data is sparse.

• custom.law: Allows you to fit a model that is not defined in the package. Accepts as input a function.

• show: Choose whether to display a progress bar during the fitting process. Logical. Default: FALSE.

• ...: Arguments to be passed to or from other methods.

Returns

The output is of the "MortalityLaw" class with the components: - input: List with arguments provided in input. Saved for convenience.

• info: Brief information about the model.

• coefficients: Estimated coefficients.

• fitted.values: Fitted values of the selected model.

• residuals: Deviance residuals.

• goodness.of.fit: List containing goodness of fit measures like AIC, BIC and log-Likelihood.

• opt.diagnosis: Resultant optimization object useful for checking the convergence etc.

Details

Depending on the complexity of the model, one of following optimization strategies is employed:

1. 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 optim;
2. PORT routines: provides unconstrained optimization and optimization subject to box constraints for complicated functions. For details check nlminb;
3. Levenberg-Marquardt algorithm: damped least-squares method. For details check nls.lm.

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

See Also

availableLaws

availableLF

LifeTable

ReadHMD

Author(s)

Marius D. Pascariu