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:

(1) \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:

(2) \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:

(3) \begin{align} \phi _k &= u_k + v_k\\[-10pt]\nonumber \end{align}

(4) \begin{align} \mathbf{u}|\tau _u &\sim \mathcal{N}(\mathbf{0},\tau _u^{-1}\mathbf{I})\\[-10pt]\nonumber \end{align}

(5) \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:

(6) \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:

(7) \begin{align} \delta _t &= b_t + c_t\\[-10pt]\nonumber \end{align}

(8) \begin{align} \mathbf{b}|\tau _b &\sim \mathcal{N}\left(\mathbf{0}, \tau _b^{-1}\mathbf{I}\right)\\[-10pt]\nonumber \end{align}

(9) \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:

(10) \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:

(11) \begin{align} \psi _a &= f_a + g_a\\[-10pt]\nonumber \end{align}

(12) \begin{align} \mathbf{f}|\tau _f &\sim \mathcal{N}\left(\mathbf{0}, \tau _b^{-1}\mathbf{I}\right)\\[-10pt]\nonumber \end{align}

(13) \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.