Mitigating Unobserved Spatial Confounding when Estimating the Effect of Supermarket Access on Cardiovascular Disease Deaths

TL;DR

This paper addresses unmeasured spatial confounding in causal inference, proposing a new affine estimator and Bayesian method to accurately estimate the effect of supermarket access on cardiovascular deaths, revealing harmful effects previously masked.

Contribution

It introduces an affine estimator and a Bayesian framework to mitigate unmeasured spatial confounding, improving causal effect estimation in spatial epidemiology.

Findings

Standard spatial models may not reduce bias from unmeasured confounders.

The proposed method uncovers a harmful effect of limited supermarket access on cardiovascular mortality.

Results suggest unobserved confounding can mask true adverse effects in spatial health studies.

Abstract

Confounding by unmeasured spatial variables has received some attention in the spatial statistics and causal inference literatures, but concepts and approaches have remained largely separated. In this paper, we aim to bridge these distinct strands of statistics by considering unmeasured spatial confounding within a causal inference framework, and estimating effects using outcome regression tools popular within the spatial literature. First, we show how using spatially correlated random effects in the outcome model, an approach common among spatial statisticians, does not necessarily mitigate bias due to spatial confounding, a previously published but not universally known result. Motivated by the bias term of commonly-used estimators, we propose an affine estimator which addresses this deficiency. We discuss how unbiased estimation of causal parameters in the presence of unmeasured spatial confounding can only be achieved under an untestable set of assumptions which will often be application-specific. We provide a set of assumptions which describe how the exposure and outcome of interest relate to the unmeasured variables, and we show that this set of assumptions is sufficient for identification of the causal effect based on the observed data when spatial dependencies can be represented by a ring graph. We implement our method using a fully Bayesian approach applicable to any type of outcome variable. This work is motivated by and used to estimate the effect of county-level limited access to supermarkets on the rate of cardiovascular disease deaths in the elderly across the whole continental United States. Even though standard approaches return null or protective effects, our approach uncovers evidence of unobserved spatial confounding, and indicates that limited supermarket access has a harmful effect on cardiovascular mortality.

Tables

Table 1. Table 1: Set of assumptions based on which the causal exposure-response curve derivative is identifiable using observed data and can be estimated using the affine estimator.

Causal Assumptions
Temporal Order	The exposure is temporally precedent to the outcome
SUTVA	No interference between units, no hidden levels of the treatment, $Y_i=Y_i\left(Z_i\right)$
No unmeasured non-spatial confounding	$Z⟂⟂Y\left(z\right)|{𝑪}^m,U,z\in 𝒵$
Positivity∗	$P\left(Z=z|{𝑪}^m,U\right)>0,z\in 𝒵$, which implies that:
Spatial scale restriction	$Z$ has variation at a smaller spatial scale than that of $U$
Structural Assumptions
Outcome Additivity	The exposure and measured covariates do not interact with the unmeasured covariates.
$E\left[𝑼|𝒁\right]$ identification
Normality	$Z,U,\epsilon$ are jointly normal, conditional on ${𝐗}_{\left(-z\right)}$
Cross-Markov	$Z_i⟂⟂{𝑼}_{-i}|U_i,{𝒁}_{-i},{𝐗}_{\left(-z\right)}$, where
	${𝑼}_{-i}={(U_1,U_2,\dots ,U_{i-1},U_{i+1},\dots ,U_n)}^T$, ${𝒁}_{-i}$ defined similarly
Conditional correlation	$\text{Cor}(U_i,Z_i|{𝑼}_{-i},{𝒁}_{-i},{𝐗}_{\left(-z\right)})$ is constant
Precision matrices	The precision matrices of $𝑼|𝒁,{𝐗}_{\left(-z\right)}$ and $𝒁|𝑼,{𝐗}_{\left(-z\right)}$ are of CAR form
Mean Specification	$E\left[𝒀|𝐗,𝑼\right]$ and $E\left[𝒁|{𝐗}_{\left(-z\right)},𝑼\right]$ are correctly specified.

Table 2. Table 2: Simulation results from 500 data sets of size n = 300 𝑛 300 n=300 . The -RS suffix indicates estimators with the restriction φ Z ≤ φ U subscript 𝜑 𝑍 subscript 𝜑 𝑈 \varphi_{Z}\leq\varphi_{U} .

	Mechanism	Estimator	Bias	Std. Err.	RMSE	95% CI Coverage
GM 1	Unconfounded	Non-spatial	$0.00$	$0.14$	$0.14$	$0.62$
		Spatial	$0.01$	$0.12$	$0.12$	$0.93$
		Spatial-RS	$0.02$	$0.12$	$0.12$	$0.93$
		Affine	$0.03$	$0.27$	$0.27$	$0.98$
		Affine-RS	$0.03$	$0.29$	$0.29$	$0.96$
GM 2	Large-scale	Non-spatial	$0.37$	$0.15$	$0.40$	$0.02$
	confounder	Spatial	$0.37$	$0.12$	$0.39$	$0.05$
		Spatial-RS	$0.37$	$0.12$	$0.39$	$0.04$
		Affine	$0.24$	$0.34$	$0.41$	$0.94$
		Affine-RS	$0.14$	$0.30$	$0.33$	$0.95$
GM 3	Large-scale	Non-spatial	$0.33$	$0.14$	$0.36$	$0.03$
	exposure	Spatial	$0.35$	$0.11$	$0.36$	$0.03$
		Spatial-RS	$0.34$	$0.10$	$0.35$	$0.04$
		Affine	$0.33$	$0.25$	$0.42$	$0.91$
		Affine-RS	$0.23$	$0.36$	$0.42$	$0.83$
GM 4	Same scales	Non-spatial	$0.34$	$0.15$	$0.38$	$0.03$
		Spatial	$0.35$	$0.11$	$0.37$	$0.05$
		Spatial-RS	$0.35$	$0.11$	$0.37$	$0.06$
		Affine	$0.29$	$0.27$	$0.39$	$0.93$
		Affine-RS	$0.18$	$0.32$	$0.36$	$0.89$
GM 5	Non-constant	Non-spatial	$0.35$	$0.11$	$0.37$	$0.00$
	conditional	Spatial	$0.37$	$0.06$	$0.37$	$0.00$
	correlation	Spatial-RS	$0.37$	$0.06$	$0.37$	$0.00$
		Affine	$0.38$	$0.27$	$0.47$	$0.57$
		Affine-RS	$0.22$	$0.20$	$0.30$	$0.75$
GM 6	Non-normal	Non-spatial	$0.24$	$0.08$	$0.26$	$0.06$
	joint	Spatial	$0.24$	$0.08$	$0.26$	$0.14$
	distribution	Spatial-RS	$0.24$	$0.08$	$0.26$	$0.15$
		Affine	$0.18$	$0.18$	$0.26$	$0.94$
		Affine-RS	$0.08$	$0.17$	$0.19$	$0.96$

Table 3. Table A.1: CDC WONDER Online Database query, January 9, 2020.

Parameter	Value
Dataset	Underlying Cause of Death, 1999–20017
ICD-10 Codes	I00–I99 (Diseases of the circulatory system)
Ten-Year Age Groups	65–74 years; 75–84 years; 85+ years
Year/Month	2007
Group By	County
Show Totals	True
Show Zero Values	True
Show Suppressed	True
Calculate Rates Per	$100,000$
Rate Options	Default intercensal populations for years 2001–2009 (except Infant Age Groups)

Table 4. Table A.2: Available covariate information. † Log-transformed in analysis

Name	Description	Q1	Median	Q3
PctUrban	Percentage of population in urban areas	$4.2$	$16$	$37$
PctWhite	Percentage of white population	$79$	$93$	$97$
PctBlack	Percentage of black population	$0.21$	$1.1$	$8.4$
PctHisp	Percentage of hispanic population	$0.73$	$1.5$	$4.2$
PctHighSchool	Percentage of population that attended high school	$31$	$36$	$41$
MedianHHInc†	Median household income ($\times 1000$ USD)	$30$	$34$	$39$
PctPoor	Percentage of impoverished population	$9.5$	$13$	$18$
PctFemale	Percentage of female population	$49$	$50$	$51$
PctMovedIn5	Percentage of population having lived in area for less than 5 years	$36$	$40$	$46$
PctOccupied	Percentage of housing units that are occupied	$80$	$87$	$91$
MedianHValue†	Median value of owner occupied housing ($\times 1000$ USD)	$54$	$72$	$95$
PopPerSQM†	Population per square mile	$20$	$58$	$240$
TotPop†	Total county population ($\times 1000$)	$12$	$31$	$94$
smokerate	Percentage of population that smokes	$23$	$26$	$29$

Table 5. Table F.1: Simulation results from 100 data sets of size n = 300 𝑛 300 n=300 . The -RS suffix indicates estimators with the restriction φ Z ≤ φ U subscript 𝜑 𝑍 subscript 𝜑 𝑈 \varphi_{Z}\leq\varphi_{U} .

Estimator	Bias	Std. Err.	RMSE	95% CI Coverage
Non-spatial	$0.36$	$0.17$	$0.40$	$0.03$
Spatial	$0.35$	$0.12$	$0.37$	$0.08$
Spatial-RS	$0.35$	$0.12$	$0.37$	$0.05$
Affine	$0.21$	$0.42$	$0.46$	$0.91$
Affine-RS	$0.11$	$0.26$	$0.28$	$0.95$

Table 6. Table G.1: Simulation results from 1000 data sets of size n = 100 𝑛 100 n=100 . The -RS suffix indicates estimators with the restriction φ Z ≤ φ U subscript 𝜑 𝑍 subscript 𝜑 𝑈 \varphi_{Z}\leq\varphi_{U} .

	Mechanism	Estimator	Bias	Std. Err.	RMSE	95% CI Coverage
GM 1	Independent	Non-spatial	$0.00$	$0.19$	$0.19$	$0.94$
		Spatial	$0.00$	$0.19$	$0.19$	$0.93$
		Spatial-RS	$0.00$	$0.19$	$0.19$	$0.93$
		Affine	$0.01$	$0.34$	$0.34$	$0.98$
		Affine-RS	$0.01$	$0.35$	$0.35$	$0.98$
GM 2	Large-scale	Non-spatial	$0.67$	$0.15$	$0.69$	$0.00$
	confounder	Spatial	$0.64$	$0.15$	$0.65$	$0.01$
		Spatial-RS	$0.63$	$0.14$	$0.64$	$0.01$
		Affine	$0.49$	$0.47$	$0.68$	$0.81$
		Affine-RS	$0.24$	$0.36$	$0.43$	$0.96$
GM 3	Large-scale	Non-spatial	$0.56$	$0.12$	$0.57$	$0.01$
	exposure	Spatial	$0.55$	$0.12$	$0.57$	$0.01$
		Spatial-RS	$0.54$	$0.13$	$0.56$	$0.02$
		Affine	$0.56$	$0.28$	$0.63$	$0.87$
		Affine-RS	$0.30$	$0.59$	$0.66$	$0.88$
GM 4	Same scales	Non-spatial	$0.60$	$0.14$	$0.62$	$0.01$
		Spatial	$0.59$	$0.14$	$0.61$	$0.01$
		Spatial-RS	$0.58$	$0.14$	$0.59$	$0.02$
		Affine	$0.55$	$0.37$	$0.66$	$0.85$
		Affine-RS	$0.31$	$0.47$	$0.56$	$0.92$
GM 5	Non-constant	Non-spatial	$0.37$	$0.09$	$0.38$	$0.03$
	conditional	Spatial	$0.36$	$0.09$	$0.37$	$0.04$
	correlation	Spatial-RS	$0.36$	$0.09$	$0.37$	$0.04$
		Affine	$0.21$	$0.32$	$0.38$	$0.95$
		Affine-RS	$0.09$	$0.27$	$0.28$	$0.98$
GM 6	Non-normal	Non-spatial	$1.15$	$0.50$	$1.25$	$0.01$
	joint	Spatial	$1.03$	$0.45$	$1.13$	$0.01$
	distribution	Spatial-RS	$1.07$	$0.48$	$1.17$	$0.01$
		Affine	$0.74$	$0.63$	$0.97$	$0.84$
		Affine-RS	$0.59$	$0.56$	$0.81$	$0.90$

Equations

$$U_i|{\mathbf{U}}_{-i}\sim \mathcal{N}\left[{\phi }_U\sum _{j\in {\partial }_i}{U}_j/|{\partial }_i|,({\tau }_U|{\partial }_i|{)}^{-1}\right],$$

$$\mathbf{U}\sim \mathcal{N}\left[0,{\tau }_U^{-1}(\mathbf{D}-{\phi }_U\mathbf{W}{)}^{-1}\right],$$

$$Y_i(z)=\eta (z,{\mathbf{C}}_i^m)+g({\mathbf{C}}_i^u)+{\epsilon }_i(z),$$

$$\mathbf{Y}={\mathbf{X}}_{(-z)}{\mathbf{\beta }}_{(-z)}+\mathbf{Z}{\beta }_z+\mathbf{U}+\mathbf{\epsilon },$$

$$\displaystyle\mathrm{E}\big{(}\widehat{\bm{\beta}}|\mathrm{\mathbf{X}}\big{)}$$

$$\tilde{\mathbf{\beta }}=\{{\mathbf{X}}^{⊺}(\mathrm{Var}[\mathbf{Y}|\mathbf{X}]{)}^{-1}\mathbf{X}{\}}^{-1}{\mathbf{X}}^{⊺}(\mathrm{Var}[\mathbf{Y}|\mathbf{X}]{)}^{-1}\mathbf{Y}.$$

$$\mathrm{E}\big{(}\widetilde{\bm{\beta}}|\mathrm{\mathbf{X}}\big{)}=\bm{\beta}+\{\mathrm{\mathbf{X}}^{\intercal}(\mathrm{Var}[\bm{Y}|\bm{Z}])^{-1}\mathrm{\mathbf{X}}\}^{-1}\mathrm{\mathbf{X}}^{\intercal}(\mathrm{Var}[\bm{Y}|\mathrm{\mathbf{X}}])^{-1}\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}].$$

$$\text{widebar}\mathbf{\beta }=\{{\mathbf{X}}^{⊺}(\mathrm{Var}[\mathbf{Y}|\mathbf{X}]{)}^{-1}\mathbf{X}{\}}^{-1}{\mathbf{X}}^{⊺}(\mathrm{Var}[\mathbf{Y}|\mathbf{X}]{)}^{-1}\{\mathbf{Y}-\mathrm{E}[\mathbf{U}|\mathbf{X}]\},$$

$$\left(\begin{array}{c}\mathbf{U}\\ \mathbf{Z}\end{array}\right)\sim \mathcal{N}\left[\left(\begin{array}{c}0\\ {\mathbf{X}}_{(-z)}\mathbf{\gamma }\end{array}\right),{\left(\begin{array}{cc}\mathbf{G}& \mathbf{Q}\\ {\mathbf{Q}}^{⊺}& \mathbf{H}\end{array}\right)}^{-1}\right],$$

$$q_{ij}$$

$$\mathbf{Y}|\mathbf{X}$$

$$\mathbf{Z}|{\mathbf{X}}_{(-z)}$$

$$\mathbf{M}$$

$$\mathbf{C}$$

$$\mathbf{\nu }$$

$$\mathbf{\theta }$$

$$f(\mathbf{Y},\mathbf{Z}|\mathbf{\beta },\mathbf{\gamma })\propto |\mathbf{M}{|}^{-\frac{1}{2}}\exp \left[-\frac{1}{2}(\mathbf{\nu }-\mathbf{C}\mathbf{\theta }{)}^{⊺}{\mathbf{M}}^{-1}(\mathbf{\nu }-\mathbf{C}\mathbf{\theta })\right],$$

$$RL$$

$$\mathbf{Y}=1{\beta }_0+\mathbf{f}(\mathbf{Z})+{\mathbf{X}}_{(-z)}{\mathbf{\beta }}_{(-z)}+\mathbf{U}+\mathbf{\epsilon },$$

$$\text{widebar}f(z)=\sum _{a=1}^A{\beta }_a{z}^a+\sum _{k=1}^K{l}_k|z-{\xi }_k{|}^A,$$

$$\mathbf{X}=\left(\begin{array}{cccc}1& z_1^1& \cdots & z_1^A\\ ⋮& ⋮& \ddots & ⋮\\ 1& z_n^1& \cdots & z_n^A\end{array}\right),\mathbf{L}=\left(\begin{array}{ccc}|z_1-{\xi }_1{|}^A& \cdots & |z_1-{\xi }_K{|}^A\\ ⋮& \ddots & ⋮\\ |z_n-{\xi }_1{|}^A& \cdots & |z_n-{\xi }_K{|}^A\end{array}\right),$$

$${\mathbf{G}}_n={\tau }_U\left(\begin{array}{cccccc}2& -{\phi }_U& & & -{\phi }_U& \\ -{\phi }_U& 2& -{\phi }_U& & & \\ & \ddots & \ddots & \ddots & & \\ & & -{\phi }_U& 2& -{\phi }_U& \\ -{\phi }_U& & & -{\phi }_U& 2& \end{array}\right)$$

$$\underset{n\to \infty }{\lim }\mathrm{Prec}[\mathbf{Z}{]}_{ij}$$

$$\underset{n\to \infty }{\lim }\mathrm{Var}[\mathbf{Y}|\mathbf{Z}{]}_{ij}={\tau }_U^{-1}\frac{1}{2\sqrt{1-{\phi }_U^2}}{\left(\frac{{\phi }_U}{1+\sqrt{1-{\phi }_U^2}}\right)}^{|i-j|},$$

$${\mathrm{E}}_{{\mathbf{X}}_{(-z)},U}\left[\exp \{{\beta }_{z}z+{\beta }_{(-z)}^{⊺}{\mathbf{X}}_{(-z)}+U\}\right]=\exp \{{\beta }_{z}z\}{\mathrm{E}}_{{\mathbf{X}}_{(-z)},U}\left[\exp \{{\beta }_{(-z)}^{⊺}{\mathbf{X}}_{(-z)}+U\}\right].$$

$$\frac{1}{2}\left[|\mathbf{P}|-{\mathbf{u}}^{⊺}\mathbf{G}\mathbf{u}-2{\mathbf{u}}^{⊺}\mathbf{Q}(\mathbf{z}-{\mathbf{X}}_{(-z)}\mathbf{\gamma })-(\mathbf{z}-{\mathbf{X}}_{(-z)}\mathbf{\gamma }{)}^{⊺}\mathbf{H}(\mathbf{z}-{\mathbf{X}}_{(-z)}\mathbf{\gamma })\right]+\log p({\tau }_U,{\tau }_Z,{\phi }_U,{\phi }_Z,\rho ),$$

$$\mathrm{E}[Z_i|{\mathbf{Z}}_{-i},U_i]={\mathbf{x}}_{(-z),i}^{⊺}\mathbf{\gamma }+\frac{{\phi }_Z}{|{\partial }_i|}\sum _{j\in {\partial }_i}(Z_j-{\mathbf{x}}_{(-z),j}^{⊺}\mathbf{\gamma })+\rho \sqrt{\frac{{\tau }_U}{{\tau }_Z}}{U}_i.$$

$$\mu (z)=E[Y(z)]=E\{E[Y(z)|W]\}=E\{E[Y(z)|Z=z,W]\}=E\{E[Y|Z=z,W]\},$$

$$\mathbf{P}=\left(\begin{array}{cc}\mathbf{G}& \mathbf{Q}\\ {\mathbf{Q}}^{⊺}& \mathbf{H}\end{array}\right)$$

$$\mathrm{Var}\left[\left(\begin{array}{c}{U}_i\\ Z_j\end{array}\right)|{\mathbf{U}}_{-i},{\mathbf{Z}}_{-j}\right]={\left(\begin{array}{cc}{p}_{u_i{u}_i}& p_{u_i{z}_j}\\ p_{z_j{u}_i}& p_{z_j{z}_j}\end{array}\right)}^{-1}=\frac{1}{p_{u_i{u}_i}{p}_{z_j{z}_j}-p_{u_i{z}_j}{p}_{z_j{u}_i}}\left(\begin{array}{cc}{p}_{z_j{z}_j}& -p_{z_j{u}_i}\\ -p_{u_i{z}_j}& p_{u_i{u}_i}\end{array}\right),$$

$$\mathrm{Var}[U_i|{\mathbf{U}}_{-i},{\mathbf{Z}}_{-j}]$$

$$\mathrm{Var}[Z_j|{\mathbf{U}}_{-i},{\mathbf{Z}}_{(-j)}]$$

$$\mathrm{Cov}[U_i,Z_j|{\mathbf{U}}_{-i},{\mathbf{Z}}_{-j}]$$

Full text

Keywords: Causal inference; Cardiovascular disease; Food access; Markov random field; Spatial confounding; Unmeasured confounding

1 Introduction

Over 17 million deaths were attributed to cardiovascular disease (CVD) worldwide in 2016, and the prevalence of CVD among people aged 20 or older in the United States that same year was 48% (American Heart Association, 2019). Poor nutrition and high body mass index are major risk factors of CVD (American Heart Association, 2019), and there is evidence that these risk factors are influenced by the availability of nearby supermarkets (Powell et al., 2007; Laraia et al., 2004), which have historically had a higher prevalence of heart-healthy foods (Sallis et al., 1986; Pearce et al., 2008). In a prospective cohort study of individuals who had undergone cardiac catheterization, living in an area (census tract) with low income and poor food access was associated with an increased risk of myocardial infarction or death (Kelli et al., 2019).

Here, our goal is to cast this question within a causal inference framework, and quantify the effect of county-level supermarket availability on the rate of CVD mortality among the elderly (65+ years old) in the United States. For that purpose, we compile a data set including mortality, store, and demographic information for all counties in the continental United States (see also Section 2). Demographic characteristics such as income have been previously associated with CVD risk factors (Kelli et al., 2017), and might be confounders of the effect of supermarket proximity and access on CVD mortality. Even though a number of demographic variables are included in our data set, the relationship of interest is possibly confounded by unobserved or difficult-to-define variables such as regional culture relating to personal vehicles, diet, and general health-consciousness, or state-level support for people with disabilities. Such variables could represent predictors of the exposure which influence where people live, whether or not they own a vehicle, or where businesses choose to locate, and could also represent predictors of the outcome such as how much people exercise, their stress levels, or what food they choose to eat regardless of supermarket availability. Furthermore, these potentially unobserved demographic variables are expected to be spatially correlated, in that nearby counties are expected to have similar values of these unmeasured covariates.

Therefore, we are faced with estimating the causal effect of an exposure on an outcome using spatially-referenced, observational data, and under the threat of unmeasured spatial confounding. To address this challenge, we combine spatial statistics tools and causal inference methodology within a common framework. Even though some attention has been given to causal inference topics in the spatial statistics literature (Paciorek, 2010; Hodges and Reich, 2010; Hughes and Haran, 2013; Hanks et al., 2015), and to spatial topics in the causal inference literature (Verbitsky-Savitz and Raudenbush, 2012; Keele et al., 2015; Papadogeorgou et al., 2018), there is a substantial gap in the intersection of the two fields.

In classic spatial statistics, regression models are often augmented to include spatially correlated random effects in order to “account” or “adjust” for the spatial dependence in the outcome model residuals. However, there is substantial confusion about what exactly these spatial models are capable of accounting for (Hanks et al., 2015). In some settings, spatial mixed models are employed to estimate the relationship between an exposure and outcome without conditioning on spatial information. In this context, Hodges and Reich (2010) and Hughes and Haran (2013) proposed including a spatial random effect that is orthogonal to the exposure of interest. Other times, it is asserted that spatial models adjust for unobserved covariates which have a spatial dependence structure (Congdon, 2013; Lee and Sarran, 2015). Nevertheless, the usual spatial models do not in general eliminate bias due to unobserved confounders, even when the residual variance components are known (Paciorek, 2010). Recently, Keller and Szpiro (2019) discuss the interpretation of estimates from regression models that progressively include spatial basis functions of higher complexity, and they conclude that increasing adjustment might even lead to bias amplification.

From a different perspective, causal inference methodology with spatial data and in the presence of unmeasured spatial confounding has been quite limited, and, to our knowledge, it has been restricted to classic causal inference tools. Within a regression discontinuity framework, Keele et al. (2015) match treated to control units separated by a boundary minimizing geographical distance of matched pairs and balancing observed covariates. Relatedly, Papadogeorgou et al. (2018) proposed matching treated to control units on a criterion incorporating both propensity scores and geographical distance. Although these approaches can, in some cases, address the problem of interest to spatial statisticians, they are not immediately compatible with models commonly used in spatial data analysis which are most often grounded in outcome regression. An exception is found in Thaden and Kneib (2018) where the authors propose a structural equation modeling approach treating the spatial variable as a confounder in a geoadditive model in order to eliminate bias from the unmeasured spatial variable.

In the causal inference literature, unmeasured confounding has been most often dealt with in the realm of sensitivity analysis. Sensitivity analysis is a powerful approach which aims to quantify the robustness of estimated effects to different amounts of unmeasured confounding (Rosenbaum and Rubin, 1983a; Rosenbaum, 2002; Imbens and Rubin, 2015; VanderWeele and Ding, 2017). However, sensitivity analysis does not directly adjust effect estimates for the presence of such confounders, which is what the methodology presented here and the works referenced above aim to achieve.

In this paper, we seek to bolster the bridge between spatial data analysis and causal inference. In order to do so, we consider unmeasured confounding within a formal causal inference framework and examine estimation approaches grounded on models and tools often employed by spatial statisticians. We start by focusing on continuous outcomes and linear models, studying the bias of commonly-used estimators. We propose a model-based approach to estimate the effect of a change in the exposure on an outcome of interest in the presence of unmeasured spatial confounding. Our approach is designed to easily incorporate popular tools in spatial statistics such as hierarchical and linear mixed models. We explain that identification and estimation of the causal parameter in the presence of unmeasured confounding requires untestable assumptions regarding the unmeasured confounders and their relationship with the treatment and outcome of interest. In general, these assumptions have to be application-specific, identification of the causal parameter needs to be evaluated separately for each set of assumptions, and the proposed estimator would have to be adapted to alternative identifying assumptions. For continuous treatments (referred to as exposures), we provide one set of assumptions that is sufficient for identification of the causal exposure-response curve, and one that is not. Importantly, our results illustrate that, when spatial dependencies can be represented using a ring graph, the components in our estimator involving the unmeasured confounder can be identified based solely on spatial dependencies in the observed data. Within the context of our study, we extend our approach to non-continuous outcomes and generalized linear mixed models employing a fully-Bayesian approach, and we carefully discuss the plausibility of the causal assumptions. While our development is in the context of areal data, refinements in the context of point-referenced data are possible and are discussed where applicable.

In Section 2 we present a detailed description of our data set and present preliminary analyses using non-spatial and commonly-used spatial regression models that yield suspect results. In Section 3, we define the causal estimand in terms of potential outcomes for continuous exposures and discuss commonly invoked identifiability assumptions when the observed covariates include a sufficient confounding adjustment set. The proposed methodology is introduced in Section 4 within the context of linear models. There, we re-derive the result by Paciorek (2010) stating that commonly-used spatial regression models do not recover the estimands of interest in the presence of unmeasured spatial confounding. Motivated by the bias of commonly used estimators, we propose the affine estimator, and we provide a set of assumptions relating the exposure and outcome to the unmeasured variables based on which the causal parameter is identifiable from observed data. In Section 5, we extend the affine estimator in the context of non-linear models and suggest using a Bayesian approach. The estimator is compared to the currently-used estimators under various generative mechanisms via simulation in Section 6. In Section 7 we discuss the plausibility of our assumptions within the context of our study, and we use the affine estimator to estimate the county-level effect of poor supermarket availability on CVD mortality. Our study illustrates the potential of the affine estimator in mitigating bias from unmeasured spatial confounders, returning effect estimates that are qualitatively different from the ones in Section 2, and more in line with subject-matter knowledge. We conclude with a discussion in Section 8.

2 County-level supermarket availability and CVD mortality

We compile a data set including mortality, store availability, demographic and behavioral data on 3,093 out of 3,109 counties and county equivalents in the 48 contiguous states and the District of Columbia. For each areal unit (county, or county equivalent), supermarket availability is defined as the proportion of housing units during 2006 that are more than 1 mile from the nearest supermarket or large grocery store and do not have a car, obtained from the Food Environment Atlas, June 2012 release (United States Department of Agriculture Economic Research Service, 2012). County-level population and cardiovascular disease mortality count (ICD-10 codes I00–I99) during 2007 for residents 65 years old and over were obtained from the United States Centers for Disease Control and Prevention (CDC) WONDER query system (Centers for Disease Control and Prevention, 2018). Due to privacy constraints, county death counts below 10 are censored by CDC WONDER. Figure 1 displays the exposure and crude outcome relative risk, without covariate-assisted estimation or smoothing. Demographic information was acquired from the 2000 Census and includes, among others, information on urbanicity, poverty, and population mobility. Covariate information also includes estimates of 2006 smoking rates derived from CDC Behavioral Risk Factor Surveillance System data (Dwyer-Lindgren et al., 2014). Appendix A provides additional information on the data sources, data collection, and data processing pipeline, including links to the publicly-available data sets, and a table including names and descriptive statistics of available covariates.

At this point, we consider two common analyses investigating the relationship between limited county-level access to supermarkets on cardiovascular deaths in the elderly. For the first analysis, we model CVD mortality counts as a Poisson-distributed outcome using a log link with the exposure and all covariates listed in LABEL:app_tab:table1 as predictors. Internal standardization was implemented by using log expected death count as the offset: the population age 65+ in each county was multiplied by the overall crude rate in the same age range. All covariates were standardized, and the model was fit within the Bayesian paradigm under Gaussian priors with mean 0 and standard deviation 10 on regression coefficients. Results are presented on original scales unless otherwise noted. Samples from the posterior distribution were obtained via a Gibbs sampler, with a Metropolis-Hastings block update for all regression coefficients. The Gibbs sampler was run for 10,000 iterations after a 1,000-iteration burn-in. Censored outcomes were imputed subject to the known privacy constraint.

The second analysis we implemented is a common analysis method for areal spatial data. We included a spatially correlated random effect for county, $\bm{U}$, in the linear predictor, and assumed it follows a conditional autoregressive (CAR) structure (Besag, 1974), i.e.,

[TABLE]

where $\partial_{i}$ and $|\partial_{i}|$ are the set and number of $i$’s neighbors, respectively. A multivariate normal representation of the distribution of $\bm{U}$ is then available as

[TABLE]

where $w_{ij}=1$ if $i$ and $j$ are neighbors, and 0 otherwise, and $\mathrm{\mathbf{D}}$ is diagonal with entries $|\partial_{i}|$ (Banerjee et al., 2004). A uniform prior on $(-1,1)$ was used for $\varphi_{U}$ and a gamma prior with shape and rate parameters equal to 5 was used for $\tau_{U}$, jointly restricted to require the precision matrix of $\bm{U}$ to be positive definite.

Based on these two models, we acquired what the estimated effect of a 1 percentage point increase in households with poor supermarket access on cardiovascular mortality would be if each model was specified correctly and was sufficiently adjusted for confounding variables. The non-spatial estimate (analysis 1) indicated that increasing poor supermarket availability is protective of CVD deaths, with an estimated relative risk of 0.968 (95% CI 0.962 to 0.973). A protective effect of poor supermarket access is not consistent with either theoretical or data-driven understanding of the phenomenon, enforcing our belief that the result is at least partially due to unobserved, or poorly adjusted-for confounders. The estimate from the spatial model (analysis 2) is effectively null with an estimated relative risk of 0.999 (95% CI 0.988 to 1.011). Although the spatial estimate differs from the non-spatial estimate in both location and credible interval width, we show in the next section that the spatial estimate does not necessarily mitigate bias from unobserved confounders.

3 Causal estimands and classic identifiability assumptions

Broadly speaking, the causal inference literature places substantial emphasis on defining target quantities of interest, referred to as estimands, and determining sufficient assumptions under which such estimands (which include unobservable quantities) are identifiable based on the observed data. We begin by defining estimands of interest following the potential outcome framework formalized by Rubin (1974) and extended to continuous exposures by Hirano and Imbens (2004). A necessary condition for an exposure $Z$ to have an effect on an outcome $Y$ is that $Z$ is temporally precedent. We make the stable unit treatment value assumption (SUTVA, Rubin (1980)) which states that there is a single version of each treatment level and there is no interference between units. Based on SUTVA, we can use $Y_{i}(z)$ to represent the value that would have been observed at location $i$ had it received exposure $z\in\mathcal{Z}$, where $\mathcal{Z}$ includes all possible values of the continuous $Z$, and $i=1,2,\dots,n$. Then, $Y_{i}(z)$ is the potential outcome for location $i$ at exposure level $z$, and unit $i$’s observed outcome $Y_{i}$ is the potential outcome for the observed level of the treatment, $Y_{i}=Y_{i}(Z_{i})$.

The most common estimands for continuous treatments are the population average exposure-response curve (PAERC) $\mu(z)=\mathrm{E}[Y_{i}(z)],z\in\mathcal{Z}$, and the expected rate of change in the outcome for an infintensimal change in the exposure around $z$, $\mu^{\prime}(z)$. Since $\mu(z)$ represents the average outcome value over the whole population had everyone experienced exposure $z$, it is clear that $\mu(z)$ includes unobserved quantities, and assumptions need to be made to ensure identifiability and to estimate it using data. The positivity and no unmeasured confounding assumptions (referred to together as the ignorability assumption) form a sufficient set of assumptions for identifiability of $\mu(z)$. Positivity states that all units can experience any $z\in\mathcal{Z}$, and the no unmeasured confounding assumption states that there exist measured covariates $\bm{C}$ which satisfy that, conditional on $\bm{C}$, the observed exposure $Z$ is independent of the potential outcomes $Y(z)$, denoted as $Z\perp\!\!\!\perp Y(z)|\bm{C},z\in\mathcal{Z}$. (See Appendix B for a discussion on identifiability of $\mu(z)$ based on these assumptions.)

Confounders $\bm{C}$ are generally thought of as temporally precedent to the exposure $Z$ and as common predictors of $Z$ and $Y$, as shown in Figure 2. Since temporal order of variables matters in drawing causal conclusions, observed data are conceived as if generated in the following order: $[\bm{C}]$, $[Z|\bm{C}]$ and $[Y|Z,\bm{C}]$. If the identifiability conditions of positivity and no unmeasured confounding are met in the observed data, estimation can proceed by imitating the data generating mechanism for the exposure $Z|\bm{C}$, known as the propensity score (Rosenbaum and Rubin, 1983b), the data generating mechanism for the outcome $Y|Z,\bm{C}$ (e.g., Hill, 2011; Hahn et al., 2018), or both (Robins and Rotnitzky, 1995; Hahn, 1998; Zigler and Dominici, 2014; Zhou et al., 2019). In order to adhere to common approaches of spatial statistics which emphasize analytical models imitating the outcome generative model, our primary focus is on modeling $Y|Z,\bm{C}$.

Even though confounding adjustment is necessary to draw causal conclusions, $\bm{C}$ might include components that are not measured, hence violating the no unmeasured confounding assumption. Denote $\bm{C}=(\bm{C}^{m},\bm{C}^{u})$ representing the measured and unmeasured components, respectively. At this point, we assume that at least some of the variables in $\bm{C}^{u}$ are spatial and refer to Section 4.5 for a further discussion on this requirement. We refer to a variable as “spatial” if the correlation of the variable for two observations depends on their geographic locations. For areal data like the ones in our study, this could refer to adjacency of counties. For point referenced data, it could refer to the geographical distance of two points.

In this section and the next, we discuss unobserved spatial confounders in the case of continuous outcomes and linear models. Focusing on this setting allows for straightforward application of theory from least squares estimation of regression coefficients and restricted maximum likelihood estimation of variance parameters. We return to non-Gaussian outcomes for the simulation study and data analysis, where we employ a fully Bayesian approach. Assume here that potential outcomes arise in the following manner:

[TABLE]

for some function $\eta$, and $\varepsilon_{i}(z)$ a mean zero random variable and independent of $\bm{C}$. In Equation 3, $\bm{C}^{u}$ is assumed to not interact with $Z$ and $\bm{C}^{m}$. We denote $U=g(\bm{C}^{u})$, representing the cumulative contribution of all unobserved covariates. Since at least some components of $\bm{C}^{u}$ are spatial, $U$ also has a spatial correlation structure. Without loss of generality, we may assume $\mathrm{E}[g(\bm{C}^{U})]=0$ by absorbing any non-zero mean into the the intercept in $\eta(z,\bm{C}^{m}$).

4 The affine estimator in linear models: Addressing omitted variable bias of classic estimators

In this section, we discuss how the classic approaches to estimation within the spatial statistics literature are biased for estimating causal quantities in the presence of unmeasured confounders, and in the context of linear regression. The bias results derived below are in line with results in Paciorek (2010), and they motivate the affine estimator, which is designed to explicitly remove the bias of the existing estimators. An extension to non-linear models within the Bayesian framework is presented in Section 5.

4.1 Omitted variable bias of ordinary and generalized least squares estimators

Let $\mathrm{\mathbf{X}}=(\bm{1}\ |\ \bm{Z}\ |\ \bm{C}^{m})$ be the design matrix containing an intercept, the exposure $\bm{Z}$, and measured covariates, and let $\mathrm{\mathbf{X}}_{(-z)}$ be the design matrix including an intercept and measured confounders, but not the exposure $\bm{Z}$. For simplicity of presentation, we assume that the causal exposure response curve is linear, $\eta(z,\bm{C}_{i}^{m})=\bm{\beta}_{(-z)}^{\intercal}\bm{x}_{i,(-z)}+\beta_{z}z$, which will be relaxed in Section 4.4.2. This simplification implies that $\mu^{\prime}(z)=\beta_{z}$, corresponding to the usual linear regression coefficient targeted in the spatial statistics literature. Using vector notation let $\bm{Y}=(Y_{1},Y_{2},\dots,Y_{n})^{\intercal}$, with $\bm{Z},\bm{U},\bm{\varepsilon}$ defined analogously. Then, if all of $\bm{Y},\bm{Z},\bm{U}$, and $\mathrm{\mathbf{X}}_{(-z)}$ were observed, estimation of the regression model

[TABLE]

would lead to consistent estimation of the causal effect through estimation of $\beta_{z}$.

However, the above model cannot be directly used in settings where $\bm{U}$ is not measured. In Equation 4, $\bm{Z}$ and $\bm{U}$ are correlated, but $\bm{\varepsilon}$ is independent of $(\bm{Z},\bm{U})$. Thus $\bm{U}$ and $\bm{\varepsilon}$ is a partition of the variability in $\bm{Y}$ not due to $\bm{Z}$ into one component ($\bm{\varepsilon}$) which is independent of $\bm{Z}$ and one ($\bm{U}$) which is not. If $\bm{U}$ is correlated with $\bm{Z}$ and is omitted from the outcome regression, the ordinary least squares (OLS) estimator of $\bm{\beta}=(\bm{\beta}_{(-z)}^{\intercal},\beta_{z})$, $\widehat{\bm{\beta}}$, will be biased. This is evident by examining the conditional expectation of $\widehat{\bm{\beta}}$:

[TABLE]

Considering \mathrm{E}\big{(}\widehat{\bm{\beta}}\big{)}=\mathrm{E}\big{[}\mathrm{E}\big{(}\widehat{\bm{\beta}}|\mathrm{\mathbf{X}}\big{)}\big{]}, we see that $\widehat{\bm{\beta}}$ will be biased for $\bm{\beta}$ since the second term will, in general, be non-zero for correlated $\bm{U},\bm{Z}$.

When $\bm{U}$ is omitted from the regression model, the component of $\bm{U}$ not attributed to $\mathrm{\mathbf{X}}$ will be incorporated in the residuals. Since $\bm{U}$ is spatially structured, residuals of the regression of $\bm{Y}$ on $\mathrm{\mathbf{X}}$ will also be spatially correlated. In an effort to account for residual spatial correlation, spatial linear mixed models are often adopted. Typically, such models represent mechanisms similar in form to Equation 4, but in which all right-hand-side variables are assumed to be independent, and some assumptions are made about the form of $\mathrm{Var}[\bm{U}]=\mathrm{Var}[\bm{U}|\mathrm{\mathbf{X}}]$. These models aim to explain the spatial correlation in the residuals and they are often effective at improving efficiency. However, they do not necessarily alleviate the omitted variable bias (Paciorek, 2010). If $\mathrm{Var}[\bm{Y}|\mathrm{\mathbf{X}}]$ (which depends on $\mathrm{Var}[\bm{U}]$) is known, the generalized least squares (GLS) estimator of $\bm{\beta}$ is

[TABLE]

with conditional expected value

[TABLE]

Therefore, even if $\mathrm{Var}[\bm{Y}|\mathrm{\mathbf{X}}]$ is known, $\widetilde{\bm{\beta}}$ remains biased. This result indicates that including a spatial random effect in the regression model does not necessarily mitigate or eliminate bias arising from unmeasured spatial confounders.

4.2 The affine estimator to account for omitted spatial variables

The results presented above establish that spatial correlation in the outcome model residuals might arise due to spatial predictors of $\bm{Y}$, and commonly used approaches to estimate $\beta_{1}=\mu^{\prime}(z)$ are biased in the presence of unmeasured confounding by a spatial variable $\bm{U}$. It is now clear that mitigating bias from unmeasured spatial variables cannot be achieved based solely on an outcome regression model without making additional assumptions, nor by harvesting spatial information found solely in the outcome model residuals.

An investigation of the formulas in Equation 5 and Equation 7 shows that bias of both least squares estimators arises from the non-zero correlation between $\bm{U}$ and $\bm{Z}$, leading to a non-zero $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$. Inspired by the form of the bias, we propose an estimator that includes a component that depends on $U$. Consider the affine estimator:

[TABLE]

which replaces $\bm{Y}$ by $\bm{Y}-\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ and is unbiased if $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ is known, or more practically, consistent if $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ is consistently estimated.

Since $\bm{U}$ is unmeasured, direct modeling of $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ based on traditional estimation methods is not possible. Hence, identifiability of this component and our ability to calculate the affine estimator require additional assumptions. In Section 4.3, we provide a set of assumptions based on Gaussian Markov random field theory which pertain to the joint distribution of $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$. Based on these assumptions, we discuss an approach to calculating the affine estimator in the context of restricted maximum likelihood in Section 4.4. Then, in Section 4.5 we show that these assumptions form a sufficient set for identification of the components of $\bm{U}$ on which the affine estimator is based, and therefore the estimation procedure is sound. The identifiability results illustrate that identification is achieved by exploiting the spatial correlation structure in the exposure and outcome model residuals which is driven by the unmeasured spatial variable.

4.3 A sufficient set of assumptions

In this section, we present a set of assumptions which pertain to both the spatial and causal aspect of the affine estimator, and are summarized in Table 1. We proceed with these for ease of illustration, and because they seem plausible within our study setting (see Section 7.1), but note that different or weaker assumptions for identification of $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ are likely possible (see Section 4.5 and Section 8).

4.3.1 A Gaussian Markov random field construction of the joint distribution

In viewing the model from a spatial perspective and to better align to the spatial modeling literature, we assume that the marginal distributions of $\bm{U}$ and $\bm{\varepsilon}$ are Gaussian with mean zero, independent of the measured covariates $\mathrm{\mathbf{X}}_{(-z)}$, and that $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$ is multivariate normal. We see the assumption that $\bm{U}$ is independent of $\mathrm{\mathbf{X}}_{(-z)}$ as without loss of generality, since the same procedure could be alternatively followed for $\bm{U}-\mathbb{P}(\bm{U}|\mathrm{\mathbf{X}}_{(-z)})$, where $\mathbb{P}(\bm{U}|\mathrm{\mathbf{X}}_{(-z)})$ is the projection of $\bm{U}$ on the column space of $\mathrm{\mathbf{X}}_{(-z)}$. This is also supported by results from simulated scenarios under which $U$ and $\mathrm{\mathbf{X}}_{(-z)}$ are correlated (see Section 6). Further, we recognize that joint normality might be a strong assumption and we discuss an approach to relaxing it in Section 8. We make the following assumptions about the joint distribution of $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$.

Cross-Markov property: $p(Z_{i}|\bm{Z}_{-i},\bm{U},\mathrm{\mathbf{X}}_{(-z)})=p(Z_{i}|\bm{Z}_{-i},U_{i},\mathrm{\mathbf{X}}_{(-z)})$, 2. 2.

Constant conditional correlation: $\mathrm{Cor}(U_{i},Z_{i}|\bm{U}_{-i},\bm{Z}_{-i},\mathrm{\mathbf{X}}_{(-z)})=\rho$.

The first assumption states that, conditional on measured covariates and the values of $Z$ at all other locations, $Z_{i}$ depends on $\bm{U}$ only through its value at location $i$, $U_{i}$. Thus, it accommodates correlation between nearby treatments, but it does not allow $U_{j}$ to directly affect the value of $Z_{i}$ for $i\neq j$. The second assumption states that the conditional correlation between $U_{i}$ and $Z_{i}$ does not vary by location. In the joint distribution of $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$, these assumptions can be incorporated in the precision matrix (see Appendix C for derivations). Specifically, if

[TABLE]

the cross-Markov assumption is equivalent to diagonal $\mathrm{\mathbf{Q}}$, and along with the constant conditional correlation assumption they imply that

[TABLE]

Given the above framework, the joint model for $\bm{U}$ and $\bm{Z}$ is completed by specifying $\mathrm{\mathbf{G}}$ and $\mathrm{\mathbf{H}}$, the precision matrices of $\bm{U}|(\bm{Z},\mathrm{\mathbf{X}}_{(-z)})$ and $\bm{Z}|(\bm{U},\mathrm{\mathbf{X}}_{(-z)})$ respectively, up to some unknown parameters that will be estimated from the data. For areal data like the ones in Section 7, we adopt conditional autoregressive structures (CAR; Besag (1974)) for $\mathrm{\mathbf{G}}$ and $\mathrm{\mathbf{H}}$, a common assumption in standard spatial analysis models. Then, the precision matrices $\mathrm{\mathbf{G}}$ and $\mathrm{\mathbf{H}}$ are assumed to share the same neighborhood structure which is encoded in the matrices $\mathrm{\mathbf{D}}$ and $\mathrm{\mathbf{W}}$ of Equation 2, but are allowed to differ by their precision and spatial dependence parameters $(\tau_{U},\varphi_{U})$ and $(\tau_{Z},\varphi_{Z})$. Based on Equation 10, the assumed CAR structure leads to $\mathrm{\mathbf{Q}}=-\rho\sqrt{\tau_{U}\tau_{Z}}\ \mathrm{\mathbf{D}}$.

In the analysis of point-referenced data, the precision matrices of $\bm{U}|(\bm{Z},\mathrm{\mathbf{X}}_{(-z)})$ and $\bm{Z}|(\bm{U},\mathrm{\mathbf{X}}_{(-z)})$ can be specified based on a correlation function decaying in geographical distance. In either case, since $\bm{U}$ includes all unmeasured spatial variables $\bm{C}^{u}$, the correct specification of its precision matrix $\mathrm{\mathbf{G}}$ becomes harder for a larger number of unmeasured spatial covariates. We note again here that, since $\bm{U}$ is unmeasured, estimating the components of the joint distribution in Equation 9 that contribute to the affine estimator cannot be based on traditional modeling approaches, and instead is based on harvesting information from the spatial structure in exposure and outcome model residuals (as we see in Section 4.4 and Section 4.5).

4.3.2 Spatial scale restriction for the unmeasured spatial confounder

In order to draw causal conclusions using our approach, the no-unmeasured and positivity assumptions still need to hold, conditional on the measured covariates $\bm{C}^{m}$, and the unmeasured $U$. The assumption of positivity implies that estimation of the causal effect of $Z$ on $Y$ in the presence of $U$ is only possible if there is variability in $Z$ within levels of $U$. If the spatial scale of $U$ is smaller than that of $Z$, the positivity assumption will be violated, since, loosely speaking, there may be “strata” of $U$ within which only one value of $Z$ is possible. Therefore, from a causal perspective, we assume that the spatial scales of the exposure and spatial confounder do not violate positivity of the treatment assignment conditional on the unmeasured spatial confounder and the measured covariates. In Section 4.6, we also discuss how the spatial scale restriction is also useful in settings where spatial variables mediate the effect of interest.

The spatial scale restriction has been studied from a spatial perspective. Paciorek (2010) shows that the bias and variance of spatial model estimators depend on the relative spatial scales of the exposure and the residual including the confounder, $\varepsilon+U$. Paciorek (2010) recommends only fitting spatial models when there is exposure variation on a spatial scale smaller than that of the unmeasured confounder, essentially ensuring positivity. From a spatial perspective, the spatial scale restriction ensures that we do not mistakenly attribute all spatial variability of the outcome residuals to the unmeasured spatial confounder when it is truly due to the exposure.

The spatial scale restriction can be enforced through the precision matrices $\mathrm{\mathbf{G}},\mathrm{\mathbf{H}}$ in Equation 9. For geostatistical data, the spatial scale of dependence is often an explicit modeling parameter, as in Paciorek (2010). For a conditional autoregressive model of areal data, the autocorrelation parameters $\varphi_{Z}$ and $\varphi_{U}$ do not have strict interpretations as spatial scale parameters, though the restriction $\varphi_{Z}<\varphi_{U}$ plays a similar role in identifying variance parameters.

4.4 The affine estimator within a restricted likelihood framework

In this section we describe estimation within a restricted likelihood framework. We do so because it allows for straightforward illustration of how model components correspond to components in the bias results of Section 4 and the affine estimator in Equation 8. Further, it allows us to easily discuss identifiability of the model parameters in Section 4.5. In Section 5, we describe a fully-Bayesian approach to estimation which is applicable for linear and non-linear models, and which we follow for the remainder of this paper.

4.4.1 Linear effect estimator

We start by assuming the linear structure in Equation 4. Using the conditional distribution $\bm{U}|\mathrm{\mathbf{X}}$ acquired from Equation 9, the joint model for the observed data (integrating $\bm{U}|\mathrm{\mathbf{X}}$ out) can be factored as

[TABLE]

where $\mathrm{\mathbf{R}}^{-1}=\mathrm{Cov}(\bm{\varepsilon})$ (see Appendix D.1 for the derivation). From Equation 11, we see that the likelihood depends on $\bm{U}$ through the components of the precision matrix in Equation 9. Note that, even though our focus is in estimating parameters of the outcome model ($\bm{\beta}$), an exposure model is also adopted to provide information on the spatial structure of $\bm{U}$. (This is related to many settings in causal inference where incorporating information from the exposure model improves estimation of causal effects (e.g. Wilson and Reich, 2014; Belloni et al., 2014; Antonelli et al., 2019).)

Following a common approach to estimation for mixed models, variance parameters are estimated based on the restricted likelihood derived from Equation 11, and the estimates are used to calculate the bias-adjusted affine estimator $\widebar{\bm{\beta}}$ in Equation 8. Defining

[TABLE]

we can write the joint distribution of $(\bm{Y},\bm{Z})$ as

[TABLE]

and the restricted likelihood as

[TABLE]

If $\widehat{\mathrm{\mathbf{M}}}$, $\widehat{\mathrm{\mathbf{C}}}$, and $\widehat{\bm{\nu}}$ are maximizers of the restricted likelihood in (14), we calculate $(\widebar{\bm{\beta}},\widebar{\bm{\gamma}})=(\widehat{\mathrm{\mathbf{C}}}^{\intercal}\widehat{\mathrm{\mathbf{M}}}^{-1}\widehat{\mathrm{\mathbf{C}}})^{-1}\widehat{\mathrm{\mathbf{C}}}^{\intercal}\widehat{\mathrm{\mathbf{M}}}^{-1}\widehat{\bm{\nu}}$. (Readers interested in the REML approach can find additional information in Appendix D.)

The restricted likelihood formulation allows us to make illuminating connections between our approach, the mixed effects models often used in spatial statistics, and the bias results of existing approaches in Section 4. If $\rho=0$, the matrix $\mathrm{\mathbf{Q}}$ is zero, the model in Equation 9 reduces to the case where $\bm{U},\bm{Z}$ are independent, and the restricted likelihood estimation method would lead to the estimator $\widetilde{\bm{\beta}}$. A non-zero correlation $\rho$ between $\bm{U}$ and $\bm{Z}$ leads to the inclusion of the $-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}$ component in the coefficient of $\bm{Z}$, corresponding to the bias correction term $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]=-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}(\bm{Z}-\mathrm{\mathbf{X}}_{(-z)}\bm{\gamma})$.

4.4.2 Semi-parametric effect estimator

To better accommodate continuous exposures, we can flexibly model the exposure-response relationship using penalized regression splines. Penalized regression splines may be represented as linear mixed models (Ruppert et al., 2003), allowing for the linear effect assumption in Equation 4 to be relaxed to

[TABLE]

where $\bm{f}(\bm{Z})=(f(Z_{1}),\ldots,f(Z_{n}))^{\intercal}$ and $f$ being a smooth function of $Z$. Our chosen radial basis penalized spline model for $f$ may then be written as

[TABLE]

where $A$ is the degree of the spline and the $\xi_{k}$ are pre-specified knots. Letting

[TABLE]

and $\mathrm{\mathbf{V}}=\psi^{-1}\mathrm{\mathbf{L}}\mathrm{\mathbf{L}}^{\intercal}+\mathrm{\mathbf{G}}^{-1}+\mathrm{\mathbf{R}}^{-1}$, where $\psi>0$ is a roughness penalty, the restricted likelihood is as in Equation 14 with updated components $\mathrm{\mathbf{C}}$ and $\mathrm{\mathbf{M}}$. Letting $\mathrm{\mathbf{T}}=(\mathrm{\mathbf{X}}\;\mathrm{\mathbf{L}})$ and $\mathrm{\mathbf{A}}$ be a diagonal matrix with the first $A+1$ elements equal to [math] and the rest equal to $1$ (corresponding to penalization of the $\bm{\beta}$ and $\bm{l}$, respectively), the estimate of $\bm{\theta}=(\beta_{0},\beta_{1},\ldots,\beta_{A},l_{1},l_{2},\ldots,l_{K})$ is $\widebar{\bm{\theta}}=\left(\mathrm{\mathbf{T}}^{\intercal}\widebar{\mathrm{\mathbf{V}}}^{-1}\mathrm{\mathbf{T}}+\psi\mathrm{\mathbf{A}}\right)^{-1}\mathrm{\mathbf{T}}^{\intercal}\widebar{\mathrm{\mathbf{V}}}^{-1}\left(\bm{Y}-\widebar{\mathrm{\mathbf{B}}}\bm{Z}\right)$.

4.5 Identifiability of parameters

It is evident from the form of the affine estimator in Equation 8 and the restricted likelihood in Equation 14 that calculating the estimator depends on being able to estimate components of the relationship of the unmeasured confounder with the exposure and outcome of interest. Therefore, it is natural to wonder whether these components can be identified, and if so, which assumptions are key in driving identifiability and which can be relaxed. We provide two examples: one in which identifiability can be proved analytically and one in which model components are not identifiable. The details of how identifiability is achieved or lost in these examples illuminate the key assumptions and their roles.

4.5.1 Identifiability of model components for a ring graph

The first example is one in which identifiability can be analytically proved. We consider the setting without measured variables since including them complicates the notation without providing any additional insight, and results extend trivially. We also assume that the spatial dependence can be represented in a ring graph of $n$ locations depicted in Figure 3, with CAR specifications for $\bm{U}|\bm{Z}$ and $\bm{Z}|\bm{U}$. This spatial structure and model yields the precision matrix

[TABLE]

for $\bm{U}|\bm{Z}$ and similar for $\mathrm{\mathbf{H}}_{n}$ for $\bm{Z}|\bm{U}$. For simplicity, we assume that $\mathrm{E}[\bm{Z}]=\bm{0}$ (marginally over $\bm{U}$) and that $\mathrm{E}[\bm{Y}|\bm{Z},\bm{U}]=\beta_{Z}\bm{Z}+\bm{U}$. We present two results: the first discussing parameter identifiability based on $\bm{Z}$ alone (Theorem 1), and the second based on $(\bm{Y},\bm{Z})$ (Theorem 2). Importantly, these identifiability results allow for $\bm{U}$ to be completely unmeasured.

Theorem 1**.**

In the scenario defined in this section, it can be determined whether or not $\rho\varphi_{U}=0$ by observing $\bm{Z}$. Further, if $\rho\varphi_{U}\neq 0$, then $(\tau_{Z},\varphi_{Z},\varphi_{U},|\rho|)$ is also identifiable by observing $\bm{Z}$.

Proof.

The proof relies on a few matrix lemmas which are stated and given proof outlines in Appendix E.1. We have that $\mathrm{Prec}[\bm{Z}]=\tau_{Z}\left[\mathrm{\mathbf{H}}_{n}-4\rho\mathrm{\mathbf{G}}_{n}^{-1}\right]$ (marginally over $\bm{U}$), and

[TABLE]

First, note that for any lag $l$, the number of pairs of locations with $|i-j|=l$ grows linearly with $n$. It can be determined whether or not $\rho\varphi_{U}=0$ because $\lim_{n\to\infty}\mathrm{Prec}[\bm{Z}]_{ij}=0$ for all $(i,j)$ such that $|i-j|>1$ if and only if $\rho\varphi_{U}=0$. If $\rho\varphi_{U}\neq 0$, then for $|i-j|>1$ and $|i^{\prime}-j^{\prime}|=|i-j|+1$ the ratio $\lim_{n\to\infty}\mathrm{Prec}[\bm{Z}]_{ij}/\mathrm{Prec}[\bm{Z}]_{i^{\prime}j^{\prime}}$ depends only on $\varphi_{U}$ and is bijective, thus $\varphi_{U}$ is identified. With $\varphi_{U}$ identified, the three cases in Equation 19 form a system of equations solvable for $(\tau_{Z},\varphi_{Z},|\rho|)$. ∎

Based on Theorem 1, if $\rho\varphi_{U}=0$ some parameters are not identifiable by only observing $\bm{Z}$. That is because, when $\rho=0$, $Z$ and $U$ are uncorrelated, hence $Z$ cannot provide any information on $U$, and when $\varphi_{U}=0$, the unmeasured variable is not spatial and the variability in $Z$ cannot be decomposed accordingly (a situation we examine closely in Section 4.5.2). In contrast, when $\rho\varphi_{U}\neq 0$, a number of spatial parameters are identified based solely on the vector of treatments $\bm{Z}$. For intuition about why the unobserved variable’s spatial dependence parameter $\varphi_{U}$ (but not its precision $\tau_{U}$) can be identified by examining the behavior of $\mathrm{Prec}[\bm{Z}]$ away from the tri-diagonal, recall that off of the tri-diagonal, $i$ and $j$ are not neighbors, and that the precision in Equation 19 tells us about the strength of dependence between $Z_{i}$ and $Z_{j}$ given the value of $\bm{Z}$ at all other locations. Our cross-Markov assumption states that, conditionally on $\bm{U}$, the values of $\bm{Z}$ at non-neighboring locations are independent given the values of $\bm{Z}$ at other locations. If $Z_{i}$ and $Z_{j}$ are not independent conditional only on $\bm{Z}$ (and not on $\bm{U}$), this dependence has to arise through paths in Figure 3 that pass through $\bm{U}$. The strength of dependence is a function of $\rho$, the strength of connection between $U$ and $Z$ at a given location (which is independent of the distance between $i$ and $j$), and $\varphi_{U}$ which determines how quickly the dependence of $U_{i},U_{j}$ attenuates with distance $|i-j|$. These dependencies are graphically represented in Figure 3 where the dependence between $Z_{i-1}$ and $Z_{i+1}$, conditional on $\bm{Z}$ at other locations but marginally over $\bm{U}$, arises from the paths $Z_{i-1}-U_{i-1}-U_{i}-U_{i+1}-Z_{i+1}$ and $Z_{i-1}-U_{i-1}-U_{i-2}-\cdots-U_{1}-U_{n}-\cdots U_{i+1}-Z_{i+1}$, and the dependency due to the latter path diminishes as $n$ becomes large. We can identify $\varphi_{U}$ off of the tri-diagonal by examining this attenuation, and $\rho$ is separable from $\tau_{Z}$ only when examining the tri-diagonal as well.

Theorem 2**.**

The parameter $\varphi_{U}$ is identifiable by observing $(\bm{Y},\bm{Z})$ Further, if $\varphi_{U}\neq 0$ then the full parameter $(\beta_{Z},\tau_{Z},\varphi_{Z},\tau_{U},\varphi_{U},\rho,\tau_{\varepsilon})$ is identifiable by observing $(\bm{Y},\bm{Z})$.

Proof.

Note that $\mathrm{E}[\bm{Y}|\bm{Z}]$ is identified irrespective of $\mathrm{Var}[\bm{Y}|\bm{Z}]$, and that by Theorem 1 we can identify whether $\rho\varphi_{U}=0$ by observing $\bm{Z}$.

We start by showing that $\varphi_{U}$ is identifiable by observing $(\bm{Y},\bm{Z})$. Since $\mathrm{Var}[\bm{Y}|\bm{Z}]=\mathrm{\mathbf{G}}_{n}^{-1}+\tau_{\varepsilon}^{-1}\mathrm{\mathbf{I}}_{n}$, by noting the similarity between the expressions for $\mathrm{Var}[\bm{Y}|\bm{Z}]$ and $\mathrm{Prec}[\bm{Z}]$ in (19), we have

[TABLE]

Therefore, $\lim_{n\to\infty}\mathrm{Var}[\bm{Y}|\bm{Z}]_{ij}=0$ for $i\neq j$ if and only if $\varphi_{U}=0$. Since we can identify $\mathrm{Var}[\bm{Y}|\bm{Z}]$, $\varphi_{U}$ is identifiable by $\lim_{n\to\infty}\mathrm{Var}[\bm{Y}|\bm{Z}]_{ij}/\mathrm{Var}[\bm{Y}|\bm{Z}]_{i^{\prime}j^{\prime}}$ for $|i-j|>0$ and $|i^{\prime}-j^{\prime}|=|i-j|+1$.

Next, assume that $\varphi_{U}\neq 0$, and we show that the remaining parameters are identifiable by observing $(\bm{Y},\bm{Z})$. We first note that since $\mathrm{E}[\bm{Y}|\bm{Z}]=\left(\beta_{Z}-\rho\sqrt{\tau_{Z}/\tau_{U}}\right)\bm{Z}$, the combined coefficient of $\bm{Z}$, $\beta_{Z}-\rho\sqrt{\tau_{Z}/\tau_{U}}$ is identifiable. Since we can identify whether $\rho\varphi_{U}\neq 0$ from observing $\bm{Z}$ (Theorem 1), and since we have here that $\varphi_{U}\neq 0$, we can identify whether $\rho=0$, which allows us to consider the cases where $\rho=0$ and $\rho\neq 0$ separately.

If $\rho=0$, the combined coefficient of $\bm{Z}$ is equal to $\beta_{Z}$, and $\beta_{Z}$ is identified (which is expected since for $\rho=0$ there is no confounding by $\bm{U}$). Additionally, we can return to (19) to identify $(\tau_{Z},\varphi_{Z})$ from the first two cases. Finally, with $\varphi_{U}\neq 0$ identified, $\tau_{U}$ can be identified from the off-diagonal elements of $\mathrm{Var}[\bm{Y}|\bm{Z}]$ in Equation 20.

If $\rho\neq 0$, we can identify $(\tau_{Z},\varphi_{Z},\varphi_{U},|\rho|)$ by Theorem 1. Recall that \mathrm{E}[\bm{Y}|\bm{Z}]=\big{(}\beta_{Z}-\mathrm{\mathbf{G}}_{n}^{-1}\mathrm{\mathbf{Q}}_{n}\big{)}\bm{Z} and $-\mathrm{\mathbf{G}}_{n}^{-1}\mathrm{\mathbf{Q}}_{n}=2\rho\sqrt{\frac{\tau_{Z}}{\tau_{U}}}\left(\tau_{U}\mathrm{\mathbf{G}}_{n}^{-1}\right)$, where $\tau_{U}\mathrm{\mathbf{G}}_{n}^{-1}$ does not depend on $\tau_{U}$. Since $\varphi_{U}\neq 0$, $\mathrm{\mathbf{G}}_{n}$ is not a scalar matrix, and the terms $\bm{Z}$ and $\tau_{U}\mathrm{\mathbf{G}}_{n}^{-1}\bm{Z}$ are known and not collinear. Thus we can separately identify their coefficients $\beta_{Z}$ and $\rho\sqrt{\tau_{Z}/\tau_{U}}$. Since $\tau_{U}$ and $\tau_{Z}$ are both positive, and $(\tau_{Z},|\rho|)$ has been previously identified, we can identify $\operatorname*{sign}(\rho)$, and $\tau_{U}$.

In both cases ($\rho=0$ or $\rho\neq 0$), $\tau_{\varepsilon}$ can be identified from the diagonal elements of $\mathrm{Var}[\bm{Y}|\bm{Z}]$.

∎

Theorems 1 and 2 combined establish that when spatial dependencies can be represented using a ring, and the unmeasured variable is spatial ($\varphi_{U}\neq 0$) and is truly a confounder ($\rho\neq 0$), the effect of the exposure $\bm{Z}$, $\beta_{Z}$, is identifiable.

The above results and the details of the proofs indicate that the cross-Markov property is the critical assumption in identifying $\beta_{Z}$. This is not to say that we can remove the more parametric or the constant conditional correlation assumptions, but rather that these assumptions are most likely not the only or weakest ones allowing for identification. In turn, the cross-Markov assumption can be viewed as a relaxation of the usual assumption of no unobserved confounders, expressing that confounding by unobserved spatial variables is local. To see this, the assumption may be rewritten as $Z_{i}\perp\!\!\!\perp\bm{U}_{-i}|(\bm{Z}_{-i},U_{i},\mathrm{\mathbf{X}}_{(-z)})$, expressing that conditional on measured covariates, the value of $U$ at location $i$, and the value of $Z$ everywhere else, $U_{k}$ with $k\neq i$ is not a predictor of $Z_{i}$ and therefore does not confound the relationship between $Z_{i}$ and $Y_{i}$.

Although we have proved these asymptotic results for a ring of locations, the key requirement on the structure of a graph with large connected components is that there are enough pairs of locations at varying lags. This requirement does not seem problematic for, e.g., counties in the United States, as in our study in Section 7.

4.5.2 A non-sufficient set of assumptions: unmeasured non-spatial confounders

Within framework Equation 3, any set of assumptions that suffice for identification of $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ would also allow for identification of the causal parameter $\mu^{\prime}(z)$. In Section 4.5.1 we showed that for adjacency structures described as a ring of growing size, identifiability is achieved when the unmeasured confounder is spatial $(\varphi_{U}\neq 0)$. Here, we establish that identifiability is lost when the unmeasured confounder is not spatial $(\varphi_{U}=0)$ under any adjacency structure.

Assume that $U_{i}$ are independent and identically distributed random variables (exhibiting no spatial structure) with $\mathrm{\mathbf{G}}=\tau_{U}\mathrm{\mathbf{I}}$. Further, assume that $\bm{Z}$ and $\bm{\varepsilon}$ are not spatially structured with precision matrices $\mathrm{\mathbf{H}}=\tau_{Z}\mathrm{\mathbf{I}}$ and $\mathrm{\mathbf{R}}=\tau_{\varepsilon}\mathrm{\mathbf{I}}$. Then, the parameter vector $(\tau_{U},\tau_{Z},\tau_{\varepsilon},\rho)$ in the restricted likelihood is reduced to $(\sigma^{2},\varphi)=(\tau_{U}^{-1}+\tau_{\varepsilon}^{-1},\tau_{Z}(1-\rho^{2}))$, and the parameters in $(\tau_{U},\tau_{\varepsilon})$ and in $(\tau_{Z},\rho)$ are not separately identifiable (the mathematical derivations are included in Appendix E.2). Thus, when there is no spatial structure, $\widebar{\bm{\beta}}=(\mathrm{\mathbf{X}}^{\intercal}\mathrm{\mathbf{X}})^{-1}\mathrm{\mathbf{X}}^{\intercal}(\bm{Y}-\rho\sqrt{\frac{\tau_{Z}}{\tau_{U}}}\bm{Z})$ is not identifiable either.

This result is intuitively obvious. If $\bm{U}$ is not spatially structured, there is no information in the observed data to differentiate outcome model residuals’ variability due to $U$ from that due to $\varepsilon$, and similarly nothing to differentiate intrinsic variability in $Z$ from variability due to $U$. In such case, $\mathrm{E}[\bm{U}|\mathrm{\mathbf{X}}]$ is not identifiable based on observed data indicating that adjustment for $U$ is not possible if the unmeasured confounders do not exhibit spatial structure. This is in line with recent work showing that latent variable approaches cannot be used to acquire identifiability of causal parameters without additional assumptions (D’Amour, 2019; Ogburn et al., 2019).

4.6 Spatially correlated mediating variables

From model Equation 11 and the form of the restricted likelihood in Equation 14, it is evident that information about the elements $(\mathrm{\mathbf{G}},\mathrm{\mathbf{Q}})$ in the bias correction term $-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}(\bm{Z}-\mathrm{\mathbf{X}}_{(-z)}\bm{\gamma})$ is found in the spatial variability of exposure and outcome models’ residuals. However, the estimation procedure cannot differentiate between spatial structure arising from spatial confounders (temporally precedent of $Z$) or variables found on the causal pathway between $Z$ and $Y$ (mediators). If $\bm{Z}$ and $\bm{Y}$ are measured within a small time window, it may be reasonable to assume that there are no spatial covariates mediating the effect of $Z$ on $Y$, and for that reason, our estimates correspond to estimates of $\beta_{z}=\mu^{\prime}(z)$.

On the other hand, in the presence of spatial intermediate variables, estimates using the affine estimator might more closely resemble the direct effect of $Z$ on $Y$, not due to changes to the spatial mediators (Baron and Kenny, 1986). In this setting, the spatial scale restriction provides some protection against adjustment for spatial mediators since variation in spatial scales smaller than that of the exposure is not adjusted for. Therefore, the spatial scale restriction allows for unmeasured spatial confounder bias mitigation while protecting us from adjusting for variables on the causal pathway between $Z$ and $Y$.

5 The affine estimator in non-linear settings: A Bayesian implementation

The REML framework has allowed us to analytically investigate bias and identifiability, by allowing us to integrate out the distribution of the unmeasured confounder from the observed data likelihood. However, such analytical approach is less applicable to non-continuous outcomes and non-linear models. Here, we extend the affine estimator to non-linear models within the Bayesian framework.

5.1 Estimation of causal parameters with non-linear models

When $\mu(z)=\mathrm{E}[Y(z)]$ is the estimand of interest, and the outcome model specifies $\mathrm{E}(Y|Z,\bm{C}^{m},U)$, $\mu(z)$ can be written as $\displaystyle\mathrm{E}_{\bm{C}^{m},U}\left[E(Y|Z=z,\bm{C}^{m},U)\right]$ under assumptions. This is often referred to as the g-formula, or g-computation (Robins, 1986). If the outcome model is linear without exposure-covariate interactions, the regression coefficient for the exposure can be directly interpreted as an estimate of the causal quantity $\mu^{\prime}(z)$. However, the exposure’s coefficient cannot be directly causally interpreted in the case of non-linear models, even when the model is correctly specified. For example, in the context of logistic regression with binary outcomes, the coefficient of the exposure $\beta_{z}$ is not equal to $\mu^{\prime}(z)$, and a linear specification of the exposure-response relationship in the linear predictor of the logistic regression does not imply a linear $\mu(z)$. Therefore, in non-linear outcome models, we need to proceed with care when translating estimated coefficients to estimates of causal quantities, and an integration step (over the distribution of confounders in the target population) needs to be employed in order to acquire estimates of $\mu(z)$ from a non-linear model fit.

Poisson models for count outcomes like the one in our study are an exception. In such models, the parameter $\beta_{z}$ (or $\exp(\beta_{z})$ which is often used in Poisson models) can be interpretable as causal, but for an estimand that is slightly different from $\mu^{\prime}(z)$. Let $P_{i}$ be the population at risk at location $i$, and $Y_{i}(z)$ be the potential outcome at location $i$ if the exposure was set to $z$. Consider the PAERC defined in terms of the standardized outcome rate as $\mathrm{E}[P^{-1}Y(z)]$. Under a structural model similar to Equation 3 for linear $\eta$, assume that [Y(z)|\bm{X}_{(-z)},U]\sim\text{Poisson}\Big{(}P\exp\{\beta_{z}z+\beta_{(-z)}^{\intercal}\bm{X}_{(-z)}+U\}\Big{)}. Then, the PAERC can be written as

[TABLE]

From Equation 21, the coefficient $\beta_{z}$ can be interpreted as the log relative standardized rate for a one-unit exposure change, \log\big{\{}E[P^{-1}Y(z+1)]\ /\ E[P^{-1}Y(z)]\big{\}}. Equivalently, $\beta_{z}$ can be interpreted as the instantaneous effect of the exposure in the relative scale as $\partial\{\log\mathrm{E}[P^{-1}Y(z)]\}/\partial z$. Note that this is substantially different from the standard interpretation of estimated coefficients in Poisson models, and $\beta_{z}$ cannot be used as an estimate of \log\big{\{}E[Y(z+1)/Y(z)]\big{\}}.

5.2 Bayesian implementation of the affine estimator

The presentation above indicates that, for non-linear models where non-collapsibility is an issue and estimated coefficients do not always estimate causal quantities, estimation might require an explicit model for $Y$ conditional on measured covariates and the unmeasured component. This is straightforwardly achieved for the affine estimator within a Bayesian implementation, for which $\bm{U}$ is viewed as a missing variable that is iteratively imputed through a Gibbs sampler. Therefore, placing the affine estimator within the Bayesian paradigm does not require marginalization over $\bm{U}$, which allows for estimation outside the realm of linear regression. This is exploited in Section 6 and Section 7 where we consider a count outcome and a Poisson model with the log link and linear predictor $\eta_{i}=o_{i}+\bm{x}_{i}^{\intercal}\bm{\beta}+u_{i}$, where $o_{i}$ is the offset in the usual sense, and sampling from the conditional posterior distribution of $\bm{\beta}$ can be performed without modifying standard algorithms for Poisson regression.

Apart from its generalizability to non-linear models, the Bayesian approach has a number of additional benefits over the REML approach, including computational gains. The Bayesian implementation is computationally more efficient thanks to the conditional nature of Gibbs sampling which allows us to take advantage of sparsity in $\mathrm{\mathbf{P}}=\mathrm{Prec}[\bm{U},\bm{Z}]$. For example, the log full conditional density of the dependence parameters is (up to an additive constant)

[TABLE]

where $\mathrm{\mathbf{P}}$, $\mathrm{\mathbf{G}}$, $\mathrm{\mathbf{H}}$, and $\mathrm{\mathbf{Q}}$ are sparse matrices depending on $(\tau_{U},\tau_{Z},\varphi_{U},\varphi_{Z},\rho)$, and $p(\tau_{U},\tau_{Z},\varphi_{U},\varphi_{Z},\rho)$ is the prior for these parameters. In contrast, the REML approach requires inverting $\mathrm{\mathbf{G}}$, $\mathrm{\mathbf{R}}$, $(\mathrm{\mathbf{H}}-\mathrm{\mathbf{Q}}^{\intercal}\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}})$, $\mathrm{\mathbf{M}}$, and $\mathrm{\mathbf{C}}^{\intercal}\mathrm{\mathbf{M}}\mathrm{\mathbf{C}}$ in Equation 12 and Equation 14 at each evaluation of the restricted likelihood.

Another advantage of the Bayesian approach is that it is easier to incorporate non-Gaussian exposures by distinguishing between $\bm{Z}$ in Equation 9 and the exposure parameterization in the outcome model. For example, in our analysis of the food access data in Section 7, we replace $Z$ with $\log Z$ in Equation 9 to make the assumption of joint normality more plausible, while using $Z$ in the outcome model to retain the desired interpretation of regression coefficients on the original percentage point scale. Of course, that comes with the caveat that our assumptions are now based on the transformation of the exposure variable. The exposure model and outcome model may be further decoupled by, e.g., assuming joint normality of $\bm{U}$ and a latent variable in a probit model of a binary exposure.

5.3 A regularization prior on the precision matrix of $(\bm{U},\bm{Z})$

The estimation of the joint precision matrix $\mathrm{\mathbf{P}}$ of $(\bm{U},\bm{Z})$ is critical in mitigating bias due to the unobserved spatial confounder $\bm{U}$. However, the present setting is a “low-information” one, as we neither observe $\bm{U}$ directly nor obtain independent replicates. In such settings, the restricted likelihood may have maxima at the boundary of allowed values. For example, Chung et al. (2013) noted that it is not unusual in random effects meta-analysis for the REML estimate of the between-study standard deviation to be zero, and suggested regularizing the REML estimate by multiplying the restricted likelihood by a weakly-informative gamma prior for the between-study variance. Along another thread, Won et al. (2013) considered estimating the covariance matrix in high-dimensional settings where maximum likelihood estimates of such covariance matrices are often ill-conditioned and cannot be inverted accurately. They propose a constrained maximum likelihood approach using the constraint $\kappa(\mathrm{\mathbf{\Sigma}})\leq\kappa_{max}$, where $\kappa(\mathrm{\mathbf{\Sigma}})$ is the condition number (the ratio of the largest to smallest eigenvalue) and $\kappa_{max}$ is pre-specified. They note that this optimization is equivalent to maximizing the likelihood times an exponential prior on $\kappa(\mathrm{\mathbf{\Sigma}})$ left-truncated at 1.

We have observed that this problem manifests in the fully Bayesian implementation as occasional failure of the MCMC sampler to converge. For our purposes, we adopt a truncated exponential prior for $\kappa(\mathrm{\mathbf{P}})$ with rate 1/10 and range $(1,\infty)$ which directly addresses the ill-conditioning problem. With this specification, the difference in log prior density between $\kappa(\mathrm{\mathbf{P}})=1$ and $\kappa(\mathrm{\mathbf{P}})=100$ is $9.9$.

6 Simulation study

6.1 Linear effect

We perform simulations to compare the affine estimator to the non-spatial and spatial random effect estimators under several generative models (GMs). Under all GMs, we consider a single measured covariate, $X$, generated uniformly on $(-1/2,1/2)$. We assume that the mean of $\bm{Z}|\mathrm{\mathbf{X}}_{(-z)}$ is $\bm{X}$. Four GMs reflect $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$ generation according to Equation 9 and Equation 10, with $\bm{U}|\bm{Z},\mathrm{\mathbf{X}}_{(-z)}$ and $\bm{Z}|\bm{U},\mathrm{\mathbf{X}}_{(-z)}$ being one-dimensional CAR models. The within-variable dependence parameters are denoted by $\varphi_{U}$ and $\varphi_{Z}$, and precision parameters by $\tau_{U}=\tau_{Z}=1$. The first model we consider is the unconfounded GM (GM 1), where $\bm{U}$ and $\bm{Z}$ are independent of each other ($\rho=0$), but still spatially structured with $\varphi_{U}=0.5$ and $\varphi_{Z}=0.2$. The unmeasured variable $\bm{U}$ is still predictive of the outcome, hence inducing spatial correlation in the observed outcomes. For the remaining three CAR models, we specify cross-variable dependence ($\rho=0.3$), and vary the within-variable dependence parameters $(\varphi_{U},\varphi_{Z})$ at $(0.5,0.2)$ for GM 2, representing a confounder at a larger spatial scale than the exposure, $(0.2,0.5)$ for GM 3, representing a confounder at a smaller spatial scale than the exposure violating our causal assumptions, and $(0.35,0.35)$ for GM 4, where confounder and exposure vary at the same spatial scale. The fifth and sixth GMs represent situations in which the analysis model is mis-specified. For GM 5, $\bm{U}$ was generated such that its marginal distribution is a one-dimensional CAR model with $\varphi_{U}=0.5$ and $\tau_{U}=1$, and $\bm{Z}|\bm{U},\mathrm{\mathbf{X}}_{(-z)}\sim\mathcal{N}[\bm{U}+\bm{X},\mathrm{\mathbf{I}}]$. Therefore, in this GM, the model Equation 9 is mis-specified in that $\mathrm{\mathbf{G}}$ does not describe the true precision matrix of $\bm{U}|\bm{Z},\mathrm{\mathbf{X}}_{(-z)}$, and the assumption of constant conditional correlation is violated. However, the precision matrix of $\bm{Z}|\bm{U},\mathrm{\mathbf{X}}_{(-z)}$ is still correctly specified, and the cross-Markov property holds. For GM 6, $\tan(U)$ takes the place of $U$ in Equation 9, so that the joint normality assumption on $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$ is violated. In all six GMs, the potential outcomes are generated as $Y_{i}(z)\overset{iid}{\sim}\mathrm{Poisson}\left[\exp\left\{z+X_{i}+U_{i}\right\}\right]$.

Under each GM we generate 500 data sets of size $n=300$ and fit the non-spatial, spatial, and affine estimators. When assumptions on the forms of variances are required, we assume CAR structures, and for the affine estimator we assume that $\mathrm{\mathbf{Q}}$ is of the Markov form Equation 10. Linear predictor models for the exposure and outcome are correctly specified. For the spatial and affine estimators, we evaluate variations with and without the restriction that $\varphi_{Z}<\varphi_{U}$ discussed in Section 4.3.2, with the restricted estimators denoted by (-RS). For regression coefficients we used Gaussian priors with mean 0 and standard deviation 10. For the Spatial-RS, Affine, and Affine-RS estimators we used the regularization prior discussed in Section 5.3. This prior couples the distributions of $\bm{U}$ and $\bm{Z}$ which is not usually a feature in spatial analyses, and so is not used for the Non-spatial and Spatial estimators. Instead, for the Non-spatial and unrestricted Spatial estimators flat priors were used for all variance parameters. In all cases the precision matrix $\mathrm{\mathbf{P}}$ was restricted to be positive definite. Due to the computational cost of computing the condition number of the precision matrix $\mathrm{\mathbf{P}}$ when evaluating the prior for variance parameters, we used an approximation to the condition number acquired by the same model and parameter values on a four-location, one-dimensional ring instead of the 300-location line (the first and last locations are also neighbors). Posterior samples were drawn using 10,000 Gibbs sampler iterations after 1000 burn-in iterations.

Table G.1 displays the simulation results in terms of bias, standard deviation, and root mean squared error (RMSE) of posterior means across data sets, and empirical coverage of 95% equal-tail credible intervals. We first note that there is minimal difference between the Spatial and Spatial-RS estimators, even when the spatial scale assumption is violated in the large-scale exposure scenario. Additionally, the Spatial and Spatial-RS estimators have similar biases to the Non-spatial estimator. For the unconfounded GM 1, all estimators are unbiased, all spatial estimators have approximately correct confidence interval coverage, but both affine estimators have much larger standard errors and therefore RMSE. As expected due to misspecification of the dependence structure, the posterior distributions from the Non-spatial models are too concentrated and therefore the credible intervals are anti-conservative. In all GMs with confounding (2–6), the Affine-RS estimator mitigates bias relative to the Non-spatial estimator, whereas the Spatial and Spatial-RS estimators do not. The unrestricted Affine estimator generally mitigates bias to a lesser extent, especially in the large-scale exposure, same-scales, and non-constant conditional correlation scenarios. When the restricted scale assumption is correct, the Affine-RS estimator has a smaller standard error than the unrestricted Affine estimator. Additionally, the Affine-RS estimator has a smaller RMSE than all other estimators except in the case of a large-scale exposure where its scale restriction is false. In the presence of unobserved confounding, both affine estimators have credible interval coverage rates far superior to the other estimators. Both have approximately nominal coverage rates in the independent, large-scale confounder, and non-normal joint distribution scenarios. Analogous simulation results for maximum a posteriori (restricted likelihood multiplied by a prior) estimation in the Gaussian outcome are available in Appendix G.

Although our model specifies that $\bm{U}$ is independent of $\bm{X}$, we conducted a smaller simulation for a scenario in which the large-scale unmeasured confounder in GM 2 is correlated with the measured confounder, by specifying that $\mathrm{E}[\bm{U}|\bm{X}]=\bm{X}$. We simulated 100 data sets and fit the same models (i.e., without a mean model for $\bm{U}$) as in the previous simulations. The results in estimating $\beta_{z}$ were similar to those under GM 2 (Appendix F) but the estimates of $\beta_{x}$ were biased upward (not shown). This is expected since the part of $U$ that is correlated with $X$ is captured and adjusted for with the inclusion of $X$ in the outcome model, and the affine estimator targets the component of $U$ that is orthogonal to $X$.

6.2 Nonlinear effect

The bias-variance trade-off observed between the spatial and affine estimators in the linear case was also observed for a non-linear effect. We generated 500 data sets of size 300 where $\bm{U},\bm{Z}$ are generated from Equation 9 with $(\tau_{U},\varphi_{U},\tau_{Z},\varphi_{Z},\rho)=(1,0.5,1,0.2,0.3)$, and $Y$ is a Poisson variable with log link and linear predictor in the form of the right-hand side of Equation 15, with $(\beta_{0},\beta_{x})=(0,0)$ and $f(z)=2/(1+e^{-6z})-1$. Therefore, the true effect curve is an anti-symmetric sigmoid with asymptotes $-1$ and $1$. We fit the restricted-scale semiparametric spatial and constrained affine estimators using a penalized cubic spline model with a radial basis and used the same priors as in the linear simulations. Inference was based on 10,000 posterior draws after 5,000 burn-in iterations.

Figure 4 displays a graphical summary of the simulation results. For the most part, both estimators capture the general shape of the mean response curve. However, the spatial estimator is biased toward more extreme estimates as the exposure deviates from 0, and this bias is mitigated by the constrained affine estimator. On the other hand, the constrained affine estimator exhibits substantially greater variability, especially for exposure ranges with limited available data (away from an exposure value of 0). There also appears to be an asymmetry in that when the true log PAERC is negative the bias of the spatial estimator is more pronounced and the constrained affine estimator is more effective at mitigating bias than when the log PAERC is positive. This asymmetry is likely due to the non-linear relationship between the model’s linear predictor and the expected outcome according to the Poisson likelihood.

7 Estimating the county-level effects of poor supermarket availability on CVD mortality

Here, we use the affine estimator in order to estimate the county-level effect of poor supermarket availability on CVD mortality. We consider the affine estimator with and without the spatial scale restriction, and we also consider the spatial random effect estimator with the spatial scale restriction, an extension to the results shown in Section 2. In the outcome model, we include the exposure on the percentage point scale to aid in interpretation of its coefficient. However, we replace $\bm{Z}$ with $\log\bm{Z}$ in the joint model Equation 9 to better satisfy the condition of joint normality. For the Spatial-RS model we restrict $\rho$ to be zero, and for the Spatial-RS and Affine-RS models we apply the constraint $\varphi_{U}\geq\varphi_{Z}$. Note that the prior distribution for the Spatial-RS model differs substantially from that of the unconstrained Spatial model. In all three models we use a similar approximation to the condition number prior of Section 5.3 that was used in the simulation study: rather than computing the condition number on the full joint precision matrix $\mathrm{\mathbf{P}}$, we use the analog of $\mathrm{\mathbf{P}}$ derived from a $4\times 4$ regular grid. Posterior distributions from all models were simulated by retaining 10,000 Gibbs sampler iterations after 1,000 burn-in iterations. The affine-RS model took approximately 1.5 hours to fit on a laptop for a sample size of $n=3,093$.

In Section 7.1 we examine the assumptions underlying the affine estimator, and in Section 7.2 we report summaries of the posterior distribution of the causal effect estimates.

7.1 Examining the plausibility of the assumptions in the context of our study

A number of assumptions, previously presented in Table 1, are necessary to identify the causal effect of interest in the presence of unmeasured spatial confounders.

Temporal ordering and SUTVA are standard assumptions in causal inference and are necessary to define our causal effect. The temporal ordering of the exposure and outcome is immediate satisfied since the exposure data were compiled from 2000 and 2006 data sets, while the mortality outcome data were compiled from 2007 reports. SUTVA is expected to hold, at least approximately, since we can assume that the county-level effect of poor supermarket access on CVD mortality is due to individual-level causal effects, and the home address county listed on death certificates in 2007 corresponds well to the deceased person’s county of residence from 2000 to 2006.

Outcome additivity and the appropriateness of the joint normality and CAR assumptions for $(\bm{U},\bm{Z})|\mathrm{\mathbf{X}}_{(-z)}$ are modeling assumptions that may be at least partially addressed via standard diagnostics. Maps of Pearson residuals based on posterior mean parameters and plots of those residuals versus linear predictors indicated no visually apparent residual spatial correlation or non-linearity in either the log exposure or outcome models in the affine-RS approach. A scatterplot of the joint distribution of the residual log exposure after adjusting for covariates versus mean imputed confounder $\bm{U}$ appeared Gaussian.

The cross-Markov and constant conditional correlation assumptions are assumptions about the relationship between unobserved confounders and the exposure, to which standard diagnostics are not applicable. The plausibility of these assumptions depends heavily on the application and hypothesized confounders. For example, consider an unmeasured variable representing cultural preference toward purchasing prepared food from restaurants versus cooking at home. This variable might act as a confounder in our study, since an increase in such a preference could both depress the demand for and availability of supermarkets, and might drive food choices independent of supermarket availability. The cross-Markov property for this variable (an illustration of which is shown in Figure 5) allows for such cultural preferences to have complex dependence structures across locations. However, grocery store accessibility within a county $i$, $Z_{i}$, is only allowed to depend on such cultural preferences only through its value within the county (conditional on the grocery accessibility in all other locations). This assumption is reasonable for large counties where the food culture in neighboring counties does not directly influence the demand for (and eventual availability of) supermarkets except through its correlation with the food culture within the county itself. Given the cross-Markov property, the constant conditional correlation assumption implies that the strength of the relationship between this aspect of food culture and supermarket availability is constant (conditional on other observed variables in $\mathrm{\mathbf{X}}_{(-z)}$). This can be seen by noting that

[TABLE]

Thus if both the cross-Markov and constant conditional correlation hold, $\sqrt{\tau_{U}/\tau_{Z}}U_{i}$ above behaves like an additive predictor of $Z_{i}$ with regression coefficient $\rho$. A plot of residuals from the regression implied by Equation 23 did not indicate any departures from linearity.

An assumption on the spatial scales of unobserved confounders is also critical for reliable identification of causal effects. We can evaluate this assumption within the model by examining the joint posterior distribution of $(\varphi_{U},\varphi_{Z})$. In our case, the posterior distribution of $\varphi_{U}-\varphi_{Z}$ from the unconstrained Affine model was approximately Gaussian with mean $0.091$, standard deviation $0.024$, and first percentile $0.030$. The posterior from the Affine-RS model was similar, indicating that the assumption is satisfied within the scope of the model.

7.2 Estimating the effect of poor supermarket availability on CVD mortality

Figure 6 displays the posterior distribution of the exponentiated exposure coefficient from the Non-spatial, Spatial, Spatial-RS, Affine, and Affine-RS models. As noted in Section 5.1, this quantity can be interpreted as the relative expected risk of CVD mortality among the population in the 65+ age range due to a one percentage point increase in poor supermarket access in a randomly-chosen county. The agreement of the posterior distribution of $\varphi_{U}-\varphi_{Z}$ in the spatially restricted and unrestricted models implies that the posterior densities from the Spatial and Spatial-RS models closely coincide (both posterior geometric means $0.999$, 95% CIs $(0.988,1.011)$), as do those from the affine and affine-RS models (posterior geometric mean $1.005$, 95% CI $(0.993,1.018)$ and $1.005$ $(0.992,1.018)$, respectively). As we saw in Section 2, the non-spatial model reports a definitive, protective effect of poor supermarket availability (posterior mean $0.968$, 95% CI $(0.962,0.973)$). In contrast, all four spatial models return a smaller and potentially null effect. The Spatial and Spatial-RS models return approximately null effects with posterior probability of a relative risk smaller than 1 equal to 0.54 and 0.53, respectively. In contrast, the Affine and Affine-RS models estimate that poor supermarket access might have a harmful effect on CVD mortality, with posterior probability of a relative risk greater than 1 equal to $0.8$ and $0.79$, respectively.

The change in the point estimate between the spatial random effect models and the affine models is largely attributable to the posterior distribution of $\rho$. This distribution is skewed slightly left, with a posterior mean of $-0.020$ and 95% CI $(-0.040,-0.001)$, indicating confounding by the latent variable $\bm{U}$, even though the posteriors of $\exp(\beta_{Z})$ from the Spatial(-RS) and Affine(-RS) models overlap. Among posterior draws from the Affine-RS model, the correlation between $\rho$ and $\exp(\beta_{Z})$ is $-0.448$, indicating that unobserved spatial confounding is likely to bias the spatial estimator downwards.

We also fit semiparametric versions of each model following the approach in Section 4.4.2. Due to the non-uniform distribution of exposures, we used a truncated cubic basis for the penalized spline. We retained 10,000 MCMC iterations after 5,000 burn-in iterations. The affine-RS model took approximately 2 hours to fit. Results are shown in Figure 7 and are in broad agreement with those from the generalized linear models. The non-spatial model indicates a protective effect of increasing the proportion of households with poor supermarket access on CVD mortality across the observed exposure range. In contrast, the spatial and affine models indicate null and weakly harmful effects, respectively. In all models, the posterior geometric mean of PAERC indicates a protective effect at extremely high levels of poor supermarket access (above 15%), though this is likely an artifact of extrapolation of the spline basis since data are sparse in that range, and the credible bands widen dramatically.

8 Discussion

By positing a joint model for the exposure and unmeasured spatially-correlated confounders, we were able to extend commonly-used spatial data analysis tools to mitigate bias due to such confounders. In contrast to existing approaches that indicate a potentially protective effect of poor supermarket access on cardiovascular disease deaths among the elderly at the county level, the proposed approach leads to estimated effects and inferences that are more plausible and in line with subject-matter knowledge, which indicate that poor supermarket access is likely to be harmful on cardiovascular health. Our study also contributes to the literature on racial and socio-economic disparities. Much of the recent research on food access and health makes reference to food deserts, census tracts with low income and poor food access (United States Department of Agriculture Economic Research Service, 2012). At the neighborhood level, supermarkets are more prevalent in high-income neighborhoods and in neighborhoods with predominantly white residents compared to those whose residents have lower incomes or are black (Morland et al., 2002).

At the same time, there exist limitations to our study that extend outside the scope of unobserved confounders. County-level effects of supermarket access can be extremely relevant for policy making and local planning, but they cannot be directly translated to individual-level effects that may also be of substantial scientific interest. Additionally, although outcomes were measured in the 65+ age range, the exposure and covariates were generally measured across all age ranges. Finally, it is unlikely that effects of poor supermarket access would manifest in differences in mortality in the short-term (e.g., the following year). Instead, effects are more likely to be cumulative over timespans on the order of many years.

From a statistical perspective, our approach to mitigating bias from unobserved spatial confounders is rooted in the causal inference framework and exploits spatial statistics tools that can be used to directly adjust for structured unmeasured confounding. We hope that our work contributes to the growing bridge between spatial data analysis and causal inference. The methodology is intended to be amenable to researchers accustomed to the usual spatial statistics literature, but could potentially be useful in situations calling for mixed models more generally, with appropriate modifications. For that reason, the proposed approach may be widely applicable to scenarios in spatial statistics, time-series analyses, and spatio-temporal settings.

One of the key assumptions in drawing causal conclusions is that of positivity. In the presence of spatial confounders, and for positivity to hold, the spatial scale of the confounder must be larger than that of the treatment. Checking the robustness of estimated effects using the affine estimator with and without the spatial scale restriction can provide intuition on the plausibility of positivity due to unmeasured spatial confounders. This should, however, be employed with care, if spatial mediating variables are expected to exist.

A natural question that arises is whether, and at what occasions, methodology that directly adjusts for unmeasured confounding should be preferred over classical sensitivity analysis. We find that, in settings with structured data, such as spatial and temporal data, unmeasured confounders will also be expected to be structured. In those situations, we find that methodology that directly adjusts for these variables can provide more accurate effect estimates, and strengthen the causal claim of scientific conclusions. An interesting line of future research could extend sensitivity analysis to settings where methods like the one presented here alleviate bias from structured confounders, and sensitivity of results in the presence of unstructured covariates is evaluated.

Based on the structural model Equation 3, we discussed a set of assumptions that allowed for identification of the bias correction term using only the observed data, while providing a simple expression for the expected value of the confounder conditional on the exposure: joint normality of the spatial confounder and exposure of interest, the cross-Markov and constant conditional correlation assumptions. There are several ways in which these assumptions could be relaxed. First, joint normality may be relaxed by assuming joint normality of an underlying random effect process, with the realizations of both the exposure and covariates arising from other distributions. For example, a latent probit model could be used to assess the effects of binary exposures. Furthermore, it may be possible to use more flexible, nonparametric alternatives to the multivariate normal form of the random effect structure, such as spatial Dirichlet processes, both in point-referenced (Gelfand et al., 2005) and areal (Kottas et al., 2008) data. However, questions of identifiability in less parametric settings will likely be difficult to answer. The cross-Markov assumption may be relaxed by, for example, treating the joint distribution of the exposure and confounder as a multivariate conditional autoregressive process (Gelfand and Vounatsou, 2003) and expanding the allowable neighbor relations. The constant conditional correlation assumption may be relaxed by allowing the conditional correlation to vary smoothly in space or based on the number of neighboring locations. In any case, formal treatment of general requirements on the (a) spatial dependence structure (such as the ring graph in Section 4.5.1), including (b) the cross-Markov structure specifying the conditional dependence between the exposure and the unmeasured confounder, and (c) the conditional correlation between the unmeasured confounder and exposure, allowing for identification of the causal effect in the presence of unmeasured confounding is an interesting topic of future research. On a more technical note, priors that place a positive probability on $\tau_{U}^{-1}=0$ could be allowed.

Furthermore, even though the structural model in Equation 3 allows for arbitrary interactions among the exposure and measured covariates (allowing, for example, for treatment effect heterogeneity), it assumes that there are no interactions between measured and unmeasured covariates. The extent to which this assumption can be relaxed is an interesting line of future work, especially in the light of recent results in causal inference for unmeasured confounding (D’Amour, 2019; Ogburn et al., 2019). Therefore, an interesting question that arises is: Can we harvest the spatial information of the data to mitigate bias from unmeasured confounders without imposing structural assumptions?

We consider this to be the most pressing topic of future study: what are the general conditions under which causal effects are identifiable in the presence of unmeasured spatial confounding, and to what extent is bias mitigation robust to model misspecification? Within the context of each study, researchers would need to verify whether the set of reasonable (within their context) assumptions suffices for identification of the bias correction term, while we have illustrated situations in which the causal effect is and is not identifiable.

Acknowledgements

We would like to thank Dr. Lucas R.F. Henneman for his contributions to the original version of the manuscript. We would also like to thank Dr. Jim Hodges for his valuable thoughts and input during various stages of this paper. The work presented here was partially supported by grants UL1TR002733 and R01MH118927 from the National Institutes of Health.

Appendices

Appendix A Additional information on the creation of the data set

We compile a county-level data set including supermarket availability, cardiovascular mortality, demographic and behavioral information for 3,093 out of 3,109 counties or county equivalents in the continental United States.

Supermarket availability

Store data were acquired from the United States Department of Agriculture (Ver Ploeg et al., 2009). For each county (or county equivalent), supermarket availability is defined as the proportion of housing units that are more than 1 mile from the nearest supermarket or large grocery store and do not have a car, obtained from the Food Environment Atlas, June 2012 release (United States Department of Agriculture Economic Research Service, 2012). The data are produced from the United States Department of Agriculture and Economic Research Service June 2009 Report to Congress on access to affordable and nutritious food (Ver Ploeg et al., 2009). The data set was compiled from a 2006 directory of supermarkets and large grocery stores within the continental United States and 2000 Census data. A supermarket or large grocery store was defined as stores that had at least 2 million USD in annual sales and contained all the major food departments found in a traditional supermarket.

Cardiovascular mortality data

County-level data on cardiovascular mortality data were obtained from the United States Centers for Disease Control and Prevention (CDC) WONDER query system (Centers for Disease Control and Prevention, 2018). We obtained county-level population and CVD death counts, both in the 65+ age range, during 2007. A CVD death was defined as a death registered with ICD-10 codes I00–I99 (diseases of the circulatory system) as the underlying cause of death. Due to privacy constraints, county death counts below 10 are censored by CDC WONDER. Death counts were internally standardized to the total CVD death rate in the 65+ population among all United States counties (including those outside of the continental US), which was not influenced by censoring.

Covariate information

Zip code level demographic information was acquired from the 2000 Census. The median household value was set to N/A for 448 (out of 40,646) zip codes for which the Census measurement was equal to zero, but the zip codes population size was positive. We identified the county (or county equivalent) to which a zip code belongs, and set county-level demographic measurements equal to the mean of the zip code-level measurements within the county, excluding missing values. Except for the Census information, we also acquired smoking rates from small-area estimates acquired using the CDC Behavioral Risk Factor Surveillance System data (Dwyer-Lindgren et al., 2014). All covariates along with short description and descriptive statistics are shown in LABEL:app_tab:table1.

Sixteen out of 3,109 counties are excluded from the analysis, leaving 3,093. Broomfield County, Colorado, having been created in 2001, is missing from the supermarket availability data. Loving County, Texas (2000 population 67) does not have a median home value recorded in the census data. Covington city, Emporia city, and Lexington city, Virginia are independent cities (not part of the surrounding counties) with populations on the order of 5,000 and do not appear in the census data set. Eleven contiguous counties in southwestern Georgia are missing demographic information in the census data set for an unknown reason.

Links to data sources and sets

Mortality data can be acquired from the CDC WONDER system at https://wonder.cdc.gov by specifying the following query:

Food access data can be aquired from USDA Food Environment Atlas at https://www.ers.usda.gov/data-products/food-environment-atlas/. We used the June 2012 version. 3. 3.

Complete datasets, including remaining county information are available at in the online data supplement.

Appendix B Identifiability of causal estimands

Here, we review causal identifiability of $\mu(z)=E[Y(z)]$, the expected value of the potential outcome for a fixed treatment $z$ over some population, for a binary treatment $z\in\mathcal{Z}=\{0,1\}$. Note that we do not observe $Y(z)$ for everyone and $\mu(z)$ is an expectation including many unobserved quantities.

A causal estimand is referred to as identifiable under a set of assumptions if it can be written as a function of observables. For $\mu(z)$, on set of assumptions is (1) consistency of potential outcomes, (2) positivity and (3) no unmeasured confounding, since

[TABLE]

where the third equation holds because of the no unmeasured confounding assumption, and the fourth equation holds because of the causal consistency assumption. So $\mu(z)$ is written as a function of the observed outcomes among those with $Z=z$, and for that reason it is identifiable.

Appendix C Implications of dependence assumptions on the joint precision matrix

Here we show how the cross-Markov and constant conditional correlation assumptions determine the matrix $\mathrm{\mathbf{Q}}$. We suppress dependence on $\mathrm{\mathbf{X}}_{(-z)}$ for simplicity. Let

[TABLE]

be the joint precision matrix of $(\bm{U},\bm{Z})$. Then, since the precision matrix of coordinates of a Gaussian variable conditional on other coordinates is obtained by dropping the rows and columns corresponding to those other coordinates,

[TABLE]

Where $i$ is may or may not be equal to $j$, $p_{u_{i}z_{j}}$ is the element of $\mathrm{\mathbf{P}}$ in the row corresponding to $u_{i}$ and column corresponding to $z_{j}$, and similar for other coordinates. In particular,

[TABLE]

If $i\neq j$, the cross-Markov assumption $p(Z_{i}|\bm{Z}_{-i},\bm{U})=p(Z_{i}|\bm{Z}_{-i},U_{i})$ implies that $\mathrm{Cor}[U_{i},Z_{j}|\bm{U}_{-i},\bm{Z}_{-j}]=0$, which in turn implies that $p_{u_{i},z_{j}}$, an arbitrary off-diagonal element of $\mathrm{\mathbf{Q}}$, is zero. Additionally, the constant conditional correlation assumption $\mathrm{Cor}(U_{i},Z_{i}|\bm{U}_{-i},\bm{Z}_{-i},\mathrm{\mathbf{X}}_{(-z)})=\rho$ implies that $p_{u_{i}z_{i}}=-\rho\sqrt{p_{u_{i}u_{i}}p_{z_{i}z_{i}}}$.

Appendix D Restricted likelihood

D.1 Formalization

For simplicity, we assume no covariates in the exposure model, but covariates may be included in the outcome model by inclusion in $\mathrm{\mathbf{X}}$. Based on Equation 9 the covariance matrix is

[TABLE]

Based on the properties of the multivariate normal distribution,

[TABLE]

where

[TABLE]

Based on the above, the marginal variance of $\bm{Z}$ is $(\mathrm{\mathbf{H}}-\mathrm{\mathbf{Q}}^{\intercal}\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}})^{-1}$, and the $\mathrm{Var}[\bm{Y}|\bm{Z}]=\mathrm{Var}[\bm{U}|\bm{Z}]+\mathrm{Var}[\bm{\varepsilon}|\bm{Z}]=\mathrm{\mathbf{G}}^{-1}+\mathrm{\mathbf{R}}^{-1}$. Further, $E[\bm{U}|\bm{Z}]=-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}\bm{Z}$ leading to the following outcome model integrating $\bm{U}|\bm{Z}$ out:

[TABLE]

The full data likelihood can be factored as $f(\bm{y},\bm{u},\bm{z}|\bm{\beta})=f(\bm{y}|\bm{u},\bm{z};\bm{\beta})f(\bm{u}|\bm{z})f(\bm{z})$. Using the outcome model in (D.1) and defining $\mathrm{\mathbf{B}}=-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}$, $\mathrm{\mathbf{A}}=\mathrm{Var}[\bm{Z}]=(\mathrm{\mathbf{H}}-\mathrm{\mathbf{Q}}^{\intercal}\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}})^{-1}$, and $\mathrm{\mathbf{V}}=\mathrm{\mathbf{G}}^{-1}+\mathrm{\mathbf{R}}^{-1}$, we have

[TABLE]

leading to the following restricted likelihood conditional on $\bm{Z}$,

[TABLE]

Since $f(\bm{Z})$ does not depend on $\bm{\beta}$, we can write the full restricted likelihood as

[TABLE]

D.2 Conservative bounds on the conditional correlation

In order to ensure positive definiteness of the precision matrix $\displaystyle\mathrm{\mathbf{P}}=\begin{pmatrix}\mathrm{\mathbf{G}}&\mathrm{\mathbf{Q}}\\ \mathrm{\mathbf{Q}}^{\intercal}&\mathrm{\mathbf{H}}\end{pmatrix}$, $\rho$ has to be constrained. Even though no convenient form of such constraint is available, a conservative one is given by

[TABLE]

where $\lambda_{\mathrm{\mathbf{G}},i}$ and $\lambda_{\mathrm{\mathbf{H}},i}$ are the $i$th eigenvalues of $\mathrm{\mathbf{G}}$ and $\mathrm{\mathbf{H}}$, respectively. To establish that, let

[TABLE]

we constraint $\mathrm{\mathbf{S}},\mathrm{\mathbf{T}}$ such that for any vector $\bm{v}\neq\bm{0}$ of length $2n$, $\bm{v}^{\intercal}\mathrm{\mathbf{P}}\bm{v}=\bm{v}^{\intercal}\mathrm{\mathbf{S}}\bm{v}+\bm{v}^{\intercal}\mathrm{\mathbf{T}}\bm{v}>0$. Let $\mathrm{C}=\{\bm{v}:|\bm{v}|=1\}$. It suffices to show that $\min_{\bm{v}\in\mathrm{C}}\bm{v}^{\intercal}\mathrm{\mathbf{S}}\bm{v}>-\min_{\bm{v}\in\mathrm{C}}\bm{v}^{\intercal}\mathrm{\mathbf{T}}\bm{v}$. Note that $\min_{\bm{v}\in\mathrm{C}}\bm{v}^{\intercal}\mathrm{\mathbf{S}}\bm{v}$ is the minimum eigenvalue of $\mathrm{\mathbf{S}}$, and similarly for $\bm{T}$. Also, since $\mathrm{\mathbf{S}}$ is block diagonal, $\min_{i}\{\lambda_{\mathrm{S}_{i}}\}=\min[\min_{i}\{\lambda_{\mathrm{G}_{i}}\},\min_{i}\{\lambda_{\mathrm{H}_{i}}\}]$. The eigenvalues of $\mathrm{\mathbf{T}}$ are the roots of $|\lambda\mathrm{\mathbf{I}}_{2n}-\mathrm{\mathbf{T}}|=|\lambda\mathrm{\mathbf{I}}_{n}|\cdot|\lambda\mathrm{\mathbf{I}}_{n}-\lambda^{-1}\mathrm{\mathbf{Q}}\mathrm{\mathbf{Q}}|=\prod_{i=1}^{n}(\lambda^{2}-\rho^{2}g_{ii}h_{ii})$, i.e., $\lambda=\pm\rho\sqrt{g_{ii}h_{ii}}$. Thus $-\min_{\bm{v}\in\mathrm{C}}\bm{v}^{\intercal}\mathrm{\mathbf{T}}\bm{v}=|\rho|\sqrt{\max_{i}\{g_{ii}h_{ii}\}}$, and

[TABLE]

guarantees positive definite $\mathrm{\mathbf{P}}$.

D.3 Approximate standard errors accounting for correlation between parameter estimates

For the spatial estimator $\widetilde{\bm{\beta}}$, approximate standard errors are often constructed assuming known variance parameters: \displaystyle\widehat{\mathrm{Var}}\big{(}\widetilde{\bm{\beta}}\big{)}\approx\big{(}\mathrm{\mathbf{X}}^{\intercal}\widehat{\mathrm{\mathbf{V}}}^{-1}\mathrm{\mathbf{X}}\big{)}^{-1}. We do not recommend applying this idea directly to $\widebar{\bm{\beta}}$ due to the fact that the estimates of $\rho$ and $\beta_{z}$ are strongly correlated. We account for such correlation with a small modification. If all variance parameters except $\rho$ are treated as known, then

[TABLE]

where $\mathrm{\mathbf{Q}}^{*}$ is diagonal with elements $q_{ii}^{*}=-\sqrt{g_{ii}h_{ii}}$ known and independent of $\rho$. Then, treating $\rho$ exclusively as a coefficient, we can write $\mathrm{\mathbf{D}}=[\mathrm{\mathbf{X}}|-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}^{*}\bm{Z}]$ via concatenation, and obtain an estimated variance \widehat{\mathrm{Var}}\Big{[}\Big{(}\widebar{\bm{\beta}},\widebar{\rho}\Big{)}\Big{]}=\Big{(}\widehat{\mathrm{\mathbf{D}}}^{\intercal}\widehat{\mathrm{\mathbf{V}}}^{-1}\widehat{\mathrm{\mathbf{D}}}\Big{)}^{-1}, from which an estimate of the variance of $\widebar{\bm{\beta}}$ can be acquired. Based on the estimated variance of $\widebar{\bm{\beta}}$, Wald-type confidence intervals can be obtained.

Approximate standard errors for the semi-parametric estimator may be obtained similarly by augmenting $\mathrm{\mathbf{M}}$ with $-\mathrm{\mathbf{G}}^{-1}\mathrm{\mathbf{Q}}^{*}\bm{Z}$ and $\bm{\theta}$ with $\rho$.

Appendix E Identifiability results

E.1 Matrix results supporting identifiability on the ring graph

We call a matrix STDC if it is symmetric, tridiagonal, and the diagonal, subdiagonal, and superdiagonal are all constant vectors, i.e.,

[TABLE]

Let $\mathrm{R}_{n}$ be the ring graph of order $n$, and $\mathrm{\mathbf{A}}_{n}(\varphi)$ be the unscaled CAR precision matrix

[TABLE]

That is, $\mathrm{\mathbf{A}}_{n}(\varphi)$ is the matrix $\mathrm{\mathbf{S}}_{n}(2,-\varphi)$ with the upper-right and lower-left entries modified to be $-\varphi$.

Lemma 1**.**

For an unscaled CAR precision matrix $\mathrm{\mathbf{A}}_{n}(\varphi)$ of a ring graph $\mathrm{R}_{n}$,

[TABLE]

Proof.

We use the Laplace expansion for computing the determinant. The first term of (E.3) comes directly from the first term of the expansion along the first row. The second and final terms (the only other non-zero terms) are equal because the minor matrices are transposes of each other, and their determinants may be computed via Laplace expansion along the first column of whichever matrix has a non-zero entry in the lower-left corner. In this latter expansion, one matrix is STDC and the other is triangular with a constant diagonal. ∎

Lemma 2**.**

For an STDC matrix $\mathrm{\mathbf{S}}_{n}(2,-\varphi)$ with $|\varphi|<1$,

[TABLE]

Proof.

The recurrence relation can be obtained by computing the Laplace expansion along the first row, yielding two non-zero terms. The minor matrix in one non-zero term is $\mathrm{\mathbf{S}}_{n-1}(2,-\varphi)$, and in the other term the Laplace expansion along the first column has one non-zero term, whose minor matrix is $\mathrm{\mathbf{S}}_{n-2}(2,-\varphi)$. Initial conditions for the recurrence relation $r(n)=2r(n-1)-\varphi^{2}r(n-2)$ can be set to $r(0)=1$ and $r(1)=2$ by letting $\mathrm{\mathbf{S}}_{0}(2,-\varphi)$ be empty and $\mathrm{\mathbf{S}}_{1}(2,-\varphi)=2$, and the characteristic roots technique yields the solution $r(n)=\frac{1+\sqrt{1-\varphi^{2}}}{2\sqrt{1-\varphi^{2}}}\left(1+\sqrt{1-\varphi^{2}}\right)^{n}-\frac{1-\sqrt{1-\varphi^{2}}}{2\sqrt{1-\varphi^{2}}}\left(1-\sqrt{1-\varphi^{2}}\right)^{n}$. ∎

Theorem 3**.**

For an unscaled CAR precision matrix $\mathrm{\mathbf{A}}_{n}(\varphi)$ of a ring graph $\mathrm{R}_{n}$,

[TABLE]

Proof.

The result follows immediately from Lemma 1 and Lemma 2. ∎

Lemma 3**.**

Let $\{C_{ij}^{(n)}\}$ be the cofactors of $\mathrm{\mathbf{A}}_{n}(\varphi)$. Then, for $n>3$,

[TABLE]

Note that the second term of $C_{1,j}^{(n)}$ converges linearly with rate $\varphi$ to 0 as $n$ increases.

Proof.

The first cofactor $C_{1,1}^{(n)}$ can be verified by inspection. In computing the other cofactors along the first row, for each corresponding minor matrix, the first column has $-\varphi$ in the first and final positions, and [math] elsewhere. We examine the minors corresponding to those non-zero positions.

In the first position, each minor matrix is an upper-triangular block matrix: the upper-left $(j-1)\times(j-1)$ block is upper-triangular with $-\varphi$ along the diagonal, the lower-left block is $\mathrm{\mathbf{S}}_{n-j}(2,-\varphi)$, and the lower-left block is $\mathrm{\mathbf{0}}$. Thus the first term in the second line of (E.6).

In the second position, each minor matrix is a lower-triangular block matrix: the upper-left block is $\mathrm{\mathbf{S}}_{j-2}(2,-\varphi)$, the lower-right $(n-j+1)\times(n-j+1)$ block is lower-triangular with $-\varphi$ along the diagonal, and the upper-right block is $\mathrm{\mathbf{0}}$. Thus the second term in the second line of (E.6). ∎

When $\mathrm{\mathbf{A}}_{n}(\varphi)$ is a CAR precision matrix, $C_{1,j}^{(n)}$ is proportional to the marginal correlation between locations $1$ and $j$ (due to Cramer’s rule). The first term in the second line of (E.6) can be thought of as representing the correlation due to the “leftward” path from location $j$ to location $1$, and the second term as that due to the “rightward” path. For a fixed $j$, the “rightward” path becomes long as $n$ increases, thus the fact that the second term’s limit is zero corresponds to the marginal correlation due to the long “rightward” path decreasing to 0.

Theorem 4**.**

For $|\varphi|<1$,

[TABLE]

Proof.

The case $i=1$ follows from Theorem 3 and Lemma 3 via Cramer’s rule. Other cases follow by noting that $\left[\mathrm{\mathbf{A}}_{n}(\varphi)^{-1}\right]_{ij}=\mathrm{Cov}[Z_{i},Z_{j}]$ and thus depends only on $|i-j|$. ∎

E.2 Lack of identifiability in the absence of spatial structure

We illustrate the way in which the proposed method fails to adjust for confounding in the absence of spatial structure. In the interest of simplicity, we do not include measured covariates and $\bm{X}=(\bm{1}\ |\ \bm{Z})$. Suppose that $\mathrm{\mathbf{G}}=\tau_{U}\mathrm{\mathbf{I}}$, $\mathrm{\mathbf{H}}=\tau_{Z}\mathrm{\mathbf{I}}$, and $\mathrm{\mathbf{R}}=\tau_{\varepsilon}\mathrm{\mathbf{I}}$ are all scalar matrices. We then have $\mathrm{\mathbf{Q}}=-\rho\sqrt{\tau_{U}\tau_{Z}}$, and in the notation of (14), $\mathrm{\mathbf{V}}=(\tau_{U}^{-1}+\tau_{\varepsilon}^{-1})\mathrm{\mathbf{I}}$, $\mathrm{\mathbf{A}}=\tau_{Z}^{-1}(1-\rho^{2})^{-1}\mathrm{\mathbf{I}}$, and $\mathrm{\mathbf{B}}=\rho\sqrt{\frac{\tau_{Z}}{\tau_{U}}}\mathrm{\mathbf{I}}$. Note that $\mathrm{\mathbf{A}}$ does not depend on $\tau_{U}$. Various simplifications to Equation 14 are then available:

[TABLE]

The log restricted likelihood is then

[TABLE]

where the last equation holds because $\bm{Z}$ is a column of $\mathrm{\mathbf{X}}$, and therefore $\bm{Z}^{\intercal}\left\{\mathrm{\mathbf{I}}-\mathrm{\mathbf{X}}(\mathrm{\mathbf{X}}^{\intercal}\mathrm{\mathbf{X}})^{-1}\mathrm{\mathbf{X}}^{\intercal}\right\}=\bm{0}$. Treating the restricted likelihood as a function of $\rho$, due to the term $+\frac{n}{2}\log(1-\rho^{2})$, the log restricted likelihood approaches $-\infty$ as $\rho\to\pm 1$, and so all maxima on $\rho\in[-1,1]$ are in the interior.

Writing $(\tau_{U}^{-1}+\tau_{\varepsilon}^{-1})=\sigma^{2}$ and $\varphi=\tau_{Z}(1-\rho^{2})$ we then have

[TABLE]

Thus $\tau_{U}$ and $\tau_{\varepsilon}$ are not identifiable, and the global maximum is at

[TABLE]

Similarly, $\tau_{Z}$ and $\rho$ are not identifiable, and the global maximum is at $\varphi=n/\bm{Z}^{\intercal}\bm{Z}$. As a final result, $\widebar{\bm{\beta}}=(\mathrm{\mathbf{X}}^{\intercal}\mathrm{\mathbf{V}}^{-1}\mathrm{\mathbf{X}})^{-1}\mathrm{\mathbf{X}}^{\intercal}\mathrm{\mathbf{V}}^{-1}(\bm{Y}-\mathrm{\mathbf{B}}\bm{Z})=(\mathrm{\mathbf{X}}^{\intercal}\mathrm{\mathbf{X}})^{-1}\mathrm{\mathbf{X}}^{\intercal}(\bm{Y}-\rho\sqrt{\frac{\tau_{Z}}{\tau_{U}}}\bm{Z})$ is undetermined.

Appendix F Simulation results when unobserved confounder is correlated with observed confounder

Appendix G Simulation results maximum a posteriori estimation

Here we present some simulation results for linear models for continuous outcomes. Estimation is based on the REML approach presented in the manuscript. The exact regularization prior is used because it does not add a substantial computational burden beyond that otherwise necessary for the rest of the MAP estimation procedure.