2.1 Dataset

Built by the HMD team, the United States Mortality Database (USMDB) United States Mortality Database (USMDB) contains a complete historical set of state-level life tables for every calendar year during 1959–2019 for all 50 U.S. states and the District of Columbia (D.C.). The USMDB covers ages 0–110+ and includes separate datasets for the Male and Female populations. The raw data contain birth and death counts from the U.S. vital statistics system and incorporate the census counts and population estimates from the U.S. Census Bureau.

Our objective is to estimate and smooth historical mortality rates and then forecast short-term calendar trends through analyzing and estimating mortality improvement factors. Throughout this paper, we thus focus on the following subset of the USMDB: (a) Calendar Years 1990–2018; (b) Males and Females, considered separately; and (c) Ages 60–84. This subset of older ages and recent years is the most relevant for actuarial applications. Omitting data from the 20th century is in line with the data-driven machine-learning framework we employ, so that long-past mortality experience is not only less relevant but potentially misleading for our model construction, in particular due to nonstationary mortality improvement patterns. We omit very old Ages since the population data underlying USMDB are only available up to an open age interval at 85+ years (Barbieri, Reference Barbieri2020).

Due to operating on raw data, nearly all of the works cited above consider aggregated mortality across Age groups, e.g., bins of 10 years (Becker et al., Reference Becker, Majmundar and Harris2021), or 5 years (Dwyer-Lindgren et al., Reference Dwyer-Lindgren, Bertozzi-Villa, Stubbs, Morozoff, Kutz, Huynh, Barber, Shackelford, Mackenbach, van Lenthe and Murray1980; Shull et al., Reference Shull, Richardson, Groendyke and Hartman2025; Vierboom et al., Reference Vierboom, Preston and Hendi2019), or look at life expectancy at a given age ( $e_0, e_{50}, e_{65}$ , etc.) (Barbieri, Reference Barbieri2020; Boing et al., Reference Boing, Boing, Cordes, Kim and Subramanian2020; Chetty et al., Reference Chetty, Stepner, Abraham, Lin, Scuderi, Turner, Bergeron and Cutler2016; Vierboom & Preston, Reference Vierboom and Preston2020; Wilmoth et al., Reference Wilmoth, Boe and Barbieri2011). This fundamentally obscures the non-constant relative impact of Age, which is well documented. Indeed, multiple studies conclude that spatial and socio-economic inequalities decline with age: “health disparities narrow with age” (Vierboom et al., Reference Vierboom, Preston and Hendi2019) and “the gap [in probabilities of dying by SIS deciles] declines progressively after age 55 years and becomes small (less than 10 percent) at ages 85 and above” (Barbieri, Reference Barbieri2020). At the same time, most references disaggregate mortality by other, say socio-economic or cause-of-death, factors, or into smaller spatial units. In contrast, we follow the USMDB to fully disaggregate into 1-year bins by Age, but otherwise consider aggregated mortality across all individuals in the state.

Let $\mathcal{S} = (s_i)_{1 \le i \le 51}$ represent the 51 U.S. states (throughout we count D.C. as the 51st “state”). The information regarding each state $s \in \mathcal{S}$ provided by the USMDB is organized as follows:

1. (i) Independent variables: calendar year $x_{t}$ and age $x_{a}$ . The pair $(x_{a},x_t)$ refers to the set of persons from state $s$ aged $a \in \{60,\ldots ,84\}$ during year $t \in \{1990,\ldots , 2018\}$ .

2. (ii) Dependent variables: The death counts $D^{(a,t)}$ along with the total number of persons lived (exposed to risk) $E^{(a,t)}$ at a given age-time cell $(x_{a},x_t)$ . We record the log-mortality rate

   (1) \begin{equation} y^{(a,t)} \;:\!=\; \log \left [ \frac {\text{$\#$ of Deaths during $(x_{a},x_t)$ age-time Interval}}{\text{$\#$ of Exposed to Risk during $(x_{a},x_t)$ age-time Interval}}\right ] \equiv \log \left [ \frac {D^{(a,t)}}{E^{(a,t)}} \right ]. \end{equation}
   To simplify notation, when the age and calendar year are clear from the context, we drop the superscript $(a,t)$ . While death counts are generally highly accurate, exposed-to-risk are based on census estimates and systematically undercount undocumented immigrants.

In summary, the USMDB provides 5 inputs, $(s,x_a,x_t,D,E) \equiv$ (state, age, year, deaths, exposed-to-risk) and the associated output $y \equiv$ log-mortality, see Table 1. The complete dataset is denoted as $\mathcal{D}$ , with $\mathcal{D}_s \subset \mathcal{D}$ representing the subset of rows associated with state $s \in \mathcal{S}$ .

Table 1. Sample rows from the USMDB dataset $\mathcal{D}$

2.2 Multi-population Gaussian process models

Multi-population modeling aims to identify and capture mortality dependence patterns among several populations in order to fuse data and achieve coherent forecasts. We follow Huynh et al. (Reference Huynh, Ludkovski and Zail2020) and Huynh & Ludkovski (Reference Huynh and Ludkovski2021) in implementing a Multi-Output Gaussian Process (MOGP) to model USMDB longevity data. The MOGP model quantifies mortality uncertainty by probabilistically smoothing raw data and simultaneously generates stochastic out-of-sample forecasts by projecting mortality surfaces across the Age and Year dimensions. For multi-population analysis, MOGP imposes a cross-population correlation structure on top of the Age-Period pattern in each population. The data fusion employed in MOGP improves model fit, reduces model risk, and provides insights into the discrepancies among mortality trends across populations.

2.2.1 Gaussian process regression

First, we describe the mechanics of the Single-Output Gaussian Process (SOGP) model; see Ludkovski et al. (Reference Ludkovski, Risk and Zail2018) for more details. Fix an arbitrary state $l \in \mathcal{S}$ . We are given a sample of $n = 29 \times 25 = 725$ observed log-mortality rates. Mortality is described as a function of age and time:

• • $\mathbf{x} \;:\!=\; (x^1, \ldots , x^n)$ where $x^i \;:\!=\; (x_a^i,x_t^i)$ for $a \in \{60, \ldots , 84\}$ and $t \in \{ 1990,\ldots ,2018\}$ .

• • $\mathbf{y}_l(\mathbf{x}) \equiv \mathbf{y}_l \;:\!=\; (y_l^1, \ldots , y_l^n)$ where $y^i_l, i=1,\ldots , n$ denotes the observed log-mortality rate provided by the USMDB for a state $l \in \mathcal{S}$ at age $x_a^i$ during year $x_t^i$ .

Observation Likelihood. We assume that the relationship between $y_l^i$ and $x^i$ can be described with a latent black box function $f_l(\cdot )$ and white noise term,

(2) \begin{equation} y_l^i = f_l(x^i) + \epsilon _l^i; \qquad i =1, \ldots , n. \end{equation}

The observation noise is Gaussian with $\boldsymbol{\epsilon }_l = (\epsilon _l^1,\ldots , \epsilon _l^n) \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma }_l \;:\!=\; \textrm {diag}(\sigma _l^2))$ . We assume that the observation likelihood is population-dependent but constant in age and year, $\sigma _l = \text{StDev}(\epsilon ^i_l) \; \forall i$ . The underlying function $f_l(x^i)$ intuitively represents the true mortality rate which would materialize in the absence of random idiosyncratic shocks. While a Poisson observation likelihood is appropriate for small populations, typical death counts for the considered age groups are at least in the high hundreds, making the Gaussian approximation (which is methodologically much preferred for GPs) highly accurate.

Distribution of the Latent Process. A priori, for any sample $\mathbf{x}$ of $m \ge 1$ observations, the finite dimensional distributions (fdds) of $f_l(\mathbf{x}) = (f_l(x^1), \ldots , f_l(x^m))$ are postulated to follow a multivariate Gaussian law $f_l \sim \mathcal{GP}\big (m_l, C_l\big )$ with prior (parametric) mean function, $m_l(\cdot )$ , and covariance matrix $C_l(\cdot ,\cdot )$ ,

\begin{equation*} m_l(\mathbf{x}) \;:\!=\; \mathbb{E}[f_l(\mathbf{x})] =\big ( \mu _l(x^1), \ldots , \mu _l(x^m) \big ) \quad \text{and} \quad C_l({x}, {x}') \;:\!=\; \mathbb{E}\Big [\big (f_l({x}) - m_l({x})\big ) \big (f_l({x}') - m_l({x}')\big ) \Big ]. \end{equation*}

Assuming that $\epsilon ^i_l$ s are independent across $x^i$ s and from $f$ , it follows that

\begin{equation*} \mathbf{y}_l \sim \mathcal{N}\left (m_l, C_l + \boldsymbol{\Sigma }_l \right ), \qquad \boldsymbol{\Sigma }_l = \sigma _l^2 \boldsymbol{I}, \end{equation*}

since Cov $(y_l^i,y_l^j) = \text{Cov}\big (f_l(x^i), f_l(x^j)\big ) + \sigma _l^2 \delta (x^i,x^j)$ , where $\delta (x^i,x^j)$ is the Kronecker delta.

Mean and Covariance Structure. We assume functional representations for the mean and covariance functions $m_l(\cdot )$ , $C_l(\cdot ,\cdot )$ which represent the prior beliefs about the dataset. Recall that all the observed properties of a stochastic process with Gaussian fdds are characterized by $m_l(\cdot )$ and $C_l(\cdot ,\cdot )$ .

1. (i) The GP Mean function describes the prior trend in log-mortality rates. We use a parametric prior mean function, $m_l(x^i) = \beta _{0,l} + \sum _{j=1}^p \beta _{j,l}h_j({x})$ , where $h_j({x})$ s are given basis functions and the $\beta _{j,l}$ s are unknown coefficients to be estimated. Letting $\boldsymbol{\beta }_l = \big ( \beta _{0,l}, \ldots , \beta _{p,l})^T$ , $\boldsymbol{h}(x) = \big ( h_1(x), \ldots , h_p(x) \big )$ , we use the shorthand $m_l(x) = \boldsymbol{h}(x) \boldsymbol{\beta }_\ell$ . Below, we postulate a linear trend in the Age dimension:

   (3) \begin{equation} m_l(x^i) = \beta _{0,l} +\beta _{1,l}^a \cdot x^i_a. \end{equation}
   The choice (3) is used to de-trend the data according to an exponential increase in mortality as a function of Age (the so-called Gompertz Law of Mortality) in our segment of interest $x_a \in \{60, \ldots , 84\}$ .

2. (ii) The GP Covariance kernel captures the dependence of the response surface $f_l$ on the varying Age and Year dimensions $x_a, x_t$ . The GP kernel characterizes the smoothing process by quantifying the influence of inputs on the likelihood of the output. Our kernels are distance-based, capturing the logic that the mortality experience should be similar at neighboring data points, and separable across the Age and Period coordinates.

   We concentrate on a common family of covariance functions known as the Matérn class, equipped with automatic relevance determination. The Matérn-5/2 kernel defines the covariance between arbitrary univariate inputs $x,x_* \in \mathbb{R}$ as:

   (4) \begin{equation} C^{(M52)}(x,x_*; \theta ) \;:\!=\; \left ( 1 + \frac {\sqrt {5}}{\theta }|x- x_{*}| + \frac {5}{3 \theta ^2}|x-x_{*}|^2\right ) \cdot \exp \left \{-\frac {\sqrt {5}}{\theta }|x-x_{*}| \right \}. \end{equation}

   This kernel is parameterized by the length scale (hyper)parameter $\theta$ , to be estimated. To construct the overall dependence structure we use a multiplicative Age-Period-Cohort (APC) structure, so that the covariance between two mortality table entries $x^i \equiv (x_a^i,x_t^i)$ , $x^{j} \equiv (x_{a}^j,x_{t}^j)$ is

   (5) \begin{align} C^{(APC)}_l(x^i,x^j) \;:\!=\; \eta ^2 \cdot C^{(M52)}(x^i_a,x^j_{a}; \theta _{l,a}) \cdot C^{(M52)}(x^i_t,x^j_{t}; \theta _{l,t}) \cdot C^{(M52)}(x^i_c,x^j_{c}; \theta _{l,c}), \end{align}
   where $x_c \;:\!=\; x_t - x_a$ is the Birth Cohort (i.e., year of birth of an individual who is $x_a$ -old in year $x_t$ ), and $\eta ^2$ is the process variance hyperparameter, scaling covariances to capture the typical amplitude of the response. Observe that there are three length scale hyperparameters $\theta _{l,a}, \theta _{l,t}, \theta _{l,c}$ , jointly estimated together. The last cohort term in (5) is essential to capture well-known generational effects, such as the special 1918 and 1939 cohorts. The product structure is analogous to the Age-times-Year terms in the classical Lee-Carter framework. See Ludkovski & Risk (Reference Ludkovski and Risk2024) for a further discussion of appropriate GP kernels for mortality.

2.2.2 GP posterior

The GP paradigm models input-output relationships by algebraically conditioning its prior distribution on the training data. The resulting posterior yields a probabilistic projection regarding the latent log-mortality surface at desired input vector $\mathbf{x}_*$ , given the information in the USMDB. Note that $\mathbf{x}_*$ can refer to in-sample cells (historical smoothing) or out-of-sample cells (future forecasts), both obtained from exactly the same formulas below.

Given a prior distribution $f_l \sim \mathcal{GP}(m_l, C_l)$ and a training set $\mathcal{T} = (\mathbf{x}, \mathbf{y}_l)$ , we calculate the posterior distribution $\mathbf{y}_{l,*}|\mathcal{T} \equiv \mathbf{y}_l(\mathbf{x}_*)|\mathcal{T}$ at predictive cells $\mathbf{x}_*$ . Observe that $(\mathbf{y}_l,\mathbf{y}_{l,*})$ follows the Multivariate Normal distribution (MVN)

(6) \begin{equation} \begin{bmatrix} \mathbf{y}_l \\ \mathbf{y}_{l,*} \end{bmatrix} \sim \mathcal{N}\left ( \begin{bmatrix} m_l(\mathbf{x}) \\ m_{l}(\mathbf{x}_*) \end{bmatrix}, \begin{bmatrix} C_l(\mathbf{x},\mathbf{x}) + \boldsymbol{\Sigma }_l & C_l(\mathbf{x},\mathbf{x}_*) \\ C_l(\mathbf{x}_*, \mathbf{x}) & C_l(\mathbf{x}_*, \mathbf{x}_*) + \boldsymbol{\Sigma }_{l,*} \end{bmatrix}\right ). \end{equation}

Applying MVN conditioning expressions, the Universal Kriging equations (Rasmussen & Williams, Reference Rasmussen and Williams2006, Section 2.7) below provide both the estimated mean-function coefficients $\boldsymbol{\beta }_l = \big ( \beta _{0,l}, \beta _{1,l})^T$ in (3) and the posterior distribution of $\mathbf{y}_{l,*}$ $ p\big (\mathbf{y}_{l,*} | \mathbf{y}_l \big ) \sim \mathcal{N}\big (m_{l, *}(\mathbf{x}_*), C_{*}(\mathbf{x}_*,\mathbf{x}_*) \big )$ with the posterior mean-variance

(7) \begin{align} m_{l,*}(\mathbf{x}_*) \;:\!=\; & \mathbb{E}[ \mathbf{y}_{l,*}| \mathcal{T}] = \boldsymbol{h}(\mathbf{x}_*)^T\boldsymbol{\hat {\beta }}_l + C_l(\mathbf{x}_*, \mathbf{x})(C_l(\mathbf{x}, \mathbf{x}) + \boldsymbol{\Sigma }_l)^{-1}(\mathbf{y}_l-\boldsymbol{H}\boldsymbol{\hat {\beta }}_l); \end{align}

(8) \begin{align} \boldsymbol{\hat {\beta }}_l \;:\!=\; & \big (\boldsymbol{H}^T\boldsymbol{D} \big )^{-1} \boldsymbol{H}^T \left (C_l(\mathbf{x}, \mathbf{x}) + \boldsymbol{\Sigma }_l \right )^{-1}\mathbf{y}_l; \end{align}

(9) \begin{align} C_{*}(\mathbf{x}_*, \mathbf{x}_*) \;:\!=\; & \text{Cov}(\mathbf{y}_{l,*} | \mathcal{T}) = C_l(\mathbf{x}^*, \mathbf{x}^*) + \boldsymbol{\Sigma }_l + \nonumber \\ & \quad +(\boldsymbol{h}(\mathbf{x}_*)^T - C_l(\mathbf{x}_*,\mathbf{x})\boldsymbol{D})^T (\boldsymbol{H}^T\boldsymbol{D})^{-1}\big ( \boldsymbol{h}(\mathbf{x}_*)^T - C_l(\mathbf{x}_*, \mathbf{x})\boldsymbol{D} \big ) \end{align}

where the matrix ${C}_l (\mathbf{x}, \mathbf{x}_*)_{i,j} = C_l(x_i,x_{j,*})$ represents the covariance between inputs in the training set and predictive locations $\mathbf{x}_*$ , $\boldsymbol{H} = \big ( \boldsymbol{h}(x^1), \ldots , \boldsymbol{h}(x^n) \big )$ and $\boldsymbol{D} \;:\!=\; \big ( C_l(\mathbf{x}, \mathbf{x}) + \boldsymbol{\Sigma }_l)^{-1}\boldsymbol{H}$ .

The predictive distribution of mortality rates at different age and time coordinates represented by $\mathbf{x}_*$ yield the estimated mortality surfaces. Note that (7) can be applied both in-sample (for $\mathbf{x}_*$ s that are in the training set) to obtain smoothed historical experience, as well as out-of-sample, to predict mortality into the future using exactly the same algebraic expressions. In parallel, the GP model also outputs confidence intervals around $m_{l,*}(\mathbf{x}_*)$ based on the posterior covariance $C_*(\cdot , \cdot )$ , depicting the confidence of the model in its own projections.

2.2.3 Shared covariance structure

Assume that we have selected a collection of $L \subset \mathcal{S}$ states and that each state-specific mortality surface $f_l$ , $1 \le l \le L$ , follows a GP $f_l \sim \mathcal{GP}(m_l,C_l)$ . We proceed to create a joint model for the vector $\boldsymbol{f}=(f_1,\ldots , f_L)$ through correlating its components. The motivation is that similar states should share alike mortality rates. Therefore, we impose a shared covariance structure which captures the dependencies between mortality rates across states.

To jointly model $L$ outputs, we need to specify the mean and covariance kernel of the joint GP $\boldsymbol{f}$ . More specifically, let ${\vec {x}}^i \equiv (x_a^i,x_t^i, x^i_{1}, \ldots , x^{i}_L)$ where $x^i_l = \mathbb{I}_{\{\text{population = }l\}}$ ; then we take

(10) \begin{equation} \boldsymbol{f} \sim \mathcal{GP}(\boldsymbol{m}, \boldsymbol{\mathcal{C}}) \quad \text{ where } \boldsymbol{f}(\mathbf{x}) = \big (f_1(\vec {\mathbf{x}}), \ldots , f_L(\vec {\mathbf{x}}) \big ), \end{equation}

$\boldsymbol{m} \in \mathbb{R}^{Ln \times 1}$ is the mean vector whose elements represent the mean functions of each population $l \in L$ , $\{ m_l(\vec {x}) \}_{l=1}^L$ , and $\boldsymbol{\mathcal{C}} \in \mathbb{R}^{Ln \times Ln}$ denotes the covariance matrix across the entire system.

Intrinsic Coregionalization Model (ICM). Directly specifying the cross-covariances of each output pair $f_l,f_{l'}$ becomes unwieldy for $L\gt 3$ , so instead we rely on coregionalized kernels (Huynh & Ludkovski, Reference Huynh and Ludkovski2021) which assume that each output $f_l$ , $1 \le l \le L$ is a linear combination of $Q$ independent latent GPs $\boldsymbol{u} = \big (u_1(\mathbf{x}), \ldots , u_Q(\mathbf{x}) \big )$ with shared covariance kernel $C^{(u)}(\cdot ,\cdot )$ . In our case, we use the $C^{(u)}= C^{(APC)}$ APC kernel from (5).

Let $\boldsymbol{a}^*_q = (a_{1,q}, \ldots , a_{L,q})^T$ , $1 \le q \le Q$ , be the vector containing the $q$ -th factor loadings across all populations $L$ . Then $\boldsymbol{f}(\mathbf{x}) = \sum _{q=1}^Q \boldsymbol{a}_q^* u_q(\mathbf{x})$ , or for the $l$ -th population,

(11) \begin{equation} f_l(\mathbf{x}) = a_{l,1}u_1(\mathbf{x}) + \ldots + a_{l,Q}u_{Q}(\mathbf{x}). \end{equation}

The switch from $L$ separate GPs to $Q$ GPs ( $u_1,\ldots , u_Q$ ) is similar to a PCA or singular value decomposition approach and allows to reduce the number of hyperparameters in the cross-population covariance matrix from $\frac {L(L-1)}{2}$ to $Q \times L$ :

\begin{align*} & \boldsymbol{\mathcal{C}}({x}, {x}') = \text{Cov}\big ( \boldsymbol{f}({x}), \boldsymbol{f}({x}') \big ) = \text{Cov} \left ( \sum _{q=1}^Q \boldsymbol{a}_q^* u_q({x}), \sum _{q=1}^Q \boldsymbol{a}_q^*u_q({x}') \right ) \\ & \hspace {2cm} = \left ( \sum _{q=1}^Q \boldsymbol{a}_q^* \boldsymbol{a}_q^{*T} \right ) \otimes \text{Cov}\big (u_q({x}), u_q({x}') \big ) \equiv B \otimes C^{(u)}({x},{x}'). \end{align*}

The $L \times L$ coregionalization matrix $B \;:\!=\; A A^T$ with entries $B_{l,k}= \sum _{q=1}^Q a_{l,q}a_{k,q}$ has rank $Q$ .

Hyperparameters. The modeling task is ultimately to learn the covariance structure, i.e., the mean and kernel functions based on the training data. The overall set of the ICM MOGP hyperparameters is $\Theta = ( (\theta _{j})_{j \in \{a,t,c\} }, (a_{l,q})_{l=1,\ldots , L, q=1,\ldots , Q}, (\beta _{0,l}, \beta _{1,l}, \sigma ^2_l)_{l=1,\ldots ,L})$ . In our results below, the R package kergp Deville et.al (2015) is used to carry out the respective Maximum Likelihood estimation through Kronecker decompositions, see Huynh & Ludkovski (Reference Huynh and Ludkovski2021, Reference Huynh and Ludkovski2024) for more details. Once $\Theta$ is estimated, ICM MOGP inference reduces to evaluating the linear-algebraic formulas in (7)–(9).

3.1 Motivation

To motivate the issue of how to group states, we briefly discuss two GP-based alternatives. First, we recall the base case of constructing a Single Output GP (SOGP) model state-by-state. This can be done using the methodology in Ludkovski, et al. (Reference Ludkovski, Risk and Zail2018) and yields independently fitted GP models. Specifically, we fit 51 SOGP models utilizing the kernel (5) with hyperparameters $\theta _a, \theta _t, \theta _c, \sigma ^2, \eta ^2, \beta _0, \beta _1$ . Second, we construct MOGP models based on a geographic partitioning, namely the nine U.S. Census regions, see Appendix A.1. These regions are 3–9 states in size and are purely geographically aligned. From a statistical perspective, the widely varying group size and the large size of some groups (e.g., the South Atlantic region includes 9 states with a total population of over 50 million) are challenging and slow down ICM MOGP performance.

Figure 1 Raw data (black circles, covering years 1990–2018) and smoothed/projected mortality curves (for years 1990–2020) for two representative states based on three different groupings: (i) a single-state SOGP model; (ii) MOGP-GEO with geographic groupings based on U.S. Census regions in Appendix A.1; (iii) proposed MOGP-PCA.

Fig. 1 displays estimated mortality experience at Age 65 in two states as a function of year for the above two choices, as well as the proposed state-characteristic grouping from Section 3. These smoothed mortality estimates are contrasted with the raw mortality rates shown as black circles. For the latter, we observe that the observation noise is much higher in the right panel of Fig. 1 compared to the left panel. This matches the intuition that the noise in mortality data is roughly proportional to the underlying exposed population. Since Arizona’s population is 7 million compared to about 20 million in New York, we expect triple the respective observation variance. The estimated standard deviation of this noise is the hyperparameter $\sigma _l = 3.42\%$ for N.Y. Females, and $\sigma _l = 5.09\%$ for Ariz. Males.

All three GP forecasts are statistically unbiased and carry out viable non-parametric penalized curve fitting. In particular, the in-sample (for years 1990–2018) estimated mortality rates closely match each other. At the same time, there are non-trivial differences among the models when it comes to out-of-sample prediction, best understood as differences in mortality improvement factors. Fig. A.1 in the Appendix highlights some of the discrepancies. In general, we find that single-state models tend to produce more extreme MI and tend to over-smooth the data, while the MOGP-GEO models occasionally overfit. For example, in Fig. A.1 Vermont Males and Minnesota Females’ geographic groupings lead to a negative MI estimate, while the PCA-based grouping projects a positive MI. The take-away is that better MOGP groupings can help maximize stability and interpretability. Moreover, the ideal group size is 3–6 states, enabling sufficient information borrowing from similar states without making the groups too big. Larger groups are challenging for modeling the MOGP cross-covariance described in Section 2.2.3 and affect computational time which is cubic in $L$ .

3.2 Covariates

As documented in prior studies and confirmed in our results below, there is a strong correlation between economic and demographic variables and observed state-level mortality discrepancies. Hence we use various state characteristics to compute a customized similarity metric that drives our group selection. We work with a diverse set $\mathcal{C}$ of 18 non-mortality state-level covariates, chosen to be a representative collection of (1) economic, (2) demographic, and (3) geographic characteristics. The respective state-level data are obtained from several sources, including (FRED https://fred.stlouisfed.org; US Census Bureau https://www.census.gov, and U.S. Bureau of Economic Analysis (BEA). https://www.bea.gov). These sources do not distinguish between the Male and Female subpopulations, hence all the covariates are shared among the genders. Table 2 lists the covariates we work with; Appendix A.3.1 provides a complete description. Our covariates overlap with and are broadly similar to the 20 used in Chetty et al. (Reference Chetty, Stepner, Abraham, Lin, Scuderi, Turner, Bergeron and Cutler2016, Table 8) (who grouped them into Health Behaviors, Healthcare, Environmental, Labor Market, Social Cohesion, and Other Factors), the six used in Couillard et al. (Reference Couillard, Foote, Gandhi, Meara and Skinner2021) and the 11 covariates in Barbieri (Reference Barbieri2020).

Table 2. The 18 selected state covariates in $\mathcal{C}$ . See Appendix A.3.1 for definitions of each covariate

Rather than manually combining the above covariates, many of which are highly correlated, to identify which states are similar, we opt for a statistical solution of principal component analysis (PCA). Thus, we apply PCA to the $51 \times 18$ matrix $\mathcal{C}$ , to identify the main sources of differences across states. Use of PCA to summarize covariates is also advocated in Barbieri (Reference Barbieri2020). Fig. 2 and Table A.1 report the PCA results, where we focus on the first three PCA components $PC1, PC2, PC3$ that together explain over 66% of the variance in $\mathcal{C}$ and allow us to succinctly summarize the drivers of heterogeneity.

The PC1 component can be interpreted as an economic or wealth factor: the most expensive states (N.Y., Cali.) display the highest PC1 loadings, while states like Miss. and La. have lower PC1 loadings. As confirmation, cf. Table A.1 in Appendix A.3.2, most of the economic covariates from $\mathcal{C}$ are positively correlated with the first PC component. The PC2 component can be interpreted as a climate-related factor, with the Southern-most warmest states having highest PC2 loadings. Lastly, PC3 component loosely corresponds to the Sun Belt states: the Southwest plus Texas, Georgia, and Florida. The population in these states is rapidly growing due to internal migration and immigration.

Table A.1 in the Appendix shows the correlation between each covariate and LE at birth, $e_0$ , as constructed by the CDC Center for Disease Control and Prevention (CDC). Similar statistical association between income (MI), poverty rate (PR), and age-specific mortality is documented in Li & Hyndman (Reference Li and Hyndman2021). Other observed correlations are harder to explain, for example between state Religious Percentage and its LE. This supports our motivation to use statistical, rather than causal, ways to capture similarity. A further alternative, which however could be deemed as circular logic, would be to use past mortality trends to identify which states are similar.

Figure 2 State-wise PCA factor loadings $\mathcal{P}_k(s)$ , $k=1,2,3$ .

3.3 Grouping by covariate similarity

The PCA factor loadings obtained above are used to construct groupings of states for the MOGP. We first define a distance metric $\mathfrak{D}(s_1,s_2)$ between any two states $s_1,s_2 \in \mathcal{S}$ . We then group $s$ with its most similar states according to $\mathfrak{D}(s,\cdot )$ , aiming for geographical closeness to help with interpretability. The groups are constructed in a stepwise agglomerative manner, with 3–6 states per group. Note that since the covariates are the same across genders, all the groupings pertain both for Males and Females.

Distance between States. We use PCA components from previous section to define a distance between states. To this end, we compare the respective factor loadings $\mathcal{P}_{k}(s)$ , generating a Euclidean metric weighted by the eigenvalues associated with each PCA component:

(12) \begin{equation} \mathfrak{D}(s_1,s_2)^2 \;:\!=\; \lambda _1\big (\mathcal{P}_{1}(s_1) - \mathcal{P}_{1}(s_2)\big )^2+ \lambda _2\big (\mathcal{P}_{2}(s_1) - \mathcal{P}_{2}(s_2)\big )^2+ \lambda _3\big (\mathcal{P}_{3}(s_1) - \mathcal{P}_{3}(s_2)\big )^2, \end{equation}

where $\lambda _k$ is the eigenvalue associated with PC component $k= 1,2,3$ , see Table 3. The intuition is to prioritize the first factor that explains the majority of the observed variability in our covariates.

Identifying Similar States. Fix a state $s \in \mathcal{S}$ . We construct the set $\mathcal{N}_{s}$ of nearest neighbors of $s$ as follows:

1. (i) Compute the geographical neighbors of $s$ , $O_{1}(s) = \{ s_* : s \text{ and }s_* \text{ are geographically}$ $\text{contiguous}\}$ . That is, $O_{1} \equiv O_1(s)$ is the collection of states which share a physical boundary with $s$ .

2. (ii) Add to $\mathcal{N}_{s}$ the state $s_{1}$ that minimizes the PCA-based distance to $s$ in $O_1$ , $s_{1} = \arg \min _{s_* \in O_{1}} \mathfrak{D}(s,s_*)$ .

3. (iii) Update the neighborhood definition to include $s_1$ :

   \begin{equation*}O_2 \equiv O_{2}(s\cup s_{1}) = \{ s_* : s_* \text{ is geographically contiguous with either } s \text{ or } s_1 \}.\end{equation*}

4. (iv) Repeat steps (ii) and (iii) with $s_n = \arg \min _{s_* \in O_{n}}\mathfrak{D}(s,s_*)$ for $n=2,\ldots , 10$ .

Table 3. Summary of principal component analysis for state covariates

Figure 3 The ten most similar states $\mathscr{N}_s$ for $s=$ Cali. (left), Idaho (middle), and Mich. (right).

Fig. 3 visualizes the results from running the above algorithm on a few representative states and until $|\mathcal{N}_s| = 10$ . The y-axis shows the distance $\mathfrak{D}(s,s_n)$ between the target state $s \in \mathcal{S}$ and the “neighbor” added in round $n$ . We observe the common "zig-zaggging" pattern, indicating that there are states which are further apart geographically, yet closer with respect to $\mathfrak{D}$ ; see, for example, Michigan and Iowa, or Idaho and North Dakota. Indeed, N.Dak. is the closest to Idaho in terms of $\mathfrak{D}(s,\cdot )$ even though they are not geographically neighboring.

Proposed Grouping $\mathscr{O}_s$ for $s \in \mathcal{S}$ . By construction, the neighborhoods $\mathcal{N}_s = \{ s, s_1, s_2, \ldots \}$ are contiguous; given the above non-monotonicity in Fig. 3, we adjust them to take into account $\mathfrak{D}$ -similarity. We start with $s, s_1, s_2 \in \mathscr{O}_s$ ; that is, the first three states from $\mathcal{N}_s$ are in the group of $s$ . Next, let $s_3' = \arg \min _{s^* \in \mathcal{N}_s \backslash (s_1,s_2)} \mathfrak{D}(s^*,s)$ be the next closest state to $s$ in terms of $\mathfrak{D}(\cdot , s)$ . We include $s_3' \in \mathscr{O}_s$ , if and only if $\mathfrak{D}(s_3',s) \lt \max _{i=1,2} \mathfrak{D}(s_i, s)$ ; so that we create a group of 4 states if $s_3'$ is a better (closer) fit in distance to $s$ than either of the two neighbors $s_1, s_2$ already in the group. The above procedure enforces a high degree of geographic contiguity (at least two states are guaranteed to be contiguous with $s$ ) but also allows to add one more state that is not geographically close but is very similar to $s$ in terms of the PCA loadings. As an example, Iowa is added to Mich.’s group (Fig. 4, right), and N. Dak. is added to Idaho’s group (middle panel), while Cali. stays in a group of 3 (no further state is closer than $s_3=$ Wash. in that case, left panel).

As a final step, we ensure that all groups have sufficient underlying population to yield credible estimation. To this end, we augment additional states (in order of their distance $\mathfrak{D}$ ) until the group $\mathscr{O}_s$ has a total population of at least 5 million. This is particularly relevant for the Mountain West region, where Idaho, Montana, Wyoming, North and South Dakota (all with populations under 1.5 M) tend to group together.

A sample of the resulting 51 groupings $\mathscr{O}_s$ (one for each $s$ ) are shown in Fig. 4; the full list of $\mathscr{O}_s$ s is in Appendix A.4. The colors correspond to $\mathfrak{D}(s_*,s)$ , i.e., how close a given state is to its selected neighbors. We note that closeness in terms of $\mathfrak{D}$ varies; there are many very similar (in terms of PCA factor loadings) states in the Midwest, while California’s neighbors are all much less similar to it.

Figure 4 Selected state groupings $\mathscr{O}_s, \, s =$ Cali. (left panel), Idaho (middle), and Mich. (right).

Fig. A.3 in Appendix A.4 visualizes the size of the constructed groups. Recall that by default $|\mathscr{O}_s| = 3$ ; groups of four arise when one of the two most similar (in terms of $\mathfrak{D}$ ) states is not contiguous with $s$ , but farther away. Most of these cases arise in the Midwest and Mid-Atlantic regions. In addition, in Mountain West we have groups of 5 or 6 due to low state population counts.

The groups $\mathscr{O}_s$ are constructed separately for each state. To get a sense of how similar are the groups for different target states, we define the concept of reciprocity. States $A$ and $B$ are in positive reciprocity (PR) if $B \in \mathscr{O}_A$ and $A \in \mathscr{O}_B$ , i.e., they are mutual members in the respective groups. The states that do not experience any PR are La., Miss., R.I., S.Dak., Texas, and W.Va.. This means that these states are quite different from all their geographical neighbors and are not getting selected for their neighbors’ groups. The triple

\begin{equation*} \big (\text{Arizona, Nevada, Utah}\big ) \end{equation*}

is the resulting $\mathscr{O}_s$ for these three states which are both contiguous geographically and are very similar across covariates, forming a mini-cluster. Similarly, the following pairs of states result in the same constructed groups $\mathscr{O}_s$ :

\begin{equation*} \big (\text{Montana, Wyoming}\big ); \quad \big (\text{Connecticut, New York}\big ). \end{equation*}

4.1 Time structure

Fig. 5 shows the bulk behavior of state mortality rates across years 1990–2020 fixing the Age, namely at age 65. As a baseline, we also compute and show the fitted national-level U.S. mortality rate. This curve (dark blue dashed lines in the plot) is generated by fitting a SOGP model to the aggregated U.S. mortality experience. As expected, it lies in the core of the state curves and is close to the population-weighted average of the state-level mortality projections. See also Fig. A.4 in Appendix A.5 for Age 75 counterpart, and the interactive RShiny widget (Ludkovski & Padilla, Reference Ludkovski and Padilla2023) where users can select any other desired Age. Notice that the mortality trend before 2010 is universally positive, while in the past decade mortality has been either stagnating or deteriorating. We also observe that although there is a strong common trend, mortality evolution has been rather variable across states. There are several outlying curves in Fig. 5, notably D.C. and Hawaii, as well as many “cross-overs” where states change relative ranks over time. This heterogeneity increased during the 2010s compared to 2000s, and there are also more cross-overs in the Male populations compared to Females. Moreover, the spread among states is very substantial: e.g., from just over 0.75% inferred mortality in 2020 for 65-year old Females in N.Dak. and Conn., to 1.45% in Miss. and Okla. – a ratio of 1.88 at the right edge of the right panel.

Figure 5 MOGP-PCA mortality rates for age 65 Males (left panel) and Females (right) and years 1990–2020. U.S. national SOGP-smoothed rate is shown as dashed blue.

Figure 6 MOGP-PCA smoothed mortality rates for age 65 Males (left panel) and Females (right) and years 1990–2021. States are ordered from left to right by their mortality in 2018.

A more in-depth visualization is provided by the heatmaps in Fig. 6 where columns denote states and rows denote years; darker gradient corresponds to lower mortality. The states are sorted in increasing order of their mortality as of 2018. The expected behavior as we move up across a column is a smooth transition from red to blue/purple corresponding to improving mortality throughout 1990–2020. While this occurs for some states (notably D.C. is displaying this pattern, having gone from one of the biggest laggards to being middle-ranked by 2020), for many columns the pattern is much more checkered. Note the colored stripes that indicate different historical mortality paths for states that nowadays have similar rates. Only a few states, such as Utah, Colo., and Hawaii show a consistent positive mortality trend throughout the past three decades.

Relative ranks of states are of great interest. By looking at the first few/last columns of the heatmaps in Fig. 6, we can read off the states with the worst/best mortality. Table 4 summarizes which states are projected to have the highest (worst) and lowest (best) mortality in 2020. Northeast and Pacific states have the overall best mortality, and the Southern states are at the bottom. There is a lot of consistency between the genders as well, with Miss., Ala., and Ky. being among the bottom-5 in all four columns on the left of Table 4. We caution against reading too much into the precise rankings: our fitted models yield predictive standard deviations of latent log-mortality on the order of $\sqrt {C_*(x_*, x_*)} \in [0.04,0.1]$ (depending on size of the state), which corresponds to standard errors of about 0.2%–0.3% on the original scale of mortality rates. For example, Colo.’s Male mortality in 2020 is projected to be 1.304% (ranked fifth) but with a 95% predictive interval of $[1.141\%,1.489\%]$ which would be anywhere within the top-20. There is a complex correlation between the projections of different states which makes relative ranks less volatile, but the upshot is that it is not statistically possible to decide whether a given state is in top-5 or top-10. Nevertheless, these rankings echo Chetty et al. (Reference Chetty, Stepner, Abraham, Lin, Scuderi, Turner, Bergeron and Cutler2016) regarding life expectancy for individuals in the bottom income quartile, who have singled out Tenn., Ark., Okla. as the worst performers and Cali. and Vt. as the top-ranked states.

Table 4. Top-5 and bottom-5 states ranked by MOGP-PCA-projected mortality rate in 2020 at ages 65 and 75. The best and worst states are in the first rows

Figure 7 State rankings in terms of age 65 mortality rate across years 2000, 2010, and 2020.

To better visualize the relative ranks of states across time, Fig. 7 shows a bump chart ranking states by their fitted/projected mortality rates in 2000, 2010, and 2020 at Age 65. We observe that all states experienced a decrease in their mortality rates between 2000 and 2010. However, for the 2010s the picture is largely reversed. In fact for Males, only Maine and New York are projected to improve between 2010 and 2020. For Females, the picture is mixed; 16 states (AL, AR, DC, FL, HI, IN, KY, LA, MO, MS, NM, OH, OK, TN, UT, WV) are projected to have a worse Female Age 65 mortality in 2020 compared to 2010, while the rest are improving. Our analysis corroborates the county-level findings in Shull et.al (2025) who report an improvement in mortality between 2000 and (approximately) 2014, followed by a recent deterioration in county-level rates popularly attributed to deaths of despair.

Moreover, Fig. 7 indicates a certain split among Female mortality: laggard states, primarily in the South, have deteriorating Female mortality in 2020 compared to 2010, while the best-performing states in the Northeast and West continue to experience improving mortality. In other words, there is a divergent pattern where states with low (Female) mortality are doing relatively better compared to states with higher mortality, with this pattern aligning with regional partitions. Similar age-aggregated conclusions appear in Fenelon (Reference Fenelon2013), Chetty et al. (Reference Chetty, Stepner, Abraham, Lin, Scuderi, Turner, Bergeron and Cutler2016).

Figure 8 MOGP-PCA projected mortality rates (on the log scale) across ages 60–85 for year 2020. The dashed blue line shows the fitted SOGP model for the U.S. nationwide mortality.

4.2 Age structure

To complement our discussion of the temporal pattern of mortality, we next discuss its structure in Age. Fig. 8 presents the MOGP-smoothed age structure of mortality for Males and Females. The online RShiny app (Ludkovski & Padilla, Reference Ludkovski and Padilla2023) shows an interactive version of this Figure for any user-specified Year. Given that mortality increases exponentially in Age, we plot log-mortality rates, which are roughly linear, matching the prior mean function in (3). We observe that in bulk the Age structure is very consistent across states, meaning that there is a high correlation between relative mortality experience at different Ages. Consequently, state-level mortality rankings are largely age invariant. For example, Connecticut Male mortality is the lowest in the nation at Age 60 and is 4th, 5th, and 4th lowest at Ages 65, 70, 75, respectively (Females similarly rank 3rd, 3rd, 4th, and 5th at those ages). At the other end of the spectrum, Mississippi has the highest Male mortality across all ages and is in the worst-5 for all ages for Females as well.

Although the overall shape in Fig. 8 is cylindrical, implying a fixed spread between log-mortality of the best and worst state, as a function of Age, the individual behavior of states relative to the U.S. average exhibits highly heterogeneous patterns. The right panel of Fig. 9 (also available interactively) shows those ratios for Males in 2018 for a selection of representative states. For some states, younger Ages are doing better than average, while older Ages are worse off, see e.g., Wisconsin where Male mortality at Age 60 is 84% of national average, while at Age 80 it is 103% of national average. For other states, younger Ages do relatively worse than old Ages; this is the pattern for Tenn. and S.C. in Fig. 9. Male Tennesseans aged 60 have mortality that is 37% above U.S. average, while those aged 80 are only 21% above average. Yet other states have no discernible trend: Delaware Males are within 8% of U.S. average rates across all ages.

In all, we observe distinct clusters of state mortality age structure: the increasing pattern of Wisc. is repeated for AK, ID, IA, ME, NE, NH, RI, UT, VT, VA; the decreasing pattern of S.C. and Tenn. is repeated in AR, DC, KY, LA, MS, NV, OK, and WV. CA, CO, CT, MA, MN, MT, NJ, NY, ND, OR, SD, WA, and WY have a largely Age-invariant negative gap to national average (doing better at all ages), while IN, KS, MO, NC, OH, PA, and TX have a largely Age-invariant positive gap to national average (doing a bit worse at all Ages). Finally, there are a couple of idiosyncratic outliers like Florida, where younger ages are worse than average but older ages are much better than average. The above reveals a novel structural affinity among state mortalities.

We moreover document a widening gap between mortality rates of best/worst states across time. Thus, the spread between states in Fig. 8 has been increasing over the past two decades, something already observed in Fig. 5 at Age 65. The right panel of Fig. 9 shows the ratio of mortality of the second best state vs. the second worst state (we remove the best and the worst to stabilize this metric, akin to winsorizing) at four representative Ages and across years. As expected, the relative spread in mortality is lower for older ages because the underlying rates themselves increase in Age (so the absolute spread is in fact growing in Age). However, the noteworthy pattern is that the dispersion dramatically increased since $\sim$ 2005. In 2020, the Female mortality in the worst states at Age 60 is estimated to be more than double compared to the best states, see the right edge of the right panel in Fig. 9, while thirty years ago it used to be only 50% higher. At Age 70 the ratio is more than 75% in 2020 compared to 45% in 2000.

Figure 9 Left: mortality rates for Males in 2020 expressed as a ratio of U.S. national average, as a function of age for 6 representative states. right: ratio of second-worst to second-best state mortality at four different ages, as a function of year.

Fig. A.5 in Appendix A.6 provides another visualization that addresses two further aspects. First, the figure shows the relative ranks of states at two different Ages, complementing the Year structure in Fig. 7. We observe that several states have drastically different ranks at Age 65 versus Age 75. For example, D.C. Female mortality in 2018 ranks 47th at Age 65 but 10th at Age 75; Ariz. moves from 15th to 7th and Cali. from 8th to 3rd at those ages. In the opposite direction, Vt. drops from 5th to 24th and Maine drops from 20th at Age 65 to 31st for Age 75. Second, Fig. A.5 compares the rankings based on the MOGP setup versus those from the single-population SOGP models. While for most states the rankings are very stable across the models (showcasing that all our results are indeed data-driven), for about a dozen states there are sizable changes. In particular, the SOGP rankings for some of the smaller states are often significantly different. For example, at Age 65, Vt. Females rank 5th best across the MOGP models, but are only 9th based on SOGP smoothing; Mont. ranks 22nd according to MOGP but 12th according to SOGP. At Age 75, SOGP suggests that R.I. is 26th rather than 15th according to MOGP. This illustrates our previous point about credibility and the need to borrow information from “neighbors” to achieve model fit. As discussed in Section 3.1, the SOGP estimates can be unstable when the raw data is very noisy.