An Autoregressive Disease Mapping Model for Spatio-Temporal Forecasting

Abstract

One of the more evident uses of spatio-temporal disease mapping is forecasting the spatial distribution of diseases for the next few years following the end of the period of study. Spatio-temporal models rely on very different modeling tools (polynomial fit, splines, time series, etc.), which could show very different forecasting properties. In this paper, we introduce an enhancement of a previous autoregressive spatio-temporal model with particularly interesting forecasting properties, given its reliance on time series modeling. We include a common spatial component in that model and show how that component improves the previous model in several ways, its predictive capabilities being one of them. In this paper, we introduce and explore the theoretical properties of this model and compare them with those of the original autoregressive model. Moreover, we illustrate the benefits of this new model with the aid of a comprehensive study on 46 different mortality data sets in the Valencian Region (Spain) where the benefits of the new proposed model become evident.

1. Introduction

Statistical techniques for disease mapping have received substantial attention during the last two decades [1,2,3]. These methods seek the smoothing of health indicators by taking into account the geographic structure of the units of study in order to obtain reliable risk estimates, regardless of the statistical size of these units. These methods consider dependence between nearby geographic units, according to their geographical proximity. In principle, that dependence improves the risk estimates of the units of study. Among all the approaches for smoothing risk estimates proposed in the literature, that by Besag et al. (BYM) [4] is particularly popular, and many applications of this model can be easily found.

Although BYM is a benchmark for disease mapping studies, this proposal ignores the temporal evolution of the risks for the units of study, assuming that these are static quantities over time. This hypothesis is not necessarily realistic, especially in those studies where large periods of analysis are considered, since therefore, temporal changes in risks are likely to have occurred somewhere. In addition, if large periods of study are aggregated, we would not retrieve an updated view of the risks at the end of the period, which would be surely our main goal, but its average during all that time period. Thus, any risk excess that could have been estimated for some geographical unit could be just the consequence of a past risk excess that does not exist anymore and on which we should not focus our attention. Additionally, Ocaña and Riola (2007) [5] showed also that the aggregation of data over long time periods in disease mapping studies produces biased risk estimates, and therefore, it would be convenient to consider the temporal variation of risks in those models for correcting that bias. For all these reasons, and evidently for forecasting, spatio-temporal models should be borne in mind in very general contexts, mainly if the available period of study is particularly long. Regretfully, despite the benefits of spatio-temporal analyses, the number of studies of this kind in the applied literature in comparison to regular spatial studies is low. Surely, the reason for this situation is the additional complexity of spatio-temporal models, which have to integrate two sources of dependence, spatial and temporal dependence, into a single proposal. Moreover, spatio-temporal disease mapping models have been mainly used for visualizing the spatio-temporal evolution of risks, and despite the evident interest, forecasting has not focused particular attention on this area.

Spatio-temporal models are useful tools to study the evolution of the geographical distribution of diseases and to predict the risks in future periods of time. Forecasting the risks of diseases could have innumerable applications of practical interest such as planning and locating health resources, the establishment of priorities, and guidance in the development and implementation of public health and social programs to alleviate the geographic inequalities found.

Many proposals can be found in the Bayesian literature for the spatio-temporal modeling of risks. These proposals have followed very different approaches, ranging from spatially dependent temporal parametric models [6,7,8] to more flexible and parameterized spatio-temporal spline models [9,10,11,12,13]. However, the ability of disease mapping models to predict future mortality or disease incidence has been hardly explored [14,15,16]. An interesting proposal among this class of models, of particular importance to this work, was introduced by Martinez-Beneito et al. (2008) [17], and we will refer to this proposal simply as the autoregressive (AR) model. In that paper, the authors proposed a spatio-temporal model for disease mapping that concatenates several BYM spatial models in time by means of first order autoregressive time series. As a consequence, a spatio-temporal structure is defined in that model where the risks for each location and moment are spatially and temporally dependent at the same time. This allows nearby places to have, in general, similar risk estimates and time trends. In addition, Martinez-Beneito et al.’s proposal does not assume any parametric family for the temporal evolution of risks that could restrict them, which seems in principle an advantage of this model. Another interesting advantage of this proposal is that it can be easily implemented in practice in regular Bayesian inference packages such as WinBUGS or INLA, given its spatio-temporal separability [17,18]. Finally, the use of time series to induce temporal dependence in this model could be an important advantage in predictive terms, given the good forecasting capabilities of these tools over other modeling alternatives such as splines, with a much more modest predictive performance.

In this work, we start by introducing Martinez-Beneito et al.’s AR model (Section 2), and afterwards, we describe a motivating application of that model to a large set of mortality causes in the Valencian Region (Section 3), one of the seventeen autonomous regions of Spain. Some limitations of the AR model that could affect its forecasting properties will be evidenced according to the results of the motivating example. For this reason, we will propose (Section 4) an enhancement of the AR model that allows solving the issues found and, therefore, improving the corresponding risk estimates. We also show and compare (Section 5) the results of this new model in the Valencia mortality data set used in Section 3. In addition, we also compare these proposals with two simpler spatio-temporal models proposed in the literature, which consider spatially varying quadratic and linear time trends [6,7], in order to check if less parameterized models could show greater ability for forecasting future out-of-sample risks. Finally, we give some conclusions (Section 6) about the results and models described in the previous sections.

2. An Autoregressive Spatio-Temporal Model for Disease Mapping

In this section, we introduce the autoregressive spatio-temporal disease mapping proposal of [17]. Let $O_{ij}$ and $E_{ij}$ denote, respectively, the number of observed and expected cases for the i-th geographical unit and j-th time interval of the period of study. The data likelihood for this model assumes:

$$O_{ij}\sim Poisson(E_{ij}·exp\left(R_{ij}\right))i=1,\cdots ,I,j=1,\cdots ,J$$

where $exp\left(R_{ij}\right)$ are the corresponding relative risks to be estimated. For the first time interval, the log-relative risks are modeled as in a BYM spatial model, that is as the sum of an intercept and two random effects, one of them inducing spatial dependence and the second inducing heterogeneous variability between spatial units. Specifically, for the first period of study ($j=1$) and $i=1,\dots ,I$, we consider:

$$\begin{array}{cc}\hfill & R_{i1}=\mu +{\alpha }_1+{\left(1-{\rho }^2\right)}^{-1/2}·({\theta }_{i1}+{\varphi }_{i1}),\hfill \\ \hfill & {\theta }_{i1}\sim \mathcal{N}(0,{\sigma }_{\theta }^2),\hfill \\ \hfill & {\mathit{\varphi }}_1=({\varphi }_{11},\dots ,{\varphi }_{I1})\sim ICAR\left({\sigma }_{\varphi }^2\right).\hfill \end{array}$$

(1)

The intercept in the linear predictor above ($\mu +{\alpha }_1$) is the sum of two terms, $\mu$ the overall intercept for all geographical units and time intervals and ${\alpha }_1$ the average deviation from that overall intercept for the first time interval. The second term in the linear predictor above depends on $\rho$, the temporal correlation parameter, which will weight the temporal dependence between time intervals. The contribution of this parameter will be clearer when we introduce the modeling of the subsequent time intervals. Nevertheless, the term ${\left(1-{\rho }^2\right)}^{-1/2}$ in the expression above basically guarantees that the variance of the log risks in the first period is equal to that of the subsequent periods (stationary variance of the autoregressive time series). This term depends on two additional vectors, $\mathit{\theta }$${}_1=({\theta }_{11},\cdots ,{\theta }_{I1})$ and $\mathit{\varphi }$${}_1=({\varphi }_{11},\dots ,{\varphi }_{I1})$, the heterogeneous and spatial random effects, respectively, for the first period. As can be observed, the spatial random effect follows an intrinsic Gaussian conditional autoregressive distribution (ICAR) and the heterogeneous random effect a conditionally independent normal distribution of mean zero and common variance. The combination of the spatial and heterogeneous random effects corresponds to the spatial modeling proposal of BYM. Nevertheless, the autoregressive proposal that we are introducing could be easily modified to include other alternative spatial processes such as Leroux et al.’s distribution [19] or just a proper CAR distribution.

For the subsequent time periods the autoregressive model assumes, for $i=1,\dots ,I$ and $j=2,\dots ,J$:

$$\begin{array}{cc}\hfill & R_{ij}=\mu +{\alpha }_j+\rho ·(R_{i(j-1)}-\mu -{\alpha }_{j-1})+{\theta }_{ij}+{\varphi }_{ij},\hfill \\ \hfill & {\theta }_{ij}\sim \mathcal{N}(0,{\sigma }_{\theta }^2),\hfill \\ \hfill & {\mathit{\varphi }}_j\sim ICAR\left({\sigma }_{\varphi }^2\right),\hfill \\ \hfill & \mathit{\alpha }=({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_J)\sim ICAR\left({\sigma }_{\alpha }^2\right).\hfill \end{array}$$

(2)

Note that the spatial and heterogeneous random effects are temporally and mutually independent for each period in both Expressions (1) and (2). However, the linear predictor in Expression (2) induces temporal dependence for each areal unit, specifically as a first order autoregressive time series. Moreover, the jumps of the time series corresponding to the different areal units are spatially correlated because of the spatial dependence of $\mathit{\varphi }$${}_2,\dots ,$$\mathit{\varphi }$${}_J$. This makes the time series corresponding to the risks of neighboring units follow similar trends. In this way, the relative risk in each geographical unit and time interval depends on the neighboring units, but also on the neighboring time intervals, making information be shared in both time and space.

As for any Bayesian hierarchical model, the autoregressive model needs prior distributions for the rest of the parameters in the model: $\mathit{\alpha }$, $\mu ,\rho ,{\sigma }_{\varphi },{\sigma }_{\theta },{\sigma }_{\alpha }$. In the original proposal, an ICAR distribution in time (this is equivalent to a random walk of first order) is also proposed for modeling $\mathit{\alpha }$, considering as neighbors all the consecutive components of this vector. This induces temporal dependence on the intercepts of all time intervals in the period of study and makes the overall time trend follow a smooth trajectory. The temporal dependence parameter $\rho$ is assumed to follow a uniform prior distribution between −1 and 1, which implicitly entails that the time series for $R_{i·}$ for each i are stationary. Finally, vague prior distributions are assumed for $\mu$, on the whole real line, and the rest of the standard deviations present in the model, on the positive real line.

The AR model just introduced has interesting links with the spatio-temporal modeling framework proposed by Knorr-Held (2000) [20], in particular the four spatio-temporal interaction types defined there. Specifically, the spatio-temporal interaction types (I and III) showing temporal independence can be reproduced by just setting $\rho =0$ in the AR model. Similarly, spatial dependence could be removed from the AR model (interaction Types I and II) by making ${\sigma }_{\varphi }=0$. The Type IV Knorr-Held’s spatio-temporal interaction, showing both spatial and temporal dependence, can be seen as a limit case of the AR model without a heterogeneous term (${\sigma }_{\theta }=0$) and with $\rho$ fixed to one, as this interaction was originally formulated as a random walk instead of an autoregressive process. However, in addition to all four Knorr-Held’s interaction types, which correspond to some particular values of ${\sigma }_{\varphi },{\sigma }_{\theta }$, and $\rho$ in the AR model, this model can reproduce far more spatio-temporal interactions by setting these values to take some more additional values. In this sense, the AR model can be seen as a comprehensive proposal able to reproduce, within a single model, all four Knorr-Held’s interaction types and much more interaction types in between.

The AR model is able to reproduce all four Knorr-Held’s interaction types; this modeling framework cannot be strictly considered a particular case of the AR model. Namely, Knorr-Held’s modeling proposal splits the spatio-temporal variability into three terms: spatial, temporal, and (pure) spatio-temporal terms. Thus, this model has a pure spatial term, common for the whole study period, that the AR model does not have. In part, the AR model ignores this term due to the computational and confounding problems between terms posed by the joint inclusion of those three components. These confounding problems make some of the proposed implementations of this model include substantial constraints in order to avoid them [21,22]. These constraints increase in a significant manner the computational burden and computing time of the model.

3. A Motivating Analysis: A Comprehensive Spatio-Temporal Mortality Study in the Valencian Region

As a motivating example for the rest of the paper, we performed the following analysis. The AR model was originally used as the spatio-temporal smoothing tool for the Spatio-temporal Mortality Atlas of the Valencian Region (1987–2006) [23]. This atlas studies the spatio-temporal evolution for 46 combinations of the cause of death and sex corresponding to 27 causes of mortality, where some particular combinations without enough deaths, such as breast cancer in men, or without biological sense, such as prostate cancer in women, have been removed. Mortality is disaggregated for this study at the municipal level for a total of 540 municipalities and 20 annual periods from 1987 to 2006.

Regarding the results derived from this atlas, those shown in

Table 1. Posterior mean and 95% credible interval for the temporal correlation coefficient ($\rho$) of the AR model for each mortality cause and sex.

	Mortality Causes	Sex
		Men	Women
1.	All tumors	0.981	0.963
		(0.966–0.990)	(0.938–0.983)
2.	Mouth	0.956	–
		(0.898–0.974)
3.	Stomach	0.960	0.890
		(0.912–0.992)	(0.603–0.985)
4.	Colorectal	0.902	0.921
		(0.802–0.962)	(0.794–0.981)
5.	Colon	0.903	0.828
		(0.808–0.971)	(0.499–0.919)
6.	Rectum	0.945	0.888
		(0.848–0.990)	(0.709–0.979)
7.	Liver	0.970	0.822
		(0.932–0.993)	(0.562–0.950)
8.	Vesicle	–	0.960
			(0.912–0.992)
9.	Pancreas	0.957	0.575
		(0.868–0.995)	(−0.239–0.940)
10.	Larynx	0.972	–
		(0.929–0.995)
11.	Lung	0.989	0.947
		(0.977–0.997)	(0.882–0.987)
12.	Breast	–	0.965
			(0.914–0.995)
13.	Uterus	–	0.893
			(0.661–0.988)
14.	Ovary	–	0.874
			(0.422–0.989)
15.	Prostate	0.960	–
		(0.908–0.989)
16.	Bladder	0.974	–
		(0.950–0.991)
17.	Lymphatic	0.934	0.842
		(0.807–0.994)	(0.607–0.961)
18.	Leukemia	0.960	0.108
		(0.912–0.992)	(−0.685–0.895)
19.	Diabetes	0.934	0.961
		(0.887–0.964)	(0.935–0.979)
20.	Hypertensive	0.941	0.949
		(0.867–0.984)	(0.917–0.973)
21.	Ischemic	0.961	0.950
		(0.949–0.972)	(0.875–0.959)
22.	Cerebrovascular	0.953	0.940
		(0.936–0.968)	(0.922–0.955)
23.	Atherosclerosis	0.954	0.942
		(0.936–0.960)	(0.928–0.955)
24.	Other Cardiovascular	0.943	0.936
		(0.919–0.963)	(0.916–0.951)
25.	Pneumonia	0.932	0.892
		(0.889–0.967)	(0.837–0.910)
26.	COPD	0.961	0.951
		(0.938–0.977)	(0.915–0.978)
27.	Cirrhosis	0.980	0.991
		(0.964–0.992)	(0.978–0.999)

are of particular interest for this work. This table shows the temporal correlation coefficient (parameter $\rho$ in the AR model) estimated for each mortality cause and sex.

As can be observed, the temporal correlation parameter estimated is very high for most of the mortality causes considered, frequently greater than 0.9. This agreement is particularly striking, and these common high values for $\rho$ suggest the possible existence of a common geographical pattern for each cause of death that is repeated, for all the time intervals considered, perhaps with slight temporal variations. However, the AR model does not explicitly include a pure spatial component, but just a pure (common) time trend for the whole Valencian Region and the spatio-temporal component already described. The high values observed for the temporal correlation parameter could be just fixing the missing spatial component of the AR model. Furthermore, this model assumes that the autoregressive time series for every municipality has a mean of $\mu$. This implicit assumption could heavily influence the forecasts of this model since, in the long rung, all the predictions will tend to this common value.

4. An Autoregressive Spatio-Temporal Model with a Specific Spatial Component

These results could be suggesting the convenience of including a new pure spatial term in the AR model that could reproduce the common spatial pattern that all time intervals could show. The presence of this new term could make the autoregressive term of the AR model be in charge of reproducing just the spatio-temporal interaction of the matrix O of observed cases instead of the spatial and spatio-temporal variability, as for the results above. In that case, the term in charge of weighting the temporal dependence ($\rho$) for each spatial unit could take values not so close to one if the data really required it. For this new model, the common spatial pattern would have been explicitly modeled, and therefore, it would not have to be reproduced by fitting a high value of $\rho$ in the autoregressive term.

4.1. Modeling Proposal

In this section, we propose a modification of the spatio-temporal autoregressive model. This modification is motivated by the high temporal correlation parameter estimated with this model for all the mortality causes studied in the Spatio-temporal Mortality Atlas of the Valencian Region. As already mentioned, this fact suggests the existence of a common spatial pattern that is repeated during the whole period of study for each cause of death. Therefore, we propose to include a new random effect in the model, which we will denote by ${\delta }_i,i=1,\cdots ,I$, independent of the others effects and common for all time periods. The goal of this new effect will be modeling the hypothetical common pattern that the risk maps could show for all time periods, that is the average risk level of each geographical unit. Therefore, from a forecasting point of view, risk predictions would eventually tend to the corresponding local value instead of to a common average value.

The new spatial random effect will be modeled once again as a BYM process, that is as a sum of a spatial (${\varphi }_i^{\delta }$) ICAR distributed vector and a heterogeneous normal random effect (${\theta }_i^{\delta }$) of mean zero and common variance. More in detail, for the first time period ($j=1$) and for $i=1,\dots ,I$, we define the log relative risks as:

$$\begin{array}{cc}\hfill & R_{i1}=\mu +{\alpha }_1+{\delta }_i+{\left(1-{\rho }^2\right)}^{-1/2}·({\theta }_{i1}+{\varphi }_{i1}),\hfill \\ \hfill & {\delta }_i={\theta }_i^{\delta }+{\varphi }_i^{\delta },\hfill \\ \hfill & {\theta }_i^{\delta }\sim \mathcal{N}(0,{\sigma }_{{\theta }^{\delta }}^2)\hfill \\ \hfill & {\mathit{\varphi }}^{\mathit{\delta }}=({\varphi }_1^{\delta },\dots ,{\varphi }_I^{\delta })\sim ICAR\left({\sigma }_{{\varphi }^{\delta }}^2\right)\hfill \\ \hfill & {\theta }_{i1}\sim \mathcal{N}(0,{\sigma }_{\theta }^2),\hfill \\ \hfill & {\mathit{\varphi }}_{\mathbf{1}}=({\varphi }_{11},\dots ,{\varphi }_{I1})\sim ICAR\left({\sigma }_{\varphi }^2\right)\hfill \end{array}$$

(3)

and, for the following time periods ($j=2,\dots ,J$):

$$\begin{array}{cc}\hfill & R_{ij}=\mu +{\alpha }_j+{\delta }_i+\rho ·(r_{i(j-1)}-\mu -{\alpha }_{j-1}-{\delta }_i)+{\theta }_{ij}+{\varphi }_{ij},\hfill \\ \hfill & {\theta }_{ij}\sim \mathcal{N}(0,{\sigma }_{\theta }^2),\hfill \\ \hfill & {\mathit{\varphi }}_j\sim ICAR\left({\sigma }_{\varphi }^2\right),\hfill \\ \hfill & \mathit{\alpha }=({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_j)\sim ICAR\left({\sigma }_{\alpha }^2\right).\hfill \end{array}$$

(4)

In order to complete the model, we set prior distributions for the rest of parameters as follows:

$$\mu \sim \mathcal{U}(-\infty ,+\infty )$$

$$\rho \sim \mathcal{U}(-1,1)$$

$${\sigma }_{\alpha },{\sigma }_{\theta },{\sigma }_{\varphi },{\sigma }_{{\theta }^{\delta }},{\sigma }_{{\varphi }^{\delta }}\sim \mathcal{U}(0,L)$$

where L is a large enough value intended to make vague the corresponding uniform prior distributions. These distributions follow the same reasoning as those used for the AR model. In particular, vague prior distributions are considered for the mean (log) risk for all geographic units and time periods ($\mu$) and also for the standard deviations of the random effects. The uniform prior distribution for $\rho$ guarantees the stationarity of the log risk time series.

4.2. Dependence Structure

We explore now some theoretical properties of the model just introduced and compare them with the corresponding properties of the AR model. These results will provide useful insight about the dependence structures, and the main differences, of both models.

In order to simplify the length of the expressions, we denote ${\mathit{R}}_j=(R_{1j},\dots ,R_{Ij})$, ${\mathit{\theta }}_j=({\theta }_{1j},\dots ,{\theta }_{Ij})$, ${\mathit{\varphi }}_j=({\varphi }_{1j},\dots ,{\varphi }_{Ij})$ for $j=1,\dots ,J$ and $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_J)$, $\mathit{\delta }={\mathit{\theta }}^{\mathit{\delta }}+{\mathit{\varphi }}^{\mathit{\delta }}$ where ${\mathit{\theta }}^{\mathit{\delta }}=({\theta }_1^{\delta },\dots ,{\theta }_I^{\delta })$ and ${\mathit{\varphi }}^{\mathit{\delta }}=({\varphi }_1^{\delta },\dots ,{\varphi }_I^{\delta })$. In addition, the expression ${\mathit{R}}_j|\dots$, for each j, stands for the distribution of the log relative risk, for the j-th time period, given the model parameters $\mu$, $\mathit{\alpha }$, ${\sigma }_{{\theta }^{\delta }}$, ${\sigma }_{{\varphi }^{\delta }}$, ${\sigma }_{\theta }$, and ${\sigma }_{\varphi }$, i.e., all parameters except the spatial ($\mathit{\delta }$) and spatio-temporal (${\mathit{\theta }}_j$, ${\mathit{\varphi }}_j$, $j=1,\dots ,J$) random effects. In addition, ${\mathbf{\Sigma }}_{\mathit{\delta }}$ denotes the variance-covariance matrix of $\mathit{\delta }|{\sigma }_{{\mathit{\theta }}^{\mathit{\delta }}},{\sigma }_{{\mathit{\varphi }}^{\mathit{\delta }}}$ and ${\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}$ the variance-covariance matrix of ${\mathit{\theta }}_j+{\mathit{\varphi }}_j|{\sigma }_{\theta },{\sigma }_{\varphi }$ for any period j.

First, we deduce the joint distribution of the log risks conditioned to the rest of parameters of the model except the spatial and spatio-temporal random effects. This distribution will be multivariate normal since each risk $R_{ij}$ is a linear combination of normal random variables. In order to determine the joint distribution of all the risks, we need to calculate the corresponding variance-covariance matrix. This is the goal of the following lines.

The variance-covariance matrix for the first time period is:

$$\begin{array}{cc}\hfill var\left({\mathit{R}}_1\right|\dots )& =var(\mathit{\delta }+{\left(1-{\rho }^2\right)}^{-1/2}·({\mathit{\theta }}_1+{\mathit{\varphi }}_1))\hfill \\ \hfill & =var\left(\mathit{\delta }\right)+{\left(1-{\rho }^2\right)}^{-1}var({\mathit{\theta }}_1+{\mathit{\varphi }}_1)\hfill \\ \hfill & ={\mathbf{\Sigma }}_{\mathit{\delta }}+{\left(1-{\rho }^2\right)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }},\hfill \end{array}$$

and for the second time period:

$$\begin{array}{cc}\hfill var\left({\mathit{R}}_2\right|\dots )& =var(\mathit{\delta }+\rho ·({\mathit{R}}_1-\mathit{\delta })+{\mathit{\theta }}_2+{\mathit{\varphi }}_2)\hfill \\ \hfill & =var\left(\mathit{\delta }\right)+{\rho }^{2}var({\mathit{R}}_1-\mathit{\delta })+var({\mathit{\theta }}_2+{\mathit{\varphi }}_2)\hfill \\ \hfill & =var\left(\mathit{\delta }\right)+{\rho }^2{(1-{\rho }^2)}^{-1}var({\mathit{\theta }}_1+{\mathit{\varphi }}_1)+var({\mathit{\theta }}_2+{\mathit{\varphi }}_2)\hfill \\ \hfill & ={\mathbf{\Sigma }}_{\mathit{\delta }}+{\rho }^2{(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}+{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}\hfill \\ \hfill & ={\mathbf{\Sigma }}_{\mathit{\delta }}+{(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}.\hfill \end{array}$$

Iteratively, it can be easily shown that the variance-covariance matrix for the subsequent time periods ($var\left({\mathit{R}}_j\right|\dots )$, $j=3,\dots ,J$) has also as the expression:

$${\mathbf{\Sigma }}_{\mathit{\delta }}+{(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}.$$

Martinez-Beneito et al. showed that the corresponding variance-covariance matrices for ${\mathit{R}}_j,j=1,\dots ,J$ conditioned to the rest of model parameters, for the AR model, have as the expression: ${(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}$. Therefore, including the spatial term $\mathit{\delta }$ in the new model induces an additive term, the spatial variance-covariance matrix ${\mathbf{\Sigma }}_{\mathit{\delta }}$, on the variance-covariance matrix of the (log) risk vector for each time period.

Regarding the covariances of the (log) risk vectors for the different time periods, for any $k>0$, we have:

$$\begin{array}{cc}\hfill cov({\mathit{R}}_j,{\mathit{R}}_{j+k}|\dots )& =cov\left({\mathit{R}}_j,\mathit{\delta }+\rho ({\mathit{R}}_{j+k-1}-\mathit{\delta })+{\mathit{\theta }}_{j+k}+{\mathit{\varphi }}_{j+k}\right)\hfill \\ \hfill & =cov\left({\mathit{R}}_j,\mathit{\delta }+\rho (\rho ({\mathit{R}}_{j+k-2}-\mathit{\delta })+{\mathit{\theta }}_{j+k-1}+{\mathit{\varphi }}_{j+k-1})+{\mathit{\theta }}_{j+k}+{\mathit{\varphi }}_{j+k}\right)\hfill \\ \hfill & =cov\left({\mathit{R}}_j,\mathit{\delta }+{\rho }^2({\mathit{R}}_{j+k-2}-\mathit{\delta })+\sum _{i=0}^1{\rho }^i({\mathit{\theta }}_{j+k-i}+{\mathit{\varphi }}_{j+k-i})\right)\hfill \\ \hfill & =cov\left({\mathit{R}}_j,\mathit{\delta }+{\rho }^k({\mathit{R}}_j-\mathit{\delta })+\sum _{i=0}^{k-1}{\rho }^i({\mathit{\theta }}_{j+k-i}+{\mathit{\varphi }}_{j+k-i})\right)\hfill \\ \hfill & =cov\left({\mathit{R}}_j,(1-{\rho }^k)\mathit{\delta }\right)+{\rho }^{k}var\left({\mathit{R}}_j\right)\hfill \\ \hfill & =(1-{\rho }^k)cov\left({\mathit{R}}_j,\mathit{\delta }\right)+{\rho }^{k}var\left({\mathit{R}}_j\right),\hfill \end{array}$$

where:

$$\begin{array}{cc}\hfill cov({\mathit{R}}_j,\mathit{\delta })& =cov(\mathit{\delta }+\rho ({\mathit{R}}_{j-1}-\mathit{\delta })+{\mathit{\theta }}_j+{\mathit{\varphi }}_j,\mathit{\delta })\hfill \\ \hfill & =cov(\mathit{\delta }+\rho (\rho ({\mathit{R}}_{j-2}-\mathit{\delta })+{\mathit{\theta }}_{j-1}+{\mathit{\varphi }}_{j-1})+{\mathit{\theta }}_j+{\mathit{\varphi }}_j,\mathit{\delta })\hfill \\ \hfill & =cov(\mathit{\delta }+{\rho }^2({\mathit{R}}_{j-2}-\mathit{\delta })+\rho ({\mathit{\theta }}_{j-1}+{\mathit{\varphi }}_{j-1})+{\mathit{\theta }}_j+{\mathit{\varphi }}_j,\mathit{\delta })\hfill \\ \hfill & =cov(\mathit{\delta }+{\rho }^2({\mathit{R}}_{j-2}-\mathit{\delta })+\sum _{i=0}^1{\rho }^i({\mathit{\theta }}_{j-i}+{\mathit{\varphi }}_{j-i}),\mathit{\delta })\hfill \\ \hfill & =cov(\mathit{\delta }+{\rho }^{j-1}({\mathit{R}}_1-\mathit{\delta })+\sum _{i=0}^{j-2}{\rho }^i({\mathit{\theta }}_{j-i}+{\mathit{\varphi }}_{j-i}),\mathit{\delta })\hfill \\ \hfill & =cov(\mathit{\delta }+{\rho }^{j-1}{(1-{\rho }^2)}^{-1/2}({\mathit{\theta }}_1+{\mathit{\varphi }}_1)+\sum _{i=0}^{j-2}{\rho }^i({\mathit{\theta }}_{j-i}+{\mathit{\varphi }}_{j-i}),\mathit{\delta })\hfill \\ \hfill & =var\left(\mathit{\delta }\right)+{\rho }^{j-1}{(1-{\rho }^2)}^{-1/2}cov({\mathit{\theta }}_1+{\mathit{\varphi }}_1,\mathit{\delta })+\sum _{i=0}^{j-2}{\rho }^{i}cov({\mathit{\theta }}_{j-i}+{\mathit{\varphi }}_{j-i},\mathit{\delta })\hfill \\ \hfill & =var(\mathit{\delta },\mathit{\delta })={\mathbf{\Sigma }}_{\mathit{\delta }}.\hfill \end{array}$$

Therefore:

$$\begin{array}{cc}\hfill cov({\mathit{R}}_j,{\mathit{R}}_{j+k}|\dots )& =(1-{\rho }^k){\mathbf{\Sigma }}_{\mathit{\delta }}+{\rho }^k{\mathbf{\Sigma }}_{\mathit{\delta }}+{\rho }^k{(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}\hfill \\ \hfill & ={\mathbf{\Sigma }}_{\mathit{\delta }}+\frac{{\rho }^k}{1-{\rho }^2}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}.\hfill \end{array}$$

As a consequence, in a more compact form, the variance-covariance matrix of ${({\mathit{R}}_1,\dots ,{\mathit{R}}_J)}^{\prime }$ can be expressed as:

$$cov\left({({\mathit{R}}_1,\dots ,{\mathit{R}}_J)}^{\prime }\right)={\mathbf{1}}_J{\mathbf{1}}_J^{{}^{\prime }}\otimes {\mathbf{\Sigma }}_{\mathit{\delta }}+\frac{1}{1-{\rho }^2}{\mathbf{\Sigma }}_{\mathit{\rho }}\otimes {\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }},$$

where ⊗ denotes the Kronecker product of two matrices, and

$${\mathbf{\Sigma }}_{\mathit{\rho }}=\left(\begin{array}{cccc}1& \rho & \dots & {\rho }^{J-1}\\ \rho & 1& \dots & {\rho }^{J-2}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }^{J-1}& {\rho }^{J-2}& \dots & 1\end{array}\right)$$

is the correlation matrix of a first order autoregressive process with correlation coefficient $\rho$.

Note that $cov({\mathit{R}}_j,{\mathit{R}}_{j+k}|\dots )$ for the original AR model is just [17]:

$$\frac{{\rho }^k}{1-{\rho }^2}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }};$$

thus, the spatial variance-covariance matrix ${\mathbf{\Sigma }}_{\mathit{\delta }}$ has the same additive effect on the covariances between different time periods as for the variances for the (log) risk vectors, in comparison to the original AR model.

In addition, the components of any ${\mathit{R}}_j$ for $j=1,\dots ,J$ (given the hyperparameters of the model) follow a normal distribution (of mean $\mu +{\alpha }_j$), since they are linear combinations of normal variables. These results yield the following joint distribution:

$$\left(\begin{array}{c}{\mathit{R}}_1\\ {\mathit{R}}_2\\ ⋮\\ {\mathit{R}}_J\end{array}\right)|\cdots \sim \mathcal{N}\left((\mu ·{\mathbf{1}}_J+\mathit{\alpha })\otimes {\mathbf{1}}_I,{\mathbf{1}}_J{\mathbf{1}}_J^{{}^{\prime }}\otimes {\mathbf{\Sigma }}_{\mathit{\delta }}+\frac{1}{1-{\rho }^2}{\mathbf{\Sigma }}_{\mathit{\rho }}\otimes {\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}\right).$$

As shown in Martinez-Beneito et al. (2008) [17], the covariance of the log risks at different time periods in that model has the form:

$$\begin{array}{c}\hfill cov({\mathit{R}}_j,{\mathit{R}}_{j+k}|\dots )={\rho }^k{(1-{\rho }^2)}^{-1}{\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }},\end{array}$$

for any $k>0$. Therefore, for the AR model, the joint distribution of the (log) risks is:

$$\left(\begin{array}{c}{\mathit{R}}_1\\ {\mathit{R}}_2\\ ⋮\\ {\mathit{R}}_J\end{array}\right)|\cdots \sim \mathcal{N}\left((\mu ·{\mathbf{1}}_J+\mathit{\alpha })\otimes {\mathbf{1}}_I,\frac{1}{1-{\rho }^2}{\mathbf{\Sigma }}_{\rho }\otimes {\mathbf{\Sigma }}_{\mathit{\theta }+\mathit{\varphi }}\right).$$

These results allow us to compare the spatial and temporal dependence structure of both models and how information is shared between the different time periods. Note that the variance-covariance matrix of the new proposal is the sum of the variance-covariance matrix of the AR model and a block matrix with all its elements equal to ${\mathbf{\Sigma }}_{\mathit{\delta }}$. This has important consequences on the correlation between log risks of different time periods. Thus, let us consider $k>0$ and $1\le i\le I$; therefore:

$$\begin{array}{cc}\hfill corr(R_{ij},R_{ij+k}|\dots )& =\frac{cov(R_{ij},R_{ij+k}|\dots )}{\sqrt{var\left(R_{ij}\right|\dots )var\left(R_{ij+k}\right|\dots )}}\hfill \\ \hfill & =\frac{{\left({\Sigma }_{\delta }\right)}_{ii}+\frac{{\rho }^k}{1-{\rho }^2}{\left({\Sigma }_{\theta +\varphi }\right)}_{ii}}{{\left({\Sigma }_{\delta }\right)}_{ii}+\frac{1}{1-{\rho }^2}{\left({\Sigma }_{\theta +\varphi }\right)}_{ii}}\hfill \\ \hfill & \stackrel{k\to \infty }{\to }\frac{{\left({\Sigma }_{\delta }\right)}_{ii}}{{\left({\Sigma }_{\delta }\right)}_{ii}+\frac{1}{1-{\rho }^2}{\left({\Sigma }_{\theta +\varphi }\right)}_{ii}}.\hfill \end{array}$$

Consequently, the correlation between spatial patterns does not vanish to zero for distant time periods, unlike the AR model where ${\left({\Sigma }_{\delta }\right)}_{ii}=0$. This correlation will be larger when the variance of the spatial term is also larger, in comparison to the variance of the spatio-temporal term. Note also that, for the AR model, the convergence rate of the correlations above to zero depends just on $\rho$, and this will be faster for $\rho \approx 0$ and slower for $\left|\rho \right|=1$. Therefore, if there was an underlying persistent spatial pattern for all time periods, which seems reasonable from an epidemiological perspective, the AR model would fit a high value of $\rho$ in order to reproduce that setting. This fact could easily explain the results obtained in the Spatio-temporal Mortality Atlas of the Valencian Region (1987–2006) described in the previous section.

5. A Re-Analysis of the Mortality Study in the Valencian Region

In this section, we return to the analysis of the spatio-temporal trends of the 27 mortality causes considered in the Spatio-temporal Mortality Atlas of the Valencian Region (1987–2006). We now implement the new modeling proposal described in the previous section and compare the results with those of the AR model. Specifically, we compare the temporal correlation coefficients, deviance information criteria (DICs) [24], and log predictive scores [25] on the models considered. The log predictive score allows us to compare the models according to their predictive performance. Furthermore, as suggested by a reviewer, we also considered the spatio-temporal model of Bernardinelli et al. (1995) [6], which considers a linear time trend for each municipality, and an extension of this model to spatially varying quadratic time trends [7]. We will compare these models, less parameterized than those discussed in this paper, in predictive terms for out-of-sample data for the period 2007–2011.

All four models, for each data set, were run in WinBUGS. The syntax of each model and the R code used to fit them can be found, as Annex Material, from the RCode.pdf supplementary document. Three chains were run for each model with 7000 iterations, whose first 2000 iterations were discarded as the burn-in period. Of these, one of every 15 iteration was saved, yielding a final sample size of 1000 iterations. Convergence was assessed by means of the Brooks–Gelman–Rubin statistic [26]; in particular, we required this parameter to be lower than 1.1 for each variable in the models, and the effective sample size [27] was required to be at least 100 for every sampled variable in the models. Convergence was assessed with the aid of the R2WinBUGS package of R.

5.1. Comparison of the Temporal Correlation Coefficients of Both Models

We begin by comparing the temporal correlation coefficients estimated with the AR model and the new model proposed in this paper for each mortality cause and sex.

Table 2. Posterior mean and 95% credible interval for the temporal correlation coefficient ($\rho$) for each model. Model 1 refers to the original AR model in Martinez Beneito et al. (2008) [17], and Model 2 refers to the new model proposed in this paper.

	Mortality Cause	Men		Women
		Model 1	Model 2	Model 1	Model 2
1.	All tumors	0.981	0.872	0.963	0.890
		[0.966, 0.990]	[0.632, 0.979]	[0.938, 0.983]	[0.687, 0.970]
2.	Mouth	0.956	0.178	–	–
		[0.898, 0.991]	[−0.707, 0.920]
3.	Stomach	0.960	0.484	0.890	−0.219
		[0.912, 0.992]	[−0.413, 0.940]	[0.603, 0.985]	[−0.883, 0.857]
4.	Colorectal	0.902	0.138	0.921	0.350
		[0.801, 0.962]	[−0.345, 0.565]	[0.794, 0.981]	[−0.538, 0.901]
5.	Colon	0.903	0.320	0.828	0.234
		[0.808, 0.971]	[−0.184, 0.829]	[0.499, 0.973]	[−0.698, 0.931]
6.	Rectum	0.945	−0.201	0.888	0.207
		[0.848, 0.990]	[−0.857, 0.685]	[0.709, 0.979]	[−0.645, 0.880]
7.	Liver	0.970	−0.077	0.822	0.250
		[0.932, 0.993]	[−0.690, 0.712]	[0.562, 0.950]	[−0.470, 0.861]
8.	Vesicle	–	–	0.887	0.037
				[0.444, 0.990]	[−0.702, 0.875]
9.	Pancreas	0.957	0.047	0.575	0.199
		[0.868, 0.995]	[−0.680, 0.827]	[−0.239, 0.940]	[−0.620, 0.870]
10.	Larynx	0.972	0.412	–	–
		[0.929, 0.995]	[−0.559, 0.961]
11.	Lung	0.988	0.741	0.947	0.516
		[0.977, 0.997]	[0.101, 0.976]	[0.882, 0.987]	[−0.597, 0.955]
12.	Breast	–	–	0.965	0.223
				[0.914, 0.995]	[−0.621, 0.904]
13.	Uterus	–	–	0.893	0.106
				[0.661, 0.989]	[−0.604, 0.819]
14.	Ovary	–	–	0.874	−0.027
				[0.421, 0.989]	[−0.882, 0.854]
15.	Prostate	0.960	0.866	–	–
		[0.908, 0.989]	[0.553, 0.978]
16.	Bladder	0.974	0.640	–	–
		[0.950, 0.991]	[−0.129, 0.976]
17.	Lymphatic	0.934	0.100	0.842	0.461
		[0.807, 0.994]	[−0.777, 0.910]	[0.607, 0.961]	[−0.652, 0.949]
18.	Leukemia	0.596	−0.018	0.108	−0.157
		[−0.531, 0.973]	[−0.773, 0.858]	[−0.685, 0.895]	[−0.781, 0.709]
19.	Diabetes	0.934	0.862	0.961	0.916
		[0.887, 0.964]	[0.699, 0.950]	[0.935, 0.979]	[0.778, 0.972]
20.	Hypertensive	0.941	0.632	0.949	0.905
		[0.867, 0.984]	[−0.419, 0.964]	[0.917, 0.973]	[0.801, 0.958]
21.	Ischemic	0.961	0.935	0.950	0.929
		[0.949, 0.972]	[0.890, 0.965]	[0.933, 0.963]	[0.879, 0.957]
22.	Cerebrovascular	0.953	0.946	0.940	0.922
		[0.936, 0.968]	[0.912, 0.967]	[0.922, 0.955]	[0.887, 0.947]
23.	Atherosclerosis	0.954	0.947	0.942	0.935
		[0.936, 0.967]	[0.916, 0.966]	[0.928, 0.955]	[0.911, 0.953]
24.	Other Cardiovascular	0.943	0.913	0.936	0.926
		[0.919, 0.963]	[0.849, 0.954]	[0.916, 0.951]	[0.895, 0.948]
25.	Pneumonia	0.931	0.834	0.893	0.875
		[0.889, 0.967]	[0.578, 0.950]	[0.837, 0.935]	[0.793, 0.934]
26.	COPD	0.961	0.936	0.951	0.887
		[0.938, 0.977]	[0.861, 0.974]	[0.915, 0.977]	[0.687, 0.968]
27.	Cirrhosis	0.980	0.939	0.991	0.320
		[0.964, 0.992]	[0.682, 0.993]	[0.978, 0.999]	[-0.630, 0.990]

shows the posterior mean and 95% credible interval for the temporal correlation coefficient (parameter $\rho$) estimated for each model.

As can be seen, the posterior mean for the temporal correlation coefficients ($\rho$) for the new model proposed in this work (Model 2) is always lower than that estimated for the AR model (Model 1). In addition, it can be observed how the decrease in the temporal correlation coefficient is much more remarkable for the tumoral mortality causes. On the contrary, for the rest of the diseases, this difference is not so important. The results of Table 2 seem to confirm the starting hypothesis: the AR model fits high values of the temporal correlation coefficient to compensate the lack of a common spatial pattern for all time periods. When such a pattern is explicitly included in the model, the estimate of the temporal correlation coefficient changes evidently. Therefore, when the common spatial pattern is included in the model, the temporal correlation coefficient is used to fit the dependence strength of the time series, which can be very different for the different data sets, instead of being used to reproduce the missing common spatial pattern of the whole spatio-temporal system as for the original AR model.

5.2. Models’ Fit Comparison According to the DICs

In this subsection, we compare the fit of the AR and the new model in terms of the DICs, for each mortality cause and sex.

Table 3. DICs for each model and data set. Model 1 refers to the original AR model and Model 2 to the new model proposed in this paper. In order to facilitate comparisons, the lowest DIC in each case is highlighted in bold.

	Mortality Cause	Men		Women
		Model 1	Model 2	Model 1	Model 2
1.	All tumors	12,128.3	12,126.5	11,495.0	11,482.6
2.	Mouth	4801.6	4800.8	−	−
3.	Stomach	6901.1	6905.0	5798.1	5796.8
4.	Colorectal	7906.1	7898.7	7652.7	7651.2
5.	Colon	7261.5	7260.1	6966.8	6961.9
6.	Rectum	4721.0	4719.7	4166.9	4162.7
7.	Liver	5949.5	5945.6	5053.5	5052.8
8.	Vesicle	−	−	3866.0	3863.4
9.	Pancreas	5251.4	5248.6	4853.4	4848.9
10.	Larynx	4897.8	4899.1	−	−
11.	Lung	9692.4	9694.5	4675.0	4677.8
12.	Breast	−	−	7712.1	7707.7
13.	Uterus	−	−	5008.8	5004.3
14.	Ovary	−	−	4648.4	4642.6
15.	Prostate	8466.1	8469.1	−	−
16.	Bladder	6548.3	6552.0	−	−
17.	Lymphatic	5140.8	5138.1	4636.0	4638.2
18.	Leukemia	4759.5	4753.4	4282.2	4275.9
19.	Diabetes	7293.5	7294.7	8727.0	8730.9
20.	Hypertensive	5224.8	5231.9	6609.3	6615.8
21.	Ischemic	11,480.4	11,477.0	10,871.8	10,882.5
22.	Cerebrovascular	11,611.7	11,608.1	12,012.9	11,986.7
23.	Atherosclerosis	6400.6	6405.0	7623.1	7600.7
24.	Other Cardiovascular	11,457.9	11,459.2	12,101.0	12,097.3
25.	Pneumonia	7425.7	7426.7	7436.5	7414.3
26.	COPD	10,180.8	10,176.2	7388.2	7396.4
27.	Cirrhosis	6946.4	6944.9	5179.2	5182.3

shows the corresponding DICs for each model and data set.

For men, the new model with a common spatial component (Model 2) has a lower DIC in 13 out of the 23 mortality causes considered, i.e., for 57% of the data sets. For women, this model has a lower DIC in 16 out of the 23 mortality causes, i.e., for 70% of the data sets. Therefore, the DIC seems to show a slight superiority of the new proposal over the original AR model. Interestingly, this superiority is more evident for the tumoral causes of death considered, Causes 2–18 of Table 3. Thus, for these causes, the DIC was lower for the new model in eight out of 13 data sets for men and 11 out of 13 data sets in women. We remind that, precisely, the tumoral causes of death were those where $\rho$ showed the sharpest decrease, so DICs are more similar for the causes of death where this parameter is, in general, similar for both models.

5.3. Prediction for the Next Five Years

Finally, in this subsection, we compare the predictive performance of:

• The original AR model of Martinez-Beneito et al.
• The autoregressive model proposed in this paper.
• The spatio-temporal model of Bernardinelli et al. (1995) [6], which would assume spatially correlated linear time trends for the municipalities of the Valencian Region.
• The spatio-temporal model of Assunção et al. (2001) [7], which would assume spatially correlated quadratic time trends for the municipalities of the Valencian Region.

We will assess the predictive performance of these models for the five years of the period 2007–2011 (out-of-sample data not used in model fitting) and each of the data sets studied in the Spatio-temporal Mortality Atlas of the Valencian Region. We will compare the predictive performance of these four models by means of log predictive scores [25]. These are the sum of the log predictive scores for the number of observed deaths, for each municipality, for the next five upcoming years after the end of the period of study (2007–2011).

Table 4. Log predictive scores for each data set and year for the 2007–2011 period. Model 1 refers to the AR model; Model 2 refers to the new model introduced in this paper; Model 3 refers to Bernardinelli et al.’s proposal; and Model 4 refers to the model proposed by Assunção et al. In order to facilitate comparisons, the largest log predictive score in each case is highlighted in bold.

			2007	2008	2009	2010	2011	Total
All tumors	Men	Model 1	−1016.3	−1060.1	−1113.5	−1099.5	−1114.8	−5404.3
		Model 2	−1019.6	−1063.1	−1119.0	−1107.6	−1117.9	−5427.1
		Model 3	−1077.8	−1112.7	−1169.8	−1155.6	−1142.5	−5658.4
		Model 4	−1046.2	−1096.1	−1171.4	−1179.0	−1182.6	−5675.2
	Women	Model 1	−807.2	−896.8	−885.3	−871.8	−927.4	−4452.0
		Model 2	−871.2	−898.9	−885.5	−872.5	−926.4	−4454.5
		Model 3	−890.9	−921.0	−906.0	−882.3	−917.6	−4517.8
		Model 4	−885.1	−918.3	−905.9	−884.6	−920.7	−4514.6
Mouth	Men	Model 1	−212.7	−200.9	−221.4	−232.0	−233.9	−1101.1
		Model 2	−210.6	−199.1	−220.9	−230.8	−232.2	−1093.6
		Model 3	−213.1	−205.8	−222.7	−234.1	−229.0	−1104.7
		Model 4	−209.9	−203.4	−222.8	−236.6	−232.0	−1104.7
Stomach	Men	Model 1	−331.8	−318.9	−336.0	−344.2	−354.8	−1685.7
		Model 2	−330.2	−318.0	−331.6	−339.9	−348.4	−1668.2
		Model 3	−370.1	−359.2	−357.4	−368.3	−367.1	−1822.1
		Model 4	−366.8	−357.8	−357.4	−369.5	−369.2	−1820.7
	Women	Model 1	−256.0	−282.7	−266.6	−282.6	−259.7	−1374.6
		Model 2	−256.3	−281.7	−265.4	−281.6	−257.9	−1342.9
		Model 3	−270.7	−294.6	−279.4	−297.2	−266.8	−1408.8
		Model 4	−271.3	−294.8	−279.1	−296.6	−266.1	−1407.9
Colorectal	Men	Model 1	−484.5	−532.3	−539.5	−539.9	−544.8	−2641.0
		Model 2	−481.3	−528.8	−539.0	−534.7	−535.0	−2618.9
		Model 3	−470.9	−526.0	−537.0	−525.2	−520.8	−2579.8
		Model 4	−472.8	−528.7	−536.9	−523.1	−518.5	−2579.9
	Women	Model 1	−411.1	−470.4	−439.0	−452.1	−452.6	−2225.3
		Model 2	−410.0	−470.1	−437.7	−451.8	−449.9	−2219.4
		Model 3	−408.2	−468.6	−436.3	−448.5	−446.7	−2208.3
		Model 4	−407.2	−468.1	−436.4	−449.1	−447.5	−2208.2
Colon	Men	Model 1	−409.8	−483.8	−461.4	−493.3	−486.3	−2334.7
		Model 2	−407.5	−480.3	−458.8	−485.8	−477.3	−2309.8
		Model 3	−399.2	−482.9	−461.5	−479.8	−468.2	−2291.5
		Model 4	−402.0	−487.6	−461.6	−478.1	−465.2	−2294.5
	Women	Model 1	−371.5	−423.6	−388.0	−405.0	−414.6	−2002.6
		Model 2	−369.9	−422.5	−386.6	−404.0	−413.2	−1996.1
		Model 3	−367.2	−419.7	−384.0	−398.8	−406.7	−1976.5
		Model 4	−366.7	−419.4	−384.0	−399.3	−407.7	−1977.1
Rectum	Men	Model 1	−243.2	−236.6	−294.4	−260.7	−266.2	−1301.2
		Model 2	−242.1	−235.9	−295.2	−258.2	−263.1	−1294.5
		Model 3	−240.3	−234.4	−293.8	−256.5	−261.3	−1286.3
		Model 4	−240.4	−234.2	−294.1	−256.5	−261.4	−1286.6
	Women	Model 1	−179.4	−198.3	−218.7	−214.3	−194.0	−1004.8
		Model 2	−177.8	−198.3	−216.3	−211.8	−191.6	−995.8
		Model 3	−178.4	−199.5	−213.3	−208.7	−189.2	−989.0
		Model 4	−179.1	−199.9	−213.4	−208.8	−188.6	−989.6
Liver	Men	Model 1	−274.0	−280.5	−275.3	−320.8	−315.6	−1466.2
		Model 2	−272.7	−279.6	−274.3	−318.7	−311.0	−1456.3
		Model 3	−279.6	−288.7	−280.8	−315.9	−307.4	−1472.4
		Model 4	−273.9	−285.2	−280.8	−317.7	−311.1	−1468.8
	Women	Model 1	−198.7	−191.9	−221.5	−210.3	−213.3	−1035.8
		Model 2	−197.3	−190.1	−220.2	−206.8	−210.8	−1025.3
		Model 3	−206.3	−197.6	−223.2	−206.2	−209.1	−1042.4
		Model 4	−206.6	−197.6	−223.2	−205.8	−208.6	−1041.9
Vesicle	Women	Model 1	−126.4	−138.9	−145.8	−138.3	−142.4	−691.7
		Model 2	−126.6	−138.0	−146.0	−137.7	−141.8	−690.1
		Model 3	−142.6	−147.7	−153.9	−145.4	−146.0	−735.7
		Model 4	−138.0	−145.8	−153.7	−146.9	−148.4	−732.8
Pancreas	Men	Model 1	−290.2	−301.3	−310.9	−334.2	−351.1	−1587.6
		Model 2	−291.4	−301.5	−309.4	−333.8	−349.0	−1585.3
		Model 3	−291.4	−301.5	−308.5	−332.7	−347.6	−1581.7
		Model 4	−292.5	−302.2	−308.3	−332.4	−346.0	−1581.4
	Women	Model 1	−247.6	−239.00	−268.3	−328.1	−309.7	−1392.7
		Model 2	−247.7	−236.9	−266.6	−324.4	−307.1	−1382.6
		Model 3	−247.3	−236.7	−267.0	−326.6	−307.5	−1385.3
		Model 4	−247.6	−236.8	−267.0	−326.6	−307.2	−1385.3
Larynx	Men	Model 1	−214.8	−220.6	−215.1	−199.8	−211.2	−1061.6
		Model 2	−214.2	−218.9	−214.2	−199.6	−210.4	−1057.4
		Model 3	−233.2	−226.7	−234.7	−214.4	−224.8	−1133.8
		Model 4	−227.1	−224.7	−234.7	−217.3	−230.4	−1134.2
Lung	Men	Model 1	−683.1	−682.9	−723.2	−707.8	−713.2	−3510.1
		Model 2	−684.1	−682.1	−724.4	−709.3	−710.7	−3510.6
		Model 3	−688.1	−679.9	−722.7	−708.5	−706.6	−3505.9
		Model 4	−683.6	−677.7	−723.1	−715.3	−719.6	−3519.2
	Women	Model 1	−307.7	−316.6	−304.3	−308.6	−357.5	−1594.7
		Model 2	−305.3	−313.9	−303.0	−305.6	−354.1	−1582.0
		Model 3	−322.7	−328.5	−310.4	−319.9	−366.8	−1648.2
		Model 4	−319.6	−327.2	−310.8	−322.1	−370.8	−1650.5
Breast	Women	Model 1	−427.1	−428.4	−436.4	−444.6	−470.6	−2207.2
		Model 2	−426.9	−428.0	−434.0	−443.2	−469.1	−2201.3
		Model 3	−438.8	−436.7	−438.9	−450.2	−468.4	−2233.0
		Model 4	−432.9	−434.2	−439.0	−452.7	−471.6	−2230.4
Uterus	Women	Model 1	−218.2	−203.0	−231.6	−225.0	−247.9	−1125.7
		Model 2	−217.7	−202.2	−230.5	−224.9	−245.6	−1120.9
		Model 3	−222.0	−204.9	−231.8	−222.5	−239.6	−1120.8
		Model 4	−222.0	−204.9	−231.6	−222.4	−239.6	−1120.6
Ovary	Women	Model 1	−230.2	−244.1	−223.4	−239.5	−258.0	−1195.1
		Model 2	−229.8	−244.5	−222.9	−239.2	−257.8	−1194.2
		Model 3	−230.2	−242.0	−219.7	−235.8	−254.3	−1181.9
		Model 4	−232.4	−241.9	−219.8	−236.3	−255.3	−1185.7
Prostate	Men	Model 1	−439.7	−429.6	−447.5	−473.9	−497.5	−2288.3
		Model 2	−439.9	−429.0	−447.5	−473.0	−495.8	−2285.3
		Model 3	−457.1	−461.4	−468.2	−486.6	−508.4	−2381.7
		Model 4	−445.9	−452.0	−468.7	−494.2	−523.8	−2384.6

shows the log predictive score for the first 15 causes of death of those considered. The full table with the log predictive scores for all the causes of death is annexed as Supplementary Material to this paper due to space limitations (Table S1). Log predictive scores are reported for each data set of the study and, separately, for each of the years of the period 2007–2011, the five subsequent years to the end of the period of study used originally for the Atlas. In principle, the model with a greater log predictive score is supposed to have a better performance in predictive terms. Moreover, considering separately each of the years of the period 2007–2011 allows us to study how the predictive performance of these models evolves as the distance with the end of the study increases. Unlike other alternative model fit assessment criteria, the log predictive scores used do not reward the overfitting of any model as, in this case, the fit is measured on new data, in contrast to those criteria that use the data used to fit those models.

As can be observed in Table 4, and more in depth in the full table of the Supplementary Material (Table S1), the new model introduced in this paper shows in general larger log predictive scores than the original AR model. Specifically, the log predictive score of the new model for the full period of study is greater that of the AR model for 18 out of the 23 (78%) mortality causes studied in men and for 21 out of the 23 (91%) mortality causes studied in women. This shows the enhanced performance, in predictive terms, of the new proposal. In addition, the predictions of the new model with a common spatial component are even better for the latest years considered, that is, as shown in Table 4, differences for the log predictive scores of Model 1 and 2 are, in general, mild for both models for 2007, but these differences become more evident for the latest predicted years (2010 or 2011). This suggests that the new proposal captures and reproduces the mid-/long-term time trends better than the original AR model. That is, the new model with a spatial term sets a mean value for each municipality where time trends should return in the long run, and this seems to be beneficial in predictive terms. Regarding the full comparison, which includes the rest of the (parametric) models considered, the new model introduced in this paper also shows, in general, higher log predictive scores than Bernardinelli et al.’s and Assunção et al.’s proposals. Thus, we find that the log predictive scores of the new model are the highest for 15 out of the 23 (65%) mortality causes studied in men and for 14 out of the 23 (61%) mortality causes studied in women.

6. Discussion

As described in this paper, the autoregressive spatio-temporal model proposed by Martinez-Beneito et al. (2008) [17] tends to estimate high values of the temporal correlation coefficient ($\rho$) for all the data sets studied in the Spatio-temporal Mortality Atlas of the Valencian Region (1987–2006). As shown, this is an artifact of this model, which is due to the lack of an explicit spatial pattern for the whole period of study. As a consequence, the term in charge of modeling temporal dependence ($\rho$) is forced to take high values in order to reproduce that missing common pattern. One of the main contributions of this work is to disclose this limitation of the autoregressive spatio-temporal model, which questions the time trend fit of this model.

In this paper, we propose a modification of the spatio-temporal original AR model that includes a new random effect, which models the common spatial pattern for the whole period of study. Under a theoretical perspective, this model generalizes Knorr-Held’s spatio-temporal modeling framework [20], allowing reproducing all four models proposed there (and many more models) in a single proposal. In addition, the new proposal introduced makes the temporal correlation coefficients $\rho$ considerably lower for many of the considered data sets. As a consequence, this parameter for the new model is used to tune the strength of the data temporal dependence instead of reproducing the common spatial pattern, explicitly introduced in this model. In addition, we observed that the decrease in $\rho$ is much more remarkable in the tumoral mortality causes; therefore, these causes seem to be benefit more from the inclusion of the common spatial component. This result is of important epidemiological interest, as it points out a significant difference between the temporal correlation of the tumor mortality causes and the other types of diseases, which would be convenient to study in depth in future studies. The study of the differential behavior of the tumoral mortality causes, with respect to the other diseases, is an interesting line of research, which has been already pointed out [28]. The results of this paper seem to confirm that particular performance of tumoral mortality causes. Research on the causes of that particular performance can make an interesting contribution to the spatial epidemiology of tumoral mortality causes.

Regarding the fit of the original AR model and the new proposal introduced in this paper, we conclude that the fit differences during the period of study are mild. This seems to be confirmed by the DICs in Table 3, which point out a better performance of the new model, mainly when the spatial term takes a more evident effect. Nevertheless, these DIC differences are generally mild, which seems to be confirmed by the similar risk estimates of both models. As a consequence, we could conclude that both models are in general flexible enough for modeling the risk time trends for the data sets considered. However, in predictive terms, things are different. As shown, the log predictive score of the new modeling proposal is higher than that of the original AR model for almost all the data sets considered, and therefore, this model has a better behavior when used for forecasting. Moreover, this new model seems to be even better for medium-/long-term forecasting, which makes it particularly interesting in practice. In addition, this new model also shows in general better predictive capabilities than Bernardinelli et al.’s and Assunção et al.’s parametric spatio-temporal models, which suggests that parametric models may be somewhat restrictive for this particular purpose.

As a summary, this new proposal makes a contribution to spatio-temporal disease mapping, of particular use/interest when forecasting is the main goal of the corresponding study, something very common in spatio-temporal modeling. The proposed model would be mainly useful for the analysis of diseases with stable spatial patterns over time, such as those corresponding to chronic diseases, in contrast to more chaotic spatio-temporal mortality patterns, such as COVID-19 mortality. The interest of this approach is its enhancement of the AR model, which already uses the time series methodology, particularly suited for predictive purposes. As shown, this makes this new proposal more suited for predictive purposes than the AR model and much more suited than other spatio-temporal approaches relying on modeling tools, such as polynomial or splines fit, with much more limited predictive skills.