The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
79,215 characters
Causal clustering: design of cluster experiments under network interference
\if00
{
\title{\bf Causal clustering: design of cluster experiments under network interference\footnote{We thank Toru Kitagawa, Michael Leung, Jesse Shapiro, Jann Spiess, Max Tabord-Meehan for helpful comments. All mistakes are our own.} }
\author{Davide Viviano\footnote{Department of Economics, Harvard University. Email: [email removed]} $\quad$ Lihua Lei\footnote{Graduate School of Business, Stanford University. Email: [email removed]} $\quad$ Guido Imbens\footnote{Graduate School of Business and Department of Economics, Stanford University. Email: [email removed]} \\ Brian Karrer\footnote{FAIR, Meta. Email: [email removed] } $\quad$ Okke Schrijvers\footnote{Central Applied Science, Meta. Email: [email removed].} $\quad$ Liang Shi\footnote{Central Applied Science, Meta. Email: [email removed].}}
\date{This version: \today \\ First version: October, 2023}
\maketitle
} \fi
\if10
{
\bigskip
\bigskip
\bigskip
\
\\
\
\\
\medskip
} \fi
\begin{abstract}
\noindent
This 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 \textit{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.
\end{abstract}
\noindent
{\it Keywords:} Experimental Design, Spillover Effects, Causal Inference, Cluster Designs. \\
{\it JEL Codes:} C10, C14, C31, C54
\vfill
\newpage
\bigskip
\spacingset{1.5}
\section{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 \cite{egger2022general} for an empirical example.}
This paper provides an econometric framework to choose when and how to design the \textit{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 \textit{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 \citep[e.g.,][]{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 \citep[e.g.][]{baird2018optimal}, or using some particular modeling assumptions. We provide a practical discussion in Section \ref{rem:psi}.
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 \citep[e.g.][]{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 \cite{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 \cite{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 \citep{ugander2013graph, eckles2017design}. Variants of these algorithms have been recently studied in \cite{leung2022rate} \citep[and more recently in contemporaneous work by][]{leung2023design} for spatial spillovers, and \cite{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 \cite{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 \textit{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 \citep[e.g.][]{brennan2022cluster, taylor2018randomized}, or on rate-optimal designs only \citep{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 \textit{absence} of cluster experiments are \cite{basse2018model}, and \cite{jagadeesan2020designs} for estimating direct instead of global average treatment effects studied here; \cite{kang2016peer} who study encouragement designs, without focusing on the problem of mean-squared error optimal designs; \cite{viviano2020experimental} for experiments on networks using information from a pilot study; \cite{basse2018limitations} discuss limitations of design-based causal inference under interference. Recent econometric literature studying experimental design focusing however on $i.i.d.$ data include \cite{cai2024performance, tabord2023stratification, armstrong2022asymptotic} among others.
None of these papers study (optimal) cluster designs.
The literature on treatment effects under network interference includes \cite{aronow2017estimating}, \cite{hudgens2008toward}, \cite{toulis2013estimation}, \cite{manski2013identification}, \cite{leung2020treatment}, \cite{athey2018exact}, \cite{goldsmith2013social}, \cite{savje2017average}, \cite{ogburn2017causal}, \cite{manresa2013estimating}, \cite{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 \cite{wooldridge2003cluster}, \cite{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 \cite{von2007tutorial}, \cite{newman2013community,newman2013spectral}, \cite{lei2019unified}, \cite{lei2020consistency}, \cite{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{sec:setup} presents the setup; Section \ref{sec:when} characterizes the bias and variance and contrast cluster and Bernoulli designs; Section \ref{sec:opt} formalizes the optimization problem for the optimal clustering; Section \ref{sec:empirics} presents two applications and numerical studies, Section \ref{sec:extensions} presents several extensions, and Section \ref{rem:psi} contains recommendations for practice.
\section{Setup} \label{sec: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{rem:unobserved_a} 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,
\begin{equation} \label{eqn:overall}
\small
\begin{aligned}
\tau_n = \frac{1}{n} \sum_{i=1}^n \Big[Y_i(\mathbf{1}) - Y_i(\mathbf{0})\Big].
\end{aligned}
\end{equation}
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., \cite{egger2022general}, \cite{muralidharan2017experimentation}, \cite{karrer2021network} (see Appendix \ref{sec:saturations} for a discussion on other estimands). Table \ref{tab:methods3} summarizes the notation.
\vspace{-4mm}
\subsection{Potential outcomes and spillover effects}
Next, we impose restrictions on the spillover effects.
\begin{ass}[First-order local interference] \label{ass:first_order} 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}$.
\end{ass}
Assumption \ref{ass:first_order} 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., \cite{cai2015social}, \cite{sinclair2012detecting}, \cite{muralidharan2017experimentation}. \cite{athey2018exact} provide a framework for testing Assumption \ref{ass:first_order}. Higher-order interference can be accommodated, although higher-order effects are often small and difficult to detect (see Remark \ref{rem:higher_order}). 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 \eqref{eqn:overall} 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 \textit{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.
\begin{ass}[Class of potential outcome functions] \label{ass:exposure_restriction} 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 \textit{all} functions $\mu_i(d, \cdot)$ satisfying the following two conditions:
\begin{itemize}
\item[(i)] $\mu_i(d, \cdot)\in [\psi_L, \psi_R]$ for some constants $-\infty < \psi_L < \psi_R < \infty$ that do not vary with $n$;
\item[(ii)] for all $\mathbf{d} \in \{0,1\}^{|\mathcal{N}_i|}$, for some $\bar{\phi}_n \le \psi_R - \psi_L$
\begin{equation} \label{eqn:interference}
\small
\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}
\end{ass}
Assumption \ref{ass:exposure_restriction} imposes two restrictions. Condition (i) states that potential outcomes are uniformly bounded.\footnote{This restriction is common in the literature, e.g., \cite{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 \textit{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{sec:random_outcomes}, we relax Assumption \ref{ass:exposure_restriction}, by imposing $(ii)$ only on the expected value of the potential outcomes (instead of the realized potential outcomes). In Appendix \ref{sec:heterogeneity}, we relax Assumption \ref{ass:exposure_restriction} 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.
\begin{ass}[Spillover local asymptotics] \label{defn:local} 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$.
\end{ass}
Assumption \ref{defn:local} 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{defn:local} can be sharpened by imposing restrictions on
the \textit{average} second-order degree instead of maximum degree; see Equation \eqref{eq:sum_square_degree} in the proof of Lemma \ref{lem:rel}.} 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 \citep[e.g.][where $\mathcal{N}_{n, \max}$ is bounded]{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 \citep[e.g.][]{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{rem:direct_asymptotics}.
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 \cite{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 \citep[see also][for a discussion]{athey2023semi}. Similar examples in digital advertisements with small effects are
\cite{athey2023digital} with $i.i.d.$ data and \cite{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 \cite{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.
\begin{exmp}[Linear exogenous peer effects] \label{exmp:1}
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{ass:exposure_restriction} holds. \qed
\end{exmp}
\begin{rem}[Direct and overall effect] \label{rem:direct_asymptotics} Our asymptotic framework also allows direct treatment effects (and global effects) to be of the \textit{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, \cite{karrer2021network}). \qed
\end{rem}
\begin{rem}[Higher order interference and endogenous peer effects] \label{rem:higher_order}
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{defn:local} 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{example:endogenous}.
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
\end{rem}
\vspace{-4mm}
\subsection{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
$$
\small
\begin{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 \text{ for } k \neq k' \right\}.
\end{aligned}
$$
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{sec:opt}.
\begin{ass}[Cluster designs] \label{ass:clusters} $\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.
\end{ass}
Assumption \ref{ass:clusters} 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{ass:clusters} 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 \citep[which are those used in most of applications on online platforms,][]{karrer2021network}; see Remark \ref{rem:estimator} 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.}
\begin{equation} \label{eqn:estimator2}
\small
\begin{aligned}
\hat{\tau}_n(\mathcal{C}_{ n}) = \frac{2}{n} \sum_{i=1}^n \Big[D_i Y_i - (1 - D_i) Y_i \Big] = \frac{2}{n}\sum_{i=1}^{n}(2D_i - 1)Y_i.
\end{aligned}
\end{equation}
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}$ \textit{only} because the distribution of the treatments depends on the clusters under Assumption \ref{ass:clusters}. Studying the estimator in Equation \eqref{eqn:estimator2} 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{rem:reg_adj}) are often used or studied in practice \citep{karrer2021network, holtz2020reducing, savje2017average}. One could also normalize each sum in Equation \eqref{eqn:estimator2} 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).
\begin{rem}[Alternative estimators] \label{rem:estimator}
Alternative estimators studied in the literature are inverse probability weights estimators \citep[e.g.,][]{aronow2017estimating, ugander2013graph}. These estimators reweight by the probability that the treatment of \textit{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 \textit{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 \textit{choice of the estimators and designs}. We study this class of estimators and designs in Appendix \ref{sec:other_estimators}.
\qed
\end{rem}
\begin{rem}[Covariate adjustment] \label{rem:reg_adj} 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} \label{eqn:reg_adj}
\small
\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{ass:exposure_restriction}. To account for this, we can apply the more general result in Appendix \ref{sec:heterogeneity}.
\qed
\end{rem}
\begin{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{sec:saturations}. \qed
\end{rem}
\vspace{-5mm}
\section{Bias, Variance and MSE of a given clustering } \label{sec:when}
As a necessary first step, this section focuses on studying the bias, variance, and mean-squared error of a \textit{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{ass:clusters}.
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
\begin{equation} \label{eqn:worst-case}
\small
\begin{aligned}
\mathcal{B}_n(\mathcal{C}_{ n}, \lambda) = \sup_{\mu \in \mathcal{M}} \left\{\mathbb{E}_\mu\Big[\Big(\hat{\tau}_n(\mathcal{C}_{ n}) - \mathbb{E}_\mu[\hat{\tau}_n(\mathcal{C}_{ n})]\Big)^2\Big] + \lambda \Big(\tau_{n,\mu} - \mathbb{E}_\mu[\hat{\tau}_n(\mathcal{C}_{ n})]\Big)^2 \right\},
\end{aligned}
\end{equation}
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{ass:exposure_restriction}. 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 \textit{length of confidence intervals} for the overall effects, which we defer to Section \ref{sec:bias_aware_inferece}.
\vspace{-3mm}
\subsection{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
\begin{equation}\label{eq:bn}
\small
\begin{aligned}
b_n(\mathcal{C}_n) = \frac{1}{n}\sum_{i=1}^n\frac{1}{|\mathcal{N}_i|} \Big|\mathcal{N}_i \bigcap \Big\{j: c(i) \neq c(j)\Big\} \Big|
\end{aligned}
\end{equation}
and $\Big\{j: c(i) \neq c(j)\Big\}$ denotes the set of units $j$ in a different cluster from unit $i$.
\begin{lem}[Worst-case bias] \label{lem:worst_case_bias} Let Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{ass:clusters} 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).
$
\end{lem}
\begin{proof} See Appendix \ref{sec:proof_bias}.
\end{proof}
Lemma \ref{lem:worst_case_bias} 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{lem:worst_case_bias} 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 \eqref{eqn:interference} holds with equality (Appendix \ref{sec:heterogeneity} studies settings where effects are heterogeneous across individuals). The expression for the bias connects to measures of cluster purity such as those in \cite{goel2009respondent}, \cite{toulis2013estimation}, motivating clustering purity based on the class of models in Assumption \ref{ass:exposure_restriction}.
To derive the worst-case variance, we define
\begin{equation}\label{eq:psibar}
\small
\begin{aligned}
\bar{\psi}^{1/2} = 2\max\{|\psi_R|, |\psi_L|\} .
\end{aligned}
\end{equation}
By Assumption \ref{ass:exposure_restriction} (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$.
\begin{lem}[Worst-case variance] \label{prop:main}
Let Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{defn:local}, \ref{ass:clusters} 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}
$$
\end{lem}
\begin{proof} See Appendix \ref{app:propmain}.
\end{proof}
Lemma \ref{prop:main} 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{defn:local}, 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.
\begin{thm} \label{thm:Bn_Bnstar} Suppose that Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{defn:local}, \ref{ass:clusters} 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}\label{eqn:bound2}
\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}
\end{thm}
\begin{proof}
See Appendix \ref{app:worst_case_mse}.
\end{proof}
\begin{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{sec:random_outcomes}, 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 \textit{moments} of potential outcomes. \qed
\end{rem}
\begin{rem}[Unobserved $\mathbf{A}$] \label{rem:unobserved_a} Suppose that $\mathbf{A}$ is unobserved or \textit{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 \citep[][]{breza2020using}. In this case, the sparsity restrictions in Assumption \ref{ass:first_order} can be replaced by assumptions on the expected squared degree over the prior -- see Equation \eqref{eq:sum_square_degree} in the proof of Lemma \ref{lem:rel}.
\qed
\end{rem}
\vspace{-3mm}
\subsection{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\}\}$.
\begin{thm}[Rule of thumb] \label{thm:thumb} Suppose that Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{ass:clusters} 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}
$$
\end{thm}
\begin{proof}[Proof of Theorem \ref{thm:thumb}] See Appendix \ref{app:thumb}.
\end{proof}
Theorem \ref{thm:thumb} 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 \textit{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{fig:min_size} illustrates the exact threshold for the spillovers that justifies a cluster experiment.
In our empirical application in Section \ref{sec:empirics} we find $\bar{\phi}_n \ge 0.003$ would motivate a cluster experiment.
\begin{table}[htb]
\caption{Minimum spillover effects' size that justifies running a cluster experiment instead of a Bernoulli design in the presence of equally sized clusters $(\gamma_k = 1$ for all $k \in \{1, \cdots, K_n\})$, and where the clustering $\mathcal{C}_n$ is given. The bias denotes $\frac{1}{n} \sum_{i=1}^n \frac{1}{|\mathcal{N}_i|} |j \in \mathcal{N}_i: c(i) \neq c(j)| \in [0,1]$. The table on the left-hand side reports the value of $\bar{\phi}_n$ for different values of the bias and number of clusters, fixing $\bar{\psi} = 4$ (e.g., valid when outcomes are bounded between $[-1,1]$). In the figure on the right-hand side, the y-axis reports in log-scale the minimum size of the spillover effects $\bar{\phi}_n/\sqrt{\bar{\psi}}$ to run a cluster experiment as in Theorem \ref{thm:thumb}. The x-axis reports the number of clusters. } \label{fig:min_size}
\begin{tabularx}{\linewidth}{*{2}{>{\centering\arraybackslash}X}}
\begin{tabular}[b]{cccc}
Bias & $K_n$ & $\bar{\psi}$ & $\bar{\phi}_n$ \\
\hline
0.25 & 100 & 4 & 0.21 \\
0.5 & 100 & 4 & 0.23 \\
0.75 & 100 & 4 & 0.3 \\
0.25 & 200 & 4 & 0.15 \\
0.5 & 200 & 4 & 0.16 \\
0.75 & 200 & 4 & 0.21 \\
0.25 & 500 & 4 & 0.09 \\
0.5 & 500 & 4 & 0.10 \\
0.75 & 500 & 4 & 0.14 \\
\hline
\end{tabular}
&
\includegraphics[scale = 0.3]{figures/min_spillovers.eps} \\
&
\end{tabularx}
\end{table}
\vspace{-15mm}
\section{Optimization over the choice of the clustering} \label{sec:opt}
In this section, we turn to designing the optimal clustering.
As a first observation, we note that
Theorem \ref{thm:Bn_Bnstar} implies that any clustering solving the original objective in Equation \eqref{eqn:worst-case} must be on the Pareto frontier for the worst-case bias and variance as derived in Lemmas \ref{lem:worst_case_bias} and \ref{prop:main}, 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
\begin{equation} \label{eqn:optimization_v2}
\small
\begin{aligned}
L_c(\xi) = \frac{\xi}{n^2} \sum_{k=1}^{K_n}n_k^2 + \frac{1}{n} \sum_{i=1}^n \frac{1}{|\mathcal{N}_i|} \Big|j \in \mathcal{N}_i: c(i) \neq c(j) \Big| .
\end{aligned}
\end{equation}
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 \eqref{eqn:worst-case} must minimize $L_c(\xi)$ for \textit{some} value of $\xi$. Once we obtain the Pareto frontier, computing the clustering that minimizes our original objective in Equation \eqref{eqn:worst-case}, becomes a simple task, as it only requires comparing the objective in Equation \eqref{eqn:worst-case} within the frontier.
Therefore, as a first step, we focus on minimizing $L_c(\xi)$ for each given $\xi$.
Equation \eqref{eqn:optimization_v2} solves a weighted min-cut problem with an additional \textit{penalization} term that depends on the variance of the clusters' size. The optimization problem in Equation \eqref{eqn:optimization_v2} differs from the minimum normalized-cut problem \citep[e.g.,][]{shi2000normalized, ling2020certifying}, or related graph clustering problems \citep[see][and references therein]{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 \textit{and variance} in the choice of the clusters.
To the best of our knowledge, the problem in Equation \eqref{eqn:optimization_v2} 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
\begin{equation}\label{eq:laplacian}
\small
\begin{aligned}
\mathbf{L} = \mathbf{V}^{-1}\mathbf{A}.
\end{aligned}
\end{equation}
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
$$
\small
\begin{aligned}
\mathbf{M}_{c, ik}(K) = 1\{c(i)=k\}, \quad i\in \{1, \cdots, n\}, \quad k\in \{1, \cdots, K\},
\end{aligned}
$$
where $\mathbf{M}_c(K)$'s dimension depends on $K$. Here, $\mathbf{M}_{c, ik}(K) = 1$ if unit $i$ is in cluster $k$.
\begin{thm} \label{thm:opt1} Let $\mathrm{tr}(\cdot)$ denote the trace operator. $c^\star$ minimizes \eqref{eqn:optimization_v2} 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}
$$
\end{thm}
\begin{proof}[Proof of Theorem \ref{thm:opt1}] See Appendix \ref{app:optimization}.
\end{proof}
Theorem \ref{thm:opt1} formalizes the optimization problem as a trace-optimization program for a fixed number of clusters $K$ and fixed $\xi$. Theorem \ref{thm:opt1} 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$):
\begin{equation}\label{eq:SDP}
\small
\begin{aligned}
\max_{\mathbf{X}(K)} \text{tr}(\mathbf{L}_{\xi} \mathbf{X}(K)), \quad \text{s.t.}\quad \text{diag}(\mathbf{X}(K)) = \mathbf{1}_{n}, \,\, \mathbf{X}(K)\succeq 0, \quad \mathbf{L}_{\xi;} = n \mathbf{L} - \xi \mathbf{1}_{n}\mathbf{1}_{n}^T.
\end{aligned}
\end{equation}
where
Equation \eqref{eq:SDP} defines a sequence of semi-definite optimization programs, each indexed by the number of clusters $K$ and multiplier $\xi$.
The main distinction between Equation \eqref{eqn:optimization_v2}, and Equation \eqref{eq:SDP} 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 \eqref{eq:SDP}, 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 \citep{von2007tutorial}. Finally, we compare solutions for different values of $K$ and choose the clustering with the largest objective.
Equation \eqref{eq:SDP} is a convex relaxation of the problem in Theorem \ref{thm:opt1} 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 \cite{hong2021optimal} for a review.
We can now describe the complete algorithm in Algorithm \ref{alg:1}. 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 \eqref{eqn:bound2} 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.
\begin{algorithm} [!ht] \caption{Causal Clustering}\label{alg:1}
\begin{algorithmic}[1]
\Require Adjacency matrix $\mathbf{A}$, $\underline{K}$, $\bar{K}$ smallest and largest number of clusters,
expected relative magnitude of spillover effects $\bar{\phi}_n^2/\bar{\psi}$, weight on the variance $\lambda$, grid of values for $\xi \in \mathcal{S} \subset \mathbb{R}_+$.
\For{$\xi \in \mathcal{S}$}
\For{$K \in \{\underline{K}, \cdots, \bar{K}\}$}
\begin{algsubstates}
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space Solve Equation \eqref{eq:SDP} and obtain $\hat{\mathbf{X}}(K)$ as the minimized of Equation \eqref{eq:SDP}
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space Retrive the clusters $c_K$ via
$K$-means algorithm on the first $K$ eigenvectors of $\hat{\mathbf{X}}(K)$
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space Compute the objective function corresponding to the chosen clustering in
\eqref{eqn:optimization_v2}
\end{algsubstates}
\EndFor
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space Save clustering with the lowest objective function in \eqref{eqn:optimization_v2}. Call this $\mathcal{C}(\xi)$.
\EndFor
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space For each $\xi \in \mathcal{S}$ compute the objective in Equation \eqref{eqn:bound2} under each $\mathcal{C}(\xi)$.
\stepcounter{algsubstate}
\Statex {\footnotesize\alph{algsubstate}:}\space Define $\mathcal{C}^\star \in \left\{\mathcal{C}(\xi), \xi \in \mathcal{S}\right\}$ the clustering with the lowest objective in \eqref{eqn:bound2}.
\Return $\mathcal{C}^\star$
\end{algorithmic}
\end{algorithm}
\begin{rem}[Spectral relaxation and comparison to spectral clustering]
Unlike the minimum normalized cut, there is no simple \textit{spectral} relaxation of the optimization problem in Theorem \ref{thm:opt1}, unless all clusters are equally sized. For the special case where all clusters are equal-sized, \eqref{eq:SDP} 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} \label{eqn:UU}
\small
\begin{aligned}
\text{tr}(\mathbf{L}_{\xi} \mathbf{U}(K) \mathbf{U}^T(K)) = \text{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 \citep{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 \citep{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 \textit{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
\end{rem}
\begin{rem}[Reducing the computational complexity with a surrogate relaxation] \label{rem:surrogate} To reduce the complexity of Algorithm \ref{alg:1} it is possible to consider a \textit{single} value of $\xi$ (instead of a grid of values as described in Algorithm \ref{alg:1}) 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
\end{rem}
\vspace{-5mm}
\section{Empirical illustration and numerical studies}\label{sec:empirics}
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 \cite{cai2015social}. We provide further evidence supporting our theoretical analysis using simulated networks in Appendix \ref{sec:numerics}.
\vspace{-3mm}
\subsection{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 \cite{karrer2021network}, outside the scope of this paper.} In Section \ref{sec:info} we turn to the implementation of Causal Clustering.
\medskip
\noindent \textbf{Bias and variance computations} We consider two clustering algorithms implemented at scale on social networks owned by Meta: Louvain algorithm \citep{blondel2008fast}, and Balanced Partitioning \citep{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{tab:meta1}.
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.
$
\medskip \noindent \textbf{Results}
Table \ref{tab:meta1} 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{fig:meta1}. Figure \ref{fig:meta1} 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.
\medskip
\noindent \textbf{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{thm:thumb}, and the dense Messaging graph, Theorem \ref{thm:thumb} 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 \cite{karrer2021network} based on AA-tests.
\definecolor{glaucous}{rgb}{0.38,0.51,0.71}
\begin{figure}[!ht]
\centering
\includegraphics[scale = 0.6]{./figures/plot_meta1}
\includegraphics[scale = 0.6]{./figures/plot_meta2}
\caption{Clusters comparisons for Louvain clustering. Different Types correspond to different numbers of clusters (with Type 1 having the largest number of clusters). Different panels correspond to different graphs where two individuals are not connected if the connection (measured with a continuous variable) is below the $5^{th}, 10^{th}, 50^{th}$ percentile (dense, moderate, and sparse graph). The two graphs in the panels are Facebook friendship and Facebook messaging. }
\label{fig:meta1}
\end{figure}
\vspace{-3mm}
\subsection{Clustering in the field} \label{sec:info}
\begin{figure}[!ht]
\centering
\includegraphics[scale = 0.5]{./figures/cai_app}
\caption{Example of the number of clusters as a function of $\xi$ (left-panel, with dotted line corresponding to 3.2) and objective function in Theorem \ref{thm:opt1} for different clustering. Algorithm \ref{alg:1} corresponds to causal clustering. Data from \cite{cai2015social}, where we report the average result across $47$ regions in the dataset. Here log indicates the natural logarithmic function.}
\label{fig:1}
\end{figure}
Next, we present an application of a field experiment where the treatment consists of informing individuals to increase insurance take-up \citep{cai2015social}.
\cite{cai2015social} collected network information in 185 villages in rural China from 48 larger regions. We use network data collected by \cite{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 \textit{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{alg:1}.
We compare the procedure (denoted as ``Causal Clustering") with the following alternative clustering algorithms: $\varepsilon$-net clustering as in \cite{eckles2017design} with $\varepsilon = 3$ as suggested in \cite{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. }
\begin{figure}[!ht]
\centering
\includegraphics[scale = 0.45]{./figures/mse_2}
\caption{Mean-squared error in calibrated simulations (in natural logs) through comparisons of $M(\lambda)$ in Equation \eqref{eqn:M_lambda} using a surrogate objective with $\xi = \lambda^{-1}$. Simulations are calibrated the model to data from \cite{cai2015social}, averaged over 47 regions. The first three plots vary the variance of the residuals in the outcome model $\sigma^2 \in \{1/4, 1/2, 1\}$, and calibrating the remaining parameters to the model in \cite{cai2015social}, Table 2, Column 4, where outcomes are functions of neighbors' treatments. The last plot calibrate the simulations to settings where the outcomes are functions of the neighbors' \textit{outcome} violating Assumption \ref{ass:first_order}, and $\sigma^2 = 1/4$. }
\label{fig:cai_mse}
\end{figure}
\subsubsection{Main results}
Next we present results for the Pareto frontier and weighted mean-squared error.
\medskip \noindent \textbf{Pareto Frontier} We contrast the objective function with alternative algorithms on the right-hand side of Figure \ref{fig:1}. The proposed method achieves the smallest objective function across all competitors over a wide range of $\xi$, providing suggestive evidence that Algorithm \ref{alg:1} (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{fig:1} (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 \cite{cai2015social} and $\bar{\phi}_n= 0.27$ as in Table 2 in \cite{cai2015social}. Following Remark \ref{rem:surrogate}, a \textit{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.
\medskip \noindent \textbf{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 \cite{cai2015social} (Table 2, Column 4), which assumes that individual outcomes depend on neighbors' treatments, consistently with our Assumption \ref{ass:first_order}. 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
\begin{equation} \label{eqn:M_lambda}
\small
\begin{aligned}
M(\lambda) = \lambda \Big(\mathbb{E}[\hat{\tau}] - \tau\Big)^2 + \mathbb{V}(\hat{\tau})
\end{aligned}
\end{equation}
as the sum of the squared bias and variance, as a function of the parameter $\lambda$. We compare the performance of Algorithm \ref{alg:1} implemented to minimize $L(\xi)$ in Equation \eqref{eqn:optimization_v2}, with $\xi = \lambda^{-1}$ as described in Remark \ref{rem:surrogate}. 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 \eqref{eqn:M_lambda} 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{fig:cai_mse} (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{fig:cai_mse}, we consider settings where the local interference assumption is violated, and we calibrate our model to the specification in \cite{cai2015social} (Column 4, Table 5), where the individual outcome depends on the neighbors' outcomes. Although this is not the main specification in \cite{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.
\vspace{-5mm}
\section{Extensions} \label{sec:extensions}
Here we present a brief description of our extensions, and defer to Appendix \ref{app:more_extensions} details.
\subsection{Experimental design for bias-aware inference} \label{sec:bias_aware_inferece}
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
\begin{equation}\label{eq:exante_CBA_CI}
\small
\begin{aligned}
\mathcal{I}_{\mathrm{ea}}(\beta, \sigma) = [\hat{\tau} \pm \chi(b, \sigma))],
\end{aligned}
\end{equation}
where
\begin{equation}
\small
\begin{aligned}
\label{eq:chi}
\chi(b, \sigma) = (1-\alpha)\text{-th quantile of }|b + \sigma Z|\text{ where }Z \sim N(0, 1).
\end{aligned}
\end{equation}
In addition, define
$$
\small
\begin{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.
\end{aligned}
$$
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{sec:when}.
\begin{lem} \label{thm:inference1} Suppose that Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{ass:clusters} 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.
$
\end{lem}
\begin{proof} See Appendix \ref{proof:lem:inference1}.
\end{proof}
Lemma \ref{thm:inference1} 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{thm:inference1} 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{prop:main}, under a local asymptotic framework with $\varepsilon = o(1)$ in Assumption \ref{defn:local} (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$.
\begin{thm} \label{thm:worst_case_inference} Suppose that Assumptions \ref{ass:first_order}, \ref{ass:exposure_restriction}, \ref{defn:local}, \ref{ass:clusters} 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}
$$
\end{thm}
\begin{proof} See Appendix \ref{proof:thm:worst_case_inference}.
\end{proof}
Theorem \ref{thm:worst_case_inference} shows that the length of the confidence intervals using the worst-case approximations for the bias and variance derived in Section \ref{sec:when} 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{alg:1} (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{alg:1}, 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{thm:inference1} for inference. Such estimator can be constructed using, for instance, the HAC estimator in \cite{gao2023causal}. We provide details in Appendix \ref{sec:estimated_variance2}.
\subsection{Additional extensions}
In Appendix \ref{app:more_extensions} 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.
\vspace{-5mm}
\section{Conclusions and recommendations for practice} \label{rem:psi}
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 \cite{cai2015social}. Our framework directly extends to settings where researchers only have some prior distribution on the network, as in \cite{breza2020using}, see Remark \ref{rem:unobserved_a}.
Similarly to standard power analysis \citep[e.g.][]{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{rem:reg_adj} and consider an estimator as in Equation \eqref{eqn:reg_adj}. 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 \citep{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 \cite{bramoulle2009identification, leung2022causal} among others. Topics of future research also include clustering in dynamic settings, or optimal clustering for a larger class of estimators.
\vspace{-5mm}
\spacingset{1.6}
\bibliography{reference}
\bibliographystyle{chicago}