1. Introduction
Life insurers, among others, are interested in modeling and predicting mortality experience for many different groups. There are many ways that the life experiences of individuals vary across a population. Some of these ways are easily quantifiable, such as age, sex, and income. Other things are harder to measure, such as happiness, health, and social connection. Even though it is not possible to measure all of the ways that these lives differ, it is nonetheless expected that life experiences will tend to be more similar in areas that are geographically closer to one another.
Before the 1990s, mortality models were largely deterministic. Dickson et al. (Reference Dickson, Hardy and Waters2020) give a good overview of the development of these methods. Now stochastic modeling is the standard practice for mortality modeling with Lee and Carter (Reference Lee and Carter1992) proposing one of the earliest stochastic models for mortality. Despite their model’s shortcomings, such as lack of smoothness across ages and a lack of spatiotemporal interactions, it is an improvement over purely deterministic methods and has good interpretability.
Many variations on the Lee–Carter model have been proposed; the Cairns-Blake-Dowd model (Cairns et al., Reference Cairns, Blake and Dowd2006) has been widely used for modeling mortality improvements. Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) provide an overview of the differences and a quantitative comparison between the Cairns–Blake–Dowd model and the Lee–Carter model. Melnikov and Romaniuk (Reference Melnikov and Romaniuk2006) and Booth and Tickle (Reference Booth and Tickle2008) discuss these different models as well as older models and how they have been used for mortality forecasting.
One thing that is missing from all of these models is spatial correlation. It seems quite reasonable that mortality would be more similar in areas that are geographically closer as there are similarities in areas that are close together and are not easily measured. Thus, including spatial effects in our modeling could improve our understanding of mortality. Many researchers have incorporated spatial and spatiotemporal effects into their mortality modeling. Clayton and Kaldor (Reference Clayton and Kaldor1987) and Manton et al. (Reference Manton, Woodbury, Stallard, Riggan, Creason and Pellom1989) were some of the first to use Empirical Bayes to account for spatial correlation. A review of Empirical Bayes and fully Bayesian approaches for modeling spatial variation in mortality rates is given by Bernardinelli et al. (Reference Bernardinelli, Clayton, Pascutto, Montomoli, Ghislandi and Songini1995). Waller et al. (Reference Waller, Carlin, Xia and Gelfand1997) used a hierarchical Bayesian approach with spatial, temporal, and spatiotemporal effects to model county-level lung cancer death rates in the state of Ohio. Xia and Carlin (Reference Xia and Carlin1998) studied the same data using a similar technique but also incorporating relevant covariates such as age and smoking prevalence.
More recently, Ayele et al. (Reference Ayele, Zewotir and Mwambi2015) used Gaussian Markov random fields to account for spatial variation in an additive logistic regression model for child mortality rates in Ethiopia. We note that Biffis and Millossovich (Reference Biffis and Millossovich2006) also developed a random field framework to study mortality and its impact on the liabilities of insurers, but used dimensions based on demographic, cohort, and other risk factors, rather than geographic locations. Dwyer-Lindgren et al. (Reference Dwyer-Lindgren, Bertozzi-Villa, Stubbs, Morozoff, Kutz, Huynh, Barber, Shackelford, Mackenbach and van Lenthe2016) used a Bayesian approach to fit a hierarchical model that looked at the relationship between the effects at adjacent counties and utilized county-level covariates. Alexander et al. (Reference Alexander, Zagheni and Barbieri2017) fit a hierarchical model to obtain subnational mortality estimates. Their model, which smoothed across space and time, turned out to be a better fit for mortality data – both simulated U.S. data and real French data divided into 19 age groups – than simpler methods. The improvement was especially noticed in areas with low population. Boing et al. (Reference Boing, Boing, Cordes, Kim and Subramanian2020) studied geographic variation in longevity at three different geographic levels: state, county, and census tract. They used these geographic levels as random effects in a linear regression model and found that the census tract accounts for the largest portion of the geographical variation in longevity, making longevity inequality in the U.S. a more local phenomenon than is often assumed. In a similar vein, Kim and Subramanian (Reference Kim and Subramanian2016) used a model with two geographic levels (county and state) to study mortality differences by location. They argued that including state in models was superior to using counties alone, due to the fact that legislation, policies, and programs that can impact health are often implemented at the state level. Li and Hyndman (Reference Li and Hyndman2021) also considered inequalities in mortality across the United States. They first used Lee–Carter models to produce independent state-level forecasts then used a forecast reconciliation approach to reconcile the state-level and national-level mortality rate forecasts. They projected mortality rates 10 years into the future and found that mortality inequality among states is likely to persist and that mortality improvement rates will slow down in the future. Haberman and Renshaw (Reference Haberman and Renshaw2011) compared several mortality modeling approaches with various versions of ways of treating random and fixed effects.
Recently, Gibbs et al. (Reference Gibbs, Groendyke, Hartman and Richardson2020) fit a spatiotemporal model to county-level mortality data from the contiguous United States using conditional auto-regressive priors and a county-varying linear time trend for each age group. Here we utilize the same data and build on this model in two significant ways. First, we model all age groups together in a multivariate spatiotemporal model. Using a multivariate approach adds significant complexity to the model but improves the model in a similar way that adding spatial or temporal correlation would, by allowing data at one age group to borrow strength from data at neighboring age groups (Royle & Berliner, Reference Royle and Berliner1999; Gelfand, Reference Gelfand2021). We address the added complexity by using strong Markov assumptions on the correlations and estimate the model parameters with integrated nested Laplace approximations (INLA) (Blangiardo & Cameletti, Reference Blangiardo and Cameletti2015). This is one of the first significant applications modeled using multivariate INLA models for spatiotemporal data (Vicente et al., Reference Vicente, Goicoa and Ugarte2020). Incorporating multivariate dependence in spatio-temporal models is rarely done due to the computational difficulty of building models with compounding sources of dependence, so adding this is a significant feature of the model we propose. While using random effects in a mortality model is not new (Biffis, Reference Biffis2005; Loisel & Serant, Reference Loisel and Serant2007), this is the first attempt we are aware of that uses all three axes of dependence as random effects.
The other novel element in our mortality model includes non-linear and spatially varying covariate effects. This allows for an unprecedented degree of flexibility in how variables such as employment, education level, and others affect mortality. Several previous works have demonstrated the relationships between various socioeconomic variables and mortality rates. Fuchs (Reference Fuchs2004) discussed several socioeconomic variables that have been shown to be correlated with health, and in particular, the difficulties that arise as a result of these variables often being interrelated and correlated, and sometimes difficult to measure. Mackenbach et al. (Reference Mackenbach, Stirbu, Roskam, Schaap, Menvielle, Leinsalu and Kunst2008) studied the mortality inequalities in 22 European countries, considering education level and occupation class. They found that mortality and morbidity are inversely correlated with socioeconomic status, but that the magnitude of the inequality varied considerably. Bassanini and Caroli (Reference Bassanini and Caroli2015) focused on the relationship between work and health, in terms of both number of hours worked as well as employment status. They found that an excessive number of hours worked was detrimental to health, but also that a reduction in hours worked due to involuntary job loss also had a deleterious impact. Boing et al. (Reference Boing, Boing, Cordes, Kim and Subramanian2020) incorporated several socioeconomic and demographic variables in their multiple geographic level model and found that education and income level are the most significant of these variables. Villegas and Haberman (Reference Villegas and Haberman2014) is one of many papers that examine mortality disparity by socio-economic group.
Recent papers by Wen et al. (Reference Wen, Cairns and Kleinow2023) and Cairns et al. (Reference Cairns, Kleinow and Wen2024) explore neighborhood effects for modeling mortality as functions of socio-economic factors. Chetty et al. (Reference Chetty, Stepner, Abraham, Lin, Scuderi, Turner, Bergeron and Cutler2016), Currie and Schwandt (Reference Currie and Schwandt2016), and Ezzati et al. (Reference Ezzati, Friedman, Kulkarni and Murray2008) model mortality at the county level in the United States, examining the associations between various socioeconomic factors and life expectancy as we do. They employ different methodologies to analyze mortality patterns and disparities at the county level. Rashid et al. (Reference Rashid, Bennett, Paciorek, Doyle, Pearson-Stuttard, Flaxman and Ezzati2021) is another work that shares many similarities with the present paper, but also some important differences. In particular, they develop a Bayesian hierarchical model that is broadly similar to ours, modeling their mortality rates with a negative binomial distribution and indexing their mortality rates along age, time, and space dimensions. Their data are from England and, as in our work, they use two levels of spatial resolution. They use as covariates various “measures of deprivation,” including poverty, unemployment, and education; they find evidence of increasing inequalities of life expectancy, and that worse mortality is associated with higher levels of deprivation. Their work differs from ours significantly, however, in their use of covariates. We incorporate covariates directly into our statistical model, whereas Rashid et al. (Reference Rashid, Bennett, Paciorek, Doyle, Pearson-Stuttard, Flaxman and Ezzati2021) does not, instead using them to conduct a post-hoc analysis of the relationship between mortality in various locations and the associated values of the measures of deprivation in these areas.
Functionally, we treat the covariates as processes on the covariate space and they are given Gaussian process priors (Shi & Choi, Reference Shi and Choi2011). By allowing these to change across space, we remove the restriction that the predictor variables affect mortality equally across the whole country. It was seen in Gibbs et al. (Reference Gibbs, Groendyke, Hartman and Richardson2020) that this improved the model and helped to clarify the importance of the covariates.
The remainder of this article is organized as follows: Section 2 discusses the data that are analyzed in this study, Section 3 introduces the statistical models, methods, and notation that are used to perform the analysis, Section 4 conveys and discusses the results of this analysis, and Section 5 concludes with discussion of model limitations and potential future work.
2. Data
In this analysis, we use data from Division of Vital Statistics of the National Center for Health Statistics, part of the Centers for Disease Control and Prevention. The data contain information about every death that occurred in the United States from 2000 to 2017 along with demographic information including age, sex, county of residence, county of death, race, and marital status. In this analysis, we only use county of residence (and not county of death), and restrict our attention to those counties of residence belonging to the contiguous United States. Our model relies on connections between nearby regions, both through an overall spatial trend and covariate effects that vary by region. As a result, we decided not to include Alaska and Hawaii because they are geographically disconnected from the contiguous U.S. Additionally, both states have unique administrative structures – Alaska uses boroughs instead of counties, and Hawaii has county-level jurisdictions. While these could be treated as counties, we were concerned about how well they would integrate with the traditional county units used across the rest of the United States.
Using census data with interpolation, we are able to obtain the age group mortality exposure for all of the counties during this time period. The census data are available with ages being placed in 18 buckets, with the first being for those 0-4 years of age and the last being those individuals who are 85 years old or older. The mortality data are divided into the same age buckets so that for each sex, county, year, and age group, we had an exposure and the number of people who died.
As the census data are only gathered every ten years, the majority of the population data we have are estimates. However, we have exact information on the number of deaths. This mismatch results in some anomalies whereby there are some sex, county, year, and age group combinations where the data indicate that there are more people who died than were living in the county at the time. There are also some counties whose populations are quite small; after dividing the counties by sex and into 18 age groups, there are some counties whose population for a given group was zero in a given year. Counties that have an exposure of zero or more deaths than exposures are combined with their neighbor with the largest population. After such modifications, we have 3092 counties that we consider in the model. A full list of counties that were combined can be found in Appendix B. (We use the term “county” to refer to the subdivisions throughout all states in the United States, including Louisiana, where they are actually parishes.)
Figure 1 shows the observed mortality for males age 55–59 and females age 85+ and provides some evidence of spatial correlation. These are shown later in the paper to compare with the corresponding fitted values. Mortality tends to be higher in the south and lower in the Midwest. Similar trends exist across both sexes and other age groups. Figure 9 in Appendix A shows the observed countrywide mortality rate for each age group and sex combination. In general, males have higher mortality than females. Mortality increases with age, with the exception that the youngest group sees a mortality improvement upon reaching age 5. Across time, there appears to be mortality improvement for those individuals 45+.
Six different covariates are used in this analysis: Unemployment, Race, Home Value, Education, Marital Status, and Household Size. These covariates are measured at different frequencies (see Table 1). Table 1 also shows the source of each variable. To clarify a few points regarding the covariates, the covariate unemployment rate denotes the proportion of the labor force actively seeking employment, calculated as the number of unemployed individuals divided by the total labor force. Head of households who declare two or more races are not counted as white head of households. The covariates of unemployment, race, and home value are divided into 20 different groups based on quantiles. The remaining covariates of education, marital status, and household size are divided into groups where each group consists of the states that have the same value for the covariate. Those covariates have respectively 44, 43, and 29 unique values and so that is the number of resulting groups. The arithmetic average of the values in each group is used as the “location” for that group for calculating a Matérn correlation between points. The covariates appear to also have some strong spatial dependence. In Appendix A, Figures 10 and 11 show unemployment and race values for 2010 plotted on the map. There is clear spatial dependence in both of these graphs. The unemployment data are complete except for a few counties in Louisiana for 2005 and 2006. (These correspond to the parts of Louisiana that were most heavily affected by hurricane Katrina.) Similarly, we are missing home value information for a few states in the early 2000s. See Appendix B for a complete list of modifications. These missing data were dealt with by simply letting those affected observations not have an unemployment and home value effect during those times.
Covariates used in the model along with what they are measuring, frequency of measurement, and the source of the data
3. Methods
We employ a multivariate spatiotemporal random effects model to analyze mortality data, a widely accepted approach for handling spatiotemporal variables. While dynamic linear models are another viable option, they primarily excel in predictive analyses. In contrast, our aim is to investigate how various covariates influence mortality trends over time and space. The random effects model offers a more intuitive and flexible framework for this analysis (Banerjee et al., Reference Banerjee, Carlin and Gelfand2003; Cressie & Wikle, Reference Cressie and Wikle2011).
For each sex, we fit distinct models, which are identical in form. Let
$y_{akt}$
be the number of deaths that occurred within age group
$a$
in county
$k$
during year
$t$
. In this application,
$a\in \{1,2,\ldots , 18\}$
,
$k\in \{1,2,\ldots , 3092\}$
, and
$t\in \{1,2,\ldots ,18\}$
. The number of people within age group
$a$
in county
$k$
during year
$t$
is known, and so the binomial likelihood is appropriate for this data:
\begin{align} y_{akt}|\pi _{akt}&\sim \text{Binomial}(n_{akt},\pi _{akt}) \end{align}
The Poisson distribution is another option that has several advantages over the binomial distribution, such as zero modification and overdispersion. Our application, however, does not need either of those features. Because we are modeling at the age-county-year level, there are a few subsets with small populations. Using the Poisson distribution could generate more deaths than people, whereas the binomial distribution guarantees that the number of deaths is at most equal to the number of people. In the larger counties, the two distributions would produce similar results.
Here,
$\pi _{akt}$
is the annual mortality probability for someone who is in age group
$a$
in county
$k$
during year
$t$
, and
$n_{akt}$
is the corresponding population count. We can then relate
$\pi _{akt}$
to the desired effects using the logit link function:
\begin{align} \ln \left (\frac {\pi _{akt}}{1-\pi _{akt}}\right ) &= \beta _0 + \sum _{i=1}^3F_i(\mathbf{x}_{k}) + \sum _{i=1}^3G_{is}(\mathbf{x}_{kt}) + \phi _k + \delta _t + \psi _a + \gamma _{akt} \end{align}
The parameter
$\beta _0$
is an overall intercept for the model. Here,
$ \mathbf{x}_k$
represents the vector of covariates at the county level, and
$\mathbf{x}_{kt}$
represents the vector of covariates at the county level that are also time-dependent. Each
$F_i$
(
$i\in \{1,2,3\}$
) represents a function of the covariate value for the covariates of education, marital status, and household size. A standard linear effect would be
$F_i(\mathbf{x}_k) = \beta _i \mathbf{x}_k$
. Instead, we induce a nonlinear covariate effect by allowing
$F_i$
to be a stochastic process indexed on the covariates
$\mathbf{x}_k$
, specifically a Gaussian process with the Matérn correlation with smoothness
$\kappa = 1.5$
. It is a Gaussian process on the covariate space, meaning the finite-dimensional distribution of the
$n$
-tuple
$(F(x_1), F(x_2), \ldots , F(x_n))$
is a multivariate normal where the covariance is based on the distances between covariate values, i.e.,
$\text{cov}(F(x_1),F(x_2)) = C(|x_2 - x_1|)$
. This is a flexible way to allow the covariates to have nonlinear effects on mortality. Similarly the
$G_{is}$
functions are the nonlinear covariate effects for unemployment, race, and home value. Once again, these are Gaussian processes with the Matérn correlation function; however, instead of simply having one effect for the entire dataset, we allow each state to have its own effect. For the purpose of these effects, we define a state to be one of the 48 contiguous states in the United States (all but Alaska and Hawaii) and Washington D.C. Because we are examining county-level data, it might make more sense to use county-varying effects instead of state-level effects. However, this would not only lead to a nearly unidentifiable model, but it would be an intractable computational challenge. We are nonetheless gaining some advantage over previous models by allowing for some regionalization beyond an overall national effect. We impose a conditional autoregressive prior on the different state effects so that, conditional on all the other states, the effect for a given state only depends upon those states with which it shares a border. For identifiability, all of the individual covariate effects have a sum-to-zero constraint.
The spatial effect is broken down into two components:
\begin{align} \phi _k &= u_k + v_k\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{u}|\tau _u &\sim \mathcal{N}(\mathbf{0},\tau _u^{-1}\mathbf{I})\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{v}|\tau _v &\sim \mathcal{N}\left(\mathbf{0},\tau _v^{-1}\mathbf{W}^{-1}\right) \end{align}
with precision parameters
$\tau _u$
and
$\tau _v$
, where
$\mathbf{u}$
is our iid spatial effect and
$\mathbf{v}$
is our structured spatial effect. For identifiability,
$\mathbf{u}$
and
$\mathbf{v}$
have a sum to zero constraint. Our structured effect follows the type of conditional autoregressive model proposed by Besag et al. (Reference Besag, York and Mollié1991). We say that county
$i$
and county
$j$
are neighbors if they share a border and denote this relationship as
$i\sim j$
. We denote the number of neighbors of county
$i$
by
$n_i$
. Note that in our specification, a county is not its own neighbor.
$\mathbf{W}$
is a county adjacency matrix where, for
$i\neq j$
,
$W_{ij}=1$
if
$i\sim j$
and
$W_{ij}=0$
if
$i\nsim j$
. The diagonal entries are
$W_{ii} = -n_i$
. This creates a conditional independence where conditioned on all other counties, the effect for a given county only depends upon those counties with which it is a neighbor. This structure is more apparent if we write the effect as:
\begin{equation} v_i|v_j \sim \mathcal{N}\left (\frac {1}{n_i}\sum _{j:i\sim j}v_j,\frac {1}{n_i}\tau _v^{-1}\right ) \quad \text{ for } j\neq i \end{equation}
However, both formulations are equivalent.
The temporal effect is also broken down into two components:
\begin{align} \delta _t &= b_t + c_t\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{b}|\tau _b &\sim \mathcal{N}\left(\mathbf{0}, \tau _b^{-1}\mathbf{I}\right)\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{c}|\tau _c &\sim \mathcal{N}\left(\mathbf{0}, \tau _c^{-1}\mathbf{R}_c^{-1}\right) \end{align}
with precision parameters
$\tau _b$
and
$\tau _c$
and structure matrix
$\mathbf{R}_c$
, where
$\mathbf{b}$
is the iid temporal effect and
$\mathbf{c}$
is the structured temporal effect. For identifiability,
$\mathbf{b}$
and
$\mathbf{c}$
have a sum to zero constraint. The structured effect follows a random walk of order 1 so that
$\mathbf{R}_c$
is an
$18\times 18$
tridiagonal matrix of the form:
\begin{align} \mathbf{R}_c &= \begin{bmatrix} -1 & \quad 1 & \quad 0 & \quad \cdots & \quad 0 & \quad 0 & \quad 0 \\[4pt] 1 & \quad -2 & \quad 1 & \quad \cdots & \quad 0 & \quad 0 & \quad 0 \\[4pt] 0 & \quad 1 & \quad -2 & \quad \cdots & \quad 0 & \quad 0 & \quad 0 \\[4pt] \vdots & \quad \vdots & \quad \vdots & \quad \ddots & \quad \vdots & \quad \vdots & \quad \vdots \\[4pt] 0 & \quad 0 & \quad 0 & \quad \cdots & \quad -2 & \quad 1 & \quad 0 \\[4pt] 0 & \quad 0 & \quad 0 & \quad \cdots & \quad 1 & \quad -2 & \quad 1 \\[4pt] 0 & \quad 0 & \quad 0 & \quad \cdots & \quad 0 & \quad 1 & \quad -1 \\ \end{bmatrix} \end{align}
This implies that, given
$c_{t-1}$
and
$c_{t+1}$
,
$c_t$
is conditionally independent of
$c_{t*}$
for all other
$t^*$
not equal to
$t-1,t,$
or
$t+1$
. More flexible prior dependence could be imposed by using an autoregressive prior instead of a random walk, but a random walk is computationally faster in the INLA algorithm and with so much data, the results would be quite similar.
We then have an age group effect again broken down into two components:
\begin{align} \psi _a &= f_a + g_a\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{f}|\tau _f &\sim \mathcal{N}\left(\mathbf{0}, \tau _b^{-1}\mathbf{I}\right)\\[-10pt]\nonumber \end{align}
\begin{align} \mathbf{g}|\tau _g &\sim \mathcal{N}\left(\mathbf{0}, \tau _c^{-1}\mathbf{R}_g^{-1}\right) \end{align}
with precision parameters
$\tau _f$
and
$\tau _g$
and structure matrix
$\mathbf{R}_g$
, where
$\mathbf{f}$
is the iid age group effect and
$\mathbf{g}$
is the structured age group effect. For identifiability,
$\mathbf{f}$
and
$\mathbf{g}$
have a sum to zero constraint. The structured effect follows a random walk of order 1 so we have that
$\mathbf{R}_g$
is an
$18\times 18$
tridiagonal matrix that coincidentally is the same as
$\mathbf{R}_c$
from Equation (10) because the number of age groups and the number of time points happen to be the same. Then, given
$g_{a-1}$
and
$g_{a+1}$
,
$g_a$
is conditionally independent of
$g_{a*}$
for all other
$a^*$
not equal to
$a-1,a,$
or
$a+1$
.
Note that treating age as an axis of dependence in the same way as time and space is a novel feature of our model. Other approaches have either built age-specific models or have added age as a covariate. What our approach allows for is that the model at one age point borrows strength from nearby ages, the same way the spatial random effect allows the model at a specific location to borrow strength from nearby locations.
Finally, we have an iid error term:
\begin{equation*} \gamma _{akt}\stackrel {iid}{\sim }\mathcal{N}\left(0,\tau _\gamma ^{-1}\right) \end{equation*}
where
$\tau _\gamma$
is a precision parameter. All precision parameters (
$\tau$
) are given identical Gamma priors with mean 2000 and standard deviation 2000. This prior implies a prior mean of the standard deviation of the process effects of
$\left (\frac {1}{\sqrt {\tau }}\right )$
, which is approximately 0.022. This value is consistent with the expected range of errors in the process effects. Additionally, large values of
$\tau$
(high precision) help stabilize the model by shrinking the variance of the random effects, which is a form of Bayesian regularization that is particularly useful when dealing with models such as ours that have many random effects. Because our data set is quite large, the exact value of
$\tau$
will not affect posterior calculations significantly.
Rue et al. (Reference Rue, Martino and Chopin2009) proposed a deterministic method for performing Bayesian inference using the INLA, which is implemented in the R-INLA package (Lindgren & Rue, Reference Lindgren and Rue2015). For a given model, the computation time in R-INLA tends to be much faster than traditional MCMC algorithms that have been used for exploring the posterior of a model. As a result, this methodology has become quite popular in recent years. Blangiardo and Cameletti (Reference Blangiardo and Cameletti2015) is an informative textbook on the theory and implementation of basic spatiotemporal inference using the package. Rue et al. (Reference Rue, Riebler, Sørbye, Illian, Simpson and Lindgren2017) and Bakka et al. (Reference Bakka, Rue, Fuglstad, Riebler, Bolin, Illian, Krainski, Simpson and Lindgren2018) provide reviews of INLA and how it works with spatial data. We use this method for our computation as we have a large problem and hence computational efficiency is important.
Several techniques exist to accelerate computation in large dependent models. Some methods leverage dimension reduction, while others capitalize on sparsity in the data structure. Approaches that circumvent dimension reduction, such as INLA, are generally considered to yield better performance (Stein, Reference Stein2014). Other methods that avoid dimension reduction, like Local Approximate Gaussian Process and Nearest Neighbor Gaussian Process, are also expected to perform comparably (Gramacy & Apley, Reference Gramacy and Apley2015; Datta et al., Reference Datta, Banerjee, Finley and Gelfand2016; Heaton et al., Reference Heaton, Datta, Finley, Furrer, Guinness, Guhaniyogi, Gerber, Gramacy, Hammerling and Katzfuss2019). Importantly, we anticipate that the choice of model fitting procedure should not significantly impact the inferential conclusions drawn from the analysis.
Parameters such as precisions of the structured and unstructured effects are estimated using the INLA procedure. Several other modeling decisions, such as the structure of the Gaussian processes for the covariate effects and the structure of the random effects were made by evaluating deviance information criterion (DIC) of the model output for different versions of the model. The most relevant of these DIC comparisons are shown in the Results section. For this paper, we did not try any other model structure besides the random effects model and we did not try other estimation procedures besides INLA. Even with the sparse modeling that INLA allows us, fitting these models was computationally intensive, especially when using the spatially varying non-linear covariate effects. For this application, the additional computational time was not overly onerous, but did constrain the amount of fine tuning that could feasibly be done on the model.
4. Results
4.1. Model comparisons
The model described in Section 3 was fit to the data described in Section 2. We tried three different versions of the model. Specifically, we used one version where we eliminated all the covariate effects (
$F_i$
’s and
$G_{is}$
’s), a second version where we used only countrywide effects and no state-specific covariate effects:
\begin{equation} \ln \left (\frac {\pi _{akt}}{1-\pi _{akt}}\right ) = \beta _0 + \sum _{i=1}^6F_i(\mathbf{x}_{k}) + \phi _k + \delta _t + \psi _a + \gamma _{akt}, \end{equation}
and a third version, which was the full model described in Section 3. Equation (14) is essentially the same as Equation (2) with all spatially varying covariates,
$G_{is}$
replaced by country-wide covariate effects,
$F_i$
. Table 2 shows differences in the DIC values obtained for each model fit on both sets of data, as compared to those for the full model. DIC is a model comparison metric that favors a good fit but also penalizes the complexity of the model (Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and Van Der Linde2002). Since we are interested in seeing if the more complex model with state-specific covariate effects produces a better fit than the simpler models, while accounting for model complexity, we use DIC in our comparison of the three models. As is the case for similar model comparison criteria such as AIC and BIC, lower values of DIC indicate a better model fit, after accounting for model complexity. From Table 2, we see that the full model produces the lowest DIC values for both the male and female models, i.e., the other two models produced positive DIC differences. These differences are quite substantial in the male model, and larger yet in the female model. Thus, we have chosen to use the full model for the remainder of the analysis.
Differences from the full model in the deviance information criterion (DIC) for the three different model versions that were fit to both the male and female data
Comparison of observed and fitted mortality rates for males aged 55–59 and females aged 85+ in 2010.
Figure 1 shows the posterior means of the mortality rates for males aged 55–59 and females aged 85+. The fitted values display similar general patterns to those seen in the observed data. One difference that is noted, however, is that the fitted mortality rates exhibit a greater degree of spatial correlation than the observed ones. This is to be expected, as the fitted model will smooth extra noise in the data.
Posterior mean of the spatial effect (
$\phi _k$
) for the model fit to the female data (left) and male data (right).
4.2. Random effects
Figure 2 shows the posterior mean of the spatial effect (
$\phi _k$
) for the counties in the United States for the female and male models. Recall that these effects have sum-to-zero constraints, meaning they are centered around 0. We see that in both cases, the spatial patterns of the data appear to be very similar. We also see the same pattern manifesting itself that we saw earlier where there tends to be higher mortality in the South and lower mortality in the upper Midwest.
One important aspect of the results to consider is how mortality is changing over time. Figure 3 shows the posterior mean and 95% credible interval for the temporal effects (
$\delta _t$
). From 2000 through roughly 2014, the trend for both males and females is of mortality improvement over time (i.e., a decrease in
$\delta _t$
), though the patterns are not smooth or monotonic. Around 2014, there begins to be a small uptick in the mortality. One potential cause for this increase could be a spike in the “deaths of despair” that are often discussed in the mortality literature (see, for example, Scutchfield and Keck (Reference Scutchfield and Keck2017)).
Posterior means and 95% credible intervals of the temporal effects (
$\delta _t$
). Male values are in blue and female values are in red.
Figure 4 shows the posterior mean and 95% credible interval for the age group effects (
$\psi _t$
). Recall that this effect is built into the model as an additional axis of dependence as opposed to building age-specific models or even adding it as a covariate. For the plots, we assign the 85+ age group a value of 87.5 and use the midpoint for all other age groups. These results are consistent with the results seen in the raw data with respect to the age group. Namely, mortality improves upon aging out of the youngest age group, but then consistently deteriorates as age increases. The pattern is seen to be quite similar between the males and females; the most significant difference between the two patterns is the more pronounced “accident hump” in the males between the ages of 15 and 29.
Posterior mean and 95% credible interval of the age group effects (
$\psi _t$
). Male values are in blue and female values are in red. The credible intervals are hardly visible because they are so tight around the estimates.
4.3. Covariate effects
4.3.1. Non-spatially varying effects
There were three variables used that were the same for every county across the United States. These were education, marital status, and household size. Note that these are still flexible covariate effects, but no advantage was seen in the model results to justify the increased complexity of making these variables spatially varying. Figure 5 displays the posterior means and 95% credible intervals for the covariate effects corresponding to education, marital status, and household size for the female and male models, respectively. In particular, the plots show the effect on log odds of mortality for a given county-level variable, as in Equation (2).
Posterior mean and 95% credible interval for the covariate effects (
$F_i(x_k)$
) for the model fit to both female and male data. The effects displayed correspond to education (left), marital status (center), and household size (right).
Although the entire credible interval contains 0 for each of these covariates, they nonetheless present useful illustrations of the information that can be obtained when the covariates are allowed to have nonlinear effects. As an example, consider the case of the education covariate. If this had been treated as a linear effect, the result would likely have been a slight negative slope. However, just having a negative slope does not communicate as much information as our nonlinear effect. It can be seen from the figure that, rather than having a monotonic linear relationship, there is a trend of slow decrease in mortality as education increases, until the percentage of people with a Bachelor’s degree hits 25%. Then follows a quick decline in mortality followed by a leveling off. This demonstrates the added flexibility and more nuanced conclusions that can be obtained by utilizing nonlinear covariate effects.
4.3.2. Unemployment effect
Figure 6 shows the posterior mean of the random effect for unemployment for three specific states. The vast majority of the states follow one of two trends. Either there is no effect for unemployment, as is seen for Colorado, or there is a downward trend for low unemployment rates that flattens out, as is seen for North Carolina. The negative trend that flattens suggests that low unemployment rates lead to higher mortality.
This effect is not intuitive, and in some ways, defies current understanding of the relationship between unemployment rates and mortality rates. Marmot and Wilkinson (Reference Marmot and Wilkinson2005) review evidence of how unemployment (as well as job insecurity) can impact mental and physical health. They cite studies finding that unemployment is associated with ill health and excess mortality. They note that becoming or remaining unemployed can cause stress and other mental health issues, which can be harmful to health; in addition, they also cite studies that associate unemployment with behavior that harms health. They point out, however, that causality is difficult to establish; one reason is that those individuals with poorer health may struggle to keep or regain employment. In addition, it can be difficult to disentangle any health effects due to unemployment from those associated with poverty. These studies look at individual-level unemployment, which likely has a more direct effect on mortality experience. Since our data is at the county level, the relationship can be harder to understand. Rashid et al. (Reference Rashid, Bennett, Paciorek, Doyle, Pearson-Stuttard, Flaxman and Ezzati2021) find that the highest levels of unemployment are associated with worse life expectancy outcomes, though their study used data from England, thus making direct comparison to our results difficult.
We find that the unemployment effect is very consistent across nearly all states. Further research into the cause(s) of this effect would seem to be warranted, though this inquiry is beyond the scope of this paper. There is in fact only one state, California, where low unemployment corresponds to lower mortality. Again, the fact that the shape of this relationship is so different in California than in the other states could merit further investigation.
Figure 6 shows similar patterns for male and female data, although female mortality tends to be lower for higher unemployment rates. The histograms of county unemployment values are shown for each state. It does not seem to be a concern here, but generally speaking, effects for regions where there is little data should be trusted less. Results for both genders and for all 48 states in the continental U.S., along with the District of Columbia, are given in Supplementary Material Section S1.
Posterior mean of the unemployment effects (
$G_{1s}(x_{kt})$
) for selected states, with male and female effects plotted together.
4.3.3. Race effect
Figure. 7 shows the posterior mean of the random effect for percentage of white households for three specific states as well as the distribution of this covariate for those states. The slopes of the curves are consistently negative and in many cases, such as Alabama, the trend is nearly linear. However, there are states such as Illinois where we see a different pattern. In particular, we see pattern where there is a negative slope in the portion of the graph corresponding to lower percents of white households, but a flattening of the slope in the portion of the graph corresponding to higher percents of white households. These differences support having state-specific effects for race. We do see that some state effects are much more sharply sloped than others. Montana, for example, has a large negative slope, while Alabama is nearly flat, suggesting very little effect.
The strong relationship between race and mortality is consistent with previous results; many studies have sought to explain this relationship. Fiscella and Sanders (Reference Fiscella and Sanders2016) emphasize racial and ethnic disparities in healthcare access and quality, potentially explaining why counties with higher white household percentages exhibit lower mortality rates due to better access to medical services. Kurian and Cardarelli’s systematic review (Reference Kurian and Cardarelli2007) underscores variations in health behaviors across racial and ethnic groups, including smoking and physical activity, and how these variations impact cardiovascular disease risk factors. They thus suggest that cultural differences in lifestyle choices may contribute to the observed mortality trends. Additionally, Gee and Payne-Sturges (Reference Gee and Payne-Sturges2004) discuss the impact of environmental factors on health disparities, implying that counties with higher percentages of white households may enjoy better environmental conditions, further contributing to reduced mortality rates. In a similar vein, Morello-Frosch and Shenassa (Reference Morello-Frosch and Shenassa2006) study the effect of environmental hazards such as pollution exposure on racial disparities in maternal and child health.
Posterior mean of the race effects (
$G_{2s}(x_{kt})$
) for selected states, with male and female effects plotted together as well as the distribution of covariates for those states.
4.3.4. Home value effect
Figure 8 shows the posterior mean of the random effect for home value for three specific states. Home value seems to have a more significant effect on female mortality than male mortality. Consider, for example, the states of Delaware, Maryland, and Nevada. While both the male and female models exhibit an inverted U-shape, the effect is more exaggerated for the females than it is for the males. In particular, between about 100,000 and 250,000, the female model mortalities are consistently higher than the males. This is significant because most often when we see separation in these models between males and females, the females have lower mortality rates. Here is a rare case where we see a scenario where female mortality is consistently higher. In contrast, the female mortality is consistently lower than male mortality when average home value exceeds 250,000. These trends are reasonably consistent among most states; this inverted U-shape is seen in the home value effect covariate plots for nearly all states in the female model, and in many states for the male model. This suggests that mortality is highest in counties where the average home price is between 150,000 and 200,000. Mortality is lower when average home value is under 150,000. There are not many counties with average home values above 300,000 during the period being analyzed, so we will not try to infer too much from the end of the curve, although female mortalities are consistently lower than male mortalities. The effects for the remainder of the states and the District of Columbia (and for the male model) can be found in Supplementary Material Section S3.
Posterior mean of the home value effects (
$G_{3s}(x_{kt})$
) for selected states, with male and female effects plotted together as well as the distribution of covariates for each state.
5. Conclusion
We have fit a multivariate spatiotemporal model to mortality data in the contiguous United States. This model has built on the existing mortality modeling literature in two significant ways. First, we model all age groups together to create a multivariate spatiotemporal model. This allows for the borrowing of information not only across space and time but also across the different age groups of the model.
The other significant contribution is the inclusion of nonlinear and spatially varying nonlinear covariate effects on mortality. These nonlinear and spatially varying covariate effects allow us to examine how factors such as education, unemployment, and race affect mortality and how those effects change over space. By including a nonlinear education effect, we were able to see that while mortality generally improves as education increases, there becomes a point where additional education seems to no longer provide additional mortality improvement. This is broadly consistent with the findings of Boing et al. (Reference Boing, Boing, Cordes, Kim and Subramanian2020), which reported a positive correlation between life expectancy and education; Grzywacz et al. (Reference Grzywacz, Almeida, Neupert and Ettner2004) also found that higher levels of education were associated with lower levels of physical and mental stress, leading to better health outcomes. Allowing the covariate effects to change over space allowed us to observe both how the magnitude of the effect changes across states, as we saw for race and home value, and how the shape of the effect changes across state values, as we saw for unemployment. This flexibility allowed us to observe some otherwise undetectable trends, such as the relationship between unemployment and mortality improvement, which had a positive relationship in some areas, a negative relationship in some areas, and no relationship in other areas.
There are ways that the model could be improved. By using the Kronecker product on the precision matrices of our structured random effects, we could create space-time, space-age, and age-time interactions. This could help to make the model more realistic, as mortality is probably not changing over time and age uniformly at all points in the country. We could also create a model in which the covariates interact with age or time so that rather than looking at how the effect of a covariate on mortality changes across space, we could instead explore how it changes across time or the age of the individuals. Another possibility to consider is to include both sexes in a single model; this would increase the dimension of the model but also allow for additional borrowing of strength. Another potential improvement would be accounting for an effect known as “spatial confounding,” where spatial relationships within predictors might be affecting spatial correlations. By properly accounting for potential spatial confounding, we could look more carefully at the precision parameters of the model and draw inference from them.
Currently, we assume the correlation between all neighboring counties to be approximately the same, which is a property of a spatial structure known as isotropy. However, there are oddly shaped counties, and large counties with disproportionate population distributions within the county. There are also counties that, even though they are miles apart from each other, due to a water mass or other feature, may be several orders of connection apart from each other in the current structure. We are currently working on an anisotropic version of this model, where correlation between counties will depend on more than just borders, but will also depend on distance and population distribution. By doing this, a large county with a disproportionate population favoring one side would be more correlated with the county closest to the population center than a county that borders less populated regions of the county.
In addition to the potential improvements to this model, we are also examining how clustering the mortality curves by county can help to provide additional insights (Madrigal et al., Reference Madrigal, Matthews, Patel, Gaches and Baxter2011). Finally, because the relationships between socioeconomic factors and health outcomes are complex and often exhibit non-linear patterns, we believe that models like ours could be useful not only in the modeling of mortality rates but could also be applied to the study of morbidity. For instance, Case and Deaton (Reference Case and Deaton2015) demonstrate the alarming trend of rising morbidity and mortality rates among middle-aged white non-Hispanic Americans in the 21st century, underscoring the importance of socioeconomic disparities in shaping health trajectories. In the context of housing policies, Keene and Geronimus (Reference Keene and Geronimus2011) emphasize the need to evaluate the population health impact of public housing demolition and displacement, highlighting the potential adverse effects of housing instability on health outcomes. Geographical factors also play a significant role in health disparities, as evidenced by Kershaw and Albrecht (Reference Kershaw and Albrecht2014), who explore the impact of metropolitan-level ethnic residential segregation on body mass index among US Hispanic adults. Furthermore, Galea and Vlahov (Reference Galea and Vlahov2005) provide a comprehensive overview of urban health, emphasizing the importance of addressing social determinants of health, including housing affordability and neighborhood characteristics. These studies collectively underscore the need for comprehensive approaches and complex models to address socioeconomic disparities and promote population health.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/S1748499525000016
Acknowledgements
The authors would like to acknowledge Zoe Gibbs for her help in the data curation process. The authors would also like to acknowledge Haavard Rue for his support of our R-INLA code. We also thank two anonymous reviewers, whose comments greatly improved the quality of this paper.
Data availability statement
The mortality data can be requested from Division of Vital Statistics of the National Center for Health Statistics. The covariate data and code that support the findings of this study are available from the corresponding author upon reasonable request.
Funding statement
This work received no specific grant from any funding agency, commercial, or not-for-profit sectors.
Author contributions
Michael Shull: Methodology, Software, Validation, Formal Analysis, Investigation, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualization. Robert Richardson: Conceptualization, Methodology, Validation, Formal Analysis, Investigation, Resources, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. Chris Groendyke: Methodology, Validation, Formal Analysis, Investigation, Writing - Review & Editing, Visualization. Brian Hartman: Conceptualization, Methodology, Validation, Formal Analysis, Investigation, Resources, Writing - Review & Editing, Supervision, Project administration, Funding acquisition.
Competing interests
The authors declare no conflicts of interest.
Appendix A. Mortality and covariate plots
See Appendix Figures 9, 10, and 11.
The countrywide mortality trends for each age group and sex. The plots on the left are for females and the plots on the right are for males. The top plots are for ages 44 and under, and the bottom plots are for ages 45+.
A map of the contiguous United States showing the unemployment rates in 2010.
A map of the contiguous United States showing the proportion of heads of household in the county which are white.
Appendix B. County adjustments
Table 3 contains the Federal Information Processing Standards (FIPS) codes for those counties where we had to adjust the counties for the purpose of the analysis as well as the reason for the adjustment.
Unemployment data was unavailable in 2005 and 2006 in the state of Louisiana for the 7 counties with the following FIPS codes: 22051, 22071, 22075, 22087, 22089, 22095, and 22103.
We were unable to acquire typical home value data for Montana during 2000–2001 and for North Dakota during 2000–2004.
Table of Federal Information Processing Standards (FIPS) adjustments and justifications
Open access
