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.
87,708 characters · 15 sections · 68 citation commands
Rate-Optimal Cluster-Randomized Designs for Spatial Interference
Consider a population of $n$ experimental units. Denote by $Y_i(\bm{d})$ the potential outcome of unit $i$ under the counterfactual that the population is assigned treatments according to the vector $\bm{d} = (d_i)_{i=1}^n \in \{0,1\}^n$, where $d_i=1$ ($d_i=0$) implies unit $i$ is assigned to treatment (control). Treatments assigned to alters can influence the ego since $Y_i(\bm{d})$ is a function of $d_j$ for $j\neq i$, what is known as {\em interference}.
An important estimand of practical interest is the {\em global average treatment effect}
where $\bm{1}_n$ ($\bm{0}_n$) is the $n$-dimensional vector of ones (zeros). This compares average outcomes under the counterfactuals that all or no units are treated. Each average can only be directly observed in the data under an extreme design that assigns all units to the same treatment arm, which would necessarily preclude observation of the other counterfactual. Common designs used in the literature, including those studied here, assign different units to different treatment arms, so neither average is directly observed in the data. Nonetheless, we show that asymptotic inference on $\theta_n$ is possible for a class of cluster-randomized designs under spatial interference where the degree of interference diminishes with distance.
Many phenomena diffuse primarily through physical interaction. The government of a large city may wish to compare the effect of two different policing strategies on crime, but more intensive policing in one neighborhood may displace crime to adjacent neighborhoods blattman2021place,verbitsky2012causal. A rideshare company may wish to compare the performance of two different pricing algorithms, but these may induce behavior that generates spatial externalities, such as congestion. Other phenomena exhibiting spatial interference include infectious diseases in animal donnelly2003impact and human miguel2004worms populations, pollution giffin2020generalized, and environmental conservation programs paler2015social.
Much of the existing literature assumes that interference is summarized by a low-dimensional exposure mapping and that units are individually randomized into treatment or control either via Bernoulli or complete randomization aronow_estimating_2017,basse2019randomization,forastiere2020identification,manski2013identification,toulis2013estimation. Jagadeesan et al.\ jagadeesan2020designs and Ugander et al.\ ugander2013graph also utilize exposure mappings but depart from unit-level randomization. They propose new designs that introduce cluster dependence in unit-level assignments in order to improve estimator precision. We build on this literature by (1) studying rate-optimal choices of both cluster-randomized designs and Horvitz-Thompson estimators, (2) avoiding exposure mapping restrictions on interference, which can be quite strong eckles2017design, and (3) developing a distributional theory for the estimator and a variance estimator.
Regarding (2), most exposure mappings used in the literature imply that only units within a small, known distance from the ego can interfere with the ego's outcome. We instead study a weaker restriction on interference similar to leung2021causal, which states that the degree of interference decreases with distance but does not necessarily zero out at any given distance. This is analogous to the widespread use of mixing-type conditions in the time series and spatial literature instead of $m$-dependence because the latter rules out interesting forms of autocorrelation, including models as basic as AR$(1)$.
Regarding (1), we study cluster-randomized designs in which units are partitioned into spatial clusters, clusters are independently randomized into treatment and control, and the assignment of a unit is dictated by the assignment of its cluster. By introducing correlation in assignments, such designs can avoid overlap problems common under Bernoulli randomization, which improves the rate of convergence. For analytical tractability, we focus on designs in which clusters are equally sized squares, each design distinguished by the number of such squares. We pair each design with a Horvitz-Thompson estimator that compares the average outcomes of units with all or no treated neighbors, where the neighborhood radius is of the same order as the cluster size dictated by the design. See (ref) for a depiction of a hypothetical design and neighborhoods used to construct the estimator.
Our results inform how the analyst should choose the number of clusters (and hence, the cluster size and neighborhood radius of the estimator) to minimize the rate of convergence of the estimator. Notably, existing work on cluster randomization with interference utilizes small clusters (those with asymptotically bounded size). We show that such designs are generally asymptotically biased under the weaker restriction on interference we impose, which motivates the large-cluster designs we study.
Finally, regarding (3), we show that the estimator is asymptotically normal and provide a variance estimator. These results appear to be novel, as no existing central limit theorems seem to apply to our setup in which treatments exhibit cluster dependence, clusters can be large, and units in different clusters are spatially dependent due to interference. As usual, the variance estimator is biased due to heterogeneity in unit-level treatment effects. However, we show that, in a superpopulation setting in which potential outcomes are weakly spatially dependent, the bias is asymptotically negligible.
Based on our theory, we provide practical recommendations for implementing cluster-randomized designs in (ref). Of course, rate-optimality results do not determine the choice of nonasymptotic constants that are often important in practice under smaller sample sizes. Still, they constitute an important first step toward designing practical procedures. Due to the generality of the setting, which imposes quite minimal assumptions on interference, it seems reasonable to first study rate-optimality, as finite-sample optimality appears to require substantial additional structure on the problem. We note that existing results on graph cluster randomization, which require stronger restrictions on interference than this paper, are nonetheless limited to rates, and how “best” to construct clusters in practice has been an open question.
Most of the literature supposes interference is mediated by a network. Studying optimal design in this setting is difficult because network clusters can be highly heterogeneous in topology, and their graph-theoretic properties can closely depend on the generative model of the network leung2021network. We study spatial interference, and to make the optimal design problem analytically tractable, we focus on a class of designs that partitions space into equally sized squares while exploring in simulations the performance of more realistic designs that partition using clustering algorithms. We discuss in (ref) the (pessimistic) prospects of extending our approach to network interference.
There is relatively little work on optimal experimental design under interference. Viviano viviano2020experimental proposes variance-minimizing two-wave experiments under network interference. Baird et al.\ baird2018optimal study the power of randomized saturation designs under partial interference.
A recent literature studies designs for interference that depart from unit-level randomization. A key paper motivating our work is ugander2013graph, who propose graph cluster randomization designs under network interference. Ugander and Yin ugander2020randomized study a new variant of these designs, and harshaw2021design consider related designs for bipartite experiments. These papers assume interference is summarized by exposure mappings, which enables the construction of unbiased estimators and use of designs in which clusters are small. Under our weaker restriction on interference, we show that large clusters are required to reduce bias, which creates a bias-variance trade-off.
Eckles et al.\ eckles2017design show that graph cluster randomization can reduce the bias of common estimators for $\theta_n$ in the absence of correctly specified exposure mappings. Pouget-Abadie et al.\ pouget2018optimizing propose two-stage cluster-randomized designs to minimize bias under a monotonicity restriction on interference. Several papers basse2018model,jagadeesan2020designs,sussman2017elements study linear potential outcome models and propose designs targeting the direct average treatment effect, rather than $\theta_n$. Under a normal-sum model, basse2018model compute the mean-squared error of the difference-in-means estimator, which they use to suggest model-assisted designs.
The aforementioned papers on cluster randomization target global effects such as $\theta_n$ chin2019regression,choi2017estimation. Much of the literature on interference considers fundamentally different estimands defined by exposure mappings. When these mappings are misspecified, the estimands are functions of assignment probabilities, in which case their interpretations can be specific to the experiments run savje2021causal,savje2017average. Hu et al.\ hu2022average (\S 5) views this as “largely unavoidable” in nonparametric settings with interference. Our results show that inference on $\theta_n$, which avoids this issue, is possible under restrictions on interference weaker than those typically used in the literature. Additionally, papers in the literature impose an overlap assumption, which implicitly restricts the estimand leung2021causal. We study cluster-randomized designs that directly satisfy overlap.
There is a large literature on cluster-randomized trials hayes2017cluster,park2020assumption. This literature predominantly studies partial interference, meaning that units are divided into clusters such that those in distinct clusters do not interfere. That is, the clusters themselves impose restrictions on interference. In our setting, clusters are determined by the design and do not restrict interference.
Finally, aronow2020design, pollmann2020causal, and zigler2021bipartite study spatial interference in a different “bipartite” setting in which treatments are assigned to units or locations that are distinct from the units whose outcomes are of interest. This shares some similarities with spatial cluster randomization, where different spatial regions are randomized into treatment, so some of the ideas here may be applicable to optimal design there.
The next section defines the model of spatial interference and the class of designs and estimators studied. In (ref), we derive the estimator's rate of convergence, discuss rate-optimal designs, and provide practical design recommendations. In (ref), we prove that the estimator is asymptotically normal, propose a variance estimator, and characterize its asymptotic properties. We report results from a simulation study in (ref), exploring the use of spectral clustering to implement the designs. Finally, (ref) concludes.
Let $\mathcal{N}_n$ be a set of $n$ units. We study experiments in which units are cluster-randomized into treatment and control, postponing to (ref) the specifics of the design. For each $i\in\mathcal{N}_n$, let $D_i$ be a binary random variable where $D_i=1$ indicates that $i$ is assigned to treatment and $D_i=0$ indicates assignment to control. Let $\bm{D} = (D_i)_{i\in\mathcal{N}_n}$ be the vector of realized treatments and $\bm{d} = (d_i)_{i\in\mathcal{N}_n} \in \{0,1\}^n$ denote a non-random vector of counterfactual treatments. Recall from (ref) that $Y_i(\bm{d})$ is the {\em potential outcome} of unit $i$ under the counterfactual that units are assigned treatments according to $\bm{d}$. Formally, for each $n\in\mathbb{N}$ and $i\in\mathcal{N}_n$, $Y_i(\cdot)$ is a non-random function from $\{0,1\}^n$ to $\mathbb{R}$. We denote $i$'s factual, or observed, outcome by $Y_i = Y_i(\bm{D})$ and maintain the standard assumption that potential outcomes are uniformly bounded.
Thus far, the model allows for unrestricted interference in the sense that $Y_i(\bm{d})$ may vary essentially arbitrarily in any component of $\bm{d}$. In order to obtain a positive result on asymptotic inference, it is necessary to impose restrictions on interference to establish some form of weak dependence across unit outcomes. The existing literature primarily focuses on restrictions captured by {\em $K$-neighborhood exposure mappings}, which imply that $D_j$ can only interfere with $Y_i(\bm{D})$ if the distance between $i,j$ is at most $K$. We will discuss how this assumption is potentially restrictive and establish results under weaker conditions.
We assume each unit is located in $\mathbb{R}^2$. Label each unit by its location, so that $\mathcal{N}_n \subset \mathbb{R}^2$, and equip this space with the sup metric $\rho(i,j) = \max_{t=1,2} \lverti_t-j_t\rvert$ for $i=(i_1,i_2)$, $j=(j_1,j_2)$, and $i,j\in\mathbb{R}^2$. Let $Q(i,r) = \{j\in\mathbb{R}^2\colon \rho(i,j) \leq r\}$, the ball of radius $r$ centered at $i$. Under the sup metric, balls are squares, and the radius is half the side length of the square. Letting $\bm{0}$ denote the origin, we consider a sequence of {\em population regions} $\{Q(\bm{0},R_n)\}_{n\in\mathbb{N}}$ such that
That is, units are located in the square $Q(\bm{0},R_n)$ with growing radius $R_n$. Combined with the next increasing domain assumption, the number of units in the region is $O(n)$, but throughout, we will simply assume the number is exactly $n$.
This allows units to be arbitrarily irregularly spaced, subject to being minimally separated by some distance $\rho_0$, a widely used sampling framework in the spatial literature jenish2009central. In contrast, “infill” asymptotic approaches that do not require minimal separation and instead assume increasingly dense sampling from a fixed region can yield nonstandard limiting behavior lahiri1996inconsistency. For some applications, the spatial distribution of units may exhibit “hotspots” with unusually high densities, perhaps making the infill approach more plausible. Some work adopts a hybrid of the two approaches lahiri2003central,lahiri2006resampling, and it may be possible to extend our results to this framework.
Let $\mathbb{R}_+$ denote the set of non-negative reals and
denote the {\em $K$-neighborhood} of $i$. We study the following model of interference similar to that proposed by leung2021causal.
To interpret this, observe that $\max\big\{\lvertY_i(\bm{d}) - Y_i(\bm{d}')\rvert\colon \bm{d},\bm{d}'\in\{0,1\}^n, d_j=d_j' \,\,\forall j\in\mathcal{N}(i,s)\big\}$ maximizes over pairs of treatment assignment vectors that fix the assignments of units in $i$'s $s$-neighborhood but allow assignments to freely vary outside of this neighborhood. It therefore measures the degree of spatial interference in terms of the maximum change to $i$'s potential outcome caused by manipulating treatments assigned to units $k$ “distant” from $i$ in the sense that $\rho(i,k) > s$. The assumption requires the degree of interference to vanish with the neighborhood radius $s$ so that treatments assigned to more distant alters interfere less with the ego. The rate at which interference vanishes is controlled by $\psi(s)$, which is required to decay at a rate faster than $s^{-2}$.
We next discuss two models of interference satisfying (ref). The first is the standard approach of specifying a $K$-neighborhood exposure mapping. Such a mapping is given by $T_i = T(i,\bm{D},\mathcal{N}_n)$ with the crucial property that its dimension does not depend on $n$, unlike that of $\bm{D}$. The approach is to assume that the low-dimensional $T_i$ summarizes interference by reparameterizing potential outcomes as
That is, once we fix $i$'s exposure mapping $T_i$, its potential outcome is fully determined. No less important, it is also typically assumed that exposure mappings are restricted to a unit's $K$-neighborhood, where $K$ is small, meaning fixed with respect to $n$. Formally, $T(i,\bm{d},\mathcal{N}_n) = T(i,\bm{d}',\mathcal{N}_n)$ for any $\bm{d},\bm{d}'$ such that $d_j=d_j'$ for all $j \in \mathcal{N}(i,K)$, which implies that the treatment assigned to a unit $j$ only interferes with $i$ if $\rho(i,j) \leq K$. In practice, choices with $K=1$ are most common, for example $T_i = (D_i,S_i)$ for $S_i = \bm{1}\{\sum_{j\in\mathcal{N}_n} G_{ij}T_j > 0\}$ or $S_i = \sum_{j\in\mathcal{N}_n} G_{ij}T_j$ where $G_{ij} = \bm{1}\{\rho(i,j)\leq 1\}$. In these examples, $D_i$ captures the direct effect of the treatment, and $S_i$ captures interference from units near $i$.
Crucially, $T_i$ and $K$ must be known to the analyst in this approach, which is often a strong requirement. In contrast, (ref) enables the analyst to impose (ref) {\em while requiring neither to be known}. Indeed, if there exists a $K$-neighborhood exposure mapping satisfying (ref), then (ref) holds with $\psi(s) = c\,\bm{1}\{s \leq K\}$ for some $c$ sufficiently large.
Furthermore, (ref) allows for more complex forms of interference ruled out by (ref) in which interference decays more smoothly with distance, rather than being truncated at some distance $K$. The former is analogous to mixing conditions, which are widespread in the time series and spatial literature, while the latter is analogous to $m$-dependence, which rules out interesting forms of autocorrelation, including models as basic as AR$(1)$.
In the spatial context, our assumption accommodates, for example, the Cliff-Ord autoregressive model cliff1973spatial,cliff1981spatial, which is a workhorse model of spatial autocorrelation used in a variety of fields, including geography getis2008history, ecology valcu2010spatial, and economics anselin2001spatial. A typical formulation of the model is
where we assume $\varepsilon_i$ is uniformly bounded to satisfy (ref). Let $\bm{W}$ be the $n\times n$ spatial weight matrix whose $ij$th entry is $W_{ij}$. These weights typically decay with distance $\rho(i,j)$ in a sense to be made precise below. While this model is highly stylized, the important aspect it captures is autocorrelation through the spatial autoregressive parameter $\lambda$. If this is nonzero, then there is no $K$-neighborhood exposure mapping for which (ref) holds, a point previously noted by eckles2017design in the context of network interference.
To see this, first note that coherency of the model requires nonsingularity of $\bm{I} - \lambda \bm{W}$, where $\bm{I}$ is the $n\times n$ identity matrix. Let $\bm{V}$ be the inverse of this matrix and $V_{ij}$ its entry corresponding to units $(i,j)$. Then the reduced form of the model is
a spatial “moving average” model with spatial weight matrix $\bm{V}$. (See (ref) for some examples of $\bm{V}$.) Noticeably, $Y_i(\bm{D})$ can potentially depend on $D_j$ for any $j\in\mathcal{N}_n$, which is ruled out if one imposes a $K$-neighborhood exposure mapping.
Outcomes satisfying (ref) are near-epoch dependent, a notion of weak spatial dependence, when the weights decay with spatial distance in the following sense:
for some $\gamma>0$ jenish2011spatial. The next result shows that this condition is sufficient for verifying (ref) if $\gamma>2$.
This result shows that, unlike the standard approach of imposing a $K$-neighborhood exposure mapping, (ref) can allow for richer forms of interference in which alters that are arbitrarily distant from the ego can interfere with the ego's response.
Much of the literature on interference considers designs in which units are individually randomized into treatment and control, either via Bernoulli or complete randomization. A common problem faced by such designs is limited overlap, meaning that some realizations of the exposure mapping occur with low probability. For example, suppose that (ref) holds with exposure mapping $T_i = \sum_{j\in\mathcal{N}_n} \bm{1}\{\rho(i,j)\leq K\}D_j$, the number of treated units in $i$'s $K$-neighborhood. Then in a Bernoulli design, for large values of $K$, ${\bf P}(T_i=0)$ is small, tending to zero with $K$ at an exponential rate. This is problematic for a Horvitz-Thompson estimator such as $n^{-1} \sum_{i\in\mathcal{N}_n} (p_i(t)^{-1} \bm{1}\{T_i=t\} - p_i(t')^{-1} \bm{1}\{T_i=t'\})Y_i$ where $p_i(t) = {\bf P}(T_i=t)$ since its variance grows rapidly with $K$ if either $t$ or $t'$ is zero. Ugander et al.\ ugander2013graph propose cluster-randomized designs, which reduce this problem by deliberately introducing dependence in treatment assignments across certain units.
We consider the following class of such designs. We assign units to mutually exclusive clusters by partitioning the population region $Q(\bm{0},R_n)$ into $m_n \leq n$ equally sized squares, assuming for simplicity that $m_n \in \{4^s\colon s\in\mathbb{N}\}$. That is, to obtain increasingly more clusters, we first divide the population region into four squares, then divide each of these squares into four squares, and so on, as in (ref). Label the $m_n$ squares $Q_1, \dots, Q_{m_n}$, and call
cluster $k$. Then the number of units in each cluster is uniformly $O(n/m_n)$ under (ref), and the radius of each cluster is
which we assume is greater than 1. We also assume there are no units on the common boundaries of different squares, so the squares partition $\mathcal{N}_n$.
A {\em cluster-randomized design} first independently assigns each cluster to treatment and control with some probability
that is fixed with respect to $n$. Then within a treated (control) cluster $C_k$, all $i\in C_k$ are assigned $D_i=1$ ($D_i=0$). In order to emphasize that we use this design in later theorems, we state it as a separate assumption.
Note that $m_n$ will be required to diverge with $n$ since a large number of clusters is needed for the estimator to concentrate. If $m_n$ is order $n$, then $r_n=O(1)$, so clusters are asymptotically bounded in size, the usual case studied in the literature, which includes unit-level Bernoulli randomization as a special case. If $m_n$ is of smaller order, then cluster sizes grow with $n$.
To construct the estimator, define the {\em neighborhood exposure indicator}
This is an indicator for whether $i$'s $\kappa_n$-neighborhood is entirely treated ($t=1$) or untreated ($t=0$). Unlike $K$-neighborhood exposure mappings, the radius $\kappa_n$ will be allowed to diverge. Let $p_{ti} = {\bf E}[T_{ti}]$. We study the Horvitz-Thompson estimator
Intuitively, $n^{-1} \sum_{i=1}^n Y_i T_{1i} p_{1i}^{-1}$ estimates $n^{-1} \sum_{i\in\mathcal{N}_n} Y_i(\bm{1}_n)$ using the outcomes of units whose neighbors within radius $\kappa_n$ are all treated. Since the radius depends on $m_n$ through $r_n$, $\hat\theta$ is a function of the number of clusters dictated by the design. (ref) depicts the relationship between the clusters and the $\kappa_n$-neighborhoods that determine exposure.
Since nontrivial designs will include both treated and untreated units, $\hat\theta$ is biased for the global average treatment effect. The choice of design can trade off the size of the bias against that of the variance. In particular, small choices of $m_n$ (few clusters, large radii) induce lower bias and higher variance. In (ref), we discuss nearly rate-optimal choices of $m_n$ for which the bias is asymptotically negligible.
We next derive the rate of convergence of $\hat\theta$ as a function of $n$, $m_n$, and $\psi(\cdot)$, which we use to obtain rate-optimal choices of $m_n$. Recall that designs are parameterized by $m_n$, which determines the number and sizes of clusters, and also that $\hat\theta$ depends on $m_n$ through the neighborhood exposure radius $\kappa_n$, so we will be optimizing over both the design and the radius that determines the estimator.
We first provide asymptotic upper bounds on the bias and variance of $\hat\theta$. For two sequences $\{a_n\}_{n\in\mathbb{N}}$ and $\{b_n\}_{n\in\mathbb{N}}$, we write $a_n \lesssim b_n$ to mean $a_n/b_n = O(1)$ and $a_n \gtrsim b_n$ to mean $b_n/a_n = O(1)$.
Our second result provides asymptotic lower bounds.
The result shows that we can construct potential outcomes satisfying the assumptions of (ref) such that the bias is at least order $\psi(1.5r_n)$ and the variance at least $m_n^{-1}$. As discussed in (ref), existing work on cluster randomization under interference assumes clusters have asymptotically bounded size, which, in our setting, implies $r_n=O(1)$. (ref) implies that the bias of the Horvitz-Thompson estimator can then be bounded away from zero, showing that existing results strongly rely on the exposure mapping assumption to obtain unbiased estimates. In the absence of this assumption, it is necessary to consider designs in which cluster sizes are large to ensure the bias vanishes with $n$.
(ref) implies the mean-squared error of $\hat\theta$ is at most of order $\psi(r_n/2)^2 + m_n^{-1}$, and (ref) provides a similar asymptotic lower bound. Under either bound, the bias increases with $m_n$ while the variance decreases, so there exists a bias-variance trade-off in the choice of design. We next derive rates for $m_n$ that minimize or nearly minimize the upper bound under different assumptions on $\psi(\cdot)$. Based on these results, we make recommendations for practical implementation in the next subsection.
{\em Oracle design.} Suppose $\psi(s)$ is known to decay with $s$ at some rate $\phi(s)$. Then by definition of $r_n$, a rate-optimal design chooses $m_n$ to minimize $\phi( 0.5R_nm_n^{-1/2} )^2 + m_n^{-1}$.
{\em Exposure mappings.} If we assume (ref) holds for some $K$-neighborhood exposure mapping, then $\psi(s) = 0$ for all $s>K$. If $K$ is known, then by choosing $m_n = R_n^2(2K)^{-2}$, we have $\kappa_n = K$ and zero bias. In this case, clusters are asymptotically bounded in size, the estimator converges at rate $n^{-1/2}$, and both the design and estimator qualitatively coincide with those of ugander2013graph.
On the other hand, if $K$ is unknown, then for a nearly rate-optimal design, we can choose $\kappa_n$ to grow at a slow rate so that it eventually exceeds any fixed $K$. This may be achieved by choosing $m_n$ to grow slightly slower than $n$, say $n/\log(n)$. Then for large enough $n$, the bias is zero, and the rate of convergence is $\sqrt{\log(n)/n}$.
{\em Exponential decay.} Common specifications of the spatial weight matrix $\bm{W}$ in the Cliff-Ord model imply that $\psi(s)$ decays exponentially with $s$, for example, the row-normalized matrix
If $\psi(s)$ is known to decay at some exponential rate but the exponent is unknown, then we may choose $m_n = n^{1-\epsilon}$ for any small $\epsilon>0$ for a nearly rate-optimal design, which yields a rate of convergence of $n^{-0.5(1-\epsilon)}$, close to an $n^{-1/2}$-rate. This shows that rates close to $n^{-1/2}$ are attainable in the absence of exposure mapping assumptions, despite targeting the global average treatment effect.
{\em Worst-case decay.} In practice, we may have little prior knowledge about $\psi(s)$. Recall that (ref) requires the rate of decay to be no slower than $s^{-2(1+\epsilon)}$ for $\epsilon>0$. As discussed in (ref), this is the slowest rate for spatial dependence that ensures a finite variance. For this rate, since $R_n$ is order $\sqrt{n}$, the bias is order $(n/m_n)^{-(1+\epsilon)}$. Without knowledge of $\epsilon$, we can settle for a nearly rate-optimal design by setting $\epsilon=0$ and choosing $m_n$ to minimize $(n/m_n)^{-2} + m_n^{-1}$, which yields $m_n = n^{2/3}$ and an $n^{-1/3}$-rate of convergence. Under this design, cluster sizes grow at the rate $r_n^2 = n^{1/3}$.
In the last three designs, the bias is $o(m_n^{-1/2})$, which is of smaller order than the variance. This makes the bias negligible from an asymptotic standpoint, but it would be useful to develop bias-reduction methods. We also reiterate that, while this analysis only provides rates, it is apparently by necessity at this level of generality. A finite-sample optimal design seems to require substantially more knowledge of the functional form of $\psi(\cdot)$.
The designs in the previous section rely on varying degrees of knowledge of $\psi(\cdot)$, the rate at which interference vanishes with distance. In practice, this is likely unknown, so we recommend operating under the worst-case rate of decay discussed in the previous subsection. The {\em default conservative choice} we recommend using in practice is the near-optimal rate described there, namely
To construct the clusters, we recommend partitioning space into $m_n$ clusters using a clustering algorithm, such as spectral clustering. A confidence interval (CI) for $\theta_n$ is given in (ref). In (ref), we explore in simulations the performance of the CI when clusters are constructed according to these recommendations.
Our large-sample theory assumes space is subdivided into evenly sized squares in order to avoid the difficult problem of optimizing over arbitrary shapes. However, since units are typically irregularly distributed in practice, division into equally sized squares may be inefficient, which is why we recommend the use of clustering algorithms. We suggest spectral clustering because it recovers, under weak conditions, low-conductance clusters peng2017partitioning, and low conductance is the key property of clusters utilized in our proofs, as discussed in (ref).
We next state results for asymptotic inference on $\theta_n$. Define $\sigma_n^2 = \text{Var}(\sqrt{m_n}\hat\theta)$.
This is a standard condition and reasonable to impose in light of the lower bound on the variance derived in (ref).
The result centers $\hat\theta$ at its expectation, not the estimand $\theta_n$. However, designs discussed in (ref) result in small bias, meaning $\lvert{\bf E}[\hat\theta] - \theta_n\rvert = o(m_n^{-1/2})$, so we can replace ${\bf E}[\hat\theta]$ with $\theta_n$ on the left-hand side. Also note that the assumption $m_n=o(n)$ implies that cluster sizes grow with $n$. If instead $m_n$ were of order $n$, then $r_n=O(1)$, so by (ref), we would additionally need to assume that there exists a $K$-neighborhood exposure mapping in the sense of (ref) in order to guarantee that the bias vanishes at all. In this case, it is straightforward to establish a normal approximation using existing results.
To our knowledge, there is no off-the-shelf central limit theorem that we can apply to $\hat\theta$. Under (ref), the outcomes appear to be near-epoch dependent on the input process $\{D_i\}_{i\in\mathcal{N}_n}$, but the treatments are cluster-dependent with growing cluster sizes, rather than $\alpha$-mixing, as required by jenish2012spatial. To prove a central limit theorem, they split the average into two parts: its expectation conditional on the dependent input process $\{D_i\}_{i\in\mathcal{N}_n}$, and a remainder that they show is small. Rather than conditioning on all treatments, we find that the following unit-specific conditioning event is more useful for proving our result.
Let $\mathcal{C}_i$ be the cluster containing unit $i$, and $\mathcal{F}_i = \{D_j\colon j \in \mathcal{C}_i \cup \mathcal{N}(i,\kappa_n)\}$. Rewrite the estimator as
We first show that the last term is relatively small, $o_p(m_n^{-1/2})$ to be precise, which means that, on average, $Z_i$ is primarily determined by $\mathcal{F}_i$. The proof of this claim is somewhat complicated jenish2012spatial, but it is similar to the argument showing $[B] \lesssim (nm_n)^{-1/2}$ in (ref). To then establish a central limit theorem for $n^{-1} \sum_{i\in\mathcal{N}_n} {\bf E}[Z_i \mid \mathcal{F}_i]$, we observe that the dependence between “observations” $\{{\bf E}[Z_i \mid \mathcal{F}_i]\}_{i\in\mathcal{N}_n}$ is characterized by the following dependency graph $\bm{A}$, which, roughly speaking, links two units only if they are dependent. Recalling the definition of $\Lambda_i$ from (ref), we connect units $i,j$ in $\bm{A}$ if and only if $j \in \Lambda_i$ (or equivalently $i \in \Lambda_j$). Then $\bm{A}$ is indeed a dependency graph because, under (ref), $j \not\in\Lambda_i$ implies that the treatment assignments that determine $\mathcal{F}_i$ are independent of those that determine $\mathcal{F}_j$. The result follows from a central limit theorem for dependency graphs.
The proof highlights two sources of dependence. The first-order source is the first term on the right-hand side of (ref). Dependence in this average is due to cluster randomization, which induces correlation in the treatments determining $\mathcal{F}_i$ across $i$. The second-order source is the second term on the right-hand side of (ref). Dependence in this average is due to interference, which decays with distance due to (ref). S\"{a}vje savje2021causal derives a similar decomposition in a different context with misspecified exposure mappings. The previous arguments show that the second-order source of dependence is small relative to the first-order source because, with large clusters, dependence induced by cluster randomization dominates dependence induced by interference. This is generally untrue with small clusters.
The proof sketch suggests that, to estimate $\sigma_n^2$, it suffices to account for dependence induced by cluster randomization. Define $A_{ij} = \bm{1}\{j \in \Lambda_i\}$, where $\Lambda_i$ is defined in (ref), and note that $A_{ii}=1$ and $A_{ij}=A_{ji}$. Let $\bar{Z} = n^{-1}\sum_{i\in\mathcal{N}_n} Z_i$, which is equivalent to $\hat\theta$. Our proposed variance estimator is
The bias term $\mathcal{R}_n$ is typically nonzero due to the unit-level heterogeneity. That is, $\lvert{\bf E}[Z_i] - {\bf E}[\bar{Z}]\rvert$ does not approach zero asymptotically, except in the special case of homogeneous treatment effects where $Y_i(\bm{1}_n)-Y_i(\bm{0}_n)$ does not vary across $i$. In the no-interference setting, it is well-known that the variance of the difference-in-means estimator is biased for the same reason and that consistent estimation of the variance is impossible. However, due to the term $m_n/n = o(1)$ in $\mathcal{R}_n$, we will argue that typically $\mathcal{R}_n = o_p(1)$, meaning that $\hat\sigma^2$ is asymptotically exact.
Let us first compare $\mathcal{R}_n$ to its formulation under no interference. In this case, $Y_i(\bm{D}) = Y_i(D_i)$, and we replace $T_{1i}$ with $D_i$ and $T_{0i}$ with $1-D_i$ to estimate the usual average treatment effect. Furthermore, we set $A_{ij}=0$ for all $i\neq j$ because units are independent and set $m_n=n$ since there is no longer a need to cluster units. With these changes, $Z_i = (D_i/p - (1-D_i)/(1-p))Y_i$, and
for $\tau_i = Y_i(1)-Y_i(0)$ and $\bar{\tau} = n^{-1} \sum_{i\in\mathcal{N}_n} \tau_i$. This is the well-known expression for the bias in the absence of interference imbens_causal_2015.
In our setting, we have additional “covariance” terms included in $\mathcal{R}_n$ due to the non-zero off-diagonals of the dependency graph $A_{ij}$. These would be problematic if they were negative and larger in magnitude than the main variance terms since that would make $\hat\sigma^2$ anti-conservative. We show that this occurs with small probability, and in fact, that $\mathcal{R}_n$ is $o_p(1)$. Observe that $m_n/n = o(1)$ and $\hat{\tilde{\sigma}}^2$ has the form of a HAC (heteroskedasticity and autocorrelation consistent) variance estimator andrews1991heteroskedasticity,conley1999gmm. Hence, under conventional regularity conditions, $\hat{\tilde{\sigma}}^2$ is consistent for a variance term $\tilde\sigma^2 \geq 0$, in which case $\mathcal{R}_n$ is non-negative in large samples, and furthermore, $o_p(1)$. To formalize this intuition, we need to specify conditions on the superpopulation from which potential outcomes are drawn. In (ref), we show that, if potential outcomes are $\alpha$-mixing, then $\hat{\tilde{\sigma}}^2$ is asymptotically unbiased for $\tilde{\sigma}^2 = \text{Var}(n^{-1/2} \sum_{i=1}^n {\bf E}[Z_i \mid \{Y_i(\bm{d})\}_{\bm{d}\in\{0,1\}^n}])$, and furthermore, $\text{Var}(\hat{\tilde{\sigma}}^2) = O(n^2/m_n^3)$. Consequently, $\text{Var}(\mathcal{R}_n) = O(m_n^{-1})$ due to the $m_n/n$ term in its expression.
We next present results from a simulation study illustrating the quality of the normal approximation in (ref) and coverage of the CI (ref) when constructing clusters using spectral clustering. To generate spatial locations, let $\{\tilde\rho_i\}_{i\in\mathcal{N}_n}$ be i.i.d.\ draws from $\mathcal{U}([-1,1]^2)$. Unit locations in $\mathbb{R}^2$ are given by $\{\rho_i\}_{i\in\mathcal{N}_n}$ for $\rho_i = R_n\tilde\rho_i$ with $R_n=\sqrt{n}$. We let $\rho(i,j) = \lVert\rho_i - \rho_j\rVert$ where $\lVert\cdot\rVert$ is the Euclidean norm.
We set the number of clusters according to (ref), rounded to the nearest integer, which corresponds to the near-optimal design under the worst-case decay discussed in (ref). To construct clusters, we apply spectral clustering to $\{\rho_i\}_{i\in\mathcal{N}_n}$ with the standard Gaussian affinity matrix whose $ij$th entry is $\text{exp}\{-\rho(i,j)^2\}$. Clusters are randomized into treatment with probability $p=0.5$. (ref) displays the clusters and treatment assignments for a typical simulation draw.
We generate outcomes from three different models. Let $\{\varepsilon_i\}_{i\in\mathcal{N}_n} \stackrel{iid}\sim \mathcal{N}(0,1)$ be drawn independently of the other primitives. The first model is Cliff-Ord:
with $(\alpha,\lambda,\delta,\beta) = (-1,0.8,1,1)$ and spatial weight matrix given by the row-normalized adjacency matrix (ref). As discussed in (ref), this model features exponentially decaying $\psi(s)$, in fact of order $\lambda^s$ leung2021causal.
We construct the second and third models to explore how our methods break down when (ref) is violated or close to violated. For this purpose, we use the “moving average” model (ref) with $(\alpha,\beta) = (-1,1)$ and $V_{ij} = \rho(i,j)^{-\eta}$ for $\eta=4,5$ for the two respective models, so that $\psi(s)$ decays at a polynomial rate. Notably, the choice of $\eta=4$ implies that the rate of decay is slow enough that (ref) can fail to hold. This is because
for some $c>0$ by Lemma A.1(iii) of jenish2009central. The right-hand side does not converge for some $\gamma>2$, as required by (ref). On the other hand, the choice of $\eta=5$ is large enough for (ref) to be satisfied since we now replace the 3 on the right-hand side of the previous display with 4. However, in smaller samples, $\eta=4$ or 5 may not be substantially different, so our methods may still break down from the assumption being “close to” violated.
(ref) displays the results of 5000 simulation draws. Row “Bias$(\hat\theta)$” displays $\lvert{\bf E}[\hat\theta - \theta_n]\rvert$, estimated by taking the average over the draws, while “Var$(\hat\theta)$” is the variance of $\hat\theta$ across the draws. The next rows display the coverage of three different confidence intervals. “Our CI” corresponds to the empirical coverage of (ref). “Naive CI” corresponds to (ref) but replaces $\hat\sigma m_n^{-1/2}$ with the i.i.d.\ standard error, so the extent to which its coverage deviates from 95 percent illustrates the degree of spatial dependence. “Oracle CI” corresponds to (ref) but replaces $\hat\sigma m_n^{-1/2}$ with the “oracle” SE, which is the standard deviation of $\hat\theta$ across the draws. Note that the oracle SE approximates $\sigma_n^2 + \mathcal{R}_n$ because the variance is taken over the randomness of the design as well as of the potential outcomes. Lastly, “SE” displays our standard error $\hat\sigma m_n^{-1/2}$.
There are at most 100 clusters in all designs, and the rate of convergence is quite slow at $n^{-1/3}$ for our choice of $m_n$. Nonetheless, across all designs, the coverage of the oracle CI is close to 95 percent or above, which illustrates the quality of the normal approximation. For the Cliff-Ord model, our CI attains at least 95 percent coverage even for small sample sizes, despite $m_n$ being chosen suboptimally for the worst-case decay. For the moving average model with $\eta=5$, we see some under-coverage in smaller samples due to the larger bias, which is unsurprising from the above discussion, but coverage is close to the nominal level for larger $n$. The results for $\eta=4$, as expected, are worse since it is deliberately constructed to violate our main assumption. Once again, our CI exhibits under-coverage due to the larger bias, but coverage improves and bias decreases as $n$ grows.
This paper studies the design of cluster-randomized experiments targeting the global average treatment effect under spatial interference. Each design is characterized by a parameter $m_n$ that determines the number and sizes of clusters. We propose a Horvitz-Thompson estimator that compares units with different neighborhood exposures to treatment, where the neighborhood radius is of the same order as clusters' sizes given by the design. We asymptotically bound the estimator's bias and variance as a function of $m_n$ and the degree of interference and derive rate-optimal choices of $m_n$. Our lower asymptotic bound shows that designs using small clusters (those with asymptotically bounded sizes) generally result in a non-negligible asymptotic bias. On the other hand, constructing large clusters reduces the total number of clusters, resulting in a bias-variance trade-off that we seek to optimize in terms of rates through the choice of design.
In the worst case where the degree of interference is substantial, the estimator has an $n^{-1/3}$-rate of convergence under a nearly rate-optimal design, whereas in the best case where interference is characterized by a $K$-neighborhood exposure mapping, the rate is $n^{-1/2}$ under a rate-optimal design. We derive the asymptotic distribution of the estimator and provide an estimate of the variance.
Important areas for future research include data-driven choices of $m_n$ and $\kappa_n$ and methods to reduce the bias of the estimator. However, a rigorous theory appears to require more substantive restrictions on interference than what we impose.
Our results focus on the canonical case of spatial data in $\mathbb{R}^2$. We conjecture that they can be extended to $\mathbb{R}^d$ for $d>2$ because our proofs fundamentally rely on the following key property of Euclidean space, which is true for any dimension: it is always possible to construct many clusters with low {\em conductance}, or boundary-to-volume ratio, for example by partitioning space into hypercubes or by spectral clustering leung2021network. This appears in our proofs through the use of Lemma A.1 of jenish2009central, which, together with (ref), is crucial to establish that spatially distant units have small covariance, despite dependence induced by cluster randomization and interference. In this sense, the technical idea behind this paper is to exploit a useful property of Euclidean space -- the existence of many low-conductance clusters -- to show that cluster-randomized designs may be fruitfully applied to the problem of spatial interference.
The story for network interference appears to be different. Existing cluster-randomized designs have theoretical guarantees under exposure mapping assumptions, but it is an open question whether such designs work under weaker restrictions on interference such as (ref). In order to directly apply our idea in the previous paragraph, the network must possess many low-conductance clusters across which we can randomize. Unfortunately, this is a strong requirement in practice because, as discussed in leung2021network, not only do some networks not possess multiple low-conductance clusters, but, of those that do, some apparently possess only a small number of such clusters. Because network “space” differs from Euclidean space in this fundamental aspect, under network interference, clusters can be strongly dependent in the absence of exposure mapping assumptions.
\ifarXiv \foreach \x in {1,...,\the\pdflastximagepages} { \includepdf[pages={\x}]{supplement.pdf} } \fi