EconBase
← Back to paper

Causal clustering: design of cluster experiments under network interference

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.

79,239 characters · 16 sections · 89 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Causal clustering: design of cluster experiments under network interference

\if00 \fi

\if10 { \ \\

\ \\

} \fi

abstractThis paper studies the design of cluster experiments to estimate the global treatment effect in the presence of network spillovers. We provide a framework to choose the clustering that minimizes the worst-case mean-squared error of the estimated global effect. We show that optimal clustering solves a novel penalized min-cut optimization problem computed via off-the-shelf semi-definite programming algorithms. Our analysis also characterizes simple conditions to choose between any two cluster designs, including choosing between a cluster or individual-level randomization. We illustrate the method's properties using unique network data from the universe of Facebook's users and existing data from a field experiment.

{\it Keywords:} Experimental Design, Spillover Effects, Causal Inference, Cluster Designs. \\ {\it JEL Codes:} C10, C14, C31, C54

\spacingset{1.5}

Introduction

Consider a (large) population of $n$ individuals connected under a single observed network. Researchers are interested in conducting an experiment to estimate the global average treatment effect, i.e., the difference between the average effect of treating all versus none of the individuals in the population. Treating an individual may generate spillovers to her friends in the network. To capture such effects, researchers conduct a cluster experiment. Individuals are first partitioned into clusters. Within a cluster, either all units are assigned to the treatment or all units are assigned to the control group. Finally, researchers estimate treatment effects by taking a difference between the average outcomes of treated and control units (possibly adjusting for baseline covariates). The cluster design does not require modeling the dependence of individual outcomes on neighbors' assignments, but it requires a choice of clusters and some assumptions on the extent of the spillovers along the network. For example, cluster experiments on online platforms require choosing a partition of the social network, and field experiments require choosing the unit of randomization, such as villages or regions. This raises the question of how many and which clusters to use in experiments.

Typical approaches in economic research assume prior knowledge of many independent clusters. There are many settings where this information is not available, and instead, units in the population have different degrees of connections.\footnote{For example, when using villages as clusters, individuals may interact in the same and nearby villages. See egger2022general for an empirical example.} This paper provides an econometric framework to choose when and how to design the clusters in cluster experiments. Different from existing clustering algorithms geared towards community detection, we motivate the choice of clusters based on the task of estimating global treatment effects. The choice of clustering must balance two competing objectives: the larger the clusters (and the smaller the number of clusters), the smaller the bias of the estimated global effect, but the larger its variance. We introduce an algorithmic procedure -- entitled Causal Clustering -- to choose the clustering that minimizes a weighted combination of the worst-case bias and variance as a function of the network and clusters. The worst-case approach encodes uncertainty over the dependence of individual outcomes on neighbors' assignments. We study the optimal clustering (and how many clusters to use) in this scenario.

We focus on a class of models where spillover effects are local (within neighbors) and small relative to the outcomes' variance but possibly non-negligible for inference. This framework encodes the researchers' uncertainty on the magnitude of spillover effects by modeling first-order neighbors' effect (and possibly also direct and overall effects) local to zero possibly at arbitrary rates.\footnote{The assumption of spillovers within first-order neighbors can be relaxed here by assuming that higher-order spillovers are an order of magnitude smaller than first-order spillovers.} From a practical perspective, we therefore focus on a class of designs which are (approximately) optimal when spillover (and direct) effects are small, common, for instance, in online experiments karrer2021network. We characterize optimal clustering as a function of the expected magnitude of the largest spillover effects that the experiment can generate. The largest size of the spillover effects is a key input in our algorithms, whose characterization, in practice, is necessary for the design of the experiment but can be challenging. This parameter can be informed by previous experiments, in the same spirit of minimum detectable effects used in power analysis baird2018optimal, or using some particular modeling assumptions. We provide a practical discussion in Section (ref).

Our analysis proceeds as follows. First, we provide a formal characterization of the worst-case bias and variance. The worst-case bias is closely related to a particular notion of between-clusters connectedness, defined as the per-individual average number of friends in other clusters, similar for instance to notions of cluster purity brennan2022cluster, toulis2013estimation. On the other hand, the worst-case variance can potentially be an arbitrary function of within-clusters covariances and between-clusters covariances: individuals in the same cluster have identical assignments and in different clusters may share common neighbors. We show that the variance only depends on the average squared cluster size, up to a small approximation error. This result formalizes the intuition that a larger number of clusters, with a small variance in cluster size, decreases the variance of the estimator. Taken together, these two results provide, as we show, a simple and intuitive expression of the worst-case mean squared error associated with a given clustering.

We draw the implications of these results for choosing between a cluster experiment (for a given clustering) or assigning treatments independently between individuals (i.e., Bernoulli design or completely randomized design). Suppose the magnitude of the spillover effects is smaller than the square root of the number of clusters. In that case, the variance component dominates the bias, and a Bernoulli design is preferred (where a Bernoulli design is a special case of a cluster design with clusters containing a single unit). Vice-versa, a cluster design is preferred if the bias dominates the variance. We provide practitioners with a simple and explicit decision rule between cluster and Bernoulli designs that can be easily implemented.

We then turn to the design of optimal clustering, i.e., which clustering and how many clusters to use. The choice of optimal clustering reduces to a novel sequence of penalized minimum cut optimization problems, with a penalty that depends both on a measure of cluster purity and (different from standard clustering algorithm) on the variance contribution of each cluster. Each problem is solved via off-the-shelf semidefinite programming algorithms.

We provide an empirical application using unique network data from the universe of Facebook users. We show that our procedure provides an explicit ranking between different clustering algorithms implemented at scale at Facebook. We also illustrate trade-offs in the choice of the clustering algorithm and properties of the graph. A second application using network data from cai2015social illustrates the method's advantages for choosing clusters in field experiments. We show that the choice of the clusters based on village identity is sub-optimal in this application because the number of village clusters is too small relative to the optimal clustering because the number of village clusters is too small relative to the optimal clustering, and provide an explicit choice of the clustering, and provide an explicit choice of the clustering. These results illustrate the applicability of the method in choosing the clustering and number of clusters.

This paper connects to the literature on experimental design, causal inference with spillover effects, and clustering. Existing methods for experimental designs with spillover effects include cluster and saturation designs, studied in taylor2018randomized, baird2018optimal, basse2016analyzing, pouget2018dealing, viviano2020policy, karrer2021network, cai2022optimizing, among others. These papers provide an analysis of the properties of particular designs for a given clustering or clustering algorithm. Different from the current paper, these references either do not study the question of the optimal clustering algorithm for experimental design or only provide heuristic comparisons of different clustering methods.

A particular class of clustering algorithms is $\varepsilon$-net clustering algorithms -- which consist of assigning sequentially individuals in the same neighborhood to the same cluster ugander2013graph, eckles2017design. Variants of these algorithms have been recently studied in leung2022rate leung2023design for spatial spillovers, and faridani2022rate in general non-Euclidean spaces. These papers provide an optimal rate of the clusters' size as a function of how interference decays in space for asymptotically unbiased estimators. Here, instead, we provide an explicit trade-off and novel characterization of the bias and variance, leveraging local asymptotics.

In particular, as leung2022rate conjectures, network dependence is difficult to handle for general cluster experiments. A key challenge is that the estimator's variance depends not only on the network structure but also on the clustering structure: outcomes in two different clusters are dependent if they share common neighbors or if they are in the cluster with some cross-connection between the two clusters. The lack of explicit and simple-to-optimize expressions for the variance have motivated designs that either focus on minimizing the bias only brennan2022cluster, taylor2018randomized, or on rate-optimal designs only leung2022rate, faridani2022rate. Our local asymptotic framework allows us to explicitly characterize the variance through a simple and intuitive expression. This result allows us to characterize the optimal clustering as the solution of a novel and simple trace-optimization program (different from the algorithms studied in the previous references).

Additional references of experiments with and without networks in the absence of cluster experiments are basse2018model, and jagadeesan2020designs for estimating direct instead of global average treatment effects studied here; kang2016peer who study encouragement designs, without focusing on the problem of mean-squared error optimal designs; viviano2020experimental for experiments on networks using information from a pilot study; basse2018limitations discuss limitations of design-based causal inference under interference. Recent econometric literature studying experimental design focusing however on $i.i.d.$ data include cai2024performance, tabord2023stratification, armstrong2022asymptotic among others. None of these papers study (optimal) cluster designs.

The literature on treatment effects under network interference includes aronow2017estimating, hudgens2008toward, toulis2013estimation, manski2013identification, leung2020treatment, athey2018exact, goldsmith2013social, savje2017average, ogburn2017causal, manresa2013estimating, li2022random, often focusing on first-order degree interference.\footnote{Our local asymptotic framework extends to settings with higher order dependence, assuming that higher order neighbors generate spillovers of smaller order compared to first order neighbors, and the network is sufficiently sparse.} None of the above references study experimental designs. Finally, we relate more broadly to the literature on clustering in economics -- see wooldridge2003cluster, abadie2017should and references therein -- with the difference that here we focus on optimal clustering, rather than taking the clusters as given -- and the literature on graph clustering, including von2007tutorial, newman2013community,newman2013spectral, lei2019unified, lei2020consistency, li2022hierarchical, and references therein. Different from this last set of papers on graph clustering, this paper focuses on treatment effect estimation instead of community detection.

The remainder of the paper is organized as follows: Section (ref) presents the setup; Section (ref) characterizes the bias and variance and contrast cluster and Bernoulli designs; Section (ref) formalizes the optimization problem for the optimal clustering; Section (ref) presents two applications and numerical studies, Section (ref) presents several extensions, and Section (ref) contains recommendations for practice.

Setup

We consider a setting with $i \in \{1, \cdots, n\}$ units. Let $Y_i \in \mathbb{R}$ denote the observed outcome of interest for individual $i$, $D_i \in \{0,1\}$ denote a binary treatment assignment, and $\mathbf{D} \in \{0,1\}^n$ the vector of treatment assignments of each unit. Define $\mathbf{A}$ a symmetric adjacency matrix with $\mathbf{A}_{i,j} \in \{0,1\}$, and $ Y_i(\mathbf{d}), \mathbf{d} \in \{0,1\}^n $ the potential outcome as a function of the entire vector of treatment assignments, with $Y_i = Y_i(\mathbf{D})$. We implicitly condition on $\mathbf{A}$, observed to the researchers, and (unobserved) potential outcomes $Y(\mathbf{d})$, unless otherwise specified (see Remark (ref) for an extension with unobserved $\mathbf{A}$). We define $\mathcal{N}_i = \{j: \mathbf{A}_{i,j} = 1\}$ the set of individuals connected to $i$. To avoid dividing by zero, we define $|\mathcal{N}_i| = \max\{1, \sum_j \mathbf{A}_{i,j}\}$. Let $\mathcal{N}_{n, \max} = \max_{i \in \{1, \cdots, n\}} |\mathcal{N}_i|$, the maximum degree. We focus on estimating the global average effect,

equation[equation omitted — 151 chars of source]

The overall effect defines the effect if all units had received the treatment compared to none of the individuals receiving the treatment. This is a target estimand in many applications when researchers are interested in implementing the policy at scale, e.g., egger2022general, muralidharan2017experimentation, karrer2021network (see Appendix (ref) for a discussion on other estimands). Table (ref) summarizes the notation.

Potential outcomes and spillover effects

Next, we impose restrictions on the spillover effects.

ass[First-order local interference] For $i \in \{1, \cdots,n\}$, $ Y_i(\mathbf{d}) = \mu_i(\mathbf{d}_i, \mathbf{d}_{\mathcal{N}_i}), \forall \mathbf{d} \in \{0,1\}^n, $ for some functions $\mu_i(1, \cdot) \in \mathcal{M}_{1,i}, \mu_i(0, \cdot) \in \mathcal{M}_{0, i}$ for some set of functions $\mathcal{M}_{1, i}, \mathcal{M}_{0, i}$.

Assumption (ref) states that spillovers occur between neighbors and allows for arbitrary dependence of potential outcomes with neighbors' assignments. One-degree neighborhood is consistent with models often used in applications, e.g., cai2015social, sinclair2012detecting, muralidharan2017experimentation. athey2018exact provide a framework for testing Assumption (ref). Higher-order interference can be accommodated, although higher-order effects are often small and difficult to detect (see Remark (ref)). With a slight abuse of notation, we occassionally write $\mu_i(\mathbf{d})$ for $\mu_i(\mathbf{d}_i, \mathbf{d}_{\mathcal{N}_i})$ for notational convenience.

We will refer to $\tau_{n, \mu}$ as the overall effect in Equation (ref) to make the dependence of $\tau_n$ on $(\mu_1, \cdots, \mu_n)$ explicit. We do not assume that we know or can estimate consistently $\mu_i$. (The functions $\mu_i$ and their classes $\mathcal{M}_{0, i}, \mathcal{M}_{1, i}$ are indexed by $i$, and therefore they cannot be consistently estimated.) Instead, we allow for arbitrary classes $\mathcal{M}_{1, i}, \mathcal{M}_{0, i}$ of potential outcome functions, as long as such classes satisfy the conditions below. These classes have to be sufficiently large to accommodate rich structures in the data, as well as satisfy some restrictions to be able to estimate features of causal effects.

ass[Class of potential outcome functions] The potential outcomes $\mu = (\mu_1, \cdots, \mu_n) \in \mathcal{M} = \otimes_{i=1}^n \mathcal{M}_i$ where $\mathcal{M}_i$ is the set such that $\mu_i(1, \cdot) \in \mathcal{M}_{1, i}, \mu_i(0, \cdot) \in \mathcal{M}_{0, i}$. For all $i \in \{1, \cdots, n\}$, $d \in \{0,1\}$, the function class $\mathcal{M}_{d,i}$ is the class of functions that include all functions $\mu_i(d, \cdot)$ satisfying the following two conditions: \begin{itemize} • $\mu_i(d, \cdot)\in [\psi_L, \psi_R]$ for some constants $-\infty < \psi_L < \psi_R < \infty$ that do not vary with $n$; • for all $\mathbf{d} \in \{0,1\}^{|\mathcal{N}_i|}$, for some $\bar{\phi}_n \le \psi_R - \psi_L$ \begin{equation} \begin{aligned} \Big|\mu_i(0, \mathbf{d}) - \mu_i(0, \mathbf{0})\Big| & \le \frac{\bar{\phi}_n}{|\mathcal{N}_i|} \sum_{k \in \mathcal{N}_i} \mathbf{d}_k, \\ \Big|\mu_i(1, \mathbf{d}) - \mu_i(1, \mathbf{1})\Big| & \le \frac{\bar{\phi}_n}{|\mathcal{N}_i|} \sum_{k \in \mathcal{N}_i} \Big(1 - \mathbf{d}_k\Big). \end{aligned} \end{equation} \end{itemize}

Assumption (ref) imposes two restrictions. Condition (i) states that potential outcomes are uniformly bounded.\footnote{This restriction is common in the literature, e.g., kitagawa2020should, and can be relaxed by assuming random sub-gaussian potential outcomes. } Condition (ii) is our main restriction on the exposure mapping. Condition (ii) is attained if $\mathcal{M}_{1, i}, \mathcal{M}_{0, i}$ are Lipschitz function classes in the share of treated neighbors. Condition (ii) states that potential outcomes vary in the share of neighbors' treatments by at most $\bar{\phi}_n$. Here $\bar{\phi}_n$ captures the magnitude of (largest) spillovers. The condition that $\bar{\phi}_n \le \psi_R - \psi_L$ is a regularity assumption that guarantees that $(ii)$ is consistent with the boundedness restriction in $(i)$. In Appendix (ref), we relax Assumption (ref), by imposing $(ii)$ only on the expected value of the potential outcomes (instead of the realized potential outcomes). In Appendix (ref), we relax Assumption (ref) by allowing for heterogeneity in the range of potential outcomes and the spillover effects.\footnote{Also, note that because $\mathcal{M}_{d,i}$ contains all functions satisfying $(i)$ and $(ii)$ it also contains constant functions as well as linear functions in the share of treated neighbors with linear coefficient bounded by $\bar{\phi}_n$.}

The component $\bar{\phi}_n$ quantifies the largest magnitude of spillover effects and will play an important role in our asymptotic analysis as $n \rightarrow \infty$. We quantify the magnitude of spillover effects through the following asymptotic framework.

ass[Spillover local asymptotics] We consider a sequence of $\Big(\mathbf{A}, (Y_{i}(\cdot))_{i=1}^n\Big)$, indexed by $n$, where $\lim \sup_n \bar{\phi}_n \min\{n, \mathcal{N}_{n, \max}^2\} \le \varepsilon$, for some finite constant $\varepsilon \ge 0$.

Assumption (ref) formalizes the local asymptotic framework considered here, where $\bar{\phi}_n \mathcal{N}_{n, \max}^2 $ denotes the product of the largest spillover effects with the squared maximum degree.\footnote{Assumption (ref) can be sharpened by imposing restrictions on the average second-order degree instead of maximum degree; see Equation (ref) in the proof of Lemma (ref).} We assume that such a product is bounded in the limit by $\varepsilon \ge 0$. For $\varepsilon = 0$, spillover effects converge to zero asymptotically. This may occur if $\bar{\phi}_n \rightarrow 0$ at an arbitrary rate for networks with bounded degree de2018identifying. In this case, the rate of convergence of $\bar{\phi}_n$ characterizes the magnitude of the spillover effects. As in a local asymptotic framework hirano2009asymptotics, rates of converges of order equal or slower than the effective number of observations (allowed here) indicates non-negligible spillover effects even with $\varepsilon = o(1)$. Although we impose no restrictions on $\varepsilon$, we will think of $\varepsilon$ as small. Specifically, we provide a framework that characterizes approximately optimal experimental design, where the quality of the approximation depends on $\varepsilon$. Also, note that overall effects may also exhibit the same (or similar) asymptotic behavior, see Remark (ref).

Scenarios where $\bar{\phi}_n$ is local to zero formalize the idea that spillover effects (and possibly but not necessarily also overall treatment effects $\tau_n$) are small but possibly non-negligible for inference. For example, in contexts with online advertising lewis2015unfavorable observe that a treatment effect of 0.005 of a standard deviation of the outcome would be highly profitable for the company given the cost of advertising athey2023semi. Similar examples in digital advertisements with small effects are athey2023digital with $i.i.d.$ data and karrer2021network with network (Facebook) data, where total effect size (sum of direct and spillover effects) are also of the order of 0.005 standard deviations. Motivated by these applications we therefore focus on settings where first-order effects are small but possibly non-negligible. As we show below, this framework allows us to provide explicit characterization of the bias and variance of the clustering, different from models of decaying interference as for example in leung2022causal, where such characterization is not possible. In other scenarios where either the network is dense, or spillover effects have large magnitudes, the quality of the approximations we provide below will deteriorate.

exmp[Linear exogenous peer effects] Consider a class of functions of the form $\mu_i(\mathbf{d}) = \mu(T_i(\mathbf{d})) + \nu_i$, with $T_i(\mathbf{d}) = \Big[\mathbf{d}_i, (1 - \mathbf{d}_i) \times \frac{\sum_{j \neq i} \mathbf{A}_{i,j} \mathbf{d}_j}{\sum_{j \neq i} \mathbf{A}_{i,j}}, \mathbf{d}_i \times \frac{\sum_{j \neq i} \mathbf{A}_{i,j} \mathbf{d}_j}{\sum_{j \neq i} \mathbf{A}_{i,j}}\Big]$, and $\mu(t) = t^\top \beta$ for $\beta \in [-\bar{\Delta}_n,\bar{\Delta}_n] \times [-\bar{\phi}_n,\bar{\phi}_n]^2,$ for some arbitrary $\bar{\Delta}_n, \bar{\phi}_n$, and $\nu_i$ that is not a function of $\mathbf{d}$. Then Assumption (ref) holds. \qed
rem[Direct and overall effect] Our asymptotic framework also allows direct treatment effects (and global effects) to be of the same magnitude of spillover effects. Specifically, it is possible that $\tau_n$ is of the same order of $\bar{\phi}_n$. Therefore, our local asymptotic assumption does not require that spillover effects are local to the global effect (since $\tau_n$ can also converge to zero). Instead, our local asymptotics formalizes settings where the noise-to-signal ratio can be small (common in experiments on online platforms, karrer2021network). \qed
rem[Higher order interference and endogenous peer effects] Our setting generalizes to higher-order interference in two scenarios. First, suppose that friends up to degree $d < \infty$ generate spillovers in magnitude similar to first-degree friends. Our results extend after we define the set of friends as the set of friends up to degree $d$, and $\bar{\phi}_n$ the largest effect that such friends generate. Sparsity restrictions on the largest degree in Assumption (ref) are with respect to the number of friends up to degree $d$. In the second scenario, suppose that the assumption of first-order effects approximates higher-order effects up to a term of smaller order than first-order effects. In this case our results can extend as we discuss in Appendix (ref). In practice, higher-order effects can be small, leading to underpowered studies, especially when individuals have many friends and first-order effects capture most of the spillovers. \qed

Experimental design and estimation

Next, we turn to the class of designs and estimators considered here. Define a clustering of size $K_n$ as a set of sets of indicators satisfying $$

aligned\mathcal{C}_{ n} = \left\{c_k \subseteq \{1, \cdots, n\}, k \in \{1, \cdots, K_n\}, \bigcup_k c_k = \{1, \cdots, n\}, c_k \bigcap c_{k'} = \emptyset for k \neq k' \right\}.

$$ Here, $\mathcal{C}_{ n}$ denotes a particular partition of the units in the population with $K_n$ exclusive clusters. With an abuse of notation, let $c(i) \subseteq \{1, \cdots, n\}$ denote the cluster assigned to unit $i$, and $|c_k| = n_k$ the number of individuals in cluster $k$. For ease of exposition, we focus our discussion and assumptions below in the presence of a given clustering $\mathcal{C}_{ n}$ and return to choosing the optimal clustering in Section (ref).

ass[Cluster designs] $\mathcal{C}_{ n}$ is a deterministic function of $\mathbf{A}$ and is such that the number of units in cluster $k$ is $n_k = \gamma_k \frac{n}{K_n}$ with $\max_k \gamma_k \le \bar{\gamma} < \infty$, for some $\bar{\gamma} < \infty$. Assume that $D_i = \tilde{D}_{c(i)}$ almost surely, where for $k \in \{1, \cdots, K_n\}$, $\tilde{D}_{k}, | \mathbf{A}, \{Y_i(\mathbf{d})\}_{i \in \{1, \cdots, n\} \mathbf{d} \in \{0,1\}^N} \sim \mathrm{Bern}(0.5)$ and $\tilde{D}_1, \cdots, \tilde{D}_{K_n} | \mathbf{A}, \{Y_i(\mathbf{d})\}_{i \in \{1, \cdots, n\} \mathbf{d} \in \{0,1\}^N}$ are independent across clusters.

Assumption (ref) states the clustering is constructed using information from the adjacency matrix only (i.e., it is independent of potential outcomes), clusters are proportional in size\footnote{This restriction is sufficient but not necessary for our analysis, and it is imposed for expositional convenience. It is possible to relax such an assumption by assuming that $\sum_{k = 1}^{K_n} \gamma_k^2/K_n = \mathcal{O}(1)$.}, and individuals in a given cluster are all assigned either treatment or control with equal probability. Assumption (ref) restricts the class of designs to cluster designs, motivated by our focus on overall treatment effects and empirical practice.

Motivated by standard practice both in industrial applications and in field experiments with clusters, we consider estimators obtained by simple differences in means between treated and control clusters karrer2021network; see Remark (ref) for comparison with other estimators studied in the literature. Because $P(D_i = 1) = 1/2$, we construct a (biased) estimator of treatment effects as\footnote{Here, for expositional convenience, we make explicit the dependence of $\hat{\tau}_n$ with $\mathcal{C}_n$, where the clustering $\mathcal{C}_n$ only characterizes the distribution of treatment assignments $D_i$ but not the structure of the estimator.}

equation[equation omitted — 212 chars of source]

The estimator $\hat{\tau}_n(\mathcal{C}_{ n})$ is a simple difference in means between treated and control units that normalizes by the probability of treatment. Therefore, $\hat{\tau}_n(\mathcal{C}_{ n})$ depends on the clustering $\mathcal{C}_{ n}$ only because the distribution of the treatments depends on the clusters under Assumption (ref). Studying the estimator in Equation (ref) is a natural starting point for the analysis of cluster experiments. Variants of difference in means estimators (possibly also with regression adjustments discussed in Remark (ref)) are often used or studied in practice karrer2021network, holtz2020reducing, savje2017average. One could also normalize each sum in Equation (ref) by the number of treated and control units (instead of using knowledge about the treatment probability). This would improve the stability of the estimator but complicate the analysis when the number of treated units is stochastic (e.g., when clusters have different sizes).

rem[Alternative estimators] Alternative estimators studied in the literature are inverse probability weights estimators aronow2017estimating, ugander2013graph. These estimators reweight by the probability that the treatment of all friends being either zero or one. Unless researchers impose additional restrictions on the exposure mapping, in the network context, these estimators are highly sensitive to instability of the propensity score because the propensity score scale to zero exponentially in the maximum degree.\footnote{Even if the maximum degree grows at a slow rate, such as $\mathcal{N}_{n, \max} \propto \log(n), \bar{\phi}_n \propto 1/\log^2(n)$, so that effects are local to zero but non negligible for inference the propensity score converges to zero at rate $1/n$.} An alternative are model-based estimators, which are subject to model misspecification. A third class of estimators are trimming estimators, which trim observations based on their position in the network and also optimize over the choice of the estimators and designs. We study this class of estimators and designs in Appendix (ref). \qed
rem[Covariate adjustment] Our framework generalizes to settings that use covariate adjustment. Denote $\bar{\mu}_i$ an arbitrary predictor for $\mu_i(\mathbf{0})$ that only uses information from some arbitrary baseline observable characteristics (i.e., does not depend on the treatments or end-line outcomes in the experiment). The estimator with such an adjustment takes the form \begin{equation} \begin{aligned} \frac{2}{n} \sum_{i=1}^n (Y_i - \bar{\mu}_i)(2 D_i - 1). \end{aligned} \end{equation} Note that $Y_1 - \bar{\mu}_1, \ldots, Y_n - \bar{\mu}_n$ typically have different ranges ex ante even under Assumption (ref). To account for this, we can apply the more general result in Appendix (ref). \qed
rem[Saturation designs] It is possible to consider designs where treatments are assigned with different cluster-level probabilities. Our main insights for the bias and variance characterization follow similarly to what we present for cluster experiments, with appropriate modifications presented in Appendix (ref). \qed

Bias, Variance and MSE of a given clustering

As a necessary first step, this section focuses on studying the bias, variance, and mean-squared error of a given clustering $\mathcal{C}_{ n}$, deferring the question of choosing an optimal clustering to the next section. Denote $\mathbb{E}_\mu[\cdot]$ the expectation for given potential outcome functions $\mu$ (conditional on $\mathbf{A}$). The components $\mathcal{O}(\cdot), o(\cdot)$ in the following lemmas and theorems hold uniformly over all cluster designs with $K_n$ many clusters that satisfy Assumption (ref).

Our goal is to design experiments that minimize a weighted combination of the variance and squared bias, worst-case over the class of potential outcomes.

Our objective reads as follows

equation[equation omitted — 362 chars of source]

where the supremum is intended over $(\mu_i)_{i=1}^n$, and $\mathcal{M}$ is the product space of potential outcome functions as in Assumption (ref). The user may specify $\lambda$.

For $\lambda = 1$, $\mathcal{B}_n(\mathcal{C}_n, \lambda)$ coincides in the limit with the worst-case mean-squared error. It is also possible to study optimal clustering to minimize the length of confidence intervals for the overall effects, which we defer to Section (ref).

Worst-case bias, variance, and MSE

As a preliminary step in the characterization of the objective function, we first characterize the worst-case bias and variance. Let

equation[equation omitted — 198 chars of source]

and $\Big\{j: c(i) \neq c(j)\Big\}$ denotes the set of units $j$ in a different cluster from unit $i$.

lem[Worst-case bias] Let Assumptions (ref), (ref), (ref) hold. Then $ \sup_{\mu \in \mathcal{M} } \Big|\tau_{n, \mu} - \mathbb{E}_\mu[\hat{\tau}_n(\mathcal{C}_{ n})]\Big| = \bar{\phi}_n b_n(\mathcal{C}_n). $
proofSee Appendix (ref).

Lemma (ref) shows that the worst-case bias can be expressed as the number of friends of $i$ in a different cluster from $i$, appropriately reweighted by the overall number of friends of $i$. The size of the worst-case bias also depends on the magnitude of spillover effects $\bar{\phi}_n$. Lemma (ref) provides a formal justification to notions of between clusters connectedness based on the worst-case bias of the treatment effect estimator. In the proof, the worst-case bias is attained as the condition in Equation (ref) holds with equality (Appendix (ref) studies settings where effects are heterogeneous across individuals). The expression for the bias connects to measures of cluster purity such as those in goel2009respondent, toulis2013estimation, motivating clustering purity based on the class of models in Assumption (ref).

To derive the worst-case variance, we define

equation[equation omitted — 119 chars of source]

By Assumption (ref) (i), $\bar{\psi}$ can also be expressed as $\max_{\mu_i\in \mathcal{M}_i}(\mu_i(\mathbf{1}) + \mu_i(\mathbf{0}))^2$ for every $i$.

lem[Worst-case variance] Let Assumptions (ref), (ref), (ref), (ref) hold. Then as $K_n \rightarrow \infty$ $$ \small \begin{aligned} & \sup_{\mu \in \mathcal{M}} \mathbb{E}_\mu\Big[\Big(\hat{\tau}_n(\mathcal{C}_{ n}) - \mathbb{E}[\hat{\tau}_n(\mathcal{C}_{ n})]\Big)^2\Big] = \frac{1}{n^2} \sum_{k=1}^{K_n} n_k^2 \times \Big[\bar{\psi} + \mathcal{O}(\varepsilon) \Big]. \end{aligned} $$
proofSee Appendix (ref).

Lemma (ref) characterizes the worst-case variance as the number of clusters $K_n$ grows (at an arbitrary rate). Such bounds depend on the variance in cluster sizes, $ \frac{1}{n^2} \sum_{k=1}^{K_n} n_k^2, $ times a positive constant $\bar{\psi}$. By leveraging the local asymptotic framework in Assumption (ref), a key insight is that the variance is (asymptotically) driven by the within-cluster correlations instead of cross-cluster connections, up-to a $\varepsilon$-approximation. To our knowledge, this characterization of the worst-case variance is novel to the literature.

thmSuppose that Assumptions (ref), (ref), (ref), (ref) hold. Then as $K_n \rightarrow \infty$ $ \mathcal{B}_n(\mathcal{C}_{ n}, \lambda) = \mathcal{B}_n^*(\mathcal{C}_{ n}, \lambda)(1 + \mathcal{O}(\varepsilon)), $ where \begin{equation} \mathcal{B}_n^*(\mathcal{C}_{ n}, \lambda) = \bar{\psi}\cdot \frac{1}{n^2} \sum_{k=1}^{K_n} n_k^2 + (\lambda \bar{\phi}_n^2) \cdot b_n^2(\mathcal{C}_n). \end{equation}
proofSee Appendix (ref).
rem[Choice and sharper bound for $\bar{\psi}$] Because here we consider a finite population with fixed potential outcomes, the worst-case variance depends on $\bar{\psi}$, a constant proportional to the largest value of the realized outcomes. In Appendix (ref), we provide similar calculations in settings where potential outcomes are random functions, in which case we can interpret $\bar{\psi}$ as a sharper bound on the moments of potential outcomes. \qed
rem[Unobserved $\mathbf{A}$] Suppose that $\mathbf{A}$ is unobserved or partially observed, and researchers have a prior over $\mathbf{A}$. In this case, the worst-case variance does not change, the worst-case bias is replaced by $\bar{\phi}_n\mathbb{E}_{\mathbf{A}}[b_n(\mathcal{C}_n)]$, and $b_n^2(\mathcal{C}_n)$ in $\mathcal{B}_{n}^{*}(\mathcal{C}_n, \lambda)$ is replaced by $\mathbb{E}_{\mathbf{A}}[b_n^2(\mathcal{C}_n)]$, where the expectation may depend on partial network information breza2020using. In this case, the sparsity restrictions in Assumption (ref) can be replaced by assumptions on the expected squared degree over the prior -- see Equation (ref) in the proof of Lemma (ref). \qed

Comparison between a given clustering and a Bernoulli design

Before discussing the optimal choice of the cluster, we explore a simpler question, that is central in many randomized trials with spillovers. Suppose we are given a clustering $\mathcal{C}_n$, but we have the option to run a Bernoulli design where treatments are assigned independently between individuals. When should be prefer the (more complex) cluster design over a Bernoulli design? Denote a Bernoulli design as $\mathcal{C}_{\mathrm{B},n} = \{c_k = \{k\}, k \in \{1, \cdots, n\}\}$.

thm[Rule of thumb] Suppose that Assumptions (ref), (ref), (ref) hold. Let $\underline{\gamma} = \frac{1}{K_n} \sum_{k=1}^{K_n} \gamma_k^2$, where $\gamma_k = n_k K_n/n$. Then $$ \small \begin{aligned} \frac{\mathcal{B}_n^{*}(\mathcal{C}_{\mathrm{B},n}, \lambda)}{\mathcal{B}_n^{*}(\mathcal{C}_{ n}, \lambda)} \ge 1 \quad \text{ if } \quad (\lambda \bar{\phi}_n^2/\bar{\psi})^{-1} \le \frac{K_n(1 - b_n(\mathcal{C}_n)^2) }{\underline{\gamma} }. \end{aligned} $$
proof[Proof of Theorem (ref)] See Appendix (ref).

Theorem (ref) provides an explicit rule for choosing a cluster design with given clusters $\mathcal{C}_{ n}$ over a Bernoulli design $\mathcal{C}_{\mathrm{B},n}$, using the approximate objective $\mathcal{B}_n^*(\mathcal{C}_n, \lambda)$. The right-hand side depends on three observables: (i) $\underline{\gamma}$, (ii) the expected bias of the clustering method (as a function of $\mathbf{A}$), and (iii) the number of clusters $K_n$. The left-hand side depends on $\bar{\phi}_n^2/\bar{\psi}$. For $\lambda = 1$, known $\bar{\psi}$, the rule of thumb provides the smallest spillover effects that would guarantee that the cluster design dominates the Bernoulli design. Figure (ref) illustrates the exact threshold for the spillovers that justifies a cluster experiment. In our empirical application in Section (ref) we find $\bar{\phi}_n \ge 0.003$ would motivate a cluster experiment.

table[table omitted — 1,443 chars of source]

Optimization over the choice of the clustering

In this section, we turn to designing the optimal clustering.

As a first observation, we note that Theorem (ref) implies that any clustering solving the original objective in Equation (ref) must be on the Pareto frontier for the worst-case bias and variance as derived in Lemmas (ref) and (ref), respectively. Therefore, we first focus on drawing the frontier between the variance and bias. The optimization problem to draw the frontier between the worst-case bias and variance can be equivalently expressed in terms of the absolute instead of quadratic bias, so that we study

equation[equation omitted — 231 chars of source]

We are interested in minimizing $L_c(\xi)$ over a grid of values for $\xi$. This effectively will allow us to draw a Pareto frontier between the worst-case variance and absolute bias. Any clustering solving the original objective in Equation (ref) must minimize $L_c(\xi)$ for some value of $\xi$. Once we obtain the Pareto frontier, computing the clustering that minimizes our original objective in Equation (ref), becomes a simple task, as it only requires comparing the objective in Equation (ref) within the frontier.

Therefore, as a first step, we focus on minimizing $L_c(\xi)$ for each given $\xi$. Equation (ref) solves a weighted min-cut problem with an additional penalization term that depends on the variance of the clusters' size. The optimization problem in Equation (ref) differs from the minimum normalized-cut problem shi2000normalized, ling2020certifying, or related graph clustering problems pokhilko2019d, parker2017optimal, zhang2023asymptotics whose objective only depend on measure of cluster purity such as $ \sum_{k=1}^K \sum_{i: c(i) = k}|j \in \mathcal{N}_i: c(i)\neq c(j) |/ \sum_{i: c(i) = k}|\mathcal{N}_i|. $ Different from our proposal, standard min-cut problems use a different denominator (since the objective is not motivated by the bias of treatment effect estimators), and, most importantly, do not account for the additional variance component. Here instead, researchers should balance the bias and variance in the choice of the clusters. To the best of our knowledge, the problem in Equation (ref) has not been studied before.

Let $\mathbf{V} = \text{diag}(\mathbf{A}\mathbf{1}_{n}) = \text{diag}(|\mathcal{N}_{1}(1)|, \ldots, |\mathcal{N}_{n}(1)|)$. Define the left-normalized Laplacian

equation[equation omitted — 112 chars of source]

Given a set of $K$ clusters (i.e., fixing the number of clusters), for any cluster mapping $c: \{1, \cdots, n\} \mapsto \{1, \cdots, K\}$, let $\mathbf{M}_c(K)\in \mathbb{R}^{n\times K}$ with $$

aligned\mathbf{M}_{c, ik}(K) = 1\{c(i)=k\}, \quad i\in \{1, \cdots, n\}, \quad k\in \{1, \cdots, K\},

$$ where $\mathbf{M}_c(K)$'s dimension depends on $K$. Here, $\mathbf{M}_{c, ik}(K) = 1$ if unit $i$ is in cluster $k$.

thmLet $\mathrm{tr}(\cdot)$ denote the trace operator. $c^\star$ minimizes (ref) over the set of clusterings with exactly $K$ clusters if and only if $$ \small \begin{aligned} c^\star \in \mathrm{arg} \max_{c: \{1, \cdots, n\} \mapsto \{1, \cdots, K\}} \mathrm{tr}\Big((n \mathbf{L} - \xi \mathbf{1}_{n}\mathbf{1}_{n}^T) \mathbf{M}_c(K) \mathbf{M}_c^T(K) \Big). \end{aligned} $$
proof[Proof of Theorem (ref)] See Appendix (ref).

Theorem (ref) formalizes the optimization problem as a trace-optimization program for a fixed number of clusters $K$ and fixed $\xi$. Theorem (ref) does not characterize a convex optimization program, but it provides a natural starting point to study convex relaxations of the proposed optimization problem. To obtain a convex relaxation, we relax the constraint on the matrix $\mathbf{X}(K) = \mathbf{M}_c(K) \mathbf{M}_c^T(K)$. We propose solving the following semi-definite programming (SDP) problem (for a given number of clusters $K$):

equation[equation omitted — 304 chars of source]

where Equation (ref) defines a sequence of semi-definite optimization programs, each indexed by the number of clusters $K$ and multiplier $\xi$. The main distinction between Equation (ref), and Equation (ref) is that the matrix $\mathbf{X}(K)$ must be positive-definite instead of a matrix with binary entries. Let $\hat{\mathbf{X}}(K)$ be the solution of (ref), for given $K$, that can be obtained using off-the-shelf optimization routines. We then apply the $K$-means algorithm on the first $K$ eigenvectors of $\hat{\mathbf{X}}(K)$ to retrieve the mapping $c$, whose properties are well studied in the literature von2007tutorial. Finally, we compare solutions for different values of $K$ and choose the clustering with the largest objective.

Equation (ref) is a convex relaxation of the problem in Theorem (ref) because it substitutes the original constraint on $\mathbf{X}$ to contain binary entries with a semi-definite constraint. Such convex relaxations, together with $K$-means algorithm to retrieve the clusters, are common in the clustering literature and have been widely studied both from theoretical and numerical perspectives; see hong2021optimal for a review.

We can now describe the complete algorithm in Algorithm (ref). We first posit the grid of values for $\xi$ that serves to approximate the Pareto frontier and a grid of values on the number of clusters $K_n$ we want to consider. A more fined grid provides a better approximation, at the cost of the optimization program scaling linearly in the size of each grid. For each value of $\xi$ and number of clusters $K_n$, we then solve a semi-definite program using off-the-shelf computational methods. For each value of $\xi$, we compute the best number of clusters $K_n$ by also minimizing $L_c(\xi)$ over the number of clusters. The resulting estimated clustering provides us with the set of clustering on a discretized Pareto frontier (approximated through the grid for $\xi$). Finally, we compute the mean-squared error optimal clustering by comparing the original objective in Equation (ref) across all such estimated clustering over the grid of values for $\xi$ and $K_n$. The resulting optimization program is a simple sequence of semi-definite programs.

algorithm[algorithm omitted — 1,837 chars of source]
rem[Spectral relaxation and comparison to spectral clustering] Unlike the minimum normalized cut, there is no simple spectral relaxation of the optimization problem in Theorem (ref), unless all clusters are equally sized. For the special case where all clusters are equal-sized, (ref) can be relaxed to $$ \small \begin{aligned} \max_{\mathbf{U}} \text{tr}(\mathbf{L}_{\xi} \mathbf{U}(K) \mathbf{U}^T(K)), \quad \text{s.t.}\,\, \mathbf{U}^T(K) \mathbf{U}(K) = \mathbf{I}_{K}, \end{aligned} $$ where $\mathbf{I}_{K}$ denotes the identity matrix of dimension $K$. We first symmetrize \begin{equation} \begin{aligned} tr(\mathbf{L}_{\xi} \mathbf{U}(K) \mathbf{U}^T(K)) = tr\left(\frac{\mathbf{L}_{\xi} +\mathbf{L}_{\xi}^T}{2} \mathbf{U} \mathbf{U}^T\right) \end{aligned} \end{equation} The solution to the above problem $\hat{\mathbf{U}}(K) \in \mathbb{R}^{n\times K}$ is given by the matrix of top-$K$ eigenvectors of $(\mathbf{L}_{\xi} + \mathbf{L}_{\xi}^T) / 2$. Then we can perform a constrained $K$-means algorithm with equal-sized clusters to recover the clusters bradley2000constrained. Finally, we iterate over the grid of values for $(\xi,K)$ to retrieve the optimal clustering. This special case allows us to illustrate sharp distinctions with standard clustering methods such as spectral clustering von2007tutorial. Once we fix the number of clusters $K$, the first distinction is which matrix we use to retrieve the clusters via K-means clustering. Spectral clustering would estimate the clustering through an eigendecomposition of the Laplacian matrix. Here instead, we use an eigendecomposition of the different matrix $\mathbf{U} \mathbf{U}'$ that is a function of a left-normalized Laplacian and the matrix capturing the variance contribution of the estimated effect (implicit in $\mathbf{L}_{\xi}$). Instead of fixing the number of clusters, we then choose the clustering (and number of clusters) based on our different objective. \qed
rem[Reducing the computational complexity with a surrogate relaxation] To reduce the complexity of Algorithm (ref) it is possible to consider a single value of $\xi$ (instead of a grid of values as described in Algorithm (ref)) to minimize the mean-squared error. In particular, we can write $ \frac{1}{n^2} \sum_{k=1}^{K_n}n_k^2 + \frac{\lambda \bar{\phi}_n^2}{\bar{\psi}} \Big(\frac{1}{n} \sum_{i=1}^n \frac{1}{|\mathcal{N}_i|} \Big|j \in \mathcal{N}_i: c(i) \neq c(j) \Big|\Big)^2 \le \frac{1}{\xi} L_c\left(\xi\right), $ for $\xi = \frac{\bar{\psi}}{\lambda \bar{\phi}_n^2}$, where $\frac{1}{\xi} L_c\left(\xi\right)$ is monotonically decreasing in $\xi$. Therefore, minimizing $\frac{1}{\xi} L_c(\xi)$ for any $\xi \le \frac{\bar{\psi}}{\lambda \bar{\phi}_n^2}$ serves as a surrogate objective for the mean-squared error objective, where such surrogate imposes a larger weights to the bias relative to the variance. \qed

Empirical illustration and numerical studies

In this section, we illustrate the properties of the procedure in two empirical applications, one using unique data from the Facebook friendship and messaging graphs and one using data from an experiment conducted in rural China by cai2015social. We provide further evidence supporting our theoretical analysis using simulated networks in Appendix (ref).

Comparisons between existing clusterings at Faceboook

First, we focus on evaluating existing clusters currently implemented on Facebook at scale.\footnote{We focus on existing clusterings (instead of implementing the causal clustering proposed here), as implementing new clustering methods requires changes to the entire engineering pipeline as described in karrer2021network, outside the scope of this paper.} In Section (ref) we turn to the implementation of Causal Clustering.

Bias and variance computations We consider two clustering algorithms implemented at scale on social networks owned by Meta: Louvain algorithm blondel2008fast, and Balanced Partitioning kabiljo2017social. Louvain and Balanced Partitioning clustering produces a hierarchical clustering structure with a growing number of clusters. We consider three “types" of such algorithms, corresponding to three levels in the clustering structure hierarchy (i.e., different numbers of clusters), defined as Type 1, 2, and 3, respectively. For Louvain, higher-order types denote more clusters, whereas it is the opposite for Balanced Partitioning. For each clustering, we also report $\log(n/K_n)$ in Table (ref).

We consider two graphs owned by Meta. In the first graph, edges capture the strength of the friendship in the Facebook graph, and in the second graph, they capture connections based on messaging on Facebook. In both cases, the data was aggregated and de-identified. In each graph, edges are continuous variables. For each graph, we construct three adjacency matrices by setting the edges to be zero if below the $5^{th}$, $10^{th}$, and $50^{th}$ percentile value of the edges. We will refer to these graphs obtained as dense, moderate (mod), and sparse after thresholding. We also note that all graphs have a bounded maximum degree.

We refer to the worst-case bias and variance weighted by $\xi$ as $ b_n(\mathcal{C}_n), \frac{\xi}{n^2} \sum_{k=1}^{K_n} n_k^2. $

Results Table (ref) in the appendix collects each of these algorithms' worst-case bias and variance and the number of clusters per individual. A larger number of clusters increases the bias but decreases the variance throughout all graphs. (The bias increases for sparser networks because the denominator $|\mathcal{N}_i|$ decreases.) Louvain algorithm dominates Balanced Partitioning uniformly in variance and bias, suggesting that practitioners may prefer Louvain.

We study different types of Louvain algorithms (each having a different number of clusters) in Figure (ref). Figure (ref) reports the worst-case mean squared error for different graphs, Louvain algorithms, and degrees of sparsity. For $\xi \approx 1$, the Louvain algorithm with the largest number of clusters (Type 3) dominates the other two algorithms in all but one case.

As $\xi$ increases, Louvain clustering with the largest number of clusters (Type 3) is optimal for a dense graph but sub-optimal as the graph becomes more sparse. This result is intuitive: for denser networks, a larger number of clusters may best control the estimator's bias relative to the bias induced by other clustering algorithms. However, we observe that Type 1 is the most robust in terms of variance in the mean-squared error across different values of spillover effects. This is suggestive that Louvain clustering with the largest number of clusters (Type 1) may be preferred in various applications.

Comparison between Bernoulli and Cluster design To conclude, we study when a cluster design with our preferred clustering (Louvain Type 1) should be preferred over a Bernoulli design taking $\lambda = 1$. Using Theorem (ref), and the dense Messaging graph, Theorem (ref) for Louvain clustering (Type 1), for $ \sqrt{\frac{\bar{\psi}}{\bar{\phi}_n^2}} \le \bar{\xi}^{1/2}, \quad \bar{\xi}^{1/2} = 706, $ we should run a cluster experiment instead of a Bernoulli design. Therefore, researchers should run a Cluster design if $ \bar{\phi}_n^2 \ge \bar{\psi}/\bar{\xi}. $ This implies that the researcher should run a cluster design if spillovers are larger than $\bar{\phi}_n \ge \bar{\psi}^{1/2}/706$. This comparison is suggestive that a cluster design with a Louvain clustering may be preferred over a Bernoulli design over a wide range of values of spillover effects (and outcomes variation $\bar{\psi}$). For binary outcomes, we should prefer cluster experiments for $\bar{\phi}_n \ge 0.003$.

In summary, our results shed light on using clustering algorithms for large-scale implementation on online platforms. These results suggest that Louvain clustering with possibly many clusters performs best in practice, and may often be preferred to a Bernoulli design. It confirms heuristic arguments in karrer2021network based on AA-tests.

\definecolor{glaucous}{rgb}{0.38,0.51,0.71}

figure[figure omitted — 649 chars of source]

Clustering in the field

figure[figure omitted — 516 chars of source]

Next, we present an application of a field experiment where the treatment consists of informing individuals to increase insurance take-up cai2015social.

cai2015social collected network information in 185 villages in rural China from 48 larger regions. We use network data collected by cai2015social to study the properties of the proposed method. The network has a total of $7649$ nodes, once we also include individuals who are friends of surveyed individuals. Individuals between villages present on average $50\%$ of their connections outside their village. On the other hand, individuals in different regions (with $48$ regions) have $99\%$ of their connections within the same region. We use an adjacency matrix a “weak ties" adjacency matrix, where two individuals are connected if either indicates the other has a friend.

We construct a clustering in each of the 48 regions (since individuals are not connected between different regions). We estimate optimal clustering over the Pareto frontier of values of $\xi$ as in Algorithm (ref).

We compare the procedure (denoted as “Causal Clustering") with the following alternative clustering algorithms: $\varepsilon$-net clustering as in eckles2017design with $\varepsilon = 3$ as suggested in eckles2017design; spectral clustering with a fixed number of clusters equal to $n/3$, where $n$ is the population size in a given region (note that spectral clustering do not optimize over the number of clusters); Louvain clustering with default parameters selected by the R-package igraph; clustering based on village identity of each individual (creating one additional cluster for those individuals whose village is missing in the data).\footnote{Note that all these algorithms are off-the-shelf and none of these clustering have been previously studied in the context of mean-squared error reduction for the overall effect as studied in this paper. Therefore, whenever possible we select the tuning parameters based on default implementations either suggested by the literature, as for $\varepsilon$-clustering, or by default in standard packages such as for Louvain clustering. }

figure[figure omitted — 860 chars of source]

Main results

Next we present results for the Pareto frontier and weighted mean-squared error.

Pareto Frontier We contrast the objective function with alternative algorithms on the right-hand side of Figure (ref). The proposed method achieves the smallest objective function across all competitors over a wide range of $\xi$, providing suggestive evidence that Algorithm (ref) (with a semi-definite relaxation) works well in practice. Interestingly, clustering based on village identity leads to the largest loss function. This is because individuals are connected within and between villages in this application.

Therefore, Figure (ref) (right-panel) shows that each of the existing clustering methods are Pareto dominated by the proposed method for the grid of values of $\xi$ considered.

In addition, our conclusions suggest that we should choose a larger number of clusters compared to the overall number of villages. To gain further intuition of its practical implication, take $\bar{\psi} = 0.24$ approximately equal to the outcomes' variance in cai2015social and $\bar{\phi}_n= 0.27$ as in Table 2 in cai2015social. Following Remark (ref), a conservative choice of $\xi$ (i.e., a choice that favors a smaller bias) is approximately $\bar{\psi}/\bar{\phi}_n^2 = 3.29$. Even with such a conservative choice, the suggested number of clusters is large relative to the population size. In this case, the optimal number of clusters is around $15\%$ of the overall population size, approximately $600$ clusters with a surveyed sample size of $3600$ individuals. This result is suggestive of the relevance of our method to optimally choose the number of clusters, taking into account the downstream causal inference objective.

Weighted mean-squared error comparisons Next, we evaluate the performance of each methods in terms of mean-squared error. We take each clustering estimated from the real-world network as described above. We then simulate treatments and outcomes following the estimated model in cai2015social (Table 2, Column 4), which assumes that individual outcomes depend on neighbors' treatments, consistently with our Assumption (ref). We vary the variance of the residuals $\sigma^2 \in \{1/4, 1/2,1\}$ to emulate settings with high, medium, and low signal-to-noise ratios. We then compute

equation[equation omitted — 161 chars of source]

as the sum of the squared bias and variance, as a function of the parameter $\lambda$. We compare the performance of Algorithm (ref) implemented to minimize $L(\xi)$ in Equation (ref), with $\xi = \lambda^{-1}$ as described in Remark (ref). This comparison is useful for two reasons. First, by selecting $\xi = \lambda^{-1}$, $L(\xi)$ serves as a surrogate for the underlying objective, where such surrogate uses an upper bound $\bar{\phi}_n^2/\bar{\psi} \le 1$. Second, because we compare the objective in Equation (ref) we contrast the squared instead of absolute bias with the variance.\footnote{More refined comparisons where we would also optimize over a grid of $\xi$ for each $\lambda$ instead of using an upper bound for $\xi$ are possible but omitted for brevity.}

In Figure (ref) (first three panels), we report $M(\lambda)$. Our method uniformly outperforms competitors regarding mean-squared error, except for a slight underperformance relative to the Louvain algorithm when the signal-to-noise ratio is particularly high. However, as the signal to noise ratio is medium or low, our methods uniformly outperforms all others.

In the last panel of Figure (ref), we consider settings where the local interference assumption is violated, and we calibrate our model to the specification in cai2015social (Column 4, Table 5), where the individual outcome depends on the neighbors' outcomes. Although this is not the main specification in cai2015social, such calibration is useful to measure robustness to misspecification of first-degree neighbors' interference. We fix the residuals' variance to be $\sigma^2 = 1/4$, whereas the results are robust as we increase $\sigma^2$. In this setup, we observe a larger mean-squared error of all methods, with the proposed method achieving the lowest mean-squared error.

Extensions

Here we present a brief description of our extensions, and defer to Appendix (ref) details.

Experimental design for bias-aware inference

In this section, we turn to settings where researchers are interested in conducting inference on the global treatment effect $\tau_n$, and focus on minimizing the length of the confidence interval. First, define

equation[equation omitted — 168 chars of source]

where

equation[equation omitted — 180 chars of source]

In addition, define $$

aligned\beta_n(\mu) = \tau_{n,\mu} - \mathbb{E}_\mu[\hat{\tau}_n(\mathcal{C}_n)], \quad \nu_n(\mu) = \sqrt{\mathbb{E}_\mu\left[\Big(\mathbb{E}_\mu[\hat{\tau}_n(\mathcal{C}_n)] - \hat{\tau}_n(\mathcal{C}_n)\Big)^2\right]}, \quad b_{*,n} = \bar{\phi}_n b(\mathcal{C}_n),\quad \sigma_{*,n}^2 = \frac{\bar{\psi}}{n^2} \sum_{k=1}^{K_n} n_k^2.

$$ Here, $\mathbb{E}_\mu[\cdot]$ indicates the expectation operator under potential outcomes $Y_i(\mathbf{d})= \mu_i(\mathbf{d}_i, \mathbf{d}_{\mathcal{N}_i})$. Here $b_{*,n}$ and $\sigma_{*,n}^2$ correspond to the worst-case bias and variance derived in Section (ref).

lemSuppose that Assumptions (ref), (ref), (ref) hold. Let $\mathcal{N}_{n, \max} = o(K_n^{1/4})$. For any $\mu$ and sequences $\bar{b}_n, \bar{\sigma}_n$ such that $|\beta_n(\mu)| \le \bar{b}_n$, $\nu_n(\mu) \le \bar{\sigma}_n$, with $\bar{\sigma}_n^2 \ge \frac{\eta + o(1)}{K_n}$ for some $\eta > 0$, it follows that $ \lim_{n \rightarrow \infty} \mathbb{P}(\tau_{\mu, n} \in \mathcal{I}_{\mathrm{ea}}(\bar{b}_n, \bar{\sigma}_n)) \ge 1 - \alpha. $
proofSee Appendix (ref).

Lemma (ref) establishes that the confidence interval $\mathcal{I}_{\mathrm{ea}}(\bar{b}_n, \bar{\sigma}_n)$ guarantees valid coverage pointwise for any $\mu$, provided that $\bar{b}_n$ and $\bar{\sigma}_n$ are valid upper bounds on the bias and standard deviation under a potential outcome model $\mu$. In addition, Lemma (ref) requires that such upper bound on the variance $\bar{\sigma}_n^2$ scales at order $1/K_n$ (or slower). This condition is attained when observations are independent or positively dependent and also attained by the worst-case variance derived in Lemma (ref), under a local asymptotic framework with $\varepsilon = o(1)$ in Assumption (ref) (i.e., spillover effects are local to zero).

We optimize over an objective of the form $ \sup_{\mu \in \mathcal{M}} |\mathcal{I}_{\mathrm{ea}}(\beta_n(\mu), \nu_n(\mu)|, $ i.e., we study the worst-case length of the confidence interval once we use the true bias and standard deviation as a function of $\mu$.

thmSuppose that Assumptions (ref), (ref), (ref), (ref) hold. Then $$ \small \begin{aligned} \lim_{n\rightarrow \infty} \frac{\sup_{\mu \in \mathcal{M}} |\mathcal{I}_{\mathrm{ea}}(\beta_n(\mu), \nu_n(\mu)|}{|\mathcal{I}_{\mathrm{ea}}^*|(1 + \mathcal{O}(\varepsilon))} = 1, \quad |\mathcal{I}_{\mathrm{ea}}^*| = 2|\chi(b_{*,n}, \sigma_{*,n})| \end{aligned} $$
proofSee Appendix (ref).

Theorem (ref) shows that the length of the confidence intervals using the worst-case approximations for the bias and variance derived in Section (ref) approximates the worst-case length of the confidence interval. Given this result, we can then derive a simple clustering algorithm to design the clusters with minimum worst-case confidence interval length.

The algorithm proceeds as follows. First, we compute the Pareto frontier between worst-case variance and bias as in Algorithm (ref) (Step 1 - 5). Clearly, the algorithm that minimizes $|\mathcal{I}_{\mathrm{ea}}^*|$ must be on such frontier. Second, instead of choosing between the clusters $\mathcal{C}(\xi), \xi \in \mathcal{S}$ on the frontiers, the one with the smallest mean-squared error, we choose the one with the smallest $|\mathcal{I}_{\mathrm{ea}}^*|$. As for Algorithm (ref), the algorithm solves a sequence of semi-define programs.

Finally, note that once the experiment is conducted, it is possible to conduct inference on $\tau_\mu$ by using an estimator for the variance (instead of a worst-case bound), and invoke Lemma (ref) for inference. Such estimator can be constructed using, for instance, the HAC estimator in gao2023causal. We provide details in Appendix (ref).

Additional extensions

In Appendix (ref) we consider a series of other extensions. These are: (i) potential outcomes are random variables and we impose restrictions only on the moments of such random variables; (ii) spillover effects may be asymmetric between different individuals; (iii) researchers consider a larger class of estimators and optimize both over designs and estimators; (iv) saturation designs with positive treatment probabilities.

Conclusions and recommendations for practice

This paper provides a formal characterization of the bias and variance of the clustering algorithm to estimate the overall treatment effect. Using such characterization, it introduces a novel clustering method that takes into account the downstream task of causal inference.

The algorithm requires a choice of the adjacency matrix and a sensible prior value of the largest magnitude of the spillover effect. The choice of the adjacency matrix must depend on researchers' prior knowledge of which dimension spillovers propagate over, as, for example, in cai2015social. Our framework directly extends to settings where researchers only have some prior distribution on the network, as in breza2020using, see Remark (ref).

Similarly to standard power analysis baird2018optimal, our method with $\lambda = 1$ also depends on the choice of $(\bar{\phi}_n^2/\bar{\psi})^{-1}$ (the size of spillover effects $\bar{\phi}_n$ relative to the outcomes' largest squared deviation $\bar{\psi}$). In our empirical applications, we observe that certain clustering algorithms uniformly outperform others for a large range of values of such values.

To choose such a range of values, our recommended choice for $\bar{\psi}$ is to consider values $\bar{\psi} = c \bar{\sigma}^2, c \in [1,4]$ and $\bar{\sigma}^2$ is the baseline variance of the residuals after removing the covariate adjustment, and $c$ is a constant between one to four.\footnote{To gain further intuition, let $\lambda = 1$. Let $\bar{\mu}_i$ be a prediction for $\mu_i(\mathbf{0})$ as in Remark (ref) and consider an estimator as in Equation (ref). When using regression adjusted estimators, we have $\bar{\psi} = \Big(\mu_i(\mathbf{1}) - \bar{\mu}_i + \mu_i(\mathbf{0}) - \bar{\mu}_i\Big)^2$. By approximating $(\mu_i(\mathbf{d}) - \bar{\mu}_i)^2 \approx \bar{\sigma}^2$, we can write $\bar{\psi} \le 4 \bar{\sigma}^2$.}

The choice of $\bar{\phi}_n$ can be based on spillover effects observed in previous experiments, in the spirit of minimum detectable effects baird2018optimal, or using a priori bound.

Finally, many open questions remain. As we discuss in the main text, local interference provides us with an explicit objective function that can be optimized (as opposed to only studying rates of convergence). Future research should study explicit optimal clustering with non-local network interference as for example in bramoulle2009identification, leung2022causal among others. Topics of future research also include clustering in dynamic settings, or optimal clustering for a larger class of estimators.

\spacingset{1.6}