2.1 Stacking sub-populations

For joint modelling purposes, datasets are stacked together. Generically, we have a total of L populations, indexed by the subscript $l=1,2,\ldots, L$ . When needed, to indicate country, cause, and gender factors, we re-index by $l=(c,s,g)$ . Throughout, the two main independent variables are age and year, $(x^i_{ag},x^i_{yr})$ , and the observed mortality rate in the l-th population is denoted by:

(1) \begin{equation} \mu^i_{l} \;:\!=\; \dfrac{\text{Death counts at }(x^i_{ag},x^i_{yr})\text{ of type }l}{\text{Exposed-to-risk counts at }(x^i_{ag},x^i_{yr})}= \dfrac{D^i_{l}}{E^i_{l}},\end{equation}

where $D^i_{l}$ is the number of deaths in the lth population for the ith observation with the corresponding age and year inputs. The denominator $E^i_{l}$ refers to the corresponding exposed-to-risk counts. Some of the $E^i_l$ ’s will be the same when we consider different causes: for a fixed country–gender combination (c, g) $E^i_{(c,s,g)}$ does not depend on s while $D^i_{(c,s,g)}$ does.

Our GP models work with log-mortality rates so that our data is

\begin{align*} y^i_{l} = \log \mu^i_{l}.\end{align*}

The overall dataset is then represented as $(x^i, y^i_l)$ , $i=1,\ldots,N_l$ , $l=1,\ldots, L$ .

3.1 Gaussian process regression for mortality tables

Consider the non-parametric additive regression model for the L populations to be jointly analysed:

(2) \begin{equation} y^i_l = f_l(x^i) + \epsilon^i_l, \qquad\qquad i=1,\ldots, N,\end{equation}

where $x^i$ represents an individual entry in the mortality table (indexed by age and year), $y^i_l$ is the observed output in the lth population, and $f_l(\!\cdot\!)$ , $l=1,\ldots,L$ is the underlying latent log-mortality surface, obscured with the observation noise $\epsilon^i_l$ .

We first summarise the spatial structure that concerns the dependence of mortality rates as a function of age and year. Momentarily focusing on a single-output Gaussian process (SOGP) regression, we put a GP prior on the latent function $f_l \sim \mathcal{GP}(m,C)$ , meaning that any finite vector $f_l(\textbf{x})=(f_l(x^1),\ldots,f_l(x^n))$ at n inputs follows the multivariate Gaussian distribution:

\begin{equation*} f_l(x^1),\ldots ,f_l(x^n) \sim \mathcal{N}\big(\textbf{m}_{\textbf{l}}{(\textbf{x})},\textbf{C}_l(\textbf{x,x})\big),\end{equation*}

where ${m}_l(x)=\mathbb{E}[f_l(x)]$ is the mean vector of size n and ${C_l(x,x')}=\mathbb{E}[(f_l({x})-{m_l(x)})(f_l({x'})-{m_l(x')})]$ is the n-by-n covariance matrix. All properties of a GP are thus completely described by its mean and covariance functions. The imposition of this Gaussian structure is purely for machine learning purposes in order to leverage the extensive theory of reproducing kernel Hilbert spaces and focus the modelling efforts on the covariance kernel as a way to describe the dependence across different mortality rates.

The functions $m_l(\!\cdot\!)$ and $C_l(\cdot, \cdot)$ characterise our prior beliefs about the response surface $f_l$ . The covariance kernel of the GP defines a similarity between pairs of data points. It characterises the smoothing process by determining the influence of observations on the distribution of the output. Data points that are close are expected to behave more similarly than data points that are farther away. In terms of spatial dependence on age and year, we concentrate on a common family of covariance functions known as the Matérn class, equipped with automatic relevance determination. Specifically, the Matérn-5/2 kernel defines the covariance between two mortality table entries x, x ^′ as follows:

(3) \begin{equation} C_l(x,x') = \prod_{k \in \{ag,yr\}} \bigg(1+\dfrac{\sqrt{5}}{\theta_{k,l}}|x_k-x'_{\!\!k}|+\dfrac{5}{3\theta_{k,l}^2}|x_k-x'_{\!\!k}|^2\bigg) \exp\bigg(-\dfrac{\sqrt{5}}{\theta_{k,l}}|x_k-x'_{\!\!k}|\bigg).\end{equation}

This kernel is parametrised by the age lengthscale $\theta_{ag,l}$ and the year lengthscale $\theta_{yr,l}$ . Our choice of (3) is driven by the popularity of Matérn-5/2 in the GP literature, which produces twice-differentiable predictive surfaces (important for stably evaluating mortality improvement factors) but is more flexible than an infinitely differentiable kernel. See Huynh & Ludkovski (Reference Huynh and Ludkovski2021) for further discussion of kernel choice. In short, this is a separate, non-trivial task that has some, but not drastic, effects on the results.

The mean function $m_l(x)$ describes the relevant trends in log-mortality rates. Typically in GP models, the mean is a constant; however, this is not an appropriate choice given the strong age dependence of mortality. To capture the long-term longevity features, we thus fit a parametric trend: $m_l(x)=\beta_{0,l}+\sum_{j=1}^p \beta_{j,l} h_j(x),$ where $h_j$ ’s are given basis functions and the $\beta_{j,l}$ ’s are unknown coefficients. Motivated by Figure A.2 in Appendix A.2 which shows that cause-specific log-mortality tends to exhibit a concave, rather than a linear pattern in age, we consider a quadratic trend in age and linear trend in the year dimension:

(4) \begin{equation} m_l(x^i) = \beta_{0,l} + \beta^{ag}_{1,l}x^i_{ag} + \beta^{ag}_{2,l}(x^i_{ag})^2 + {\beta^{yr}_{l,1}} x^i_{yr}.\end{equation}

We are interested in the posterior distribution for $\textbf{f}_{\ast} \equiv f_l(\textbf{x}_{\ast})$ at specified inputs $\textbf{x}_{\ast}$ given the dataset $\textbf{y}_l = (y^{1}_l, \ldots, y^N_l)$ , $p(\textbf{f}_{\ast}|\textbf{y}_{\textbf{l}})$ , in other words the likelihood of the true response surface being $\textbf{f}_{\ast}$ given what we have observed. Using ${\textrm{Cov}}(y^i_l,y^j_l)={\textrm{Cov}}(f_l(x^i),f_l(x^j))+\sigma^2_l\delta(x^i,x^j)$ where $\delta(x^i,x^j)$ is the Kronecker delta, we have the Gaussian observation likelihood $\textbf{y}_{\textbf{l}} \sim \mathcal{N}(\textbf{m}_{\textbf{l}}\textbf{(x)},\textbf{C}_{\textbf{l}}\textbf{(x,x)}+\boldsymbol{\Sigma}_{l})$ where the error terms are assumed to be independent and Gaussian-distributed, $(\epsilon^1_l,\ldots,\epsilon^N_l) \sim \mathcal{N}(\textbf{0}, \boldsymbol{\Sigma}_l = \text{diag}(\sigma_l^2))$ . Incorporating the evidence from observations, as reflected in the likelihood function and the prior, we obtain the Gaussian posterior:

\begin{equation*} p(\textbf{f}_{\ast}|\textbf{y}_l) \sim \mathcal{N}(\textbf{m}_{\ast}(\textbf{x}_{\ast},\textbf{l}),\textbf{C}_{\ast,\textbf{l}}(\textbf{x}_{\ast},\textbf{x}_{\ast})).\end{equation*}

The Universal Kriging equations (5)–(7) (Rasmussen & Williams, Reference Rasmussen and Williams2005, Ch 2.7) below provide not only the posterior mean $m_*(\cdot,l)$ and posterior variance $s^2_*(\cdot,l)$ but also the estimated coefficients $\boldsymbol{\beta}_l=(\beta_{1,l},\ldots,\beta_{p,l})^T$ . Let $\textbf{h}(x)=\big(h_1(x), \ldots,h_p(x)\big)$ , $\textbf{H}=\big(\textbf{h}(x^1),\ldots,\textbf{h}(x^N)\big)$ and $\textbf{D}=(\textbf{C}+\boldsymbol\Sigma)^{-\textbf{1}}\textbf{H}$ where $\textbf{C}$ is the covariance matrix $C_l(x^i,x^j)_{i,j=1}^N$ . The posterior mean of $\boldsymbol{\beta}_l$ along with the predicted posterior mean $m_*(x_*,l)$ and respective variance $s_*^2(x_*,l)=C_{*,l}(x_*,x_*)$ for any input $x_*$ are as follows:

(5) \begin{align} {\hat{\boldsymbol\beta}}_l =(\textbf{H}^T\textbf{D})^{-1}\textbf{H}^T(\textbf{C}+\boldsymbol\Sigma)^{-1}\textbf{y}; \end{align}

(6) \begin{align} m_*(x_*,l) = \textbf{h}(x_*)^T {\hat{\boldsymbol\beta}}+\textbf{c}(x_*)^{T}(\textbf{C}+\boldsymbol\Sigma)^{-1}(\textbf{y}-\textbf{H}{\hat{\boldsymbol\beta}}); \end{align}

(7) \begin{align} s^2_*(x_*,l) &\ = (\textbf{h}(x_*)^T -\textbf{c}(x_*)^T\textbf{D})^T(\textbf{H}^T\textbf{D})^{-1} (\textbf{h}(x_*)^T-\textbf{c}(x_*)^T\textbf{D})\end{align}

where $\textbf{c}(x_*)=(C_l(x^1,x_*), \ldots, C_l(x^N, x_*))$ is the vector of covariances between inputs in the training set and desired test input $x_*$ . Note that the predictive distribution of observation ${y_l}$ at $x_*$ is similarly obtained as ${y_l} \sim \mathcal{N}(m_{*}(x_*,l),s^2_*(x_*,l)+\sigma^2_l)$ .

Below we use $m_*(x_*,l)$ as the model prediction for the respective (log)-mortality rate of the lth population in cell $x_*$ , and $s_*(x_*, l)$ as the corresponding posterior uncertainty which is used to obtain predictive quantiles around the former prediction.

3.2 Semi-parametric latent factor model

The vector-valued latent response variable over the age–year input space is defined as $\textbf{f(x)} = (f_1(\textbf{x}),\ldots,f_L(\textbf{x}))$ , where the functions $\{f_l({x})\}_{l=1}^L$ are the log-mortality surfaces for the corresponding lth population. Similar to single-population GP above, we place a GP prior over the latent function $\textbf{f}$ such that:

\begin{align*} \textbf{f} \sim \mathcal{GP}(\textbf{m, C}),\end{align*}

where $\textbf{m}$ is the mean vector function whose elements $\{m_l({x})\}_{l=1}^L$ are the mean functions of each output and $\textbf{C}$ is the fused (meaning combining both spatial and cross-population terms) covariance matrix. This implies that we separate the “spatial” dependence, encoded in x from the inter-task dependence encoded in l. While x has a natural Euclidean metric that is used in (3), the population indices l are generally of factor type. Hence, to describe the dependence across L outputs, we need $L(L-1)/2$ hyperparameters $C^{(f)}_{lk}, 1 \le l,k \le L$ which becomes inefficient and unstable beyond 3-4 populations. Thus, with an eye towards reducing the number of hyperparameters, additional structure is needed for $C^{(f)}$ .

In the SLFM dating back to Teh et al. (Reference Teh, Seeger and Jordan2005), each output $f_l({x})$ is assumed to be a linear combination of Q latent functions:

(8) \begin{equation} f_l({x}) = \sum_{q=1}^Q a_{l,q} u_q({x}), \end{equation}

where $u_q({x})$ ’s are independent realisations from GP priors with distinct covariances $C^{(u)}_q({x,x'})$ and $a_{l,q} \in \mathbb{R}$ ’s are the factor loadings ( $q=1,\ldots,Q$ ), considered part of the kernel hyperparameter space $\Theta$ . Thus, the semiparametric name of the model comes from the combination of a nonparametric component (several GPs) and a parametric linear mixing of the functions $u_q({x})$ .

The role of $Q \le L$ is to achieve dimension reduction for the correlation structure across the $f_l$ ’s. Let $\textbf{a}_q = (a_{1,q}, \ldots, a_{L,q})^T$ be the vector representing the collection of coefficients associated with the qth latent function across the L outputs. Then, the covariance of the vector-valued function $\textbf{f(x)} = \sum_{q=1}^Q\textbf{a}_qu_q(\textbf{x})$ is as follows:

(9) \begin{align} \notag {\textrm{Cov}}(\textbf{f}(x),\textbf{f}(x')) & = \sum_{q=1}^Q (\textbf{a}_q\textbf{a}_q^T) \otimes C^{(u)}_q({x,x'}) \\[5pt] & = \sum_{q=1}^Q (A_qA_q^T) \otimes C^{(u)}_q({x,x'}) = \sum_{q=1}^Q B_q \otimes C^{(u)}_q({x,x'}) \end{align}

where $\otimes$ symbolises the Kronecker product, $A_q = \textbf{a}_q = (a_{1,q},\ldots, a_{L,q})^T$ and each $B_q$ has rank one.

3.3 Intrinsic coregionalisation model

Similar to SLFM, ICM assumes each output function $f_l({x})$ is generated from a common pool of Q latent functions, cf. (8). However, the latent $u_q({x})$ all share the same GP prior with the covariance kernel $C^{(u)}({x,x'})$ . Then, the covariance for $\textbf{f}(x)$ is

(10) \begin{align} {\textrm{Cov}}(\textbf{f}(x),\textbf{f}(x')) & = \bigg(\sum_{q=1}^Q \textbf{a}_q\textbf{a}_q^T\bigg) \otimes {\textrm{Cov}}(u_q({x}),u_q({x'})) = B \otimes C^{(u)}({x,x'}), \end{align}

where $B = AA^T \in \mathbb{R}^{L \times L}$ has rank Q. In other words, the cross-output correlation is of rank $Q \le L$ , while the spatial covariance of each output is the same. The lth element in the diagonal of the cross-covariance matrix B (or $\sum_{q=1}^Q B_q$ in SLFM) represents the process variance of $f_l(\!\cdot\!)$ . Since $B = \sum_{q=1}^Q \textbf{a}_q\textbf{a}_q^T$ , the individual entries are $B_{l,k} = \sum_{q=1}^Q a_{l,q}a_{k,q}$ and the diagonals are $B_{l,l} = \eta^2_{l} = \sum_{q=1}^Q a^2_{l,q}$ , $1 \leq l \leq L$ . We can similarly infer the correlation matrix $R=(r_{l_1,l_2})$ between population $l_1$ and $l_2$ ( $1 \leq l_1,l_2 \leq L$ ); for ICM, it is

(11) \begin{align}r_{l_1,l_2}\;:\!=\;\frac{{\textrm{Cov}}(f_{l_1}({x}),f_{l_2}({x}))}{\sqrt{{\textrm{Var}}(f_{l_1}({x}))\times{\textrm{Var}}(f_{l_2}({x}))}} = \frac{B_{l_1,l_2}}{\sqrt{ B_{l_1,l_1} B_{l_2, l_2}}} = \frac{\sum_{q=1}^Q a_{l_1,q}a_{l_2,q}}{\sqrt{\big(\sum_{q=1}^Qa^2_{l_1,q}\big)\big(\sum_{q=1}^Qa^2_{l_2,q}\big)}}.\end{align}

In both SLFM and ICM, the fused covariance kernel belongs to the class of separable kernels (Alvarez et al. Reference Alvarez, Rosasco and Lawrence2012) and decouples using the Kronecker product into: (1) the coregionalisation matrices $B_q$ that measure the interaction between different outputs and (2) the spatial covariance over age–year dimensions $C^{(u)}_q({x,x'})$ , cf. Equations (9) & (10). In ICM, all L populations share the same spatial covariance kernel $C^{(u)}({x,x'})$ . Such assumption of a common spatial covariance over age–year inputs links to the concept of commonality in the mortality structure (but not levels) across populations. Compared to ICM, SLFM offers more flexibility at the cost of adding more hyperparameters: the total number of kernel hyperparameters in the fused covariance matrix $\textbf{C}$ of SLFM is $QL+2Q$ vis-a-vis $QL+2$ hyperparameters in the ICM.

Selecting rank Q . As Q is not one of the hyperparameters to be optimised, ad hoc ways are needed to pick it. We use the Bayesian information criterion (BIC) to select rank Q that produces the most parsimonious model; see Williams et al. (Reference Williams, Klanke, Vijayakumar, Chai, Koller, Schuurmans, Bengio and Bottou2009) and Huynh & Ludkovski (Reference Huynh and Ludkovski2021). As discussed in Bonilla et al. (Reference Bonilla, Chai and Williams2008), taking $Q<L$ in ICM corresponds to finding a rank-Q approximation (based on an incomplete Cholesky decomposition) to the full-rank $C^{(f)}$ . A similar interpretation holds for SLFM and the respective $B_q$ ’s. While attractive for dimension reduction and computational speed up, low Q may not be adequate to describe the overall dependence structure and hence clashes with the original goal of capturing the variability present in the fused mortality dataset. In particular, we observe that BIC tends to select $Q \in \{2,3\}$ which may be too small for $L \ge 5$ . Based on our case studies, we recommend $Q \in \{3, 4, 5\}$ for maximising predictive performance.

3.4 Multi-level ICM for scalable GPs

In the situation when we have multi-dimensional factor inputs (e.g. cause and gender together), one approach is to combine all factor inputs into a single covariate with L distinct outputs prior to applying ICM or SLFM. When L grows large, ICM becomes less feasible due to its time complexity $\mathcal{O}(N^3 L Q^2)$ (Bonilla et al., Reference Bonilla, Chai and Williams2008). In this section, we develop the structured Kronecker product kernel (multi-level ICM in Liu et al., Reference Liu, Ong, Shen and Cai2020; Zhe et al., Reference Zhe, Xing and Kirby2019) to mitigate this scalability issue in GP. The structured covariance kernel exploits the fact that mortality tables are gridded along each factor dimension.

We express the total number of outputs L as the product across P types of categorical inputs, $L =\prod_{p=1}^P L_p$ , where $L_p$ is the number of levels within the pth categorical input. We then decompose the cross-population covariance $\tilde{B}$ as the Kronecker product:

(12) \begin{align} \tilde{B} = \bigotimes_{p=1}^P \tilde{B}_p \end{align}

where $\tilde{B}_p,\;1\leq p \leq P$ refers to the cross-covariance matrix between sub-populations within the pth categorical input, taken to have rank $Q_p \le L_p$ . Directly marginalising the cross-covariance matrices yields a convenient interpretation of the correlation between sub-populations within a factor input and moreover allows for separate estimation of each cross-covariance sub-matrix $\tilde{B}_p$ , $1\leq p \leq P$ .

The multi-level ICM set-up implies that each output $f_l({x})$ is the weighted combination of $Q_1 \times \ldots \times Q_P$ independent latent functions, all with the spatial covariance kernel $C^{(u)}({x,x'})$ :

(13) \begin{equation} f_l({x}) = \sum_{j=1}^{Q_1\times \ldots \times Q_P} a_{l,j} u_j({x}).\end{equation}

The improvement in scalability of the multi-level ICM can be analysed via the ranks $Q_p$ of the cross-covariance sub-matrices. Thanks to the Kronecker product’s property, $\text{rank}(\tilde{B}) = \prod_{p=1}^P \text{rank}(\tilde{B}_p) = \prod_{p=1}^P Q_p$ and using Cholesky decomposition, Equation (12) can be rewritten as:

(14) \begin{align} \tilde{B} = \tilde{A}\tilde{A}^T = \bigotimes_{p=1}^P \tilde{B}_p = \bigotimes_{p=1}^P \big(\tilde{A}_p \tilde{A}_p^T \big) = \bigg(\bigotimes_{p=1}^P \tilde{A}_p \bigg)\bigg(\bigotimes_{p=1}^P \tilde{A}_p \bigg)^T, \end{align}

where $\tilde{A}_p=(\textbf{a}_1^p,\ldots,\textbf{a}_{Q_p}^p)$ , $1 \leq p \leq P$ , and each vector $\textbf{a}_k^p=(a_{1,k}^p,\ldots,a_{L_p,k}^p)^T$ , ( $1 \leq k \leq Q_p$ ) represents the collection of scalar coefficients associated with the kth latent function across $L_p$ sub-populations in the pth categorical input. Thus, the number of hyperparameters required to estimate the cross-covariance $\tilde{B}$ is $\sum_{p=1}^P Q_p L_p$ , which can be much lower than for single-level ICM when L is large. Note that when $\text{rank}(B) < \text{rank}(\tilde{B})$ , the multi-level ICM utilises more latent functions to generate the model outputs, compensating for the imposed structure in $\tilde{B}$ in (14). In terms of overall complexity, multi-level ICM requires $\mathcal{O}(N^3 (\sum_p L_p Q_p^2) )$ time compared to $\mathcal{O}( N^3 L Q^2)$ for single-level ICM.

3.5 MOGP hyperparameters

To implement a GP model requires specifying its hyperparameters. Note that actual inference reduces to linear algebraic formulas in (6)–(7), and the modeling task is to learn the spatial covariance, namely the mean and kernel functions.

Mean function: We make the prior $m_l(x)$ to be population-specific in order to maximise model flexibility in describing the mortality trend of each population. Thus, we have $3 L+1$ coefficients $\boldsymbol{\beta}=(\beta_{0},\beta^{ag}_{1,l}, \beta^{ag}_{2,l}, \beta^{yr}_{l,1} \;:\; l=1,\ldots, L)$ , cf. (4).

Observation Likelihood: We assume a constant observation noise within each population $\sigma_l = \textrm{StDev}(\epsilon^i_l)$ . This accounts for heterogeneous characteristics when observations from multiple populations are combined; in particular, $\sigma_l$ is smaller for larger populations and for more prevalent causes (Huynh et al., Reference Huynh, Ludkovski and Zail2020). The $\sigma_l$ ’s are estimated via maximum likelihood along with all other hyperparameters. More advanced GP models that either employ Poisson likelihood or infer input-dependent non-parametric $\sigma_l(x^i)$ are possible but require additional coding and are beyond the scope of this work.

Estimating Hyperparameters: For (multi-level) ICM, the set of hyperparameters is $\Theta = \left(\theta_{ag}, \theta_{yr}, (a_{l,q}), (\sigma^2_l), \boldsymbol{\beta} \right)$ ; for SLFM, it is $\Theta = \left((\theta_{ag,q}), (\theta_{yr,q}), (a_{l,q}), (\sigma^2_l), \boldsymbol{\beta} \right)$ . We use the R package kergp (Deville et al., Reference Deville, Ginsbourger, Roustant and Durrande2019) to carry out the respective maximum likelihood estimation via Kronecker decompositions. Alternatively, to account for model risk and offer more robust results, one could employ a fully Bayesian hyperparameter inference. This could be done with the Stan software (Carpenter et al., 2017) but is beyond the scope of this work.

3.6 Performance metrics

Given a test set of observed $y_*(x_*,l)$ ’s, we evaluate the effectiveness of different models using two metrics. First, we employ the mean absolute percentage error (APE) to examine the discrepancy between the observed and predicted outputs:

(15) \begin{equation} \text{APE}(y_*, m_*) \;:\!=\; \bigg|\dfrac{y_*-m_*(x_*)}{y_*}\bigg|\end{equation}

where $y_*(x_*,l)$ is the observed value at test input $x_*$ in the lth population and $m_*(x_*,l)$ is the predicted log-mortality rate. Note that APE is scale-independent, enabling us to compare model performance across populations with different exposures.

We also use the continuous ranked probability score (CRPS) metric to assess the quality of the probabilistic forecasts produced by a MOGP. Indeed, one of the major benefits of GP-based mortality models is a full distribution for future observations $y_*(x_*,l)$ which allows a more detailed uncertainty quantification beyond just looking at the predictive mean $m_*(x_*,l)$ . CRPS assesses the closeness of the realised outcome $y_*(x_*,l)$ relative to its predictive distribution $F_*(\cdot;\; x_*)$ which in the GP contest is Gaussian and leads to

(16) \begin{align} \notag \text{CRPS}( F_*, y_*) & \;:\!=\; \int_{\mathbb{R}}\left[F_*(z)-\mathbb{1}_{\{z \geq y_* \}} \right]^2 \, dz \\[5pt] & = \sqrt{s_*^2(x_*)+\sigma^2_l}\big[ \tilde{y}_* (2\Phi(\tilde{y}_*)-1) + 2 \phi(\tilde{y}_*) - \frac{1}{\sqrt{\pi}}\big], \quad \tilde{y}_* \;:\!=\; \frac{y_*-m_*(x_*)}{s_*^2(x_*)+\sigma^2_l} \end{align}

where $\phi(\!\cdot\!), \Phi(\!\cdot\!)$ are the standard Gaussian density and cdf. Observe how CRPS penalises both bias ( $2\Phi(\tilde{y}^*)-1$ ) and excessive predictive variance.

We average both APE and CRPS across a test set of age–year–population inputs. The resulting mean APE is interpreted as a normalised relative predictive error, and mean CRPS as the squared difference between the forecasted and the empirical cumulative distribution functions. Models with lower mean APE/CRPS are judged to have a better fit.

Mortality Improvement Factors. A common way to interpret a mortality surface is via the (annual) mortality improvement factors which measure longevity changes year over year. The raw annual percentage mortality improvement is $MI^{obs}_l(x)\;:\!=\; 1-\frac{\exp\big(y_l(x_{ag};x_{yr})\big)}{\exp\big(y_l(x_{ag};x_{yr}-1)\big)}$ . The smoothed improvement factor is obtained by substituting in the GP posterior means $m_{*}$ ’s:

(17) \begin{equation} MI^{GP}_l(x)\;:\!=\;\Bigg[1-\frac{\exp(m_{*}(x_{ag};\;x_{yr},l))}{\exp(m_{*}(x_{ag};\;x_{yr}-1,l))}\Bigg].\end{equation}

4.1 Commonalities in cause-specific mortality surfaces

The main difference between ICM and SLFM is the underlying assumption about the latent factors $u_q(\!\cdot\!)$ . ICM assumes that all factors have the same lengthscales $\theta_{ag},\theta_{yr}$ and is therefore appropriate for modelling homogenous mortality surfaces. SLFM is more general and fits distinct $\theta_{ag,q}, \theta_{yr,q}$ ; it is expected to perform better when the different surfaces exhibit heterogeneity (e.g. different degree of correlation across age). We examine this commonality assumption in Figure 2 by displaying side by side the mortality improvement factors of the various cancers derived from multi-cause ICM and SLFM. The BIC selection criterion yields $Q=2$ for both ICM and SLFM. In this case study, the results from SLFM are almost identical to ICM predictions. Therefore, the assumption of sharing the spatial kernel over age–year inputs across the considered cancer types is plausible. This conclusion is reinforced by Table A.1 in Appendix A which shows that the inferred lengthscales $\theta_{ag,q},\theta_{yr,q}$ for SLFM are very similar for $q=1,2$ across all three countries. In other words, both of the latent factors learned by SLFM have similar age–year spatial dependence, and so there is little loss of fidelity in a priori forcing them to be equal, as is done in ICM. Indeed, the lengthscales in SLFM are close to the ICM ones.

We can also inspect Figure 2 for insights about the relative mortality improvement trends of different cancers. Stomach cancer has the largest improvement rates for most age groups in all three countries. Decline in stomach cancer incidence tends to be associated with economic improvements resulting in healthier diet, better food preservation, clean water supply, etc. We also observe the increasing improvement trend of lung cancer among age groups below 60 years in Czech Rep. and Poland, reflecting lower smoking rates in birth cohorts after WWII. Czech males experienced large increase in the improvement rates in most cancers. In Germany, the improvement rates have been rising among age groups below 60 years but slowing down among older age groups. In Polish males, except stomach cancer, the improvement trends increased until early 2010s and then significantly declined, displaying the impact of an ageing population and an increase in lifestyle exposure to risk factors for cancers (Wojtys et al., Reference Wojtys, Antczak and Godlewski2014). The incidence of stomach cancer has been flat over time, consistent with the finding in Arash et al. (Reference Arash, Ramin, Saeid and Sepanlou2020).

Figure 3. visualises the inferred cross-cause correlation matrices $(r_{l_1,l_2} : 1 \le l_1,l_2 \le 5)$ . Both ICM and SLFM document a global positive association among mortality rates from these cancers in each country. It is consistent with strong resemblance in the mortality improvement trends between cancers in Figure 2. Since cancers within a given population share common risk factors, innovations in early detection and advancements in treatment for one cancer are likely to have positive effects on other cancer variations. Negative correlation $r_{l_1,l_2} < 0$ reflects opposite trends, for example, stomach and pancreas cancers in Poland; this may be due to a competing risks context.

Figure 3 Correlation matrices R derived from multi-cause GP models; rank $Q=2$ is chosen for both ICM and SLFM. In each country, the MOGP model is fitted on males aged 50–84 years (seven groups) during years 1998–2016, over five cancer variations: stomach, colorectal, pancreatic, lung and remaining types.

For a different take on cause-of-death commonality, we applied the multi-cause GP ICM model to jointly model $L=8$ top causes in the HCD US dataset, separately for each gender. Using BIC as the criterion, models with rank $Q=6$ (for eight populations) yield the largest BICs for both males and females. This indicates a higher degree of heterogeneity in these larger cause groupings, compared to $Q=2$ across five cancer causes in the previous study. Thus, the models employ more latent functions $u_q(\!\cdot\!)$ to adequately capture the total variability in the joint datasets. The inferred cross-cause correlation matrices are displayed in Figure B.2 in Appendix B. Overall, our results are in agreement with the SOA study (Boumezoued et al., Reference Boumezoued, Coulomb, Klein, Louvet and Titon2019), for example, confirming that there are moderate positive associations between most causes. The strongest correlations are found between heart disease and stroke, and between heart and drug overdose. As might be expected, there is little correlation between remaining causes (RMN) and most other categories. Some of the correlations vary between genders, possibly due to strong observation noise in less common causes.

A complementary way to examine commonalities across causes is to inspect the inferred factor loadings $a_{l,q}$ . Populations that have similar loadings will be highly correlated. Figure A.3 in Appendix A displays the factor loadings in a country–cause SLFM with both $Q=2$ and $Q=3$ latent factors $u_q$ . We observe that primary clustering is by causes rather than by countries. Some outliers, such as STM in Czech R., can also be seen and suggest idiosyncratic behaviour of the respective mortality surface. Less separation (and limited interpretation of $a_{l,i}$ ) is observed when using only $Q=2$ latent factors, which indirectly suggests that $Q=3$ is preferable.

4.2 Aggregating by-cause models

An important motivation for our work was to use by-cause analysis to make more precise conclusions about all-cause mortality. For example, the mortality trends of individual cancers give insights into the respective all-cancer mortality trends. Figure 4 shows the results from a multi-level country–cause GP ICM. It visualises the aggregated predictive distribution of all-cancer log-mortality observations, $y_l(x_*)$ , for male populations by country and age group, using shading to denote predictive quantiles. Note that since the model did not use all-cancer mortality during training, the fact that the predictive in-sample bands closely match the historical movement of all-cancer mortality data is a validation of a successful by-cause analysis.

Figure 4 Predictive distribution of all-cancer log-mortality rates for different age groups in three countries via multi-level country–cause ICM. The joint model is fitted on males, age groups between 50 and 84 years, years 1998–2016 across three countries ( $Q_{ctry} = 3$ ) and five cancers ( $Q_{caus}=5$ ). Shading indicates the 60%, 80%, 95% and 99% predictive quantile bands.

The effort in fighting cancer has transpired in all three countries, but the improvement is not uniform. Despite Czech Republic and Poland being socio-economically similar, Czech Rep. has a faster improvement pace than Poland. Although Germany continues to have the lowest log-mortality rates across all age groups, the Czech Rep. has drastically closed this gap in the last decade; see especially the left panel of Figure 4 (age group 55–59 years). The main driver is the rapid improvements in all common cancers in Czech Rep., for example, the mortality improvement factor for LUN being recently more than double compared to Germany, cf. Figure 2.

4.3 Trends in US top-level causes

Figure 5 presents the in-sample posterior distribution of log-mortality rates $f_l(x_*)$ during 1999–2018 along with the projected trends up to 2025 for two age groups, in both male and female US populations. The left panels display the predictive trends for the top six causes in each age group. We observe improvements in most causes and age groups; the largest improvement being in CANL. In contrast, drug abuse deaths are rising rapidly among young age groups (e.g. age 40–44 years); see the two upper-left panels of Figure 5. In the older age groups, mortality from heart disease experienced large declines in early 2000s but essentially flattened out after 2015. Boumezoued et al. (Reference Boumezoued, Coulomb, Klein, Louvet and Titon2019) emphasised the need for break-point detection to improve forecasts of such causes whose trends change over time. In the MOGP framework, the forecasts are driven by the most recent data and get automatically gradually adjusted if trends shift. So, for example, we do not need to do anything further to achieve the slow pace of future HEA improvement as shown in Figure 5.

Figure 5 Posterior distribution of true log-mortality rates for US males and females for the 40–44 and 65–69 years of age groups. In each row, left panel shows the posterior quantiles for top six individual causes (via multi-cause GP), middle panel shows aggregate log-mortality trends via by-cause model (multi-cause GP) and right panel shows aggregate log-mortality trends via all-cause model (single-output GP). All models are fitted on ages 40–69 years and years 1999–2018. The vertical lines indicate the boundary between in-sample (1999–2018) and out-of-sample forecast (2019–2025). Shading indicates the 60%, 80%, 95% and 99% predictive quantile bands. We further overlay the 50% quantile (predictive median) of by-cause and all-cause models for the out-of-sample period for convenient comparison.

The middle and right columns in Figure 5 compare the aggregate all-cause projections for the US population. We witness the pessimism of the aggregate projected trends based on the by-cause models compared to an all-cause SOGP (right column), especially in the younger ages. This discrepancy, driven by the growing importance of causes with increasing trends, such as drug overdose, highlights the additional insights from by-cause modelling. The pessimism of by-cause analysis was first mentioned in Wilmoth (Reference Wilmoth1995) and re-iterated in Boumezoued et al. (Reference Boumezoued, Coulomb, Klein, Louvet and Titon2019). For older age groups, the underlying dynamics among common causes are more stable and all-cause and by-cause forecasts are broadly similar.

Note that compared to Figure 4, the GP posterior uncertainty bands widen dramatically as we go from in-sample (up to year 2018) to extrapolating for years 2019–2025. This reflects the data-driven nature of GP forecasts which intrinsically leads to low uncertainty for in-sample smoothing and widening uncertainty as predictions are made further into the future. This phenomenon gets amplified as we add up the by-cause forecasts to obtain all-cause predictions and witness the wider band in the middle panels of Figure 5 relative to the right panel based on a SOGP.

As discussed in Huynh & Ludkovski (Reference Huynh and Ludkovski2021), MOGPs are well suited to generate expert-based projections. This is a useful feature to have given that by default, projections are driven by the historical trends that might not continue in the future. For instance, the increasing trend of drug abuse is largely fuelled by the opioid epidemic. Assuming this crisis is addressed in the future, and the projected DRU mortality should be adjusted downwards. In MOGP, this can be achieved by modifying the year trend in $m(\!\cdot\!)$ . Figure B.1 in the Appendix displays an illustration where the trend of drug abuse is reduced by one-third (through lowering the year effect $\beta^{yr}$ by one-third) of the original pace for both the male and female populations. The resulting adjusted forecast for aggregate mortality gets closer to that from the all-cause model, reducing the level of pessimism we have observed earlier among the US young population. Another potential adjustment could be for smoking-induced cancer, where the MOGP models extrapolate the historical trend of rapid longevity gains. However, the SOA report (Boumezoued et al., Reference Boumezoued, Coulomb, Klein, Louvet and Titon2019) suggests that this pace might not take place in the intermediate term, lowering aggregate mortality gains.

5.1 Multi-level versus single-level correlation structure

Figure 6 displays the inferred correlation structure $r_{l_1,l_2} \;:\; 1 \le l_1,l_2 \le 30$ across the above 30 populations. On the left, we fit a multi-level ICM ( $Q_{ctry}=3,\;Q_{caus}=5,\;Q_{gndr}=2$ ) and on the right a single-level ICM ( $Q=2$ ). Both models are fitted on age groups 50–84 years, years 1998–2016 for three countries, five cancers and two genders. Observe that the multi-level model has $\sum_p Q_p L_p + 2 = 3 \times 3 + 5 \times 5 + 2 \times 2 + 2 = 38$ hyperparameters compared to $Q L + 2 = 62$ in the single-level model. The right panel does not display any recognisable structure in the inferred correlations because the model is not aware of the different factor dimensions; after dimension reduction, the marginal associations between sub-populations within the original factor inputs are no longer accessible. In contrast, the multi-level model enforces a block structure in the correlation matrix R; see Figure 6(a). Recall that the correlation sub-matrices for each factor input are estimated separately and the Kronecker product structure implies that we can read off the correlation among any combination of factors. Figure 7 displays the derived sub cross-correlation matrices from the above three-level country–cause–gender GP model. One can then multiply these factor-based correlations to get the total correlation $r_{(c_1,s_1,g_1),(c_2,s_2,g_2)} = \prod_{p \in \{c,s,g\}} r_{p_1,p_2}$ in Figure 6(a). For example, the correlation of cancer mortality between Czech-lung-male and Polish-pancreas-male is $r_{CZE,\;POL} \times r_{LUN,\;PAN} \times r_{MAL,\;MAL} = 0.88 \times 0.91 \times 1 \approx 0.80$ . This is in fact also the correlation between Czech-lung-female and Polish-pancreas-female ( $r_{CZE,\;POL} \times r_{LUN,\;PAN} \times r_{FEM,\;FEM}$ ) mortality. The limited number of deaths in some causes explains differences in the correlation matrices between single- and multi-level ICM. For example, the correlations between Polish-stomach-female and other cancer types are mostly negatively correlated in ICM but (mildly) positively correlated in the multi-level ICM.

Figure 6 Cross-correlation matrices of MOGP models that incorporate country, cause and gender as categorical inputs. Thirty total populations.

Figure 7 Cross-correlation matrices derived from a three-level country–cause–gender multi-level ICM model ( $Q_{ctry}=3, Q_{caus}=5,\;Q_{gndr}=2$ ), fitted on age groups 50–84 years and years 1998–2016.

5.2 Model selection

Next, we compare the predictive performance of multi-level versus single-level MOGP-ICM. Following the same set-up as Table 1, Table 2 shows the 3-year median improvement in APE and CRPS for different joint models relative to SOGP models in a multinational context. Joint models produce more accurate mean forecasts with higher credibility, but the predictive gains are not uniform across countries. The largest improvements from multi-level ICM are in Czech Republic; Czech raw data have lower credibility due to its smaller exposures; thus, there is more opportunity for data fusion. The results further validate our approach of modelling cause-specific mortality across populations: models that incorporate information from foreign countries (e.g. country–cause GP) have larger predictive improvements compared to Table 1. In Poland due to the recent structural break encountered in mid-2010s (see Figure 4), the performance of aggregated country–cause MOGP is consistently worse than that of all-cause SOGP; this issue is rectified for country–cause–gender MOGPs. We expect this issue to self-correct in a couple of years if the new trend persists.

Table 2. Comparison between MOGP models with different ranks Q in terms of APE and CRPS metrics. The reported values are median relative improvements of MOGP versus SOGP of 1-year-out aggregated all-cancer forecasts for age groups 50–84 years based on three training periods: 1998–2013 (predict 2014), 1999–2014 (predict 2015) and 2000–2015 (predict 2016). Top half: MOGP models fitted on five cancer types and three countries. Bottom half: MOGP models fitted on five cancer types, three countries and two genders.

In contrast to section 4 where the performance of ICM and SLFM was almost identical, with multiple input factors ICM usually outperforms SLFM. Given the strong commonality in mortality trends across cancer types, a shared age–year covariance kernel appears preferable for information fusion. The comparison between single- and multi-level ICM depends on the number of populations to model. In country–cause–gender setting with $L=30$ populations, the multi-level ICM ( $Q_{ctry}=3,\;Q_{caus}=5,\;Q_{gndr}=2$ ) yields better mean APE and CRPS than single-level ICM and moreover uses fewer hyperparameters. For country–cause setting, the performance is comparable. Note that to make an apples-to-apples comparison between single- and multi-level ICM, one ought to equalise their number of hyperparameters, rather than the Q’s.

Table 2 also shows the impact of the latent ranks Q and $Q_p$ ’s. For country–cause, $Q=4$ tends to yield the best results in single-level models, but for country–cause–gender, ICM with $Q=3$ performs consistently worse than $Q=2$ , presumably due to unstable estimates of the more than 90 underlying hyperparameters. For multi-level ICM, we generally find that full-rank $Q_p = L_p$ works best, although low-rank set-ups $Q_p < L_p$ also yield good predictive performance, indicating the opportunity to shrink even further the number of kernel hyperparameters.

Figure 8 shows the predicted log-mortality rates for individual cancers and age group 55–59 years via full-rank multi-level ICM and single-level ICM with $Q=2$ . Both models are fitted on age groups 50–84 years and years 1998–2016 before we perform out-of-sample forecasts for the next 3 years (2014–2016). The single-level ICM produces over-smoothed forecasts $m_*(\!\cdot\!)$ for several cancers, especially cancers with large observation noise like stomach and pancreas; this problem is mitigated by the shorter lengthscale in year in multi-level ICM ( $\theta_{yr} \approx 8.8$ versus $\theta_{yr} \approx 14.8$ in single-level ICM).

Figure 8 Predictive log-mortality distributions from single-level ICM (listed as ICM, $Q=2$ ) and multi-level ICM (listed as Hier ICM, $Q_{ctry}=3,\;Q_{caus}=5,\;Q_{gndr}=2$ ) models with three factor inputs: country, cause and gender (30 populations). All models are fitted on age groups 50–84 years and years 1998–2013 and applied for 3-year-out forecasts up to 2016. The shadings indicate 95% predictive bands; the vertical lines mark the edge of training data. Note different y-axes in each panel.

Coherence in cause-specific trends: Figure 8 demonstrates that males and females do not always share similar progress in mortality reduction. While the trends in stomach, colorectal and pancreatic cancers are consistent for both genders, for lung cancer the male–female gap is diminishing rapidly, especially in Germany and Poland. This is driven by a decrease in cigarette consumption among men, while women are more likely to develop lung cancers that are not associated with smoking. Thus, the concept of forecast coherence (namely extrapolating a stable male–female spread, as observed historically) is not always well suited for cause-of-death analysis.

5.3 Borrowing the latest datasets from other populations

The period coverage in HCD varies by country as datasets for countries are uploaded asynchronously. This implies that some countries have more up-to-date datasets than others; see Figure A.1 in Appendix A.1 and offers opportunities to fuse data from other countries to update domestic forecasts. For Age-Period-Cohort models, such as Li & Lee (Reference Li and Lee2005), fusion is generally challenging as model fitting relies on having a rectangular dataset. Our MOGP framework can straightforwardly handle such “notched” datasets to take full advantage of additional observations in different countries.

As an illustration, we consider joint modelling of male mortality in France and Germany. These are two developed Western European countries with similar demographics. Relative to Germany which has data for 1998–2016, France at the time of writing covers only 2000–2015. We choose French males in 2015 as the target population and examine its 1-year-out prediction quality for several models. We take French observations for 2000–2014 and borrow the more up-to-date dataset from Germany (1998–2015) to implement the notched set-up.

We proceed to compare six different models; see Table A.2 in the Appendix. Our comparison covers single-output models for France, multi-cause models for France and multi-cause models joint for France and Germany. Figure 9 visualises relative performance by comparing two components, both for French males in 2015 and across three different age bands: (a) predictive standard deviation $s_*(x_*)$ (top panels) and (b) prediction errors $m_*(x_*)-y_*(x_*)$ relative to the realised 2015 all-cancer log-mortality rates (bottom). Our benchmark is French males SOGP fitted on all-cancer log-mortality rates from 2000 to 2015, that is, the ideal case where latest domestic data are already available. As expected, with access to only 2000–2014 French data, predictions have less credibility (higher $s_*(x_*)$ ) than the benchmark model that performs in-sample smoothing. However, a country–cause MOGP that ingests both French and German data yields lower $s_*(x_*)$ , even without seeing German 2015 mortality. Similarly, multi-cause models boost credibility compared to all-cancer analysis.

Figure 9 Comparison of the prediction accuracy for 2015 all-cancer log-mortality of French males for indicated age groups between different models with “notched” set-up. Top row: predictive standard deviation $s_*(x_*)$ ; bottom row: discrepancy between the predictive mean $m_*(x_*)$ and the observed value $y_*(x_*)$ .

For the forecasts, we observe that fusing German data yields prediction errors $|m_*(x_*)-y_*(x_*)|$ that are competitive to that from the benchmark. In fact, the prediction quality of the notched country–cause MOGP (Germany ‘98-‘15, France ‘00-‘14) is as good as the multi-cause MOGP that simultaneously models the log-mortality rates of all five cancer types in French males with 2015 observations available. Throughout, there are the intuitive gains from having 2015 rather than only 2014 experience. In particular, we observe material improvement from access to 2015 German data (right-most two points in Figure 9) which lowers prediction errors.