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.
142,202 characters · 40 sections · 85 citation commands
Cluster-Randomized Trials with Cross-Cluster Interference
\onehalfspacing
\addcontentsline{toc}{part}{Main Paper}
The literature on cluster-randomized trials (CRTs) predominantly assumes partial interference, which allows for interference within but not across clusters. However, researchers often conduct CRTs in what hayes2017cluster refer to as “arbitrary geographical zones,” where clusters are not given by nature since “the population is widely scattered and not divided up into clearly defined and well separated communities.” A well-known concern in these settings is {\em cross-cluster interference}, where units within a cluster may respond to treatments assigned to units in adjacent clusters. Some references refer to this phenomenon as {\em contamination} hudgens2008toward,staples2015incorporating, while others define contamination to mean that control units procure treatment from treated neighbors, which is suggestive of interference hayes2017cluster. In either case, a simple comparison of treatment and control clusters may be biased.
This paper proposes methods for reducing bias due to cross-cluster interference in the analysis and design of cluster-randomized trials. Instead of partial interference, we consider a spatial interference model proposed by leung2022rate which posits that interference decays with geographic distance. This captures the central concern expressed in the previous examples, namely potential interference between geographically proximate units.
At the analysis stage, a common bias-reduction strategy is the “fried-egg design” hayes2017cluster. This is not an experimental design but rather entails excluding from estimation units in the trial located near cluster boundaries, the “whites” of the “fried eggs” that are the clusters. An open question is whether excluding observations in such a fashion is worth the loss of efficiency. We show that conventional difference-in-means estimators can have large biases due to interference near cluster boundaries. We then propose to improve the efficiency of the fried-egg design by excluding not all units near cluster boundaries but rather only the subset not surrounded by clusters assigned to the same treatment arm. We prove that this can substantially reduce the asymptotic order of the bias relative to difference in means with no increase in the asymptotic order of the variance.
We then turn to optimal design of clusters. Partitioning the region into a small number of large clusters reduces bias since fewer units are located near cluster boundaries where they are most prone to cross-cluster interference. However, this comes at the cost of higher variance because the sample size in a CRT is the number of clusters. Unlike the conventional partial interference framework, our spatial interference model induces a bias-variance trade-off, which enables the study of cluster design for optimizing the trade-off.
We prove that the number of clusters $k$ that balances the asymptotic bias and variance of our estimators depends on a parameter $\gamma$ that measures the speed at which interference decays with distance. This formally characterizes how domain knowledge of interference informs optimal cluster construction. Such knowledge is implicitly used in practice when researchers define fried-egg boundaries or construct clusters with “buffer zones” to minimize interaction between units hayes2017cluster.
In practice, cluster construction is {\em ad hoc}. homan2016effect utilize a “traveling salesman algorithm,” essentially an unsupervised learning method. More commonly, researchers use administrative divisions egger2022general or manually partition the study region moulton2001design. binka1998impact study a malaria trial conducted in Northern Ghana, noting that, “Where possible, small paths or road were used to delineate the clusters. However, in most cases, the cluster boundaries did not correspond to natural barriers.” We contribute to the literature by providing formal justification for constructing clusters using $k$-medoids, a well-known unsupervised learning algorithm.
We study CRTs under a design-based framework, as in hudgens2008toward, imai2009essential, and schochet2022design, among others. Unlike these papers, we do not assume partial interference and therefore require a different formulation of standard estimands.
Our work is most closely related to leung2022rate. He studies designs targeting the global average treatment effect (GATE) in which clusters are squares in $\mathbb{R}^2$ with identical areas. We consider a more general set of estimands and richer designs that may utilize unsupervised learning to construct clusters. Also, to accommodate denser spatial settings, we consider “infill-increasing” asymptotics in which the number of units is of larger asymptotic order than the volume of the study region. Finally, we propose new standard errors that are asymptotically conservative without restrictions on the superpopulation.
Our analysis of $k$-medoids builds on cao2024inference. We extend a key result of theirs to the case where $k$ diverges and characterize other $k$-medoid properties in this regime which may be of independent interest. We generalize their Ahlfours-regularity condition on the metric space to allow for infill-increasing asymptotics, though unlike them, we also require an additional boundary condition.
faridani2023rate substantially generalize the theoretical results in leung2022rate while retaining his focus on the GATE. Their results hold for general spaces defined using topological conditions that differ from Ahlfours regularity. For this reason, their proofs differ substantially from ours. They study designs that essentially cover the study region with non-intersecting balls of approximate radius $g_n$, where the radius is chosen to grow at an optimal rate. Our design specifies an optimal choice for the number of clusters and uses unsupervised learning to construct clusters.
The paper is organized as follows. The next section defines the model and estimands. In (ref), we discuss the disadvantages of existing estimators and propose an alternative. We study the theoretical properties of the estimators in (ref) and propose an optimal design in (ref). We present simulation results in (ref) and an empirical application in (ref) using data from the egger2022general trial. Finally, (ref) concludes and summarizes the practical outputs of our analysis.
We will use the following asymptotic order notation. For two sequences of random variables $\{X_n\}_{n\in\mathbb{N}}$ and $\{Y_n\}_{n\in\mathbb{N}}$, we write $X_n \precsim Y_n$ if $\lvertX_n/Y_n\rvert = O_p(1)$, $X_n \prec Y_n$ if $\lvertX_n/Y_n\rvert = o_p(1)$, $X_n \succsim Y_n$ if $\lvertY_n/X_n\rvert = O_p(1)$, $X_n \succ Y_n$ if $\lvertY_n/X_n\rvert = o_p(1)$, and $X_n \sim Y_n$ if both $X_n \precsim Y_n$ and $X_n \succsim Y_n$. For two sequences of constants, we use the same notation for analogous notions of asymptotic boundedness and domination.
Let $(\mathcal{X},\rho)$ be a metric space where $\mathcal{X}$ is the set of spatial locations and $\rho$ the metric. We observe a set of $n$ units $\mathcal{N}_n \subseteq \mathcal{X}$, so $\rho(i,j)$ is the spatial distance between units $i,j\in\mathcal{N}_n$. Let $D_i$ denote unit $i$'s binary treatment assignment and $\bm{D} = (D_i)_{i\in\mathcal{N}_n}$ the vector of observed assignments, which is the only random quantity in our analysis. Let $\{Y_i(\cdot)\}_{i\in\mathcal{N}_n}$ be a set of functions with domain $\{0,1\}^n$ and range $\mathbb{R}$. For $\bm{d} = (d_i)_{i\in\mathcal{N}_n} \in \{0,1\}^n$, $Y_i(\bm{d})$ denotes the potential outcome of unit $i$ under the counterfactual that all units are assigned treatments according to $\bm{d}$. Unit $i$'s observed outcome is $Y_i = Y_i(\bm{D})$. We maintain the following standard assumption, that potential outcomes are uniformly asymptotically bounded. All asymptotic statements are with respect to a sequence indexed by $n\in\mathbb{N}$ unless otherwise indicated.
The primary metric space of interest is $\mathbb{R}^2$, but our results apply to a more general “Ahlfors-regular” space, augmented with a boundary condition. Define the {\em $r$-neighborhood} of unit $i$
The constant $d$ represents the spatial dimension and $\xi_n$ a measure of density (units per unit volume). The dependence on $\xi_n$ is new to this paper relative to prior work on spatial interference and accommodates applications with denser regions. To understand the assumption, consider the standard {\em increasing-domain} case in which $\xi_n = 1$ for all $n$. Part (a) says that the number of units in an $r$-neighborhood is the same order $r^d$ as the neighborhood's volume. This defines a $(C,d)$-finite Ahlfors-regular space, the same metric space studied by cao2024inference. The boundary condition in part (b) is new relative to their setup, but both (b) and the upper bound in (a) are satisfied if $\mathcal{X}=\mathbb{R}^d$ under the usual increasing-domain assumption that units are minimally separated in space jenish2009central.
(ref) also accommodates the {\em infill-increasing} case in which $\xi_n$ diverges with $n$. Here the number of units in any neighborhood is of larger order than the neighborhood volume by a factor of $\xi_n$, corresponding to a denser region. Because we require $\xi_n \prec n$, this is a hybrid of infill and increasing-domain asymptotics in which the volume of the study region $\mathcal{N}_n$ grows but at a slower rate $n/\xi_n$ than the population size $n$ lahiri2006resampling.
Under (ref)(a), the number of units in any $r$-ball is proportional to $\xi_n r^d$, which allows for some variation in density across the study region. The assumption is violated if two subregions of similar volume exist but the number of units in one is a large multiple of the other. This could be the case in practice if the study region encompasses urban and rural areas. In (ref), we discuss how our proposed methodology can be modified to accommodate large variations in density.
This corresponds to Assumption 3 of leung2022rate. Unlike partial interference, it is formulated independently of clusters, enabling us to develop a theory of optimal cluster construction. To interpret the condition, consider a unit $i$ centered at a “cluster” $\mathcal{N}(i,r)$ of radius $r$. The quantity $\lvertY_i(\bm{d}) - Y_i(\bm{d}')\rvert$ measures interference induced by manipulating the treatment assignments of units outside the cluster, and the assumption requires this to decay like $r^{-\gamma}$ or faster. Hence, the larger the minimum distance $r$ between $i$ and the units with manipulated treatments, the smaller the spillover effect. We thus interpret $1/\gamma$ as (an upper bound on) the {\em degree of interference}. In (ref), we discuss the optimal design of clusters using knowledge of $\gamma$.
(ref) further requires $\gamma$ to decay fast enough relative to the spatial dimension $d$. This coincides with the requirement imposed by leung2022rate for $d=2$. In the increasing-domain case, the intuition is that $r$-neighborhood sizes grow like $r^d$ under (ref), so weak dependence requires interference to decay faster than the rate at which neighborhoods densify with $r$. Central limit theorems for spatial processes impose analogous requirements on mixing coefficients which control the degree of spatial dependence, as $\gamma$ does in our framework jenish2009central.
We interchangeably use $k$ or $k_n$ to denote the number of clusters in the design. For a given population $\mathcal{N}_n$, denote by $\mathcal{C}_n = \{C_j\}_{j=1}^k$ the set of clusters, which is a partition of $\mathcal{N}_n$. We study standard two-stage randomized-saturation designs.
Under this design, $k$ clusters are independently randomized into treatment with probability $q$ with $W_j$ denoting cluster $j$'s treatment assignment. If $W_j=t$, all units in cluster $j$ are randomized into treatment with probability $p_t$. The literature often refers to $p_1,p_0$ as {\em saturation levels}.
Let $\xi_n$ and $d$ be given from (ref). We consider (sequences of) clusters satisfying the following properties.
This requires clusters to be globular in that they contain and are contained by balls with radii of the same order. Under (ref)(a), part (b) implies that the number of units in any cluster is order $n/k_n$, which implies that clusters are comparable in size. We show below that clusters generated by $k$-medoids satisfy these requirements. Note that even under partial interference, restrictions on cluster size heterogeneity are required for inference using conventional clustered standard errors hansen2019asymptotic, although our requirements are stronger. These may be possible to relax, but we leave this to future research.
(ref) allows clusters to have quite heterogeneous shapes, as can be seen in (ref). It depicts $k$-medoid clusters, which satisfy (ref) by (ref) below. However, the assumption is violated if clusters are highly size-imbalanced, such as elongated clusters that narrowly encompass a geographical feature such as a river or road. This may be addressed by subdividing large or abnormally shaped clusters or grouping adjacent small clusters. The assumption can also be violated if clusters are not constructed based on spatial proximity. For instance, if units are connected through an online social network connecting units far apart in space, then clusters based on network connectivity would likely violate the assumption.
(ref) is satisfied if the researcher partitions the space into $k_n$ identically-sized cubes, but such clusters do not adapt to the spatial distribution of units. We therefore suggest using unsupervised learning algorithms. The next result provides theoretical guarantees for the well-known $k$-medoids algorithm stated in (ref).
In practice, $k$-means can also be used since it typically delivers clusters similar to those of $k$-medoids. The globular clusters produced by both algorithms are sometimes viewed unfavorably relative to algorithms such as spectral clustering for certain unsupervised learning tasks. However, for our purposes, globular clusters are preferable because they minimize the number of units near cluster boundaries, which are the primary source of bias due to cross-cluster interference. More broadly, (ref) provides general design principles for constructing clusters, namely to aim for balance and globularity, which may be achieved either using these algorithms or manually.
Because we allow for cross-cluster interference, we need to redefine conventional estimands in a manner free of this source of bias. To this end, define for $t \in \{0,1\}$ the $p_t$-{\em counterfactual design}, which sets $k=1$ and $q=t$ in (ref). This design groups the population into a single cluster and assigns treatment with probability $p_t$. Throughout the paper, let ${\bf E}[\cdot]$ denote the expectation taken with respect to the observed design in (ref) and ${\bf E}^*_{p_t}[\cdot]$ the expectation taken with respect to the $p_t$-counterfactual design.
Let $\bm{D}_{-i}$ denote the assignment vector excluding the $i$th component and $Y_i(d, \bm{D}_{-i})$ denote $i$'s potential outcome under the counterfactual that $i$'s observed assignment is $d \in \{0,1\}$, holding fixed the realized assignments of other units. We study estimands of the form
for $d_1,d_0 \in \{\emptyset,0,1\}$, where we define $Y_i(\emptyset,\bm{D}_{-i}) \equiv Y_i(\bm{D})$. The following special cases are analogous to estimands defined by hudgens2008toward and hayes2017cluster. The {\em direct effect} compares treated and untreated units under the $p_1$-counterfactual design:
This is directly identified if the $p_1$-counterfactual design were implemented in practice. The estimands that follow, however, require multiple clusters assigned to different arms because they compare different counterfactual designs. Under our framework, this is the primary motivation for cluster randomization.
The {\em indirect effect} compares outcomes of untreated units under counterfactual designs with different saturation levels:
Interest often centers on the “pure control” baseline of $p_0=0$, so that ${\bf E}^*_{p_0}[Y_i(0,\bm{D}_{-i})] = Y_i(\bm{0})$. The {\em total effect} is the sum of the direct and indirect effects, equal to $\theta_T^* = \theta^*(1, 0; p_1, p_0)$. Finally, the {\em overall effect} compares outcomes under different saturation levels: $\theta_O^* = \theta^*(\emptyset, \emptyset, p_1, p_0)$. leung2024causal provides conditions under which these have causal interpretations.
Let $c(i)\in\{1,\ldots,k\}$ denote the index of the cluster containing unit $i$. A common strategy for estimating $\theta_I^*$ is to compute the difference in means between control units in clusters assigned saturation level $p_1$ and those in clusters assigned level $p_0$:
This is the sample analog of
which may be quite far from the target $\theta_I^*$. Units in control clusters near cluster boundaries may be spatially proximate to units in treated clusters and therefore at greater risk of contamination. For such units $i$, ${\bf E}[Y_i \mid W_{c(i)}=0]$ may be substantially different from ${\bf E}^*_{p_0}[Y_i(0,\bm{D}_{-i})]$.
The fried-egg design attempts to reduce bias by restricting the comparison in (ref) to the subset of units deemed sufficiently far from cluster boundaries. Unfortunately, this has two problems. First, the resulting estimator is in fact asymptotically biased because boundary units are excluded with probability one, so it only estimates an average effect for the subpopulation of units in cluster interiors. As noted by mccann2018reducing, this differs from the target estimand since units in cluster interiors may be systematically different from those near the boundaries due to spatial heterogeneity. Second, it is inefficient. The purpose of only including units in the interiors of control clusters is that such units are “well surrounded” by control clusters, or equivalently, relatively far from treated clusters. However, boundary units may be well surrounded in the same fashion if all neighboring clusters are assigned to control, so it would be just as useful to include these units.
Fix a neighborhood radius $r_n$ to be defined in (ref) below. Call a unit $i$ {\em well-surrounded} if
that is, if a unit $i$'s $r_n$-neighborhood only intersects clusters assigned to the same treatment arm. (ref) depicts $k$-medoid clusters, marking with an “X” units that are not well surrounded. Our strategy is to only exclude from estimation units that are not well surrounded. If $r_n=0$, then all units are well surrounded, and our estimators reduce to difference in means. Choosing a larger radius $r_n$ is analogous to choosing a larger boundary region for exclusion in a fried-egg design.
Under our assumptions, the number of well-surrounded units is of asymptotic order equal to the population size $n$ because any unit has nontrivial probability of being well-surrounded ((ref)). As a consequence, our strategy solves the fried-egg design's boundary bias problem and typically excludes strictly fewer units.
Let $d_1,d_0$ be given from the estimand $\theta^*$. For any $t \in \{0,1\}$ and $i \in \mathcal{N}_n$, define $T_{ti} = \bm{1}\{D_i=d_t, W_{c(i)}=t\}S_i$ if $d_t \in \{0,1\}$ and $T_{ti} = \bm{1}\{W_{c(i)}=t\}S_i$ if $d_t = \emptyset$. Let $p_{ti} = {\bf E}[T_{ti}]$ be the propensity score, which has a closed-form expression given in (ref) below. We propose the following H\'{a}jek estimator for $\theta^*$:
When $r_n=0$, $S_i=1$ for all $i$, so propensity scores are homogeneous across $i$, and $\hat\theta$ reduces to difference in means. For $r_n>0$, the scores are generally spatially heterogeneous. For instance, boundary units are less likely to be well surrounded than interior units because that would require more clusters to be assigned the same saturation level.
For any cluster $C_j$, let $m_j$ be its centroid from (ref) and $R_j = \max_{i \in C_j} \rho(i,m_j)$, the “radius” of the cluster. We propose setting
In the case where clusters are equally sized squares, this coincides with the radius suggested by leung2022rate.
Lastly, we propose a variance estimator for $\hat\theta$. Let
the set of units $\ell$ for which some cluster intersects the $r_n$-neighborhood of $\ell$ and $i$. These can be thought of as the units “most potentially correlated” with $i$. Define $A_{ij}(1) = \bm{1}\{j \in \Lambda_i\}$, $A_{ij}(2) = \bm{1}\{j \in C_{c(i)}\}$, and $\hat{Z}_i = (T_{1i}(Y_i - \hat\mu_1))/p_{1i} - (T_{0i}(Y_i - \hat\mu_0))/p_{0i}$. The variance estimator is
Notice that $\hat\sigma^2(2)$ is the conventional cluster-robust variance estimator baird2018optimal, which only accounts for within-cluster dependence. The estimator $\hat\sigma^2(1)$ is analogous to that of leung2022rate and additionally accounts for cross-cluster dependence since $C_{c(i)} \subseteq \Lambda_i$. In (ref), we discuss the advantages of taking the larger of the two.
We next derive bounds on the rates of convergence of our estimator and difference in means. We then characterize the asymptotic distribution of our estimator and prove that the variance estimator is asymptotically conservative.
As discussed in (ref), this ensures overlap or positivity. For the next two theorems, abbreviate the estimand $\theta^*(d_1,d_0,p_1,p_0)$ as $\theta^*$. Recall the asymptotic order notation from the end of (ref).
The result establishes a bias-variance trade-off in $k_n$. The $k_n^{-1/2}$ term is the contribution of (the square root of) the variance since the effective sample size in a CRT is the number of clusters $k_n$. The asymptotic bias is order $r_n^{-\gamma}$. Intuitively, if $r_n$ is large, then the expected outcome of a unit $i$ assigned treatment $d_t$ and well surrounded by clusters assigned to saturation level $p_t$ should well approximate ${\bf E}^*_{p_t}[Y_i(d_t,\bm{D}_{-i})]$, corresponding to lower bias. The bias decreases at a faster rate for larger $\gamma$ since this corresponds to a lower degree of spatial interference.
The requirement $k_n \precsim n/\xi_n$ is mild, and typically we would have $k_n \prec n/\xi_n$. In the increasing-domain case where $\xi_n=1$, necessarily $k_n \precsim n/\xi_n$ since $k_n \leq n$, and usually the number of clusters is of smaller order than the population size. In the infill-increasing case, if $k_n$ grows faster than $n/\xi_n$, which is the volume of the study region under (ref), then we would be increasingly subdividing the region into smaller clusters of shrinking volume, analogous to having more clusters than units.
Denote the difference-in-means estimator of $\theta^*$ by $\hat\theta^+$, which corresponds to setting $r_n=0$ in the definition of $\hat\theta$. Let $\hat\theta_D^+$, $\hat\theta_T^+$, and $\hat\theta_O^+$ be the difference-in-means estimates of $\theta_D^*$, $\theta_T^*$, and $\theta_O^*$ defined in (ref).
Part (a) provides an upper bound on the rate. Part (b) shows that the rate is tight for several illustrative cases. Recall that setting $p_0=0$ corresponds to the “pure control” baseline used in the estimands of hayes2017cluster.
The $k_n^{-1/2}$ term in the bound is (the square root of) the variance contribution, as in (ref), while $(k_n\xi_n/n)^{1/d}$ is the bias contribution. The latter is notably worse than that of $\hat\theta$ because a reduction in the degree of interference $1/\gamma$ has no effect on the rate. The reason is that units situated at cluster boundaries may be directly proximate to clusters assigned to different treatment arms, and as shown in the proof, the share of such units can be of order $(k_n\xi_n/n)^{1/d}$.
These results demonstrate that excluding units in the manner of $\hat\theta$ can significantly reduce the asymptotic order of the bias relative to difference in means. While excluding units may come at the cost of efficiency, there is no increase to the asymptotic order of the variance because the variances of both estimators scale not with the number of units but with the number of clusters. Hence, the efficiency loss is second-order relative to the potential reduction in bias.
Define $\mu_t = n^{-1} \sum_{i\in\mathcal{N}_n} {\bf E}[Y_i \mid T_{ti}=1]$, $\bar{\theta} = \mu_1 - \mu_0$, and
Note that $\sigma_n^2 \precsim 1$ by the proof of (ref).
The first result (ref) centers the estimator at $\bar{\theta}$, the probability limit of $\hat\theta$. This is not the estimand of interest since there is an asymptotic bias $\lvert\bar{\theta} - \theta^*\rvert$ by Theorem (ref). To use the normal limit to justify the validity of conventional CIs for $\theta^*$, the second result (ref) requires “undersmoothed designs” in which the number of clusters is of smaller order than the optimal rate discussed in the next section. This ensures that the asymptotic bias is small. It is analogous to nonparametric regression where rate-optimal tuning parameter choices result in asymptotic bias, so conventional CIs require undersmoothing. We discuss practical choices of $k_n$ in (ref). In (ref), we discuss “bias-aware” inference, an alternative to undersmoothing.
The cluster-robust variance estimator $\hat\sigma^2(2)$ has the advantage of being non-negative in finite sample, unlike $\hat\sigma^2(1)$ which is a truncation estimator andrews1991heteroskedasticity. On the other hand, $\hat\sigma^2(2)$ only accounts for within-cluster dependence, whereas $\hat\sigma^2(1)$ also accounts for cross-cluster dependence in the definition $A_{ij}(1)$. As shown in the proof, this dependence vanishes, so $\hat\sigma^2(1)$ is a valid estimator. However, the dependence vanishes at the slow rate $(k_n\xi_n/n)^{1/d}$, the same order as the bias of difference in means. Thus in smaller samples, ignoring second-order terms may result in anti-conservativeness. By taking the larger of the two estimators, we obtain the benefits of both. See (ref) for a comparison of $\hat\sigma^2$ with other variance estimators in the literature.
Recall that $d$ is the dimension of the spatial region, $\xi_n$ is the density of the region (number of units per unit volume/area), and $\gamma$ is a lower bound on the speed at which interference decays with distance. By (ref), choosing
optimizes the rate of convergence of $\hat\theta$. This formalizes how domain knowledge of interference $\gamma$ informs the design of clusters. The right-hand side is increasing in $\gamma$ since less interference means the same level of bias reduction can be achieved with more clusters. It is decreasing in density $\xi_n$ because having more units in a given area effectively corresponds to greater interference (an infectious disease may spread more easily).
Choosing the number of clusters according to (ref) balances the asymptotic orders of the bias and variance of $\hat\theta$. However, as discussed in (ref), validity of the CI
requires the bias to be of smaller order than the variance. This requires an {\em undersmoothed design} where $k$ is of smaller order than (ref).
We next recommend an undersmoothed choice of $k$ given domain knowledge of $\gamma$. Let $\mathcal{V}$ denote the area or volume of the study region containing $\mathcal{N}_n$. As discussed below, $\mathcal{V}$ is of order $n/\xi_n$, so the rate-optimal formula (ref) can be rewritten as $k \sim \mathcal{V}^{\frac{2\gamma}{2\gamma+d}}$. While one could choose $k$ equal to the right-hand side, this choice would be extremely sensitive to the unit of length used to measure geographic distance. Switching from square kilometers to square millimeters would dramatically increase $\mathcal{V}$.
Our observation is that both the unit of length and $\gamma$ determine the speed at which interference decays with distance. For any given choice of $\gamma$, the rate of decay is significantly faster if distance is measured in millimeters compared to kilometers, for example. Therefore, domain knowledge of interference informs both $\gamma$ and the unit of length, and one can determine $k$ from these ingredients as follows.
The next subsection provides empirical examples of calibrating $k$ under the conservative choice $\tilde\gamma=d$. The last subsection discusses how to determine $\gamma$ to obtain less conservative estimates.
sur2009cluster conduct a CRT in an urban slum in India spanning about 1.2 by 0.7 km with 38k participants. To compute our suggested number of clusters (ref), we need to select $\tilde\gamma$ and the unit of length. Suppose we choose the conservative bound $\tilde\gamma = d = 2$, so that interference is assumed to decay like $r^{-2}$ or faster with each unit of length $r$. If we take 35m to be the unit of length, then for units $i,j,k$ such that $\rho(i,j) = 35$m and $\rho(i,k) = 70$m, the extent to which $k$'s treatment affects $i$ is less than $1/4$th ($2^{-\tilde\gamma}=0.25$) as small as the extent to which $j$'s affects $i$. This is with only a 35m difference in distance. For this unit of length, (ref) yields $k = 78$. In other words, these are the assumptions on interference that justify the authors' choice of $k=80$, originally determined by a conventional power analysis based on a partial interference model.
The previous rate of decay may be over-optimistic, so suppose the relevant unit of length is in 100m increments. If $\rho(i,j)=100$m and $\rho(i,k) = 200$m, the extent to which $k$'s treatment affects $i$ is less than $1/4$th as small as the extent to which $j$'s affects $i$, now with a 100m difference in distance. With a slower rate of decay, bias is higher, which requires constructing fewer, larger clusters. As a result, (ref) yields only $k = 19$.
Next consider homan2016effect whose trial region is substantially larger at roughly 12 by 4 km with nearly the same number of units (34k). Due to the lower density, there is less bias from interference, so a choice of $k=81$ can be justified under weaker assumptions on interference. Specifically, if $\tilde\gamma=2$ but the unit of length is now in 250m increments, then (ref) results in $k=84$. For context, the {\em Anopheles} mosquito that is the subject of their trial typically does not fly more than 2 km from their breeding grounds anopheles.
These examples illustrate what domain knowledge of the degree of interference entails. They also show how our proposed choice of $k$ accounts for regional density, unlike the standard power analysis hemming2011sample. In both examples, weaker assumptions on interference require choosing $k$ smaller than the standard analysis to better balance bias and variance. To justify choosing larger $k$, the researcher must either collect data over a larger area or be willing to entertain stronger assumptions on interference.
Optimal design generally requires prior information on certain population parameters. The standard power analysis for CRTs assumes partial interference and requires knowledge of the intracluster correlation coefficient, a measure of within-cluster dependence in potential outcomes baird2018optimal,hemming2017design. Our results suggest that, if cross-cluster interference is of first-order importance, the focus of attention should instead be $\gamma$. As illustrated in the previous subsection, the conservative choice $\gamma=d$ can result in relatively small values of $k$ because it allows for a greater degree of interference, so to the extent that one can justify stronger assumptions on interference, that is, values of $\gamma$ that are larger than $d$ for a given unit of length, this would improve asymptotic power.
In the context of infectious disease trials, halloran2017simulations argue that CRT design should be informed by simulating models of disease transmission. Several papers utilize parametric models and simulation methods to estimate or bound contamination bias. alexander2020spatial and jarvis2019spatial use spatial models to provide evidence of contamination in prior CRTs. multerer2021analysis employ models of disease transmission for a similar purpose. Our theory provides a precise way in which modeling can inform design, namely by providing plausible values of $\gamma$. This relates to staples2015incorporating who show how to estimate a different measure of cross-cluster interference to assess the degree to which the conventional power analysis overstates trial power.
To estimate $\gamma$, models may be combined with external data sources such as data from pilot studies. In the context of malaria vector control, mosquito mark-release-recapture experiments guerra2014global provide data on geographic dispersion of malaria carriers, which is informative of spatial interference. We provide additional suggestions in (ref).
We conduct a simulation study to illustrate the finite-sample properties of our estimator and difference in means under our proposed design. We randomly draw unit locations from the square $[-(n\alpha_n)^{1/2}, (n\alpha_n)^{1/2}]^2$ with $\alpha_n = 0.8, 0.7, 0.6$, respectively. This corresponds to the infill-increasing case where the regional volume shrinks with the population size. We create clusters using $k$-medoids with $k$ given by (ref) using the conservative choice $\tilde\gamma=d=2$ and set the assignment probabilities in (ref) to $(q,p_1,p_0) = (0.7,0.5,0)$.
Let $\{\tilde\varepsilon_i\}_{i\in\mathcal{N}_n} \stackrel{iid}\sim \mathcal{N}(-0.5,1)$, $\{\beta_i\}_{i \in \mathcal{N}_n} \stackrel{iid} \sim \mathcal{N}(2,1)$, and $\{\gamma_i\}_{i \in \mathcal{N}_n} \stackrel{iid} \sim \mathcal{N}(1,1)$ be independent and drawn independently of locations. We generate spatially autocorrelated errors $\varepsilon_i = \tilde\varepsilon_i + \sum_{j\in\mathcal{N}_n} \bm{1}\{\rho(i,j) \leq 1\}\tilde\varepsilon_j / \sum_{k\in\mathcal{N}_n} \bm{1}\{\rho(i,k) \leq 1\}$. For $w_{ij} = \min\{\rho(i,j)^{-5}, 1\}$, we generate outcomes according to
Under this model, the unit-level direct and indirect effects are respectively given by
Due to the choice of $-5$ in the spatial weights $w_{ij}$, (ref) holds for $\gamma = 3$ leung2022rate, which is a fairly slow rate of decay given that (ref) requires $\gamma>2$.
We present results for the indirect and total effects using 5000 simulation draws where within each draw, we redraw potential outcomes and recompute the design-based estimand. In (ref), the “Spatial Interference” columns correspond to the outcome model described above, whereas the “Partial Interference” columns redefine $w_{ij}=0$ if $i,j$ lie in different clusters to eliminate cross-cluster interference. The “CI” rows report the coverage of 95-percent CIs using the indicated standard errors. For our estimator, the “SE” row corresponds to standard errors obtained from our variance estimator (ref), while for difference in means, it corresponds to conventional cluster-robust SEs. The “SE$^*$” rows are the true superpopulation standard errors obtained by taking the standard deviation of the estimator across the simulation draws. As such, SE should be consistent for SE$^*$ rather than conservative, while the CIs have asymptotically conservative coverage for the design-based estimand. Finally, row “% Excl” is the percentage of units that are not well surrounded.
The results are consistent with the theory. As $n$ grows, the biases of our estimators shrink, while coverage tends to or exceeds the nominal level. The bias of difference in means is more than twice that of our estimators under spatial interference, resulting in severe undercoverage even with the true superpopulation SEs. While the SEs are smaller than those of our estimators, this is not by a significant amount, and the efficiency advantage comes at a large cost to bias under spatial interference.
(ref) presents results for our estimator under spatial interference, but we multiply (ref) by $c \in \{0.8,1.2\}$. This is to assess robustness and illustrate a bias-variance trade-off in the choice of $r_n$. The $c=0.8$ columns show that this results in about half the proportion of excluded units relative to (ref). As a result, the bias is higher, resulting in undercoverage. The $c=1.2$ columns show that the proportion of excluded units a little less than doubles. The bias is lower, and as a result, the probability of coverage is higher. On the other hand, the variance increases, as can be seen in the standard error columns.
We apply our estimator to data from the unconditional cash transfer experiment mentioned in (ref). In the experiment, households eligible for the transfer (the treatment) live in homes with thatched roofs, which is a proxy for poverty. Houses are grouped in villages, which are grouped in “sublocations.” See Figure A.2 of egger2022general for a map of the study area. The CRT randomizes sublocations (the clusters) into treatment with probability $q=0.5$. Within treatment (control) sublocations, villages are randomized into treatment with probability $p_1=2/3$ ($p_0=1/3$). Within treated villages, all eligible households receive cash transfers totaling 1000 USD, which is about 75 percent of average annual household spending. Sublocations contain 7.8 villages on average (SD 3.9).
egger2022general study the effect of the cash transfers on the following household-level outcomes, which can be grouped into three categories: (1) annualized expenditures on consumption (of food and other purchases described in their footnote 31), non-durables, food alone, temptation goods, and durables; (2) asset stocks, housing value, and land value; and (3) annualized household income, transfers, taxes paid, profits, and wage earnings. The authors report model-assisted estimates of causal effects at the household level. Our analysis will be entirely design-based but at the village level. We aggregate household outcomes to the village level by averaging.
egger2022general estimate the effects of the transfers on the population of eligible households using two main specifications described in their \S 3.2. Their “RF” (reduced form) specification is an OLS regression of an outcome on village- and sublocation-level treatment indicators and covariates. They report the coefficient on the village-level indicator, which corresponds to a model-assisted estimate of a village-level direct effect $\theta^*(1,0,2/3,2/3)$ among eligible households.
However, the authors note the potential for interference across sublocations (their quote in (ref)). For this reason, their preferred specification is the following “IV” (instrumental variables) regression. The main regressors are the amount of cash transferred per capita to the household's village and the amount transferred to neighboring villages within a band of $r-2$ to $r$ km from the ego's village for a range of $r$ values. The corresponding instruments are respectively an indicator for the ego's village being treated and the share of eligible households in the band assigned to treatment. They use a BIC criterion to select the maximum range of $r$, which is 2 km. Using these estimates, the authors compute a model-assisted estimate of the total effect $\theta^*(1,0,2/3,0)$. This compares saturation levels of two-thirds and zero, which is nonparametrically unidentified since control villages have a saturation level of one-third. We instead report our estimates of the total effect $\theta^*(1,0,2/3,1/3)$.
(ref) reports the results for a subset of the outcomes, and (ref) in the supplementary appendix reports the remainder. Both tables choose $r_n$ in $\hat\theta$ according to (ref), which results in $r_n=1.6$ and 39.66 percent of units not well surrounded. In (ref) we describe how we construct cluster radii $R_j$ used in this formula. The $\hat\theta^+$ columns correspond to difference-in-means estimates with clustered standard errors.
We find that the difference-in-means estimates of the direct and total effects are comparable to the RF and IV estimates, respectively, despite the distinctions outlined above. Our estimators find larger direct and total effects. Decomposing the total effect into direct and indirect effects, we find that the latter are substantially smaller in magnitude with large standard errors. Compared to difference in means, our estimates tend to be larger in magnitude but with larger standard errors due to the restriction to well-surrounded units.
(ref) in the supplementary appendix reports results for $r_n=2$, resulting 57.89 percent of units not being well surrounded. This choice of $r_n$ coincides with the largest range of $r$ selected by the BIC procedure of egger2022general. Our estimates and standard errors become larger still in magnitude relative to difference in means, but the results are qualitatively similar.
To estimate spillover effects, egger2022general rerun their IV specification using only non-eligible households, which did not receive any transfers (their \S 3.3). In (ref) of the supplementary appendix, we compare their results with design-based estimates of the overall effect $\theta^*(\emptyset, \emptyset, 2/3, 1/3)$ on the population of non-eligibles. The effect sizes of our estimators and those of difference in means are fairly similar in magnitude to their IV results, but the standard errors are large. Combined with the results in (ref), we ultimately find strong direct effects of the cash transfers but weaker evidence for spillover effects compared to egger2022general.
When interference occurs across clusters, conventional analyses of CRTs suffer from bias induced by units near cluster boundaries. To reduce bias at the analysis stage, we provide in (ref) an estimator $\hat\theta$ that improves upon the fried-egg design by excluding from estimation units that are not surrounded by clusters assigned to the same treatment arm. This may be employed as a robustness check for difference in means. To reduce bias at the design stage, we propose a rate-optimal formula for the number of clusters $k$ in (ref). Unlike the standard power analysis that assumes partial interference, our choice balances power against the need to reduce bias and accounts for the density of units in the spatial region. Given $k$, we suggest automating cluster construction using $k$-medoids and prove that the resulting clusters are balanced and globular, thereby approximately minimizing the number of units near boundaries. We also provide valid design-based standard errors.
Under the conventional superpopulation, partial interference framework, power calculations for choosing $k$ require prior knowledge of the intracluster correlation coefficient (ICC). Under our design-based framework, the optimal choice of $k$ requires knowledge of $\gamma$, the speed at which interference decays with distance, rather than the ICC. Absent domain knowledge, one can conservatively set $\gamma$ to the dimension of the spatial region. We discuss in (ref) how to obtain less conservative estimates via modeling or prior data.
\part{Supplementary Appendix}
\makeatletter \@addtoreset{section}{part} \makeatother \setcounter{section}{0} \numberwithin{equation}{section} \numberwithin{table}{section}
The $k$-medoids algorithm selects a set of $k$ units in $\mathcal{N}_n$ -- the {\em medoids} or cluster centroids -- to minimize the total distance between all units and their nearest medoids. This creates medoids that are spatially well separated, as shown in (ref). Units are then grouped into clusters based on their closest medoids.
Whereas $k$-means allows centroids to be any element of $\mathcal{X}$, $k$-medoids constrains the centroids to the data $\mathcal{N}_n \subseteq \mathcal{X}$, but in practice the algorithms tend to produce similar output. Our theoretical results pertain to the implementation in (ref), the “partitioning around medoids” algorithm. This is simple to understand, though it has complexity $O(k(n-k)^2)$ compared to an $O(n^2)$ runtime under the fastest known implementation schubert2021fast.
For any set of candidate medoids $\mathcal{M} \subseteq \mathcal{N}_n$, let $\text{cost}(\mathcal{M}) = \sum_{i \in \mathcal{N}_n} \min_{m \in \mathcal{M}} \rho(i,m)$, the total distance between units and their nearest medoids.
In other words, this chooses the set of medoids $\mathcal{M}$ that minimizes cost by iteratively swapping out a candidate medoid with a better unit that reduces cost.
As discussed in (ref), (ref)(a) restricts the extent to which density can vary across the study region. To accommodate larger variation in density, such as when the study region encompasses both urban and rural areas, we suggest first partitioning the region $\mathcal{N}_n$ into subregions that are relatively homogeneous in density. This may be done manually, for instance using cartographic and demographic information. It can also be done using a density-based clustering algorithm; see bhattacharjee2021survey and especially kriegel2011density for surveys of this literature.
Call the density-homogeneous subregions $\mathcal{S}_1, \ldots, \mathcal{S}_m$. For each $\mathcal{S}_j$ the researcher can compute the optimal number of clusters $k_j$ given in (ref), subdivide $\mathcal{S}_j$ into $k_j$ clusters $C_{j1}, \ldots, C_{jk_j}$, say using $k$-medoids, and cluster-randomize across the collection of all clusters $\{C_{j\ell}\colon \ell=1,\ldots k_j, j=1,\ldots, m\}$. Lastly, we recommend modifying the formula for $r_n$ since cluster radii may vary substantially across subregions with differing densities. In the definition of $T_{ti}$, replace $r_n$ with half the median cluster radius among clusters in the subregion containing $i$.
Estimating the following spatial moving average model with prior data can be a starting point for determining $\gamma$:
This satisfies (ref) with $\gamma = \eta-2$ leung2022rate, so $\gamma$ can be backed out from an estimate of $\eta$.
It would also be useful to develop designs for nonparametrically estimating $\gamma$. A preliminary proposal is the following. Given $k$ clusters, say constructed using $k$-medoids, denote by $m_j$ the centroid of a cluster $C_j$. Randomize clusters to the following $T$ treatment arms. In the 0th arm, we assign all units to control. In the $t$th treatment arm for $t\geq 1$, we only assign units $i$ for which $\rho(i,m) \in (t, t+1]$ to treatment, where $m$ is the centroid of the cluster in question. That is, we treat only units in a ring at a certain distance from the centroid. Let $\hat\theta_t$ denote the difference-in-means estimate comparing average outcomes of units near centroids of clusters assigned to arm $t\geq 1$ with those assigned to arm 0. Then under (ref), $\hat\theta_t$ should decay like $t^{-\gamma}$, producing a sort of “causal covariogram.” We may then regress $\log \hat\theta_t$ on $\log t$ to estimate $-\gamma$.
Recall from (ref) that the “undersmoothed” design chooses $k_n$ smaller than the optimal rate, which is why (ref) uses a strict lower bound $\tilde\gamma$ in place of $\gamma$. By analogy to nonparametric regression, an alternative is to choose $k$ rate-optimally using (ref) with $\gamma$ in place of $\tilde\gamma$ and use the following “bias-aware” confidence interval (CI) in place of (ref):
where $c$ is given in (ref). faridani2023rate propose a similar bias-aware CI for the GATE. Compared to (ref), (ref) should have better finite-sample coverage, although its implementation requires prior knowledge of both $\gamma$ and $c$.
In addition to knowledge of these parameters, suppose the researcher has preliminary consistent estimates of $\hat\theta$ and $\hat\sigma^2$, say from a pilot study (which is more plausible in a superpopulation setup). Then following armstrong2018optimal, they can pick $k$ to minimize the length of the bias-aware CI (ref). Note that $r_n$ is an implicit function of $k$, which can be traced out using grid search.
The idea behind (ref) is that, by (ref),
As shown in the proof of (ref), specifically the argument preceding (ref),
This provides a worst-case bound on the bias that (ref) incorporates.
We next discuss how $\hat\sigma^2$ relates to variance estimators in the literature. Define $\bm{A}(1)$ ($\bm{A}(2)$) as the symmetric, $n\times n$ matrix with $ij$th entry $A_{ij}(1)$ ($A_{ij}(2)$), and recall that $\hat\sigma^2(2)$ is the cluster-robust variance estimator, while $\hat\sigma^2(1)$ is analogous to the leung2022rate variance estimator. The proof of (ref) shows that the cluster-robust variance estimator can be asymptotically decomposed as $\hat\sigma^2(2) = \sigma_n^2 + \mathcal{B}_n + o_p(1)$ where
and $\bm{\mu}$ is the $n$-dimensional vector with $i$th component $(\mu_{1i}-\mu_{0i}) - (\mu_1 - \mu_0)$ for $\mu_{ti} = {\bf E}[Y_i \mid T_{ti}=1]$. Because $\bm{A}(2)$ is block-diagonal, it is positive semidefinite, so $\mathcal{B}_n \geq 0$ for any $n$. Hence, $\hat\sigma^2(2)$ is asymptotically conservative.
The proof further shows that $\hat\sigma^2(1) = \sigma_n^2 + \mathcal{B}_n(1) + o_p(1)$ where $\mathcal{B}_n(1)$ is obtained by replacing $\bm{A}(2)$ with $\bm{A}(1)$ in (ref). This replacement adds nonzero off-diagonal elements to $\bm{A}(2)$, so $\bm{A}(1)$ is not block-diagonal and hence not guaranteed to be positive semidefinite. Accordingly, Leung imposes additional conditions to show that $\mathcal{B}_n(1) \stackrel{p}\longrightarrow c \geq 0$ in the superpopulation. We find, however, that $\mathcal{B}_n(1) = \mathcal{B}_n + o_p(1)$ without any conditions on the superpopulation, so in fact $\hat\sigma^2(1)$ is asymptotically conservative in a purely design-based setup. faridani2023rate provide a different approach to conservative inference for the GATE based on bounding $\mathcal{B}_n$.
The positive-semidefiniteness of the “kernel” $\bm{A}(2)$ results in both a non-negative variance estimator for any $n$ and asymptotic conservativeness. This mirrors the corresponding insight for HAC variance estimators, that positive semidefinite kernels ensure conservativeness in finite-population models. This fact was first pointed out by leung2019causal and has since been exploited by other papers. (In his case of network-dependent data, the difficulty with using positive-semidefinite HAC kernels is that they are sloped and tend to severely over-reject in finite samples, which is why leung2022causal recommends use of the uniform kernel even though it is not positive semidefinite. For spatial data, a variety of positive semidefinite HAC kernels exist, and these tend to have less severe issues with over-rejection relative to the network case.)
hudgens2008toward propose variance estimators for difference-in-means type estimators under partial interference. Their theory relies on an additional stratified interference assumption, which says that potential outcomes only depend on the ego's treatment assignment and the proportion of treated units in the ego's cluster, in which case $Y_i = Y_i(D_i, p)$ (since they assume complete randomization). Because $p$ is a constant, the dependence structure is the same as under no interference since outcomes within a cluster are only correlated if treatment assignments are. The variance estimators proposed by hudgens2008toward heavily rely on this structure. In contrast, we do not impose partial, let alone stratified, interference, so our setting features both within- and cross-cluster dependence in outcomes. It is therefore critical to account for additional covariance terms that are absent in the hudgens2008toward variance formula to avoid anti-conservativeness. Our estimator does so through the terms involving $A_{ij}(u)$ for $i\neq j$.
In the econometric literature, the standard approach under partial interference is clustering standard errors baird2018optimal. These allow for arbitrary within-cluster dependence and are valid even in the absence of stratified interference. A new finding of our paper is that clustered standard errors are also valid in a design-based setting with cross-cluster interference. However, we suggest combining them with $\hat\sigma^2(1)$ to better capture second-order covariance terms, as discussed in (ref).
We define $R_j$ from (ref) as half the largest distance between households in any pair of villages within cluster $j$. To compute this, we utilize supplemental data provided by Dennis Egger on distances between village centroids and publicly available data from egger2022general to estimate village radii. We first compute the largest distance between the centroids of any pair of villages $a,b$ within a given cluster, which we denote by $\Delta_{a,b}$. Since this does not account for village size, we estimate for each village $a$ its radius $\delta_a$, which is the furthest distance between a household in the village and its centroid. This data is publicly available in 1 km increments. For each cluster $C_j$, we define $R_j = \max_{a,b \in C_j} (\Delta_{a,b} + \delta_a + \delta_b)/2$. The average radius across clusters is 3.12 with a standard deviation of 0.81.
The first two lemmas establish properties of $k$-medoids clusters for $k=k_n$ possibly diverging.
Let $T_{ti}^+$ equal $T_{ti}$ with $r_n$ set to zero, $\hat{p}_t^+ = n^{-1} \sum_{i\in\mathcal{N}_n} T_{ti}^+$, and $p_t^+ = {\bf E}[T_{ti}^+]$, which does not depend on $i$ under (ref).
Let each cluster $C_n$'s “centroid” $i_n$ be its medoid and $U_n \equiv \max_j R_j$ for $R_j$ defined prior to (ref), so $C_n \subseteq \mathcal{N}(i_n, U_n)$. By (ref), $U_n \precsim (n/(k_n\xi_n))^{1/d}$. By (ref), there exists $\Delta_n \succsim (n/(k_n\xi_n))^{1/d}$ such that any pair of medoids is physically separated by at least $\Delta_n$. Then by definition of $k$-medoid clusters, for $L_n = \Delta_n/2$, $\mathcal{N}(i_n,L_n) \subseteq C_n$. \rule{0.08in}{0.08in}
Define the following Horvitz-Thompson analog of $\hat\theta$:
This is unbiased for $\bar{\theta}$ defined prior to (ref). Let $\hat\kappa_t = n^{-1} \sum_{i\in\mathcal{N}_n} T_{ti}/p_{ti}$. By (ref) and (ref), $\tilde Z_i$ is uniformly asymptotically bounded, so
by (ref). It remains to bound the bias and variance of $\tilde\theta$.
For any $i\in\mathcal{N}_n$, $d_t \in \{0,1\}$, $\bm{d} \in \{0,1\}^n$, and $S \subseteq \mathcal{N}_n$ containing $i$, we use the notation $(d_t, \bm{D}_{S\backslash\{i\}}, \bm{d}_{-S})$ to mean that we take the observed treatment vector $\bm{D}$, replace entry $D_i$ with $d_t$, and replace the subvector $\bm{D}_{\mathcal{N}_n\backslash S} = (D_j\colon j\in\mathcal{N}_n\backslash S)$ with the corresponding entries of $\bm{d}_{\mathcal{N}_n\backslash S}$.
The basic idea is that if $T_{ti}=1$, then $D_i=d_t$ and $i$ is well surrounded by units belonging to clusters with saturation level $p_t$, so ${\bf E}[Y_i \mid T_{ti}=1]$ is a good approximation of ${\bf E}^*_{p_t}[Y_i(d_t,\bm{D}_{-i})]$. Formally,
because
by (ref) and
by (ref). Therefore,
The last part of (ref) holds because $r_n$ is the median cluster radius by (ref), and cluster radii are uniformly of order $(n/(k_n\xi_n))^{1/d}$ by (ref).
Recalling the definition of $\Lambda_i$ from (ref),
By (ref), $\max_i \lvert\Lambda_i\rvert \precsim n/k_n$, so $[P1] \precsim k_n^{-1}$. It remains to show that $[P2] \precsim k_n^{-1}$. Intuitively, $[P2]$ should be well controlled since units not in $\Lambda_i$ are “far” from and therefore less correlated with $i$.
{\bf Step 1.} Recall that $c(i)$ is the index of the cluster containing unit $i$. Define $X_i^r = {\bf E}[\tilde Z_i \mid \mathcal{F}_i(r)]$ for
We bound the discrepancy between $\tilde Z_i$ and $X_i^r$. For any $r\geq r_n$ and $t\in\{0,1\}$, $T_{ti}$ is measurable with respect to the $\sigma$-algebra generated by $\mathcal{F}_i(r)$. Then by (ref), for any $q>0$, $r\geq r_n$, $n$ sufficiently large, and $t \in \{0,1\}$,
using an argument similar to (ref). By (ref), $\max_i p_{ti}^{-1}\precsim 1$, so by the law of total probability, there exists $c_n \precsim 1$ such that
{\bf Step 2.} Fix $i,j\in\mathcal{N}_n$ such that $j\not\in\Lambda_i$. We use (ref) to bound $\text{Cov}(\tilde Z_i, \tilde Z_j)$. Since the $r_n$-neighborhoods of $i$ and $j$ do not intersect a common cluster, $X_i^{r_n} \perp\!\!\!\perp X_j^{r_n}$ by (ref). Applying the Cauchy-Schwarz inequality, (ref) for $q=2$, and (ref) and (ref), there exists $c_n' \precsim 1$ such that for $n$ sufficiently large and all $i,j$ such that $j\not\in\Lambda_i$,
Now suppose additionally that $\rho(i,j) > 4\bar{R}$ for
where $R_j$ is defined prior to (ref). In this case, we derive a different covariance bound. Since $\mathcal{N}(i,\rho(i,j)/2-\bar{R})$ and $\mathcal{N}(j,\rho(i,j)/2-\bar{R}))$ are separated by a distance of at least $2\bar{R}$, which upper bounds the “diameter” of any cluster, we have $X_i^{\rho(i,j)/2-\bar{R}} \perp\!\!\!\perp X_j^{\rho(i,j)/2-\bar{R}}$. By a derivation similar to (ref) using $\rho(i,j)/2-\bar{R}$ in place of $r_n$,
{\bf Step 3.} Let $\lceil c \rceil$ ($\lfloor c \rfloor$) denote $c$ rounded up (down) to the nearest integer. Using the covariance bounds derived in step 2,
where $[P2.1]$ takes the part involving $s \leq 4\bar{R}$ and $[P2.2]$ the part involving $s > 4\bar{R}$. The sum over $s$ starts at $\lfloor 2r_n \rfloor$ because $j\not\in\Lambda_i$ implies that the $r_n$ neighborhoods of $i,j$ do not intersect.
We next show that $\eqref{covzz} \precsim \xi_n/n$. By Assumptions (ref)(a) and (ref), $\max_\ell \lvertC_\ell\rvert \precsim n/k_n$. By (ref)(b), $\sum_{j\not\in\Lambda_i} \bm{1}\{\rho(i,j)\in [s,s+1)\} \leq \lvert\mathcal{N}(i,s+1)\backslash\mathcal{N}(i,s)\rvert \leq C \xi_n \max\{s^{d-1},1\}$. Then
because
given $\bar{R} \geq r_n$ by definition. By (ref),
and $\bar{R}^d \precsim n/(k_n\xi_n)$, so
since $\gamma>d$ by (ref). Similarly, using the fact that $d \geq 1$ from (ref)(b),
Combining (ref), (ref), (ref), and (ref) yields $[P2] \precsim \xi_n/n \precsim k_n^{-1}$ since $k_n\xi_n/n \precsim 1$ by assumption. \rule{0.08in}{0.08in}
Let $T_{ti}^+ = \bm{1}\{D_i=d_t, W_{c(i)}=t\}$ for $d_t \in \{0,1\}$ and $T_{ti}^+ = \bm{1}\{W_{c(i)}=t\}$ for $d_t = \emptyset$. Define $p_t^+ = {\bf E}[T_{ti}^+]$, $\hat p_t^+ = n^{-1} \sum_{i\in\mathcal{N}_n} T_{ti}^+$, and
the Horvitz-Thompson analog of the difference in means $\hat\theta^+$. We have
because $T_{ti}^+Y_i/p_t^+$ is uniformly asymptotically bounded by (ref) and (ref), and $p_t^+ - \hat{p}_t^+ \precsim k_n^{-1/2}$ by (ref). It suffices to derive the rate of convergence of $\tilde\theta^+$.
By the argument used to bound (ref) in the proof of (ref), $\text{Var}(\tilde\theta^+) \precsim k_n^{-1}$. It remains to bound the bias, that is, to show that
Recall the definitions in (ref), letting $m_j$ denote the “centroid” of cluster $C_j$. Define $j$'s “boundary” as
By (ref), there exists $\alpha_n \precsim 1$ such that $U_n - L_n = \alpha_n (n/(k_n\xi_n))^{1/d}$, so by (ref)(b),
For $r = 0, \ldots, L_n$, define the “contour sets”
where $\mathcal{N}(i,-1) \equiv \emptyset$. As $r$ increases, $J(r,C_j)$ moves away from the boundary and towards the interior. If $i \in J(r,C_j)$, then $\mathcal{N}(i,r) \subseteq C_j$, so for such $i$, by Assumptions (ref) and (ref),
for any $t\in\{0,1\}$. Then
The second line uses (ref) and (ref). It converts the sum over all units to a sum over all clusters followed by sums over units in each contour set through the boundary. The third line bound on $\lvert\mathcal{B}(C_j)\rvert$ follows from (ref), while the bound on $\lvertJ(r,C_j)\rvert$ uses (ref)(b). This establishes (ref).
Let $d\geq 1$, and fix any sequence $\{\xi_n\}_{n\in\mathbb{N}}$ such that $1 \precsim \xi_n \prec n$. Let $\rho$ be the sup norm and $\mathcal{N}_n = \{\xi_n^{-1/d}x\colon x \in \mathbb{Z}^d\} \cap B(\bm{0},\mathcal{R}_n)$ where $B(\bm{0},\mathcal{R}_n)$ is a (hyper)cube with radius $\mathcal{R}_n \in \mathbb{Z}$ and $\mathcal{R}_n \sim (n/\xi_n)^{1/d}$. When $\xi_n\sim 1$, the observed units are positioned within a cube of radius $\sim n^{1/d}$ containing $\sim n$ units positioned on the integer lattice. When $\xi_n \succ 1$, we shrink this region towards the origin by a factor $\xi_n^{-1/d}$. This results in the same asymptotic order of number of units $\lvert\mathcal{N}_n\rvert \sim n$, but the volume of the cube is reduced to $\sim \mathcal{R}_n^d \sim n/\xi_n$, resulting in a density of $\xi_n$.
{\bf Verifying (ref).} By construction of $\mathcal{N}_n$, there exist $C_0,C_1>0$ such that for $r\geq 0$ and $i\in\mathcal{N}_n$,
Substituting $\xi_n^{1/d}r$ for $r$ in these expressions,
Letting $[c]$ denote rounding $c$ to the nearest integer, this implies
which is bounded by a constant times $\xi_n r^{d-1}$. This verifies (ref).
{\bf Verifying (ref).} Let $1 \prec k_n \prec n/\xi_n$. For the remainder of the proof, suppose the clusters are generated by $k_n$-medoids which by (ref) satisfy (ref). Also let $n$ be sufficiently large that $L_n$ in the assumption exceeds two.
{\bf Verifying Assumptions (ref) and (ref).} Set
which is unit $i$'s treatment times the fraction of treated units in $i$'s {\em 2-neighborhood}. The sum in the denominator is always at least one since it includes $i$ itself, so (ref) holds. (ref) holds because potential outcomes only depend on 2-neighborhood treatments.
{\bf Reduction to HT.} Notice $Y_i(\bm{0}) = 0$. Since $Q \in \{D,T,O\}$ and $p_0=0$, $T_{0i}^+Y_i = 0$. Then
By (ref), $p_1^+/\hat{p}_1^+ \stackrel{p}\longrightarrow 1$, so it remains to lower bound the bias and variance of $\tilde\theta^+$.
{\bf Bias.} Let $\beta = (p_1-p_1q) \bm{1}\{d_1=1\} + p_1 (p_1-p_1q) \bm{1}\{d_1=\emptyset\}$, which is strictly positive by assumption. We have
Define the boundary of a set $C \subseteq \mathcal{N}_n$ as
the set of units whose $1$-neighborhoods intersect a point outside the set. Then
By construction of $\mathcal{N}_n$,
Therefore,
By (ref),
Let $m_j$ denote the medoid of a cluster $C_j$. By (ref)(a), units in $\tilde{\mathcal{B}}(C_j)$ are at least distance $\lfloor L_n \rfloor$ from $m_j$. Define the “contour set”
By construction of $\mathcal{N}_n$, $\lvert\tilde J(r,C_j)\rvert$ is increasing in $r$, so
uniformly in $j$ by (ref) and (ref)(b). Hence
Combined with (ref), $\eqref{lbias2.5} \succsim (k_n\xi_n/n)^{1/d}$.
{\bf Variance.} For $Z_i^+ = T_{1i}^+ Y_i/p_1^+$, $\text{Cov}(Z_i^+,Z_j^+)$ equals
Since $Q\in\{D,T,O\}$, either $T_{1i}^+ = D_iW_{c(i)}$ in the case of $Q \in \{D,T\}$ or $T_{1i}^+ = W_{c(i)}$ in the case of $Q=O$, so the covariance term equals $\text{Cov}(W_{c(i)}D_iD_\ell, W_{c(j)}D_jD_m) \geq 0$. Furthermore, for $j \in C_{c(i)}$, $W_{c(i)}=W_{c(j)}$, so
If additionally $\ell \in C_{c(i)}\backslash\{i\}$ and $m\in C_{c(i)}\backslash\{j\}$, the covariance term on the right-hand side equals $\alpha \equiv qp_1^4 - (qp_1^2)^2 > 0$. Then by construction of $\mathcal{N}_n$,
By (ref) and (ref)(a), $\min_j \lvertC_j\rvert \geq \lvert\mathcal{N}(m_j,L_n)\rvert \succsim n/k_n$. Hence, the right-hand side of the above display is at least order $k_n^{-1}$. \rule{0.08in}{0.08in}
The first three steps establish (ref), and step four proves (ref).
{\bf Step 1.} We first derive an asymptotically linear representation. For $\hat\kappa_t = n^{-1} \sum_{i\in\mathcal{N}_n} T_{ti} / p_{ti}$ and $\mu_t = n^{-1} \sum_{i\in\mathcal{N}_n} {\bf E}[Y_i \mid T_{ti}=1]$,
By (ref), $\hat\kappa_t - 1 \precsim k_n^{-1/2}$, so
since $n^{-1} \sum_{i\in\mathcal{N}_n} p_{ti}^{-1} T_{ti}(Y_i - \mu_t) \prec 1$ by the variance calculation in the proof of (ref).
{\bf Step 2.} Recall the definition of $\mathcal{F}_i(r)$ from (ref) and $\Lambda_i$ from (ref). Define $B_i = Z_i - {\bf E}[Z_i \mid \mathcal{F}_i(r_n)]$. We next establish that
where the last asymptotic inequality follows from the assumption $k_n \prec n/\xi_n$. Expanding the square yields
For all $r\geq r_n$ and $t\in\{0,1\}$, $T_{ti}$ is measurable with respect to $\mathcal{F}_i(r)$, so
By (ref),
since ${\bf E}[Y_i(\bm{D}_{\mathcal{N}(i,r)}, \bm{0}_{-\mathcal{N}(i,r)}) \mid \mathcal{F}_i(r)] = Y_i(\bm{D}_{\mathcal{N}(i,r)}, \bm{0}_{-\mathcal{N}(i,r)})$. Then by (ref), there exists $c_n \precsim 1$ such that for all $i$,
By (ref), $\max_i \lvert\Lambda_i\rvert \precsim n/k_n$, so (ref) and (ref) imply
which is $\precsim k_n\xi_n/n$ since $\gamma>d$ by (ref).
Turning to $\lvert[P2]\rvert$, fix $i,j$ such that $j\not\in\Lambda_i$. By (ref), there exists $c_n \precsim 1$ such that for all such $i,j$,
We will use a different bound when additionally $\rho(i,j) > 4\bar{R}$ for $\bar{R} = \max_j R_j$. In this case we have
since $\mathcal{N}(i,\rho(i,j)/2-\bar{R}) \supseteq \mathcal{N}(i,r_n)$. Notice ${\bf E}[X_i]=0$; $\max_i \lvertX_i\rvert \precsim 1$ by (ref) and (ref); and $X_i \perp\!\!\!\perp X_j$ by (ref) since $\mathcal{N}(i,\rho(i,j)/2-\bar{R})$ and $\mathcal{N}(j,\rho(i,j)/2-\bar{R})) > 2\bar{R}$ are separated by a distance of at least $2\bar{R}$, which is an upper bound on the “diameter” of any cluster. Applying (ref) with $\rho(i,j)/2-\bar{R}$ in place of $r_n$, there exists $c_n' \precsim 1$ such that for any $i,j$,
These bounds yield
The right-hand side equals (ref) multiplied by $k_n$. Since $\eqref{covzz} \precsim \xi_n/n$ as shown in the proof of (ref), this is $\precsim k_n\xi_n/n$, as desired.
{\bf Step 3.} Letting $\tilde\sigma_n^2 = \text{Var}(\sqrt{k_n}n^{-1} \sum_{i\in\mathcal{N}_n} {\bf E}[Z_i \mid \mathcal{F}_i(r_n)])$,
by Minkowski's inequality and step 2. Since $\sigma_n^2 \succsim 1$ by assumption, it suffices to show
to establish (ref). We apply Theorem 3.6 of ross2011fundamentals, defining his $X_i$ as $n^{-1}k_n^{1/2} ({\bf E}[Z_i \mid \mathcal{F}_i(r_n)] - {\bf E}[Z_i])$ and his dependency graph $\bm{A}$ by connecting units $i,j$ in $\bm{A}$ if and only if $j \in \Lambda_i$. This is a dependency graph because $j \not\in\Lambda_i$ implies that the treatment assignments determining $\mathcal{F}_i(r_n)$ are independent of those determining $\mathcal{F}_j(r_n)$ under (ref). The maximum degree of $\bm{A}$ is at most $\max_i \lvert\Lambda_i\rvert \precsim n/k_n$ by (ref). Therefore, by (ref), the right-hand side of (3.8) in ross2011fundamentals is asymptotically bounded above by
so (ref) follows from his (3.8). This completes the proof of (ref)
{\bf Step 4.} Note that $\bar{\theta} = {\bf E}[\tilde\theta]$, where $\tilde\theta$ is the Horvitz-Thompson analog of $\hat\theta$ defined in (ref). The bias of $\tilde\theta$ is $\lvert\bar{\theta} - \theta^*\rvert \precsim (k_n\xi_n/n)^{\gamma/d}$ by (ref). Given that $k_n \prec (n/\xi_n)^{\frac{2\gamma}{2\gamma+d}}$,
in which case
so (ref) follows from (ref). \rule{0.08in}{0.08in}
In what follows, steps 1--4 concern case $\hat\sigma^2(1)$, and step 5 concerns $\hat\sigma^2(2)$. Define $Z_i$ as in (ref) and $\hat{Z}_i$ as in (ref). Let
This equals (ref), which is non-negative for any $n$.
By (ref), $\lvert\sigma_n - \tilde\sigma_n\rvert \prec 1$, where
by (ref). Thus, to show that $\hat\sigma^2(1) = \sigma_n^2 + \mathcal{B}_n + o_p(1)$, it suffices to prove
{\bf Step 1.} By definition, $\hat\sigma^2(1) = n^{-2}k_n \sum_{i\in\mathcal{N}_n} \sum_{j\in\mathcal{N}_n} \hat{Z}_i \hat{Z}_j A_{ij}(1)$. We prove that
In the formula for $\hat\sigma^2(1)$, replace $\hat{Z}_i$ with $\hat{Z}_i \pm Z_i$ to obtain
By (ref), $\max_i \sum_j A_{ij}(1) \precsim n/k_n$. By (ref) and (ref), $\max_i \lvertZ_i\rvert \precsim 1$, and
which is $\prec 1$ by the proof of (ref). Hence, $\eqref{2904yuws} \prec 1$.
{\bf Step 2.} We prove that
where
In the formula of $\check\sigma^2(1)$, replace $Z_i$ with $Z_i \pm {\bf E}[Z_i]$ to obtain
We need to show that the second term on the right is $\prec 1$. For $V_i = \sum_{j\in\mathcal{N}_n} {\bf E}[Z_j] A_{ij}(1)$,
By (ref) and (ref), $\max_i\lvert{\bf E}[Z_i]\rvert \precsim 1$. Since $\max_i \sum_{j\in\mathcal{N}_n} A_{ij}(1) \precsim n/k_n$ by (ref), $\max_i\lvertV_i\rvert \precsim n/k_n$. Therefore,
The argument in (ref)--(ref) can be applied to show that this is $\prec 1$. The one distinction is that the above expression has $Z_i$ in place of $\tilde Z_i$. These only differ because $Y_i$ in the expression of $Z_i$ is centered by $\mu_t$, unlike the expression of $\tilde Z_i$, but the centering is immaterial since it cancels out when deriving the analogs of the covariance bounds (ref) and (ref) using (ref).
{\bf Step 3.} We prove that $\mathcal{B}_n(1) = \mathcal{B}_n + o_p(1)$. Noting that $C_{c(i)} \subseteq \Lambda_i$, decompose
By (ref)(a), $\Lambda_i \subseteq \mathcal{N}(i,2(r_n+U_n)) \subseteq \mathcal{N}(i,3U_n)$, and $\Lambda_i\backslash C_{c(i)} \subseteq \mathcal{N}(i,3U_n)\backslash \mathcal{N}(i,L_n)$. By (ref)(b), there exists a non-negative sequence $\alpha_n \precsim 1$ such that $3U_n - L_n = \alpha_n (n/(k_n\xi_n))^{1/d}$. Then by (ref)(b), following (ref),
By (ref) and (ref), $\max_i \lvert{\bf E}[Z_i]\rvert \precsim 1$, so
which is $\prec 1$ since $k_n \prec n/\xi_n$ by assumption.
{\bf Step 4.} Abbreviate $\mathcal{F}_i \equiv \mathcal{F}_i(r_n)$. At the top of the proof, we noted that
By steps 1--3,
It therefore remains to show that the following is $\prec 1$:
As previously argued, $\max_i \sum_{j\in\mathcal{N}_n} {\bf E}[Z_j] A_{ij}(1) \precsim n/k_n$, and $\lvertn^{-1} \sum_{i\in\mathcal{N}_n} (Z_i-{\bf E}[Z_i])\rvert \prec 1$ by the proof of (ref), so $[P2] \prec 1$, while
Using (ref), (ref), (ref), and the assumption that $k_n \prec n/\xi_n$,
Observe that $[P1.3]$ has expectation zero and variance equal to
The covariance term is zero if $k,\ell \not\in \Lambda_i \cup \Lambda_j$, so by (ref) and (ref),
for $\Xi_{ij} = \{(k,\ell)\colon \ell\in\Lambda_k \text{ and } \{k,\ell\} \cap (\Lambda_i \cup \Lambda_j) \neq \emptyset\}$. Since $\max_i \lvert\Lambda_i\rvert \precsim n/k_n$ and $\ell\in\Lambda_k$ implies $k \in \Lambda_\ell$, we have $\max_{i,j\in\mathcal{N}_n} \lvert\Xi_{ij}\rvert \precsim (n/k_n)^2$. Therefore,
{\bf Step 5.} Recall the definition of $\tilde\sigma_n^2$ from (ref). We next prove that
Then the claimed result for $\hat\sigma^2(2)$ follows from steps 1, 2, and 4 by replacing, for all $i,j$, every occurrence of “$A_{ij}(1)$” and “$\Lambda_i$” with “$A_{ij}(2)$” and “$C_{c(i)}$”, respectively. Write
By (ref) and (ref), $\max_{i,j} \lvert\text{Cov}({\bf E}[Z_i \mid \mathcal{F}_i], {\bf E}[Z_j \mid \mathcal{F}_j])\rvert \precsim 1$, so using the argument in step 3,
\rule{0.08in}{0.08in}
We draw on the clever argument used in the proof of Theorem 3 of cao2024inference. To obtain a contradiction, suppose there is a positive sequence $\ell_n \prec 1$ such that, for infinitely many $n$, there are two clusters $C_1$ and $C_2$ with medoids $i_1,i_2$ such that $\rho(i_1,i_2) < \ell_n (n/(k_n\xi_n))^{1/d}$, and no other pair of medoids is closer in distance. There are two cases to consider along the subsequence of such $n$'s, and all asymptotic statements that follow are with respect to this subsequence.
{\bf Case 1.} $\min\{\lvertC_1\rvert, \lvertC_2\rvert\} \precsim n/k_n$. Then the argument largely proceeds as in the proof of Theorem 3 of cao2024inference but with additional adjustments to account for the divergence of $k_n,\xi_n$. Since clusters partition $\mathcal{N}_n$, for every $n$ there must exist some cluster $C_3$ of size at least $n/k_n$. Let $i_3$ denote its medoid and $R_3$ be the largest distance between $i_3$ and an element in $C_3$. By (ref)(a),
Let $i_3' \in C_3$ be such that $\rho(i_3',j) \geq 0.5 R_3$ for any medoid $j$. For instance, the unit $i_3^* \in C_3$ for which $\rho(i_3^*,i_3)=R_3$ is one such candidate. For any other medoid $j\neq i_3$, $\rho(i_3',j) \geq 0.5R_3$ because otherwise step 1 of (ref) would have assigned $i_3'$ to a closer medoid.
Consider a hypothetical update in step 2 of (ref) that replaces $i_2$ with $i_3'$. Then all units in $\mathcal{N}(i_3',R_3/8)$ are optimally reassigned to the cluster with medoid $i_3'$ because they are all by construction at least distance $3/8\cdot R_3$ away from any other medoid. This reassignment reduces the total cost by at least $\lvert\mathcal{N}(i_3',R_3/8)\rvert R_3/4 \succsim n/k_n \cdot (n/(k_n\xi_n))^{1/d}$ by (ref)(a) and (ref). On the other hand, in the worst case, all other units in $C_2$ are reassigned to medoid $i_1$, which increases the total cost by at most $\lvertC_2\rvert \ell_n (n/(k_n\xi_n))^{1/d} \prec n/k_n \cdot (n/(k_n\xi_n))^{1/d}$. Hence, the update is overall cost-reducing for $n$ sufficiently large, which contradicts for such $n$ the supposition that the clusters are the output of (ref), in particular violating step 2.
{\bf Case 2.} $\min\{\lvertC_1\rvert, \lvertC_2\rvert\} \succ n/k_n$. Without loss of generality, suppose $\lvertC_2\rvert \succ n/k_n$. Then the argument proceeds similarly to case 1 but using $C_2$ in place of $C_3$. In particular, let $R_2$ be the largest distance between $i_2$ and an element of $C_2$. By an argument similar to (ref), $R_2 \succ (n/(k_n\xi_n))^{1/d}$.
Let $i_2' \in C_2$ be such that $\rho(i_2',j) \geq 0.5 R_2$ for any medoid $j$. Consider a hypothetical update in step 2 of (ref) that replaces $i_2$ with $i_2'$. Then all units in $\mathcal{N}(i_2',R_2/8)$ are optimally reassigned to the cluster with medoid $i_2'$. This reassignment reduces the total cost by at least $\lvert\mathcal{N}(i_2',R_2/8)\rvert R_2/4 \succ \lvert\mathcal{N}(i_2',R_2/8)\rvert (n/(k_n\xi_n))^{1/d}$. On the other hand, in the worst case, all other units in $C_2$ are reassigned to medoid $i_1$, which increases the total cost by at most $\lvertC_2\rvert \ell_n (n/(k_n\xi_n))^{1/d}$. Since $\lvertC_2\rvert \leq \lvert\mathcal{N}(i_2,R_2)\rvert \sim \lvert\mathcal{N}(i_2',R_2/8)\rvert$ by (ref)(a), the update is overall cost-reducing for $n$ sufficiently large, which contradicts the supposition that the clusters are the output of (ref). \rule{0.08in}{0.08in}
{\bf Step 1.} We first prove that, for any sequence $\{C_n\}_{n\in\mathbb{N}}$ with $C_n \in \mathcal{C}_n$ for each $n$, we have $\lvertC_n\rvert \sim n/k_n$. Let $\ell_n$ be given as in (ref) and $i_n$ the medoid of cluster $C_n$. By (ref)(a), $\lvert\mathcal{N}(i_n,0.5 \ell_n (n/(k_n\xi_n))^{1/d})\rvert \succsim n/k_n$. All units in this neighborhood must be elements of $C_n$ by (ref) and step 1 of (ref), so $\lvertC_n\rvert \succsim n/k_n$. Since this must be true for all sequences of clusters, it cannot be the case that $\lvertC_n\rvert \succ n/k_n$, given that the total number of units is $n$. Hence $\lvertC_n\rvert \sim n/k_n$.
{\bf Step 2.} We prove the direction $R_n^* \precsim (n/(k_n\xi_n))^{1/d}$. Suppose to the contrary that $R_n^* \succ (n/(k_n\xi_n))^{1/d}$. Let $i_n' \in C_n$ be such that $\rho(i_n',j_n) \geq 0.5 R_n^*$ for any medoid $j_n$. For instance, the unit $i_n^* \in C_n$ for which $\rho(i_n^*,i_n)=R_n^*$ is one such candidate. For any other medoid $j_n\neq i_n$, $\rho(i_n',j_n) \geq 0.5R_n^*$ because otherwise it would be cost-reducing for (ref) to assign $i_n'$ to a closer medoid.
Consider a hypothetical update in step 2 (ref) that replaces $i_n$ with $i_n'$. Then all units in $\mathcal{N}(i_n',R_n^*/8)$ are optimally reassigned to the cluster with medoid $i_n'$. This reassignment reduces the total cost by at least $\lvert\mathcal{N}(i_n',R_n^*/8)\rvert R_n^*/4 \succ n/k_n (n/(k_n\xi_n))^{1/d}$ by (ref)(a). On the other hand, in the worst case, all other units in $C_n$ are reassigned to the nearest existing medoid.
Observe that the distance between $i_n$ and the nearest other medoid must be $\precsim (n/(k_n\xi_n))^{1/d}$. If instead that distance were $\delta_n \succ (n/(k_n\xi_n))^{1/d}$, then $C_n$ would contain $\mathcal{N}(i_n, 0.5\delta_n)$ by step 1 of (ref), in which case it would have size $\succ n/k_n$ by (ref)(a), which contradicts step 1.
Therefore, the worst-case increase in total cost from the medoid replacement is $\precsim \lvertC_n\rvert (n/(k_n\xi_n))^{1/d} \precsim n/k_n \cdot (n/(k_n\xi_n))^{1/d}$ by step 1. The update is overall cost-reducing for $n$ sufficiently large, which contradicts the supposition that the clusters are the output of (ref).
{\bf Step 3.} We prove the other direction. By (ref) and step 1 of (ref), $\mathcal{N}(i_n, \ell_n (n/(k_n\xi_n))^{1/d}/2) \subseteq C_n$. Hence, $R_n^* \geq \ell_n (n/(k_n\xi_n))^{1/d}/2$, so $R_n^* \succsim (n/(k_n\xi_n))^{1/d}$. \rule{0.08in}{0.08in}
Recall from (ref) the definition of cluster centroids, $U_n$, and $L_n$. For any $i\in\mathcal{N}_n$, all clusters intersecting $\mathcal{N}(i,r_n)$ must be subsets of $\mathcal{N}(i, r_n + 2U_n)$ by the assumption. Together with (ref), we have that there exists $\alpha_n \sim 1$ such that $\mathcal{N}(i, r_n + 2U_n) \subseteq \mathcal{N}(i,\alpha_nL_n)$ for any $n$ and $i\in\mathcal{N}_n$. Invoking (ref) once more, it follows that $\phi_i$ is at most the number of $L_n$-balls that fill $\mathcal{N}(i, \alpha_nL_n)$ without intersecting.
We seek to bound the $2L_n$-packing number wainwright2019high of $\mathcal{N}(i,\alpha_nL_n)$. By Lemma 5.5 of wainwright2019high, this is at most the $L_n$-covering number of the ball wainwright2019high, and we denote this number by $N_n(i)$. We bound this following the proof of Lemma 5.7(b) in wainwright2019high.
Construct a maximal $L_n/2$-packing of $\mathcal{N}(i,\alpha_n L_n)$ with cardinality $M_i$ and centroids $\{\theta_m\}_{m=1}^{M_i}$. This is also an $L_n$-covering of $\mathcal{N}(i,\alpha_n L_n)$, so $N_n(i) \leq M_i$. The balls of the packing $\{\mathcal{N}(\theta_m, L_n/2)\}_{m=1}^{M_i}$ are disjoint and contained in $\mathcal{N}(i,\alpha_n L_n + L_n/2)$, so
By (ref)(a), there exists $C>0$ independent of $i$ and $n$ such that
Hence, $M_i \leq C^2 (1+2\alpha_n)^d \precsim 1$, so $\max_i N_n(i) \precsim 1$, which proves the first claim of the lemma.
The second claim follows from the expression in (ref), the first claim, and (ref). \rule{0.08in}{0.08in}
Let $\bar{R} = \max_j R_j$, the latter defined prior to (ref). Observe that $\Lambda_i \subseteq \mathcal{N}(i,r_n + 2\bar{R} + r_n)$. By (ref) and (ref), $r_n \precsim \bar{R} \precsim (n/(k_n\xi_n))^{1/d}$, so by (ref)(a), $\max_i \lvert\Lambda_i\rvert \precsim n/k_n$. \rule{0.08in}{0.08in}
Since $\hat\kappa_t$ has mean one, it remains to compute the variance. By (ref),
By (ref), $\min_i p_{ti} \succsim 1$, and by Assumptions (ref)(a) and (ref), $\max_j \lvertC_j\rvert \precsim n/k_n$, so the right-hand side is at most of order
\rule{0.08in}{0.08in}
By (ref), $p_t^+ \in (0,1)$. Since
it is enough to show that $\hat{p}_t^+ - p_t^+ \precsim k_n^{-1/2}$. Clearly ${\bf E}[\hat{p}_t^+] = p_t^+$. By (ref),
By Assumptions (ref)(a) and (ref), $\max_j \lvertC_j\rvert \precsim n/k_n$, so the right-hand side of the above display is $\precsim n^{-2} k_n (n/k_n)^2 = k_n^{-1}$. \rule{0.08in}{0.08in}
\FloatBarrier \phantomsection \addcontentsline{toc}{section}{References}