Restricted Spatial Regression Methods: Implications for Inference

TL;DR

This paper critically examines Restricted Spatial Regression (RSR) methods, revealing that they often worsen inference in spatial data analysis and cannot be remedied by simple basis vector selection, challenging current practices.

Contribution

The paper provides a unified mathematical framework for RSR methods, analyzes their impact on inference, and demonstrates their counterproductive effects through simulations in both linear and generalized linear models.

Findings

RSR methods can impair inference for regression coefficients.

Counterintuitive results show RSR often performs worse than non-spatial methods.

Problems with RSR cannot be fixed by choosing better spatial basis vectors.

Abstract

The issue of spatial confounding between the spatial random effect and the fixed effects in regression analyses has been identified as a concern in the statistical literature. Multiple authors have offered perspectives and potential solutions. In this paper, for the areal spatial data setting, we show that many of the methods designed to alleviate spatial confounding can be viewed as special cases of a general class of models. We refer to this class as Restricted Spatial Regression (RSR) models, extending terminology currently in use. We offer a mathematically based exploration of the impact that RSR methods have on inference for regression coefficients for the linear model. We then explore whether these results hold in the generalized linear model setting for count data using simulations. We show that the use of these methods have counterintuitive consequences which defy the general expectations in the literature. In particular, our results and the accompanying simulations suggest that RSR methods will typically perform worse than non-spatial methods. These results have important implications for dimension reduction strategies in spatial regression modeling. Specifically, we demonstrate that the problems with RSR models cannot be fixed with a selection of "better" spatial basis vectors or dimension reduction techniques.

Equations

$$\mathbf{Z}=1{\beta }_0+\mathbf{X}{\mathbf{\beta }}^X+\mathbf{\eta }.$$

$$\mathbf{Y}=\mathbf{X}{\mathbf{\beta }}^X+\mathbf{W}\mathbf{\delta }+\mathbf{ϵ}$$

$$p(\mathbf{\delta }|{\tau }_s)\propto {\tau }_s^{\text{rank}(\mathbf{F})/2}\exp \{\frac{-{\tau }_s}{2}{\mathbf{\delta }}^T\mathbf{F}\mathbf{\delta }\},$$

$$\mathbf{Y}=\mathbf{X}{\mathbf{\beta }}^X+\mathbf{V}\mathbf{\delta }+\mathbf{ϵ},$$

$$p(\mathbf{\delta }|{\tau }_s)\propto {\tau }_s^{\kappa }\exp \left\{\frac{-{\tau }_s}{2}{\mathbf{\delta }}^T\mathbf{\Lambda }\mathbf{\delta }\right\},$$

$$\mathbf{Y}=\mathbf{X}{\mathbf{\beta }}^X+\mathbf{L}{\mathbf{\delta }}_1+\mathbf{ϵ},$$

$$p({\mathbf{\delta }}_1|{\tau }_s)\propto {\tau }_s^{\kappa }\exp \{\frac{-{\tau }_s}{2}{\mathbf{\delta }}_1^T{\mathbf{Q}}_{\mathbf{s}}{\mathbf{\delta }}_1\},$$

$$\text{E}[{\mathbf{\beta }}^{\mathbf{X}}|\mathbf{Y}]$$

$$\text{Var}[{\mathbf{\beta }}^{\mathbf{X}}|\mathbf{Y}]$$

$$\frac{\text{E}[r|\mathbf{Y}]}{E[{\tau }_s]}<E\left[\frac{1}{{{\tau }_{ϵ}}^{NS}}|\mathbf{Y}\right]=E\left[{\sigma }_{ϵ}^{NS}|\mathbf{Y}\right],$$

$$\widehat{\text{Var}}({\widehat{{\mathbf{\beta }}^X}}_i)=({\mathbf{X}}^T\mathbf{X}{)}_{ii}^{-1}{\hat{\sigma }}_{ϵ}^{NS}=({\mathbf{X}}^T\mathbf{X}{)}_{ii}^{-1}\frac{||{{\mathbf{P}}_X}^{\perp }\mathbf{Y}|{|}^2}{n-p}.$$

$$\text{Var}[{\mathbf{\beta }}_{NS,i}^X|\mathbf{Y}]=({\mathbf{X}}^T\mathbf{X}{)}_{ii}^{-1}\text{E}({\sigma }_{ϵ}^{NS}|Y)=({\mathbf{X}}^T\mathbf{X}{)}_{i,i}^{-1}\frac{b_{ϵ}^{-1}+.5||{{\mathbf{P}}_X}^{\perp }\mathbf{Y}|{|}^2}{(a_{ϵ}-1)+.5(n-p)}.$$

$$\widehat{\text{Var}}({\widehat{{\mathbf{\beta }}^X}}_i)=({\mathbf{X}}^T\mathbf{X}{)}_{ii}^{-1}\widehat{{\sigma }_{\mathbf{T}}}=({\mathbf{X}}^T\mathbf{X}{)}_{ii}^{-1}\frac{||{{\mathbf{P}}_T}^{\perp }\mathbf{Y}|{|}^2}{n-p-q}.$$

$$\frac{\frac{1}{b_{ϵ}}+\frac{r}{b_s}+.5{\mathbf{Y}}^T({{\mathbf{P}}_X}^{\perp }-\mathbf{W}(\mathbf{I}+r\mathbf{F}{)}^{-1}{\mathbf{W}}^T)\mathbf{Y}}{a_{ϵ}+a_s+.5(n-p)-1}.$$

$$\underset{r\to 0}{\lim }E[{\sigma }_{ϵ}^{RSR}|\mathbf{Y},r]=\frac{1}{b_{ϵ}(a_{ϵ}+a_s+.5(n-p)-1)}.$$

$$s_l\sqrt{n-1}\mathbf{V}{\mathbf{\rho }}_{lV}+1\bar{l}$$

$$f(y_i|{\theta }_i)\exp \{y_i{\theta }_i-b({\theta }_i)+c(y_i)\}$$

$$g({\mu }_i)={\theta }_i={\mathbf{X}}_i{\mathbf{\beta }}^X+(\mathbf{V}\mathbf{\delta }{)}_i+{\mathbf{ϵ}}_i$$

$$\mathbf{\theta }$$

$${\mathbf{\delta }}_1$$

$${\mathbf{ϵ}}_1$$

$${{{\mathbf{P}}_X}^{\perp }}^{\ast }=\mathbf{I}-{\hat{\mathbf{H}}}^{1/2}\mathbf{X}({\mathbf{X}}^T\hat{\mathbf{H}}\mathbf{X}{)}^{-1}{\mathbf{X}}^T{\hat{\mathbf{H}}}^{1/2}$$

$$\mathbf{Z}={\beta }_0+{\mathbf{X}}_1{\beta }_1^X+{\mathbf{X}}_2{\beta }_2^X+{\hat{\mathbf{H}}}^{-1/2}\mathbf{L}\mathbf{\delta }$$

$$\mathbf{\delta }\sim N(0,{\tau }_s{\mathbf{L}}^T{\hat{\mathbf{H}}}^{1/2}\mathbf{Q}{\hat{\mathbf{H}}}^{1/2}\mathbf{L})$$

$$\mathbf{Z}={\beta }_0+{\mathbf{X}}_1{\beta }_1^X+{\mathbf{X}}_2{\beta }_2^X+{\mathbf{M}}_q^{\ast }\mathbf{\delta }$$

$$\mathbf{\delta }\sim N(0,{\tau }_s{{\mathbf{M}}^{\ast }}^T\mathbf{Q}{\mathbf{M}}^{\ast })$$

$$f({\mathbf{\beta }}^X|\delta ,{\tau }_s,{\tau }_{ϵ},\mathbf{Y})$$

$$\text{E}[{\mathbf{\beta }}^X|\mathbf{Y}]$$

$$=\text{E}\left[({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{Y}|\mathbf{Y}\right]$$

$$=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{Y},$$

$$\text{Var}[{\mathbf{\beta }}^X|\mathbf{Y}]$$

$$=\text{E}\left[({\tau }_{ϵ}{\mathbf{X}}^T\mathbf{X}{)}^{-1}|\mathbf{Y}\right]$$

$$=({\mathbf{X}}^T\mathbf{X}{)}^{-1}\text{E}\left[\frac{1}{{\tau }_{ϵ}}|\mathbf{Y}\right]$$

$$\displaystyle\Big{|}\iiiint_{\bm{\beta},\bm{\delta},\tau_{s},\tau_{\epsilon}}\frac{1}{\tau_{\epsilon}}f(\bm{\beta},\bm{\delta},\tau_{s},\tau_{\epsilon}|\bm{Y})d\tau_{\epsilon}d\tau_{s}d\bm{\delta}d\bm{\beta}\Big{|}$$

$$={\text{iiiint}}_{\mathbf{\beta },\mathbf{\delta },{\tau }_s,{\tau }_{ϵ}}\frac{1}{{\tau }_{ϵ}}f(\mathbf{\beta },\mathbf{\delta },{\tau }_s,{\tau }_{ϵ}|\mathbf{Y})d{\tau }_{ϵ}d{\tau }_{s}d\mathbf{\delta }d\mathbf{\beta }$$

$$\propto {\iint }_{{\tau }_{ϵ},{\tau }_s}\exp \left\{\frac{-{\tau }_{ϵ}}{b_{ϵ}}\right\}\exp \left\{\frac{-{\tau }_s}{b_s}\right\}{{\tau }_s}^{A_{{\tau }_s}}{\tau }_{ϵ}^{A_{{\tau }_{ϵ}}}\frac{1}{\sqrt{|{\mathbf{W}}^T\mathbf{W}+\frac{{\tau }_s}{{\tau }_{ϵ}}\mathbf{F}|}}\exp \left\{\frac{-{\tau }_{ϵ}}{2}{\mathbf{Y}}^T{\mathbf{\Sigma }}_{{\tau }_{ϵ},{\tau }_s}\mathbf{Y}\right\}d{\tau }_{ϵ}d{\tau }_s$$

$$\left|{\mathbf{W}}^T\mathbf{W}+\frac{{\tau }_s}{{\tau }_{ϵ}}\mathbf{F}\right|$$

Full text

Keywords: Confounding; Spatial Statistics; Bayesian; Dimension Reduction

Kori Khan111Email: khan.746@osu.edu, Department of Statistics, The Ohio State University

and

Catherine A. Calder, Department of Statistics and Data Sciences,

University of Texas at Austin

1. Introduction

In our increasingly data rich world, large spatial data sets are becoming abundant. A booming field of research involves incorporating spatial dependence into models in a computationally efficient way. As is often the case with spatial statistics, much of this research has initially focused on developing and analyzing methods for spatial process models: models developed to make predictions at unobserved locations (e.g., Banerjee et al.,, 2008; Wikle,, 2010; Fuentes,, 2007; Stein,, 2014). There have been few attempts to understand how these methods impact inference on regression coefficients.

Recently, however, there has been a line of work designed to efficiently incorporate spatial dependence into regression models when the primary interest is on inference for the regression coefficients. In many ways, this set of work mirrors the work done in the spatial process model research. For some examples, in the context of areal data, Hughes and Haran, (2013) presented a reduced-rank approach, Prates et al., (2019) developed a sparse approximation technique, and Burden et al., (2015) introduced an approximate likelihood method. Bradley et al., (2015) extended ideas of basis selection techniques to multivariate spatio-temporal mixed effects models, while Murakami and Griffith, (2015) incorporate spatial dependence through eigen-vector spatial filtering. More recently, Thaden and Kneib, (2018) introduced a structural equation approach for estimating regression coefficients when there is an unobserved spatially dependent covariate.

In contrast to the work in the spatial process model realm, these methods are either designed or inspired by methods designed to alleviate “spatial confounding.” The first explicit reference to spatial confounding is often attributed to Clayton et al., (1993), who observed what they referred to as “confounding by location”: the situation where estimates of a regression coefficient associated with a spatially-structured covariate are affected by the presence of a spatial random effect in the model. In recent work, the phenomenon is often explained as the presence of multicollinearity between the covariates, $\bm{X}$, and the spatial random effect (Prates et al.,, 2019; Hanks et al.,, 2015; Thaden and Kneib,, 2018; Hefley et al.,, 2017).

Reich et al., (2006) and Hodges and Reich, (2010) brought recent attention to the issue (see also, Paciorek,, 2010, for a study on the effects of spatial confounding on inference for regression coefficients). These works highlighted for the first time that in the presence of spatial confounding, not only can the estimates for regression coefficients change, but the uncertainty associated with these estimates can be “overinflated.” To address these dual concerns, Reich et al., (2006) proposed a method employing synthetic predictors to smooth orthogonally to the fixed effects. Hughes and Haran, (2013) then extended this work by suggesting the use of eigenvectors of the Moran operator as synthetic predictors. They argued these basis functions would allow for dimension reduction through the selection of only those synthetic predictors associated with “attractive” spatial dependence (as opposed to “repulsive” spatial dependence). The intuitive appeal of this set of work inspired a series of follow-up investigations (Prates et al.,, 2019; Burden et al.,, 2015; Bradley et al.,, 2015; Thaden and Kneib,, 2018). However, it has been observed that these methods can lead to elevated levels of Type-S errors, the Bayesian analogue of Type I error (Prates et al.,, 2019; Hanks et al.,, 2015, the latter paper suggested a posterior predictive approach to address this concern). Recently, Hanks et al., (2015) and Prates et al., (2019) have noted that there is a need to better understand when it is appropriate to utilize such methods. Indeed, with the exception of Reich et al., (2006), there has been a lack of mathematical formalism to assess the impact of smoothing orthogonally to the fixed effects has on inference for regression coefficients.

In this work, we propose to fill that void. To do so, we consider a Bayesian analysis of Gaussian areal spatial data. We show that with respect to inference on regression coefficients, the current methods proposed to smooth orthogonally to the fixed effects can be thought of as a subset of a larger class of models. Extending the terminology currently in use for these proposed methods, we will refer to this larger class as Restricted Spatial Regression (RSR) models. We find that RSR models transform a mixed-effects model into an over-fit linear model. Specifically, any of these models will produce a posterior mean for the regression coefficients which is equivalent to the posterior mean obtained in the corresponding non-spatial model. The various approaches to smoothing orthogonally to fixed effects were designed to ensure that all the marginal posterior variances of the regression coefficients are greater than the corresponding posterior variances of the non-spatial model. However, we show that the exact opposite is often true. The reduction in the posterior variances is caused solely by the addition of the spatial basis functions. Furthermore, our analytic results and the included simulation studies indicate that sufficiently large credible intervals for the regression coefficients will generally be nested within the corresponding credible intervals from the non-spatial model. Importantly, these results are invariant to the spatial structure of the covariates and any spatial structure in the residuals. In short, our results indicate that with respect to coverage and Type-S error, one would be better off fitting a non-spatial model than any RSR model. Furthermore, simulations indicate that when there is spatial dependence unexplained by the covariates, RSR models for Gaussian data exasperate the decrease in coverage and increase in Type-S rates - even if the true unexplained spatial dependence in the response variable is generated orthogonally to the covariates. A simulation study and example application of RSR models suggest similar results might hold for count data.

The rest of the paper proceeds as follows: Section 2 provides an overview of deriving inference on regression coefficients in a spatial random effects model and a discussion of the methods developed to alleviate spatial confounding. Section 3 contains the results of this paper. Section 4 includes simulation studies, and Section 5 compares various approaches on the Slovenia stomach cancer data set. The proofs for all theorems are contained in Appendix B.

2. Background

2.1. Inference on Regression Coefficients in Random Effects Models

The spatial generalized linear mixed model (SGLMM), popularized by Diggle et al., (1998), assumes that $\bm{Y}=\{y_{i},\ldots,y_{n}\}^{T}$ is a realization from a random field where $y_{i}$ is observed at spatial location $s_{i},~{}i=1\ldots n$. The vector of transformed conditional means $\bm{Z}=\{g(E(y_{1}|\eta(s_{1}))),\ldots,g(E(y_{n}|\eta(s_{n}))\}^{T}$ for a given link function $g$ is then related to fixed effects and a spatial random effect:

[TABLE]

Here, $\bm{1}$ is the $n\times 1$ column vector of 1’s, $\bm{X}=[\bm{X}_{1},\ldots,\bm{X}_{p}]$ is the $n\times p$ design matrix whose $i$th row consists of the $p$ covariates associated with $\bm{Z}_{i}$ , $\bm{\beta}^{*}=(\beta_{0},~{}\bm{\beta}^{X})\in\mathbb{R}^{p+1}$ is a vector of regression parameters, and $\bm{\eta}=(\eta_{1},\ldots,\eta_{n})^{T}$ is a zero-mean random effect with a spatial covariance matrix $\bm{\Sigma}(\bm{\theta})$ parameterized by $\bm{\theta}\in\mathbb{R}^{m}$. In this paper, the notation denoting dependence on location will be dropped when the meaning is unambiguous. We define ${\bm{X}}^{*}\equiv[\bm{1}~{}\bm{X}]$. The distinction between ${\bm{X}}$ and ${\bm{X}}^{*}$ is important when we compare different models. However, unless otherwise stated, the points about one in this section are true for both.

Historically, much of the statistical community’s intuition for the behavior of estimates of regression coefficients is based on work for the linear model when $g$ is the identity mapping. A non-spatial (NS) analysis in this case would correspond to an assumption that $\bm{\Sigma}(\theta)\propto\bm{I}$. In this setting, the best linear unbiased estimator for the regression coefficients is the well-known ordinary least squares (OLS) estimator $\bm{\hat{\beta}}_{NS}^{*}=({{\bm{X}}^{*}}^{T}{\bm{X}}^{*})^{-1}{{\bm{X}}^{*}}^{T}\bm{Y}$. Classical statistics textbooks remind us that if the assumption of i.i.d error structure is not met, these OLS estimates will be inefficient. Instead, if $\bm{\theta}$ were known, the most efficient estimator would be the generalized least squares estimator (GLS) $\bm{\hat{\beta}}_{GLS}^{*}=({\bm{X}^{*}}^{T}\bm{\Sigma}^{-1}(\bm{\theta})\bm{X}^{*})^{-1}{\bm{X}^{*}}^{T}\bm{\Sigma}(\bm{\theta})^{-1}\bm{Y}$.

The two estimates will generally not be the same. However, in spatial statistics there is often an assumption that the addition of a spatial random effect should not change the point estimates of the regression coefficients. This expectation seems to be driven by the fact that geostatistical software developed to facilitate spatial process modeling often only implemented ordinary kriging and not universal kriging, the former allowing for an unknown, but constant mean and the latter allowing the mean to be an unknown linear combination of covariates associated with points in space (Waller and Gotway,, 2004, pg. 344). Conventional wisdom in the spatial statistics literature suggests that if there is spatial dependence unexplained by the covariates (i.e, processes in which things near in space are assumed to be more similar than things far away in space), the naive use of $\bm{\hat{\beta}}_{NS}^{*}$ will underestimate the variance of this estimator. This intuition is driven by well-understood examples which illustrate that ignoring spatial dependence will result in an underestimate of the variance of the mean of a spatial process (see, e.g., Cressie,, 1993, pages 13-15). Paciorek, (2010) pointed out that there is a lack of formal quantification of this belief for regression coefficients in general. Working in a setting where $\bm{X}$ is stochastic, Paciorek, (2010) showed that generally the naive variance estimator $\textrm{Var}(\hat{\bm{{\beta}}}^{X}_{OLS})$ will underestimate the uncertainty estimate of the correct variance estimator $\textrm{Var}(\hat{\bm{{\beta}}}^{X}_{GLS})$, thereby lending support to the common expectation in spatial statistics. In any case, the expectation that the estimates of the variance associated with regression coefficient estimators will increase in models accounting for the presence of residual spatial dependence extends beyond the linear model. In practice this expectation is often true, but there are examples in the literature where it does not hold (for an example, see Banerjee et al.,, 2003).

2.2. Spatial Confounding

The term spatial confounding is used to describe the effect that multicollinearity between the fixed covariates $\bm{X}$ and the spatial random effect $\bm{\eta}$ has on the point estimates and variance of the regression coefficient estimates. Most of the work in spatial confounding involves a fully Bayesian analysis where the point estimates for the regression coefficients are defined to be the posterior mean $E(\bm{\beta}^{X}|\bm{Y})$ and the variances of interest are the diagonal entries of $\textrm{Var}(\bm{\beta}^{X}|\bm{Y})$ (Reich et al.,, 2006; Hughes and Haran,, 2013; Hefley et al.,, 2017; Hanks et al.,, 2015; Prates et al.,, 2019). Efforts have been made to propose statistics designed to detect the presence of spatial confounding (Reich et al.,, 2006; Hefley et al.,, 2017; Prates et al.,, 2019). Thaden and Kneib, (2018) provided a formalization of spatial confounding that assumes the multicollinearity between $\bm{X}$ and the spatial random effect occurs because of the absence of another unobserved spatially varying covariate. However, there is not currently a formal definition of spatial confounding in a more general setting, in part because it is a difficult concept to quantify. It is easier to define what spatial confounding is not, and so that is what we do in 1. As will be seen shortly, all the methods designed to address spatial confounding are designed in the hopes of achieving the properties in 1.

Definition 1**.**

A method which results in posterior mean $E(\bm{\beta}^{X}|\bm{Y})$ and marginal posterior variances $\textrm{Var}(\bm{\beta}^{X}_{i}|\bm{Y})$, $i=1,\ldots,p$ alleviates spatial confounding if the following conditions are met:

1. (1)

$E(\bm{\beta}^{X}|\bm{Y})=E(\bm{\beta}_{NS}^{X}|\bm{Y})$, and 2. (2)

$\textrm{Var}({\bm{\beta}_{NS,}^{X}}_{i}|\bm{Y})\leq\textrm{Var}({\bm{\beta}}_{i}^{X}|\bm{Y})\leq\textrm{Var}({\bm{\beta}_{Spatial,}^{X}}_{i}|\bm{Y})$ for $i=1,\ldots,p$,

where $\bm{\beta}_{NS}^{X}$ are the regression coefficients of the corresponding non-spatial model and $\bm{\beta}_{Spatial}^{X}$ are the regression coefficients from an unrestricted spatial random effect.

Hodges and Reich, (2010) argued that spatial confounding can be a concern whenever a spatial random effect is included in a model. However, most of the proposed methods to address spatial confounding are developed in the context of Gaussian areal spatial data (Reich et al.,, 2006; Hodges and Reich,, 2010; Hughes and Haran,, 2013; Prates et al.,, 2019). In the areal data setting, spatial dependence is described by the introduction of an underlying, undirected graph $G=(V,E)$. Non-overlapping spatial regions that partition the study area are represented by vertices, $V=\{1,\ldots,n\}$, and edges $E$ defined so that each pair $(i,j)$ represents the proximity between region $i$ and region $j$. We represent $G$ by its $n\times n$ binary adjacency matrix $\bm{A}$ with entries defined such that $\textrm{diag}(\bm{A})=0$ and $\bm{A}_{i,j}=\mathbbm{1}_{(i,j)\in E,i\neq j}$. In models designed to address spatial confounding, the spatial random effect is typically assumed to follow the intrinsic conditional autoregressive (ICAR) model/prior (Besag et al.,, 1991). We will refer to models utilizing this prior as ICAR models.

For Gaussian areal data, we observe that the ICAR model, the non-spatial model, and models designed to alleviate spatial confounding such as those proposed by Reich et al., (2006) (RHZ), Hughes and Haran, (2013) (HH), and Prates et al., (2019) (PAR) are all special cases of the following more general form:

[TABLE]

where $\bm{\epsilon}\sim N(\bm{0},\tau_{\epsilon}\bm{I})$, and $\tau_{\epsilon}$ and $\tau_{s}$ are precision parameters. $\bm{W}$ is a $n\times q$ set of basis vectors such that $q\leq n-p$, and $\bm{F}$ is a symmetric, non-negative definite matrix. Table 1 illustrates how the three proposed RSR methods can be considered special cases of the form (2).

[TABLE]

The distinction between the use of $\bm{X}$ and $\bm{X}^{*}$ in Table 1 is a direct consequence of the impropriety of the ICAR prior. The ICAR model captures spatial dependence through a Gaussian Markov Random Field (GMRF). The precision matrix is the graph Laplacian $\bm{Q}=\textrm{diag}(\bm{A}\bm{1})-\bm{A}$. Because the adjacency matrix $\bm{A}$ is defined to be the zero/one adjacency matrix on $G$, $\bm{Q}$ is the graph Laplacian of a simple graph. A well-known fact from spectral graph theory is that the kernel of $\bm{Q}$ will be constant on each connected component of $G$ (Von Luxburg,, 2007). This fact means that the prior $p(\bm{\delta}|\tau_{s})$ will include an implicit intercept for each connected component of the graph in Bayesian analysis. For simplicity, this paper will assume that the $G$ is connected and therefore there is a single global intercept included in the ICAR prior.

Reich et al., (2006) investigated in depth the effect multicollinearity between $\bm{X}$ and the spatial random effect induced by the ICAR prior has on $E(\bm{\beta}^{X}|\bm{Y})$ and $\textrm{Var}(\bm{\beta}^{X}|\bm{Y})$ using a re-parameterization of the ICAR model:

[TABLE]

where $\bm{V}\bm{\Lambda}\bm{V}^{T}$ is the eigendecomposition of $\bm{Q}$, and $\bm{\Lambda}=\textrm{diag}(\lambda_{1},\ldots,\lambda_{n})$ where $\lambda_{1}\geq\lambda_{2},\ldots\geq\lambda_{n}$. Under the assumptions of this paper (i.e., that $G$ is connected), $\kappa=\frac{n-1}{2}$. Using a Bayesian approach with flat priors on $\bm{\bm{\beta}^{X}}$ and independent gamma priors on $\tau_{s},~{}\tau_{s}$, Reich et al., (2006) and Hodges and Reich, (2010) illustrated that $E[\bm{\bm{\beta}^{X}}|Y,\tau_{s},\tau_{\epsilon}]\neq E[\bm{\beta}_{NS}^{X}|\bm{Y}]$ in the presence of multicollinearity between $\bm{X}$ and $\bm{V}$. When regression coefficients are of primary interest, they argued that spatial random effects are added merely to account for spatial correlation in residuals when computing the posterior variance. Under this reasoning, the spatial random effect should not change the point estimates of the fixed effects. Furthermore, Reich et al., (2006) suggested that this multicollinearity causes “overinflation” of the posterior variance of the regression coefficients. This overinflation was noted to be exacerbated when $\bm{X}$ was correlated with low frequency eigenvectors of the graph Laplacian. In spatial regression, Reich et al., (2006) argued that two factors can lead to an increase in the posterior variance of regression coefficients as compared to the NS model: 1) collinearity with the spatial random effects and 2) a reduction in the effective number of observations because of spatial clustering. They proposed a model that they claimed would eliminate the former of these factors while preserving the latter. In other words, the RHZ model was designed to alleviate spatial confounding according to 1:

1. (1)

$E(\bm{\beta}_{RHZ}^{X}|\bm{Y})=E(\bm{\beta}_{NS}^{X}|\bm{Y})$ 2. (2)

$\textrm{Var}({\bm{\beta}_{NS,}^{X}}_{i}|\bm{Y})\leq\textrm{Var}({\bm{\beta}_{RHZ,}^{X}}_{i}|\bm{Y})\leq\textrm{Var}({\bm{\beta}_{ICAR,}^{X}}_{i}|\bm{Y})$ for $i=1,\ldots,p$.

The authors provided theoretical support for first inequality of the second claim.222Section 3 of Reich et al., (2006) provides a mathematical argument that this should be true. However, the argument appears to rely on the stochastic ordering of a gamma distribution. It appears the authors assume they are using a gamma distribution with a scale parameter- which is stochastically increasing in the scale. Their work suggests they are instead working with a rate parameterization- which is stochastically decreasing in the rate parameter. To accomplish these goals, the authors proposed using synthetic predictors which are orthogonal to the column space of $\bm{X}$. Define $\bm{P}_{X}=\bm{X}({\bm{X}}^{T}\bm{X})^{-1}\bm{X}^{T}$ to be the projection matrix onto the column space of $\bm{X}$, and ${\bm{P}_{X}}^{\perp}=\bm{I}-\bm{P}_{X}$ to be its orthogonal complement. Let $\bm{L}$ be the $n\times(n-p)$ matrix composed of the eigenvectors of ${\bm{P}_{X}}^{\perp}$ associated with an eigenvalue of $1$ and let $\bm{K}$ be the $n\times p$ matrix composed of the eigenvectors associated with an eigenvalues of [math]. Model (3) can be re-written as a function of $\bm{\delta}_{1}=\bm{L}^{T}\bm{V}\bm{\delta}$ and $\bm{\delta}_{2}=\bm{K}^{T}\bm{V}\bm{\delta}$. The RSR method proposed by Reich et al., (2006) involves setting $\bm{\delta}_{2}=\bm{0}$. Using this notation, the RHZ model can be expressed as:

[TABLE]

where $\bm{Q_{s}}=\bm{L}^{T}\bm{Q}\bm{L}$. Hodges and Reich, (2010) explained that this solution assigns all variability explained by both $\bm{V}$ and $\bm{X}$ to $\bm{X}$ by restricting spatial smoothing to the orthogonal complement of the fixed effects.

Hughes and Haran, (2013) also argued that variance inflation caused by spatial confounding is a concern in the SGLMM. However, they stated that the RHZ model (4) fails to account for the underlying graph in the construction of $\bm{L}$, thereby permitting structure in the spatial random effects that corresponds to negative spatial dependence. Hughes and Haran, (2013) suggested such negative spatial dependence should not be expected in the context in which spatial models are generally fit. They developed a model (hereafter the HH model) that they stated would eliminate negative dependence while alleviating the impact of spatial confounding on the marginal posterior variances. Utilizing a Bayesian framework with proper priors on the regression coefficients, they suggested that rather than using $\bm{L}$, one should instead use $\bm{M}_{q}$, where $\bm{M}_{q}$ is the $n\times q$ matrix composed of the eigenvectors of what they referred to as the Moran operator ${\bm{P}_{X}}^{\perp}\bm{A}{\bm{P}_{X}}^{\perp}$. In support of the argument that $\bm{M}_{q}$ should be preferred to $\bm{L}$, the authors pointed out that the column vectors of $\bm{M}_{q}$ had more spatial structure than the column vectors of $\bm{L}$.

Hughes and Haran, (2013) also stated that this model naturally lent itself to dimension reduction. The Moran operator is the numerator of a generalization of the Moran’s $I$ statistic, a popular non-parametric measure of spatial dependence introduced by Moran, (1950). We will refer to this generalization of the Moran’s $I$ as $I_{\bm{X}}(\bm{A})$. Relying on work on the original Moran’s $I$ statistic in the context of bounded regular tessellations by Boots and Tiefelsdorf, (2000), Hughes and Haran, (2013) noted the (standardized) spectrum of the Moran operator comprises all possible values of $I_{\bm{X}}(\bm{A})$ and its eigenvectors comprise of all possible mutually distinct patterns of clustering after accounting for $\bm{X}$ and the graph. Importantly, they argued that choosing only “attractive” eigenvectors (those associated with positive eigenvalues) of the Moran operator would improve the inference on regression coefficients by eliminating patterns of “repulsive” spatial dependence (eigenvectors associated with negative eigenvalues). Furthermore, they suggested that for most graphs only selecting the first $.1n$ eigenvectors should be sufficient to perform well for regression. The work in Hughes and Haran, (2013) relied on simulations to illustrate the points made.

A downside to both the RHS and HH models is a loss of some of the computational efficiency inherent in the ICAR prior utilized in (3). ICAR models capture spatial dependence through a Gaussian Markov Random Field (GMRF). This is appealing because the precision matrix for the random field is very sparse, which allows the use of sparse matrix routines to facilitate efficient computations (Rue and Held,, 2005; Paciorek,, 2009). The matrix $\bm{Q}_{s}$ need not be, and generally is not, sparse. Furthermore, the RHZ and HH methods do not extend naturally to spatial models that do not utilize the ICAR prior. Recent work by Prates et al., (2019) attempts to to overcome the computational limitations of RSR methods. Their work suggests projecting the vertices of the graph $G$ onto the orthogonal space of the design matrix $\bm{X}$. Using the projected vertices, they then construct a new sparse precision matrix $\bm{Q}^{\perp}$ for use in analysis. This model is referenced as the PAR model in Table 1. Guan and Haran, (2018) also suggested a projection based approach to approximate covariance matrices of spatial random effects to improve the computational efficiency of RSR methods. This method can also be extended to spatial models beyond the ICAR model.

RSR methods are not universally accepted; some find the assumption that the spatial random effect operates orthogonally to the fixed effects to be too strong (See e.g., Paciorek,, 2010). Furthermore, recent work has highlighted that RSR methods can suffer from elevated rates of Type-S errors under model misspecification (Hanks et al.,, 2015; Prates et al.,, 2019). Type-S error is the Bayesian analogue to Type I error, and we define it to occur when the equal-tailed 95% credible interval for a regression coefficient that is truly 0 does not include 0. The model misspecification studied has involved a generating model which include spatial random effects that do not operate orthogonally to the fixed effects. Currently, the prevalent belief is that the elevated Type-S error is caused by such misspecification. However, in this paper we will show that the phenomena occurs even without it.

It remains an open question when RSR methods should be preferred over traditional spatial models. Understanding when RSR methods perform well is important because RSR methods can lead to a very different interpretations of the effect of regression coefficients (Hanks et al.,, 2015). Attempts to understand when RSR methods are appropriate have led to suggestions that when $\bm{X}$ is correlated with low frequency eigenvectors of the graph Laplacian, RSR methods should be preferred (Reich et al., (2006); Prates et al., (2019), Hefley et al., (2017)).

3. Consequences of Orthogonal Smoothing

In this section, we investigate the impact that methods designed to alleviate spatial confounding have on inference for regression coefficients. Formally, we state the following definition of a RSR model:

Definition 2**.**

A Restricted Spatial Regression Model is any model of the form (2) with $\bm{W}$ chosen such that $\mathcal{C}(\bm{W})\perp\mathcal{C}(\bm{X})$.

Here and throughout the paper $\mathcal{C}(\cdot)$ denotes the column space of a matrix. Both the RHZ and HH models are special cases of this class of models. Note that the ICAR model is not a special case of a RSR model. For RSR models, any results stated for $\bm{\beta}^{X}$ (relying on $\bm{X}$) are also true for $\bm{\beta}^{*}$ (by replacing $\bm{X}$ with $\bm{X}^{*}$). Since the PAR model is an attempt to approximate these methods, we will not spend time directly investigating its performance. All of the following results assume a Bayesian analysis that uses flat priors on the regression coefficients and independent gamma priors on the precision parameters $\tau_{s}$ and $\tau_{\epsilon}$ with respective shape parameters $a_{s},~{}a_{\epsilon}>0$ and scale parameters $b_{s},~{}b_{\epsilon}>0$.

3.1. Properties of Marginal Posterior Distribution

To begin to investigate the impact RSR models have on inference for regression coefficients, we consider the mean and variance of the marginal posterior distribution of $\bm{\beta}^{X}$. As previously noted, RSR models are developed in the areal data setting with the expectation that the marginal posterior mean of $\bm{\beta}^{X}$ will be $({\bm{X}}^{T}\bm{X})^{-1}\bm{X}^{T}\bm{Y}$. Working in a geostatistical setting and utilizing an empirical Bayes approach, Hanks et al., (2015) observed this would be the case for a RSR method in that context as well. However, all of the previous work has relied either on observations regarding distributions of $\bm{\beta}^{X}$ conditional on the precision parameters (for example, Reich et al., (2006); Hodges and Reich, (2010)) or on the full conditional distribution of $\bm{\beta}^{X}$ (Hanks et al., (2015)). 1 offers the first rigorous proof of the existence of the mean and variance of the marginal posterior distributions of $\bm{\beta}^{X}$ in RSR models, as well as the forms these moments take. It illustrates that the point estimates obtained from RSR models will be the same as the non-spatial point estimates. It also shows the variance of the distribution looks functionally similar to the variance we would have obtained in the non-spatial model. For the rest of the paper, we adopt the following notational shortcut: $\sigma_{\epsilon}=\frac{1}{\tau_{\epsilon}}$.

Theorem 1**.**

Under conditions A.1[1]- A.1[5], a RSR model with $\bm{W}$ a $n\times q$ matrix will give rise to the marginal posterior distributions of $\bm{\beta}^{X}$ such that:

[TABLE]

Proof: See Appendix Section B.1.

An immediate consequence of 1 is that the dimension of $\bm{W}$ will not impact the point estimates for the regression coefficients. Thus, any choice of $q$ linearly independent column vectors which are orthogonal to $\mathcal{C}(\bm{X})$ will give the same point estimates as a HH model using $\bm{M}_{q}$. Of course, incorporating spatial dependence in a model is a computationally expensive and often theoretically complicated endeavor. If the concern were simply to obtain point estimates from a non-spatial model, there would be no need to use a spatial model. As noted previously, the focus in fitting a spatial model is often to account for spatial correlation in residuals when computing the posterior variance. Hence, our interest lies in $\textrm{Var}[\bm{\bm{\beta}^{X}}|\bm{Y}]$, which 1 illustrates is a function of $E[{\sigma}_{\epsilon}^{RSR}|\bm{Y}]$.

The prevalent expectation is that the marginal posterior variances of the regression coefficients in RHZ and HH models will be greater than the marginal posterior variances obtained in the non-spatial model. In settings where RSR models are currently recommended for use, 2 illustrates that for a broader class of RSR models including the RHZ and HH model, the opposite is true. In the statement of 2, $r=\frac{\tau_{s}}{\tau_{\epsilon}}$ is the inverse of the signal-to-noise ratio.

Theorem 2**.**

Under conditions A.1[1]-A.1[4], a RSR model with $\bm{W}$ a $n\times q$ matrix with orthonormal columns and $\bm{F}$ a positive definite and symmetric matrix will always result in a marginal posterior variance for ${\bm{\beta}^{X}}_{i},\,i=1,\ldots,p$ that is less than or equal to that of the posterior variance that would have been obtained in the non-spatial model whenever the following holds:

[TABLE]

where $E[{\tau_{s}}]$ is the prior expectation of $\tau_{s}$ under the RSR model.

Proof: See Appendix Section B.2.

As a practical matter, the conditions set forth in 2 will often be met. To see this, consider that spatial regression models, including RSR models, are typically used in settings in which the $\bm{X}\bm{\beta}^{X}$ in the NS model does not explain nearly all the variability in $\bm{Y}$ (i.e., when $E\left[\sigma_{\epsilon}^{NS}|\bm{Y}\right]$ is not nearly zero). Spatial confounding is thought to be problematic in settings where $r$ is “small” for traditional spatial regression models (see e.g., Reich et al., 2006 and Hodges and Reich, 2010). “Small” is usually taken to mean substantially less than $1$ (see e.g., Reich et al., 2006 and Prates et al., 2019 who consider values of $r=.1$ and $.2$ respectively). This corresponds to situations where there is residual spatial dependence unexplained by $\bm{X}$ (and hence larger $E\left[\sigma_{\epsilon}^{NS}|\bm{Y}\right]$ ). The better a RSR model is at explaining this residual spatial dependence, the smaller $\textrm{E}(r|\bm{Y})$ will be. Finally, in the literature, $E[{\tau_{s}}]$ ranges in value from $1$ (see e.g., Hodges and Reich,, 2010) to $1,000$ (see e.g., Hughes and Haran,, 2013). Thus, the conditions in 2 will typically be satisfied for prominent examples of RSR models.

2 offers the first fully rigorous proof detailing the bounds on the variance of the marginal posterior distribution of $\bm{\beta}^{X}$ for RSR models. Working in a geostatistical setting with an empirical Bayesian analysis, Hanks et al., (2015) stated without proof that the inequality related to posterior variances in 2 should hold for a particular RSR model. However, perhaps because of the difference in settings, this observation has not been perceived as contradicting the expectations for the more commonly studied RHZ and HH models, which are precisely the opposite. Interestingly, the work in Hanks et al., (2015) is often reported as illustrating that RSR methods can have elevated Type-S error rates under model misspecification (Prates et al.,, 2019; Page et al.,, 2017). Prates et al., (2019) further noted that these elevated Type-S error rates were exacerbated in settings in which $r$ is small. Yet, 2 suggests that RSR models will have elevated Type-S error rates in the absence of model misspecification. In a frequentist setting, the OLS estimators would be normally distributed and an ordering of variance with the same mean would suggest that the $100(1-\alpha)\%$ confidence intervals would be nested (with the interval of the distribution with less variance completely contained in the other). In that setting, 2 would indicate that any RSR model would have higher rates of Type I error and lower rates of Type II error than the non-spatial model. Since the non-spatial model yields estimators with sampling distributions whose coverage probability matches the nominal coverage probability, this indicates that the RSR coverage probability would always be less than the nominal coverage probability under the conditions of 2. Hence, 2 would explain the elevated levels of Type-S error observed for RSR models, and it would mean these elevated levels could exist even without model misspecification.

In Bayesian analysis, when inference on regression coefficients is of primary interest, practitioners typically utilize an equal-tail credible interval. This is especially true in settings like those of this paper where improper priors preclude the use of Bayes factors. Because the marginal posterior distribution $f(\bm{\beta}^{X}|\bm{Y})$ is not normal for any of the models discussed, the variance ordering in 2 need not necessarily indicate a relationship between the credible intervals.

Investigating the relationship between equal-tailed credible intervals for the non-spatial model and a RSR model for all possible choices of graph, data, and spatial basis vectors is difficult because these distributions will not generally be available in closed form. However, for the case that $\beta^{X}\in\mathbb{R}$, it is possible to make some observations that shed light on the relationship between the credible intervals of a RSR model and the NS model. This is formalized in 3.

Theorem 3**.**

Assume $\beta^{X}\in\mathbb{R}$, A.1[3] holds, and A.1[5] holds. Let $g(\beta^{X}|\bm{Y})$ be the marginal posterior probability distribution (pdf) from a RSR model with choice of $\bm{F}$ which is a symmetric and positive definite $q\times q$ matrix. Let $h(\beta^{X}|\bm{Y})$ be the marginal posterior pdf from the non-spatial model. Then, $g(\beta^{X}|\bm{Y})=\textrm{O}(h(\beta^{X}|\bm{Y}))$ as $\beta^{X}\to\infty$ and $\beta^{X}\to-\infty$.

Proof: See Appendix Section B.3.

An immediate corollary to this result is the following.

Corollary 1**.**

Define $G$ and $H$ to be the respective cumulative distribution functions of $g(\beta^{X}|\bm{Y})$ and $h(\beta^{X}|\bm{Y})$. $\exists~{}C>0$ such that:

1. (1)

$\limsup\limits_{\beta^{X}\to-\infty}\frac{G(\beta^{X}|\bm{Y})}{H(\beta^{X}|\bm{Y})}\leq C,$ 2. (2)

$\limsup\limits_{\beta^{X}\to\infty}\frac{1-G(\beta^{X}|\bm{Y})}{1-H(\beta^{X}|\bm{Y})}\leq C$

Furthermore, for $D_{g}^{*}$ and $K(C_{q},D_{h})$ as defined in the proof in Appendix B, if $D_{g}^{*}>K(C_{q},D_{h})$, then $\exists~{}\beta^{*}>0,~{}\beta_{*}<0$ such that:

1. (1)

$G(\beta|\bm{Y})\leq H(\beta|\bm{Y})~{}\forall\beta<\beta_{*}$ 2. (2)

$1-G(\beta|\bm{Y})\leq 1-H(\beta|\bm{Y})~{}\forall\beta>\beta^{*}$

Proof: See Appendix Section B.3.

Generally, what these results indicate is that the tails of the marginal posterior distribution for the non-spatial model and the tails of the marginal posterior distribution for RSR models decay roughly at the same rate. The combination of 1 and 3 suggests that that RSR models offer inference very similar to a non-spatial model with respect to regression coefficients. Importantly, these results are invariant to the choice of graph, data, and spatial basis vectors. To see how the inference will be similar, note the posterior distribution $f(\beta^{X}_{NS}|\bm{Y})$ is symmetric about its mean. While, $f(\beta^{X}_{RSR}|\bm{Y})$ is not symmetric about its mean, two applications of Jensen’s inequality can bound the distance between a mean and median by the standard deviation. In practice, the standard deviation for a RSR model is always quite small. Thus, RSR and NS models produce marginal posterior distributions with a similar medians and similar tail decay. This result is important because RSR models are not offering much beyond what an NS model offers, yet they are much more computationally expensive. For example, Hughes and Haran, (2013) reported that a NS model took under a minute, while their dimension reductions technique took over 6 hours.

If the condition put forth in Corollary (1) holds, it also means that for sufficiently large credible intervals, the credible interval for the RSR model will be completely contained in the corresponding credible interval for the non-spatial model. Importantly, these results require no assumptions about $\textrm{E}(r|\bm{Y})$. Just as in the frequentist setting, this means that RSR methods will only capture the regression coefficient if the non-spatial model does. Furthermore, the RSR methods will always have higher rates of Type-S error than the non-spatial model - even if there is no spatial confounding. However, in practice determining whether the condition in Corollary (1) is satisfied will be impractical. Therefore, we further the investigation with simulation studies in Section 4.

3.2. Implications for Inference

The premise of RSR models is that incorporating spatial dependence will yield better inference on the regression coefficients. In particular, the HH model and those influenced by it have emphasized choosing spatial basis vectors which are spatially smooth. This emphasis is inherited from work on dimension reduction techniques developed for spatial process modelling where this technique is quite useful (See e.g., Cressie and Wikle,, 2015, pages 353-354). Recall our definition of RSR models includes both models designed to capture spatial dependence as well as models with arbitrary choices of $\bm{W}$ and $\bm{F}$. Yet, the results of 1 – 3 hold whether the models are designed to incorporate spatial dependence or not. A natural question then is: For RSR models, do choices of $\bm{W}$ (and perhaps $\bm{F}$) which exhibit spatial patterns behave any differently than those which do not?

4 begins to explore the answer to this question. Recall that one of the criticisms of the RHZ model (4) is that the basis vectors $\bm{L}$ do not appear to be spatially smooth. 4 states that when $\bm{F}$ is of the form $\bm{W}^{T}\bm{Q}\bm{W}$, the inference on the regression coefficients will be invariant to $\mathcal{C}(\bm{W})$. This result suggests that the choice of basis vectors which exhibit spatial patterns (such as those from the Moran operator) does not affect inference on regression coefficients in the ways expected. As an example, for moderate $q$, given a selection of spatially smooth basis vectors, we could pick a different basis of the column space (which may not appear spatially smooth), and inference for $\beta^{X}$ would be equivalent.

Theorem 4**.**

A RSR model with ${\bm{W}}_{1}$ a $n\times q$ matrix with orthonormal columns and ${\bm{F}}_{1}={\bm{W}}_{1}^{T}\bm{B}{\bm{W}}_{1}$ for arbitrary non-null symmetric $\bm{B}$ will yield the same marginal posterior distribution $f(\bm{\beta}^{X}|\bm{Y})$ as any other choice ${\bm{W}}_{2}$ with orthonormal columns such that $\mathcal{C}({\bm{W}}_{1})=\mathcal{C}({\bm{W}}_{2})$ and ${\bm{F}}_{2}={\bm{W}}_{2}^{T}\bm{B}{\bm{W}}_{2}$.

Proof: See Appendix Section B.4.

4 suggests that the assumptions underlying the development of RSR models may be worth re-evaluating. The addition of the spatial basis vectors in RSR models is assumed to make the regression coefficient estimates less biased and more precise. However, further investigation into RSR models indicates that by choosing these spatial basis vectors to be orthogonal to the column space of $\bm{X}$, RSR models are effectively transforming a mixed-effects model into an overfit fixed effects model. Thus, RSR models may not be aiding in inference for regression coefficients as expected.

To understand the intuition behind this last claim, it is instructive to take a moment to review the non-spatial Gaussian model. The method of ordinary least squares makes the assumption that all the variability in the direction of $\bm{X}$ can be described by the linear combination of $\bm{X}\bm{\beta}^{X}$. This assumption results in fitted values which are simply a projection of $\bm{Y}$ onto the column space of $\bm{X}$. In other words, the point estimates are $\hat{\bm{\beta}^{X}}=({\bm{X}}^{T}{\bm{X}})^{-1}{\bm{X}}^{T}\bm{Y}$. Importantly, this assumption is also reflected in the estimated variance of these point estimators. This can be seen in the familiar form for the estimated variance of the OLS estimators, where the $i^{th}$ regression coefficient $\hat{{\bm{\beta}^{X}}}_{i}$ is simply the $(i,i)$ element of $({\bm{X}}^{T}{\bm{X}})^{-1}\hat{\sigma}_{\epsilon}^{NS}$:

[TABLE]

The first component reflects the relationship between the different covariates, which is not of interest in this paper. The second term, however, is a consequence of the assumption that the variation in the direction of $\bm{X}$ can be explained by $\bm{X}\bm{\beta}^{X}$. Here, $\hat{\sigma}_{\epsilon}^{NS}$ is simply a function of the magnitude of the component of $\bm{Y}$ unexplained by any linear combination of the column space of $\bm{X}$.

In Bayesian analysis, the assumption that all the variability in the direction of $\bm{X}$ can be explained by a linear combination of the columns of $\bm{X}$ has analogous implications. For the non-spatial model, the point estimates for the regression coefficients will simply be those of the OLS model. The associated posterior variance for the $i^{th}$ regression coefficient ${\bm{\beta}^{X}}_{i}$ will be the $(i,i)$ element of $({\bm{X}}^{T}{\bm{X}})^{-1}\textrm{E}(\sigma_{\epsilon}^{NS}|Y)$:

[TABLE]

Although a bit messier than the OLS setting, this relationship again shows that the uncertainty associated with the regression coefficients is a function of the magnitude of the component of $\bm{Y}$ unexplained by the the column space of $\bm{X}$. For $n\gg p$ and the typical choices of hyper-parameters ($a_{\epsilon},b_{\epsilon}^{-1}<1$), this will in fact be quite similar to the estimate obtained via OLS.

In a fixed effects model, researchers are well aware that overfitting the model with the addition of synthetic predictors can distort the true relationship between $\bm{Y}$ and $\bm{X}$. For example, consider the non-spatial model with new design matrix $\bm{T}=[\bm{X}~{}\bm{W}]$, with $\mathcal{C}(\bm{X})\perp\mathcal{C}(\bm{W})$. Because $\bm{W}$ is constructed to belong to the orthogonal complement of the columns space of $\bm{X}$, the point estimates remain unchanged. If $\bm{W}$ was a $n\times(n-p)$ matrix of rank $n-p$, $\bm{T}$ would be a basis of $\mathbb{R}^{n}$. Thus, by construction, ${\bm{P}_{T}}^{\perp}\bm{Y}=\bm{0}$. More generally, because $({\bm{T}}^{T}{\bm{T}})^{-1}$ is a block diagonal matrix with diagonals $({\bm{X}}^{T}{\bm{X}})^{-1}$ and $({\bm{W}}^{T}{\bm{W}})^{-1}$, the following will be true:

[TABLE]

For fixed $\bm{X}$, $||{\bm{P}_{T}}^{\perp}\bm{Y}||^{2}$ is monotonically decreasing as a function of the number of columns of $\bm{W}$. If we add enough (linearly independent) columns to $\bm{W}$, we will eventually explain the variation in $\bm{Y}$. Of course, the $n-p-q$ in the denominator ensures that $\hat{\sigma_{\bm{T}}}$ is not monotonically decreasing as the number of columns of $\bm{W}$ increases. However, a basis for the orthogonal complement to $\mathcal{C}(\bm{X})$ can still be added in such a way as to ensure that $\hat{\sigma_{\bm{T}}}$ will be very, very close to zero (see Section C.1 for an example). Thus, even if there is truly no relationship between $\bm{Y}$ and $\bm{X}$, with enough synthetic covariates this approach can essentially ensure that the covariates are deemed “significant” (assuming that the columns $\bm{X}$ are not very strongly correlated). Of course, as noted this problem is well-understood in fixed effect models, and this approach is not considered acceptable in the practice of regression analysis.

Returning to RSR models, 1 illustrates that like the NS model, the point estimates for these models assume that all the variation in the direction of $\bm{X}$ can be explained by the linear combination $\bm{X}\bm{\beta}^{X}$. However, 2 suggests that when it comes to the posterior variance, RSR models are going beyond this assumption. In RSR models, $\bm{F}$ and $\tau_{s}$ act as regularization parameters that can smooth some of the elements of $\bm{\delta}$ to 0. These terms make it difficult to find a closed form expression for $E[\sigma_{\epsilon}^{RSR}|\bm{Y}]$. However, we can still gain intuition by looking at the behavior of $E[\sigma_{\epsilon}^{RSR}|\bm{Y},r]$, where $r=\frac{\tau_{s}}{\tau_{\epsilon}}$ is again the inverse of the signal-to-noise ratio. Straightforward calculations show that $E[\sigma_{\epsilon}^{RSR}|\bm{Y},r]$ can be expressed as:

[TABLE]

Note that the $(n-p)$ term in the denominator will remain unchanged no matter how many columns are added to the spatial basis matrix $\bm{W}$. If we consider the behavior of $E[\sigma_{\epsilon}^{RSR}|\bm{Y},r]$ as the signal-to-noise ratio increases (i.e., $r\to 0$) for $\bm{W}$ a $n\times(n-p)$ matrix of full rank, then

[TABLE]

For moderate $p$ (relative to $n$) and the typical choices of $a_{\epsilon}$ and $b_{\epsilon}$, this value will be quite small. Note that this is slightly smaller than the value that would be obtained in a Bayesian NS model if the fixed effects $\bm{X}$ alone completely explained the variation in $\bm{Y}$ (see Eq. 5 for the case when $||{\bm{P}_{X}}^{\perp}\bm{Y}||^{2}=0$). However, this reduction is a function of how well the additional columns of $\bm{W}$ explain the variability in $\bm{Y}$, and the signal-to-noise ratio will tend to increase as “useful” columns are added to $\bm{W}$. In other words, if $\bm{X}$ did a poor job of explaining the variability of $\bm{Y}$, the addition of columns to $\bm{W}$ will lead to a more dramatic decrease in $E[\sigma_{\epsilon}^{RSR}|\bm{Y}]$ and hence $\textrm{Var}[\bm{\beta}^{X}|\bm{Y}]$. Because of the regularization induced by $\bm{F}$ and $\tau_{s}$, this decrease will not necessarily be monotonic.

To understand the implications of this discussion as well as 4 on inference, we analyze a real world dataset with a series of RSR models. We consider the SAT data originally analyzed by Wall, (2004). We use the data set provided as an example in Bivand et al., (2008). The data include statewide averaged verbal scores on the SAT and the percent of students eligible to take the test in 1999 for the 48 contiguous states of the United States. Let $y_{i}$ denote the statewide averaged SAT verbal score for the $i$th state and $\bm{Y}=(y_{i},\ldots,y_{48})$. If $\bm{X}_{1}$ = Percent of students eligible to take the exam, then we consider a model with covariates $\bm{X}^{*}=[\bm{1}~{}\bm{X}_{1}~{}\bm{X}_{1}^{2}]$ with associated regression coefficients $\beta_{0},~{}\beta_{1},~{}\beta_{2}$. We fit a series of HH models for various $M_{q},~{}q=1,\ldots,45$. For all these analyses, $a_{s}=.5,~{}b_{s}=2000$ as suggested by Hughes and Haran, (2013).

In Fig. 1, we plot the marginal posterior variance for $\beta_{0}$,  $\beta_{1}$, and $\beta_{2}$ for each of the various HH models. The green in the graphs represents vectors which are “attractive” and the red represents vectors which are “repulsive.” All RSR models result in approximately the same point estimates as those obtained as in the non-spatial model $(590.5,-2.84,.022)$ and have marginal posterior variances less than those obtained in the non-spatial model. Based on the discussion above, we expect that as the number of columns of $\bm{W}$ increases, the posterior variances should tend to decrease. And in fact, that seems to be the general pattern. Since the “attractive” vectors are added to the model first, these tend to have a larger variance than models that include both “attractive” vectors and “repulsive” vectors. Note that this same pattern can be obtained by adding arbitrary synthetic covariates to a non-spatial model (see, Section C.1).

However, note, that the distinction between “attractive” and “repulsive” eigenvectors is not directly linked to the posterior variances otherwise. For example, the posterior variance obtained for the set of “attractive” vectors proposed by Hughes and Haran, (2013) is achieved again for choices that include up to $14$ “repulsive” vectors, but is not obtained with the recommended dimension reduction to $.1n$. A restriction to “attractive” eigenvectors was originally proposed as a way to approximate the inference from the RHZ model (See, Section 7 of Hughes and Haran,, 2013). However, there is no clear association between the variances that are obtained via the RHZ model and those obtained restricting the spatial basis vectors to “attractive” ones.

It has been claimed that RSR models assume that all the variability in the direction of $\bm{X}$ can be described by the linear combination $\bm{X}\bm{\beta}^{X}$ (Hodges and Reich,, 2010; Hanks et al.,, 2015). In this section, we point out that the NS already makes this assumption. In order to gain an intuitive understanding of why the posterior variances are RHZ models are always bounded above by the posterior variance of the corresponding NS model, we offer an alternative interpretation which indicates the the reduction in the posterior variance for RHZ models is akin to what would be observed an overfit fixed effects model. When this intuition holds, we would expect that RHZ models would always have a higher Type-S error rate than NS models even in the absence of model misspecification. Importantly, these problems cannot be resolved by selecting different spatial basis vectors within the class of RSR models, and dimension reduction techniques designed to mimic smoothing orthogonally to the fixed effects will likely suffer from the same problems. However, just as is the case for the implications of 3, exploring this rigorously is impractical. Therefore, in the next section we use simulations to further our investigation.

4. Simulation Studies

In this set of simulations, we investigate the relative performances of the ICAR model, the RHZ model, and the NS model for continuous areal data and count areal data. In particular, we explore how and when the RHZ model offers inference different than the inference which would have been obtained with a NS model. We also explore whether RHZ models improve the inference for regression coefficients relative to the ICAR model in the presence of spatial confounding.

Currently, spatial confounding is thought to be a concern when $\bm{X}$ is spatially dependent and there is residual spatial dependence unexplained by the covariates. Therefore, all of the following simulations involve covariates which are spatially dependent and residual spatial dependence (either with a spatial random effect or with “unobserved” covariate). To simulate the covariates, we rely on previous work in spatial confounding. In the context of areal data, statistics developed by Reich et al., (2006) and Prates et al., (2019) essentially define spatial dependence in $\bm{X}$ to be the existence of correlation between $\bm{X}$ and low frequency eigenvectors of the graph Laplacian. These statistics reflect a common assumption in the literature that the eigenvectors of the graph Laplacian $\bm{Q}$ associated with low eigenvalues exhibit spatially smooth patterns while those associated with high eigenvalues oscillate rapidly. In the continuous case, this phenomenon can be formalized. For the low frequency eigenfunctions, Bernstein estimates show smoothness (see e.g., Zelditch, (2017) Theorem 5.17); whereas for high frequency eigenfunctions, rapid oscillation can be shown using upper bounds on the size of nodal domains (see e.g., Zelditch, (2017) Theorem 13.1). We are not aware of any formalization of the phenomenon in the discrete case. However it does seem true in practice, and we will utilize this assumption in the following work. Therefore, in the following simulations we generate covariates that are correlated with low frequency eigenvectors of the graph Laplacian. To do so, we utilize the fact that on a connected graph, any column vector $\bm{l}\in\mathbb{R}^{n}$ can be written as follows:

[TABLE]

where $s_{l}$ and $\bar{l}$ are the sample standard deviation and sample mean of $\bm{l}$, $\bm{\rho}_{lV}=(\rho_{l,V_{i}},\ldots,\rho_{l,V_{n-1}},0)^{T}$ is a column vector where $\rho_{l,V_{i}}$ is the correlation between $\bm{l}$ and the $i^{th}$ column of $\bm{V}$.

In the following simulation studies, all models are fit using Markov Chain Monte Carlo algorithms with Gibbs updates and/or Metropolis-Hastings random-walk updates. Gibbs updates are used for all parameters in the Gaussian models and for $\tau_{s}$ in the Poisson model. For the Poisson model, the $\bm{\beta}^{X}$ are updated as in Hughes and Haran, (2013) using a random walk with proposal ${\bm{\beta}^{X}}^{[j+1]}\sim N({\bm{\beta}^{X}}^{[j]},\bm{U}^{-1})$ where $\bm{U}$ is the estimated asymptotic covariance matrix from the non-spatial model, and $\bm{\delta}$ is updated with multivariate random walks with a spherical normal proposal. For all models, where applicable, $a_{\epsilon}=a_{s}=.01$ and $b_{\epsilon}=b_{s}=100$. For the Gaussian models, $\bm{\beta}^{X}$ has a flat prior and for the Poisson model $\bm{\beta}^{X}$ has a normal prior with standard deviation of 1000.

To ensure that the Monte Carlo standard errors were sufficiently small, trial simulations were run until a sample path length was found that ensured all models had Monte Carlo standard error $<.01$. For the Gaussian models, a sample path of $80,000$ was sufficient; while for the Poisson models, a sample path of $1,000,000$ was sufficient. Monte Carlo standard errors were calculated using batch means (Flegal et al.,, 2008, 2012).

4.1. Simulation 1: Gaussian Model

For this set of simulations we consider Gaussian models of the form (2) for the graph of the 48 contiguous states. We generate data from three models. All generating models are of the form $\bm{Y}=\bm{1}\beta_{0}+\bm{X}_{1}\beta^{X}+\bm{\nu}+\bm{\epsilon}$ where $\bm{\epsilon}\sim N(\bm{0},\bm{I})$, $\beta_{0}=1,~{}\beta^{X}=2,$ and $\bm{X}_{1}$ is randomly generated to be correlated with $.2n$ of the eigenvectors associated with the lowest non-zero eigenvalues. The vector $\bm{\nu}$ is one of three options: $\bm{0}$ (i.e., the model is a non-spatial linear model), a realization from $\bm{L}\bm{\delta}$ with $\bm{\delta}\sim N(\bm{0},\bm{L}^{T}\bm{Q}\bm{L})$ (the RHZ model (4)), or a realization from the ICAR prior with $\tau_{s}=1$ (realizations from the ICAR random effect are generated using Algorithm 2.6 in Rue and Held, (2005)). For each of the choices of $\bm{\nu}$ we simulated 1000 realizations of $\bm{Y}$ and fit the NS model, RHZ model, and the ICAR model.

Table 2 lists these coverage probabilities based on 95% credible intervals with equal tail weights. 98% of the credible intervals from an analysis using the RHZ model were nested within the corresponding credible intervals derived from the NS model. To explore how the RHZ analysis model differs from the NS analysis model, we consider how the coverage differed for each data set in Table 3. We consider three outcomes:

1. (1)

The credible intervals for both models either both contained the true value of $\beta^{X}$ or both failed to contain it (Agree) 2. (2)

the credible interval of the RHZ model contained the true value of $\beta^{X}$ and the NS model did not (RHZ +) 3. (3)

The credible interval of the NS model contained the true value of $\beta^{X}$ and the RHZ model did not (NS +).

With respect to coverage, both the ICAR model and the NS model obtain higher coverage rates than the RHZ model, even when the data is generated from the RHZ model. When the data is generated from a NS or RHZ model, the differences between the three models are not substantial. These differences are driven primarily by the widths of the credible intervals. The ICAR model tends to have the largest width and the RHZ model tends to have the narrowest credible intervals.

Arguably, the setting of primary interest for this and the simulation studies that follow is when the data is generated with an ICAR prior. This setting is often thought to be most closely related to what would be seen in the real world. When the generating model is the ICAR model, both the RHZ model and the NS model achieve their lowest coverage rates. This suggests, as observed in Prates et al., (2019) and Hanks et al., (2015), that the RHZ model may suffer poor coverage in the presence of spatial random effects which do not operate orthogonally to the fixed effects. In investigating the cause of this poor performance, we found that it was a combination of bias and narrow credible intervals. Here, we define bias as $\beta^{X}$ - $E[\beta^{X}|\bm{Y}]$. In Section C.2.1, we find that the bias for the point estimates for the RHZ and NS models is at its highest when the generating model is an ICAR model. Since by 1, the point estimates of the RHZ and NS models are the same, the further reduction in coverage rates for the RHZ are caused by narrower credible intervals. There is not a single case in which the credible interval of the RHZ model obtained coverage when the NS model did not.

[TABLE]

[TABLE]

4.2. Simulation 2: Type-S Error in the Gaussian Model

For this set of simulations we consider Gaussian models of the form (2) for the graph of the 48 contiguous states. We again generate data from three models. All models are of the form $\bm{Y}=\bm{1}\beta_{0}+\bm{X}_{1}\beta^{X}_{1}+\bm{X}_{2}\beta^{X}_{2}+\bm{\nu}+\bm{\epsilon}$ where $\bm{\epsilon}\sim N(\bm{0},\bm{I})$, $\beta_{0}=1,~{}\beta^{X}_{1}=2,$ and $\beta^{X}_{2}=0$. $\bm{X}_{1}$ is again randomly generated by the process described above to be correlated with $.2n$ of the eigenvectors associated with the lowest non-zero eigenvalues, and $\bm{X}_{2}$ is independently generated to be correlated with the $.5n$ eigenvectors with the lowest non-zero eigenvalues. Once again, $\bm{\nu}$ is one of the three options listed in Simulation 1.

With respect to inference on $\bm{\beta}^{X}_{1}$, we again consider coverage rates. For inference on $\beta^{X}_{2}$ we now consider the Type-S error rate: the proportion of times a credible interval does not contain zero. As in Simulation 1, we use 95% equal-tailed credible intervals. We are interested in Type-S errors in the context of data which have spatial variation unexplained by the fixed effects. When the generating model is the RHZ or the ICAR model, $\bm{\nu}$ itself provides such variation. In the context of a non-spatial model, the only way to achieve such dependence is by omitting a spatially varying covariate. Thus, when the generating model is the RHZ or ICAR, we use a design matrix of $\bm{X}^{*}=[\bm{1}~{}\bm{X}_{1}~{}\bm{X}_{2}]$. When the generating model is the non-spatial model, we use $\bm{X}^{*}=[\bm{1}~{}\bm{X}_{2}]$, omitting the spatially varying $\bm{X}_{1}$. Recall that the associated $\bm{X}$ will always be used when fitting the ICAR model.

The results are contained in Table 4- Table 7. Both the ICAR model and the NS model have higher coverage rates than the RHZ model in every scenario. As before, there is not a single synthetic data set in which the resulting credible interval obtained coverage for $\beta^{X}_{1}$ under the RHZ model and the NS model did not. The RHZ model once again suffers from poor coverage when the generating model is the ICAR model. The differences in the coverage rates are driven by the same factors as discussed in Section 4.1.

[TABLE]

[TABLE]

In Table 6, the RHZ model has extremely inflated Type-S error rates - even when the data is generated from a RHZ model. Notably, the RHZ model performs the worst with respect to Type-S errors when a spatially varying covariate is omitted from the design matrix. In Section C.2.2, we find that the bias for the NS and RHZ point estimates is at its highest in this setting. This bias leads to the elevated Type-S error rates for the NS model. The extreme inflation from the RHZ model relative to the NS model is caused by narrow credible intervals. We also compare the performance of the RHZ and NS analysis models for each synthetic data set with respect to Type-S error in Table 7. Every data set for which the NS model resulted in a Type-S error, the RHZ model did as well. However, there were a a number of data sets for which the RHZ model resulted in a Type-S error, and the NS model did not.

Table 6. Type-S Error of $\beta^{X}_{2}$

Analysis Model

Generating Model

NS

RHZ

ICAR

NS

13.9%

8.5%

10.3%

RHZ

80.2%

18.0%

21.7%

ICAR

3.7%

4.8%

6.5%

[TABLE]

4.3. Simulation Study 3: Investigating Type-S Error in Poisson Model

Much of the research in the literature as well as the results of this paper focus on Gaussian models. However, in practice most of the models discussed are typically used for areal count data. It is not altogether clear whether work in the linear case can be translated into results for the non-linear case. In this set of simulations, we attempt to investigate whether the results for the linear case are relevant for count data. To our knowledge, this is the first extensive simulation study of count data for RSR methods.

For this set of simulations, we use the graph of 194 Slovenia municipalities considered in Reich et al., (2006) and Prates et al., (2019). We consider count data generated from a Poisson model of the form $\bm{Z}=\beta_{0}+\bm{X}_{1}\beta^{X}_{1}+\bm{X}_{2}\beta^{X}_{2}+\bm{\nu}$ using the log link function. The coefficients are defined as follows: $\beta_{0}=1,~{}\beta^{X}_{1}=1$, and $\beta^{X}_{2}=0$. We generate 100 datasets for each scenario. The set-up of the simulation study is otherwise the same as in Section 4.2. Fitting the RHZ model for the Poisson model is different from that of the Gaussian model. The details of the RHZ model in the non-linear context are given in Section A.2.

It should be noted that the credible intervals from the RHZ model are nearly always the same as the credible intervals from the NS model. As an example, Fig. 2 shows the credible intervals associated with 25 of the 300 simulated datasets. In Table 8 and Table 9 we see that that the RHZ model performs essentially the same as the NS model. There are three cases in which the RHZ credible intervals lead to different conclusions than the NS credible intervals. In these cases, the respective bounds for the NS and RHZ credible intervals are within .05 of one another.

Unlike in the Gaussian case, the performance of the NS model suffers greatly in the presence of spatial variation unexplained by the covariates. This poor coverage seems to be driven by both bias and very narrow credible intervals. In Section C.2.3, we find that the bias for point estimates obtained by the RHZ and NS models is always greater than that of the ICAR model, regardless of the generating model. This phenomena is distinct from the patterns observed in the Gaussian case. With respect to coverage rates and Type-S error, the bias is compounded with very narrow credible intervals for the RHZ and NS models, as illustrated by Table 9. Although the Poisson model has not been heavily studied in simulations, we note that Hughes and Haran, (2013) observed a similar phenomenon with a single simulated dataset. In that setting, the data was simulated from a HH model. The widths of the 95% credible intervals obtained from a NS, HH and ICAR models were, respectively, $.24$, $.24$, and $1.84$. It is not clear why the credible intervals of the NS model are so narrow.

[TABLE]

[TABLE]

5. Application to Slovenia Data

In light of our findings, we now revisit the dataset which motivated the creation of RSR models. Reich et al., (2006) and Hodges and Reich, (2010) used the Slovenia stomach cancer data as an example of how spatial confounding could distort inference on regression coefficients when they observed that the 95% credible interval for the NS model did not contain 0 and the respective credible interval for the ICAR model did contain 0. They believed that this discrepancy occurred because $\textrm{E}[\bm{\beta}_{ICAR}^{X}|\bm{Y}]$ is biased and $\textrm{Var}[\bm{\beta}_{ICAR}^{X}|\bm{Y}]$ is overinflated when $\bm{X}$ displays spatial dependence. This perspective and the Slovenia dataset has continued to influence subsequent research (See e.g., Prates et al., (2019)). However, the results of this paper offer an alternative interpretation of why there is a discrepancy between the NS and ICAR models. To explore this, we now introduce the data and corresponding models in more detail.

The data was collected from 1995 to 2001 for each of the 194 municipalities in Slovenia. A full description of the data is available in Hodges, (2016). The response variable $\bm{Y}$ is defined so that $y_{i}$ is the observed count of stomach cancer cases for the $ith$ municipality. We will let $E_{i}$ be the expected count of stomach cancer cases and $\bm{X}=[x_{1},\ldots,x_{194}]$ be the vector of standardized socioeconomic scores for the municipalities (see, Hodges and Reich,, 2010, for a description). We consider five Poisson models for this model: 1) the non-spatial, 2) the RHZ model, 3) the HH model with all the attractive eigenvectors chosen, 4) the HH model with $q=19$, and 5) the ICAR model. In other words, we assume that the $y_{i}$ are conditionally independent Poisson random variables with rate parameter $\lambda_{i}$. For the non-spatial model, $\log(\lambda_{i})=\log(E_{i})+\beta_{0}+x_{i}\beta^{X}$. The RHZ model and the HH models are of the form: $\textrm{log}(\lambda_{i})=\log(E_{i})+\beta_{0}+x_{i}\beta^{X}+\left(\bm{W}\bm{\delta}\right)_{i}$. The details of the choices of $\bm{W}$ and the prior for $\bm{\delta}$ are given in Section A.2. The ICAR model takes the form $\textrm{log}(\lambda_{i})=\log(E_{i})+x_{i}\beta^{X}+\bm{\delta}_{i}$ where $\bm{\delta}$ is assumed to have the ICAR prior. For all models, a normal prior with mean $\bm{0}$ and standard deviation 1,000 is given to the regression coefficients and $\tau_{s}$ is given a gamma prior with shape and scale, respectively, .01 and 100.

[TABLE]

In Table 10, the ICAR model’s credible interval is the only one to include zero. In applications, researchers would typically conclude that the socioeconomic scores were not significant for the ICAR model but were significant in the other models. The prevailing belief in spatial confounding literature is that the RSR and NS models offer better inference than the ICAR model for this dataset. However, in Section 4.3, we observed that if there is residual spatial dependence (either through a missing spatially varying covariate or through a spatial random effect), the NS and RSR models tend to have more biased point estimates than the ICAR models. Similarly, the NS and RSR models suffer from low coverage rates in this setting. Thus, if there is residual spatial dependence unexplained by $\bm{X}$, the NS and RSR models will likely offer poor inference for $\bm{\beta}^{X}$.

Unfortunately, the ICAR model also has drawbacks in this setting. Note unlike in Section 4.3, the point estimates in Table 10 are quite small. In this setting, the ICAR model will still have smaller bias and higher rates of coverage than the NS and RHZ models (see Section C.6). However, the Bayesian analogue of power, which we define as percentage of time that the 95% credible interval does not contain 0, will be very low for the ICAR model (see again Section C.6). Pratically speaking, 1 - 3 suggest that RSR models give essentially the same inference as the NS model in the Gaussian case, and Section 4.3 seems to support this insight extends to count data. In particular, we note that for each of the RSR models, the 95% credible intervals are nested within (or equivalent to) the credible interval obtained with the NS model. While it may not be clear when the NS or ICAR model should be preferred, it does seem likely that RSR models will typically be a more computationally intensive method of arriving at the same inference as the NS model with respect to bias and coverage rates.

6. Discussion

In this paper we examine the impact that methods designed to alleviate spatial confounding have on inference for regression coefficients of areal SGLMMs. We emphasis that our inferential focus is limited to assessing whether there is a linear association between the response and the covariates in settings where an assumption of (spatial) independence for the errors is apparently invalid (e.g., evidence of residual spatial dependence in the non-spatial model). We note that the distinction between inference related to a linear association and causal inference is often blurred in the literature related to spatial confounding. For instance, the work in Paciorek, (2010) and Thaden and Kneib, (2018) is arguably more closely related to concerns in causal inference. Teasing apart whether common expectations for regression coefficients in spatial models come from an interest in inference for causal relationships or a linear association will likely be important in future related work. However, the implications of RSR methods for formal causal inference are beyond the scope of the current paper.

The results of Section 3 reexamine the intuition that prompted the development of RSR methods. Despite the expectation that RSR methods will result in marginal posterior variances that are greater than would have been obtained in a non-spatial model, we have shown the opposite is usually true. Recently, there has been a call for more research to understand when the RSR methods should be used. This interest is driven in part by the fact that some RSR methods have high rates of Type-S error, particularly when $\frac{\tau_{s}}{\tau_{\epsilon}}$ is small. Currently, RSR methods are recommended when the covariates $\bm{X}$ are correlated with low-frequency eigenvectors of the graph Laplacian. Section 3 and Section 4 illustrate how and why the elevated Type-S error rates are occurring for Gaussian RSR models. In particular, these results indicate that it is not the spatial structure of the covariates that drives the performance of RSR models. Rather, the posterior distribution of the regression coefficients for RSR models is primarily a function of how well the spatially varying $\bm{W}$ (in RSR models) explains the residual spatial dependence. The question of when the RSR models should be preferred is typically asked with the implicit comparison is to that of a traditional spatial model. Surprisingly, our results indicate that in all cases studied, one would always be better off using a Gaussian non-spatial model than a Gaussian RSR model. For the count data, the RHZ model and non-spatial model result in almost equivalent inference for the regression coefficients and both perform poorly in the presence of spatial variation unexplained by the covariates. This work also shows that the ICAR model can provide poor inference as well, particularly in the context of small effect sizes.

In spatial statistics, expectations regarding dimension reduction strategies for areal data models have largely been shaped by work in the spatial process modeling realm. For example, in the process model setting, choosing “attractive” spatial vectors of the Moran operator is analogous to modeling low frequency components of spatial variation. This is the basic premise of approaches such as reduced rank approaches in spatial process modeling. It is quite likely that this method performs well for spatial process modeling. However, this paper shows that insights from spatial process modeling can lead to misleading and counterintuitive results when the interest is no longer on predicting the spatial process. Thus, this work has important implications for dimension reduction when the primary interest is deriving inference about covariate effects.

Acknowledgements

Gabriel Khan is credited for suggesting the clever use of Hölder’s Inequality in 1, which resulted in a relatively cleaner proof. The authors were partially supported by grant NIH NICHD R01-HD088545. They also acknowledge the support of The Ohio State University’s Mathematical Biosciences Institute (NSF DMS-1440386) and the Institute for Population Research (NIH NICHD P2C-HD058484-10).

Appendix A

A.1. Conditions

1. (1)

$\bm{F}$ is a $q\times q$ non-negative definite and symmetric matrix. 2. (2)

Either $\textrm{rank}(\bm{F})\geq 2$ or $\textrm{rank}(\bm{F})=1$ and $a_{\epsilon}/2+n/2-p/2-q/2>1/2$. 3. (3)

$\textrm{rank}(\bm{W})=q\leq n-p$. 4. (4)

$\textrm{rank}(\bm{X})=p<n$. 5. (5)

$\bm{F}$ and ${\bm{W}}^{T}\bm{W}$ commute.

A.2. Details of RHZ and HH Models for Count Data

Section 4 of the Reich et al., (2006) details the construction of a model which smooths orthogonally to the fixed effects for non-linear link functions. In particular, they considered data distributed such that $y_{i}$ are independent given $\theta_{i}$ and the probability distribution function is of the form:

[TABLE]

where $E(y_{i}|\theta_{i})=\mu_{i}$, and $g$ is a link function. Using this observation and the logic from the Gaussian case, they suggested the following adjustment to address spatial confounding.

[TABLE]

where now $\bm{L}$ are the eigenvectors associated with eigenvalues of 1 from ${{\bm{P}_{X}}^{\perp}}^{*}$

[TABLE]

and $\hat{\bm{H}}$ is a diagonal matrix with the $(ii)$th element defined to be $\textrm{Var}(y_{i}|\bm{\beta}^{X},\bm{\delta},\bm{\epsilon})$ evaluated at the mode of $p(\bm{\beta}^{X},\bm{\delta},\bm{\epsilon}|\tau_{\epsilon},\tau_{s},\bm{Y})$. In all relevant code, following the choice in Hughes and Cui, (2018), $\hat{\bm{H}}$ is calculated using the values for the parameters in the last iteration of the Iterated Reweighted Least Squares (IWLS) Algorithm for the non-spatial model. The model detailed in Section 4.3 is simply a reduced case of this model excluding the heterogeneity effect $\bm{\epsilon}$. In particular the model is:

[TABLE]

[TABLE]

In the application section of Hughes and Haran, (2013), the authors defined the Moran operator to be ${\bm{P}_{R}}^{\perp}\bm{A}{\bm{P}_{R}}^{\perp}$ and $\bm{M}_{q}^{*}$ to be the eigenvectors of this Moran operator, where ${\bm{P}_{R}}^{\perp}=\bm{I}-\bm{R}^{1/2}\bm{X}(\bm{X}^{T}\bm{R}\bm{X})^{-1}\bm{X}^{T}\bm{R}^{1/2}$. In Hughes and Cui, (2018), the matrix $\bm{R}$ is defined to be a diagonal matrix with $(i,i)$th entry given to be $1/w_{i}$ where $w_{i}$ are the weights from the last iteration of the IWLS algorithm in the non-spatial model. Although not explicitly stated it appears that Hughes and Haran, (2013) then assumed a model of the following form:

[TABLE]

[TABLE]

Since this also appears to be the model used in Hughes and Cui, (2018) and studied in Prates et al., (2019), we utilize this in Section 5.

Appendix B Proofs of Theorems

B.1. Proof of 1

Proof.

Under these assumptions, straightforward calculations give the following,

[TABLE]

where $(\tau_{\epsilon}{\bm{X}}^{T}\bm{X})^{-1}$ is the covariance rather than the precision. This allows us to conclude the following:

[TABLE]

[TABLE]

The integrals computed above rely on two assumptions:

1. (1)

\Big{|}\iiint_{\bm{\delta},\tau_{s},\tau_{\epsilon}}f(\bm{\delta},\tau_{s},\tau_{\epsilon}|\bm{Y})d\tau_{\epsilon}d\tau_{s}d\bm{\delta}\Big{|}<\infty 2. (2)

\Big{|}\iiiint_{\bm{\delta},\bm{\beta},\tau_{s},\tau_{\epsilon}}\frac{1}{\tau_{\epsilon}}f(\bm{\beta},\bm{\delta},\tau_{s},\tau_{\epsilon}|\bm{Y})d\tau_{\epsilon}d\tau_{s}d\bm{\delta}d\bm{\beta}\Big{|}<\infty

To show that these conditions both hold, it is sufficient to note the following:

[TABLE]

where

1. (1)

$\bm{\Sigma}_{\tau_{\epsilon},\tau_{s}}={\bm{P}_{X}}^{\perp}-{\bm{W}}\left({\bm{W}}^{T}\bm{W}+\frac{\tau_{s}}{\tau_{\epsilon}}\bm{F}\right)^{-1}{\bm{W}}^{T}$ 2. (2)

$A_{\tau_{s}}=\textrm{rank}(\bm{F})/2+a_{s}-1$ 3. (3)

$A_{\tau_{\epsilon}}=a_{\epsilon}+n/2-p/2-q/2-2$

The final integral is finite, which verifies the two assumptions above. To show this, the general strategy is the following. The final exponential term can be bounded by a constant using the non-negative definiteness of $\bm{\Sigma_{\tau_{\epsilon},\tau_{s}}}$. We first fix $\tau_{\epsilon}$, then for $\tau_{s}>1$, convergence occurs due to the exponential decay of the Gamma priors. For $\tau_{s}<1$, we can ignore the exponential term in $\tau_{s}$. Doing so, we can use the determinant estimate listed below, Holder’s inequality, and Fubini’s theorem to directly integrate with respect to $\tau_{s}$. After applying some standard estimates, we bound the remaining integral with respect to $\tau_{\epsilon}$ by separating it into two cases and bounding each by the sum of independent convergent gamma kernels.

Lemma 1** (Determinant Estimate).**

Let $\bm{F}$ be a non-negative-definite and symmetric $q\times q$ matrix and ${\bm{W}}^{T}\bm{W}$ be a positive definite $q\times q$ matrix which commutes with $\bm{F}$. Then for $\frac{\tau_{s}}{\tau_{\epsilon}}>0$ and $j=1,\ldots,\textrm{rank}(\bm{F})$, we have the following estimate on $\left|{\bm{W}}^{T}\bm{W}+\frac{\tau_{s}}{\tau_{\epsilon}}\bm{F}\right|$:

[TABLE]

Here $C_{j}>0$ is a constant which depends on $j$ (and $\bm{F}$).

This lemma can be proven by using the following observation. Since ${\bm{W}}^{T}\bm{W}$ and $\bm{F}$ commute, there exists an invertible matrix $\bm{P}$ such that ${\bm{W}}^{T}\bm{W}=\bm{P}\bm{\Omega}\bm{P}^{-1}$ and $\bm{F}=\bm{P}\bm{\Psi}\bm{P}^{-1}$, where both $\bm{\Omega}$ and $\bm{\Psi}$ are diagonal matrices. Let $\omega_{1},\ldots,\omega_{q}$ be the eigenvalues of ${\bm{W}}^{T}\bm{W}$ (i.e., the diagonal of $\bm{\Omega})$ and let $\psi_{1},\ldots,\psi_{\textrm{rank}(\bm{F})}$ be the non-zero eigenvalues of $\bm{F}$ (i.e., the non-zero diagonal entries of $\bm{\Psi}$). For $j=1,\ldots,\textrm{rank}(\bm{F})$, we observe the following:

[TABLE]

Lemma 2** (Non-negative Definiteness of $\bm{\Sigma}_{\tau_{\epsilon},\tau_{s}}$).**

$\bm{\Sigma}_{\tau_{\epsilon},\tau_{s}}={\bm{P}_{X}}^{\perp}-\bm{W}({\bm{W}}^{T}\bm{W}+\frac{\tau_{s}}{\tau_{\epsilon}}\bm{F})^{-1}{\bm{W}}^{T}$ is non-negative definite.

This lemma can be proven using the following observations.

1. (1)

Because ${\bm{W}}^{T}{\bm{W}}$ and $\bm{F}$ are simultaneously diagonalizable and ${\bm{W}}^{T}\bm{W}$ is positive definite, $({\bm{W}}^{T}\bm{W})^{-1}-({\bm{W}}^{T}\bm{W}+\lambda\bm{F})^{-1}$ is non-negative definite for all $\lambda\geq 0$. 2. (2)

For any two matrices $\bm{A}$ and $\bm{B}$, if $\mathcal{C}(\bm{B})\subset\mathcal{C}(\bm{A})$, then $||\bm{P}_{\bm{A}}\bm{Y}||^{2}-||\bm{P}_{\bm{B}}\bm{Y}||^{2}\geq 0$. 3. (3)

If $q=n-p$, then ${\bm{P}_{X}}^{\perp}=\bm{W}({\bm{W}}^{T}{\bm{W}})^{-1}{\bm{W}}^{T}$. Hence by observation 1, ${\bm{P}_{X}}^{\perp}-\bm{W}({\bm{W}}^{T}\bm{W}+\lambda\bm{F})^{-1}{\bm{W}}^{T}$ is non-negative definite for $\lambda\geq 0$. 4. (4)

If $q<n-p$, then ${\bm{P}_{X}}^{\perp}-\bm{W}({\bm{W}}^{T}\bm{W}+\lambda\bm{F})^{-1}{\bm{W}}^{T}={\bm{P}_{X}}^{\perp}-\bm{P}_{\bm{W}}+\bm{P}_{\bm{W}}-\bm{W}({\bm{W}}^{T}\bm{W}+\lambda\bm{F})^{-1}{\bm{W}}^{T}$. By observations 2 and 3, this is the sum of two non-negative definite matrices and hence non-negative definite.

With these lemmas in hand, we now return to the problem of showing (8) is finite. Using the non-negative definiteness of $\bm{\Sigma_{\tau_{\epsilon},\tau_{s}}}$, we note the following:

[TABLE]

In the following steps, we will use $\lessapprox$ to mean there exists $D>0$ such that $\textrm{LHS}\leq D~{}\textrm{RHS}$.

Case 1: $\tau_{s}>1$

By the Lemma 1 above we can conclude the following

[TABLE]

Note that the integrand is the product of two gamma kernels: one with shape parameter $A_{\tau_{s}}+1=\textrm{rank}(\bm{F})/2+a_{s}>0$ and scale parameter $b_{s}>0$; the second with shape parameter $A_{\tau_{\epsilon}}+\textrm{rank}(\bm{F})/2+1=a_{\epsilon}+n/2-p/2-q/2-1+\textrm{rank}(\bm{F})/2$ and scale parameter $b_{\epsilon}>0$.

Note that if $\textrm{rank}(\bm{F})\geq 2$, then $A_{\tau_{\epsilon}}+\textrm{rank}(\bm{F})/2>0$. If $\textrm{rank}(\bm{F})=1$ and $a_{\epsilon}+n/2-p/2-q/2>1/2$, then $A_{\tau_{\epsilon}}+\textrm{rank}(\bm{F})/2>0$. Thus, under our assumptions this integral is finite.

Case 2: $\tau_{s}<1$

Again utilizing the determinant estimate:

[TABLE]

Then using Fubini’s theorem, reorganizing terms, and substituting the value of $A_{\tau_{s}}$ back in:

[TABLE]

Then using Hölder’s inequality,

[TABLE]

Now we consider two cases, the first being that $\tau_{\epsilon}<1$. In this case, note that:

[TABLE]

Note that there exists $C_{1}(a_{\epsilon})\geq 0$ such that $-\log\left(\frac{\tau_{\epsilon}}{2^{\frac{1}{\textrm{rank}(\bm{F})}}}\right)<\tau_{\epsilon}^{-a_{\epsilon}/2}+C_{1}(a_{\epsilon})$ $\forall~{}\tau_{\epsilon}<1$. Using this inequality,

[TABLE]

Note the integrand is the sum of the kernel of two gammas: one with shape parameter $A_{\tau_{\epsilon}}+\textrm{rank}(\bm{F})/2-a_{\epsilon}/2+1=a_{\epsilon}/2+n/2-p/2-q/2+\textrm{rank}(\bm{F})/2-1$ and scale parameter $b_{\epsilon}$; the second has shape parameter $A_{\tau_{\epsilon}}+\textrm{rank}(\bm{F})/2+1=a_{\epsilon}+n/2-p/2-q/2-1+\textrm{rank}(\bm{F})/2$ with the same scale parameter. Under our assumptions, the shape parameters are positive and hence the integral is finite.

If $\tau_{\epsilon}>1$, then by the concavity of $\log$, $\log\left(1+\frac{1}{\tau_{\epsilon}^{\textrm{rank}(\bm{F})}}\right)\leq\frac{1}{\tau_{\epsilon}^{\textrm{rank}(\bm{F})}}$ and:

[TABLE]

Note the integrand is the kernel of a gamma with shape parameter $A_{\tau_{\epsilon}}+2=a_{\epsilon}+n/2-p/2-q/2>0$ and scale parameter $b_{\epsilon}$. Thus the integral is finite.

∎

B.2. Proof of 2

Proof.

Because the marginal posterior variance of the regression coefficients for a model of form 2 under the conditions of 1 is of the form $\textrm{Var}(\bm{\beta}^{X})=({\bm{X}}^{T}{\bm{X}})^{-1}\textrm{E}(\sigma_{\epsilon}^{RSR}|\bm{Y})$, it is sufficient to show that $\textrm{E}(\sigma_{\epsilon}^{RSR}|\bm{Y})\leq\textrm{E}(\sigma_{\epsilon}^{NS}|\bm{Y})$. In the following proof, we offer a sufficient condition for this inequality to hold. Note that, as in the paper, $r=\frac{\tau_{s}}{\tau_{\epsilon}}$.

Basic calculations lead to the following:

[TABLE]

Above, we assume a shape-scale parameterization. Note, then, that by properties of the inverse-gamma distribution, we know that:

[TABLE]

By iterated expectations, we can conclude:

[TABLE]

In using iterated expectations, we rely on the assumption that the moments exist. To see this assumption holds note that it can be shown that $\textrm{E}(r|\bm{Y})$ exists by adding a power of $1$ to $\tau_{s}$ in the integral in line Eq. 8 in the proof of 1 and then following the same argument used there. We adopt the notation $A=\textrm{E}({\bm{Y}}^{T}\bm{W}(\bm{I}+r\bm{F})^{-1}{\bm{W}}^{T}\bm{Y})|\bm{Y})$ and note that $A\geq 0$. This expectation can be written as the finite sum of terms of $r$ each of which is bounded between 0 and 1, and therefore the integral converges. To show that 2 holds, we then use the following auxiliary function $M(x)$:

[TABLE]

where $k>0$. We consider $M(x)$ on the interval $[0,\frac{1}{b_{\epsilon}}]$. We note that $M(0)=\textrm{E}(\sigma_{\epsilon}^{NS}|\bm{Y})$. For the case that $k=a_{s}b_{s}$, $M(\frac{1}{b_{s}})=\textrm{E}(\sigma_{\epsilon}^{RSR}|\bm{Y})+\frac{A}{a_{\epsilon}+a_{s}+.5(n-p)-1}$. As previously noted, $A\geq 0$, and conditions A.1[1]-A.1[4] ensure that $a_{\epsilon}+a_{s}+.5(n-p)-1>0$. We now consider the first derivative of $M(x)$:

[TABLE]

Note that $M(x)$ is monotonically decreasing on $[0,\frac{1}{b_{s}}]$ whenever:

[TABLE]

Now, in the case that $k=a_{s}b_{s}$, the result follows. As a final note, we observe that $\textrm{E}(\sigma_{\epsilon}^{NS}|\bm{Y})$ is a non-trivial upper bound (i.e., the moment exists) only in the case that $a_{\epsilon}+.5(n-p)>1$.

We offer a final comment that the result will hold under less restrictive assumptions. We focus on a weaker result here because it is more intuitive and relevant to the settings in which RSR models are fit. Intuitively, the more general result suggests that $\textrm{E}(\sigma_{\epsilon}^{RSR}|\bm{Y})\leq\textrm{E}(\sigma_{\epsilon}^{NS}|\bm{Y})$ whenever a spatial precision matrix $\bm{F}$ is chosen such that components of $\bm{P}_{W}\bm{Y}$ which are large are not smoothed to zero by the prior on $\bm{\delta}$. In other words, if the spatial basis vectors do a good job of explaining the variability in $\bm{Y}$ and the prior on $\bm{\delta}$ does not assume this cannot be true, $\textrm{E}(\sigma_{\epsilon}^{RSR}|\bm{Y})\leq\textrm{E}(\sigma_{\epsilon}^{NS}|\bm{Y})$. To prove this more general result, we combine the above estimates with a lower bound on the term denoted $A$ using Jensen’s inequality. The latter term, in turn, is a function of the decomposition of $\bm{P}_{W}\bm{Y}$ in terms of the eigenvectors and eigenvalues of $\bm{F}$ and the posterior first moment of $r$, which makes the result more technical.

∎

B.3. Proof of 3 and 1

Proof.

The general strategy is as follows: we can directly calculate the form of $h(\beta^{X}|\bm{Y})$. We find an upper bound for $g(\beta^{X}|\bm{Y})$ by bounding an exponential of a non-negative definite matrix by 1, using Lemma 1, and basic inequalities. Doing so allows us to bound the ratio of $\frac{g(\beta^{X}|\bm{Y})}{h(\beta^{X}|\bm{Y})}$ by the ratio of polynomials of the same (even) order. The use of L’Hospital’s Rule then allows us to draw the conclusions about the tail behavior of the two distributions stated in the theorem.

Form of $h(\beta^{X}|\bm{Y})$

In the univariate case, with the help of Mathematica, it is possible to express $h(\beta^{X}|\bm{Y})$ in closed form.

[TABLE]

where $D_{h}=\left(\frac{{\bm{X}}^{T}\bm{X}}{2}\right)^{-a_{\epsilon}-n/2}\frac{\Gamma(d-1/2)}{\Gamma(d)}c^{1/2-d}\sqrt{\pi}$ is defined such that

1. (1)

$c=\frac{2}{{\bm{X}}^{T}\bm{X}}\left[\frac{{\bm{Y}}^{T}\bm{Y}}{2}-\frac{({\bm{X}}^{T}\bm{Y})^{2}}{2{\bm{X}}^{T}\bm{X}}+\frac{1}{b_{\epsilon}}\right]$ 2. (2)

$d=a_{\epsilon}+n/2$

Note further that:

[TABLE]

Form of $g(\beta^{X}|\bm{Y})$

[TABLE]

For the moment we will focus on the numerator of (10).

[TABLE]

We first integrate out $\bm{\delta}$. To do so, with some manipulation, we can recognize a Gaussian kernel and complete the square with:

[TABLE]

After which we are left with:

[TABLE]

Note that $\bm{I}-\bm{W}({\bm{W}}^{T}\bm{W}+\frac{\tau_{s}}{\tau_{\epsilon}}\bm{F})^{-1}{\bm{W}}^{T}$ is non-negative definite. To see this observe that

[TABLE]

${\bm{W}}^{T}\bm{W}$ and $\bm{F}$ commute, and we can use the set of observations in the proof of Lemma 2 to conclude that this is the sum of two non-negative definite matrices and hence non-negative definite. Therefore, we can bound the exponential term with this matrix by 1.

[TABLE]

Again relying on Lemma 1 we can conclude that $\exists~{}C_{q}>0$ such that $|{\bm{W}}^{T}\bm{W}+\frac{\tau_{s}}{\tau_{\epsilon}}\bm{F}|>C_{q}\left(1+\frac{\tau_{s}^{q}}{\tau_{\epsilon}^{q}}\right)$. Therefore, we have the following bound:

[TABLE]

Now let $m(\tau_{\epsilon})$ be the pdf of a gamma distribution with shape parameter $n/2+a_{\epsilon}$ and scale parameter $\left[.5({\bm{\beta}^{X}}^{T}\bm{X}^{T}\bm{X}\beta^{X}-{\beta^{X}}^{T}\bm{X}^{T}\bm{Y}-\bm{Y}^{T}\bm{X}\beta^{X})+\frac{1}{b_{\epsilon}}\right]^{-1}$. Denote $\mu(\beta^{X})=({\beta^{X}}^{T}\bm{X}^{T}\bm{X}\beta^{X}-{\beta^{X}}^{T}\bm{X}^{T}\bm{Y}-\bm{Y}^{T}\bm{X}\beta^{X})$. Substituting in $m(\tau_{\epsilon})$, we find

[TABLE]

Note that $\tau_{\epsilon}^{q}+\tau_{s}^{q}\geq\tau_{s}^{q}$, therefore we can construct another upper bound:

[TABLE]

Considering Ratios of pdfs

Consider the following, if we let $D_{g}^{*}=\iiint f(\beta^{X}|\bm{Y},\bm{\delta},\tau_{\epsilon})f(\bm{\delta}|\tau_{s})f(\tau_{\epsilon})f(\tau_{s})d\bm{\delta}d\tau_{\epsilon}d\tau_{s}d\beta^{X}$

[TABLE]

Because $\mu(\beta^{X})$ is a polynomial of even order, this allows us to conclude the following:

[TABLE]

Implications for Tail Behavior

Note then by the use of L’Hospital’s Rule:

[TABLE]

Note that if $K(C_{q},D_{h})=\frac{(2\pi)^{-n/2}\frac{\Gamma(n/2+a_{\epsilon})}{\Gamma(a_{\epsilon})b_{\epsilon}^{a_{\epsilon}}}D_{h}}{\sqrt{C_{q}}}<D_{g}^{*}$, then:

[TABLE]

∎

B.4. Proof of 4

Proof.

In the framework of model (2), Bayes theorem indicates that the marginal posterior distribution for $\bm{\beta}^{X}$ will be the same for choices of $f,g$ such that:

[TABLE]

To see this note, that by Bayes Theorem:

[TABLE]

So, integrating over $\bm{\delta}$ implies that:

[TABLE]

By the assumption of a priori independent prior distributions, then:

[TABLE]

Hence, substituting $f(\bm{Y}|\bm{\beta}^{X},\tau_{s},\tau_{\epsilon})$ with $g(\bm{Y}|\bm{\beta}^{X},\tau_{s},\tau_{\epsilon})$ will yield an equivalent $f(\bm{\beta}^{X}|\bm{Y})$. Straightforward calculations indicate that $f(\bm{Y}|\bm{\beta}^{X},\tau_{s},\tau_{\epsilon})$ will be proportional to:

[TABLE]

Assuming that $\bm{X}$ and $\bm{Y}$ are fixed, this term will be equivalent for the stated choices of ${\bm{W}}_{i},{\bm{F}}_{i},i=1,2$. To see this it is sufficient to show the following:

1. (1)

$\textrm{rank}({\bm{F}}_{1})=\textrm{rank}({\bm{F}}_{2})$ 2. (2)

$\left|{\bm{W}}_{1}^{T}{\bm{W}}_{1}+\frac{\tau_{s}}{\tau_{\epsilon}}{\bm{F}}_{1}\right|=\left|{\bm{W}}_{2}^{T}{\bm{W}}_{2}+\frac{\tau_{s}}{\tau_{\epsilon}}{\bm{F}}_{2}\right|$ 3. (3)

${\bm{W}}_{1}\left[{\bm{W}}_{1}^{T}{\bm{W}}_{1}+\frac{\tau_{s}}{\tau_{\epsilon}}{\bm{F}}_{1}\right]^{-1}{\bm{W}}_{1}^{T}={\bm{W}}_{2}\left[{\bm{W}}_{2}^{T}{\bm{W}}_{2}+\frac{\tau_{s}}{\tau_{\epsilon}}{\bm{F}}_{2}\right]^{-1}{\bm{W}}_{2}^{T}$

To see $1$, note: $\textrm{rank}({\bm{F}}_{1})=\textrm{rank}({\bm{W}}_{1})=\textrm{rank}({\bm{W}}_{2})=\textrm{rank}({\bm{F}}_{2})$. Conditions 2 and 3 are a consequence of the fact that there exists $\bm{O}\in\textrm{O}(q)$ (i.e., an orthogonal matrix) such that ${\bm{W}}_{1}={\bm{W}}_{2}\bm{O}$ and ${\bm{W}}_{2}={\bm{W}}_{1}\bm{O}^{T}$. This fact results in the following straightforward observation.

[TABLE]

The third condition can be shown with the application of two variants of the Woodbury identity given in lines 156 and 157 of Petersen and Pedersen, (2012).

[TABLE]

∎

Appendix C Additional Summaries and Simulations

C.1. An Overfit NS model for the SAT dataset

In this subsection, we return to the SAT dataset introduced in Section 3.2. We now fit a series of fixed effect models. We start with a model that includes just $\bm{X}^{*}$ (the NS model). We then sequentially add additional synthetic covariates which are columns chosen from an arbitrary basis of ${\bm{P}_{X}}^{\perp}$. To do so, we order (in decreasing value) the columns by their correlation with ${\bm{P}_{X}}^{\perp}\bm{Y}$. We add the synthetic covariates in this order to dominate the effect of the $n-p-q$ in the denominator that would be present in Eq. 5. We do this until there is insufficient data to update the prior on $\tau_{\epsilon}$ for a final model including 43 synthetic covariates and $\bm{X}^{*}$. As before, $a_{s}=.5,~{}b_{s}=2000$.

In Fig. 3, we plot the marginal posterior variance for for $\beta_{0}$,  $\beta_{1}$, and $\beta_{2}$ for each of the series of fixed effect models. As we expect from the discussion in Section 3.2, all models result in approximately the same point estimates of $(590.5,-2.84,.022)$ as additional synthetic covariates and the posterior variances tend to decrease as additional synthetic covariates are added. The relative smoothness of this decrease is due to the fact that we have reduced the impact of the $n-p-q$ term in the denominator of Eq. 5. For RSR models, such a tactic is no longer necessary because irregardless of the dimension of $\bm{W}$, the corresponding term will be $n-p$.

C.2. Bias and MSE Results for Simulations 1 to 3

In the spatial confounding literature, there is an belief that using the ICAR model results in “poor” point estimates of the regression coefficients (See e.g., Reich et al.,, 2006; Hodges and Reich,, 2010; Prates et al.,, 2019). In looking at the performance of the point estimates, most of the work falls into two groups. First, some work considers the point estimates of the ICAR model to be performing poorly when they differ from the point estimates that would have been obtained in the NS model (Reich et al.,, 2006; Hodges and Reich,, 2010). Other work investigates the poor performance with simulation studies. In this work, poor performance is typically assessed by the difference of the point estimates obtained by the ICAR model and the “true” (generating) value of the regression coefficients (Hughes and Haran,, 2013).

Almost all of the previous work has focused on the Gaussian case. Although seldom remarked upon, the magnitude of the difference between the point estimates obtained from an ICAR model and the OLS estimates is often very small. For instance, the Slovenia stomach cancer data set has influenced much of the concern with the poor performance of point estimates of the regression coefficients (See e.g., Reich et al.,, 2006; Hodges and Reich,, 2010; Thaden and Kneib,, 2018; Prates et al.,, 2019). As shown in Section 5, the absolute value of the difference between the point estimate obtained in the ICAR model and the point estimate obtained in the NS model is .119. Similar results are found when assessing the performance of point estimates with simulations (Hughes and Haran,, 2013; Prates et al.,, 2019).

In the following subsections, we investigate the performance of the ICAR point estimates in more detail. We find that the difference between the point estimates obtained in the ICAR model and the generating value tends to have small magnitude. We also find that the point estimates of the ICAR model perform better than that of the NS and RHZ models for count data. Although this latter finding contradicts the expectations of the spatial confounding literature, it is not unprecedented. Hughes and Haran, (2013) observed something similar in the context of a simulated example.

Although not previously noted, the simulation results for Section 4.1-Section 4.3 support 2: in the Gaussian case, the posterior variance for a given regression coefficient for the RHZ model is always less than or equal to that of the posterior variance for the NS model (we considered the results rounded to the 2$nd$ decimal place). As 1 indicates we should expect, the point estimates for regression coefficients for the NS and RHZ model are the same for the Guassian models (we defined same as the absolute difference being less than .01). Because 1 and 2 were only proven for the Gaussian case, we look to the simulation studies to assess whether similar results hold for the Poisson case. Unlike in the Gaussian case, the marginal posterior means of $\beta^{X}$ are not always exactly the same for the NS and RHZ models. However, they are always quite close. The maximum difference is approximately .06 across all simulations. In all but one case, the variance of the marginal posterior distribution for the RHZ model was less than (or equal to) the variance of the marginal posterior distribution for the NS model. This lends support that the findings in 2 may extend in some fort to non-Gaussian cases. The point estimates of the NS and RHZ models are similar enough that bias and mean standard error (MSE) results are the same in almost all cases. Thus, we omit separate categories for the RHZ and NS models in the following discussions of bias and MSE summaries.

C.2.1. Simulation 1

In this section, we consider the bias where we define bias as $\beta^{X}$ - $E[\beta^{X}|\bm{Y}]$. To give a sense of the spread of bias in each setting, we calculate the bias for each of the 1000 simulated data sets. Table 11 then illustrates the 10$th$ and 90$th$ percentiles of the resulting bias measurements. In Table 12, we consider the average MSE, where the average is taken across the 1,000 simulated data sets.

The results indicate that the bias for the ICAR model tends to be larger than the NS/RHZ model when the data is generated from the RHZ model. If the data is generated with an ICAR random effect, then the bias for the ICAR analysis model tends to be less than the the bias for the NS and RHZ analysis models. The results for the average MSE indicate that there are not large differences in the performance of the various models with respect to this measure. This is not surprising given the relatively small magnitude of the bias for all models. The ICAR analysis model performs marginally better with respect to the MSE than the NS/RHZ model when the data is generated with an ICAR random effect.

[TABLE]

[TABLE]

C.2.2. Simulation 2

We again consider the same summaries of bias and MSE as in Section C.2.1. As before, Table 13 shows the NS/RHZ model performs better with respect to bias for $\beta^{X}_{1}$ when the data is generated under an RHZ model. Similarly, the ICAR model performs better when the data is generated with an ICAR random effect. Again, Table 14 demonstrates that the average MSEs are relatively small.

Note, that under a RHZ generating model the 10$th$ percentile of the bias differed slightly for the NS (.195) and RHZ (.194) models. However, for convenience we still rounded both to .20 in Table 13.

[TABLE]

[TABLE]

With one exception, the bias and MSE summaries for $\beta^{X}_{2}$ are quite similar to those observed for $\beta^{X}_{1}$, and therefore omitted. The exception occurs when the generating model is the NS model. In this setting the 10$th$ and 90$th$ percentiles for the bias for the NS/RHZ model are $(-.56,.54)$ and the average MSE is .19. The results for the ICAR model remain similar seen for the settings explored in Table 13 and Table 14.

C.2.3. Simulation 3

In Table 15 and Table 16, we consider the same summaries of bias and MSE as in Section C.2.2 for $\beta^{X}_{1}$. Note that under the RHZ generating model and with rounding, the 10 $th$ percentile in Table 15 is $-.26$ for the NS model and $-.27$ for the RHZ model. For brevity, we just include one case still. Unlike in the Gaussian case, the ICAR model tends to have smaller bias than NS/RHZ models regardless of the generating model. This observation is inconsistent with the general beliefs in the spatial confounding literature (See e.g., Hodges and Reich,, 2010). However, a similar result was seen in an analysis of a single simulated dataset in Hughes and Haran, (2013).

[TABLE]

[TABLE]

The bias and MSE summaries for $\beta^{X}_{2}$ are similar to those observed for $\beta^{X}_{1}$, and therefore omitted.

C.3. Impact of Small Effect Sizes

As the previous subsections indicate, the magnitude of the bias is relatively small for all analysis models. In the context of Section 4.1-Section 4.3, where the effect sizes are all modest (1 or 2), the bias may not seem terribly concerning. That said, it should be noted that given the very narrow credible intervals of the NS and RHZ models for count data, even a small magnitude of bias can help explain the poor coverage rates and elevated Type-S errors.

There are many areas of spatial statistics, such as in environmental epidemiology, where small effect sizes are important. Here, to be clear, we define small effect sizes to be associated with regression coefficients of absolute value $\ll 1$. To investigate the impact of spatial confounding in a setting with small effect sizes, we repeat Simulations 1-3. Now, though, we simulate the response ($\bm{Y}$ or $\bm{Z}$) by replacing $\beta_{0}=1$ and $\beta^{X}=\beta^{X}_{1}=2$ with $\beta_{0}=.1$ and $\beta^{X}=\beta^{X}_{1}=.2$.

For the Gaussian cases, the summaries of bias and MSE typically remain very similar to the results observed in Section C.2.1 and Section C.2.2. Given the smaller effect sizes, now the magnitude of the bias is more of a concern. Despite this, the coverage rates observed in Section 4.1-Section 4.2 remain quite similar. In the following discussions, we omit repeating the coverage, bias, and MSE results for the Gaussian simulations absent a difference from prior results. For both the Gaussian and Poisson models, the results for the Type-S error rates change from Section 4.1-Section 4.3. The impact of the smaller effect size can also be seen in the Bayesian analogue of power, which we define as percentage of time that the 95% credible interval does not contain 0. In Section 4.1-Section 4.3, the power was 100% for all analysis models in every setting. Unsurprisingly, the power suffers with for all models in the context of the smaller effect sizes. We now discuss the results in more detail.

C.4. Simulation 1

Table 17 indicates the RHZ model always has the highest power. This difference is most pronounced when the data is generated with an ICAR random effect. The ICAR model always has the lowest power. Our investigations indicate these differences are explained primarily by the widths of the 95% credible intervals (as opposed to the impact of bias pushing point estimates away from 0).

[TABLE]

C.5. Simulation 2

As in Section C.4, Table 18 indicates that the RHZ model always has the highest power and the ICAR model always has the lowest.

[TABLE]

Table 19 indicates that when the data is generated from a NS model, the Type-S error is much lower than it was in Table 6 for both the RHZ and NS analysis models. This change is primarily driven by a reduction in the bias for $\beta^{X}_{2}$. For the larger effect sizes, Section C.2.2 indicated the 10$th$ and 90$th$ percentiles of the bias were $(-.56,.54)$. Here, the corresponding percentiles are $(-.2,.2)$. The bias for other settings remained largely unchanged.

[TABLE]

C.6. Simulation 3

For the count data, the power depends heavily on how the data was generated. As Table 20 indicates, when the data is generated with an IID random variable, all models have exceptionally poor power. The RHZ and NS models have relatively high power (compared to the Gaussian model) when the data is generated with a spatial random effect.

[TABLE]

For the count data, the effect size appears to have a more direct impact on the coverage rates, bias, and MSE summaries than it did for the Gaussian data. Table 21 indicates that the coverage rates of $\beta^{X}_{1}$ are greatly improved compared to when the generating value of $\beta^{X}_{1}$ was larger (see Table 8).

[TABLE]

The driving factor in this reduction appears to be a reduction in the bias for the NS and RHZ models. Table 22 indicates that the bias for the NS and RHZ models have decreased significantly compared to when the data was generated with larger effect sizes (see Table 15). Table 23 shows that a similar reduction in the average MSE for these models. The summaries for the bias and MSE measures for $\beta^{X}_{2}$ are once again quite similar to those for $\beta^{X}_{1}$ and thus omitted.

[TABLE]

[TABLE]

As with the Gaussian data, a reduction in the bias for the NS/RHZ model corresponds to a reduction in the Type-S error rates, as indicated in Table 24.

[TABLE]