2. Background and theory

In this section, we set up the notation and describe the theoretical background of our modeling and inference procedure. We assume to have data covering ages $x_{1},\ldots ,x_{n_X}$ and periods $t_1,\ldots ,t_{n_T}$ , where $n_X$ and $n_T$ are the total number of ages and periods, respectively. The number of cohorts is denoted by $n_C = n_X + n_T - 1$ . We will denote by $d_{x,t}$ and $E_{x,t}$ the number of deaths and exposure for age $x$ and period $t$ , respectively. We model deaths using the Poisson framework (Brouhns et al., Reference Brouhns, Denuit and Vermunt2002), $d_{x,t} \sim \text{Poisson}(E_{x,t} m_{x,t})$ , where $m_{x,t}$ is the mortality rate for the $x$ -th age and the $t$ -th period. In the case of APP data, we add the subscript $p$ to all the previous quantities to denote product, for example, $d_{x,t,p}$ denotes the number of deaths for the $x$ -th age, $t$ -th period and $p$ -th and product for APP data. The total number of products is denoted by $n_P$ .

We perform model selection in a Bayesian framework. To set up some notations, we denote the data as $d$ , the parameters as $\theta$ , and the model as $M$ , and we denote the posterior distribution of the model $M$ and the parameters $\theta$ conditional on the data $d$ by $p(\theta , M | d)$ , which is proportional to the likelihood, $p(d | \theta , M)$ , and the prior, $p(\theta , M)$ . Since we are performing model selection, we assume to have a set of models $(M_1,\ldots ,M_K)$ , each with a set of parameter $\theta _k$ , and we denote by $n_k$ the dimension of model $M_k$ . We denote by $p(\theta _k, M_k) = p(\theta _k | M_k) p(M_k)$ the joint prior on the model and the corresponding model parameters. We note that choosing a prior on the model of the form $p(M_k) \propto \frac {1}{n_k}$ corresponds to a posterior distribution with the same form as BIC, while choosing $p(M_k) \propto \text{exp}(-n_k)$ corresponds to AIC, since the log-posterior distribution $\ell + \log(p(M_k))$ is equal to $\ell - \log(n_k)$ in the first case and $\ell - n_k$ in the second case (with $\ell$ being the log-likelihood), which corresponds to the expression of the BIC and AIC, respectively. In all the following, we have assumed a prior $p(M_k)$ corresponding to the AIC penalty.

When proposing a move between models of different dimensions, the standard Metropolis-Hastings algorithm cannot be used anymore because of the difference in dimensions between models, and instead, we resort to the RJMCMC algorithm (Green, Reference Green1995), which we briefly review. The RJMCMC relies on defining two proposal distributions:

• • The move from model $M_k$ to $M_{k'}$ happens with probability $q(k \rightarrow k')$ .

• • Once a move is proposed, to define the parameter of the new proposed model, for each proposed model transition from $k$ to $k'$ , a latent variable $u$ of dimension $\zeta _{k \rightarrow k'}$ is sampled from a distribution $q_{k \rightarrow k'}(u)$ , such that the new set of parameters $\theta _{k'}$ of the proposed model $k'$ can be obtained as $\theta _{k'} = g_{k \rightarrow k'}(\theta _{k},u)$ .

Finally, the new model is accepted with probability

(1) \begin{equation} \alpha [(k, \theta _k),(k', \theta _{k'})] = \text{min}\left \{ 1, \frac { p(d | \theta _{k'}, M_{k'}) p(\theta _{k'} | M_{k'}) p(M_{k'}) q(k' \rightarrow k) q_{k' \rightarrow k}(u')}{ p(d | \theta _{k}, M_{k}) p(\theta _k | M_k) p(M_k) q(k \rightarrow k') q_{k \rightarrow k'}(u)} \left | \frac {\partial g_{k \rightarrow k'}(\theta _k, u)}{\partial (\theta _k, u)} \right | \right \} \end{equation}

where $q(k' \rightarrow k)$ and $q_{k' \rightarrow k}(u')$ are the transition probabilities of the reverse move and $u'$ is defined such that $\theta _k^{\prime } = g_{k \rightarrow k'}(g_{k' \rightarrow k}(\theta _k^{\prime },u'),u)$ . A requirement of RJMCMC is the dimension matching condition: $n_k + \zeta _{k \rightarrow k'} = n_{k'} + \zeta _{k' \rightarrow k}$ . It follows that in the case where $n_k + \zeta _{k \rightarrow k'} = n_{k'}$ , we have $\zeta _{k' \rightarrow k} = 0$ , with the last expression meaning that the reverse move to the smaller model is deterministic. For example, when models are nested, if we can find a proposal (from the smaller to the larger model) with a dimension equal to the difference in the model dimensions, the proposal to the smaller model will be deterministic.

Generally, the closer the proposal distribution $q_{k \rightarrow k'}(u)$ to the posterior distribution of the variables $u$ , the more efficient the sampling scheme will be. To achieve this, we choose as proposal distributions $q(u)$ a Laplace approximation of the posterior distribution. The Laplace approximation is given by a normal distribution $\text{N}(\mu ^{\star }, \Sigma ^{\star })$ , where $\mu ^{\star }$ is the set of a parameters corresponding to the maximum of the posterior distribution, and $\Sigma ^{\star }$ is the inverse of the negative Hessian of the posterior distribution computed in $\mu _{\star }$ . We note that since this corresponds to solving an optimization tasks for each proposal, it is preferred to choose moves with the smallest possible dimension changes $\zeta _{k \rightarrow k'}$ and $\zeta _{k' \rightarrow k}$ , since this will require performing optimization tasks over smaller number of variables.

We give a simple example of how RJMCMC works in practice. Let us consider two models for age-only data: model $1$ , having mortality as a constant $\text{log}(m_{x}) = a$ , and model $2$ , having a age-varying mortality $\text{log}(m_{x}) = a_x$ . We could propose the move in two different ways. In the first case, when going from model $1$ to model $2$ , we can propose $a_x$ from a proposal distribution $q_{1\rightarrow 2}(a_x)$ (of dimension $n_X$ ), and for the reverse move (from model $2$ to model $1$ ), we propose $a$ from a proposal distribution $q_{2\rightarrow 1}(a)$ (of dimension $1$ ). The dimension matching condition is satisfied as $n_1 (= 1) + \zeta _{1 \rightarrow 2} (= n_X) = n_2 (= n_X) + \zeta _{2 \rightarrow 1} (= 1)$ . Alternatively, we could design the move from model $1$ to model $2$ by only proposing the increments $\epsilon _x$ from $a$ to $a_x$ , such that $a_x = a + \epsilon _x$ , where $\sum _x \epsilon _x = 0$ (and therefore, $\epsilon _x$ is of dimension $n_X - 1$ ). The dimension matching condition is again satisfied since $n_1 (= 1) + \zeta _{1 \rightarrow 2} (= n_X - 1) = n_2 (= n_X)$ , and the reverse move is deterministic, where the proposed $a$ is the mean of the current $a_x$ . As mentioned, the second strategy is likely to be more efficient, since it only requires finding the maximum of one posterior distribution on a smaller dimensional space ( $X -1$ ) against finding the maximum of two posterior distributions (of dimension $1$ and $n_X$ ).

3. Mortality model selection for age-period data

In this section, we will define a modeling framework for building mortality models for age and period data and will next use this framework to perform model selection. We define our framework by separating the log-mortality rate into three terms: age-effect, period and age-period effect, and cohort effect, which have already been considered in a wide variety of models (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009).

We assume that the log-mortality rate can be expressed as

(2) \begin{equation} \text{log}(m_{x,t}) = f_1(x) + f_2(t, f_3(x)) + f_4(t-x) \end{equation}

where $f_1(x)$ is the age effect, $f_2(t, f_3(x))$ encompasses the period and the age-period effects, and $f_4(t-x)$ is the cohort effect. We model each term in the following way:

• • Age effect: we assume that age effect can be modeled linearly ( $f_1(x) = a + b (x - \bar {x})$ ), where $\bar {x} = \frac {1}{n_X} \sum _{i=1}^{n_X} x_i$ , or nonlinearly ( $f_1(x) = a_x$ ). We introduce the variable $\delta _1$ , which has value $1$ if the age effect is not linear and $2$ if the effect is linear.

• • Period effect: we assume that the period effect can be modeled in three ways and introduce the variable $\delta _2$ to represent each outcome: independent from age $(f_2(t,f_3(x)) = k^1_t)$ ( $\delta _2 = 1$ ), modeled through an interaction with age $(f_2(t,f_3(x)) = k^2_t (1 + f_3(x))$ ( $\delta _2 = 2$ ), and a sum of the previous two $(f_2(t,f_3(x)) = k^1_t + k^2_t f_3(x))$ , ( $\delta _2 = 3$ ).

• • Age-period effect: if we have chosen an interaction between age and period ( $\delta _2 = 2$ or $3$ ), we introduce the variable $\delta _3$ that is $1$ if the age term in the interaction is modeled linearly with respect to age ( $f_3(x) = \bar {b}(x - \bar {x})$ , where $\bar {b}$ is the slope) and $2$ if modeled nonlinearly ( $f_3(x) = b_x$ ).

• • Cohort effect: we introduce the variable $\delta _4$ that is $2$ if there is a cohort effect $\gamma _{c}$ in the model ( $f_4(c) = \gamma _c$ ), where $c = t - x$ , and $1$ otherwise ( $f_4(c) = 0$ ).

After reparametrization, the terms in the log-mortality rate in (2) can be expressed as

\begin{equation*} f_1(x) = \begin{cases} a_x \hspace {1.81cm} \text{if } \delta _1 = 1 \\[5pt] a + b (x - \bar {x}) \quad \text{if } \delta _1 = 2 \end{cases} \hspace {1cm} f_2(t, f_3(x)) = \begin{cases} k^1_t \hspace {1.9cm} \text{if } \delta _2 = 1 \\[5pt] k^2_t (1 + f_3(x)) \hspace {.4cm} \text{if } \delta _2 = 2 \\[5pt] k^1_t + k^2_t f_3(x) \hspace {.5cm} \text{if } \delta _2 = 3 \end{cases} \end{equation*}

\begin{equation*} f_3(x) = \begin{cases} \bar {b} (x - \bar {x}) \hspace {.42cm} \text{if } \delta _3 = 1 \\[5pt] b_x \hspace {1.3cm} \text{if } \delta _3 = 2 \end{cases} \hspace {1cm} f_4(t - x) = \begin{cases} 0 \hspace {.9cm} \text{if } \delta _4 = 1 \\[5pt] \gamma _{t-x} \hspace {.42cm} \text{if } \delta _4 = 2 \end{cases} \end{equation*}

The complete list of models that can be obtained as a result of the different choices is summarized in Table 1. It can be seen that models such as the APC, LC, and CBD, are also included. However, as mentioned, we are also able to achieve intermediate models between the ones already considered.

Table 1.

Complete set of models that can be obtained in the framework for each choice of $\delta _1$ , $\delta _2$ , $\delta _3$ , and $\delta _4$ . We have also denoted in brackets the models defined in Table 1 of Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009)

The set of identifiability constraints (ICs) are:

• • $\sum _t k^1_t = 0$ , since any constant added to $k^1_t$ can be absorbed into $f_1(x)$ .

• • $\sum _t k^2_t = 0$ , since we can remove a constant $C$ from $k^2_t$ and absorb $C f_3(x)$ into $f_1(x)$ . However, if $\delta _2 = \delta _3 = 2$ , this constraint is not strictly necessary since the nonlinear term $C b_x$ cannot be absorbed into the linear term $a + b (x - \bar {x})$ . For reasons described later, we assume the constraint throughout anyway.

• • $\sum _{x} b_x = 0$ . If $\delta _2 = 3$ , we also need on $b_x$ the additional constraint $b_{x_1} = -1$ , since we can otherwise replace $b_x$ with $C b_x$ and absorb the rest of the constants into the other parameters.

• • $\sum _{c} \gamma _{c} = 0$ , since any additive constant can be absorbed into $f_1(x)$ .

• • If $\delta _2 = 1$ or $\delta _2 = 3$ , following Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009), another constraint is needed since we can remove a term $C ( (t - \bar {t}) - (x - \bar {x}))$ from $\gamma _{t-x}$ , for any real constant $C$ , where $\bar {t} = \frac {1}{n_T} \sum _{i=1}^{n_T} t_i$ , and add $- C (x - \bar {x})$ to $a_x$ and $ C (t - \bar {t})$ to $k_t$ . We, therefore, choose to assume no linear trend by enforcing $\sum c \gamma _c = 0$ . We note that although not immediately apparent, it is necessary to assume another constraint on $\gamma _c$ also if only a term of the form $k_t b_x$ is present ( $\delta _2 = 2$ ) and therefore we choose to impose the previous constraint also for $\delta _2 = 2$ .

• • When $\delta _1=1,\delta _2=3,\delta _3=1$ , we also need another constraint on the cohort effect, since we can add a quadratic term $a_1 + a_2 ( (t - \bar {t}) - (x - \bar {x})) + a_2 ((t - \bar {t}) - (x - \bar {x}))^2$ (but still with mean zero and no linear trend) to $\gamma _{t-x}$ , and remove it from $f_1(x)$ , $k^1_t$ and $k^2_t$ , obtaining an equivalent model. We therefore impose $\sum c^2 \gamma _{c} = 0$ .

We note that our constraints are slightly different than those usually used in the literature because of the different parametrizations we have used.

3.1 Inference

We sample from the posterior distribution of $(\theta _1,\delta _1,\theta _2,\delta _2,\theta _3,\delta _3,\theta _4,\delta _4)$ , where $\theta _i$ is the set of parameters of $f_i$ , using an MCMC algorithm. The algorithm will alternate between sampling each set of parameters $\theta _i$ and the model configuration $\delta _i$ of each term in (2), conditional on the rest of the parameters and the current model configuration. We start the MCMC from the simplest model $(2,1,1,1)$ and choose as starting value of $(a,b)$ the MLE using a model of the form $\text{log}(m_{x,t}) = a + b (x - \bar {x})$ , and $k^1_t = 0$ .

3.1.1 Parameters update

For the parameter $\theta _i$ , we sample from the posterior distribution of the parameters with a Metropolis-Hastings algorithm using as proposal the Laplace approximation of the posterior distribution. To find the Laplace approximation, we maximize the log-posterior distribution $\text{log}(p(\theta _i | \theta _{-i}, d))$ to find the optimal $\hat {\theta _i}$ and the Hessian $\hat {H}$ at the maximum. We use an explicit expression of the solution $\hat {\theta _i}$ and of the Hessian $\hat {H}$ when analytically available, and resort to using an optimization procedure such as the optim function in R when analytical expressions are not available. To achieve better performances, we also provide optim with the explicit expression of the gradient $\nabla _{\theta _i} \text{log}(p(\theta _i | \theta _{-i}, d))$ . We have used a Laplace approximation using a Student’s t distribution in place of a normal distribution, since acceptance rate using a normal distribution can in some cases be too small when the current points is very distant from the mean of the Laplace approximation, because of the light tails of the normal. The Student’s t approximation has been chosen by matching the mean and the scale matrix and using $\nu = 10$ degrees of freedom. The proposed value $\theta _i^*$ from $t_{\nu }(\hat {\theta _i}, \hat {H})$ is accepted with probability

\begin{equation*} \min \left \{1, \frac {p(d | \theta _i^*, \theta _{-i}) t_{\nu }(\theta _i | \hat {\theta _i}, \hat {H})}{p(d | \theta _i, \theta _{-i}) t_{\nu }(\theta _i^{\star } | \hat {\theta _i}, \hat {H})}\right \}. \end{equation*}

If constraints are not present on the parameters, we update each parameter individually; otherwise, we update parameters in blocks. To handle constraints in the updates, we reparametrize each set of parameters in terms of a set of free variables, ensuring that all constraints are automatically satisfied. We detail how constraints are enforced in each model in Section $1$ of the Supplementary Material.

3.1.2 Model space update

To explore the model space, we propose to change each variable $\delta _i$ at a time. For $\delta _1$ , $\delta _3$ , and $\delta _4$ , we simply propose the value $2$ if we are currently in the value $1$ and vice versa. For $\delta _2$ , we propose the following transitions: $(\delta _2 = 1) \Rightarrow (\delta _2^{\star } = 2, \delta _3^{\star } = 1)$ , $(\delta _2 = 2) \Rightarrow (\delta _2^{\star } = 1)$ , $(\delta _2 = 2) \Rightarrow (\delta _2^{\star } = 3)$ , and $(\delta _2 = 3) \Rightarrow (\delta _2^{\star } = 2)$ . If more than one move is available from the current configuration, we choose each of them with equal probability.

To formally describe the updates, we set up some notation. We assume that the current model $\delta$ has parameter $(\eta _1,\eta _2)$ and the new proposed model $\delta ^{\star }$ has parameter $(\eta _1,\eta _2^{\star })$ (noting that $\eta _2$ or $\eta _2^{\star }$ can also be empty). For example, if the current model is nested within the proposed model, we have $\eta _2 = \emptyset$ .

The moves we will consider fall into one of the following cases:

1. 1. $\eta _2 = \emptyset$ , $\eta _2^* \neq \emptyset$ . This is the case where we are proposing to add new parameters to the model. In this case, we can simply choose as auxiliary variables $u$ the new set of parameters $\eta _2^{\star }$ . Specifically, we generate $u = \eta _2^{\star }$ from the normal Laplace approximation to the posterior $p(\eta _2 | \eta _1, \delta ^*, d)$ , and accept the proposed model with the new set of parameters, $((\eta _1,\eta _2^*),\delta ^*)$ with probability

   \begin{equation*} \text{min}\left \{ 1, \frac {p(d | (\eta _1, \eta _2^*),\delta ^*) p(\eta _1, \eta _2 | \delta ^*) p(\delta ^*) q(\delta ^* \rightarrow \delta ) }{p(d | \eta _1, \delta ) p(\eta _1 | \delta ) p(\delta ) q(\delta \rightarrow \delta ^*) \text{N}(\eta _2^* | \hat {\eta _2}, \hat {H})} \right \}. \end{equation*}

2. 2. $\eta _2^* = \emptyset$ , $\eta _2 \neq \emptyset$ . This is the case where we are proposing to remove parameters of the model. In this case, the acceptance ratio is the inverse of the ratio in the previous move.

3. 3. $\eta _2 \neq \eta _2^*$ both not empty. This is the case where we are proposing to change the parameters of the current model. In this case, we choose as auxiliary variables $u = \eta _2^*$ and $u' = \eta _2$ . The new model is accepted with probability as in (1).

We will need to pay particular attention to the cases where the move will cause changes in the ICs on the set of parameters that do not change ( $\eta _1$ ). When this happens, we cannot straightforwardly use the updates defined above because the dimension matching condition is not satisfied anymore. There are three possible ways to work around this. One option is to propose jointly the new set of parameters $\eta _2^*$ and the subset of parameters of $\eta _1$ for which the ICs change. The drawback is that this would require an optimization over a large number of parameters, which is costly. Another option is to simply impose the ICs also for the model where these are not theoretically needed. The preferred option is to reparametrize the model such that the reduction in ICs is replaced with the addition of new variables, which are then proposed together with the variables $\eta ^*$ . We can now describe each different update of $\delta _i$ :

• • Update $\delta _1$ To propose the move from $\delta _1$ to $\delta ^*_1$ , we use the third type of update by choosing as auxiliary variables $u = a_x$ (if $\delta ^{*}_1 = 1$ ) and $u = (a,b)$ (if $\delta ^{*}_1 = 2$ ). In the case $\delta _2 = 2$ or $\delta _2 = 3$ , changing $\delta _1$ causes a change in the ICs on $k^2_t$ , since $\sum _t k^2_t = 0$ is required when $\delta _1 = 1$ but it is not in the case $\delta _1 = 2$ . We therefore choose to keep the IC $\sum _t k^2_t = 0$ throughout to avoid complicating the RJMCMC. In the case where $(\delta _2, \delta _3, \delta _4) = (3,1,2)$ , the move from $\delta _1 = 1$ to $\delta _2 = 2$ causes the constraints $\sum c^2 \gamma _c = 0$ on $\gamma _c$ to drop. As mentioned above, we aim to reparametrize the model with fewer constraint ( $\delta _1 = 1$ ) such that the new constraint $\sum c^2 \gamma _c = 0$ is added (and therefore, a degree of freedom is lost) by adding an additional variable into the model. To find the parametrization, we note that the additional constraint $\sum c^2 \gamma _c = 0$ removes quadratic trend from $\gamma _{t-x}$ , hence, if this constraint is added, we need to introduce a quadratic term with zero mean and no linear trend into the model. Therefore, we propose to reparametrize the model

  \begin{equation*} a + b(x - \bar {x}) + k^1_t + k^2_t \bar {b} (x - \bar {x}) + \gamma _{t-x} \qquad \end{equation*}
  with constraint $\sum \gamma _c = \sum c \gamma _c = 0$ , to the form
  \begin{equation*} a + b(x - \bar {x}) + k^1_t + k^2_t \bar {b} (x - \bar {x}) + h_{\zeta }(t-x) + \gamma _{t-x} \qquad \end{equation*}
  with constraint $\sum \gamma _c = \sum c \gamma _c = \sum c^2 \gamma _c = 0$ , where $h_{\zeta }(c) = \zeta \frac {(n_C + 1)(n_C + 2)}{6}- \zeta n_C (c + n_X) + \zeta (c + n_X)^2$ , and $\zeta$ is a free parameter. $h_{\zeta }(c)$ is a quadratic polynomial in $c$ that follows that specific parametrization to satisfy the existing constraints $\sum _c h_{\zeta }(c) = \sum _c c h_{\zeta }(c) = 0$ . The degree of freedom lost because of the constraint has been absorbed in the parameter $\zeta$ . Having the two models the same constraints, we can propose $\delta ^*_1 = 2$ from $\delta _1 = 1$ by proposing $u = (a,b, \eta )$ and in the opposite move, we propose $u = a_x$ .

• • Update $\delta _2$ We separately consider each move:

  1. – $(\delta _2 = 1) \leftrightarrow (\delta _2 = 2, \delta _3 = 1)$ :

     To perform the move $(\delta _2 = 1) \rightarrow (\delta ^*_2 = 2, \delta ^*_3 = 1)$ , we use the first type of update by proposing $\bar {b}$ . The opposite move is performed deterministically, by setting $\bar {b} = 0$ .

  2. – $(\delta _2 = 2) \leftrightarrow (\delta _2 = 3)$

     In this update, we need to pay particular attention to the cases where the ICs of other parameters change as a consequence of the change in $\delta _2$ . As it can be seen in Table 1, this happens where $(\delta _1, \delta _3, \delta _4) = (1,1,2)$ (since $\sum c^2 \gamma _c = 0$ applies only for $\delta _2 = 3$ ), and when $\delta _3 = 2$ (since $b_{x_1} = -1$ applies only for $\delta _2 = 3$ ). We discuss each case separately.

     1. * $\delta _3 = 1$ and, $\delta _1 = 2$ or $\delta _4 = 1$ . This is the case in which no ICs change on the rest of the parameters. In this case, we use the third update by using $u = k^2_t$ when moving to $\delta ^*_2 = 2$ and $u = (k^1_t, k^2_t)$ when moving to $\delta ^*_2 = 3$ .

     2. * $\delta _3 = 2$ . In this case, the ICs on $b_x$ change. We propose to reparameterize the model for $\delta _2 = 2$ :

        \begin{equation*} f_1(x) + k^1_t + k^1_t b_x + f_4(t-x) \qquad \end{equation*}
        in the form
        \begin{equation*} f_1(x) + k^1_t + k^1_t \bar {b} b_x + f_4(t-x) \qquad \end{equation*}
        with the additional constraint $b_{x_1} = -1$ . The degree of freedom lost because of the constraint has been absorbed in the parameter $\bar {b}$ . Now that the model has the same constraints as the model of the case $\delta _2 = 3$
        \begin{equation*} f_1(x) + k^1_t + k^2_t b_x + f_4(t-x) \qquad \end{equation*}
        we can move to $\delta ^*_2 = 3$ by proposing $u = (k^1_t, k^2_t)$ and we move to $\delta ^*_2 = 2$ by proposing $u = (k^2_t, \bar {b})$ .

     3. * $(\delta _1, \delta _3, \delta _4) = (1,1,2)$ . In this case, the constraint $\sum c^2 \gamma _c = 0$ does not hold for $\delta _2 = 2$ . Similarly to what we did in the $\delta _1$ case, we propose to reparametrize the model for $\delta _2 = 2$ :

        \begin{equation*} a_x + k^2_t + k^2_t \bar {b} (x - \bar {x}) + \gamma _{t-x} \qquad \end{equation*}
        with constraint $\sum \gamma _c = \sum c \gamma _c = 0$ , to
        \begin{equation*} a_x + k^2_t + k^2_t \bar {b} (x - \bar {x}) + h_{\zeta }(t-x) + \gamma _{t-x} \qquad \end{equation*}
        with constraints $\sum \gamma _c = \sum c \gamma _c = \sum c^2 \gamma _c = 0$ , Having the two models the same constraints, we can propose $\delta ^*_2 = 3$ from $\delta _2 = 2$ by proposing $u = (k^1_t,k^2_t)$ and in the opposite move we can propose $u = (k^2_t, \bar {b}, \zeta )$ .

• • Update $\delta _3$ We distinguish two cases.

  1. – $\delta _2 = 2$ : in this case we are moving between the models $k^2_t (1 + \bar {b}(x- \bar {x}))$ and $f_2(t,f_3(x)) = k^2_t (1 + b_x)$ . We use the third type of update, using the auxiliary variables $u = b_x$ in the move from $\delta _3 = 1$ to $\delta _3^* = 2$ and $u = \bar {b}$ in the reverse move.

  2. – $\delta _2 = 3$ : in this case we are moving between the models $f_2(t,f_3(x)) = k_t^1 + k^2_t (x - \bar {x})$ and $k^1_t + k^2_t b_x$ . Since one model is nested within the other, we can use the first type of update with $u = \bar {b}_x$ when moving from $\delta _3 = 1$ to $\delta _3^* = 2$ and the second type of update in the reverse move.

     We separately consider the case $(\delta _1, \delta _2, \delta _4) = (1,3,2)$ , since the change from $\delta _3 = 1$ to $\delta ^*_3 = 2$ causes the constraint $\sum c^2 \gamma _c = 0$ to drop.

     Similarly to before, we reparametrize the model $a_x + k^1_t + k^2_t (x - \bar {x}) + \gamma _{t-x}$ , with constraints $\sum \gamma _c = \sum c \gamma _c = 0$ to the model $a_x + k^1_t + k^2_t (x - \bar {x}) + h_{\zeta }(t - x) + \gamma _{t-x}$ , with constraints $\sum \gamma _c = \sum c \gamma _c = \sum c^2 \gamma _c = 0$ . Next, to change between the two models, we use the first type of update with $u = (\zeta , b_x)$ . The reverse move is deterministic.

• • Update $\delta _4$ The move from $\delta _4 = 1$ to $\delta _4^* = 2$ falls into the first scheme and is performed by choosing $u = \gamma _{c}$ , while the move from $\delta _4 = 2$ to $\delta _4^* = 1$ falls into the second scheme and is a deterministic move.

We note that even in cases where models are nested, we have not used the first/second type of update as in our experiments, the acceptance ratio was too low. The scheme is summarized in Figure 1.

Figure 1

Summary of the model moves $\delta \rightarrow \delta ^{\star }$ . We use $\tilde {m}_{x,t}$ to denote all the residual terms.

3.2 Simulation study

We have performed a simulation study to assess the ability of the model selection procedure to identify the true model. For each of the $20$ models, we have simulated $4$ different datasets and for each dataset, and we run the model with $2$ chains of $5000$ burn-in iterations and $5000$ iterations. We have used flat priors on all the model parameters. We have used the following settings to perform the simulations:

• • $f_1(x)$ : if $\delta _1 = 2$ , we set $f_1(x) = a + b (x - \bar {x})$ , with $a = -3.5, b = 0.105$ , if $\delta _1 = 1$ , we set $f_1(x) = a + b (x - \bar {x}) + \epsilon _x$ , with $a,b$ as before and $\epsilon _x \sim \text{N}(0, 0.3^2)$ .

• • $f_2(t)$ , If $\delta _2 = 1$ or $\delta _2 = 2$ , we let $k^1_t$ (or $k^2_t$ ) vary linearly from $0.675$ to $-0.675$ . If $\delta _2 = 3$ , we let $k^1_t$ vary on a grid from $-0.675$ to $0.675$ and assume $k^2_t = k^1_t + \epsilon _t$ , with $\epsilon _t \sim \text{N}(0, 0.3^2)$ (and then centering $\epsilon _t$ ).

• • $f_3(x)$ . if $\delta _3 = 1$ , we assume $\bar {b} = 0.0818$ , if $\delta _3 = 2$ , we set $b_x = \bar {b} (x - \bar {x}) + \epsilon _x$ , with $\epsilon _x \sim \text{N}(0, 0.3^2)$ . If $b_{x_1} = -1$ , we have renormalized to account for it.

• • $f_4(t-x)$ : if $\delta _4 = 2$ , we have simulated $\gamma _{t-x} = A \eta$ , where $\eta \sim N(0, 0.3^2)$ is a $(n_C-2)$ or $(n_C-3)$ dimensional vector, and $A$ is a $n_C \times (n_C-2)$ (or $n_C \times (n_C-3)$ ) projection matrix such that $\gamma _{t-x}$ will satisfy the relevant constraints.

The previous values have been chosen to reflect values typically obtained from real data. The results are presented in Table 2. The inference algorithm can generally recover the true model, except for a few instances. When the true model was not selected, very similar models were selected instead. An issue we have noticed is that the model $(1,3,1,2)$ is often selected even when it is not the correct model. In our experiments, we noticed that this was due to issues with slow mixing related to this model, since when $b_x \approx k (x - \bar {x})$ the model approaches unidentifiability (as another constraint on the cohort effect would be needed).

Table 2.

Results of the simulation study. In each cell, we report the posterior probability of each model averaged across the $4$ simulations. The true model used to simulate the data is indicated by the row, while the model for which the posterior probability is reported is indicated by the column. The index of the model is the same as the order used in Table 1

Figure 2

HMD data results: 95% PCI of the log-mortality rates $m_{x,t}$ . The dots represent the crude rates.

3.3 Case study: Human mortality database

We apply the models defined in the previous section to mortality data of $40$ countries as part of the Human Mortality Database (HMD, 2024). For each country, we selected a subset of data comprising individuals aged between $60$ and $90$ , and covering the period from $1990$ to $2022$ . We note that some data are available on a smaller timescale, and in that case, we use the subset of data available within that period. For each country, we ran the MCMC for $2$ chains of $10,000$ iterations, with $20,000$ burn-in iterations. We have used flat priors on all the model parameters. The posterior credible intervals (PCI) of the estimated log-mortality rates are presented in Figure 2. For simplicity, we only show here the results for a few countries (the full results are presented in the Supplementary Material). We also present the posterior probabilities of the model selection procedure in Table 3. In many countries, it can be seen that the evidence for the selection of one term is strong enough that the posterior probability is $1$ , while for other countries, multiple models are selected. The countries with the largest exposure, such as USA, Russia, Japan, Germany, and the UK, select the most complex model, $(1,3,2,2)$ , $a_x + k^1_t + k^2_t b_x + \gamma _{t-x}$ . Countries with lower exposures, such as Croatia, Norway, and Slovakia are the ones that result in higher uncertainty in the model selection. The rest of the mortality rates are presented in Section  $2.1$ of the Supplementary Material, while the MCMC mixing summaries are reported in Section  $2.2$ of the Supplementary Material. To evaluate parameter mixing, we have reported summaries of the effective sample size and the Gelman–Rubin diagnostics. To evaluate the model choice mixing, we have reported two quantities. First, we have reported the posterior model probabilities across different chains. Excessive variation between these two probabilities shows poor mixing. In addition, we have looked at adjusted acceptance rates toward different model configurations. The rationale behind the proposed efficiency rate is the following. If the MCMC is perfectly efficient, we expect the probability of transition from model $i$ to model $j$ , $p(\delta ^{(t)} = j | \delta ^{(t-1)} = i)$ , to be close to the posterior probability of model $j$ , $p(\delta = j| d)$ . Therefore, can look at the ratio $\frac {p(\delta ^{(t)} = j | \delta ^{(t-1)} = i)}{p(\delta = j| d)}$ as a measure of efficiency, and we term this as “efficiency rate” in the Supplementary Material. Excessive variation from $1$ can be seen as a measure of poor mixing in the RJMCMC. As shown in the summaries, aside from a few countries such as Russia, Japan, and Germany, parameter mixing is generally acceptable, and even when individual model terms do not exhibit adequate effective sample sizes, the mortality rate tends to mix well. This suggests that the slow mixing is likely due to autocorrelation in the MCMC between the different model terms. Model mixing can only be evaluated for countries that select more than one model. For some of these countries, such as Canada and Slovakia, the diagnostics indicate that model mixing is less satisfactory. Although it is unrealistic to expect efficiency rates close to $1$ , our results show that, in some cases, the rates are close to $0$ , highlighting how the RJMCMC can struggle to transition between models. We note that, to the best of our knowledge, the efficiency rate has not been used in the literature, so there are no standard benchmarks to contextualize the efficiency of our procedure. However, we can conclude that efficiency rates below 1% are strong evidence of poor RJMCMC performance. The within-chain model probabilities support a similar conclusion, as the model proportions can for some countries differ noticeably between chains.

Table 3.

HMD data results: posterior probabilities of each model $(\delta _1,\delta _2,\delta _3,\delta _4)$

4. APCI model stratified by product

In this section, we demonstrate how to define mortality models for data stratified by an additional variable, that is, the APP data. For example, the additional stratifying variable could represent different insurance products, sexes, or countries. For simplicity, we will refer to this additional stratifying variable as “product” throughout this section. Note that the notation used here is different from the one used in the previous section.

We build the model by extending the APCI model (Richards et al., Reference Richards, Currie, Kleinow and Ritchie2019). As mentioned in the introduction, we have chosen this base model for two main reasons. First, the APCI model is one of the most popular models used to fit mortality data and is currently used by the CMI (Continuous Mortality Investigation Bureau, 2024). Second, we aim to define an extension of the base model that involves fewer changes in the ICs when each model changes according to product, which, similar to the previous section, will help the design of the RJMCMC framework, and the APCI serves this purpose well.

According to the APCI model, log-mortality is modeled as

\begin{equation*} \text{log}(m_{x,t}) = a_x + b_x (t - \bar {t}) + k_t + \gamma _{t-x} \ , \end{equation*}

where $a_x$ represents the age effect, $b_x$ represents the change in period effect by age, $k_t$ represents the period effect and $\gamma _{t-x}$ represents the cohort effect. The ICs that we use are $\sum _x b_x = \sum _t k_t = \sum _{c} \gamma _{c} = \sum _c c \gamma _c = \sum _c c^2 \gamma _c = 0$ . We then extend the base APCI model to allow for product effect by proposing different functional forms for the terms $a_x$ , $b_x$ and $k_t$ .

• • For $a_x$ , we propose the product forms: (a) $a_x$ (no change), (b) $a_x + c^1_p$ , and (c) $a_{x,p}$ , and we denote each choice with the latent variable $\delta _1 = 1$ , $2$ and $3$ , respectively. The ICs are $\sum c^1_p = 0$ .

• • For $b_x$ , we propose the product forms: (a) $b_x$ (no change) and (b) $b_x c^2_p$ , modeled with the latent variable $\delta _2 = 1$ and $2$ , respectively. The ICs are $\sum c^2_p = n_P$ , since $c^2_p$ appears in a product with $b_x$ .

• • For $k_t$ , we propose the forms: (a) $k_t$ (no change) and (b) $k_t + k^2_t c^3_p$ , modeled with the latent variable $\delta _3 = 1$ and $2$ , respectively. The ICs are the following. $\sum k^2_t = 0$ , since any constant $K c^3_p$ can be incorporated into $ c^1_p$ (or into $a_{x,p}$ ), $\sum c^3_p = 0$ since any additive constant can be incorporated into $k_t$ . $c^3_1 = 1$ , Since $c_p$ can still be multiplied by a constant which can be incorporated into $k^2_t$ .

Similarly to the previous section, we note that some constraints depend on the modeling choices. The third constraint $\sum c^2 \gamma _c = 0$ is not necessary when $\delta _2 = 2$ . $\sum k^2_t = 0$ only holds if $\delta _1 \neq 1$ , and $\sum _{x} b_x = 0$ only holds in the case $\delta _2 = 1$ . We note that we still assume $\sum _x b_x = 0$ also in the case $\delta _2 = 2$ since, in our experiments, the model was very close to being unidentifiable without the constraint $\sum _x b_x = 0$ . Moreover, similarly to the previous section, to help the design of the RJMCMC, we choose to still impose the constraint $\sum k^2_t = 0$ even when $\delta _1 = 1$ .

The model for age-period-product data can be summarized as

\begin{equation*} \text{log}(m_{x,t,p}) = a_{\delta _1}(x,p) + b_{\delta _2}(x,p) (t - \bar {t}) + k_{\delta _3}(t,p) + \gamma _{t- x} \end{equation*}

where

\begin{equation*}a_{\delta _1}(x,p) = \begin{cases} a_x \hspace {.9cm} \text{ if } \delta _1 = 1 \\[5pt] a_x + c^1_p \hspace {.2cm}\text{ if } \delta _1 = 2 \\[5pt] a_{x,p} \hspace {.7cm} \text{ if } \delta _1 = 3 \end{cases}\!\!\!\!\!, b_{\delta _2}(x,p) = \begin{cases} b_{x} \hspace {.5cm} \text{ if } \delta _2 = 1 \\[5pt] b_x c^2_p \hspace {.2cm }\text{ if } \delta _2 = 2 \\[5pt] \end{cases}\!\!\!\!\!, k_{\delta _3}(t,p) = \begin{cases} k_{t} \hspace {1.3cm} \text{ if } \delta _3 = 1 \\[5pt] k_t + k^2_t c^3_p \hspace {.3cm} \text{ if } \delta _3 = 2 \\[5pt] \end{cases}\!\!\!\!\!. \end{equation*}

The full set of models is presented in Table 4.

Table 4.

Complete set of models for the APP data

4.2 Case study: Insurance products

We apply the previously described framework on mortality data stratified by different products. The data were sourced by the CMI and were provided by different life insurance providers. The available products are pension annuities in payment (Annuities), term assurances (TA), accelerated critical illness (ACI), and standalone critical illness (SCI). The data range from $2016$ to $2020$ and cover both males and females and across varying commencement years, durations, and sum assured bands. In our example, we have used ages from $50$ to $80$ . We note that the data come with several limitations, such as very low exposure for some combination of age, year, and product as shown in Figure 3. However, working in a Bayesian framework naturally allows the model to take the low sample size into account in the estimation. For more details on the data, we refer the reader to the CMI Working Papers 117 and 162 (Continuous Mortality Investigation Bureau, 2019, 2022). We have run $2$ chains of $50,000$ iterations with $5000$ burn-in iterations. Because of the lack of data, we have used priors $a_{x,p} \sim \text{N}(-3,3^2)$ for $\delta _1 = 3$ . We have used flat priors elsewhere.

Figure 3

Insurance product case study: log(Exposures).

The credible intervals of the log-mortality curves are shown in Figure 4. It is clear that individuals associated with term assurance products exhibit lower mortality than those associated with other products. Surprisingly, mortality rates for the TA product are lower than the ones for the Annuities product, which is counter-intuitive since individuals signing up to the Annuities product are generally expected to live longer. For older ages, the conclusion for the products ACI and SCI is less clear since the variability of their estimates considerably increases in older ages, which results from the very low exposure as visible in Figure 3.

Figure 4

Insurance product case study: $95$ % PCI of the log-mortality rates. The crude mortality rates are represented by the dots. Data points where the crude mortality rates are $0$ are represented by a triangle.

The posterior probabilities of the model choice variables are shown in Table 5. The most supported model has the form

\begin{equation*}a_{x,p} + b_x (t - \bar {t}) + (k_t + k^2_t c^3_p) + \gamma _{t-x}\end{equation*}

with a posterior probability of $0.99871$ , followed by the model

\begin{equation*}a_{x,p} + b_x c^2_p (t - \bar {t}) + (k_t + k^2_t c^3_p) + \gamma _{t-x}\end{equation*}

with a lower posterior probability of $0.00129$ . We note that this result highlights the advantage of using our model selection framework, where the preferred model shares (with high probability) the term $b_x (t- \bar {t})$ across all the different products. By contrast, fitting the model separately for each product would involve specifying a separate term $b_{x,p} (t - \bar {t})$ for each $p$ , which may not be the most ideal choice (as in the case here).

We also present in Figure 5 the posterior distributions for the parameters of the model. In Figure 5c, we report the PCI of $k(t,p)$ . It can be seen that the period effect is generally negative, since mortality tends to decrease with time, but for some products, the trend is not always present, potentially as a result of the small sample size (only $5$ years of data). In Figure 5b, we present the PCIs of $b_x$ , which represents the interaction between age and period effect, after accounting for age and period effect individually. Lower values of $b_x$ represent larger improvement in mortality as the period increases, and vice versa. In our results, it can be seen that for older ages, the improvement in mortality is lower. The result seems to suggest that ages between $60$ and $75$ are the ones who are likely to experience the largest improvement in mortality in the long run. We note that, due to the different parametrization used, the interpretation of $b_x$ here differs from the interpretation under the standard constraints. However, we suggest analyzing the results with care because of the very small number of years used (only $5$ years). Some MCMC diagnostics are reported in Section $3$ of the Supplementary Material. Similarly to the previous case study, parameter mixing is satisfactory for the mortality rate $m_{x,t,p}$ , while mixing worsens for each model term. Results for the model mixing are encouraging, although we noticed one model is selected very rarely, making it difficult to properly assess the efficiency of the RJMCMC.

Table 5.

Product data: posterior probabilities on the model choice variables $\delta _i$