Bayesian regression with spatio-temporal varying coefficients

TL;DR

This paper introduces a Bayesian spatio-temporal regression model with a novel dependent process prior to analyze how climate variables affect disease morbidity across different regions and times in Mexico.

Contribution

It proposes a new Bayesian spatio-temporal dependent process prior and demonstrates its ability to model varying climate effects on disease incidence.

Findings

Climate effects vary across space and time.

The model effectively captures spatio-temporal dependence.

Results quantify changing climate impacts on diseases.

Abstract

To study the impact of climate variables on morbidity of some diseases in Mexico, we propose a spatio-temporal varying coefficients regression model. For that we introduce a new spatio-temporal dependent process prior, in a Bayesian context, with identically distributed normal marginal distributions and joint multivariate normal distribution. We study its properties and characterise the dependence induced. Our results show that the effect of climate variables, on the incidence of specific diseases, is not constant across space and time and our proposed model is able to capture and quantify those changes.

Figures

Tables

Table 1. Table 1: Simulated data: Goodness of fit measures for different models defined by neighbourhood type and temporal dependence q 𝑞 q . Bold numbers correspond to the best fit.

		Scenario 1		Scenario 2
Type	$q$	LPML	DIC	LPML	DIC
$s$	$-$	$-42.77$	84.02	$-47.38$	92.27
$t_q$	1	$-43.82$	83.23	$-47.52$	92.55
$t_q$	2	$-40.43$	78.48	$-48.40$	93.69
$s+t_q$	1	$-43.65$	81.93	$-46.70$	90.58
$s+t_q$	2	$-43.76$	81.60	$-48.35$	95.06
$s\times t_q$	1	$-43.10$	82.17	$-47.13$	91.79
$s\times t_q$	2	$-44.86$	85.54	$-47.81$	93.99

Table 2. Table 2: Mexican states and adjacencies.

ID	State	Adjacent states	ID	State	Adjacent states
1	Aguascalientes	14,32	17	Morelos	9,12,15,21
2	Baja California	3,26	18	Nayarit	10,14,25,32
3	Baja California Sur	2	19	Nuevo Leon	5,24,28,32
4	Campeche	23,27,31	20	Oaxaca	7,12,21,30
5	Coahuila	8,10,19,32	21	Puebla	12,13,15,17,20,29,30
6	Colima	14,16	22	Queretaro	11,13,15,16,24
7	Chiapas	20,27,30	23	Quintana Roo	4,31
8	Chihuahua	5,10,25,26	24	San Luis Potosi	11,13,14,19,22,28,30,32
9	Distrito Federal	15,17	25	Sinaloa	8,10,18,26
10	Durango	5,8,18,25,32	26	Sonora	2,8,25
11	Guanajuato	14,16,22,24,32	27	Tabasco	4,7,30
12	Guerrero	15,16,17,20,21	28	Tamaulipas	19,24,30
13	Hidalgo	15,21,22,24,29,30	29	Tlaxcala	13,15,21
14	Jalisco	1,6,11,16,18,24,32	30	Veracruz	7,13,20,21,24,27,28
15	Mexico	9,12,13,16,17,21,22,29	31	Yucatan	4,23
16	Michoacan	6,11,12,14,15,22	32	Zacatecas	1,5,10,11,14,18,19,24

Table 3. Table 3: Real data: Goodness of fit measures for different models defined by neighbourhood type, temporal dependence q 𝑞 q and association parameter c i , t subscript 𝑐 𝑖 𝑡 c_{i,t} . The median neighbourhood size and competing models fittings are also reported. Gastrointestinal data (gastro) and respiratory data (resp). Bold numbers correspond to the best fit in each neighbourhood definition.

				gastro		resp
Type	$q$	$\text{med}\left|{\partial }_{i,t}\right|$	$c_{i,t}$	LPML	DIC	LPML	DIC
$s$	–	5	10	$-88$	$-85$	$-95$	$17$
$s$	–	5	50	$\mathrm{𝟓𝟒𝟗}$	$-\mathrm{𝟏𝟑𝟑𝟐}$	$\mathrm{𝟒𝟗𝟕}$	$-\mathrm{𝟏𝟏𝟒𝟎}$
$s$	–	5	100	$514$	$-1304$	$442$	$-1102$
$t_q$	3	4	10	$-202$	$257$	$-161$	$184$
$t_q$	6	7	10	$172$	$-432$	$188$	$-461$
$t_q$	3	4	50	$634$	$-1399$	$569$	$-1308$
$t_q$	6	7	50	$744$	$-1716$	$\mathrm{𝟕𝟔𝟗}$	$-\mathrm{𝟏𝟕𝟕𝟎}$
$t_q$	3	4	100	$\mathrm{𝟕𝟕𝟔}$	$-\mathrm{𝟏𝟕𝟒𝟖}$	$678$	$-1627$
$t_q$	6	7	100	$673$	$-1668$	$689$	$-1715$
$s+t_q$	3	8	10	$208$	$-552$	$170$	$-492$
$s+t_q$	6	11	10	$370$	$-883$	$394$	$-954$
$s+t_q$	3	8	50	$810$	$-1867$	$738$	$-1703$
$s+t_q$	6	11	50	$\mathrm{𝟗𝟕𝟒}$	$-\mathrm{𝟐𝟏𝟗𝟗}$	$\mathrm{𝟗𝟖𝟐}$	$-\mathrm{𝟐𝟏𝟒𝟖}$
$s+t_q$	3	8	100	$740$	$-1757$	$690$	$-1587$
$s+t_q$	6	11	100	$817$	$-1955$	$935$	$-2061$
$s\times t_q$	3	20	10	$662$	$-1490$	$581$	$-1342$
$s\times t_q$	6	35	10	$\mathrm{𝟕𝟑𝟒}$	$-\mathrm{𝟏𝟕𝟒𝟎}$	$\mathrm{𝟕𝟑𝟒}$	$-\mathrm{𝟏𝟔𝟔𝟓}$
$s\times t_q$	3	20	50	$504$	$-1304$	$379$	$-1043$
$s\times t_q$	6	35	50	$224$	$-631$	$285$	$-752$
$s\times t_q$	3	20	100	$205$	$-595$	$178$	$-551$
$s\times t_q$	6	35	100	$1$	$-88$	$66$	$-206$
CAR, $\rho =0.99$, ${\tau }_c=0.1$				$-2652$	$6042$	$-2620$	$5967$
AR, $\phi =0.99$, ${\tau }_{ar}=0.1$				$-799$	$1627$	$-806$	$1803$
MCAR, $\rho =\phi =0.99$, ${\tau }_{ar}=0.1$				$-97$	$127$	$-96$	$102$
Model 0				$211$	$-417$	$513$	$-1021$
Model 00				$-408$	$819$	$-485$	$970$

Equations

$${\beta }_{i,t}\mid \mathbf{\gamma }$$

$${\gamma }_{i,t}\mid \omega$$

$$\omega$$

$$\text{mbox}Corr({\beta }_{i,t},{\beta }_{j,s})=\frac{c_0\left({\sum }_{(k,r)\in {\partial }_{i,t}\cap {\partial }_{j,s}}{c}_{k,r}\right)+\left({\sum }_{(k,r)\in {\partial }_{i,t}}{c}_{k,r}\right)\left({\sum }_{(k,r)\in {\partial }_{j,s}}{c}_{k,r}\right)}{\left(c_0+{\sum }_{(k,r)\in {\partial }_{i,t}}{c}_{k,r}\right)\left(c_0+{\sum }_{(k,r)\in {\partial }_{j,s}}{c}_{k,r}\right)},$$

$$\boxed{{\beta }_1\mid {\beta }_2\mid {\beta }_3}.$$

$${\eta }_{i,t}=\alpha +{\mathbf{\beta }}_{i,t}^{\prime }{\mathbf{x}}_{i,t}.$$

$$f(\alpha ,\mathbf{\beta },\mathbf{\gamma },\omega )=\text{mbox}N(\alpha \mid m_{\alpha },c_{\alpha })\left\{\prod _k\text{mbox}N\left({\omega }_k\mid m_0,c_0\right)\right\}$$

$$\times \left\{\prod _i\prod _t\prod _k\text{mbox}N\left({\beta }_{i,t,k}|\frac{c_0{m}_0+{\sum }_{(j,s)\in {\partial }_{i,t}}{c}_{j,s}{\gamma }_{j,s,k}}{c_0+{\sum }_{(j,s)\in {\partial }_{i,t}}{c}_{j,s}},c_0+\sum _{(j,s)\in {\partial }_{i,t}}{c}_{j,s}\right)\text{mbox}N\left({\gamma }_{i,t,k}\mid {\omega }_k,c_{i,t}\right)\right\}.$$

$$Y_{i,t}\mid {\mathbf{x}}_{i,t}\sim \text{mbox}N(\mu ({\mathbf{x}}_{i,t}),{\tau }_{i,t})$$

$${\mu }_{\alpha }=\frac{c_{\alpha }{m}_{\alpha }+{\sum }_i{\sum }_t{\tau }_{i,t}\left(y_{i,t}-{\sum }_k{\beta }_{i,t,k}{x}_{i,t,k}\right)}{c_{\alpha }+{\sum }_i{\sum }_t{\tau }_{i,t}}\text{mbox}and{\tau }_{\alpha }=c_{\alpha }+\sum _i\sum _t{\tau }_{i,t}$$

$${\mu }_{\beta }=\frac{c_0{m}_0+{\sum }_{(j,s)\in {\partial }_{i,t}}{c}_{j,s}{\gamma }_{j,s,k}+{\tau }_{i,t}{x}_{i,t,k}(y_{i,t}-\alpha -{\sum }_{l\ne k}{\beta }_{i,t,l}{x}_{i,t,l})}{c_0+{\sum }_{(j,s)\in {\partial }_{i,t}}{c}_{j,s}+{\tau }_{i,t}{x}_{i,t,k}^2}$$

$${\tau }_{\beta }=c_0+\sum _{(j,s)\in {\partial }_{i,t}}{c}_{j,s}+{\tau }_{i,t}{x}_{i,t,k}^2$$

$${\mu }_{\gamma }=\frac{{\omega }_k+{\sum }_{(j,s)\in {\varrho }_{i,t}}\left\{{\beta }_{j,s,k}-\frac{1}{D_{j,s}}\left(c_0{m}_0+{\sum }_{(l,u)\ne (i,t)\in {\partial }_{j,s}}{c}_{l,u}{\gamma }_{l,u,k}\right)\right\}}{1+c_{i,t}{\sum }_{(j,s)\in {\varrho }_{i,t}}\frac{1}{D_{j,s}}}$$

$${\tau }_{\gamma }=c_{i,t}\left(1+c_{i,t}\sum _{(j,s)\in {\varrho }_{i,t}}\frac{1}{D_{j,s}}\right),\text{mbox}with{D}_{j,s}=c_0+\sum _{(l,u)\in {\partial }_{j,s}}{c}_{l,u},$$

$${\mu }_{\omega }=\frac{c_0{m}_0+{\sum }_i{\sum }_t{c}_{i,t}{\gamma }_{i,t,k}}{c_0+{\sum }_i{\sum }_t{c}_{i,t}}\text{mbox}and{\tau }_{\omega }=c_0+\sum _i\sum _t{c}_{i,t}$$

$$a_{\tau }=a_0+\frac{T}{2}\text{mbox}and{b}_{\tau }=b_0+\frac{1}{2}\sum _t{\left(y_{i,t}-\alpha -\sum _k{\beta }_{i,t,k}{x}_{i,t,k}\right)}^2$$

$${\widehat{\text{mbox}CPO}}_{i,t}={\left(\frac{1}{R}\sum _{r=1}^R\frac{1}{\text{mbox}N\left(y_{it}\mid {\alpha }^{(r)}+{\mathbf{\beta }}_{i,t}^{(r{)}^{\prime }}{\mathbf{x}}_{i,t},{\tau }_i^{(r)}\right)}\right)}^{-1},$$

$${\beta }_i\mid \mathbf{\gamma }$$

$${\gamma }_j^i\mid \omega$$

$$\omega$$

$$\text{mbox}Corr({\beta }_i,{\beta }_j)=\frac{c_0{c}_j^{i}I(i\sim j)+\left({\sum }_{{\partial }_i}{c}_k^i\right)\left({\sum }_{{\partial }_j}{c}_k^j\right)}{\left(c_0+{\sum }_{{\partial }_i}{c}_k^i\right)\left(c_0+{\sum }_{{\partial }_j}{c}_k^j\right)},$$

Full text

Keywords: Autoregressive processes, climate analysis, disease mapping, latent variables, stationary processes.

1 Introduction

The Mexican National Institute for Ecology and Climate Change (INECC) wants to study the impact of climate variables in the health sector in Mexico. In particular, they want to characterise the effect of pluvial precipitation and temperature, on the morbidity of gastrointestinal and respiratory diseases. The study is for each of the 32 states of Mexico and for each month in a window of five years. Previous studies that relate climate variables to diseases are those of Morral-Puigmala, (2018) and D’Amato et al., (2014) for gastrointestinal and respiratory diseases, respectively.

Such a study lies with in the scope of disease mapping methods (Lawson,, 2009). These methods aim to describe and quantify spatial, and sometimes temporal, variation in disease risk, including the identification of patterns of possible association between neighbours in a certain region. As a result of a disease mapping study, it is possible to determine areas, and times, of high and low risk in order to contribute to the disease aetiology.

Most disease mapping models usually belong to the family of generalized linear mixed model, being the Poisson likelihood the most popular choice. If we denote by $\eta_{i,t}$ an appropriate transformation of the disease (mortality or morbidity) rate for area or location $i$ at time $t$, then the typical specification has the form (e.g. Torabi and Rosychuk,, 2010) $\eta_{i,t}=\boldsymbol{\beta}^{\prime}{\bf x}_{i,t}+\theta_{i,t}$, where the first part corresponds to the fixed effects, driven by covariates ${\bf x}_{i,t}$, and the second to the random effects. The random effects are, in turn, expressed in terms of spatial effects, temporal effects and sometimes interaction (spatio-temporal) effects. To be specific $\theta_{i,t}=\alpha_{i}+\delta_{t}+\xi_{i,t}$. Spatial dependencies $\{\alpha_{i}\}$ are routinely incorporated into the covariance structure through normal conditionally auto-regressive (CAR) specifications (Besag et al.,, 1991; Banerjee et al.,, 2003). Temporal dependencies $\{\delta_{t}\}$ are captured by first order autoregressive (AR), or dynamic, normal models (Waller et al.,, 1997). For the interaction effect, proposals are $\xi_{i,t}=\xi_{i}t$ or $\xi_{i,t}=\xi_{i}(t)$, where $\xi_{i}$ is a spatial (CAR) model and $\xi_{i}(t)$ is a temporal spline function for each area $i$ (e.g. MacNab and Dean,, 2001).

Random effects are also assumed to be contaminated by observational noise (Best et al.,, 2005), or to include seasonal effects (Torabi and Rosychuk,, 2010). Multivariate CAR (MCAR) models for the random effects have also been proposed in survival data (Jin and Carlin,, 2005). Posterior inference of disease mapping models under a Bayesian approach were optimised using integrated nested Laplace approximations (Schrödle and Held,, 2010). Additionally, Goicoa et al., (2018) studied identifiability constraints in these kinds of models.

Instead of placing temporal and spatial associations into random effects, they can also be incorporated into the regression coefficients. In other words, the linear predictor would be $\eta_{i,t}=\boldsymbol{\beta}_{i,t}^{\prime}{\bf x}_{i,t}$, where $\boldsymbol{\beta}_{i,t}$ is a set of spatio-temporal varying coefficients. In a continuous setting, Gelfand et al., (2003) used Gaussian processes to define spatial and spatio-temporal varying coefficients with correlations defined via the Matèrn function of the spatial and temporal distances. Gelfand and Vounatsou, (2003) defined multivariate conditional autoregressive processes to cope with spatial association and dependence across multiple regression coefficients. Choi et al., (2012) assumed spatial clusters in which each has a set of homogeneous time-varying coefficients. In contrast, Cai et al., (2013) generalised parametric to non parametric spatial specifications by considering an area-specific Dirichlet process prior. Under a Bayesian framework, the spatial or spatio-temporal formulations provide easy borrowing of information across the whole study region and across time to provide efficient smooth estimates of the overall spatio-temporal risk patterns as well as variance reduction through the use of shrinkage estimators.

In this work we introduce a novel spatio-temporal dependent process through a hierarchical model that relies on latent variables. The process has a joint multivariate normal distributions with very flexible associations that can accommodate spatial, temporal and spatio-temporal interaction dependences. We then use it as a prior distribution for spatio-temporal varying coefficients $\boldsymbol{\beta}_{i,t}$ in a generalised linear regression framework. We use the proposed model to identify space and time variations in our motivating study of the impact of climate variables in the morbidity of gastrointestinal and respiratory diseases.

The contents of the rest of the paper is as follows: In Section 2 we introduce our space-time dependent process and study its association properties. In Section 3 we define the framework of generalised linear regression models with time-varying coefficients and obtain the corresponding posterior distributions when we use our space-time process as prior distribution. Section 4 contains a detail study of the climate variables impact on the morbidity of some diseases in Mexico. A comparison with alternative priors is also included. We conclude with some remarks in Section 5.

Before proceeding we introduce notation. $\mbox{N}(m,c)$ denotes a univariate normal density with mean $m$ and precision $c$ ; $\mbox{N}_{p}({\bf M},{\bf C})$ denotes a $p$-variate normal density with mean vector ${\bf M}_{p\times 1}$ and precision matrix ${\bf C}_{p\times p}$; $\mbox{Ga}(a,b)$ denotes a gamma density with mean $a/b$. The density evaluated at a specific point $x$, will be denoted, for instance for the normal case, as $\mbox{N}(x\mid m,c)$.

2 A novel space-time process

Let $\beta_{i,t}$ be a parameter of interest for area $i$ at time $t$, for $i=1,\ldots,n$ and $t=1,2,\ldots,T$. The idea is to define a dependence structure among the $\boldsymbol{\beta}=\{\beta_{i,t}\}$ in space and time. Since we further want to use this construction as prior distributions for regression coefficients, we also want $\beta_{i,t}$ to have a normal marginal distribution. To achieve this we follow ideas from Jara et al., (2013), Nieto-Barajas and Bandyopadhyay, (2013) and Nieto-Barajas and Huerta, (2017).

Let $\partial_{i,t}$ be the set of “neighbours”, in the broad sense, of area $i$ at time $t$. For spatial dependence we consider actual neighbours, i.e. areas that share a border, plus the current index $i$, whereas for temporal dependence we consider lagged times up to time $t$. Combinations of spatial and temporal dependences are also possible. To be specific,

• $(s)$

Spatial neighbours of areas at the same time, that is

$\partial_{i,t}^{(s)}=\{(j,t):\{j\sim i\}\,\cup\{j=i\}\}$, where “$\sim$” denotes spatial neighbour;

• ($t_{q}$)

Temporal neighbours of order $q>0$, that is

$\partial_{i,t}^{(t_{q})}=\{(i,s):s\in\{t-q,\ldots,t-1,t\}\}$;

• ($s+t_{q}$)

Spatial plus order $q$ temporal neighbours, that is

$\partial_{i,t}^{(s+t_{q})}=\partial_{i,t}^{(s)}\cup\partial_{i,t}^{(t_{q})}$; and

• ($s\times t_{q}$)

Interaction between spatial and temporal dependence or order $q$, i.e., the spatial neighbours are neighbours for all lagged times, that is

$\partial_{i,t}^{(s\times t_{q})}=\{(j,s):\{j\sim i\}\,\cup\{j=i\};s\in\{t-q,\ldots,t-1,t\}\}$.

Alternative definitions of the set $\partial_{i,t}$ are also possible, for instance seasonal or periodic temporal dependencies (e.g. Jara et al.,, 2013) or second order spatial dependencies. Note that spatial neighbours are reciprocal, that is, $i\sim j$ iff $j\sim i$, whereas temporal neighbours are directed, that is, $s\rightarrow t$ does not imply that $t\rightarrow s$, because we may expect a variable at time $t$ to depend on past values (lagged times $s$), but not the other way around.

For each area $i$ at time $t$ we require a latent parameter $\gamma_{i,t}$, plus a common parameter $\omega$. Let $\boldsymbol{\gamma}=\{\gamma_{i,t}\}$, then the proposed model for $\boldsymbol{\beta}=\{\beta_{i,t}\}$ is defined through a three level hierarchical model of the form

[TABLE]

where $m_{0}\in{\rm I}\!{\rm R}$, $c_{0}>0$ and $c_{i,t}>0$ for $i=1,\ldots,n$ and $t=1,2,\ldots,T$ are known hyper-parameters, and $\partial_{i,t}$ is a set of neighbours which could be any of those defined above. We denote construction (1) as $\mbox{STN}(m_{0},c_{0},{\bf c})$, where ${\bf c}=\{c_{i,t}\}$ parameters determine the importance of location $(i,t)$ in the definition of the net, and thus the degree of dependence among $\{\beta_{i,t}\}$. Properties of this construction are given in the following proposition.

Proposition 1

Let $\boldsymbol{\beta}\sim\mbox{STN}(m_{0},c_{0},{\bf c})$, that is $\boldsymbol{\beta}=\{\beta_{i,t}\}$ for $i=1,\ldots,n$ and $t=1,2,\ldots,T$ is a set of parameters whose joint distribution is defined by (1). Then, $\beta_{i,t}\sim\mbox{N}(m_{0},c_{0})$ marginally for all $i$ and $t$. Moreover, $\boldsymbol{\beta}\sim\mbox{N}_{nT}({\bf M}_{0},{\bf C})$, that is, $\boldsymbol{\beta}$ has a multivariate normal distribution of dimension $nT$ with mean vector ${\bf M}_{0}=m_{0}{\bf 1}_{nT}$, and ${\bf 1}$ a vector of 1’s, and covariance matrix ${\bf C}^{-1}$ of dimension $nT\times nT$ with diagonal elements $1/c_{0}$ and off-diagonal elements appropriately defined by correlations between any two $\beta_{i,t}$ and $\beta_{j,s}$ given by

[TABLE]

where $\partial_{i,t}$ is a set of neighbours like those defined at the beginning of this Section.

Proof To prove the marginal distribution we note that, conditionally on $\omega$, dropping the summation indexes, $\sum c_{j,s}\gamma_{j,s}\mid\omega\sim\mbox{N}\left(\omega\sum c_{j,s}\,,1/(\sum c_{j,s})\right)$ and integrating with respect to $\omega$ we get $\sum c_{j,s}\gamma_{j,s}\sim\mbox{N}\left(m_{0}\sum c_{j,s}\,,c_{0}/\{(\sum c_{j,s})(c_{0}+\sum c_{j,s})\}\right)$. Now relying on conjugacy properties of the normal model, or simply integrating with respect to $\boldsymbol{\gamma}$ in the first equation of (1), we obtain the result. For the multivariate normality, we note that $\gamma_{i,t}$ given $\omega$ are conditionally independent, then jointly $\boldsymbol{\gamma}$, after marginalizing $\omega$, are multivariate normal with common covariance $1/c_{0}$. Finally, using the same argument, $\beta_{i,t}$’s are conditionally independent given $\boldsymbol{\gamma}$, after marginalizing $\boldsymbol{\gamma}$ we obtain that $\boldsymbol{\beta}$ are jointly multivariate normal. To obtain the correlation, we first obtain the covariance using conditional expectation properties and then standardise it using the marginal variance. $\diamond$

The correlation expression in Proposition 1 has a nice interpretation, the first term in the numerator is a function of the common parameters $c_{k,r}$ that appear in the definition of both $\beta_{i,t}$ and $\beta_{j,s}$ and the second is a function of all dependence parameters in each $\beta_{i,t}$ and $\beta_{j,s}$. In other words, two locations that share the same neighbours will have a higher correlation, even though they are not direct neighbours, and two locations that do not share any neighbour will still be correlated. Additionally, a larger value of $c_{i,t}$ for location $(i,t)$ will make the correlation larger for all pairs of locations that share the same parameter $c_{i,t}$, therefore $c_{i,t}$ can also be seen as a measure of the importance of location $(i,t)$ in the whole network.

Let us consider a single location so we can concentrate on temporal associations. Let $\boldsymbol{\beta}=\{\beta_{t}\}$ for times $t=1,\ldots,T$ and assume $q=1$. From Proposition 1, our STN model implies $\mbox{Var}(\beta_{t})=1/c_{0}$ for all $t$, and correlations for time $t=2$, $\mbox{Corr}(\beta_{2},\beta_{3})=\frac{c_{0}c_{2}+(c_{1}+c_{2})(c_{2}+c_{3})}{(c_{0}+c_{1}+c_{2})(c_{0}+c_{2}+c_{3})}$ and $\mbox{Corr}(\beta_{2},\beta_{s})=\frac{(c_{1}+c_{2})(c_{s-1}+c_{s})}{(c_{0}+c_{1}+c_{2})(c_{0}+c_{s-1}+c_{s})}$ for $s\geq 4$. If all the $c_{t}$’s had the same value, the correlation between $\beta_{2}$ and $\beta_{3}$ would be higher than that between $\beta_{2}$ and any $\beta_{s}$ for $s\geq 4$. In general, the correlation can remain high between times $t$ and $s$ even if $|t-s|>q$. This is in contrast to temporal AR models where the correlation decays exponentially to zero as $|t-s|$ increases.

Now, if we only consider a single time, we can concentrate in spatial associations. Let $\boldsymbol{\beta}=\{\beta_{i}\}$ for areas $i=1,\ldots,n$, then STN model (1) can be considered as a flexible alternative to the conditional autoregressive (CAR) model. To see this we recall the definition of a CAR model. This is defined through a multivariate normal distribution with mean vector zero and precision matrix $\tau(D_{W}-\rho W)$, where $W=(w_{ij})$ is the spatial neighbourhood matrix such that $w_{ij}=I(i\sim j)$ and $D_{W}=\mbox{diag}(w_{1+},\ldots,w_{n+})$ with $w_{i+}=\sum_{j=1}^{n}w_{ij}$. To better illustrate the differences between STN and CAR models, let us consider a toy example that consists of a country with only three areas $n=3$, with a neighbourhood structure given by

[TABLE]

From Proposition 1, STN model implies a common variance, $\mbox{Var}(\beta_{i})=1/c_{0}$ for $i=1,2,3$, and correlations given by $\mbox{Corr}(\beta_{1},\beta_{2})=\frac{c_{0}(c_{1}+c_{2})+(c_{1}+c_{2})(c_{1}+c_{2}+c_{3})}{(c_{0}+c_{1}+c_{2})(c_{0}+c_{1}+c_{2}+c_{3})}$, $\mbox{Corr}(\beta_{1},\beta_{3})=\frac{c_{0}c_{2}+(c_{1}+c_{2})(c_{2}+c_{3})}{(c_{0}+c_{1}+c_{2})(c_{0}+c_{2}+c_{3})}$ and $\mbox{Corr}(\beta_{2},\beta_{3})=\frac{c_{0}(c_{2}+c_{3})+(c_{1}+c_{2}+c_{3})(c_{2}+c_{3})}{(c_{0}+c_{1}+c_{2}+c_{3})(c_{0}+c_{2}+c_{3})}$. On the other hand, the CAR model implies different variances, $\mbox{Var}(\beta_{1})=\mbox{Var}(\beta_{3})=\frac{2-\rho^{2}}{2\tau(1-\rho^{2})}$ and $\mbox{Var}(\beta_{2})=\frac{1}{2\tau(1-\rho^{2})}$, and correlations given by $\mbox{Corr}(\beta_{1},\beta_{2})=\mbox{Corr}(\beta_{2},\beta_{3})=\frac{\rho}{\sqrt{2-\rho^{2}}}$ and $\mbox{Corr}(\beta_{1},\beta_{3})=\frac{\rho^{2}}{2-\rho^{2}}$. In our proposed model, spatial dependence is controlled by the set of parameters $\{c_{0},c_{1},c_{2},c_{3}\}$, whereas in the CAR model it is controlled by a single parameter $\rho$. In both cases, the correlation between two non-neighbouring areas remains positive, however in the CAR model the correlation becomes smaller for further apart areas, whereas in our model the correlation could remain relatively high.

These features become relevant when we use any of these spatial and /or temporal models as prior distribution for space-time varying coefficients. STN borrows more strength from further apart neighbours, whereas CAR and AR models models mainly borrow strength from close neighbours. Additionally, it might not be justifiable to impose a different prior variance for the coefficients, as the CAR model does.

3 Regression models

Let $Y_{i,t}$ be a response variable and ${\bf X}_{i,t}^{\prime}=(X_{i,t,1},\ldots,X_{i,t,p})$ a set of $p$ covariates for area $i$ at time $t$, for $i=1,\ldots,n$ and $t=1,2,\ldots,T$. We consider a generalised linear model framework (e.g. McCullagh and Nelder,, 1989), so that the response variable $Y_{i,t}$, conditionally on the explanatory variables ${\bf X}_{i,t}$, has a density $f(y_{i,t}\mid{\bf x}_{it})$ which is a member of the exponential family.

We model the conditional expectation in terms of the explanatory variables as $\mbox{E}(Y_{i,t}\mid{\bf X}_{it})=g^{-1}(\eta_{i,t})$, where $g(\cdot)$ is an appropriate link function and $\eta_{i,t}$ is a linear predictor, which for static coefficients has the form $\eta_{i,t}=\alpha+\boldsymbol{\beta}^{\prime}{\bf x}_{i,t}$. We introduce a space-time dynamic in the regression coefficients and define the linear predictor as

[TABLE]

To state the prior for the model parameters we expand the inner product in the linear predictor and write $\eta_{i,t}=\alpha+\beta_{i,t,1}x_{i,t,1}+\cdots+\beta_{i,t,p}x_{i,t,p}$ and take $\alpha\sim\mbox{N}(m_{\alpha},c_{\alpha})$ for the intercept and for the regression coefficients of each explanatory variable $k$ we take a normal space-time processes, that is, $\boldsymbol{\beta}_{k}=\{\beta_{i,t,k}\}\sim\mbox{STN}(m_{0},c_{0},{\bf c})$ defined by (1), independently for $k=1,\ldots,p$.

The likelihood for the regression model is simply $f({\bf y}\mid{\bf x})=\prod_{i}\prod_{t}f(y_{i,t}\mid{\bf x}_{i,t})$, and the joint prior for all parameters, including the latent ones, has the form

[TABLE]

[TABLE]

Alternatively, after integrating with respect to the latent parameters $(\boldsymbol{\gamma},\omega)$, the prior becomes $f(\alpha,\boldsymbol{\beta})=\mbox{N}(\alpha\mid m_{\alpha},c_{\alpha})\prod_{k=1}^{p}\mbox{N}_{nT}(\boldsymbol{\beta}_{k}\mid{\bf M}_{0},{\bf C})$.

To specify the posterior distributions induced, let us consider a normal likelihood for the regression model of the form

[TABLE]

for $i=1,\ldots,n$ and $t=1,\ldots,T$, where $\mu({\bf x}_{i,t})=\eta_{i,t}$ is the mean defined by an identity function of the linear predictor (3), and $\tau_{i,t}$ is the precision parameter. In this case, the conditional posterior distribution for each $\boldsymbol{\beta}_{k}$ becomes a multivariate normal of dimension $nT$. For large $n$ and $T$, which is the case of our motivating example where $nT=32\times 60=1920$, dealing with a multivariate normal of such dimension is not computationally feasible. Alternatively, we suggest to sample from the univariate conditional posterior distributions for each $\beta_{i,t,k}$ given the latent parameters and the data. The full conditional distributions for all model parameters are:

1. i)

$f(\alpha\mid{\bf y},\mbox{rest})=\mbox{N}(\alpha\mid\mu_{\alpha},\tau_{\alpha})$, where

[TABLE] 2. ii)

$f(\beta_{i,t,k}\mid{\bf y},\mbox{rest})=\mbox{N}(\beta_{i,t,k}\mid\mu_{\beta},\tau_{\beta})$, where

[TABLE]

and

[TABLE]

for $i=1,\ldots,n$, $t=1,\ldots,T$ and $k=1,\ldots,p$. 3. iii)

$f(\gamma_{i,t,k}\mid{\bf y},\mbox{rest})=\mbox{N}(\gamma_{i,t,k}\mid\mu_{\gamma},\tau_{\gamma})$, where

[TABLE]

and

[TABLE]

where $\varrho_{i,t}$ is the set of reversed neighbours, that is, the set of pairs $(j,s)$ such that $(i,t)\in\partial_{j,s}$, for $i=1,\ldots,n$, $t=1,\ldots,T$ and $k=1,\ldots,p$. 4. iv)

$f(\omega_{k}\mid{\bf y},\mbox{rest})=\mbox{N}(\omega_{k}\mid\mu_{\omega},\tau_{\omega})$, where

[TABLE]

for $k=1,\ldots,p$.

Finally, we assume $\tau_{i,t}=\tau_{i}$ to have a common precision along time for each area $i$. If we further take $\tau_{i}\sim\mbox{Ga}(a_{0},b_{0})$ a-priori then its conditional posterior distribution has the form

1. v)

$f(\tau_{i}\mid{\bf y},\mbox{rest})=\mbox{Ga}(\tau_{i}\mid a_{\tau},b_{\tau})$, where

[TABLE]

for $i=1,\ldots,n$.

Posterior inference of the model parameters relies on a Gibbs sampler (Smith and Roberts,, 1993). Conditional posterior distributions i)–v) are of standard form because we are dealing with a normal regression model for which our normal and space-time normal prior distributions are conditionally conjugate. For other likelihoods different to the normal, say Bernoulli, Poisson or gamma, distributions i) and ii) must be appropriately modified and a Metropolis-Hastings step (Tierney,, 1994) would be required to sample from them.

4 Numerical analyses

4.1 Simulation studies

We carry out two simulation studies to test the performance of our posterior simulation procedure. We consider a simple spatial setting with $n=3$ regions as in (2), plus the temporal dimension. The first scenario consists in fixing the regression coefficients $\beta_{i,t}$ in a deterministic way as follows: $\beta_{1,t}=(1.1)^{t}$, $\beta_{2,t}=1+(1.01)^{t}$ and $\beta_{3,t}=(0.9)^{t}$, for $t=1,\ldots,T$ with $T=10$. We then sample response variables from model (4) with $\alpha=-1$ and $\tau_{i,t}=\tau_{i}=1$ and a single covariate $x_{i,t}\sim\mbox{Un}(0,1)$ to form a sample of size $nT=30$. This experiment was repeated 100 times. The prior distribution was defined by $m_{\alpha}=0$, $c_{\alpha}=0.01$, $m_{0}=0$, $c_{0}=0.01$, $c_{i,t}=0.1$ and $a_{0}=b_{0}=0.01$. The sets of neighbours were those defined in Section 3 with $q=1,2$. Posterior inference was obtained through a Gibbs sample with two parallel chains of 21,000 iterations with a burn-in of 1,000 and keeping one of every 10th iteration.

For model comparison we compute two goodness of fit (gof) measures, the logarithm of the pseudo marginal likelihood (LPML) and the deviance information criterion (DIC). LPML is a summary measure of conditional predictive ordinates (CPO), commonly used for model comparison and introduced by Geisser and Eddy, (1979). Given posterior samples of model parameters, $\alpha^{(r)}$, $\beta_{i,t,k}^{(r)}$ and $\tau_{i}^{(r)}$ for $r=1,\ldots,R$, a Monte Carlo estimate $\widehat{\mbox{CPO}}_{i,t}$, for each data point $(i,t)$, is obtained as

[TABLE]

for $i=1,\ldots,n$ and $t=1,\ldots,T$. Finally these values are summarised to define $\mbox{LPML}=\sum_{i=1}^{n}\sum_{t=1}^{T}\log(\mbox{CPO}_{i,t})$. On the other hand, DIC was introduced by Spiegelhalter et al., (2002) and is a model selection criterion that penalises for model complexity. Larger/smaller values of LPML/DIC measures are preferable.

Table 1 reports the average LPML and DIC across the 100 repetitions of the experiment. Looking at the third and fourth columns we observe that the best fitting is obtained with a temporal neighbourhood structure $(t_{q})$ with $q=2$. This makes sense since there is no relation among spatial neighbours in the definition of the coefficients. For this wining model, the average coverage of 95% credible intervals (CI) for all $\beta_{i,t}$ is, in average, 93.3%. Figure 1 includes posterior estimates (averaged across the experiments). Our estimates capture the temporal trends in the three regions.

The second scenario consists in random coefficients simulated from model (1) with parameters $m_{0}=0$, $c_{0}=1$, $c_{i,t}=1$ and neighbourhood structure $(s+t_{q})$ with $q=1$. The data were generated as in the first scenario and 100 repetitions of the experiment were also obtained. The objective of this study is to see whether the gof measures are able to detect the correct neighbourhood structure. Therefore we took the same specifications in the prior as those used to simulate the coefficients $\beta_{i,t}$, plus $m_{\alpha}=0$, $c_{\alpha}=0.01$ for $\alpha$, and $a_{0}=b_{0}=0.01$ for $\tau_{i}$. We tried different sets of neighbours, as those defined in Section 3 with $q=1,2$, and the same Gibbs sampler specifications as above.

Results are shown in the fifth and sixth columns in Table 1. The best model chosen by both the LPML and the DIC is that with neighbourhood $(s+t_{q})$ and $q=1$ as it should be.

4.2 Climate impact on morbidity

In general, the climate in a region is characterised by the atmospheric conditions of temperature and pluvial precipitation, among other factors. The Mexican Autonomous National University (UNAM) through the group of Climate and Society of the Center for Atmospheric Sciences has created a monthly database of climate records for all 32 states of Mexico from 1901 to the date. Measured variables are average and maximum temperature, both in centigrades degrees, and average pluvial precipitation in millimetres.

Public health records in Mexico are obtained through the National System for Epidemiology Surveillance (SINAVE). This system gathers information from the whole country and all institutions of the national system of health. Information is available from 1985 to date for all 32 states of Mexico. We will concentrate in two big groups of diseases, gastrointestinal and respiratory. Each of this groups is formed by several specific diseases, say: amebiasis, cholera, typhoid, intestinal infection, intoxication, salmonellosis and shigellosis, for the gastrointestinal group; and asthma, pharyngitis, acute respiratory infections, influenza and pneumonia, for the respiratory group.

Since the information of all specific diseases is not available for the same years, we will concentrate on a window of five years, from 2011 to 2015 (60 months) for the two disease groups. Additionally, we will require two more variables, the percentage of illiterate people as a marginality indicator, to capture the trend, and the population size (number of inhabitants) as an offset.

We now apply our Bayesian regression model, with spatio-temporal varying coefficients and STN prior, to the motivating study of climate impact on morbidity of gastrointestinal and respiratory diseases in Mexico, mentioned in Section 1.

To state the model we define variables: $Y_{i,t}=\log\left(N_{i,t}/P_{i,t}\right)$ is the morbidity rate (in log scale), where $N_{i,t}$ is number of disease cases and $P_{i,t}$ is the population size; $X_{i,t,1}$ is the average pluvial precipitation (in log scale); $X_{i,t,2}$ is the average temperature (in log scale); and $X_{i,t,3}$ is the percentage of illiterate people (in log scale), for state $i=1,\ldots,n$ and month $t=1,\ldots,T$, with $n=32$ and $T=60$. The maximum temperature will not be used in the model due to a high correlation with the average temperature. The model is therefore as in equations (3) and (4) with $p=3$ and $\tau_{i,t}=\tau_{i}$.

Considering the 32 states and the 60 months we have around 28 million cases registered in the gastrointestinal group and 133 million cases in the respiratory group. From this total, the percentage of cases for each state is presented in Figure 3. We produced this graph for both, the gastrointestinal and the respiratory groups, however they were almost identical, confirming that the number of cases is a function of the number of people at risk (population size). Therefore we only present the graph for the gastrointestinal group. The state with the largest proportion of cases, is the State of Mexico ($i=15$) with around 11% of the cases, followed by Mexico City ($i=9$) with around 7% of the cases. The smallest states, in terms of population size are Baja California Sur ($i=3$) and Colima ($i=6$), which present the smallest proportion with less than 1% of the cases.

The variables involved in the model are presented in Figure 4 as time series. In the panels we show $Y_{i,t}$ for the gastrointestinal and respiratory groups in the top row, and the three explanatory variables $X_{i,t,k}$ for $k=1,2,3$ in the bottom row. Apart from the illiteracy proportion (bottom right panel) which shows a decreasing trend, the rest of the variables present a 12 months seasonal pattern. We also note that the log rates for gastrointestinal and respiratory groups are shifted a period of 6 months. This is explained by the fact that respiratory infections have a peak in the winter, whereas gastrointestinal cases occur more often in the summer.

To specify the model, we consider the four sets of neighbours defined in Section 2, $(s)$, $(t_{q})$, $(s+t_{q})$ and $(s\times t_{q})$. In particular we took temporal neighbours that define a short term dependence, lags up to a quarter year, or a medium term dependence, lags up to a semester, i.e. $q\in\{3,6\}$ to compare; and spatial neighbours given by the geographical adjacencies of the 32 states of Mexico, which are included in Table 2. The number of adjacencies per state goes from 1 to 8 with a median value of 4.

Prior distributions for our model parameters were defined by: $m_{\alpha}=0$ and $c_{\alpha}=0.01$, for $\alpha$; $m_{0}=0$, $c_{0}=0.01$ and constant $c_{i,t}$ $\forall i,t$ with values in $\{10,50,100\}$ to define low, medium and high prior dependence, for $\boldsymbol{\beta}_{k}$; and $a_{0}=b_{0}=0.01$, for $\tau_{i}$. Posterior inference requires the implementation of a Gibbs sampler with a very large number of parameters, therefore we expect to experience high autocorrelation in the chains. We ran two parallel chains for $115,000$ iterations with a burn-in of $15,000$ and kept one of every $50^{th}$ iteration. Convergence of the chains was assessed by monitoring the trace plots, the ergodic means, as well as the autocorrelation function. The code is written in Fortran and is available as Supplementary Material.

For comparison purposes, we considered three competing priors: $\boldsymbol{\beta}_{t,k}=\{\beta_{i,t,k},\,i=1,\ldots,n\}\sim\mbox{N}_{n}({\bf 0},\tau_{c}(D_{W}-\rho W))$, i.e. CAR priors to account for spatial dependence; $\boldsymbol{\beta}_{i,k}=\{\beta_{i,t,k}\,t=1,\ldots,T\}\sim\mbox{N}_{T}({\bf 0},\Sigma^{-1})$, i.e. stationary first order AR processes, to account for temporal dependence, defined by covariance elements $\sigma_{t,s}=\varphi^{|t-s|}/\tau_{ar}$; and $\boldsymbol{\beta}_{k}=\{\beta_{i,t,k},\,i=1,\ldots,n,\,t=1,\ldots,T\}\sim\mbox{N}_{nT}({\bf 0},(D_{W}-\rho W)\otimes\Sigma^{-1})$, i.e. MCAR priors (Gelfand and Vounatsou,, 2003) to account for space and time dependence, with covariance matrix $\Sigma$ defined as in the first order AR process. These priors were specified by taking $\tau_{c}=\tau_{ar}=0.1$ and $\rho=\varphi=0.99$.

Additionally, we also considered two baseline models: Model 0, which assumes a different regression coefficient for each individual $i$ but common for all times $t$. In notation, the linear predictor for this model is $\eta_{i,t}=\alpha+\boldsymbol{\beta}_{i}^{\prime}{\bf x}_{i,t}$ and the precision is $\tau_{i}$; and model 00 that assumes a common regression coefficient for all individuals $i$ and for all times $t$. This is obtained with linear predictor $\eta_{i,t}=\alpha+\boldsymbol{\beta}^{\prime}{\bf x}_{i,t}$ and precision $\tau$. For these two models we took $\mbox{N}(0,0.01)$ independent prior distributions for all parameters $\alpha$, $\beta_{i,k}$ and $\beta_{k}$, respectively, and $\mbox{Ga}(0.01,0.01)$ priors for the precisions $\tau_{i}$ and $\tau$, respectively.

We assess model fit by computing the two gof measures, LPML and DIC. These are reported in Table 3 for the 21 different versions of our model, the three competitors, plus models 0 and 00, all fitted to both datasets, gastrointestinal and respiratory. Additionally, Table 3 reports the median size of the neighbourhoods, $\mbox{med}|\partial_{i,t}|$, that define each of our models. For spatial models $(s)$ the neighbourhood size ranges from 2 to 9, so the median size reported is 5. For temporal models $(t_{q})$ the exact neighbourhood size is $q+1$.

In all cases, the two gof measures are in agreement in selecting the best model. For both datasets, the worst model is model 00, however model 0 is better than some versions of our space-time model. In general we note that models with larger dependence parameters $c_{i,t}$ equal to 50 or 100 achieve better fitting in models with small neighbourhood size, which is the case for spatial $(s)$, temporal $(t_{q})$, and spatial plus temporal $(s+t_{q})$ models. However, for models with space-time interaction $(s\times t_{q})$, which have larger neighbourhood size, the smaller value $c_{i,t}=10$ is preferred. In other words, as the number of neighbours increases, both datasets prefer models with smaller values of the dependence parameters.

Considering the best spatial model $(s)$ with $c_{i,t,}=50$, we obtain a better fit than model 0, for gastrointestinal data, but a worst fit for respiratory data. This tell us that, somehow, gastrointestinal diseases are more spatially dependent than respiratory diseases. On the other hand, the best temporal model $(t_{q})$ for gastrointestinal data is obtained with a lag of order $q=3$ and a dependence parameter $c_{i,t}=100$, whereas for respiratory data it is obtained with $q=6$ and $c_{i,t}=50$. In both cases, the best temporal model is better than model 0.

Now, comparing the space and time models $(s+t_{q})$ and $(s\times t_{q})$, the former achieves better fit. In fact, for both datasets the best model, overall, is obtained with a neighbourhood that considers spatial plus medium term neighbours in time (order $q=6$), and with a medium strength prior dependence ($c_{i,t}=50$). In other words, the best model is choosing a neighbourhood structure for the pair $(i,t)$ that depends on the same state $i$ in the previous 6 months and on the geographical adjacent states for the same time $t$. This model has a median neighbourhood size of 11 and obtains gof statistics of, LPML$=974$ and DIC$=-2199$, for the gastrointestinal data, and LPML$=982$ and DIC$=-2148$, for the respiratory data.

An alternative way of assessing the gain in the fitting, we computed the residuals $y_{i,t}-E(y_{i,t}\mid\mbox{data})$ for model 00 and our best fitting model STN. Box plots of them are presented in Figure 2. For both datasets the variance not explained by the model is highly reduced when using our STN model.

Analysing the gof values of the three competing priors, we note that none of them achieve a good fitting. AR prior behaves better than the CAR prior, and both priors combined in the MCAR is the best among these three. The performance of MCAR is compared to that of our model with only spatial dependence with association parameters $c_{i,t}=10$, for both datasets.

For reference purposes, posterior 95% credible intervals for the model parameters under model 00 are: $\alpha\in(-7.91,-7.59)$, $\beta_{1}\in(-0.018,-0.006)$, $\beta_{2}\in(0.79,0.91)$, $\beta_{3}\in(-0.18,-0.13)$, for gastrointestinal data; and $\alpha\in(-2.44,-2.04)$, $\beta_{1}\in(-0.04,-0.03)$, $\beta_{2}\in(-0.57,-0.44)$, $\beta_{3}\in(-0.04,0.01)$, for respiratory data. From these numbers we can say that pluvial precipitation (${\bf X}_{1}$) has a negative effect in both disease groups, whereas temperature (${\bf X}_{2}$) has a different effect, it is positive for gastrointestinal data and negative for respiratory data, which makes sense. Finally, the percentage of illiterate people (${\bf X}_{3}$), has a negative effect for gastrointestinal disease, and shows no effect for respiratory disease. Since the percentage of illiterate people shows a decreasing trend in time, a negative effect means that gastrointestinal cases have a positive trend in time. We emphasize that the illiteracy might not have a direct relationship with the diseases, its inclusion in the analysis is merely to capture trends.

Interpreting the coefficients for the best models in both datasets, posterior 95% CI for the intercept are: $\alpha\in(-7.78,-7.42)$ for the gastrointestinal data; and $\alpha\in(-2.19,-1.69)$ for the respiratory data. These values are in accordance with those obtained from model 00. Figure 5 contains posterior estimates for the precision parameters $\tau_{i}$, $i=1,\ldots,32$, where the dots correspond to the mean and the vertical lines to 95% CI. The left panel corresponds to gastrointestinal data and the right panel to respiratory data. In both diseases, point estimates for the precisions lie between 100 and 250 with some exceptions, Jalisco ($i=14$) in the gastrointestinal case, and Baja California ($i=2$) in the respiratory case, whose precisions are lower. There is not particular reason for this to happen, other than the fact that the variability not explained by the model, in these two specific states, is larger than for the other states.

For the gastrointestinal data, posterior estimates of $\boldsymbol{\beta}_{1}$ are shown in Figure 6, where we include 32 panels, one for each state, with time series for $t=1,\ldots,60$. Reported are point estimates (solid line) and 95% CI (shadows). We can see that the estimates are not constant and vary across $i$ and $t$. However, the effect of pluvial precipitation is mainly no significant since most CI’s contain the value of zero, with some few exceptions for specific states and specific times. This is the case of Aguascalientes ($i=1$) whose effect shows an increasing trend, but only becomes significant in the rainy season (June, July and August) of 2015, where a 10% increment in pluvial precipitation produced around 1% increment in the disease rate. This is in contrast to Veracruz ($i=30$) where a 10% increment in pluvial precipitation produced around 1% reduction in the disease rates in the last quarter of the years 2013, 2014 and 2015.

On the other hand, Figure 7 reports posterior estimates of $\boldsymbol{\beta}_{2}$. Again, the estimates are not constant, and for all states and times, posterior CI’s only contain positive values, implying a positive significant effect of temperature on the gastrointestinal rates. For all states, posterior means have values between 0.7 and 0.9, which is consistent with the common effect of model 00. For some states, the effect fluctuates more, as is the case of Sonora ($i=26$), and for some others the effect is more steady, as is the case of Tabasco ($i=27$). However, temperature effects are not the same in all states. A 10% increment in the temperature produced around a 9% increment in the gastrointestinal rates in Aguscalientes ($i=1$), and around a 7% increment in Baja California ($i=2$).

Finally, Figure 8 includes posterior estimates of $\boldsymbol{\beta}_{3}$, loosely speaking we can say that the effect of the illiteracy indicator is negative significant for all states, perhaps for some states the credibility should be 90% instead of 95%. The point estimates lie between $-0.1$ and $-0.2$, which is also consistent with the common effect of model 00. This implies that there is a slight positive trend in the gastrointestinal rates in all states.

For the respiratory data, posterior estimates of $\boldsymbol{\beta}_{1}$ are included in Figure 9. There are only three states with a negative significant effect of pluvial precipitation for all months, these are Chiapas ($i=7$), Puebla ($i=21$) and Veracruz ($i=30$), where a 10% increment in the pluvial precipitation produced a 1% decrement in the disease rates. For other states like the State of Mexico ($i=15$) and San Luis Potosi ($i=24$), only the rain season (June, July and August) shows a negative significant effect of the same amount as in the former three states. For the rest of the states there is no effect.

Posterior estimates of $\boldsymbol{\beta}_{2}$ are shown in Figure 10. It is clear that the temperature has a negative significant effect for all states and all months in the respiratory rates. Posterior means lie between $-0.6$ and $-0.7$, consistent with the common effect estimate of model 00. In particular, a 10% increment in the temperature produced around 6% decrement in Aguascalientes ($i=1$) and around a 7.5% decrement in Baja California ($i=2$), in the respiratory rates.

Finally, posterior estimates of $\boldsymbol{\beta}_{3}$ are given in Figure 11. For most states there is no significant effect of the illiteracy indicator in the respiratory rates, perhaps the exceptions are Hidalgo ($i=13$), Sinaloa ($i=25$) and Zacatecas ($i=32$) which show a positive significant effect, implying that there is a negative trend in the respiratory rates in these states.

5 Concluding remarks

We have proposed a space-time dependent process, STN, which has a multivariate normal distribution and identically distributed marginal distributions. We use our process as a prior distribution for the coefficients in a regression model. The prior is very flexible so it allows us to identify areas and times where explanatory variables show differentiated effects, as was the case in the study of gastrointestinal and respiratory diseases in Mexico.

When comparing our model with competitors previously proposed in the literature, our model outperforms all of them. We believe that none of the typical priors (CAR, AR and MCAR), borrow enough strength from spatial and temporal neighbours. The reason is that there is a single data point ($(y_{i,t},x_{i,t})$) to update the coefficients $\alpha$ and $\beta_{i,t,k}$, $k=1,\ldots,p$. A more structured prior as the one we are proposing is preferred in this varying coefficients regression setting.

Other use of our space-time process is to model responses with spatial and temporal dependence. In this case the model would be used as a sampling distribution (likelihood) instead of a prior. Additionally, to make our model even more flexible, instead of setting a value for each $c_{i,t}$ as we did here, we could put a hyper-prior (hierarchical) distribution on all of them so the data can help us determine their best value.

An alternative construction to define a spatial process that gives a different strength to each connection is the following. Instead of defining a latent parameter for each location $i$, we define a latent parameter for each connection. That is, if location $i$ is neighbour of location $j$, in notation $i\sim j$, then we define a latent $\gamma_{j}^{i}$. Of course we do require a symmetry condition such that $\gamma_{j}^{i}\equiv\gamma_{i}^{j}$. If we now let $\boldsymbol{\gamma}=\{\gamma_{j}^{i}\}$, a different model for the distribution of spatial $\boldsymbol{\beta}=\{\beta_{i}\}$ would be

[TABLE]

where $m_{0}\in{\rm I}\!{\rm R}$, $c_{0}>0$ and $c_{i}^{j}>0$ if $i\sim j$, and $c_{i}^{j}=0$ otherwise, similarly $\gamma_{i}^{j}=0$ with probability one if $i\not\sim j$, for $i,j=1,\ldots,n$, and $\partial_{i}$ a set of spatial neighbours, as that given in (s) at the beginning of Section 3. Now, parameters ${\bf c}=\{c_{i}^{j}\}$ determine the strength of dependence between locations $i$ and $j$. It can be proven that the joint distribution of $\boldsymbol{\beta}$ is a multivariate normal with mean $m_{0}{\bf 1}_{n}$ and correlation between any two $\beta_{i}$ and $\beta_{j}$ for $i\neq j$ given by

[TABLE]

which clarifies that a larger value of $c_{j}^{i}$ will induce a larger association between neighbours $\beta_{i}$ and $\beta_{j}$. Although the pairwise dependence is more flexible due to the existence of a linking parameter $c_{j}^{i}$, the correlation between two regions that are second order neighbours (i.e., are not neighbours but share a common neighbour) will be a lot less than that induce by the first model (1). Moreover, extending this spatial process to a spatio-temporal process is not trivial due to the lack of symmetry in the temporal dependence.

Acknowledgements

This research was supported by Asociación Mexicana de Cultura, A.C.–Mexico. The author is grateful to the United Nations Development Programme in Mexico and to two anonymous referees for their insightful comments.