morthump: Isolate the young adult mortality hump from a set of...
In MortHump: Measure the Young Adult Mortality Hump

Description Usage Arguments Details Value References See Also Examples

This function estimates a model of mortality on the provided set of age-specific death rates. Both parametric and non-parametric models are available. They are all designed to estimate the size of the young adult mortality hump, i.e. the deviation in the force of mortality often observed during adolescence and early adulthood.

1 2	morthump(data, model, method = "port", w = 1/data$m, start = NULL, maxit = 500, x1 = 30, x2 = 50, lambda.sen = 100, lambda.hump = 1)

`data`	data frame produced with `HMD2MH` or similarly structured
`model`	the name of the mortality model to be fitted
`method`	the optimization algorithm to be used
`w`	weights for the weighted least square estimation (parametric models only)
`start`	optional list of three elements giving the starting, lower and upper values of the parameters (parametric models only)
`maxit`	maximum number of iterations (sse.fit)
`x1`	upper age limit for the hump component (sse.fit)
`x2`	lower age limit for the senescence component (sse.fit)
`lambda.sen`	smoothing parameter for the senescence part (sse.fit)
`lambda.hump`	smoothing parameter for the hump component (sse.fit)

Mortality models can be broadly divided into parametric and non-parametric types.

Parametric models

Dozens of parametric models have been published (see Wunsch et al. 2002 for a review), but only three models are implemented here. They are all based on the structure of the Heligman-Pollard model and are conceived as nested models that offer a range of ways to address the specificity of young adult mortality. The "hps" model is the simplest and does not take the hump into account. The "hp" model assumes a symetrical hump with varying height (D), spread (E) and location (F). The "hpk" model relaxes the assumption of symetry by introducing a different spread before and after the peak of the hump. Algebraically, "hp" is equal to "hps" iff E = 0, and "hpk" is equal to "hp" iff E1 = E2.

"hp" : Published by Heligman and Pollard (1980), this model contains eight parameters in three additive terms.

A ^ ((x + B) ^ C) + D * exp(-E * (log(x) - log(F)) ^ 2) + (G * (H ^ x))/(1 + G * (H ^ x))

where (1) A, B and C describe the infant mortality component, (2) D, E and F describe the hump component, and (3) G and H describe the senescence component.

"hpk" : Published by Kostaki (1992), this extension of "hp" contains nine parameters in three additive terms.

A ^ ((x + B) ^ C) + D * exp(-E1 * (log(x) - log(F)) ^ 2) + (G * (H ^ x))/(1 + G * (H ^ x)) for x ≤ F

A ^ ((x + B) ^ C) + D * exp(-E2 * (log(x) - log(F)) ^ 2) + (G * (H ^ x))/(1 + G * (H ^ x)) for x > F

In contrast to the "hp" model, it allows the hump to be asymetrical by differenciating the spread of the hump before (E1) and after (E2) its peak.

"hps" : This model is an adaptation of the "hp" model that does not account for the young adult hump. It resembles in this respect the model published by Siler (1979), but is nested in the "hp" model.

A ^ ((x + B) ^ C) + D + (G * (H ^ x))/(1 + G * (H ^ x))

In contrast to the "hp" model, it only has one parameter for the hump component (D) that has the same function as the constant term of the Makeham model (1860).

Non-parametric model

The age-specific mortality rates can alternatively be fitted with non-parametric models, but very few of them adopt the multiple-component approach that is often found in parametric models. In order to isolate the young adult mortality hump from the rest of the force of mortality, a non-parametric model should also include the notion of additive components. This is what makes the so-called Sum of Smooth Exponentials model an attractive alternative to parametric models.

Briefly, the SSE model describes the observed force of mortality over age μ as the sum of three vectors <γ_1:γ_2:γ_3> over age. The subscripts (1 to 3) denote which mortality component each γ_j refers to: infant, midlife and old-ages mortality, respectively. In other words it assumes that deaths are realizations from Poisson distributions with mean composed of three parts:

y ~ P(e μ = C γ)

Given m age groups, the matrix C is a m by 3 matrix repeating the exposures in each column, which, when multiplied with the γ_j, sums the components and transforms the result into death counts. In matrix notation, C is given by

C = 1_{1,3} x diag(e)

where 1_{1,3} is a 1 by 3 matrix of ones, diag(e) is the diagonal matrix of the exposure population and x denotes the Kronecker product.

Unlike parametric models, in the SSE model there is no need to make strong assumptions about the functional form of each component. For each component a discrete sequence is assumed and the exponential function is used to ensure non-negative elements:

γ_j = exp(X_j β_j), j = 1,2,3

In other words, each component has to be described by a linear combination of a model matrix X_j and associated coefficients β_j. The design matrices X_j can represent parametric or non-parametric structures. In this way the composite force of mortality μ can be viewed as sum of 3 exponential components, which potentially can be smooth. Furthermore, the SSE model allows incorporating shape constraints to enforce senescent and young-adult components to be monotonically increasing and log-concave, respectively. A Penalized Composite Link Model (PCLM) has been proposed to estimate the SSE model. For more details, see Camarda et al. (2016).

Further Notes

All models use the age-specific mortality rates (mx) as the response variable in order to preserve comparability. This is despite the original definitions that often use the mortality quotients (qx) or one of its transforms.

By default, the parametric models are estimated with the "port" algorithm from the nls function. Is some rare cases, it may become stuck into local minimums. If this happens, try switching to the Levenberg-Marcquart algorithm. The non-parametric models use a tailored type of optimization algorithm (see sse.fit).

The optimization algorithm uses starting values, as well as lower and upper boundaries for each parameter. Most of the time, these values will be good enough to get a convergence. If you are not satisfied with the fitted values (for instance, if one of the fitted parameter is a round number there is a good chance that the algorithm encountered a boundary value), you can specify your own starting, lower and upper values. To change these values, use the start argument and construct it as a list of three objects: start, lower and upper. The first argument, start, defines the starting values for each of the parameters and is itself a named list (see example). The other two elements defined the lower and upper boundaries of the parameters and are defined as vectors.

Returns an object of class morthump containing the arguments used to fit the model as well as the estimated coefficients.

Camarda, C. G., Eilers, P. H. C., & Gampe, J. (2016). Sums of smooth exponentials to decompose complex series of counts. Statistical Modelling.

Wunsch, G., Mouchart, M., & Duchene, J. (2002). The life table : modelling survival and death. Dordrecht ; London: Kluwer Academic.

Heligman, L., & Pollard, J. H. (1980). The age pattern of mortality. Journal of the Institute of Actuaries, 107.

Kostaki, A. (1992). A nine-parameter version of the Heligman-Pollard formula. Mathematical Population Studies, 3(4), 277-288.

Siler, W. (1979). A Competing-Risk Model for Animal Mortality. Ecology, 60(4), 750-757.

sse.fit, summary.morthump, plot.morthump

data("CHE2010m")

# fits the Heligman-Pollard model (parametrical)
fit <- morthump(data = CHE2010m, model = "hp")

# use custom starting, lower and upper values for the parameters
sv <- list(start = list(A=0.001, B=0.005, C=0.11, D=0.0015, E=8, F=20, G=3e-5, H=1.105),
           lower = c(0.0001, 0.000001, 0.0001, 0, 1, 16, 0.0000001, 0.5),
           upper = c(0.1, 0.5, 1, 0.01, 50, 30, 0.01, 1.5))

fit <- morthump(data = CHE2010m, model = "hp", start = sv)

# fits a Sum of Smooth Exponentials (non-parametrical)
fit <- morthump(data = CHE2010m, model = "sse")