C[1]¿p#1 R[1]¿p#1 L[1]¿p#1

# Joint Models for Cause-of-Death Mortality in Multiple Populations

Nhan Huynh and Mike Ludkovski Department of Statistics and Applied
Probability, University of California at Santa BarbaraDepartment of Statistics
and Applied Probability, University of California at Santa Barbara,
<EMAIL_ADDRESS>
###### Abstract

We investigate jointly modeling Age-specific rates of various causes of death
in a multinational setting. We apply Multi-Output Gaussian Processes (MOGP), a
spatial machine learning method, to smooth and extrapolate multiple cause-of-
death mortality rates across several countries and both genders. To maintain
flexibility and scalability, we investigate MOGPs with Kronecker-structured
kernels and latent factors. In particular, we develop a custom multi-level
MOGP that leverages the gridded structure of mortality tables to efficiently
capture heterogeneity and dependence across different factor inputs. Results
are illustrated with datasets from the Human Cause-of-Death Database (HCD). We
discuss a case study involving cancer variations in three European nations,
and a US-based study that considers eight top-level causes and includes
comparison to all-cause analysis. Our models provide insights into the
commonality of cause-specific mortality trends and demonstrate the
opportunities for respective data fusion.

##  1 Background and motivation

In-depth modeling of the evolution of human mortality necessitates analysis of
the prevalent causes of death. This is doubly so for making mortality
forecasts into the future across different age groups, populations and
genders. In this article we develop a methodology for probabilistic
forecasting of cause-specific mortality in a multi-population (primarily
interpreted as a multi-national) context. Thus, we simultaneously fit multiple
cause-specific longevity surfaces via a spatio-temporal model that accounts
for the complex dependencies across causes and countries and across the Age-
Year dimensions.

While there have been many works on modeling mortality across several
populations (Dong et al. 2020, Enchev et al. 2017, Guibert et al. 2019,
Hyndman et al. 2013, Kleinow 2015, Li and Lu 2017, Tsai and Zhang 2019), as
well as an active literature on cause-of-death mortality, there are very few
that do both simultaneously. As we detail below, there are many natural
reasons for building such a joint model, and this gap is arguably driven by
the underlying “Big Data” methodological challenge. Indeed, with dozens of
mortality datasets that are indexed by countries, causes-of-death, genders,
etc., developing a scalable approach is daunting. We demonstrate that this
issue may be overcome by adapting machine learning approaches, specifically
techniques from multi-task learning (Bonilla et al. 2008, Caruana 1997, Letham
and Bakshy 2019, Williams et al. 2009). To this end, we employ multi-output
Gaussian Processes (MOGP) combined with linear coregionalization. GPs are a
kernel-based data-driven regression framework that translates mortality
modeling into smoothing and extrapolating an input-output response surface
based on noisy observations. It yields a full uncertainty quantification for
mortality rates and mortality improvement factors. Coregionalization is a
dimension reduction technique that enables efficiently handling many
correlated outputs.

This work is a continuation of our series of articles Ludkovski et al. (2018),
Huynh et al. (2020) and Huynh and Ludkovski (2021) that discussed the
application of GPs to model all-cause mortality in the single-population and
multi-population contexts, respectively. Unlike all-cause mortality in
different geographic regions, which tends to exhibit strong correlation and
long-term coherence, different causes have less commonality, and thus require
a more flexible structure for the respective cross-dependence. Moreover, joint
analysis of 10+ mortality surfaces brings computational scalability challenges
of potentially hundreds of model parameters to calibrate. Thus, in order to
carry out cause-specific mortality analysis, we make methodological
innovations along two directions. First, we compare two distinct versions of
the MOGP that implement dimension reduction by fusing outputs through linear
combinations of latent functions: the Semiparametric Latent Factor Model
(SLFM) and the Intrinsic Coregionalization Model (ICM). ICM assumes a fixed
spatial kernel in Age-Year, while SLFM does not. This distinction corresponds
to different assumptions about the structural commonality in the modeled
mortality surfaces. We conduct sensitivity analysis to assess these two
choices for the tasks of in-sample fitting and of forecasting cause-specific
mortality. Second, we implement a multi-level MOGP-ICM model that separates
latent factors across the different types (countries, causes, genders, etc.)
of categorical inputs describing the populations. The separability assumption
in the joint covariance kernel yields a product-type cross-population
correlation structure, providing insights about the relationships between the
mortality improvement trends.

Our models are driven by the scalability issue which has been a critical
obstacle to analyze _in bulk_ the large-volume mortality datasets that have
become available recently. Thus, to cope with many mortality surfaces, we
leverage the twin pillars of multi-output models (Teh et al. 2005, Letham and
Bakshy 2019, Williams et al. 2009) and the structured Kronecker covariance
that mitigates the typical cubic computational complexity of GPs (Flaxman et
al. 2015, Gilboa et al. 2015, Saatçi 2011, Zhe et al. 2019).

Decomposing all-cause mortality rates leads to reduced signal-to-noise ratio
since less common causes intrinsically have limited death counts.
Consequently, cause-specific analysis must contend with much noisier data. One
motivation for the MOGP approach is to explicitly enable data fusion across
populations, sharing information to improve model fitting. We demonstrate
significant de-noising of mortality experience that successfully captures
cause-specific trends, including for causes with low data credibility. We also
document the benefit of joint models to reduce model risk, i.e. improved
inference of model hyper-parameters. As another feature, our framework can
handle non-rectangular datasets, e.g. countries with different period
coverages. In one of our case studies we exploit this to borrow the most
recent data from other countries to update predicted domestic mortality rates.

Literature Review

One of the few works approaching cause-specific mortality modeling within a
multi-population context is the recent study in Lyu et al. (2021). The authors
introduced a nested model in the spirit of Li and Lee (2005) to jointly model
major causes from three European countries by capturing the cross-cause and
cross-country dependencies through common factors.

Since a given death is associated with a single cause-of-death, direct
dependence among causes is not observable. As such, many researchers choose to
model and forecast each cause in isolation. A variety of forecasting methods
for individual causes are employed: univariate time-series, such as ARIMA
methods in Caselli (1996), Knudsen and McNown (1993), McNown and Rogers
(1992), dynamic parametrization in Tabeau et al. (1999), least squares methods
and variations of the Lee-Carter model in Caselli et al. (2019). To capture
dependence between causes, a common approach is via copulas within the
framework of dependent competing risks. The main effort is to characterize the
joint distribution of survival times in terms of un-observed cause-specific
mortality rates, see Dimitrova et al. (2013), Lo and Wilke (2010) and Li and
Lu (2019). Alternatively, Alai et al. (2018) utilized multinomial logistic
regression. Outside the competing risks framework, Arnold and Sherris(2013)
(-Gaille) proposed a multivariate Vector Error Correction time-series model to
examine the existing cointegration relationships. Another approach is to link
multiple causes through a list of clinical factors, then applying a stochastic
model to forecast these factors (Foreman et al. 2018).

Several studies relied on compositional data analysis to achieve coherence in
the sense that cause-specific forecasts sum to the all-cause forecast. This
entails modeling the by-cause distribution of deaths, directly incorporating
dependence between causes, see Bergeron-Boucher et al. (2017), Kjaergaard et
al. (2019). We refer to Wilmoth (1995), Caselli et al. (2019), Tabeau et al.
(1999) for discussions on the usefulness and limitations in using cause-
specific models to project all-cause mortality.

The remainder of this paper is organized as follows. Section 2 introduces
cause-of-death mortality datasets from the HCD. Section 3 describes the MOGP-
ICM framework and its extensions within multinational context. Section 4
focuses on how MOGPs can maximize predictive gains over single-population
models and provide insights about the projected trends of aggregate mortality.
Section 5 compares the results from Multi- and Single-level ICM. Finally,
Section 6 concludes with main findings and directions for further analysis.

##  2 The Human Cause-of-death Database

The Human Cause-of-death Database (HCD) (HCD 2021) provides detailed cause-
specific mortality data for more than a dozen developed countries. The HCD
offers three levels for the classification of causes. The short list (with 16
broad categories) and the intermediate list (103 categories) are the same
across countries, while the full list is country-specific. The data for each
country is organized by calendar years, age groups, gender, and causes.
Datasets for different countries do not line up, both due to historical
availability and different timelines for updates, see Figure 10 in Appendix A.

As part of their extensive and well-documented post-processing the HCD
conducted a series of bridge-coding studies to reduce disruptions in mortality
trends due to changes in the International Classification of Disease (ICD).
The current 10th Revision of the ICD is far more detailed with the addition of
8,000 categories compared to the 9th Revision (Anderson et al. 2001). As
cause-of-death records switched to the 10th Revision, death counts were
shifted among some categories. To minimize such jumps, the HCD reconstructed
death assignments between the old and the new ICD Revisions; this
reconciliation is already part of the HCD datasets we used.

In our first case study we analyze sub-categories of Cancer, available from
the 103-category intermediate list. Being the leading cause of death worldwide
(WHO 2021), it is useful to understand the trends in cancer mortality for
different age groups and explore the dependence between its common variations.
Moreover, cancer types generally do not feature any substantial trend
discontinuity due to ICD revisions (Anderson et al. 2001), allowing a better
assessment of model performance in terms of in-sample smoothing and out-of-
sample forecasting.

Figure 1: Log-mortality rates of selected Cancers by Country & Age groups
among Males. Note the different $y$-axes in each panel. Source: HCD.

For our testbed we select 5 variations of cancer: lung and bronchus (LUN),
colon and rectum (COL), pancreas (PAN), stomach (STM), and all other cancers
(RMN) from 3 countries in the HCD: Czech Republic (CZE), Germany (DEU), and
Poland (POL). The chosen Cancer causes-of-death are common in both Male and
Female populations, enabling us to jointly model mortality rates across Cause,
Country, and Gender factors. For example, we exclude breast cancer from the
study as the respective male counts are very limited. Age in the HCD is
formatted as discrete five-year age groups; we treat it as a continuous
covariate encoded via the respective group average (52 for Age group 50–54, 57
for group 55–59, etc). In Figure 1, we visualize the raw log-mortality rates
across different cancer types in Male populations from Czech Rep., Germany,
and Poland. Lung and Colorectal are the leading causes of cancer deaths, while
Pancreatic and Stomach are less prevalent and exhibit more volatility. Figure
11 in the Appendix visualizes the respective patterns in Age; we observe a
strong convexity in Age, especially for LUN.

Our second case study uses a subset of the HCD data for United States, based
on the CDC’s National Center for Health Statistics, where we decompose all-
cause mortality into its major categories. This analysis is inspired by a
similar study conducted by the SOA (Boumezoued et al. 2019) that employed a
multivariate Lee-Carter model. The SOA study looked at 11 major causes, all
Age groups between 0 and 95+ in 1999–2016 with data combined from the HCD and
the Global Burden of Disease project. We restrict analysis to Ages 40–69 and
years 1999–2018, and moreover reduce to the eight most common causes for the
examined age groups: Heart (HEA), Stroke (STK), Diabetes (DIA), Cancers not
induced by smoking (CAN), (lung) Cancers induced by smoking (CANL),
Respiratory (RES), Drug abuse (DRU), and all other remaining causes (RMN). The
cause mapping for this case study closely follows the SOA guidance.

###  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:

$\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}},$ (1)

where $D^{i}_{l}$ is the number of deaths in the $l$th population for the
$i$th observation with the corresponding Age & Year inputs. The denominator
$E^{i}_{l}$ refers to the corresponding exposed-to-risk counts. Note that
different causes in the same country-gender combination with have the same
$E^{i}_{l}$, but different $D^{i}_{l}$’s. Our GP models work with log-
mortality rates so that our data is

$y^{i}_{l}=\log\mu^{i}_{l}.$

The overall dataset is then represented as $(x^{i},y^{i}_{l})$,
$i=1,\ldots,N_{l}$, $l=1,\ldots,L$.

##  3 Methodology

###  3.1 Gaussian Process Regression for Mortality Tables

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

$y^{i}_{l}=f_{l}(x^{i})+\epsilon^{i}_{l},\qquad\qquad i=1,\ldots,N,$ (2)

where $x^{i}$ represents an individual entry in the mortality table (indexed
by Age, Year), $y^{i}_{l}$ is the observed output in the $l$th 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 summarize 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}(\mathbf{x})=(f_{l}(x^{1}),\ldots,f_{l}(x^{n}))$ at $n$ inputs follows
the multivariate Gaussian distribution

$f_{l}(x^{1}),\ldots,f_{l}(x^{n})\sim\mathcal{N}\big{(}\mathbf{m_{l}(x)},\mathbf{C}_{l}(\mathbf{x,x})\big{)},$

where ${m}_{l}(x)=\mathbb{E}[f_{l}(x)]$ is the mean vector of size $n$ and
${C_{l}(x,x^{\prime})}=\mathbb{E}[(f_{l}({x})-{m_{l}(x)})(f_{l}({x^{\prime}})-{m_{l}(x^{\prime})})]$
is the $n$-by-$n$ covariance matrix. All properties of a GP are thus
completely described by its mean and covariance functions.

The functions $m_{l}(\cdot),C_{l}(\cdot,\cdot)$ characterize our prior beliefs
about the response surface $f_{l}$. The covariance kernel of the GP defines a
similarity between pairs of data points. It characterizes 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, 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^{\prime}$ as follows:

$C_{l}(x,x^{\prime})=\prod_{k\in\\{ag,yr\\}}\bigg{(}1+\dfrac{\sqrt{5}}{\theta_{k,l}}|x_{k}-x_{k}^{\prime}|+\dfrac{5}{3\theta_{k,l}^{2}}|x_{k}-x_{k}^{\prime}|^{2}\bigg{)}\exp\bigg{(}-\dfrac{\sqrt{5}}{\theta_{k,l}}|x_{k}-x_{k}^{\prime}|\bigg{)}.$
(3)

This kernel is parametrized by the Age lengthscale $\theta_{ag,l}$ and the
Year lengthscale $\theta_{yr,l}$.

The mean function $m_{l}(x)$ describes the relevant trends in log-mortality
rates. We fit a parametric prior mean to capture the long-term longevity
features: $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. For example, to reflect a quadratic trend in Age (Figure 11 in
Appendix A) and linear trend in Year dimension we use:

$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}x^{i}_{yr}.$
(4)

We are interested in the posterior distribution for $\mathbf{f_{*}}\equiv
f_{l}(\mathbf{x_{*}})$ at specified inputs $\mathbf{x_{*}}$ given the dataset
$\mathbf{y}_{l}=(y^{1}_{l},\ldots,y^{N}_{l})$, $p(\mathbf{f_{*}|y}_{l})$, in
other words the likelihood of the true response surface being $\mathbf{f_{*}}$
given what we have observed. Using
$\mathrm{Cov}(y^{i}_{l},y^{j}_{l})=\mathrm{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
$\mathbf{y}_{l}\sim\mathcal{N}(\mathbf{m_{l}(x)},\mathbf{C_{l}(x,x)}+\mathbf{\Sigma}_{l})$
where the error terms are assumed to be independent and Gaussian-distributed,
$(\epsilon^{1}_{l},\ldots,\epsilon^{N}_{l})\sim\mathcal{N}(\mathbf{0},\mathbf{\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

$p(\mathbf{f_{*}|y}_{l})\sim\mathcal{N}(\mathbf{m_{*}(x_{*},l),C_{*,l}(x_{*},x_{*})}).$

The Universal Kriging equations below provide not only the posterior mean
$m_{*}(\cdot,l)$ and posterior variance $s^{2}_{*}(\cdot,l)$ but also the
estimated coefficients $\bm{\beta}_{l}=(\beta_{1,l},\ldots,\beta_{p,l})^{T}$.
Let $\mathbf{h}(x)=\big{(}h_{1}(x),\ldots,h_{p}(x)\big{)}$,
$\mathbf{H}=\big{(}\mathbf{h}(x^{1}),\ldots,\mathbf{h}(x^{N})\big{)}$, and
$\mathbf{D=(C+\Sigma)^{-1}H}$ where $\mathbf{C}$ is the covariance matrix
$C_{l}(x^{i},x^{j})_{i,j=1}^{N}$. The posterior mean of $\bm{\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:

$\displaystyle\bm{\hat{\beta}}_{l}$ $\displaystyle\
=(\mathbf{H}^{T}\mathbf{D})^{-1}\mathbf{H}^{T}\mathbf{(C+\Sigma)}^{-1}\mathbf{y};$
(5) $\displaystyle m_{*}(x_{*},l)$ $\displaystyle\
=\mathbf{h}(x_{*})^{T}\bm{\hat{\beta}}+\mathbf{c}(x_{*})^{T}\mathbf{(C+\Sigma)}^{-1}(\mathbf{y-H}\bm{\hat{\beta}});$
(6) $\displaystyle s^{2}_{*}(x_{*},l)$ $\displaystyle\
=(\mathbf{h}(x_{*})^{T}-\mathbf{c}(x_{*})^{T}\mathbf{D})^{T}(\mathbf{H}^{T}\mathbf{D})^{-1}(\mathbf{h}(x_{*})^{T}-\mathbf{c}(x_{*})^{T}\mathbf{D})$
(7)

where $\mathbf{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 $l$th 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 $\mathbf{f(x)}=(f_{1}(\mathbf{x}),\ldots,f_{L}(\mathbf{x}))$, where
the functions $\\{f_{l}({x})\\}_{l=1}^{L}$ are the log-mortality surfaces for
the corresponding $l$th population. Similar to single-population GP above, we
place a GP prior over the latent function $\mathbf{f}$ such that:

$\mathbf{f}\sim\mathcal{GP}(\mathbf{m,C}),$

where $\mathbf{m}$ is the mean vector function whose elements
$\\{m_{l}({x})\\}_{l=1}^{L}$ are the mean functions of each output, and
$\mathbf{C}$ is the fused 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\leq l,k\leq 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 Semiparametric Latent Factor Model (SLFM) dating back to Teh et al.
(2005), each output $f_{l}({x})$ is assumed to be a linear combination of $Q$
latent functions:

$f_{l}({x})=\sum_{q=1}^{Q}a_{l,q}u_{q}({x}),$ (8)

where $u_{q}({x})$’s are independent realizations from GP priors with distinct
covariances $C^{(u)}_{q}({x,x^{\prime}})$, 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\leq L$ is to achieve dimension reduction for the correlation
structure across the $f_{l}$’s. Let
$\mathbf{a}_{q}=(a_{1,q},\ldots,a_{L,q})^{T}$ be the vector representing the
collection of coefficients associated with the $q$th latent function across
the $L$ outputs. Then the covariance of the vector-valued function
$\mathbf{f(x)}=\sum_{q=1}^{Q}\mathbf{a}_{q}u_{q}(\mathbf{x})$ is:

$\displaystyle\mathrm{Cov}(\mathbf{f}(x),\mathbf{f}(x^{\prime}))$
$\displaystyle=\sum_{q=1}^{Q}(\mathbf{a}_{q}\mathbf{a}_{q}^{T})\otimes
C^{(u)}_{q}({x,x^{\prime}})$
$\displaystyle=\sum_{q=1}^{Q}(A_{q}A_{q}^{T})\otimes
C^{(u)}_{q}({x,x^{\prime}})=\sum_{q=1}^{Q}B_{q}\otimes
C^{(u)}_{q}({x,x^{\prime}})$ (9)

where $\otimes$ symbolizes the Kronecker product,
$A_{q}=\mathbf{a}_{q}=(a_{1,q},\ldots,a_{L,q})^{T}$ and each $B_{q}$ has rank
one.

###  3.3 Intrinsic Coregionalization 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^{\prime}})$. Then, the covariance for $\mathbf{f}(x)$ is:

$\displaystyle\mathrm{Cov}(\mathbf{f}(x),\mathbf{f}(x^{\prime}))$
$\displaystyle=\bigg{(}\sum_{q=1}^{Q}\mathbf{a}_{q}\mathbf{a}_{q}^{T}\bigg{)}\otimes\mathrm{Cov}(u_{q}({x}),u_{q}({x^{\prime}}))=B\otimes
C^{(u)}({x,x^{\prime}}),$ (10)

where $B=AA^{T}\in\mathbb{R}^{L\times L}$ has rank $Q$. In other words, the
cross-output correlation is of rank $Q\leq L$, while the spatial covariance of
each output is the same. The $l$th 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}\mathbf{a}_{q}\mathbf{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

$\displaystyle
r_{l_{1},l_{2}}:=\frac{\mathrm{Cov}(f_{l_{1}}({x}),f_{l_{2}}({x}))}{\sqrt{\mathrm{Var}(f_{l_{1}}({x}))\times\mathrm{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}^{Q}a^{2}_{l_{1},q}\big{)}\big{(}\sum_{q=1}^{Q}a^{2}_{l_{2},q}\big{)}}}.$
(11)

In both SLFM and ICM, the fused covariance kernel belongs to the class of
separable kernels (Alvarez et al. 2012) and decouples using the Kronecker
product into: (1) the coregionalization matrices $B_{q}$ that measure the
interaction between different outputs and (2) the spatial covariance over Age-
Year dimensions $C^{(u)}_{q}({x,x^{\prime}})$, cf. Eqs (9) & (10). In ICM, all
$L$ populations share the same spatial covariance kernel
$C^{(u)}({x,x^{\prime}})$. 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 $\mathbf{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 optimized,
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. (2009) and Huynh and Ludkovski (2021). As discussed in Bonilla
et al. (2008), 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\geq 5$. Based on our case
studies, we recommend $Q\in\\{3,4,5\\}$ for maximizing 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}LQ^{2})$ Bonilla et al. (2008). In this section, we develop
the structured Kronecker product kernel (multi-level ICM in Liu et al. (2019),
Zhe et al. (2019)) 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 $p$th categorical input. We then decompose the cross-
population covariance $\tilde{B}$ as the Kronecker product:

$\displaystyle\tilde{B}=\bigotimes_{p=1}^{P}\tilde{B}_{p}$ (12)

where $\tilde{B}_{p},~{}1\leq p\leq P$ refers to the cross-covariance matrix
between sub-populations within the $p$th categorical input, taken to have rank
$Q_{p}\leq L_{p}$. Directly marginalizing 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 setup 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^{\prime}})$:

$f_{l}({x})=\sum_{j=1}^{Q_{1}\times\ldots\times Q_{P}}a_{l,j}u_{j}({x}).$ (13)

The improvement in scalability of the multi-level ICM can be analyzed 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, Eq. (12) can be rewritten as:

$\displaystyle\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},$
(14)

where $\tilde{A}_{p}=(\mathbf{a}_{1}^{p},\ldots,\mathbf{a}_{Q_{p}}^{p})$,
$1\leq p\leq P$, and each vector
$\mathbf{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
$k$th latent function across $L_{p}$ sub-populations in the $p$th 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}LQ^{2})$ for single-level ICM.

Remark: In the typical situation, $L_{p}$’s are small and so it is feasible to
consider full-rank multi-level ICM, i.e. $Q_{p}=L_{p}$. Otherwise, (14) allows
to exploit simultaneously the Kronecker product structure, as well as the low-
rank approximation. See Table 2 for results on the impact of $Q_{p}$ values in
multi-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 maximize model flexibility in describing the mortality trend of each
population. Thus, we have $3L+1$ coefficients
$\bm{\beta}=(\beta_{0},\beta^{ag}_{1,l},\beta^{ag}_{2,l},\beta_{yr,l}:l=1,\ldots,L)$,
cf. (4).

Observation Likelihood: We assume a constant observation noise within each
population $\sigma_{l}=\mathrm{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. 2020). 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 hyper-parameters
is
$\Theta=\left(\theta_{ag},\theta_{yr},(a_{l,q}),(\sigma^{2}_{l}),\bm{\beta}\right)$;
for SLFM it is
$\Theta=\left((\theta_{ag,q}),(\theta_{yr,q}),(a_{l,q}),(\sigma^{2}_{l}),\bm{\beta}\right)$.
We use the R package kergp (Deville et al. 2019) 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:

$\text{APE}(y_{*},m_{*}):=\bigg{|}\dfrac{y_{*}-m_{*}(x_{*})}{y_{*}}\bigg{|}$
(15)

where $y_{*}(x_{*},l)$ is the observed value at test input $x_{*}$ in the
$l$th 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 realized outcome $y_{*}(x_{*},l)$ relative
to its predictive distribution $F_{*}(\cdot;x_{*})$ which in the GP contest is
Gaussian and leads to

$\displaystyle\text{CRPS}(F_{*},y_{*})$
$\displaystyle:=\int_{\mathbb{R}}\left[F_{*}(z)-\mathbbm{1}_{\\{z\geq
y_{*}\\}}\right]^{2}\,dz$
$\displaystyle=\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}}$
(16)

where $\phi(\cdot),\Phi(\cdot)$ are the standard Gaussian density and cdf.
Observe how CRPS penalizes 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 normalized 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:

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

##  4 Modeling Multiple Causes of Death

To understand the behavior of Age-specific mortality across different causes
of death, we begin by generating MOGP models for Cancer variants. Using the
HCD database we fit both ICM and SLFM and assess their performance in 3
European countries. We test the resulting predictive performance by computing
APE and CRPS for one-year-ahead mortality forecasts in three separate test
sets, using SOGP as the baseline. All models, including SOGP models, are
trained on the same Ages from 50–84 (7 age groups) and three overlapping
periods: 1998–2013 for the 2014 prediction, 1999–2014 for 2015 prediction, and
lastly 2000–2015 for 2016 prediction. We report APE and CRPS for all-Cancer
log-mortality rates since some of the cancer variations, such as Stomach and
Pancreatic have relatively few (and therefore more noisy) recorded deaths. The
differences in APE and CRPS between MOGP and SOGP models are expressed as the
3-year median percentage improvement over SOGP models. Positive improvement
means joint models have smaller mean APE/ CRPS values. To compute all-Cancer
CRPS, we leverage the closed-form expression of the MOGP multivariate
predictive distribution in (6)-(7) to simulate the forecast distribution of
all-Cancer log-mortality for each Age group in the data. To do so, we first
draw ($5\times 10^{5}$) stochastic samples of the joint
$\mathbf{f_{*}}(x_{*})$ across all the Cancer types. After exponentiating and
summing, we then obtain corresponding samples of (unlogged) all-Cancer
mortality rates.

Table 1: Comparison between MOGP ICM/SLFM with different ranks $Q$, reported as the relative improvement in APE and CRPS of MOGP vis-a-vis SOGP. The reported values are medians of one-year-out aggregated all-Cancer forecasts for age groups 50–84 based on 3 training periods: 1998–2013 (predict 2014), 1999–2014 (predict 2015), and 2000–2015 (predict 2016). Both ICM and SLFM are fitted on 5 Cancer types, Male populations in each selected country. Cause ($L=5$) |  # Kernel Hyperparameters | Czech Rep. | Germany | Poland
---|---|---|---|---
APE | CRPS | APE | CRPS | APE | CRPS
ICM | $Q=2$ | 12 | 19.16 | 30.18 | 3.86 | 15.05 | 12.85 | 17.72
$Q=3$ | 17 | 17.69 | 19.08 | 6.53 | 15.16 | 13.46 | 18.64
$Q=4$ | 22 | 22.31 | 21.43 | 6.30 | 14.58 | 17.02 | 11.01
SLFM | $Q=2$ | 14 | 18.04 | 29.49 | 5.09 | 13.51 | 10.82 | 16.05
$Q=3$ | 21 | 22.87 | 23.72 | 2.78 | 13.55 | 4.11 | 10.65
$Q=4$ | 28 | 16.97 | 18.24 | 8.75 | 16.73 | 5.13 | 3.81

Table 1 shows the 3-year median improvement in MAPE and CRPS for Multi-cause
ICM and SLFM over SOGP models. Overall, joint models produce more accurate
mean forecasts (positive improvement in APE) with higher credibility (positive
improvement in CRPS). We observe the opportunity to borrow information across
different cancers to better estimate the kernel hyperparameters. As a result,
joint models can describe important trends for individual cancers, leading to
the reduction in disparity between the predicted values and the observed all-
Cancer mortality in the test sets.

Figure 2: The predicted YoY improvement rates derived from multi-cause GP
models by Country and Age groups among Male observations. In each country,
MOGP-ICM ($Q=2$) and SLFM ($Q=2$) are fitted on Ages 50–84 (7 groups), Years
1998–2016 over 5 Cancer variations: Stomach, Colorectal, Pancreatic, Lung and
Remaining types.

###  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
modeling 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 (for example 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 3 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 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 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 aging
population and an increase in lifestyle exposure to risk factors for cancers
(Wojtyś et al. 2014). The incidence of Stomach cancer has been flat over time,
consistent with the finding in Arash et al. (2020).

(a) Multi-cause ICM in Czech R. (b) Multi-cause ICM in Germany (c) Multi-cause
ICM in Poland

(d) Multi-cause SLFM in Czech R. (e) Multi-cause SLFM in Germany (f) Multi-
cause SLFM in Poland

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 (7 groups) during Years 1998–2016, over 5 Cancer
variations: Stomach, Colorectal, Pancreatic, Lung and Remaining types.

Figure 3 visualizes the inferred cross-cause correlation matrices
$(r_{l_{1},l_{2}}:1\leq l_{1},l_{2}\leq 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,
e.g. Stomach and Pancreas cancers in Poland; this may be due to a competing
risks context.

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 8 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 14 in Appendix B. Overall, our results are in
agreement with the SOA study (Boumezoued et al. 2019), 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 12 in Appendix A displays the factor
loadings in a Country-Cause SLFM with $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 behavior
of the respective mortality surface.

###  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 visualizes 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.

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

Figure 4: Predictive distribution of all-Cancer log-mortality rates for
different age groups in 3 countries via multi-level Country-Cause ICM. The
joint model is fitted on Males, Age groups between 50–84, Years 1998–2016
across 3 Countries ($Q_{ctry}=3$) and 5 Cancers ($Q_{caus}=5$). Shading
indicates the 60%, 80%, 95% and 99% predictive quantile bands.

###  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 6 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), see the 2 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.
(2019) emphasized 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 are thus automatically
adjusted if trends shift. So for example, we do not need to do any break-point
analysis to achieve the slow pace of future HEA improvement shown in Figure 5.

AGE 40-44, US MALES

AGE 40–44, US FEMALES

AGE 65–69, US MALES

AGE 65–69, US FEMALES

Figure 5: Posterior distribution of true log-mortality rates for US Males and
Females for the 40-44 and 65-69 Age groups. In each row: left panel shows the
posterior quantiles for top 6 individual causes (via multi-cause GP); middle
panel shows aggregate log-mortality trends via by-cause model (multi-cause
GP); right panel shows aggregate log-mortality trends via all-cause model
(single-output GP). All models are fitted on Ages 40–69 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 modeling. The
pessimism of by-cause analysis was first mentioned in Wilmoth (1995) and re-
iterated in Boumezoued et al. (2019). 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, witness the wider band in the middle panels of Figure 4 relative
to the right panel based on a SOGP.

As discussed in Huynh and Ludkovski (2021), 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 fueled by the opioid epidemic. Assuming this crisis is addressed in
the future, the projected DRU mortality should be adjusted downward. In MOGP
this can be achieved by modifying the Year trend in $m(\cdot)$. Figure 13 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. 2019)
suggests that this pace might not take place in the intermediate term,
lowering aggregate mortality gains.

##  5 Cause-of-Death Joint Modeling in a Multinational Context

We proceed to consider simultaneously 30 populations in the Cancer variations
case study, arranged by the three factor inputs of Cause $(L_{caus}=5$),
Country ($L_{ctry}=3$) and Gender ($L_{gndr}=2$).

###  5.1 Multi-level vs Single-level Correlation Structure

(a) Multi-level ICM ($Q_{ctry}=3$, $Q_{caus}=5$, $Q_{gndr}=2$) with 40
hyperparameters (b) Single-level ICM ($Q=2$) with 62 hyperparameters

Figure 6: Cross-correlation matrices of MOGP models that incorporate Country,
Cause, and Gender as categorical inputs. Thirty total populations.

Figure 6 displays the inferred correlation structure $r_{l_{1},l_{2}}:1\leq
l_{1},l_{2}\leq 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
1998–2016 for 3 Countries, 5 Cancers, and 2 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 $QL+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.

(a) Cross-country correlation ($R_{C}$) (b) Cross-cause correlation ($R_{S}$)
(c) Cross-gender correlation ($R_{G}$)

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 and Years 1998–2016.

###  5.2 Model Selection

Next, we compare the predictive performance of multi-level vs single-level
MOGP-ICM. Following the same setup 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 has 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
structural break encountered in mid-2010s (see Fig. 4), the performance of
aggregated Country-Cause MOGP is consistently worse than those of all-cause
SOGP; this issue is rectified for Country-Cause-Gender MOGPs.

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 vs SOGP of one-year-out aggregated all-Cancer forecasts for age groups
50–84 based on 3 training periods: 1998–2013 (predict 2014), 1999–2014
(predict 2015), and 2000–2015 (predict 2016). Top half: MOGP models fitted on
5 cancer types and 3 countries. Bottom half: MOGP models fitted on 5 cancer
types, 3 countries and 2 genders.

Country$+$Cause ($L=15$) |  # Kernel Hyperparameters | Czech Rep. | Germany | Poland
---|---|---|---|---
APE | CRPS | APE | CRPS | APE | CRPS
ICM | $Q=2$ | 32 | 35.54 | 30.74 | 12.02 | 21.06 |  $-$40.44 |  $-$59.47
$Q=3$ | 47 | 45.51 | 50.07 | 8.50 | 20.35 |  $-$37.09 |  $-$38.48
$Q=4$ | 62 | 31.64 | 36.56 | 14.08 | 35.66 |  $-$13.32 | 9.56
SLFM | $Q=2$ | 34 | 31.51 | 24.57 | 4.44 | 15.87 |  $-$47.66 |  $-$45.98
$Q=3$ | 51 | 33.78 | 41.59 | 0.88 | 8.62 |  $-$37.25 |  $-$24.54
$Q=4$ | 68 | 44.79 | 43.29 | 8.93 | 20.49 |  $-$45.63 |  $-$37.66
Multi-level ICM | $Q_{ctry}=2,~{}Q_{caus}=2$ | 18 | 26.42 | 24.24 | 4.55 | 6.84 |  $-$48.46 |  $-$61.53
$Q_{ctry}=2,~{}Q_{caus}=4$ | 28 | 25.74 | 32.26 | 5.04 | 15.75 |  $-$35.94 |  $-$44.24
$Q_{ctry}=3,~{}Q_{caus}=4$ | 36 | 42.03 | 36.80 | 3.07 | 11.12 |  $-$44.45 |  $-$53.69
Country$+$Cause$+$Gender ($L=30$) |  # Kernel Hyperparameters | Czech Rep. | Germany | Poland
APE | CRPS | APE | CRPS | APE | CRPS
ICM | $Q=2$ | 62 | 9.99 | 11.75 |  $-$5.40 |  $-$13.40 | 12.21 | 6.55
$Q=3$ | 92 | 0.52 | 3.07 |  $-$10.78 |  $-$20.35 |  $-$15.31 | 4.09
Multi-level ICM | $Q_{ctry}=2,~{}Q_{caus}=4,~{}Q_{gndr}=2$ | 32 | 14.85 | 24.91 |  $-$68.45 |  $-$15.22 |  $-$25.11 | 7.75
$Q_{ctry}=3,~{}Q_{caus}=5,~{}Q_{gndr}=2$ | 40 | 21.34 | 27.22 |  $-$40.45 |  $-$10.22 | 5.12 | 8.04

STOMACH CANCER

COLORECTAL CANCER

PANCREATIC CANCER

LUNG CANCER

Figure 8: Predictive log-mortality distributions from single-level ICM ($Q=2$)
and multi-level ICM ($Q_{ctry}=3,~{}Q_{caus}=5,~{}Q_{gndr}=2$) models with 3
factor inputs: Country, Cause, and Gender (30 populations). All models are
fitted on Age groups 50–84 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.

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 upon 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.

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 setups $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 via full-rank multi-level ICM and single-level ICM with $Q=2$.
Both models are fitted on Age groups 50–84 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 length-scale in Year in multi-level
ICM ($\theta_{yr}\approx 8.8$ versus $\theta_{yr}\approx 14.8$ in single-level
ICM).

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 10. This offers opportunities to fuse data
from other countries to update domestic forecasts. For Age-Period-Cohort
models, such as Li and Lee (2005), this is generally challenging as the 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 choose Czech Males as a reference population and
examine one-year-out prediction quality for several models. We take Czech
observations for 1998–2015 and borrow the more up-to-date dataset from Germany
(1998–2016) to implement the notched setup. Table 4 in the Appendix describes
the setup of the different models we consider in this experiment.

Figure 9: Comparison of the prediction accuracy for 2016 all-Cancer log-
mortality of Czech Males for indicated Age groups between different models
with ‘notched’ setup. Top row: Predictive standard deviation $s_{*}(x_{*})$;
bottom row: discrepancy between the predictive mean $m_{*}(x_{*})$ and the
observed value $y_{*}(x_{*})$.

Figure 9 visualizes relative performance by comparing two components: (a)
predictive standard deviation of the underlying 2016 $s_{*}(x_{*})$ (top
panels), and (b) prediction errors $m_{*}(x_{*})-y_{*}(x_{*})$ relative to the
realized 2016 all-Cancer log-mortality rates of Czech Males (bottom). Our
benchmark is Czech Males SOGP fitted on all-Cancer log-mortality rates from
1998–2016. We see that a Country-Cause MOGP yields more credible predictions
(lower $s_{*}(x_{*})$) than the benchmark model that performs in-sampled
smoothing of the 2016 Czech all-Cancer experience. The forecasts from this
joint model are nearly as good as the benchmark as both models produce the
smallest prediction errors. In fact, the prediction quality of Country-Cause
MOGP is as competitive as the Multi-cause MOGP that simultaneously models the
log-mortality rates of all 5 cancer types in Czech Males with 2016
observations available. Thus, borrowing the latest information from Germany
improves the prediction quality of the recent all-Cancer mortality rates of
Czech Males. The results are as good as when we have domestic data available
in Czech Republic.

Remark: At the moment, the HCD offers cause-of-death data for less than 15
countries. Many countries do not have recent data available yet (e.g. only up
to 2014), leading to limited options in terms of the selection and the number
of countries we can incorporate with Czech Males in this experiment. Based on
our analysis in Huynh and Ludkovski (2021), choosing countries that are highly
correlated with Czech Males is essential to maximize the prediction quality.

##  6 Conclusion

In this article, we develop multi-output GP models for cause-specific
mortality modeling within a multi-population context. With the MOGP mechanism,
we are able to capture the cross-cause dependencies that allow joint models to
gain more predictive power over single-output models that treat each
population independently. Among MOGP variants, SLFM offers more flexibility
and is recommended for modeling heterogenous causes, such as top-level ICD
categories. Multi-level ICM is demonstrated to work well for interpretable
modeling across multiple factor inputs. Our case studies show the
applicability of MOGPs to understand cause-specific and aggregate mortality
trends, both within a country and across nations, whereby our framework is
convenient for information fusion and credibility boosting.

The current work focuses on exploiting the structured Kronecker covariance to
efficiently learn the joint covariance kernel. This is sufficient for handling
a moderately large number of sub-populations (up to 30 in our case studies);
additional techniques would be needed to handle larger datasets, for example
across more causes-of-death or across all the countries in HCD. There is an
active research area looking at alternative methods (local kernel
interpolation, inducing points, block structures, etc., see e.g., Flaxman et
al. (2015)) for massive scalable GP well-suited for gridded mortality
datasets. Another methodological extension would be to consider a Linear
Coregionalization Model (LCM) for mortality modeling. LCM generalizes ICM and
SLFM and allows multiple latent functions from GP priors with different
covariance kernels. A third direction would be to investigate other kernel
families, such as composite kernels or kernels that can incorporate structural
changes. The latter would be useful to model (sub-)causes that exhibit strong
disruptions over time in their mortality trends.

##  Appendix A Additional Plots for the Case Study on Cancer Sub-Types

### Year ranges available in the HCD by country

Figure 10: Countries in the HCD and their historical data coverage.

### Age patterns of log-mortality rates of different cancers

Figure 11: Raw log mortality as a function of Age for five Cancer variations
in Male populations. We show 19 faint curves for each year in 1998–2016 range
used for training, as well as the respective average historical log-mortality
rate in bold.

### Kernel Hyperparameter Learning

Searching for optimal hyperparameters in GP can be challenging if the marginal
likelihood features multiple local maxima or flattens around its global
maximum (Rasmussen and Williams 2005). When the optimizer fails to find the
global maximum, unsuitable length-scales in Age and Year sometimes result,
leading to a poor fit of the data. Table 3 reports the inferred length-scale
in Age ($\theta_{ag}$) and Year ($\theta_{yr}$) of the SOGP models fitted on
the 5 Cancer types for the Male populations in the three considered countries.
Many SOGP models have $\theta_{yr}$ being too short (less than 5 or so),
resulting in oscillatory fitted $m_{*}(\cdot)$’s. Such models lack the ability
to distinguish between true signal and the inherent randomness in the data.
Similarly, when the estimated lengthscales are too large, the fitted GP
surfaces are over-smoothed. Joint models tend to better learn the hyper-
parameters by enabling data fusion across multiple populations and utilizing
more observations. In Table 3 we show that when we fit multi-cause MOGP (both
ICM and SLFM) on all 5 cancer types, the resulting lengthscales are all well
calibrated.

Table 3: GP lengthscales $\theta_{ag},\theta_{yr}$ in Age and Year for
different models. All models are fitted on Age groups 50–84 and Years
1998–2016.

Czech Rep. | SOGP on each cancer type | Multi-cause ICM | Multi-cause SLFM
---|---|---|---
Stomach | Colorectal | Pancreas | Lung | Others | ($Q=2$) | ($Q=2$)
$\theta_{ag}$ | 159.92 | 16.32 | 6.21 | 17.69 | 12.52 | 16.90 | 21.26 | 19.67
$\theta_{yr}$ | 38.49 | 9.03 | 7.82 | 13.87 | 3.69 | 10.12 | 13.05 | 11.42

Germany | SOGP on each cancer type | Multi-cause ICM | Multi-cause SLFM
---|---|---|---
Stomach | Colorectal | Pancreas | Lung | Others | ($Q=2$) | ($Q=2$)
$\theta_{ag}$ | 8.65 | 3.72 | 8.99 | 0.00 | 2.79 | 8.59 | 11.75 | 9.46
$\theta_{yr}$ | 7.35 | 4.90 | 7.69 | 4.82 | 4.44 | 10.26 | 9.92 | 6.77

Poland | SOGP on each cancer type | Multi-cause ICM | Multi-cause SLFM
---|---|---|---
Stomach | Colorectal | Pancreas | Lung | Others | ($Q=2$) | ($Q=2$)
$\theta_{ag}$ | 23.96 | 16.28 | 5.72 | 15.25 | 22.80 | 16.98 | 16.20 | 21.89
$\theta_{yr}$ | 253.34 | 3.69 | 6.19 | 9.04 | 6.01 | 12.83 | 11.19 | 12.11

### Training Designs in Notched Datasets

Table 4: Descriptions for models being applied to forecast 2016 all-Cancer
log-mortality rates of Czech Males. All the models are fitted on Ages 50–84 (7
age groups) and Year coverage listed below.

GP models | Outcome variable | Country | Year | Pred. Type | Abbrev.
---|---|---|---|---|---
SOGP | All-cancer log-mortality | Czech R. | 1998–2016 | In-sample | SOGP (’98-’16)
1998–2015 | Out-of-sample | SOGP (’98-’15)
Multi-cause GP |  By-cancer log-mortality (5 variations) | Czech R. | 1998–2016 | In-sample | Multi-cause (’98-’16)
1998–2015 | Out-of-sample | Multi-cause (’98-’15)
Country-Cause GP |  By-cancer log-mortality (5 variations) | Czech R. | 1998-2015 | Out-of-sample | Country-Cause
Germany | 1998-2016

### Illustrating Latent Factor Loadings in SLFM

Figure 12: Factor loadings $\mathbf{a}_{q}=(a_{1,q},\ldots,a_{L,q})^{T}$ in
the Country-Cause SLFM with $Q=3$. The model is fitted on Ages 50–84, Years
1990–2016, over 3 Countries and 5 Cancer types. For each of the 15
populations, we plot $(a_{l,1},a_{l,2},a_{l,3})$ as a point in 3-space.

##  Appendix B Additional Plots for the US Top-Level-Cause Study

### Adjusting Drug Overdose Trend

Age Group 40–44. Top row: US Males $\uparrow$. Bottom row: US Females
$\downarrow$

Figure 13: Comparison of original and adjusted predicted log-mortality for Age
group 40–44, US Males and Females. Left: same model after reducing Drug Abuse
yearly trend $\beta_{yr}$ by one-third. Right: corresponding aggregated all-
cause trends. Vertical lines indicate the edge of training data (1999–2018).

### Cross-correlation Matrices for the US All-Cause Analysis

(a) US Males (b) US Females

Figure 14: Cross-cause correlation matrices derived from multi-cause single-
level ($Q=6$) ICM MOGP, fitted on Ages 40–69 (6 age groups) and Years
1999–2018, separately for US Males and US Females.

## References

  * Alai et al. [2018] Daniel H. Alai, Séverine Arnold (-Gaille), Madhavi Bajekal, and Andrés M. Villegas.  Mind the gap: a study of cause-specific mortality by socioeconomic circumstances.  _North American Actuarial Journal_ , 22(2):161–181, 2018.
  * Alvarez et al. [2012] Mauricio A. Alvarez, Lorenzo Rosasco, and Neil D. Lawrence.  Kernels for vector-valued functions: a review.  _Foundations and Trends in Machine Learning_ , 4(3):195–266, 2012.
  * Anderson et al. [2001] Robert N Anderson, Arialdi M Miniño, Donna L Hoyert, and Harry M Rosenberg.  Comparability of cause of death between ICD-9 and ICD-10: preliminary estimates.  _National Vital Statistics Report_ , 49(2):1–32, May 2001.
  * Arash et al. [2020] Etemadi Arash, Shakeri Ramin, Safiri Saeid, and Sadaf G Sepanlou.  The global, regional, and national burden of stomach cancer in 195 countries, 1990–2017: a systematic analysis for the global burden of disease study 2017.  _The Lancet Gastroenterology and Hepatology_ , 5(1):42–54, 2020.
  * Arnold and Sherris(2013) [-Gaille] Séverine Arnold (-Gaille) and Michael Sherris.  Forecasting mortality trends allowing for cause-of-death mortality dependence.  _North American Actuarial Journal_ , 17(4):273–282, 2013.
  * Bergeron-Boucher et al. [2017] Marie-Pier Bergeron-Boucher, Vladimir Canudas-Romo, James E. Oeppen, and James Vaupel.  Coherent forecasts of mortality with compositional data analysis.  _Demographic Research_ , 37(17):527–566, 2017\.
  * Bonilla et al. [2008] Edwin V Bonilla, Kian M. Chai, and Christopher Williams.  Multi-task Gaussian Process prediction.  In _Advances in Neural Information Processing Systems 20_ , pages 153–160. Curran Associates, Inc., 2008.
  * Boumezoued et al. [2019] Alexandre Boumezoued, Jean-Baptiste Coulomb, Al Klein, Damien Louvet, and Eve Titon.  Modeling and forecasting cause-of-death mortality.  Technical report, Society of Actuaries, 2019.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matthew D.Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: a probabilistic programming language.  _Journal of Statistical Software_ , 76(1), 2017.
  * Caruana [1997] Rich Caruana.  Multi-task learning.  _Machine Learning_ , 28(1):41–75, Jul 1997.
  * Caselli [1996] Graziella Caselli.  Future longevity among elderly populations.  In _Health and mortality among elderly populations_ , pages 235––265. Oxford University Press, 1996.
  * Caselli et al. [2019] Graziella Caselli, Jacques Vallin, and Marco Marsili.  How useful are the causes of death when extrapolating mortality trends. an update.  In Tommy Bengtsson and Nico Keilman, editors, _Old and New Perspectives on Mortality Forecasting_ , pages 237–259. Springer International Publishing, 2019.
  * Deville et al. [2019] Yves Deville, David Ginsbourger, and Olivier Roustant. Contributors: Nicolas Durrande.  _kergp: Gaussian Process Laboratory_ , 2019.  URL https://CRAN.R-project.org/package=kergp.  R package version 0.5.0.
  * Dimitrova et al. [2013] Dimitrina S. Dimitrova, Steven Haberman, and Vladimir K. Kaishev.  Dependent competing risks: cause elimination and its impact on survival.  _Insurance: Mathematics and Economics_ , 53(2):464–477, 2013.
  * Dong et al. [2020] Yumo Dong, Fei Huang, Honglin Yu, and Steven Haberman.  Multi-population mortality forecasting using tensor decomposition.  _Scandinavian Actuarial Journal_ , 2020(8):754–775, 2020.
  * Enchev et al. [2017] Vasil Enchev, Torsten Kleinow, and Andrew J. G. Cairns.  Multi-population mortality models: fitting, forecasting and comparisons.  _Scandinavian Actuarial Journal_ , 2017(4):319–342, 2017.
  * Flaxman et al. [2015] Seth Flaxman, Andrew Wilson, Daniel Neill, Hannes Nickisch, and Alex Smola.  Fast Kronecker inference in Gaussian processes with Non-Gaussian likelihoods.  In Francis Bach and David Blei, editors, _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pages 607–616, Lille, France, 07–09 Jul 2015. PMLR.
  * Foreman et al. [2018] Kyle Foreman, N. Marquez, A. Dolgert, K. Fukutaki, N. Fullman, M. McGaughey, Martin A Pletcher, A. Smith, Kendrick Tang, Chun-Wei Yuan, J. Brown, Joseph Friedman, J. He, Kyle Heuton, Mollie Holmberg, Disha J Patel, P. Reidy, Austin Carter, Kelly M Cercy, Abigail Chapin, D. Douwes-Schultz, Tahvi D Frank, F. Goettsch, P. Liu, Vishnu Nandakumar, M. Reitsma, Vincent Reuter, Nafis Sadat, Reed J. D. Sorensen, V. Srinivasan, R. Updike, Hunter York, Alan D. Lopez, R. Lozano, Stephen S. Lim, A. Mokdad, S. Vollset, and C. Murray.  Forecasting life expectancy, years of life lost, and all-cause and cause-specific mortality for 250 causes of death: reference and alternative scenarios for 2016–40 for 195 countries and territories.  _Lancet (London, England)_ , 392:2052 – 2090, 2018.
  * Gilboa et al. [2015] Elad Gilboa, Yunus Saatçi, and John P. Cunningham.  Scaling multidimensional inference for structured Gaussian processes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 37(2):424–436, 2015.
  * Guibert et al. [2019] Quentin Guibert, Olivier Lopez, and Pierrick Piette.  Forecasting mortality rate improvements with a high-dimensional VAR.  _Insurance: Mathematics and Economics_ , 88:255–272, 2019\.
  * HCD [2021] HCD.  The Human Cause-of-Death Database. French Institute for Demographic Studies (France) and Max Planck Institute for Demographic Research (Germany), 2021.  URL www.causeofdeath.org.
  * Huynh and Ludkovski [2021] Nhan Huynh and Mike Ludkovski.  Multi-output Gaussian processes for multi-population longevity modelling.  _Annals of Actuarial Science_ , 15(2):318–345, 2021.
  * Huynh et al. [2020] Nhan Huynh, Mike Ludkovski, and Howard Zail.  Multi-population longevity models: a spatial random field approach.  _Proceedings of the Society of Actuaries 2020 Living to 100 Symposium_ , 2020.
  * Hyndman et al. [2013] Rob J Hyndman, Heather Booth, and Farah Yasmeen.  Coherent mortality forecasting: the product-ratio method with functional time series models.  _Demography_ , 50(1):261–283, 2013.
  * Kjaergaard et al. [2019] Søren Kjaergaard, Yunus Emre Ergemen, Malene Kallestrup-Lamb, Jim Oeppen, and Rune Lindahl-Jacobsen.  Forecasting causes of death by using compositional data analysis: the case of cancer deaths.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 68(5):1351–1370, 2019.
  * Kleinow [2015] Torsten Kleinow.  A common age effect model for the mortality of multiple populations.  _Insurance: Mathematics and Economics_ , 63:147–152, 2015\.
  * Knudsen and McNown [1993] Christin Knudsen and Robert McNown.  Changing causes of death and the sex differential in the USA: Recent trends and projections.  _Population Research and Policy Review_ , 12(1):27–41, 1993.
  * Letham and Bakshy [2019] Benjamin Letham and Eytan Bakshy.  Bayesian optimization for policy search via online-offline experimentation.  _Journal of Machine Learning Research_ , 20(145):1–30, 2019.
  * Li and Lu [2017] Hong Li and Yang Lu.  Coherent forecasting of mortality rates: A sparse vector-autoregression approach.  _ASTIN Bulletin: The Journal of the IAA_ , 47(2):563–600, 2017.
  * Li and Lu [2019] Hong Li and Yang Lu.  Modeling cause-of-death mortality using hierarchical Archimedean copula.  _Scandinavian Actuarial Journal_ , 2019(3):247–272, 2019.
  * Li and Lee [2005] Nan Li and Ronald Lee.  Coherent mortality forecasts for a group of populations: An extension of the Lee-Carter method.  _Demography_ , 42(3):575–594, 2005.
  * Liu et al. [2019] Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai.  When Gaussian process meets big data: A review of scalable GPs, 2019\.  Arxiv 1807.01065.
  * Lo and Wilke [2010] Simon M. S. Lo and Ralf A. Wilke.  A copula model for dependent competing risks.  _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , 59(2):359–376, 2010.
  * Ludkovski et al. [2018] Mike Ludkovski, Jimmy Risk, and Howard Zail.  Gaussian process models for mortality rates and improvement factors.  _ASTIN Bulletin: The Journal of the IAA_ , 48(3):1307–1347, 2018.
  * Lyu et al. [2021] Pintao Lyu, Anja De Waegenaere, and Bertrand Melenberg.  A multi-population approach to forecasting all-cause mortality using cause-of-death mortality data.  _North American Actuarial Journal_ , 25(1):S421–S456, 2021.
  * McNown and Rogers [1992] Robert McNown and Andrei Rogers.  Forecasting cause-specific mortality using time series methods.  _International Journal of Forecasting_ , 8(3):413–432, 1992.
  * Rasmussen and Williams [2005] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)_.  The MIT Press, 2005.
  * Saatçi [2011] Yunus Saatçi.  _Scalable Inference for Structured Gaussian Process Models_.  PhD thesis, University of Cambridge, 2011.
  * Tabeau et al. [1999] Ewa Tabeau, Peter Ekamper, Corina Huisman, and Alinda Bosch.  Improving overall mortality forecasts by analysing cause-of-death, period and cohort effects in trends.  _European Journal of Population / Revue Européenne de Démographie_ , 15(2):153–183, 1999.
  * Teh et al. [2005] Y. Teh, M. Seeger, and Michael I. Jordan.  Semiparametric latent factor models.  In _International Workshop on Artificial Intelligence and Statistics (AISTATS)_ , pages 333–340. PMLR, 2005.
  * Tsai and Zhang [2019] Cary Chi-Liang Tsai and Ying Zhang.  A multi-dimensional Bühlmann credibility approach to modeling multi-population mortality rates.  _Scandinavian Actuarial Journal_ , 2019(5):406–431, 2019.
  * WHO [2021] WHO.  Cancer, 2021.  URL https://www.who.int/news-room/fact-sheets/detail/cancer.
  * Williams et al. [2009] Christopher Williams, Stefan Klanke, Sethu Vijayakumar, and Kian Chai.  Multi-task Gaussian process learning of robot inverse dynamics.  In D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, editors, _Advances in Neural Information Processing Systems_ , volume 21. Curran Associates, Inc., 2009.
  * Wilmoth [1995] John R. Wilmoth.  Are mortality projections always more pessimistic when disaggregated by cause of death?  _Mathematical Population Studies_ , 5(4):293–319, 1995.
  * Wojtyś et al. [2014] Piotr Wojtyś, Andrzej Antczak, and Dariusz Godlewski.  Predictions of cancer mortality in Poland in 2020.  _Central European Journal of Medicine_ , 9:667–675, 2014\.
  * Zhe et al. [2019] Shandian Zhe, Wei Xing, and Robert M. Kirby.  Scalable high-order Gaussian process regression.  In _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , volume 89 of _Proceedings of Machine Learning Research_ , pages 2611–2620. PMLR, 16–18 Apr 2019.