Rate-Optimal Cluster-Randomized Designs for Spatial Interference

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 [32]. 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 § 6 the (pessimistic) prospects of extending our approach to network interference.

There is relatively little work on optimal experimental design under interference. Viviano [50] proposes variance-minimizing two-wave experiments under network interference. Baird et al. [5] 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 [46], who propose graph cluster randomization designs under network interference. Ugander and Yin [47] study a new variant of these designs, and [19] 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. [15] 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. [40] propose two-stage cluster-randomized designs to minimize bias under a monotonicity restriction on interference. Several papers [6, 23, 44] study linear potential outcome models and propose designs targeting the direct average treatment effect, rather than $\theta_{n}$. Under a normal-sum model, [6] 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}$ [also see 9, 10]. 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 [42, 43]. Hu et al. [21] (§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 [31]. We study cluster-randomized designs that directly satisfy overlap.

There is a large literature on cluster-randomized trials [e.g. 20, 37]. 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, [4], [39], and [51] 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 § 3, we derive the estimator’s rate of convergence, discuss rate-optimal designs, and provide practical design recommendations. In § 4, 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 § 5, exploring the use of spectral clustering to implement the designs. Finally, § 6 concludes.

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 $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}\lvert i_{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 population regions $\{Q(\bm{0},R_{n})\}_{n\in\mathbb{N}}$ such that

Let $\mathbb{R}_{+}$ denote the set of non-negative reals and

We next discuss two models of interference satisfying Assumption 3. 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

Crucially, $T_{i}$ and $K$ must be known to the analyst in this approach, which is often a strong requirement. In contrast, Assumption 3 enables the analyst to impose (1) while requiring neither to be known. Indeed, if there exists a $K$-neighborhood exposure mapping satisfying (1), then Assumption 3 holds with $\psi(s)=c\,\bm{1}\{s\leq K\}$ for some $c$ sufficiently large.

Furthermore, Assumption 3 allows for more complex forms of interference ruled out by (1) 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 [11, 12], which is a workhorse model of spatial autocorrelation used in a variety of fields, including geography [17], ecology [48], and economics [2]. A typical formulation of the model is

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

Outcomes satisfying (2) are near-epoch dependent, a notion of weak spatial dependence, when the weights decay with spatial distance in the following sense:

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 (1) 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^{\prime})^{-1}\bm{1}\{T_{i}=t^{\prime}\})Y_{i}$ where $p_{i}(t)={\bf P}(T_{i}=t)$ since its variance grows rapidly with $K$ if either $t$ or $t^{\prime}$ is zero. Ugander et al. [46] 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 Figure 2. Label the $m_{n}$ squares $Q_{1},\dots,Q_{m_{n}}$, and call

A cluster-randomized design first independently assigns each cluster to treatment and control with some probability

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 neighborhood exposure indicator

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 § 3, we discuss nearly rate-optimal choices of $m_{n}$ for which the bias is asymptotically negligible.

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.

Theorem 1 implies the mean-squared error of $\hat{\theta}$ is at most of order $\psi(r_{n}/2)^{2}+m_{n}^{-1}$, and Theorem 2 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.

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_{n}m_{n}^{-1/2})^{2}+m_{n}^{-1}$.

Exposure mappings. If we assume (1) 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 [46].

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}$.

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

Worst-case decay. In practice, we may have little prior knowledge about $\psi(s)$. Recall that Assumption 3 requires the rate of decay to be no slower than $s^{-2(1+\epsilon)}$ for $\epsilon>0$. As discussed in Remark 1, 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 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 (11). In § 5, 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 [38], and low conductance is the key property of clusters utilized in our proofs, as discussed in § 6.

To our knowledge, there is no off-the-shelf central limit theorem that we can apply to $\hat{\theta}$. Under Assumption 3, 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 [27]. 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

The proof highlights two sources of dependence. The first-order source is the first term on the right-hand side of (8). 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 (8). Dependence in this average is due to interference, which decays with distance due to Assumption 3. Sävje [42] 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 (4), 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

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

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 [1, 13]. 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 § A, 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.