Alleviating confounding in spatio-temporal areal models with an application on crimes against women in India

Abstract

Assessing associations between a response of interest and a set of covariates in spatial areal models is the leitmotiv of ecological regression. However, the presence of spatially correlated random effects can mask or even bias estimates of such associations due to confounding effects if they are not carefully handled. Though potentially harmful, confounding issues have often been ignored in practice leading to wrong conclusions about the underlying associations between the response and the covariates. In spatio-temporal areal models, the temporal dimension may emerge as a new source of confounding, and the problem may be even worse. In this work, we propose two approaches to deal with confounding of fixed effects by spatial and temporal random effects, while obtaining good model predictions. In particular, restricted regression and an apparently—though in fact not—equivalent procedure using constraints are proposed within both fully Bayes and empirical Bayes approaches. The methods are compared in terms of fixed-effect estimates and model selection criteria. The techniques are used to assess the association between dowry deaths and certain socio-demographic covariates in the districts of Uttar Pradesh, India.

1 Introduction

Spatial and spatio-temporal disease mapping techniques have been widely used in epidemiology and public health. Though these analyses are somewhat descriptive, they have undoubted value as they provide information about the geographical pattern of the disease, how this pattern evolves in time, and where regions with extreme risk (high or low) are located. An overall spatio-temporal view of the phenomenon under study is extremely useful for generating hypotheses about factors that may be associated with the disease. Throughout this article and to avoid any misleading conclusion, a risk factor is simply a predictor of one outcome (i.e., cancer, or crimes against women in this article), which can be useful for identifying those areas with high risk requiring special actions or intervention programmes.

When covariates related to the study question are unknown, spatio-temporal areal models are a first step in developing an understanding of the disease or crime under study. Incorporating potential risk factors in a model is usually known as ecological regression, and it confers an inferential perspective on spatio-temporal areal models as it quantifies the relationship between a response and covariates (see, e.g., Martínez-Beneito and Botella-Rocamora (2019)).

Spatio-temporal areal models are practical and valuable tools, but they are not free from inconveniences. Goicoa et al.(2018)Goicoa, Adin, Ugarte, and Hodges highlight some identifiability problems involving the intercept and the spatial and temporal main effects, and involving the main effects and the interaction. To overcome these problems, these authors propose to reparameterize the models or to use constraints, though other approaches exist (see, e.g., Fahrmeir and Lang (2000); Kneib et al.(2019)Kneib, Klein, Lang, and Umlauf for identifiability issues in a class of generalized additive models). Another key issue in spatial and spatio-temporal areal models is potential confounding between the fixed effects and random effects. Spatial confounding occurs when the covariates have a spatial pattern and are collinear with the spatial random effects. Although some authors warn about its effects (see, e.g., Clayton et al.(1993)Clayton, Bernardinelli, and Montomoli; Zadnik and Reich (2006)), it has been and still is often ignored in practice. Though we restrict our attention to models for areal data, the effects of spatial confounding have been studied in other areas such as causal inference (see for example Papadogeorgou et al.(2019)Papadogeorgou, Choirat, and Zigler), or interpolation/prediction (Page et al.(2017)Page, Liu, He, and Sun). Reich et al.(2006)Reich, Hodges, and Zadnik show that adding a conditional autoregressive spatial random effect (CAR) to a fixed effects model can lead to a great change in the posterior mean or a great increase in the posterior variance of the fixed effects, compared to the non-spatial regression model. More precisely, these authors reconsider the Slovenia stomach cancer data, and observe that when the CAR random effect is included in the model, the relationship between stomach cancer and the covariate is diluted. To overcome this problem, they proposed to specify the random effects as orthogonal to the fixed effects. Later, Hodges and Reich (2010) explain why such spatial confounding occurs, and show that adding spatially correlated random effects does not adjust fixed effects for spatially structured missing covariates, as has been generally understood, though they do smooth fitted values.

Following the approach of reparameterizing random effects, Schnell and Bose (2019) also examine the mechanisms of confounding of fixed effects by random effects but they focus on diagnostics to evaluate the effects of confounding rather than on methods to overcome it. Additional work includes Hughes and Haran (2013), Hanks et al.(2015)Hanks, Schliep, Hooten, and Hoeting and Prates et al.(2019)Prates, Assunção, Rodrigues, et al.. Recently, Khan and Calder (2020) observe that restricted regression offers inference similar to models without random effects. Khan and Calder, however, did not discuss model fit or prediction.

The present article deals with spatio-temporal models with covariates and spatial, temporal, and spatio-temporal random effects, pursuing two main goals, namely, correct estimation of linear association between socio-demographic covariates and dowry death (a crime against women very specific to India), and good model predictions. Because we have data and covariates in space and time, confounding may be spatial, temporal or both, and it may mask the association of the outcome with the covariates.

This article proposes different methods to deal with confounding but also with identifiability, as identifying spatial, temporal and spatio-temporal random effects is important for interpretation. On one hand we consider restricted spatial regression (Reich et al.(2006)Reich, Hodges, and Zadnik) and study its application to a spatio-temporal setting. On the other hand, we examine use of constraints to deal with both identifiability and confounding issues. While both methods seem to solve the confounding issue, they lead to notably different results in terms of model fit. Here, we try to disentangle why this happens.

The rest of the article is organized as follows. Section 2 poses the spatial and the spatio-temporal models and briefly revisits identifiability and confounding issues. Section 3 proposes two procedures to alleviate confounding and identify the model: a model reparameterization followed by restricted regression, and the use of orthogonality constraints. In Section 4, both techniques are used to assess the association between dowry death and socio-demographic covariates, and to predict spatio-temporal patterns of risk in Uttar Pradesh, the most populated state in India.

2 Pitfalls in spatial and spatio-temporal areal models

Spatial and spatio-temporal models for areal data have been and still are valuable tools to give a complete picture of the status of a disease, crime or other variable of interest measured using areal counts. Although the benefit and soundness of these models are beyond any doubt, they are not free from inconveniences that should be conveniently addressed. In this article, we revisit spatial and spatio-temporal models with intrinsic conditional autoregressive (ICAR) priors for space and random walks priors for time, and focus on two issues: model identifiability and confounding of fixed effects by random effects. The first usually arises because the spatial and temporal random effects implicitly include an intercept, and the interaction term and the main effects overlap. The second arises from collinearity between fixed and random effects, which may lead to bias and variance inflation of the fixed effects and hence erroneous inference.

2.1 Identifiability and confounding in spatial models

This section focuses on a spatial model for areal count data that includes an intrinsic conditional autoregressive (ICAR) prior for space and highlights the identifiability and confounding issues.

Suppose the area under study (a country, a state) is divided into small areas (counties, districts) denoted by $i=1,\dots S$ , and that $O_i$ stands for the number of observed cases (death or incident cases, number of crimes) in the $i$ th small area. Conditional on the relative risk $r_i$ , a Poisson distribution with mean ${\mu }_i=e_i{r}_i$ is assumed for $O_i$ , where $e_i$ is the number of expected cases computed using, for example, internal standardization. That is

$$O_i|r_i\sim \mathit{Pois}({\mu }_i=e_i{r}_i),\text{log}{\mu }_i=\text{log}{e}_i+\text{log}{r}_i,$$

where $\text{log}(e_i)$ is an offset and $\text{log}{r}_i$ is modelled as

$$\text{log}{r}_i={\beta }_0+{\text{x}}_i^{\prime }\beta +{\xi }_i.$$

(2.1)

Here ${\beta }_0$ is the intercept, ${\text{x}}_i^{\prime }=(x_{i1},\dots ,x_{\mathit{ip}})$ is a $p$ -vector of standardized covariates in the $i$ th area, $\beta =({\beta }_1,\dots ,{\beta }_p{)}^{\prime }$ is the $p$ -vector of fixed effects coefficients, and ${\xi }_i$ is the spatial random effect with an ICAR prior (Besag (1974)). Then, the vector of spatial effects $\xi =({\xi }_1,\dots ,{\xi }_S)\prime$ follows the improper distribution with Gaussian kernel, $p(\xi )\propto \text{exp}\left(-\frac{1}{2{\sigma }_{\xi }^2}\xi \prime {\text{Q}}_{\xi }\xi \right)$ , where ${\text{Q}}_{\xi }$ is the $S\times S$ spatial neighbourhood matrix with $(i,j)$ element ${\text{Q}}_{\xi (\mathit{ij})}=-1$ if areas $i$ and $j$ are neighbours and 0 otherwise, and the $i\mathit{th}$ diagonal element ${\text{Q}}_{\xi (\mathit{ii})}$ is the number of neighbours of the $i\mathrm{th}$ area. In disease mapping studies, typically two regions are neighbours if they share a common border. Because ${\text{Q}}_{\xi }$ ’s rows sum to zero, ${\text{Q}}_{\xi }{1}_S=0$ where $1_S$ is a vector of ones of length $S$ , so an intercept is implicit in the ICAR specification, leading to an identifiability problem with the model intercept. We will assume the spatial map is connected, so ${\text{Q}}_{\xi }$ ’s 0 eigenvalue has multiplicity 1; more general cases are easily accommodated and omitted for simplicity. Goicoa et al.(2018) Goicoa, Adin, Ugarte, and Hodges use this spectral decomposition of the precision matrix of the random effects to reveal the identifiability issue:

$${\text{Q}}_{\xi }={\text{U}}_{\xi }{\Sigma }_{\xi }{\text{U}}_{\xi }^{\prime }=\left[{\text{U}}_{\xi n}:{\text{U}}_{\xi r}\right]\left(\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\tilde{\Sigma }}_{\xi }\hfill \end{array}\right)\left[\begin{array}{c}\hfill {\text{U}}_{\xi n}^{\prime }\hfill \\ \hfill {\text{U}}_{\xi r}^{\prime }\hfill \end{array}\right],$$

where ${\tilde{\Sigma }}_{\xi }$ is a diagonal matrix with the non-null eigenvalues of ${\text{Q}}_{\xi }$ in the main diagonal, and ${\text{U}}_{\xi }=[{\text{U}}_{\xi n}:{\text{U}}_{\xi r}]$ is an orthogonal matrix with columns the eigenvectors of ${\text{Q}}_{\xi }$ . The matrix ${\text{U}}_{\xi }$ is split into the matrix of eigenvectors having null eigenvalues, ${\text{U}}_{\xi n}$ , and the matrix of eigenvectors having non-null eigenvalues, ${\text{U}}_{\xi r}$ . Here the identifiability issue is clearly revealed as ${\text{U}}_{\xi n}$ equals the vector of ones $1_S$ divided by a normalizing constant. Consequently, the spatial Model (2.1) can be expressed in matrix form as

$$\text{log}\text{r}=1_S{\beta }_0+\text{X}\beta +\xi =1_S{\beta }_0+\text{X}\beta +1_S{\beta }_{\xi }+{\text{U}}_{\xi r}{\alpha }_{\xi },$$

(2.2)

where $\text{r}=(r_1,\dots ,r_S{)}^{\prime }$ , $\text{X}=({\text{X}}_1,\dots ,{\text{X}}_p)$ is the $S\times p$ fixed effects design matrix (excluding the intercept) with ${\text{X}}_j=(X_{1j},\dots ,X_{\mathit{Sj}}{)}^{\prime },j=1,\dots ,p$ , ${\beta }_{\xi }=1_S^{\prime }\xi$ , and ${\alpha }_{\xi }={\text{U}}_{\xi r}^{\prime }\xi$ , ${\alpha }_{\xi }\sim \mathcal{N}\left(0,{\sigma }_{\xi }^2{\tilde{\Sigma }}_{\xi }^{-1}\right)$ . Unlike ${\alpha }_{\xi }$ , which has a proper prior, ${\beta }_{\xi }$ has prior precision zero, leading to the identifiability issue: two intercepts are present, the model's and the one implicit in the ICAR. Removing one redundant intercept from the linear predictor, Model (2.2) can be written as

$$\text{log}\text{r}=1_S{\beta }_0+\text{X}\beta +{\text{U}}_{\xi r}{\alpha }_{\xi },$$

resolving the identifiability issue. Alternatively, a sum-to-zero constraint ${\sum }_{i=1}^S{\xi }_i=0$ can be considered.

This reparameterization does not preclude the other potential pitfall of spatial models: spatial confounding. Spatial confounding can be briefly defined as the impossibility of dissociating covariate effects from spatial random effects. As far as we know, Reich et al.(2006)Reich, Hodges, and Zadnik is the first article describing how a CAR random effect can produce changes in the estimates of the fixed effects and inflate the variance compared to a non-spatial model. Later, Hodges and Reich (2010) study the effect of spatial confounding more deeply and show that the usual belief that random effects adjust fixed-effects estimates for missing confounders cannot be sustained. They show that the variance inflation is large if the correlation is large between the covariate ${\text{X}}_j$ and the eigenvector of the spatial matrix ${\text{Q}}_{\xi }$ having the smallest non-null eigenvalue, that is, there is a collinearity problem.

To identify situations where confounding may be a serious issue, these authors hypothesize that the random effects will mask the association between the response and the covariate if the latter exhibits a trend in the long axis of the map. To overcome this confounding, they propose to retain in the model only the part of the random effects lying in the space orthogonal to the fixed effects; for a CAR model with normal response $\text{Y}$ , this is the model

$$\text{Y}=1_S{\beta }_0+\text{X}\beta +\text{L}{\text{L}}^{\prime }\xi ,$$

or its reparameterized version

$$\text{Y}=1_S{\beta }_0+\text{X}\beta +\text{L}{\text{L}}^{\prime }{\text{U}}_{\xi r}{\alpha }_{\xi },$$

where the columns of $\text{L}$ are eigenvectors having non-null eigenvalues (which in fact are all 1) of the projection matrix ${\text{I}}_S-{\text{X}}_{*}({\text{X}}_{*}^{\prime }{\text{X}}_{*}{)}^{-1}{\text{X}}_{*}^{\prime }$ onto the orthogonal space of the fixed effects, ${\text{I}}_S$ is the $S\times S$ identity matrix, and ${\text{X}}_{*}=[1_S:\text{X}]$ . According to Hodges and Reich (2010), this restricted spatial regression takes account of the spatial correlation without changing the estimates of the fixed effects. However, with non-normal responses, for example, a Poisson model, this method requires adjustments. In particular, to deal with collinearity between the fixed and random effects we would use the linear predictor

$$\text{log}\text{r}=1_S{\beta }_0+\text{X}\beta +{\widehat{\text{W}}}^{-1/2}\text{L}{\text{L}}^{\prime }{\widehat{\text{W}}}^{1/2}\xi ,$$

where $\text{W}$ is a diagonal matrix of weights with diagonal elements $W_{\mathit{ii}}=\mathit{Var}(O_i|{\beta }_0,\beta ,\xi )={\mu }_i$ . $\text{W}$ is the weight matrix in the iteratively reweighted least squares algorithm, and to remove collinearity between the fixed and random effects we should delete the combinations of ${\widehat{\text{W}}}^{1/2}\xi$ in the span of ${\widehat{\text{W}}}^{1/2}\text{X}$ (Reich et al.(2006)Reich, Hodges, and Zadnik). Accordingly, $\text{L}$ is now the matrix whose columns are the eigenvectors with non-zero eigenvalues of the orthogonal projection matrix ${\text{I}}_s-{\widehat{\text{W}}}^{1/2}{\text{X}}_{*}({\text{X}}_{*}^{\prime }\widehat{\text{W}}{\text{X}}_{*}{)}^{-1}{\text{X}}_{*}^{\prime }{\widehat{\text{W}}}^{1/2}$ onto the orthogonal space of ${\widehat{\text{W}}}^{1/2}{\text{X}}_{*}$ . It is also possible to reparameterize the model to remove identifiability issues and specify the random effects as orthogonal to the fixed effects:

$$\text{log}\text{r}=1_S{\beta }_0+\text{X}\beta +{\widehat{\text{W}}}^{-1/2}\text{L}{\text{L}}^{\prime }{\widehat{\text{W}}}^{1/2}{\text{U}}_{\xi r}{\alpha }_{\xi },$$

Note that in practice, $\widehat{\text{W}}$ is obtained by fitting the spatial Model (2.1).

2.2 Identifiability and confounding in spatio-temporal models

Now suppose that for each small area or district $i$ , we have data for time periods denoted by $t=1,\dots ,T$ . Similar to the spatial case, and conditional on the spatio-temporal relative risk $r_{\mathit{it}}$ , assume the number of observed cases $O_{\mathit{it}}$ in area $i$ and time $t$ follows a Poisson distribution with mean ${\mu }_{\mathit{it}}=e_{\mathit{it}}{r}_{\mathit{it}}$ where $e_{\mathit{it}}$ is the number of expected cases, that is,

$$O_{\mathit{it}}|r_{\mathit{it}}\sim \mathit{Pois}({\mu }_{\mathit{it}}=e_{\mathit{it}}{r}_{\mathit{it}}),\text{log}{\mu }_{\mathit{it}}=\text{log}{e}_{\mathit{it}}+\text{log}{r}_{\mathit{it}},$$

where the log relative risk is now modelled as

$$\text{log}{r}_{\mathit{it}}={\beta }_0+{\text{x}}_{\mathit{it}}^{\prime }\beta +{\xi }_i+{\gamma }_t+{\delta }_{\mathit{it}}.$$

(2.3)

Here, ${\text{x}}_{\mathit{it}}^{\prime }=(x_{\mathit{it}1},\dots ,x_{\mathit{itp}})$ is a $p$ -vector of standardized spatio-temporal covariates in area $i$ and time $t$ , $\beta =({\beta }_1,\dots ,{\beta }_p{)}^{\prime }$ is the $p$ -vector of fixed effect coefficients, ${\gamma }_t$ is the temporal main effect and ${\delta }_{\mathit{it}}$ is the spatio-temporal interaction term. In matrix form, Model (2.3) is

$$\text{log}\text{r}=1_{\mathit{TS}}{\beta }_0+\text{X}\beta +(1_T\otimes {\text{I}}_S)\xi +({\text{I}}_T\otimes 1_S)\gamma +{\text{I}}_{\mathit{TS}}\delta ,$$

(2.4)

where $\text{r}=(r_{11},\dots ,r_{S1},\dots ,r_{1T},\dots ,r_{\mathit{ST}}{)}^{\prime }$ , $1_{\mathit{TS}}$ and $1_T$ are columns of ones of length $\mathit{TS}$ and $T$ respectively, $\text{X}=({\text{X}}_1,\dots ,{\text{X}}_p)$ is the $\mathit{TS}\times p$ matrix of standardized spatio-temporal covariates with ${\text{X}}_j=(X_{11j},\dots ,X_{S1j},\dots ,X_{1\mathit{Tj}},\dots ,X_{\mathit{STj}}{)}^{\prime }$ , $j=1,\dots ,p$ , and ${\text{I}}_T$ and ${\text{I}}_{\mathit{TS}}$ are $T\times T$ and $\mathit{TS}\times \mathit{TS}$ identity matrices respectively. As before, we consider an ICAR prior for the spatial random effect. For the vector of temporal random effects $\gamma =({\gamma }_1,\dots ,{\gamma }_T{)}^{\prime }$ , we use a first-order random walk (RW1), that is $p(\gamma )\propto \text{exp}\left(-\frac{1}{2{\sigma }_{\gamma }^2}\gamma \prime {\text{Q}}_{\gamma }\gamma \right)$ , where ${\text{Q}}_{\gamma }$ is the RW1’s structure matrix (p. 95 Rue and Held (2005)). The vector of interaction random effects $\delta =({\delta }_{11},\dots ,{\delta }_{S1},\dots ,{\delta }_{1T},\dots ,{\delta }_{\mathit{ST}}{)}^{\prime }$ is assumed to follow a distribution with the Gaussian kernel $p(\delta )\propto \text{exp}\left(-\frac{1}{2{\sigma }_{\delta }^2}\delta \prime {\text{Q}}_{\delta }\delta \right)$ , and ${\text{Q}}_{\mathrm{\delta }}=({\text{Q}}_{\gamma }\otimes {\text{Q}}_{\xi })$ . This interaction term corresponds to the Type IV interaction of Knorr-Held (2000). We consider this type of interaction as it includes spatial and temporal dependence and consequently can produce spatial or temporal confounding. Now consider the spectral decomposition of ${\text{Q}}_{\mathrm{\gamma }}$ :

$${\text{Q}}_{\gamma }={\text{U}}_{\gamma }{\Sigma }_{\gamma }{\text{U}}_{\gamma }^{\prime }=\left[{\text{U}}_{\gamma n}:{\text{U}}_{\gamma r}\right]\left(\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\tilde{\Sigma }}_{\gamma }\hfill \end{array}\right)\left[\begin{array}{c}\hfill {\text{U}}_{\gamma n}^{\prime }\hfill \\ \hfill {\text{U}}_{\gamma r}^{\prime }\hfill \end{array}\right],$$

where ${\text{U}}_{\gamma }=\left[{\text{U}}_{\gamma n}:{\text{U}}_{\gamma r}\right]$ is the $T\times T$ matrix of eigenvectors, ${\text{U}}_{\gamma n}=1_T$ (up to a constant) is the eigenvector with null eigenvalue, ${\text{U}}_{\gamma r}$ is the $T\times (T-1)$ matrix of eigenvectors with non-null eigenvalues, and ${\tilde{\Sigma }}_{\gamma }$ is a diagonal matrix with the non-null eigenvalues in the main diagonal. The spectral decomposition of ${\text{Q}}_{\delta }={\text{Q}}_{\gamma }\otimes {\text{Q}}_{\xi }$ can be expressed as ${\text{Q}}_{\delta }={\text{Q}}_{\gamma }\otimes {\text{Q}}_{\xi }={\text{U}}_{\delta }{\Sigma }_{\delta }{\text{U}}_{\delta }^{\prime }=\left[{\text{U}}_{\delta n}:{\text{U}}_{\delta r}\right]\left(\begin{array}{cc}0& 0\\ 0& {\tilde{\Sigma }}_{\delta }\end{array}\right)\left[\begin{array}{c}{\text{U}}_{\delta n}^{\prime }\\ {\text{U}}_{\delta r}^{\prime }\end{array}\right]$ , where, as in the previous cases, ${\text{U}}_{\delta n}$ is the matrix of eigenvectors having null eigenvalues, ${\text{U}}_{\delta r}$ is the matrix of eigenvectors having non-null eigenvalues, and ${\tilde{\Sigma }}_{\delta }={\tilde{\Sigma }}_{\gamma }\otimes {\tilde{\Sigma }}_{\xi }$ is a diagonal matrix with the non-null eigenvalues in the main diagonal. It can be shown easily that ${\text{U}}_{\delta n}=[{\text{U}}_{\gamma n}\otimes {\text{U}}_{\xi n}:{\text{U}}_{\gamma n}\otimes {\text{U}}_{\xi r}:{\text{U}}_{\gamma r}\otimes {\text{U}}_{\xi n}],{\text{U}}_{\delta r}=[{\text{U}}_{\gamma r}\otimes {\text{U}}_{\xi r}]$ .

Similar to the spatial case, the spatio-temporal Model (2.4) can be expressed in matrix form as

$$\begin{array}{ccc}\hfill \text{log}\text{r}& \hfill =\hfill & 1_{\mathit{TS}}{\beta }_0+\text{X}\beta +1_{\mathit{TS}}{\beta }_{\xi }+(1_T\otimes {\text{U}}_{\xi r}){\alpha }_{\xi }+1_{\mathit{TS}}{\beta }_{\gamma }+({\text{U}}_{\gamma r}\otimes 1_S){\alpha }_{\gamma }\hfill \\ \hfill & \hfill +\hfill & [1_{\mathit{TS}}:1_T\otimes {\text{U}}_{\xi r}:{\text{U}}_{\gamma r}\otimes 1_S]{\beta }_{\delta }+({\text{U}}_{\gamma r}\otimes {\text{U}}_{\xi r}){\alpha }_{\delta },\hfill \end{array}$$

(2.5)

where ${\beta }_{\gamma }={\text{U}}_{\gamma n}^{\prime }\gamma =1_T^{\prime }\gamma$ , ${\alpha }_{\gamma }={\text{U}}_{\gamma r}^{\prime }\gamma \sim \mathcal{N}\left(0,{\sigma }_{\gamma }^2{\tilde{\Sigma }}_{\gamma }^{-1}\right)$ , ${\beta }_{\delta }={\text{U}}_{\delta n}^{\prime }\delta$ , and ${\alpha }_{\delta }={\text{U}}_{\delta r}^{\prime }\delta \sim \mathcal{N}\left(0,{\sigma }_{\delta }^2{\tilde{\Sigma }}_{\delta }^{-1}\right)$ . The reparameterized form of Model (2.5) sheds light on the identifiability issues in spatio-temporal models as it lays bare repeated terms. Removing those superfluous terms gives the following model:

$$\text{log}\text{r}=1_{\mathit{TS}}{\beta }_0+\text{X}\beta +(1_T\otimes {\text{U}}_{\xi r}){\alpha }_{\xi }+({\text{U}}_{\gamma r}\otimes 1_S){\alpha }_{\gamma }+({\text{U}}_{\gamma r}\otimes {\text{U}}_{\xi r}){\alpha }_{\delta },$$

(2.6)

which overcomes identifiability issues. For more details about this reparameterization, a generalization to RW2 priors for time, and the other interaction types described by Knorr-Held (2000), the reader is referred to Goicoa et al.(2018)Goicoa, Adin, Ugarte, and Hodges, where sum-to-zero constraints are alternatively derived to achieve model identifiability.

Confounding issues are more challenging in spatio-temporal settings than in spatial settings. In a spatio-temporal model, the covariates can exhibit spatial patterns each year, or temporal patterns in each area. As the model includes both spatial and temporal random effects along with interaction terms, the source of confounding can be spatial, temporal, or both. Note that the reparameterized Model (2.6) may present confounding problems as the covariates may be collinear with the design matrix of the spatial, temporal, or spatio-temporal random effects. The following section proposes two methods to alleviate confounding in spatio-temporal models, restricted spatial regression and constraints that make the estimated random effects orthogonal to the fixed effects.

3 Alleviating confounding in spatio-temporal models

Constraints can be used to make the random effects orthogonal to the fixed effects and, thus, alleviate confounding. The idea of inducing orthogonality between the fixed and random effects is similar to restricted regression, but they have some differences that can lead to notably distinct results. This section shows that constraining the random effects to be orthogonal to the fixed effects is not equivalent to removing from the linear predictor the component of the random effects in the span of the fixed effects, as might reasonably be assumed. In some cases, these differences are hardly noticeable in the spatial case but they become important in spatio-temporal settings.

3.1 Model reparameterization and restricted regression

Consider Model (2.6). This reparameterized spatio-temporal model is convenient as it solves the identifiability problems by removing repeated terms in the spatial and temporal main effects and the interaction random effects. The confounding issues are more challenging now because covariates are spatio-temporal, so we can have collinearity between the covariates and the spatial term, the covariates and the temporal term, or the covariates and the spatio-temporal interaction. Another concern is how to assess these correlations as the covariates have $T\times S$ entries while the spatial and temporal random effects have $S$ and $T$ elements. This section considers the matrix $\text{L}$ with columns that are the eigenvectors having non-null eigenvalues (all equal 1) of the projection matrix $\text{P}c={\text{I}}_{\mathit{TS}}-\widehat{\text{W}}1/2{\text{X}}_{*}({\text{X}}_{*}^{\prime }\widehat{\text{W}}{\text{X}}_{*})-1{\text{X}}_{*}^{\prime }\widehat{\text{W}}1/2$ , where now ${\text{X}}_{*}=[1_{\mathit{TS}}:\text{X}]$ , and $\text{X}$ is the design matrix of covariates. The matrix $\text{P}c$ projects onto the space orthogonal to the (scaled) fixed effects $\widehat{\text{W}}1/2{\text{X}}_{*}$ . Consider also the matrix $\text{K}$ , with columns that are the eigenvectors having eigenvalue 1 of the projection matrix $\text{P}={\widehat{\text{W}}}^{1/2}{\text{X}}_{*}({\text{X}}_{*}^{\prime }\widehat{\text{W}}{\text{X}}_{*}{)}^{-1}{\text{X}}_{*}^{\prime }\widehat{\text{W}}1/2$ . Note that $\text{K}=\widehat{\text{W}}1/2{\text{X}}_{*}$ and that $(\text{K}\text{K}\prime +\text{L}\text{L}\prime )={\text{I}}_{\mathit{TS}}$ ; we therefore propose the following spatio-temporal model:

$$\begin{array}{rr}\hfill \text{log}r=& \hfill 1_{TS}{\beta }_0+X\beta +{\widehat{W}}^{{}^{-1/2}}\left(K{K}^{\prime }+L{L}^{\prime }\right){\widehat{W}}^{{}^{1/2}}\left(1_T\otimes U_{{\xi }_r}\right){\alpha }_{\xi }\\ \hfill & \hfill +{\widehat{W}}^{{}^{-1/2}}\left(K{K}^{\prime }+L{L}^{\prime }\right){\widehat{W}}^{{}^{1/2}}\left(U_{{\gamma }_r}\otimes 1_S\right){\alpha }_{\gamma }\\ \hfill & \hfill +{\widehat{W}}^{{}^{-1/2}}\left(K{K}^{\prime }+L{L}^{\prime }\right){\widehat{W}}^{{}^{1/2}}\left(U_{{\gamma }_r}\otimes U_{{\xi }_r}\right){\alpha }_{\delta }\end{array}$$

As the terms involving the matrix $\text{K}$ are in the span of the fixed effects ${\text{X}}_{*}$ , they are removed and the model becomes

$$\begin{array}{c}\log r=1_{TS}{\beta }_0+X\beta +{\widehat{W}}^{-1/2}L{L}^{\prime }{\widehat{W}}^{1/2}\left(1_T\otimes U_{{\xi }_r}\right){\alpha }_{\xi }+{\widehat{W}}^{-1/2}L{L}^{\prime }{\widehat{W}}^{1/2}\left(U_{{\gamma }_r}\otimes 1_S\right){\alpha }_{\gamma }\\ +{\widehat{W}}^{-1/2}L{L}^{\prime }{\widehat{W}}^{1/2}\left(U_{{\gamma }_r}\otimes U_{{\xi }_r}\right){\alpha }_{\delta }\end{array}$$

(3.1)

Model (3.1) deserves comment. First, by using the matrix $\text{L}$ , all random effects have been restricted to be orthogonal to the fixed effects ( $\text{L}$ is orthogonal to ${\widehat{W}}^{1/2}{X}_{*}$ by construction, that is ${X^{\prime }}_{*}{\widehat{W}}^{1/2}L=0$ ), but we could also restrict only some of the random effects. Restricting all the random effects may not be necessary if only some of the random effects confound a fixed effect. Thus, it is possible to orthogonalize only the spatial, or temporal, or interaction random effects. The other important issue is the final estimation of the spatial and temporal random effects. Compared to Model (2.6), the spatial and temporal random effects in Model (3.1) undergo a substantial change. In Model (2.6), the spatial and temporal main effects are $(1_T\otimes {\text{U}}_{\xi r}){\alpha }_{\xi }$ and $({\text{U}}_{\gamma r}\otimes 1_S){\alpha }_{\gamma }$ respectively; clearly $S$ spatial effects are repeated in time, while $T$ temporal effects are repeated for the $S$ small areas. In the orthogonalized Model (3.1), however, the spatial and temporal main effects are ${\widehat{W}}^{-1/2}L{L}^{\prime }{\widehat{W}}^{1/2}(1_T\otimes U_{\xi r}){\alpha }_{\xi }$ and ${\widehat{W}}^{-1/2}L{L}^{\prime }{\widehat{W}}^{1/2}(U_{\gamma r}\otimes 1_S){\alpha }_{\gamma }$ respectively. Because these terms include the matrix $\text{L}$ , which has $T\times S$ rows (the covariates are spatio-temporal, in general), the spatial effect associated with the $i\mathrm{th}$ area is different in each time period as it depends on the value of the covariates in that period. Similarly, the effect of time period $t$ is also different for each area as it depends on the covariates in that area. Consequently, Model (3.1) has time-varying spatial effects and space-varying temporal effects.

3.2 Constraints to alleviate confounding

Placing constraints is a way to achieve identifiability in spatio-temporal disease mapping models (Goicoa et al.(2018)Goicoa, Adin, Ugarte, and Hodges). Constraints can be used not just to identify models but also to alleviate spatial confounding, by constraining the random effects to be orthogonal to the fixed effects. This section considers such constraints. Specifically, consider Model (2.4) and consider the linear predictor

$${\widehat{W}}^{1/2}{1}_{\text{TS}}{\beta }_0+{\widehat{W}}^{1/2}X\beta +{\widehat{W}}^{1/2}\left(1_T\otimes I_S\right)\xi +{\widehat{W}}^{1/2}\left(I_T\otimes 1_S\right)\gamma +{\widehat{W}}^{1/2}{I}_{TS}\delta .$$

To make the random effects orthogonal to the fixed effects, these constraints are required:

$$\begin{array}{l}{\left[1_{TS}:X\right]}^{\prime }\widehat{W}\text{(}{1}_T\otimes I_s\text{)}\xi =0\iff {\displaystyle \sum _{i=1}^S{\xi }_i{\widehat{w}}_{i.}=0,}{\displaystyle \sum _{i=1}^S{\xi }_i(\widehat{w}{x}_j)i.=0,\forall j=1,\mathrm{...},p,}\\ {\left[1_{TS}:X\right]}^{\prime }\widehat{W}\text{(}{1}_T\otimes I_s\text{)}\gamma =0\iff {\displaystyle \sum _{i=1}^S{\gamma }_i\widehat{w}{.}_i=0,}{\displaystyle \sum _{t=1}^T{\gamma }_t(\widehat{w}{x}_j).t=0,\forall j=1,\mathrm{...},p,}\\ {\left[1_{TS}:X\right]}^{\prime }\widehat{W}\delta =0\iff {\displaystyle \sum _{i=1}^S{\displaystyle \sum _{i=1}^T{\widehat{w}}_{it}{\delta }_i{}_t=0,}}{\displaystyle \sum _{i=1}^S{\displaystyle \sum _{t=1}^T{\widehat{w}}_{it}{x}_{it}{\delta }_{it}=0,\forall j=1,\mathrm{...},p,}}\end{array}$$

(3.2)

where ${\widehat{w}}_{i·}={\sum }_{t=1}^T{\widehat{w}}_{\mathit{it}}$ , ${\widehat{w}}_{·t}={\sum }_{i=1}^S{\widehat{w}}_{\mathit{it}}$ , $(\widehat{w}{x}_j{)}_{i·}={\sum }_{t=1}^T{\widehat{w}}_{\mathit{it}}{x}_{\mathit{jit}}$ , and $(\widehat{w}{x}_j{)}_{·t}={\sum }_{i=1}^S{\widehat{w}}_{\mathit{it}}{x}_{\mathit{jit}}$ . That is, the spatial random effect is constrained to be orthogonal to the time-averaged covariates at each location, the temporal random effect is constrained to be orthogonal to the space-averaged covariates at each time, and the interaction random effect is constrained to be orthogonal to the full fixed-effects design matrix. However, the interaction term $\delta$ is confounded with the spatial and temporal main (random) effects and these additional constraints are required:

$$[(1_T\otimes {\text{I}}_S):({\text{I}}_T\otimes 1_S)]\prime \widehat{\text{W}}\delta =0\iff \begin{array}{c}{\displaystyle \sum }_{i=1}^S{\hat{w}}_{it}{\delta }_{it}=0\forall t=1,\dots ,T.\\ {\displaystyle \sum }_{t=1}^T{\hat{w}}_{it}{\delta }_{it}=0\forall i=1,\dots ,S.\end{array}$$

(3.3)

making the constraint ${\sum }_{i=1}^S{\sum }_{t=1}^T{\widehat{w}}_{\mathit{it}}{\delta }_{\mathit{it}}=0$ redundant. Avoiding confounding between the interaction and the spatial and temporal main effect terms is crucial for model interpretation. The spatial main effects capture spatial variation that is not accounted for by the covariates so they may help identify spatial risk factors that have not been included in the model. Similarly, the temporal main effects may help identify risk factors associated with time. The interaction random effects usually explain a small portion of variability and capture deviations from the main effects. Given that the interaction term's null space equals the design matrix of the spatial and temporal main effects (see Model (2.5)), confounding between main effect and interaction random effects may make it impossible to achieve one of the goals of disease mapping: smoothing in space and time to uncover spatial and temporal patterns.

Though the restricted regression and the constraints approaches would seem to be equivalent, they are in fact different. The restricted regression approach focuses on the reduced Model (3.1), where the part of the random effects in the span of the fixed effects has been removed, and orthogonality is achieved through the matrix $\text{L}$ . The constraints approach, by contrast, starts with the full Model (2.4) and Equation (3.2)’s constraints force the random effects to be orthogonal to the fixed effects, which is a way to remove collinearities between them. A key distinction between the restricted regression and constraints approaches is that in the former, the spatial and temporal effects change in time and space respectively due to the spatio-temporal nature of the matrix $\text{L}$ , while in the constraints approach, the spatial effects remain constant in time and the temporal effects do not change in space.

One may think the methods are equivalent because removing the part of the random effects in the span of the fixed effects means $\text{K}\prime {\widehat{\text{W}}}^{1/2}(1_{\text{T}}\otimes {\text{I}}_S)\xi =0$ , $\text{K}\prime {\widehat{\text{W}}}^{1/2}({\text{I}}_T\otimes 1_S)\gamma =0$ , and $\text{K}\prime {\widehat{\text{W}}}^{1/2}{\text{I}}_{TS}\delta =0$ , which match the constraints in Equation (3.2), taking into account that $\text{K}={\widehat{\text{W}}}^{1/2}[1:\text{X}]$ . However, placing constraints is in fact equivalent to oblique projections of the random effects. That is, the constraints in Equations (3.2) and (3.3) correspond to this model:

$$\begin{array}{ccc}\hfill \text{log}\text{r}& \hfill =\hfill & 1_{\mathit{TS}}{\beta }_0+\text{X}\beta +(1_T\otimes {\text{I}}_S){\text{P}}_{\xi }\xi +({\text{I}}_T\otimes 1_S){\text{P}}_{\gamma }\gamma +{\text{P}}_{\delta }\delta \hfill \\ \hfill & \hfill =\hfill & 1_{\mathit{TS}}{\beta }_0+\text{X}\beta +(1_T\otimes {\text{P}}_{\xi })\xi +({\text{P}}_{\gamma }\otimes 1_S)\gamma +{\text{P}}_{\delta }\delta ,\hfill \end{array}$$

where

$${\text{P}}_{\xi }={\text{L}}_{\xi }{[{\text{L}}_{\xi }{}^{\prime }{\text{Q}}_{\xi }{\text{L}}_{\xi }]}^{-1}{\text{L}}_{\xi }{}^{\prime }{\text{Q}}_{\xi },{\text{P}}_{\gamma }={\text{L}}_{\gamma }{[{\text{L}}_{\gamma }{}^{\prime }{\text{Q}}_{\gamma }{\text{L}}_{\gamma }]}^{-1}{\text{L}}_{\gamma }{}^{\prime }{\text{Q}}_{\gamma },{\text{P}}_{\delta }={\text{L}}_{\delta }{[{\text{L}}_{\delta }{}^{\prime }{\text{Q}}_{\delta }{\text{L}}_{\delta }]}^{-1}{\text{L}}_{\delta }{}^{\prime }{\text{Q}}_{\delta }$$

(3.4)

are made to be, respectively, oblique projections onto the orthogonal complements of the row spaces of the matrices

$${\text{B}}_{\xi }={\text{X}}_{*}^{\text{'}}\widehat{\text{W}}\left(1_T\otimes {\text{I}}_{\text{S}}\right),\hspace{1em}\text{B}\gamma ={\text{X}}_{*}^{\text{'}}\widehat{\text{W}}\left({\text{I}}_T\otimes {\text{1}}_S\right),\hspace{1em}{\text{B}}_{\delta }={\left[\left(1_T\otimes {\text{I}}_{\text{S}}\right):\left({\text{I}}_T\otimes {\text{1}}_S\right):\text{X}\right]}^{\text{'}}\widehat{\text{W}}$$

(3.5)

by setting L_ξ to be the matrix whose columns are the eigenvectors with non-zero eigenvalues of the projection matrix ${\text{I}}_S-{\text{B}}_{\xi }{}^{\prime }{({\text{B}}_{\xi }{\text{B}}_{\xi }{}^{\prime })}^{-1}{\text{B}}_{\xi }$ , and similarly for ${\text{L}}_{\gamma }$ and ${\text{L}}_{\delta }$ . See Supplementary Material A for the equivalence of the oblique projections and the constraints and detailed expressions of matrices ${\text{B}}_{\xi }$ , ${\text{B}}_{\gamma }$ , and ${\text{B}}_{\delta }$ .

3.3 Model fitting and inference

Two main methods have been used to fit spatial and spatio-temporal disease mapping models: a fully Bayesian approach and an empirical Bayes approach. The latter provides point estimates of quantities of interest, traditionally using penalized quasi-likelihood (PQL; Breslow and Clayton (1993)). It has proven to be interesting because it is relatively simple and has few convergence problems, and it has been used to fit different models such as the ICAR or P-splines (Dean et al.(2001)Dean, Ugarte, and Militino; Dean et al.(2004)Dean, Ugarte, and Militino; Ugarte et al.(2010)Ugarte, Goicoa, and Militino; Ugarte et al.(2012)Ugarte, Goicoa, Etxeberria, and Militino). However, PQL automatically places sum-to-zero constraints due to the rank deficiency of the random effects covariance matrices, and placing additional constraints is not so straightforward (see the Appendix for full details).

The fully Bayesian approach is probably the most-used technique for model fitting and inference because it provides a full posterior distribution for quantities of interest. Though it has traditionally relied on Markov chain Monte Carlo (MCMC), the computational burden of MCMC prompted development of other attractive procedures, including INLA (integrated nested Laplace approximations; Rue et al.(2009)Rue, Martino, and Chopin). INLA's main advantage is that it provides approximate Bayesian inference without using MCMC, leading to substantial reduction in computational cost. INLA is ready to use in the free software R using the package R-INLA, which has implemented general models that can be adapted to disease mapping. Also, imposing constraints in INLA is relatively simple.

4 Data analysis: The association between dowry deaths and socio-demographic covariates in Uttar Pradesh, India

In this section, the two approaches to alleviate confounding are used to assess the potential association between dowry deaths and some socio-demographic covariates in Uttar Pradesh, the most populated state in India. Dowry deaths is a cruel form of violence against women deep-rooted in India. It is strongly related to dowry, which can be defined as the amount of money, property, or goods that the bride's family gives to the groom or his relatives for the marriage. Although dowry was originally designed to protect women from unfair traditions, such as the impossibility of women owning immovable property (see Banerjee (2014)), it has become a means by which the husband or husband's relatives extort higher dowries under the threat of physical violence against the wife. This form of violence can be extended over time, ending in what is known as a dowry death. If a woman commits suicide because she has experienced mental or physical violence related to the dowry, this is also considered a dowry death.

We have data on dowry deaths in 70 districts of Uttar Pradesh, the Indian state with the highest rate of dowry deaths, during the period 2001–2014. In 2014, the last year of the study period, 8 455 dowry deaths were registered representing 29.2% of all dowry deaths in India. One of the difficulties of combatting this crime against women is the lack of knowledge about potential risk factors that might be associated with dowry deaths, and thus help to predict this crime. One hypothesized risk factor is the sex ratio, that is, the number of females per 1 000 males. The literature has contradictory results about the sex ratio: some authors find a negative association between dowry deaths and sex ratio (Mukherjee et al.(2001)Mukherjee, Rustagi, and Krishnaji), while others find a positive association (Dang et al.(2018)Dang, Kulkarni, and Gaiha). In a more in-depth study of dowry deaths in Uttar Pradesh, Vicente et al.(2020)Vicente, Goicoa, Fernandez-Rasines, and Ugarte consider spatio-temporal models and include some potential risk factors as covariates to assess their association with dowry deaths. However, these authors do not address confounding issues. We now focus on some of those covariates, namely sex ratio ( $x_1)$ , population density ( $x_2)$ , female literacy rate ( $x_3)$ , per capita income ( $x_4)$ , murder rate ( $x_5)$ , and burglary rate ( $x_6)$ , estimating their association with dowry deaths accounting for confounding. The goal is to see the effect of confounding on the estimates and standard errors of the fixed effects. We also compare restricted regression and constraints in terms of model fit and complexity using DIC (Deviance Information Criterion, Spiegelhalter et al.(2002)Spiegelhalter, Best, Carlin, and Van Der Linde) and AIC (Akaike Information Criterion, Akaike (1974)) in a Bayesian and frequentist approach respectively.

4.1 Data analysis

For purely spatial models, Reich et al.(2006)Reich, Hodges, and Zadnik and Hodges and Reich (2010) argued that spatial confounding is created by a high correlation between a covariate and the eigenvector of the spatial precision matrix having the smallest non-null eigenvalue. In the spatio-temporal setting, this suggests examining analogous correlations for the spatial, temporal, and spatio-temporal random effects. Therefore, before fitting any models, we examine the data and compute some correlations. Specifically, we split the covariates into spatial vectors for each year of the period, and computed Pearson's correlation between those spatial vectors and ${\text{U}}_{{\xi }_{69}}$ , the eigenvector of the spatial precision matrix with the smallest non-null eigenvalue, so for each covariate, we compute fourteen spatial correlations. Similarly, we split the covariates into temporal vectors for each district, and evaluate the correlations between those temporal vectors and ${\text{U}}_{{\gamma }_{13}}$ , the eigenvector of the temporal precision matrix with the smallest non-null eigenvalue, seventy temporal correlations for each covariate. Boxplots of correlations between the covariates and the spatial eigenvector ${\text{U}}_{{\xi }_{69}}$ for each year (left) and correlations between the covariates and the temporal eigenvector ${\text{U}}_{{\gamma }_{13}}$ for each area (right) are displayed in Figure B.1 in Supplementary Material B, together with some explanation.

Now we consider four spatio-temporal models.

Model ST1: Simple spatio-temporal Poisson model (Equation (2.4) without random effects)

Model ST2: Spatio-temporal model without accounting for confounding (Equations (2.4) and (2.6) for INLA and PQL respectively).

Model ST3: Spatio-temporal model with restricted regression (Equation (3.1)).

Model ST4: Spatio-temporal model with orthogonality constraints (Equations (3.2) and (3.3))

Within the full Bayesian approach, we have used Normal prior distributions with mean 0 and variance equal to 1 000 for fixed effects. Regarding random effects, we have used uniform prior distributions on the positive real line for the standard deviations. We have also considered logGamma(1,0.00005) priors for the log-precisions and PC priors for the standard deviations (see, e.g., Simpson et al.(2017)Simpson, Rue, Riebler, Martins, and Sørbye), but these gave similar results and hence they are not shown here to save space.

Table 1:

Posterior mean, posterior standard deviation, and 95% credible intervals of the fixed effects for models fitted with INLA (left), and point estimates, standard errors and 95% confidence intervals obtained with PQL (right). The results for Model ST2 are in bold when the covariate is not significantly associated with dowry death

		INLA (Simplified Laplace)				PQL (tol=1e-5)
		Sex ratio
	Model	Mean	SD	$q_{0.025}$	$q_{0.975}$	Estimate	SE	$q_{0.025}$	$q_{0.975}$
${\beta }_1$	ST1	$-$ 0.2366	0.0085	$-$ 0.2532	$-$ 0.2200	$-$ 0.2366	0.0085	$-$ 0.2532	$-$ 0.2200
	ST2	$-$ 0.0920	0.0470	$-$ 0.1836	0.0015	$-$ 0.0924	0.0459	$-$ 0.1825	$-$ 0.0024
	ST3	$-$ 0.2303	0.0085	$-$ 0.2470	$-$ 0.2135	$-$ 0.2299	0.0082	$-$ 0.2459	$-$ 0.2139
	ST4	$-$ 0.2284	0.0082	$-$ 0.2445	$-$ 0.2125	$-$ 0.2281	0.0081	$-$ 0.2441	$-$ 0.2121
		Population Density
${\beta }_2$	ST1	$-$ 0.0917	0.0065	$-$ 0.1044	$-$ 0.0791	$-$ 0.0917	0.0065	$-$ 0.1044	$-$ 0.0792
	ST2	$-$ 0.0069	0.0318	$-$ 0.0682	0.0568	$-$ 0.0081	0.0304	$-$ 0.0677	0.0515
	ST3	$-$ 0.0901	0.0066	$-$ 0.1031	$-$ 0.0773	$-$ 0.0904	0.0061	$-$ 0.1024	$-$ 0.0785
	ST4	$-$ 0.0946	0.0063	$-$ 0.1071	$-$ 0.0824	$-$ 0.0949	0.0063	$-$ 0.1073	$-$ 0.0826
		Female literacy rate
${\beta }_3$	ST1	0.0992	0.0076	0.0843	0.1141	0.0992	0.0076	0.0843	0.1141
	ST2	$-$ 0.0478	0.0501	$-$ 0.1482	0.0487	$-$ 0.0469	0.0474	$-$ 0.1398	0.0460
	ST3	0.0946	0.0080	0.0789	0.1104	0.0946	0.0077	0.0796	0.1096
	ST4	0.0975	0.0077	0.0823	0.1126	0.0975	0.0077	0.0823	0.1126
		Per capita income
${\beta }_4$	ST1	$-$ 0.0661	0.0084	$-$ 0.0827	$-$ 0.0498	$-$ 0.0661	0.0084	$-$ 0.0827	$-$ 0.0499
	ST2	$-$ 0.0196	0.0296	$-$ 0.0776	0.0385	$-$ 0.0198	0.0288	$-$ 0.0763	0.0368
	ST3	$-$ 0.0651	0.0085	$-$ 0.0819	$-$ 0.0486	$-$ 0.0650	0.0081	$-$ 0.0809	$-$ 0.0491
	ST4	$-$ 0.0680	0.0080	$-$ 0.0837	$-$ 0.0525	$-$ 0.0678	0.0079	$-$ 0.0833	$-$ 0.0522
		Murder rate
${\beta }_5$	ST1	0.0833	0.0076	0.0682	0.0982	0.0833	0.0076	0.0682	0.0982
	ST2	0.0846	0.0203	0.0446	0.1244	0.0845	0.0202	0.0449	0.1241
	ST3	0.0906	0.0081	0.0748	0.1064	0.0907	0.0077	0.0756	0.1057
	ST4	0.0881	0.0079	0.0726	0.1034	0.0881	0.0079	0.0727	0.1035
		Burglary rate
${\beta }_6$	ST1	0.0419	0.0062	0.0297	0.0541	0.0419	0.0062	0.0297	0.0540
	ST2	0.0535	0.0161	0.0220	0.0850	0.0535	0.0158	0.0226	0.0844
	ST3	0.0424	0.0068	0.0291	0.0557	0.0423	0.0063	0.0300	0.0547
	ST4	0.0431	0.0063	0.0307	0.0554	0.0429	0.0063	0.0306	0.0552

Table 1 shows the posterior mean and standard deviation of the fixed effects with a 95% credible interval, computed using INLA (simplified Laplace strategy). Point estimates obtained with PQL are also displayed with their standard error and 95% confidence interval. For sex ratio, the estimate with Model ST2 and INLA is about 40% (in absolute value) of the estimates obtained with Models ST1, ST3, and ST4, and the posterior SD is 5.5 times higher. For Model ST2, the 95% credible interval includes 0, while the intervals from the other models are far from zero. The results with PQL are similar although Model ST2’s confidence interval for sex ratio barely excludes 0. For population density, the estimate from Model ST2 is less than 10% of the estimates obtained with Models ST1, ST3, and ST4, with INLA or PQL, and again the posterior SD with Model ST2 is about 5 times higher than with the other models. The consequence is that using Model ST2, the association between population density and dowry deaths is not significant. Per capita income has similar results. The effect of confounding on the estimated association with female literacy rate is also noteworthy: with Model ST2, the estimated effect is negative (though not significant) whereas with the rest of models is positive and significant. The estimated associations with murder rate and burglary rate are similar for the four models, though the posterior SD is clearly larger for Model ST2. These results are revealing and illustrate the potential harmful consequences of ignoring the effects of confounding: the estimated association between the response and the covariate may be diluted or dramatically changed. Posterior estimates of the hyperparameters (standard deviations) are displayed in Table B.1 in Supplementary Material B. Similar estimates for the standard deviations are obtained with Models ST2, ST3, and ST4.

Table 2:

Model selection criteria and computing time (seconds) for spatio-temporal models fit with INLA and PQL. Computations were made on a twin superserver with four processors, Intel Xeon 6C and 96GB RAM, using the R-INLA (stable) version 19.09.03.

	INLA (simplified Laplace)				PQL (tol=1e-5)
	$\bar{D}$	$p_D$	DIC	Time	Deviance	$\mathit{Df}$	AIC	Time
Model ST1	8471.72	7.18	8478.90	3	8466.20	7.00	8480.20	1
Model ST2	5962.60	239.26	6201.86	18	5727.89	236.48	6200.85	74
Model ST3	5962.59	239.27	6201.86	192	5727.85	236.46	6200.77	86
Model ST4	6492.87	232.01	6724.88	22	6269.53	227.64	6724.81	180

Spatio-temporal models accounting for confounding (Models ST3 and ST4) lead to similar estimates of the fixed effects and their posterior standard deviations but they differ in terms of model selection criteria. Table 2 displays the mean deviance ( $\bar{D}$ ), the effective number of parameters $p_D$ , DIC and computing time for the INLA fits. For the PQL fit, the deviance, the number of parameters ( $\mathit{Df}$ ) and AIC are provided. As expected, the model without random effects has the worst fit; the six covariates are not enough to explain the data's variation. For Models ST3 and ST4, the difference in fit is remarkable. Clearly, Model ST3 provides a much better fit: the differences in $\bar{D}$ and DIC are about 500 points. The INLA fit with a simplified Laplace strategy and using constraints (Model ST4) is much faster than restricted regression (22 and 192 seconds respectively). Computing time for Model ST3 in INLA has been reduced (about one half at least) by plugging in the posterior modes of the hyperparameters obtained from Model ST2 as initial values. Note that the total time required to fit Model ST3 should include the computing time of Model ST2 (18 seconds in our data analysis). Similarly, computing times for Models ST3 and ST4 in PQL have been reduced about 23% using the variance components estimates obtained in Model ST2 as initial values in the estimation algorithm. In this case, the computing times shown in Table 2 correspond to the total time needed to fit the corresponding models. Additionally, Supplementary Material B provides scatter plots of the estimated relative risks from Models ST2, ST3, and ST4 fitted with INLA and PQL (Figure B.2); posterior spatial patterns (top row), posterior temporal patterns (middle row) obtained from Models ST2, ST3, and ST4 fitted with INLA, and posterior spatio-temporal patterns (see Adin et al.(2017)Adin, Martí nez-Beneito, Botella-Rocamora, Goicoa, and Ugarte) for three districts, Agra, Balrampur, and Gautam Buddha Nagar (bottom row) in Figure B.3. Finally, INLA relative risk estimates (posterior means) obtained with models ST3 and ST4 are shown in Figure B.4 for the same three districts.

Why do Models ST3 and ST4 fit so differently? First, the spatial and temporal terms in Models ST4 are $(1_T\otimes {\text{I}}_S)\xi$ and $({\text{I}}_T\otimes 1_S)\gamma$ respectively. Consequently the spatial effects are repeated in every year and the temporal effects are repeated in each area. In Model ST3, however, the spatial and temporal terms $(1_T\otimes {\text{U}}_{\xi r})\xi$ and $({\text{U}}_{\gamma r}\otimes 1_S)\gamma$ are premultiplied by ${\widehat{\text{W}}}^{-1/2}\text{LL}{\widehat{\text{W}}}^{-1/2}$ . Because the matrix $\text{L}$ contains spatio-temporal information — it depends on the spatio-temporal covariate $\text{X}$ — the spatial effect for the $i$ th area is different in each time period and the temporal effect in year $t$ is also different for each area. Second, no substantial change in relation to Model ST2 is made in Model ST3. Basically, the random effects in this latter model are split into two pieces, one in the span of the covariates and one orthogonal to the covariates. When Model ST3 removes the part in the span of the covariates, it is simply discarding redundant information. By contrast, Model ST4 changes the model by forcing the random effects to be orthogonal to the fixed effects. Moreover, the spatial random effects lie in the space orthogonal to the time-weighted-added covariates, and the temporal random effects are orthogonal to the spatial-weighted-added covariates. This is equivalent to an oblique projection onto the orthogonal subspace of the fixed effects, unlike Model ST3, where the projection is orthogonal. Because the orthogonal projection minimizes the distance between the original random effects and the projection, this could explain the improvements in fit over the oblique projection (constraints).

Finally, note that we also fit Model ST3 without premultiplying the temporal and spatio-temporal effects by ${\widehat{\text{W}}}^{-1/2}{\text{LL}}^{\text{'}}{\widehat{\text{W}}}^{-1/2}$ and the results are nearly identical, indicating that all confounding arises from the spatial term. We also fit Model ST3 premultiplying by ${\widehat{\text{W}}}^{-1/2}{\text{LL}}^{\text{'}}{\widehat{\text{W}}}^{-1/2}$ only to the temporal and the spatio-temporal effects, and the confounding effects are not avoided.

5 Discussion

Including spatially correlated random effects in a model can seriously affect inference about fixed effects due to confounding. This is particularly dramatic in ecological spatial regression where the main objective is to estimate associations between the response variable and certain covariates. These relationships can be masked due to bias and variance inflation of the fixed effects caused by confounding. Though documented in the literature, spatial confounding has generally been ignored in applications and this practice has carried over to spatio-temporal settings. Here we study confounding in spatio-temporal ecological models in which including temporally correlated random effects and space-time interaction random effects (spatially and temporally correlated) can exacerbate confounding problems. We have considered two procedures to remedy the potentially harmful effects of confounding, restricted regression and constraints.

In light of this article's results, we would like to emphasize some points and provide some guidelines to practitioners. First, the relative risk estimates are not affected by confounding, so if the relative risks are of primary interest, ignoring confounding is not a problem. Second, both restricted regression and orthogonal constraints alleviate confounding and provide rather similar estimates of the fixed effects and their standard errors. However, the two approaches differ importantly in terms of model selection criteria and computing time. Although the constraints approach is computationally more efficient than restricted regression (in INLA), the latter gives clearly better fits. Consequently, if the target of the analysis is to establish associations between risk factors and the phenomenon under study, along with studying spatio-temporal patterns of risk, we recommend using restricted regression. We have also observed that differences between the approaches become more evident as the number of covariates increases. The difference in fits may arise because the deviance, and hence DIC (and AIC), are directly related to the orthogonal projection, the projection used in restricted regression. Because the orthogonal projection minimizes least squares, and $\bar{D}$ (mean deviance) is a transform of least squares, any deviation from the orthogonal projection gives worse mean deviance and worse DIC as long as the effective number of parameters $p_D$ does not change much. Consequently, the more the oblique projection (constrained approach) deviates from the orthogonal projection, the worse the mean deviance and thus DIC. In the data analysis considered here, the constrained approach has the smallest value of $p_D$ because some degrees of freedom are removed due to the constraints, but DIC is worse because the reduction in $p_D$ does not compensate for the increase in mean deviance.

A desirable property of the constrained approach is that it keeps the (random) spatial and temporal main effects constant in time and space respectively while restricting all random effects to be orthogonal to the fixed effects. However, at least in the dataset analysed here, this behaviour comes at an important price in terms of model fit, unlike other proposals such as, for example, Hughes and Haran (2013). It is not easy to guess when the constrained and restricted regression approaches will lead to similar fits. We suspect that if the covariates do not have a substantial spatio-temporal interaction the constrained approached could work well. This is consistent with the observation that in a spatial analysis of our data for the year 2011 (not shown), both procedures are nearly equivalent. Moreover, fitting spatio-temporal models including only female literacy rate (a covariate with scarcely any spatio-temporal interaction), the difference in fit between both approaches is reduced considerably.

Both procedures have been fitted using a fully Bayesian and a classical approach. Though INLA provides the posterior distributions of all quantities and hence the maximum information, PQL is still a valuable tool that allows model fitting in a reasonable time providing essential information to understand the phenomenon under study. Computing times for restricted regression can be reduced in INLA by plugging in as initial values the posterior modes of the hyperparameters obtained from the spatio-temporal model with confounding. Similarly, restricted regression can be sped up in PQL by using as initial values the variance parameter estimates obtained from the spatio-temporal model with confounding.

To sum up, models with spatial, temporal, and spatio-temporal random effects lead to a good fit of the spatio-temporal patterns of risk, but they fail to account for the correct linear association between the response and the covariates. Models including only covariates provide unbiased estimates of the fixed effects coefficients but they give a poor fit as they do not capture the variability unexplained by the covariates. The restricted regression considered in this article offers in one single model the best of the two approaches: similar to the model without random effects, it provides correct estimates of the fixed effects but substantially improves model fit and prediction of spatio-temporal patterns of risk in line with the model with random effects. This is an advantage of the restricted regression approach over the models without random effects and classical spatio-temporal models. Recent research (Khan and Calder (2020)) concludes that for Gaussian models credible intervals for fixed effects obtained with restricted regression are nested within credible intervals obtained with models without spatial random effects. These authors acknowledge that restricted regression inference on the coefficients is very similar to the one from non-spatial models. In the Poisson spatio-temporal case study analysed here, credible intervals for some of the coefficients in the restricted regression approach are not nested within those of the model without random effects, but differences are small. Further research is needed to investigate this issue more in depth. The use of constraints also provides correct estimates of the fixed effects and improves model fit, but to a lesser extent than the restricted regression approach, at least in the dataset analysed here. We would like to highlight that in this work, associations between the covariates and the response should be understood as correlations, with the final objective of identifying good predictors of the response, and not as causal relationships within a causal-inference framework. This latter approach is beyond the scope of this article.

Finally, we want to emphasize the consequences of ignoring confounding for the data analysed here. Dowry death in India, particularly in Uttar Pradesh, is a complex problem for which risk factors (socio-demographic, economic, cultural or religious) are not yet clearly identified. Ignoring confounding may lead researchers to discard some potential risk factors and wrongly estimate their associations with dowry death. In this article, sex ratio, population density, female literacy rate, per capita income, murder rate, and burglary rate were found to be associated with dowry deaths when confounding was taken into account. Ignoring such effects masks the association between dowry deaths and some of those risk factors, which obscures understanding of this atrocious practice that takes the lives of thousands of women in India.

Appendix

This subsection is about adding orthogonality constraints in the PQL approach, in particular about modifying the algorithm to replace the sum-to-zero identifiability constraints with new constraints that identify the model and also achieve the desired orthogonality. Here we use ‘conditioning by kriging’ (pp. 37, 93 Rue and Held (2005)) where the covariance matrix of the random effects conditional on the constraints is used. Consider the covariance matrices

$\begin{array}{l}\text{Cov}\left(\xi \left|{\text{B}}_{\xi }\right.\xi =0\right)={\text{Q}}_{\xi }^{-}-{\text{Q}}_{\xi }^{-}{\text{B}}_{\xi }^{\text{'}}{\left({\text{B}}_{\xi }{\text{Q}}_{\xi }^{-}{\text{B}}_{\xi }^{\text{'}}\right)}^{-1}{\text{B}}_{\xi }{\text{Q}}_{\xi }^{-}\\ \text{Cov}\left(\gamma \left|{\text{B}}_{\gamma }\right.\gamma =0\right)={\text{Q}}_{\gamma }^{-}-{\text{Q}}_{\gamma }^{-}{\text{B}}_{\gamma }^{\text{'}}{\left({\text{B}}_{\gamma }{\text{Q}}_{\gamma }^{-}{\text{B}}_{\gamma }^{\text{'}}\right)}^{-1}{\text{B}}_{\gamma }{\text{Q}}_{\gamma }^{-}\\ \text{Cov}\left(\delta \left|{\text{B}}_{\delta }\right.\delta =0\right)={\text{Q}}_{\delta }^{-}-{\text{Q}}_{\delta }^{-}{\text{B}}_{\delta }^{\text{'}}{\left({\text{B}}_{\delta }{\text{Q}}_{\delta }^{-}{\text{B}}_{\delta }^{\text{'}}\right)}^{-1}{\text{B}}_{\delta }{\text{Q}}_{\delta }^{-}\end{array}$

where ‘ ${}^{-}$ ’ denotes the Moore-Penrose generalized inverse, and the matrices ${\text{B}}_{\xi }$ , ${\text{B}}_{\gamma }$ , and ${\text{B}}_{\delta }$ capture the orthogonality constraints of the previous section and are defined in Equation (3.5). Using these covariance matrices, the PQL algorithm would automatically place the desired orthogonality constraints. However, due to the rank deficiency of the matrices ${\text{Q}}_{\xi }$ , ${\text{Q}}_{\gamma }$ , and ${\text{Q}}_{\delta }$ , the usual sum-to-zero constraints are also imposed in addition to the weighted sum-to-zero constraints. To overcome that problem, consider instead these covariance matrices (see Theorem 1 in Supplementary Material A for details):

$$\begin{array}{l}\text{Cov}\left(\xi \left|{\text{B}}_{\xi }\right.\xi =0\right)={\text{V}}_{\xi }={\text{L}}_{\xi }{\left[{\text{L}}_{\xi }^{\text{'}}{\text{Q}}_{\xi }{\text{L}}_{\xi }\right]}^{-1}{\text{L}}_{\xi }^{\text{'}}\\ \text{Cov}\left(\gamma \left|{\text{B}}_{\gamma }\right.\gamma =0\right)={\text{V}}_{\gamma }={\text{L}}_{\gamma }{\left[{\text{L}}_{\gamma }^{\text{'}}{\text{Q}}_{\gamma }{\text{L}}_{\gamma }\right]}^{-1}{\text{L}}_{\gamma }^{\text{'}}\\ \text{Cov}\left(\delta \left|{\text{B}}_{\delta }\right.\delta =0\right)={\text{V}}_{\delta }={\text{L}}_{\delta }{\left[{\text{L}}_{\delta }^{\text{'}}{\text{Q}}_{\delta }{\text{L}}_{\delta }\right]}^{-1}{\text{L}}_{\delta }^{\text{'}}\end{array}$$

where ${\text{L}}_{\xi }$ , ${\text{L}}_{\gamma }$ , and ${\text{L}}_{\delta }$ are as in Equation (3.4). Doing this, the null spaces of the covariance matrices are now spanned by the vectors of constraints, that is the rows of the matrices ${\text{B}}_{\xi }$ , ${\text{B}}_{\gamma }$ , and ${\text{B}}_{\delta }$ , and the PQL algorithm automatically circumvents identifiability issues and provides estimates satisfying the orthogonality requirements. To see this briefly, the PQL algorithm requires a working vector

$$O^{*}={\text{X}}_{*}\beta +{\text{Z}}_{\xi }\xi +{\text{Z}}_{\gamma }\gamma +{\text{Z}}_{\delta }\delta +(O-\mu )g\prime (\mu ),$$

where ${\text{Z}}_{\xi }=1_T\otimes {\text{I}}_S$ , ${\text{Z}}_{\gamma }={\text{I}}_T\otimes 1_S$ , and ${\text{Z}}_{\delta }={\text{I}}_{\mathit{TS}}$ are the design matrices of the spatial and temporal main effects and the interaction effect respectively; $g\prime (\mu )=1/\mu$ is the derivative of the link function g, which here is the logarithmic function, $\in =\left(O-\mu \right)g\text{'}\left(\mu \right)\sim N\left(0,{\widehat{\text{W}}}^{-1}\right)$ , and $\widehat{W}=diag\left({\mu }_{it}\right)$ . Then the fixed effect estimator (including the intercept) is $$\widehat{\beta }={({\text{X}}_{*}{}^{\prime }{\widehat{\text{V}}}^{-1}{\text{X}}_{*})}^{-1}{\text{X}}_{*}{}^{\prime }{\widehat{\text{V}}}^{-1}{O}^{*}$$ , where $\text{V}={\text{W}}^{-1}+{\text{Z}}_{\xi }{\text{V}}_{\xi }{\text{Z}}_{\xi }{}^{\prime }+{\text{Z}}_{\gamma }{\text{V}}_{\gamma }{\text{Z}}_{\gamma }{}^{\prime }+{\text{Z}}_{\delta }{\text{V}}_{\delta }{\text{Z}}_{\delta }{}^{\prime }$ and the random effects are estimated as

$$\hat{\xi }={\hat{\text{V}}}_{\xi }{\text{Z}}_{\xi }{}^{\prime }{\hat{\text{V}}}^{-1}(O^{*}-{\text{X}}_{*}\hat{\beta }),\hat{\gamma }={\hat{\text{V}}}_{\gamma }{\text{Z}}_{\gamma }{}^{\prime }{\hat{\text{V}}}^{-1}(O^{*}-{\text{X}}_{*}\hat{\beta }),\hat{\delta }={\hat{\text{V}}}_{\delta }{\text{Z}}_{\delta }{}^{\prime }{\hat{\text{V}}}^{-1}(O^{*}-{\text{X}}_{*}\hat{\beta }).$$

Clearly, ${\text{B}}_{\xi }\widehat{\xi }=0$ , ${\text{B}}_{\gamma }\widehat{\gamma }=0$ , and ${\text{B}}_{\delta }\widehat{\delta }=0$ , as the rows of ${\text{B}}_{\xi }$ , ${\text{B}}_{\gamma }$ , and ${\text{B}}_{\delta }$ , span the null spaces of ${\text{V}}_{\xi }$ , ${\text{V}}_{\gamma }$ , and ${\text{V}}_{\delta }$ respectively.