Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,168 characters · 14 sections · 30 citation commands
A Design-Based Approach to Spatial Correlation
A common concern among empirical researchers is spatial correlation, as policy or event treatments are often spatially correlated. A merger between two gasoline companies, for example, could affect the local gasoline retail market spatially and disproportionately in the neighborhoods close to the rebranded stations houde2012spatial. Closure and demolition of public housing affects the areas closest to the projects more than those farther away aliprantis2015blowing. Along with spatially correlated treatments, there can also be spillover effects to adjacent entities that diminish with distance. One common point of confusion is whether standard errors should be adjusted for spatial correlation. In this paper, we answer this question by disentangling sampling scheme, assignment design, and model specification under a finite population framework.
Adopting the finite population paradigm when studying spatial data has clear advantages. For one, any cases where spatial effects are of interest, the unit of observations is determined by geography or generalized distance measure, and we often sample all units in the population. For example, typically we can collect information on all counties in the United States. As pointed out by pinkse2007central, “with spatial data, it is common for the sample and the population to be the same (e.g., the set of all firms in a market)" (p. 216). The sampling uncertainty underlying the superpopulation approach is unnatural for thinking about uncertainty in such settings; see also abadie2020sampling.
Second, with a well-defined finite population, we can explicitly introduce different sampling schemes after imposing spatial correlations on the population units. Rather than fixing the sites and treating different realizations of the data at all sites as random draws from a superpopulation, it is more practical to treat the lattice as a finite population and draw different collections of sites from it. The latter is something we can actually do in real life, but the former is merely a thought experiment. Along with sampling, we observe only one potential outcome at each location. Therefore, our framework combines sampling- and design-based uncertainty. We can also differentiate between spatial correlation among all neighbors in the population and correlation within a subset of neighbors observed in the sample, which has a major impact on statistical inference. With that said, because we allow the sampling probability to shrink to zero, our asymptotic theory can also accommodate sampling from superpopulations as a special case.
By means of a set of newly developed limit theorems, we derive the finite population spatial heteroskedasticity and autocorrelation consistent (SHAC) variance-covariance matrix of M-estimators. We show that, in the case where a nontrivial fraction of the population is sampled, the new asymptotic variance matrix is smaller than the superpopulation SHAC variance matrix. This finding generalizes neyman1923application's result on conservativeness of the finite population variance estimation and the extension to regression settings in abadie2020sampling. Based on the alternate finite population variance-covariance matrix, we provide guidelines on whether standard errors should be adjusted for spatial correlation instead of relying on heuristic arguments of spatial correlation of unobserved characteristics.
To summarize, whenever there is spatial assignment, meaning that the “treatment” variables are spatially correlated, or a spillover effect is estimated, even in the absence of spillover, one should adjust standard errors for spatial correlation, for instance, by reporting the SHAC standard errors. However, there are a few exceptions. When we independently sample a negligible portion of the population, the Eicker-Huber-White (EHW) standard errors would suffice irrespective of the existence of spatial correlation. Similarly, if we sample a small fraction of clusters from all the clusters in the population, the cluster-robust standard errors would suffice.
Our paper contributes to three strands of literature. Our first contribution is to the literature on limit theorems for random fields. There is a very general asymptotic theory for spatial processes under either the mixing condition or near-epoch dependence in jenish2009central and jenish2012spatial. Unfortunately, due to the introduction of sampling indicators necessary for our paper, their theory cannot incorporate sampling from a superpopulation. Recently, bradley2017central develop a central limit theorem for nonstationary random fields using strong mixing conditions with certain restrictions. By borrowing the techniques in the aforementioned articles, we derive new laws of large numbers and a central limit theorem for near-epoch dependent (NED) processes. We apply these limit theorems to finite population asymptotic theory but they also accommodate superpopulations by including zero sampling probability in the limit. Meanwhile, we allow for nonstationary processes with unbounded moments on irregularly spaced lattices. We also introduce cluster correlation on top of spatial correlation by explicitly including the sampling indicators and allowing for cluster sampling.
Second, we contribute to the literature on finite population inference. Several articles under the finite population framework study the EHW standard error and the cluster-robust standard error for both linear and nonlinear estimators. See, for instance, abadie2020sampling, abadie2017should, xu2021potential, and xu2021asymptotic. bojinov2021panel extend the finite population framework to panel experiments. They allow for spillovers across time periods but maintian the assumption of no spillover across units. savje2021average, savje2021causal, and leung2022causal consider network/spatial interference in randomized experiments with access to the entire population. Their work focuses on estimating the treatment effect or exposure effect consistently. Our paper is the first to examine the necessity of spatial-correlation robust inference under different scenarios for a general class of estimators. Lastly, we contribute to the literature on spatial econometrics. A summary of recent developments can be found in xu2019theoretical.
The remainder of the paper is organized as follows. We derive the asymptotic distribution for M-estimators under finite populations with spatial correlation in Section 2. In Section 3, the asymptotic theory is extended to functions of M-estimators. Among the leading examples are average partial effect (APE) estimators resulting from nonlinear models. Using simulation studies in Section 4, we compare the small sample performance of the EHW, cluster-robust, and SHAC standard errors. The research is concluded in Section 5. Proofs are collected in the appendix.
Let $D\subseteq \mathbb{R}^d$, $d\geq 1$, be a lattice of (possibly) unevenly placed locations in $\mathbb{R}^d$. Consider a sequence of finite subsets of $D$, $\{D_M\}$, where $M$ indexes the sequence of finite populations. $|D_M|$ diverges to infinity in deriving the asymptotic properties, where $|V|$ denotes the cardinality of a finite subset $V\subseteq D$. We adopt the metric $\nu(i,j)=\max_{1\leq l\leq d}|j_l-i_l|$ in space $\mathbb{R}^d$, where $i_l$ is the $l$-th component of $i$. The distance between any subsets $K,V\subseteq D$ is defined as $\nu(K,V)=\inf\{\nu(i,j): i\in K\mbox{ and } j\in V\}$. For any random vector W, $\left\lVertW\right\rVert_p=(E\left\lVertW\right\rVert^p)^{1/p}$, $p\geq 1$, denotes its $L_p$-norm. Let $\mathcal{F}_{iM}(s)=\sigma(U_{jM}; j\in T_M: \nu(i,j)\leq s)$ be the $\sigma$-field generated by the random vectors $U_{jM}$ located in the $s$-neighborhood of location $i$. Lastly, $C$ denotes a generic positive constant that may be different in different circumstances.
Let $(X, z, Y)=\{(X_{iM}, z_{iM}, Y_{iM}), i\in D_M, M\geq 1\}$ and $U=\{U_{iM}, i\in T_M, M\geq 1\}$ be triangular arrays of random fields defined on a probability space $(\Omega, \mathcal{F}, P)$. For each unit $i$, we observe $(X_{iM}, z_{iM}, Y_{iM})$, where $X_{iM}$ is the vector of assignment variables, $z_{iM}$ is a set of attributes, and $Y_{iM}$ is the realized outcome. The categorization of assignments and attributes is based on the posed empirical question. Typically, the key variables of interest in an empirical study could be viewed as assignment variables, and the remaining covariates as attribute variables. There is no restriction in terms of the nature of the triple above: they can be discrete, continuous, or mixed.
We relax the stable unit treatment value assumption in the standard potential outcome framework by allowing for interference among individuals. There exists a mapping, denoted by the potential outcome function $y_{iM}(\bm{x}_M)$, from a vector of assignment variables of all population units to the potential outcomes, where $\bm{x}_M=\{x_{iM}, i\in D_M, M\geq 1\}$.\footnote{I emphasize $\bm{x}_M$ as the argument of the potential outcome function because its realization, $\bm{X}_M$, is the only stochastic vector in the function.} The potential outcome function, $y_{iM}(\cdot)$, along with the observed attributes $z_{iM}$, are non-stochastic. By contrast, the assignment vector $\bm{X}_M$ is random, where $\bm{X}_M=\{X_{iM}, i\in D_M\}$ is the realization of $\bm{x}_M$. As a result, the realized potential outcome, $Y_{iM}=y_{iM}(\bm{X}_M)$, is random. Alternatively, the finite population setting can be understood as conditioned on the potential outcomes and attributes of the $|D_M|$ units in the population. For the most part, we denote $W_{iM}=(\bm{X}_M, Y_{iM})$ for brevity.
We do not take a stance on whether an underlying model is correctly specified. There currently exist discussions on incorporating interference into causal inference and how to consistently estimate the exposure effect with misspecified exposure mappings; see, for instance, hudgens2008toward and savje2021causal. With misspecified interference, we estimate some approximation of the spillover effects in the sample. manski2013identification discusses the identification of potential outcome distributions with social interactions, which is out of the scope of the current paper. Throughout, we assume that the finite population parameters are identified whether or not any feature of the model is correctly specified.
According to the sampling scheme, the population $M$ can be partitioned into $G_M$ mutually exclusive clusters, $\{D_{gM}: g=1,2,\dots,G_M\}$, based on the primary sampling units. $C_{iM}\in\{1,2,\dots,G_M\}$ denotes the cluster that unit $i$ belongs to. Each cluster size is denoted by $|D_{gM}|$ with $|D_M|=\sum_{g=1}^{G_M}|D_{gM}|$, $g\in\{1,2,\dots,G_M\}$. When the primary sampling units are individual entities, random sampling is included as a special case.
As is the starting point in the superpopulation paradigm, we study solutions to a population minimization problem, where the estimand of interest is a $k\times 1$ vector denoted by $\theta^*_M$:
The interpretation of the finite population parameter is up to the researcher; see for instance, rambachan2020design, for discussion on when the finite population estimand carries a causal interpretation. Notice that the expectation $\mathbb{E}$ in equation ((ref)) is taken over the distribution of $X_M$ since $X_M$ is the source of randomness here. The function $q_{iM}(\cdot, \cdot)$ is the objective function for a single unit. The subscripts of the objective function indicate its dependence on the non-stochastic attribute variables $\{z_{iM},i\in D_M\}$.
Let $R_{iM}$ denote the binary sampling indicator, which is equal to one if unit $i$ is sampled. Hence, the sample size is $|D_N|=\sum^{G_M}_{g=1}\sum_{i\in D_{gM}}R_{iM}=\sum_{i\in D_M}R_{iM}$. Below we will be precise about the nature of the sampling scheme but, unless the sample equals the population, the sample size is random. The spatial M-estimator of $\theta^*_M$ is denoted by $\hat{\theta}_N$, which solves the minimization problem in the sample.
In order to establish desirable asymptotic properties of the spatial M-estimator, we need to impose restrictions on the spatial dependence of the stochastic components in the estimation problem. We adopt the definition of mixing coefficients in bradley2017central, near-epoch dependent (NED) random fields in jenish2012spatial, and $m$-dependent random fields in moricz2008strong.
We make the following assumptions. Detailed regularity assumptions are listed in Appendix A.
Assumption 1 is taken from jenish2012spatial. Consistent with the increasing domain asymptotics, the assumption of the minimum distance ensures the expansion of the sample region. In Assumption 2, we introduce the possibility of cluster correlation via the two-stage sampling scheme. It is helpful to go through various sampling schemes resulted from different values of the sampling probabilities. With $\rho_{cM}=\rho_{uM}=1$, we observe the entire population; $\rho_{cM}=1$ and $\rho_{uM}<1$ means random sampling; $\rho_{cM}<1$ and $\rho_{uM}\leq 1$ implies cluster sampling. In particular, sampling from a superpopulation is nested in our unified theory since the sampling probabilities in both steps are allowed to be zero in the limit. When $\rho_c=0$, only a negligible fraction of clusters are sampled from a population of a large number of clusters; while when $\rho_c=1$ and $\rho_u=0$, a negligible portion of units are randomly drawn from a large population. Regardless, the expected sample size, $|D_M|\rho_{uM}\rho_{cM}$, diverges to infinity.
Assumption 3 imposes boundedness of the cluster sizes whenever the population units are partitioned into clusters because of either cluster assignment or cluster sampling. As a result, the number of clusters in the population diverges to infinity along with the population size. Assumption 4 implies that the sampling process and the assignment process are independent of each other, which rules out sample selection due to assignment status.
Assumption 5 is the key weak dependence assumption in the proof of the asymptotic properties. Assignment variables are allowed to be spatially correlated as long as the spatial correlation dies out along with distance. In addition to the mixing condition, we also require the maximal correlation coefficient to be less than one in the limit. With the extension to NED processes in Assumption 6, we can study a more generalized class of random fields.\footnote{Assumption 6 is more than sufficient because in some cases, we only need $L_1$-NED, which can be implied by $L_2$-NED. For certain functions, the NED coefficients are only required to be of size $-d$.} Here, we impose NED conditions on the objective function, score function, and the Hessian matrix directly so that both continuous and discrete variables are allowed. See, for example, Chapter 4 in gallant1988unified for primitive conditions that ensure preservation of the NED property under transformations. Assumption 7 is required to establish consistency of the SHAC variance estimator. The last part of Assumption 7 is a high-level condition, which requires that the kernel weights $\omega\Big(\frac{\nu(i,j)}{b_M}\Big)$ converge to one sufficiently fast as $M\to\infty$.
We introduce the following notation for the variance-covariance matrix of M-estimators. Define
and
where
and
Theorem (ref) is derived using a set of new limit theorems given in the Appendix. It shows that the M-estimators are asymptotically normal with an alternative finite population SHAC variance-covariance matrix. The result holds for spatial assignments both at the individual level or the cluster level. For the latter assignment design, there is cluster assignment in the first place, where assignments within clusters can be strongly correlated. On top of that, the assignments across clusters are allowed to be spatially correlated as long as the key mixing condition in Assumption (ref) is satisfied. With bounded cluster sizes, units far away from each other fall into different clusters and the spatial correlation eventually can die out. For instance, while studying individual outcomes, policies are imposed at the school district level and the school districts nearby may coordinate in certain ways.
Let us first focus on positive sampling probabilities. Based on the variance-covariance matrix in (3) and (4), we can see that sampling scheme plays a role in cluster correlation but not in spatial correlation. The two terms involving spatial correlation, $\Delta_{spatial,M}(\theta^*_M)$ and $\Delta_{ES,M}$, cancel out whenever the assignment variables are independent. Therefore, it is only necessary to adjust the standard errors for spatial correlation if the assignment variables are spatially correlated either because of the assignment design or the inclusion of the spillover effects. Also notice that, even in the absence of spillover effects in the potential outcome function, we manually introduce spatial correlation among assignment variables when we explicitly specify spillover effects in the model. Here, assignment variables include your own “treatment" and your neighbors' “treatments." Another implication of Remark (ref) is that we can test the necessity of spatial-correlation robust inference since assignment variables are observable to us.
We report the SHAC standard errors in the simulation below as one way to adjust standard errors for spatial correlation, as this is the dominant approach in the literature so far. Nevertheless, there are alternative ways to construct standard errors and confidence intervals robust to spatial correlation; see for instance, muller2022spatiala and muller2022spatialb. Regardless, Remark (ref) goes through in the population variance-covariance matrix.
When we switch to zero sampling probabilities, there are exceptions. In theory, when we observe the entire population or sample a large portion from a finite population, the spatial-correlation robust standard errors should be reported to account for spatial assignments. However, if we sample a minimal amount from the population, either the cluster-robust or the EHW standard errors would suffice under spatial assignments. Intuitively, when we independently sample a small proportion of units from the population, the majority of your neighbors would not be observed. Hence, there would be little difference between the EHW, cluster-robust, and the SHAC standard errors. In another case, when we randomly draw a small fraction of clusters from the population, most of your neighbors would be contained in the clusters drawn. Therefore, the cluster-robust standard errors take the majority of the spatial correlation into account and would deviate little from the SHAC standard errors with a sufficient bandwidth. In addition, when either of the sampling probability is zero, we are essentially sampling from a superpopulation, so the usual EHW standard error or the usual cluster-robust standard error would no longer be conservative.
Define
where
and
In Section (ref), we have summarized when you should and should not adjust the standard errors for spatial correlation. In terms of the estimation of the variance matrix, Theorem (ref) shows that when spatial correlation needs to be accounted for, there is an upward bias of the usual SHAC variance estimator. leung2022causal reaches a similar conclusion for a weighted difference-in-means estimator for network data with interference, assuming the entire population is observed. We additionally introduce sampling probabilities and study both linear and nonlinear estimators.
We do not discuss the superpopulation limit of the extra term, $S_E$, as this requires weak dependence assumptions on the superpopulation, which we have not touched upon yet. However, since the usual SHAC variance estimator is consistent for the superpopulation variance-covariance matrix, the finite population SHAC variance matrix should be smaller than the superpopulation version of the variance matrix, in the matrix sense.
Given that we study M-estimators, the APE estimator from nonlinear models would be of great interest. As a result, we study the asymptotic distribution of a generic function of M-estimators to complete the discussion.
Let $f_{iM}(W_{iM},\theta^*_M)$ be a $q \times 1$ function of $W_{iM}$ and $\theta^*_M$. We wish to estimate $\gamma^*_M=\frac{1}{|D_M|}\sum_{i\in D_M}\mathbb{E}\big[f_{iM}(W_{iM},\theta^*_M)\big]$. Let $\hat{\gamma}_N=\frac{1}{|D_N|}\sum_{i\in D_M}R_{iM}f_{iM}(W_{iM},\hat{\theta}_N)$ be the estimator of $\gamma^*_M$. Denote the finite population variance matrix by
And the usual SHAC variance estimator is denoted by $\hat{V}_{f,SN}$. The detailed definition of each term can be found in Appendix (ref).
We see that the same results for M-estimators carry over to functions of M-estimators. For the APE estimator, we should only report the spatial-correlation robust standard errors if assignments are spatially correlated or spillover effects are estimated. Additionally, the usual SHAC standard errors are supposed to be conservative. However, as a well-known fact, the SHAC standard errors suffer from downward bias when the spatial correlation is high. Hence, the actual finite sample performance of the usual SHAC standard errors remains unclear.
We consider an uneven lattice. The population units lie within a square of dimension $\sqrt{M}\times \sqrt{M}$, where M is the population size. The locations $(s_{1,iM}, s_{2,iM})$ are drawn once and kept fixed across designs. $s_{1,iM}\sim \mathcal{U}(0,\sqrt{M})$, $s_{2,iM}\sim \mathcal{U}(0,\sqrt{M})$, and they are independent of each other. The distance between units $i$ and $j$ is measured by $\nu(i,j)=\max(|s_{1,iM}-s_{1,jM}|,|s_{2,iM}-s_{2,jM}|)$. We have four major designs involving spatial assignments at the individual level, spatial assignments at the cluster level, spatial assignments allowing for spillover effects, and spatial assignments in nonlinear models.
The potential outcome function is given below:
where half of $\beta_{ig}$ are equal to 1 and the other half -1. $a$ is a constant scalar. The group unobserved heterogeneity $c_g$ is independently drawn from the standard normal distribution and remains fixed across replications. Without cluster sampling, we can ignore the $g$ subscript. The individual unobservables are generated from the linear spatial autoregressive model below,
where $u_M$ and $\epsilon_M$ are $M\times 1$ vectors. $\epsilon_M$ are i.i.d. draws from a standard normal distribution and kept fixed. $W_u$ is a contiguity matrix and units $i$ and $j$ are neighbors if $\nu(i,j)\leq \sqrt{2}$. It is row-standardized with the diagonal elements being zero.
There are two sub-designs. In the first design, the assignment variables are i.i.d. draws from a Bernoulli distribution with a probability of 0.5 for each replication. The spatial correlation of the individual unobservables (in the sense of superpopulation) depends on $p_u$, which varies from 0 to 0.9 with an increment of 0.1. In the second sub-design, the assignment variables are binary variables equal to one when the input value $\xi_{iM}$ is greater than or equal to its population average $\sum^M_{i=1}\xi_{iM}/M$, where $\xi_M$ is an $M\times 1$ vector drawn from a multivariate normal distribution with mean zero and a variance-covariance matrix equal to $p_x$ raised to the power of the distance. $p_u$ is fixed at 0.3, while $p_x$ takes value from 0 to 0.9. Therefore, individual assignments can be spatially correlated.
The expected size of each dimension of the lattice is 18, unless otherwise noted, which leads to an expected sample size of 324. The clusters in the sampling scheme are formed by grouping the consecutive three units by order, resulting in an expected number of 108 clusters in the sample. There are five sampling schemes: (\romannumeral 1) we observe the entire population; (\romannumeral 2) we independently sample clusters from all the clusters in the population with a probability of 0.25; (\romannumeral 3) we independently draw units from the entire population with a probability of 0.25; (\romannumeral 4) we independently sample clusters from all the clusters in the population with a probability of 0.01; (\romannumeral 5) we independently draw units from the entire population with a probability of 0.01. The last two sampling schemes mimic cluster sampling and independent sampling from the infinite population, respectively. For the last two sampling schemes, the expected size of each dimension of the lattice is decreased to 12 to reduce the computational burden. The constant $a=2$ for the first three sampling schemes and $a=1$ for the last two.
We first regress $Y_i$ on 1 and $X_i$ and report different standard errors of the slope coefficient estimator. For cluster sampling with probability 0.25, we also report results for fixed effect by demeaning variables within clusters. The standard errors among comparison are the EHW standard errors, the cluster-robust standard errors, and the SHAC standard errors. We use the Parzen kernel to estimate the SHAC standard errors. The SHAC standard errors are highly sensitive to the choice of bandwidth. In the first sub-design, where assignments are independent, we report the SHAC standard errors with bandwidth $d^*\in \{1, 2, 3\}$. In the second sub-design, we report the SHAC standard errors with two bandwidths chosen in the following way. The SHAC standard errors are estimated using bandwidth up to 20 with a distance increment of one. Among the 20 bandwidths, We choose either the one that minimizes the mean square error or the one that minimizes the bias of the SHAC standard error with respect to the Monte Carlo standard deviation of the slope coefficient estimator. The number of iterations is 1,000.
When the assignments are independent with the entire population observed, we see from Figure 1 that all standard errors are larger than the Monte Carlo standard deviation of the slope coefficient estimator and the corresponding coverage rates of the 95% confidence interval are almost all above the benchmark line of 0.95.\footnote{In the top panel of all figures below, “oracle" means the Monte Carlo standard deviation of the coefficient estimator or the APE estimator. While in the bottom panels, “oracle" stands for the benchmark coverage rate of the 95% confidence interval, 0.95.} This is expected because the superpopulation standard errors are supposed to be conservative. Among all the standard errors reported, the EHW standard errors are the closest to the Monte Carlo standard deviation, regardless of the spatial correlation of the unobservables in the superpopulation. As well, the coverage rate of the 95% confidence interval based on the EHW standard errors is closest to the theoretical level. The cluster-robust standard errors are generally too large and the SHAC standard errors increase along with the bandwidth. With independent assignments, the latter two standard errors are unnecessarily conservative.
When we sample clusters from the population, things become slightly different. Although assignments are still independent, cluster sampling introduces cluster correlation. As a result, the cluster robust standard error is the one closest to the Monte Carlo standard deviation and the coverage rate of the corresponding confidence interval hovers around the benchmark line. The EHW standard errors are too small, as shown in Figures 2 and 5. We observe the same pattern for the standard errors of the fixed effect estimators. However, the comparison between the cluster robust standard errors and the SHAC standard errors depends on the sampling probability. When the sampling probability of each cluster is nonnegligible, the SHAC standard errors grow along with the bandwidth and eventually become too conservative. On the other hand, when the sampling probability is as small as 0.01, resembling cluster sampling from an infinite population, the difference between the cluster-robust standard errors and the SHAC standard errors narrows down. This is especially true when the chosen bandwidth contains enough nearby units.
When we independently draw units from the population, the EHW standard errors again turn out to be the appropriate one to report. Similarly, when the sampling probability is 0.25, you see a clear discrepancy among the standard errors. When the sampling probability is 0.01, which resembles the case of independent sampling from an infinite population, the EHW, cluster, and SHAC standard errors are almost identical to each other.
In the second sub-design, spatial assignments are introduced through the spatial correlation parameter, $p_x$. When $p_x=0$, we are back to the case of independent assignments with $p_u=0.3$. As we can see from Figures 7 and 10, the EHW standard errors are the ones with the best performance when we either observe the entire population or independently sample units from the population. With cluster sampling, the EHW standard errors underestimate the standard deviation and the SHAC standard errors are quite similar to the cluster-robust standard errors. Nevertheless, with spatial assignments we see that both the EHW and the cluster-robust standard errors gradually become too small when the spatial correlation among the assignment variables increases unless the sampling probability is as small as 0.01. The same observation holds for both pooled OLS and fixed effect estimators when there is cluster partition.
Though not perfect, we see patterns along the lines of the theoretic prediction in Remark (ref) when the sampling probability is 0.01. With cluster sampling, the gap between the SHAC standard errors and the cluster-robust standard errors is much smaller even with high spatial correlation compared with the case of a larger sampling probability. Similarly, with independent sampling, the difference among the EHW, cluster, and the SHAC standard errors is relatively small compared with Figures 7-10. All standard errors suffer from some downward bias. However, keep in mind that the sample size is only over 100 and we allow the correlation parameter, $p_x$, to be as large as 0.9, which results in the correlation between the assignment variables within a distance of one averaging at 0.85. I expect that the discrepancy among the different standard errors would gradually disappear when the sampling probability becomes even lower, especially when the spatial correlation is not too high. However, this would drastically increase the computation burden by enlarging the size of the lattice to a great extent. Hence, results are not reported here due to computation restrictions.
When the ratio of the sample to the population size is small, the intuition of reporting either the EHW or the cluster-robust standard errors remains the same no matter whether we introduce spatial correlation at the cluster level, spillover effects, or nonlinearity later on. Hence, for the simulation designs below, we omit the results for sampling probabilities of 0.01.
In this design, assignments are imposed at the cluster level, and we introduce spatial correlation across clusters. The potential outcome function and the individual unobservables are the same as in equations ((ref)) and ((ref)) with $a=1$. $p_u$ is also fixed at 0.3. On the contrary, the assignments are generated differently. We construct a contiguity matrix among cluster pairs, where the distance between clusters is measured as the minimum distance of units in the cluster pair. Hence, the cluster contiguity matrix $W_G$ is a $G\times G$ matrix, where $G$ is the number of clusters in the population and clusters $l$ and $m$ are neighbors if $\nu(l,m)\leq 2$. We adopt the first three sampling schemes as in the baseline design.
The assignment variables are generated in the following way:
where $\tilde{X}_G$ and $\xi_G$ are $G\times 1$ vectors, and $\xi_g \overset{i.i.d.}\sim \mathcal{N}(0,1), g=1, \dots,G$. The cluster assignment $X_g=\mathbbm{1}\{\tilde{X}_g>\frac{1}{G}\sum^G_{l=1}\tilde{X}_l\}$, $\forall\ g=1,2,\dots, G$. Units within the same cluster receive the same assignment.
When the spatial correlation is imposed at the cluster assignment variables rather than the individual assignment variables, the comparison of the standard errors in the baseline design carries over. In summary, we should report the SHAC standard errors under spatial assignments. In the absence of spatial assignments, we should report the cluster-robust standard errors because cluster assignments always occur regardless of sampling scheme. However, to account for the cluster correlation in addition to the spatial correlation across clusters, I introduce a slightly different distance measure based on the distance between cluster pairs. Using this cluster distance measure, units within the same cluster receive the same weights when estimating the SHAC standard errors. The SHAC standard errors based on the cluster distance slightly outperform the ones based on the unit distance. This is because the latter typically suffer from downward bias even when the spatial correlation is low or moderate. However, the former can also become too conservative when there is a lot of heterogeneity across units.
In this simulation design, we allow for spillover effects in addition to spatial assignments by setting up a linear-in-means type of model. The expected sample size of each dimension of the lattice is 36. The potential outcome function is given below:
where $\bm{x}_M$ is the $M\times 1$ vector containing assignments for all units in the population. The realized assignments $\bm{X}_M$ follows a multivariate normal distribution with mean zero and a variance-covariance matrix equal to $p_x$ raised to the power of the distance. The contiguity matrix, $W_x$, is constructed in the same way as $W_u$. Elements in $W_x$ equal one if the corresponding distance is less than or equal to 0.5. $\beta_{ig}$ are defined in the same way as in Section 4.1. The unobservables $\epsilon_{ig}$ are independently drawn from a standard normal distribution and kept fixed.
We adopt the first three sampling schemes in the baseline design. Within each sampling scheme, we have three assignment and spillover combinations: (\romannumeral 1) $p_x=0$ and $\gamma=0$; (\romannumeral 2) $p_x=0$ and $\gamma=1$; (\romannumeral 3) $p_x=0.1$ and $\gamma=1$. Each represents the case of no spatial assignments and no spillover effects, no spatial assignments with spillover effects, and spatial assignments combined with spillover effects, respectively. Regardless, we estimate spillover effects by regressing $Y_i$ on 1, $X_i$, and $W_s\bm{X}_s$, where $\bm{X}_s$ is the collection of assignments of all units in the sample and $W_s$ is a $N\times N$ contiguity matrix with a distance cutoff of 0.5.
Table 1 reports the coefficient estimates and the standard errors of the spillover effects coefficient. Even though all standard errors suffer from downward bias in all designs, the SHAC standard errors perform the best among the three classes of standard errors. When we include spillover effects in our models, regardless of the sampling scheme we use or whether there are spillover effects in the potential outcome function or not, we should always report the spatial-correlation robust standard errors.
It is worth mentioning that without observing the entire population, what we are identifying is the spillover effects in the sample because we can only observe our neighbors who happen to be selected into the sample, which is an imperfect measure of the spillover in the population. As a result, the coefficient estimator on the spillover effects is biased in the latter two sampling schemes. With cluster sampling, the spillover effects estimator seems to be less biased given that more neighbors are included in the cluster sample.
The coefficient estimates on the individual assignment variable are also reported in Table 2, along with the standard errors. Depending on the sampling scheme, reporting the EHW or cluster-robust standard errors is sufficient. Here, the individual assignment variables and the spillover effects are independent of each other because of independent assignments. In light of this observation, it appears unnecessary to make inference robust to spatial correlation for both coefficient estimators on two assignment variables when one is assigned independently, the other is spatially correlated, and the two assignment variables are independent of each other. We can easily verify the simulation results using the Frisch-Waugh theorem in linear models.
In Table 3, we change $p_x$ to 0.3 for the case of spatial assignments, and omit spillover effects in the estimation by regressing $Y_i$ on 1 and $X_i$ only. Due to the omitted variable bias, the coefficient on the individual assignment variable is biased upward when the assignments are spatially correlated. Nevertheless, we must make inference robust to spatial correlation because of the assignment design. With independent assignments, however, the individual assignment variable and the constructed spillover effect are independent. Therefore, not only is the slope coefficient estimator unbiased, but the standard errors need not be adjusted for spatial correlation either. In other words, even with the spillover effect omitted from the estimation, correct inference depends solely on the nature of the individual assignment variable.
In the last design, we consider spatial assignments in a nonlinear model.
where $X_M$ is an $M\times 1$ continuous random vector following a multivariate normal distribution with mean zero and a variance-covariance matrix equal to $p_x$ raised to the power of the distance. Other aspects of the population generating process are the same as the baseline design. Without spatial assignments, $p_x=0$ and $p_u$ ranges from 0 to 0.9; with spatial assignments, $p_u$ is fixed at 0.3 while $p_x$ varies from 0 to 0.9.
We run probit regressions of $Y_i$ on 1 and $X_i$ and report standard errors of the APE estimator of $X_i$. We report the results when the entire population is observed, which are identical to the linear case with or without spatial assignments. Simulation results from other sampling schemes are also similar to the linear case and hence are omitted.
Using a design-based approach, we identify the sources of uncertainty underlying spatial data. Whenever there are spatial assignments or when spillover effects are estimated, we must make inference robust to spatial correlation, unless the sampling probability is negligible. Since we allow for cluster partition, our results can be applied to short panel data as well with both spatial and temporal dependence.