A black box approach to fitting smooth models of mortality

Abstract

Actuaries have long been interested in the forecasting of mortality for the purpose of the pricing and reserving of pensions and annuities. Most models of mortality in age and year of death, and often year of birth, are not identifiable so actuaries worried about what constraints should be used to give sensible estimates of the age and year of death parameters, and, if required, the year of birth parameters. These parameters were then forecast with an ARIMA model to give the required forecasts of mortality. A recent article showed that, while the fitted parameters were not identifiable, both the fitted and forecast mortalities were. This result holds if the age term is smoothed with P-splines. The present article deals with generalized linear models with a rank deficient regression matrix. We have two aims. First, we investigate the effect that different constraints have on the estimated regression coefficients. We show that it is possible to fit the model under different constraints in R without imposing any explicit constraints. R does all the necessary booking-keeping ‘under the bonnet’. The estimated regression coefficients under a particular set of constraints can then be recovered from the invariant fitted values. We have a black box approach to fitting the model subject to any set of constraints.

Keywords

constraints identifiability invariance mortality

1 Introduction

Actuaries have long been interested in the forecasting of human mortality for the purpose of the pricing and reserving of pensions and annuities. Gompertz (1825), in a landmark paper, gave the first mathematical model of mortality. He only had death rates by age at ten year intervals but he observed that these rates were linear on a log scale. This enabled him to interpolate for all ages for the range of ages he had available. This was the celebrated Gompertz law of mortality.

At the present time extensive mortality data are available not only by single ages but also by year of death. The Human Mortality Database is a comprehensive resource of mortality data and gives mortality data (a) by age and year of death in a calendar year, and (b) an estimate of the mid-year number of lives exposed to the risk of death in that year. These data are available on over 40 countries for the male and female populations of these countries.

These richer datasets gave rise to the possibility not only of modelling mortality by age and year but also the forecasting of mortality. Such models were in terms of age and year of death, and often year of birth. The idea was to estimate the parameters in the models and then forecast the year of death parameters and, if required, the year of birth parameters; see Cairns et al. (2009) for a comprehensive treatment of such models. However, there was a problem—most models were not identifiable so it was not obvious how to obtain parameter estimates that would give sensible forecasts. Various sets of linear constraints on the parameters were proposed that would solve this problem.

Most models were generalized linear models or GLMs. Currie (2020) showed that, while the fitted parameters were not identifiable, the fitted mortalities were identifiable (a fundamental result in regression modelling) but crucially the forecast values of mortality were also identifiable when forecasting was done with an ARIMA model. This article illustrated this result by considering two sets of constraints, one a standard set of constraints found in the literature (see Cairns et al., 2009, for example) and the other a set of random constraints. A general proof of the result was also provided.

The purpose of the present article is twofold: first, to show that we can ignore the explicit use of constraints completely and use R’s lm and glm commands to fit the models, and second, to use the invariance of the fitted mortalities to explore the relationship between the estimates of the parameters under different sets of constraints. R does all the necessary booking-keeping ‘under the bonnet’. We have a black box approach to the fitting and forecasting of mortality.

We do not present any general theory but rather present our method through a single model, the Age-Period-Cohort model. We consider this model in a basic form and also in a smooth form where we use P-splines (Eilers an Marx 1996) to smooth some terms in the model. It should be clear that our method applies more generally.

2 Data and notation

We illustrate our ideas with data from the Human Mortality Database on UK males (data downloaded February 2022). We adopt the convention that matrices are denoted by upper case bold font as in M and column vectors by lower case bold font as in m . The data comprise two matrices: $D_{o b s} = (d_{x, y})$ , the matrix of the number of deaths at age $x$ in year $y$ , and $E_{o b s} = (e_{x, y})$ , the total time lived or central exposure at age $x$ in year $y$ . (We will use $D$ later to denote a difference matrix.) In general we suppose we have $n_{a}$ ages and $n_{y}$ years. We denote the age indices by $x_{a} = (1, \dots, n_{a})^{'}$ and the year indices by $x_{y} = (1, \dots, n_{y})^{'}$ where the $^{'}$ indicates the transpose of a vector (or matrix). We index the rows and columns of $D_{o b s}$ and $E_{o b s}$ by $x_{a}$ and $x_{y}$ respectively. Thus $D_{o b s}$ and $E_{o b s}$ are both $n_{a} \times n_{y}$ . Further we index the cohorts or years of birth with $x_{c} = (1, \dots, n_{c})^{'}$ where $n_{c} = n_{a} + n_{y} - 1$ is the number of distinct cohorts. We adopt the convention that the oldest cohort, ie, for age $n_{a}$ in year 1, is indexed 1; hence, the cohort index for cell $(x, y)$ is $c (x, y) = n_{a} - x + y$ . In our case we have data on ages 50 to 90 and years 1970 to 2018. Hence $n_{a} = 41$ , $n_{y} = 49$ and $n_{c} = 89$ .

In addition, the following notation is useful: $1_{n}$ denotes a column vector of 1s of length $n$ and $I_{n}$ denotes the identity matrix of size $n$ . Further we let $0_{n}$ stand for either a column vector of 0s of length $n$ or an $n \times n$ zero matrix; the context should make clear which applies. We may omit the suffix $n$ if no confusion results. We will make frequent use of the Kronecker product which we denote $\otimes$ . The function $v e c (M)$ stacks the columns of a matrix $M$ on top of each other in column order.

3 The Age-Period-Cohort model

Let $D_{x, y}$ be the random variable corresponding to the observed deaths $d_{x, y}$ . We suppose that $D_{x, y}$ has the Poisson distribution with mean $e_{x, y} λ_{x, y}$ where $λ_{x, y}$ is the force of mortality or instantaneous death rate at age $x$ in year $y$ , ie, $D_{x, y} \sim P (e_{x, y} λ_{x, y})$ . There is a slight abuse of notation here; strictly we should write $λ_{x + 1 / 2, y + 1 / 2}$ since $e_{x, y}$ is the central exposure; we will use the simpler notation throughout for clarity. Let $d = v e c (D_{o b s})$ and $e = v e c (E_{o b s})$ be the stacked vectors of deaths and exposures; let $λ$ be the corresponding vector of forces of mortality.

The Age-Period-Cohort or APC model is the simplest model which involves the age and year of death, and the year of birth. We define

\log λ_{x, y} = α_{​ x} + κ_{​ y} + γ_{c (x, y)}, x = 1, \dots, n_{a}, y = 1, \dots, n_{y},

(3.1)

where $c (x, y) = n_{a} - x + y$ is the cohort index for age $x$ in year $y$ . Let $θ = (α^{'}, κ^{'}, γ^{'})^{'}$ . The model matrix is

X = [X_{a} : X_{y} : X_{c}], n_{a} n_{y} \times (2 n_{a} + 2 n_{y} - 1),

(3.2)

where $X_{a} = 1_{n_{y}} \otimes I_{n_{a}}$ , $X_{y} = I_{n_{y}} \otimes 1_{n_{a}}$ and the row of $X_{c}$ corresponding to age $x$ and year $y$ contains a one in column $n_{a} - x + y$ and zeros elsewhere. With this setup we have $\log λ = X θ$ . We have defined a generalized linear model or GLM with dependent variable $d$ , model matrix $X$ , Poisson errors and log link; the exposures $e$ enter into the model as an offset. Now $X$ is $n_{a} n_{y} \times (2 n_{a} + 2 n_{y} - 1)$ but its rank is $2 n_{a} + 2 n_{y} - 4$ .

The usual actuarial approach is to choose three suitable constraints and then forecast the resulting estimates of $κ$ and $γ$ . A common choice is

\sum_{1}^{n_{y}} κ_{y} = \sum_{1}^{n_{c}} γ_{c} = \sum_{1}^{n_{c}} c γ_{c} = 0

(3.3)

see Cairns et al. (2009). We refer to (3.3) as the standard constraints. We can fit the APC model in R with the glm command. Here Dth and Exp contain the R variates holding the deaths $d$ and the exposures $e$ respectively. We have

Fit.APC <- glm(Dth \sim -1+X+offset(log(Exp)),family=poisson)

(3.4)

Theta <- Fit.APC$coef

(3.5)

Theta[is.na(Theta)] <- 0

(3.6)

Inspection of Fit.APC$coef shows that the estimates of $κ_{n_{y}}$ , $γ_{n_{a} + n_{y} - 2}$ and $γ_{n_{a} + n_{y} - 1}$ are returned as NA.

What is the connection between the approach to identifiability used by R and that often used by actuaries? We can mimic R’s method with the constraints

κ_{y} = γ_{n_{c} - 1} = γ_{n_{c}} = 0,

(3.7)

ie, we constrain the final period parameter and the final two cohort parameters to be zero. We refer to these constraints, implicit in R’s approach, as the R constraints. We refit the APC model with the constraints (3.7). Currie (2013) gave the following generalization of the GLM scoring algorithm to fit a GLM with constraints $H θ = 0$ but note that this algorithm plays no part in the treatment presented in this article

(\begin{matrix} X^{'} \tilde{W} X & : & H^{'} \\ H & : & 0 \end{matrix}) (\begin{matrix} \hat{θ} \\ \hat{ω} \end{matrix}) = (\begin{matrix} X^{'} \tilde{W} \tilde{z} \\ 0 \end{matrix}) .

(3.8)

Here $\tilde{W}$ is the diagonal matrix of weights, $\tilde{z}$ is the working variable, $\sim$ indicates a current estimate and ^ denotes an updated estimate. The fitted coefficients are precisely those in Fit.APC$coef but with the NA values replaced by zeros. Equation (3.8) is used simply to verify this equality.

We have fitted the model without the need to specify constraints. If we wish for some reason to find a fit with specified constraints we proceed as follows. We seek a solution of the GLM scoring algorithm

(X^{'} \tilde{W} X) \hat{θ} = X^{'} \tilde{W} \tilde{z} .

(3.9)

Let ${\hat{θ}}_{1}$ and ${\hat{θ}}_{2}$ be two solutions of (3.9) where ${\hat{θ}}_{1}$ satisfies the known constraint $H {\hat{θ}}_{1} = 0$ and ${\hat{θ}}_{2}$ satisfies some unspecified constraint. We know

(X^{'} \tilde{W} X) {\hat{θ}}_{1} = (X^{'} \tilde{W} X) {\hat{θ}}_{2}

\Rightarrow (X^{'} \tilde{W} X + H^{'} H) {\hat{θ}}_{1} = (X^{'} \tilde{W} X) {\hat{θ}}_{2} .

Now $X^{'} \tilde{W} X + H^{'} H$ is positive definite so if we can find a ${\hat{θ}}_{2}$ we have the solution which satisfies $H θ = 0$ as

{\hat{θ}}_{1} = (X^{'} \tilde{W} X + H^{'} H)^{- 1} X^{'} \tilde{W} X {\hat{θ}}_{2} .

(3.10)

This expression exists in another form which is particularly illuminating. We observe that $X {\hat{θ}}_{2} = \log \hat{λ}$ and so we can write (3.10) as

{\hat{θ}}_{1} = (X^{'} \tilde{W} X + H^{'} H)^{- 1} X^{'} \tilde{W} \log \hat{λ} .

(3.11)

In other words, we can compute the estimated coefficients subject to any particular set of constraints directly from the invariant fitted values. Invariance is a fundamental and powerful idea and (11) demonstrates just one result from applying it.

Figure 1 emphasizes the effect that constraints can have on parameter estimates. The effect on the estimates of the time parameters is particularly striking. Under the standard constraints one may be tempted to interpret the parameter estimates as indicating improving mortality. The estimates under the R constraints indicate exactly the opposite. What can be said is that mortality does improve over time but that can only be concluded securely from the invariant fitted mortality shown in the lower right panel.

Figure 1

Estimates of $α$ , $κ$ and $γ$ under standard constraints (solid line), R constraints (dashed line) and the observed and invariant fitted log mortality for age 70.

We examine the differences in the parameter estimates. We define $Δ α = {\hat{α}}_{s} - {\hat{α}}_{r}$ where the suffixes indicate the estimates under the standard and R constraints. We have corresponding definitions for $Δ κ$ and $Δ γ$ . Currie (2020) showed that these diferences obeyed a very simple relation. We have

Δ α = (A + B) 1_{n_{a}} + C x_{a}

(3.12)

Δ κ = - A 1_{n_{y}} - C x_{y}

(3.13)

Δ γ = - (B + n_{a} C) 1_{n_{c}} + C x_{c}

(3.14)

for some constants A, B and C. We conclude that the three sets of differences are linearly related with slopes C, -C and C. With our data and the standard and R constraints we found A=-1.163 B=1.465 and C=0.03199. Despite appearances to the contrary in Figure 1 the two sets of estimates are intimately related.

4 The smooth Age-Period-Cohort model

Gompertz smoothed the age parameters with a simple linear function. Here we use the method of P-splines to smooth the age parameters, $α$ , in the APC model. As well as the general consideration that we expect the age effect to be a smooth monotonic function, there are good actuarial reasons for smoothing the age parameters. The resulting estimates in time by age will be more regular and hence the pricing of insurance products will also be more regular.

Let $B_{a}, n_{a} \times c_{a},$ be a regression matrix of B-splines along age. The model matrix is as in (2) but with $X_{a} = 1_{n_{y}} \otimes I_{n_{a}}$ replaced by $1_{n_{y}} \otimes B_{a}$ . We define the penalty matrix $P = b l o c k d i a g {τ D_{2}^{'} D_{2}, 0_{n_{y}}, 0_{n_{c}}}$ . Here $D_{2}$ is a second order difference matrix and $τ$ is a tuning constant. See Eilers and Marx (2021) for a thorough introduction to the method of P-splines. We seek solutions of the penalized scoring algorithm

(X^{'} \tilde{W} X + P) \hat{θ} = X^{'} \tilde{W} \tilde{z}

(4.1)

where

\tilde{W} = diag {\tilde{d}}, \tilde{z} = X \tilde{θ} + d / \tilde{d} - 1 .

(4.2)

for a GLM with Poisson errors and a log link.

We could use a generalization of (3.8) where the upper left block becomes $X^{'} \tilde{W} X + P + H^{'} H$ . However, our purpose here is to use R to fit the model and then use invariance to calculate the estimates of the regression coefficients under the standard and R constraints. We will see that many of the results from the unsmoothed model carry over in a natural way.

We can avoid the use of algorithms like (3.8) with data augmentation; see Eilers and Marx (2021), p35. Data augmentation enables us to fit the smooth APC model with code similar to (3.4), (3.5) and (3.6) for the original APC model discussed in section 3. We define the augmented model matrix and augmented working variable as follows:

X_{Aug} = [\begin{matrix} {\tilde{W}}^{1 / 2} X_{a} & : & {\tilde{W}}^{1 / 2} X_{y} & : & {\tilde{W}}^{1 / 2} X_{c} \\ \sqrt{τ} D_{2} & : & 0_{(c_{a} - 2) \times n_{y}} & : & 0_{(c_{a} - 2) \times n_{c}} \end{matrix}], z_{Aug} = [\begin{matrix} {\tilde{W}}^{1 / 2} z \\ 0_{c_{a} - 2} \end{matrix}] .

(4.3)

It follows immediately that

X_{Aug}^{'} \tilde{W} X_{Aug} = X^{'} \tilde{W} X + P, X_{Aug}^{'} \tilde{W} z_{Aug} = X^{'} \tilde{W} \tilde{z} .

(4.4)

Thus solutions of the penalized scoring algorithm (4.1) and

X_{Aug}^{'} \tilde{W} X_{Aug} \hat{θ} = X_{Aug}^{'} \tilde{W} {\tilde{z}}_{Aug}

(4.5)

are equal. Equation (4.5) is weighted least squares so we can fit the model with R’s lm function. We note that we must iterate with (4.2) until convergence. Good initial values are the observed deaths for $d$ and the observed log mortality rates for $z$ . Convergence is very fast since we have essentially implemented the GLM scoring algorithm via the lm function.

We choose $τ$ by minimizing AIC. The AIC curve is quite flat in the region of the optimal value and a short grid search gives $τ = 100$ as a suitable value; this value is used in what follows.

We use the following extension of (3.11)

{\hat{θ}}_{1} = (X^{'} \tilde{W} X + H^{'} H + P)^{- 1} X^{'} \tilde{W} \log \hat{λ}

(4.6)

to investigate the relation between the two sets of estimates obtained under the standard and R constraints. We put $H$ equal to $H_{s}$ and $H_{r}$ in turn where $H_{s}$ and $H_{r}$ are the constraints matrices for the standard and R constraints.

As in the unsmoothed case, we define $Δ_{a}$ , $Δ_{κ}$ and $Δ_{γ}$ as the differences in the two sets of estimates but note here that $Δ_{a}$ refers to the coefficients of the B-spline basis. We do not have the general theory that led to equations (3.12), (3.13) and (3.14) when $τ > 0$ ; Currie (2020) gives the theory when $τ = 0$ . We proceed numerically. Plots of $Δ_{a}$ , $Δ_{κ}$ and $Δ_{γ}$ suggest that each is linear and this is confirmed by simple regression. We find in our example that the slopes are 0.17915, $- 0.03583$ and 0.03583 respectively. We note that $0.17915 / 0.03583 = 5$ . We have used a knot spacing of five which suggests that the slopes are $Δ C$ , $- C$ and $C$ where $Δ$ is the knot spacing. Support for this conjecture comes from Currie (2020) where it is shown that the slopes do indeed satisfy $Δ C$ , $- C$ and $C$ in the case $τ = 0$ .

Figure 2 shows the observed and fitted log mortalities for age 50 to 90 at ten year intervals. Evidently, there has been a steady improvement in mortality over the observation period. This improvement corresponds to an approximate increase in life expectancy of two years every ten years. For example, life expectancy in the UK in 1970 was 72.3 years and this has risen to 81.2 years in 2019. This has major consequences for the provision of both state and private pensions.

Figure 2

Observed and fitted log $λ$ for ages 50 to 90 at ten year intervals for the smooth APC model.

5 Discussion

This short article has highlighted the effect that different constraint systems can have on parameter estimates. Figure 1 is a dramatic example of such effects. We have also highlighted the role of the invariant fitted values. Equation (3.11) shows that knowledge of these invariant values yields parameter estimates under any particular constraint system. We have used the basic lm function in R as a black box to estimate the invariant fitted values and from these values diagrams like Figure 1 can be obtained and equation (3.11) can be applied.

Footnotes

Acknowledgements

I am grateful for the invitation to submit an article for this special issue. I met Brian on many occasions: it would be hard to imagine anyone with a greater zest for life in general and statistics in particular. I am also grateful for his great paper with Paul Eilers which introduced the world to P-splines. It has dominated my working life for over twenty years. I would also like to thank Paul Eilers for his encouragement to write this short article.

Declaration of Conflicting Interests

The author declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author received no financial support for the research, authorship and/or publication of this article.

References

Currie

(2013) Smoothing constrained generalized linear models with an application to the Lee-Carter model. Statistical Modelling , 13, 69–93.

Cairns

AJG

, Blake

, Dowd

, Coughlan

, . (2009). A quantitative comparison of stochastic mortality models using data from England and Wales and the United States. North American Actuarial Journal , 13, 1–35.

Currie

(2020) Constraints, the identifiability problem and the forecasting of mortality. Annals of Actuarial Science , 14, 537–566.

Eilers

PHC

and Marx

(1996) Flexible smoothing with B-splines and penalties. Statistical Science , 11, 89–121.

Eilers

PHC

and Marx

(2021) Practical smoothing: the joys of P-splines . Cambridge: Cambridge University Press.

Gomperts

(1825) On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society of London , 115, 513–583.

Human Mortality Database (2023) University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany). Available at www.mortality.org.

R Core Team (2019) R: A language and environment for statistical computing . R Foundation for Statistical Computing, Vienna, Austria. URL https://www.Rproject.org/.