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.
152,097 characters
Simple Standard Errors for Complicated (Quasi-)Experimental Designs
\maketitle
\begin{abstract}
We develop standard errors for estimates of causal effects in settings where the distribution of treatment exposures is known but potentially exhibits complicated dependence patterns across units. Such dependence arises, for example, in settings with spillovers/interference across units or when the treatment assignment is spatially correlated. Our variance estimator generalizes several properties of the Neyman variance estimator for completely randomized experiments to general designs: it is unbiased for an upper bound on the variance from the design-based perspective; the conservativeness of this upper bound depends on the heterogeneity of treatment effects; and it is also valid under certain forms of super-population sampling of the potential outcomes. Our estimator reduces to standard approaches in simple designs, such as the Neyman variance estimator in completely randomized experiments and the cluster-robust variance estimator in cluster-randomized experiments. We show that the proposed estimator is admissible among the class of quadratic-form variance estimators, and derive asymptotic conditions under which it is consistent for an upper bound on the asymptotic variance of the treatment effect estimator. In simulations and applications based on two empirical studies, one with peer effects and another with a spatially correlated treatment, we find that the proposed variance estimator performs well, and in some cases yields standard errors that are 2-3 times shorter than those of existing alternatives.
\end{abstract}
\section{Introduction}
Social scientists often study settings where the distribution of treatment assignments is assumed to be known, either because it was explicitly randomized or because it was determined by some quasi-random process resembling randomization. In such cases it is natural to leverage knowledge of the assignment scheme to obtain valid inference on causal effects without imposing strong assumptions about how the potential outcomes were generated. This motivates so-called \emph{design-based} inference, which conditions on the potential outcomes and leverages only the randomness of treatment assignment, and therefore is valid regardless of how the potential outcomes were generated. One approach to design-based inference is to conduct Fisher randomization tests (a.k.a. permutation tests), which are valid under arbitrary treatment assignment schemes for testing sharp null hypotheses (e.g. no effects of treatment for any unit). Yet, researchers often wish to obtain standard errors that allow them to test hypotheses and form confidence intervals for (possibly-weighted) average treatment effects. A challenge in the design-based setting is that the variance of the treatment effect estimator typically depends on products of potential outcomes that are not simultaneously observable, and therefore it is not possible to unbiasedly (or consistently) estimate its variance. It is often possible, however, to obtain valid standard errors that estimate an upper bound on the variance of the estimator.
For simple treatment assignment schemes, there are well-accepted approaches for obtaining valid design-based standard errors. In a completely randomized experiment, for example, the most common approach to obtain standard errors for the difference-in-means estimator is to use the \citet{neyman_application_1923} variance estimator, $\frac{1}{N_1} \hat{s}_1^2 + \frac{1}{N_0} \hat{s}_0^2$, where $N_d$ denotes the number of observations in arm $d$ and $\hat{s}_d^2$ is the sample variance of the outcome in arm $d$. This is equivalent to the heteroskedasticity-robust (HC2) variance estimator applied to a regression of the outcome on a treatment dummy. The Neyman variance estimator has three well-known properties (see, e.g., \citet{imbens_causal_2015}, for a textbook treatment). First, from the design-based perspective, it is unbiased for an upper bound on the estimator’s variance. Second, the conservativeness of the upper bound is proportional to the variance of treatment effects, and hence the Neyman variance is exact under treatment effect homogeneity. And third, the Neyman variance corresponds to the standard variance estimator if, in contrast to the design-based paradigm, one imagines the units as having been sampled iid from an infinite-superpopulation. The Neyman variance is thus valid from both the design-based and super-population sampling perspectives. This approach to inference extends naturally to a cluster-randomized experiment, where it is common to use the Neyman variance applied to cluster-level sums (or equivalently, cluster-robust standard errors). The resulting variance estimator is valid from the design-based perspective, with its conservativeness depending on the heterogeneity of treatment effects across clusters; it is also valid from the sampling perspective if we imagine the clusters as drawn independently from a super-population of clusters.
However, in many social-science applications, treatment assignments are correlated across units in more complex ways than in these simple designs. One example where treatment assignment exhibits complex dependence patterns is when treatments occur at a point in space. For example, the researcher may be interested in the effect of having a nearby grocery store on foot traffic to restaurants \citep{pollmann_causal_2026}, or the impact of having a nearby factory on economic outcomes \citep{dell_development_2020}. In such settings, nearby locations will tend to have correlated treatment assignments. A second example is settings in which there are spillovers across units. For example, assigning one student to an anti-bullying intervention as in \citet{paluck_changing_2016} may have impacts on that student's friends. In such settings, it is common to specify an exposure mapping that formalizes how one unit's outcome depends on other units' treatments \citep{aronow_estimating_2017}; e.g. we may specify that $i$'s outcome depends both on whether $i$ is treated and on whether any of $i$'s friends are treated. Even if the original treatment is assigned independently across units, the distribution of treatment \emph{exposures} will then exhibit dependence among units who are connected in the network. Correlated treatment assignments arise in a variety of other settings, such as when the treatment is determined by a centralized assignment mechanism (e.g. assignment of students to schools), when experiments are conducted in two-sided markets, or when an experiment uses a complex assignment mechanism involving re-randomization, stratification, or sequential treatment assignment. More generally, researchers frequently conduct permutation tests that simulate counterfactual treatment assignments by permuting observed treatments along a particular dimension, inducing a known but dependent assignment scheme.
This raises the question: can we construct a variance estimator for settings with complicated assignment schemes that has similar properties to the Neyman variance estimator in completely randomized experiments? In this paper, we propose a generalization of the Neyman variance estimator for completely randomized experiments that is applicable under complex assignment schemes. The estimator reduces to the Neyman variance estimator when applied to a completely randomized experiment and to a cluster-robust variance estimator when applied to a cluster-randomized experiment. Moreover, the estimator possesses similar properties to the Neyman variance estimator under more complex designs: it is valid from the design-based perspective, its conservativeness depends on treatment-effect heterogeneity, and it is also valid under certain superpopulation sampling schemes.
We now describe the setting and proposed variance estimator in more detail. At a high level, we suppose that there is a finite population of $K$ units each with an attribute $Y_k$, and we are interested in estimating the average $\tau = \frac{1}{K} \sum_k Y_k$ from a sample wherein only the attributes for individuals with $W_k = 1$ are observed. As described in \Cref{sec:setup}, this deceptively simple setup nests a large number of causal settings by setting $k$ to index the relevant potential outcomes, setting $Y_k$ to be an appropriate transformation of the $k$th potential outcome, and setting $W_k$ to be an indicator for whether the relevant potential outcome is observed. For example, if the treatment exposure $D_i$ can take values $\{0,\ldots, Q\}$, then we can let $k=(i,d)$ index the $K=N(Q+1)$ combinations of $(i,d)$. Setting $Y_k = \omega_{i,d} Y_i(d)$ for some chosen weights $\omega_{i,d}$ then yields the causal contrast $\tau = \frac{1}{K} \sum_i \sum_d \omega_{i,d} Y_i(d)$. We then observe $Y_k$ if and only if $W_k = 1\{D_i =d\}$ is equal to 1. The framework thus nests estimation of average causal contrasts in settings where the treatment exposure takes on a finite number of values. Beyond causal settings, this set-up also nests the problem of learning about the mean of an attribute in a finite population under a complicated survey sampling scheme.
Within this framework, we focus on the problem of estimating the variance of the Horvitz-Thompson (inverse-probability weighting) estimator, $\hat\tau = \frac{1}{K} \sum_k \tfrac{W_k}{\pi_k} Y_k$, where $\pi_k = P(W_k =1)$. (As we describe below, a similar approach extends to other estimators that have an asymptotic linear representation, such as those obtained by linear regressions.) We adopt a design-based approach wherein the $Y_k$ are viewed as fixed and the variation in the data comes from the stochastic nature of the assignment $W_k$. We refer to the distribution of $W_k$ as the \emph{assignment} or \emph{design} distribution, and we assume that it is known, so that $\pi_k$ and $\pi_{kk'} = P(W_k=1, W_{k'}=1)$ are both known. As first shown in \citet{horvitz_generalization_1952}, the variance of $\hat\tau$ over the design distribution is given by $\operatorname{Var}(\hat\tau) = K^{-2} \sum_{k,k'} v_{kk'} Y_kY_{k'}$, where $v_{kk'} = \operatorname{Cov}(\tfrac{W_k}{\pi_k}, \tfrac{W_{k'}}{\pi_{k'}} )$ is a function of the distribution of $W_k$ and hence is known. If $Y_k$ and $Y_{k'}$ are both observed with positive probability, so that $\pi_{kk'}= P(W_k=1,W_{k'}=1)>0$, then it is straightforward to unbiasedly estimate $Y_kY_{k'}$ by inverse probability weighting, $\frac{W_k W_{k'}}{\pi_{kk'}} Y_k Y_{k'}$. The challenge for variance estimation is that certain pairs of potential outcomes are not jointly observable, so that $\pi_{kk'}=0$ for some pairs $k,k'$. For example, if $k$ and $k'$ correspond to different potential outcomes for the same unit, say $Y_i(0)$ and $Y_i(1)$, then it is impossible to observe both $W_k=1$ and $W_{k'}=1$ simultaneously. Under complex assignment schemes, other pairs of potential outcomes may also be jointly unobservable. For example, if the treatment has spillovers and units $i$ and $j$ are neighbors, then we cannot observe that unit $i$ is treated while unit $j$ has no treated neighbors. To derive valid standard errors, we therefore need to obtain some upper bound on the terms involving products $Y_k Y_{k'}$ for pairs of potential outcomes that are jointly unobservable.
Our approach to constructing an upper bound on the variance is based on what we call the \emph{joint unobservability graph}. This is defined to be the graph among the nodes $k=1,\ldots,K$, where $k$ and $k'$ are adjacent in the graph if $Y_{k}$ and $Y_{k'}$ cannot be jointly observed, i.e. $\pi_{kk'}=0$. We can then partition the joint unobservability graph into its connected components $\mathcal{K}_1,\ldots,\mathcal{K}_C$. This construction implies that if $\pi_{kk'}=0$, then $k,k'$ are in the same $\mathcal{K}_c$ for some $c$. For our main approach, we assume that there are at least two connected components, $C \geq 2$, although we consider an extension of the main idea to settings with $C=1$ in \Cref{sec:arbitrary-clusters}. In a completely randomized experiment, each connected component contains the two potential outcomes for a given unit $i$. Likewise, in a cluster-randomized experiment, each connected component collects all the potential outcomes for units in a particular cluster. We thus see that under these two simple assignment schemes, the standard approaches to variance estimation have conservativeness that is proportional to the variance of treatment effect sums across the connected components.
Our proposed estimator generalizes this property to more general designs: its conservativeness is proportional to the variance of treatment effect sums across connected components. Specifically, we consider the variance estimator
$$\widehat{V}_{\mathrm{CC}} = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} (v_{kk'} + s_{kk'}) \frac{W_k W_{k'}}{ \pi_{kk'} } Y_k Y_{k'}, \text{ where } s_{kk'} = \begin{cases}1 & k,k' \in \mathcal{K}_c \text{ for some } c \\
-1/(C-1) & \text{o.w.}
\end{cases} $$
which we refer to as the connected components (CC) estimator. The estimator $\widehat{V}_{\mathrm{CC}}$ does inverse-probability weighting only over pairs $k,k'$ that are jointly observable. However, rather than using the weights $v_{kk'}$ found in the expression for $\operatorname{Var}(\hat\tau)$, it uses $v_{kk'}+s_{kk'}$, where $s_{kk'}$ is 1 if individuals are in the same connected component and $-1/(C-1)$ otherwise. The CC estimator thus compensates for the jointly unobservable pairs by adding the slack term $s_{kk'}$, which depends on whether $k$ and $k'$ are in the same component, to the observable pairs. We show that $\widehat{V}_{\mathrm{CC}}$ matches the Neyman variance in a completely randomized experiment, and a cluster-robust variance estimator in a cluster-randomized completely randomized experiment. More generally, our first main result (\Cref{prop:cc-valid}) shows that $\widehat{V}_{\mathrm{CC}}$ is weakly upwards-biased for $\operatorname{Var}(\hat\tau)$, with conservativeness proportional to the variance of treatment effects across clusters, $\operatorname{Var}_c(\tau_c)$, where $\tau_c = \sum_{k \in \mathcal{K}_c} Y_k$. It follows that $\widehat{V}_{\mathrm{CC}}$ is unbiased if treatment effect sums are homogeneous across connected components (which necessarily occurs under the sharp null when all effects are zero). In addition to its design-based interpretation, we further show that $\widehat{V}_{\mathrm{CC}}$ can be viewed as an unbiased estimator of $\operatorname{Var}(\hat\tau)$ under a sampling scheme in which the component sums $\tau_c$ are sampled iid across components. The CC estimator thus possesses the three main properties of the Neyman variance estimator for completely randomized experiments in more general designs.
We further establish that $\widehat{V}_{\mathrm{CC}}$ is admissible in the class of quadratic-form variance estimators (defined formally below), i.e. there is no quadratic-form variance estimator that is upwards-biased for $\operatorname{Var}(\hat\tau)$ and that dominates $\widehat{V}_{\mathrm{CC}}$ in the sense that it is always weakly less conservative and strictly so for some choice of potential outcomes. Our proof of admissibility builds on the results of \citet{harshaw_optimized_2026}, who show that any admissible quadratic-form estimator corresponds to the solution of a class of optimization problems with an objective over matrices that is increasing in the Loewner order. We show that $\widehat{V}_{\mathrm{CC}}$ is a solution to this optimization for a particular choice of objective, and hence is admissible. The objective that $\widehat{V}_{\mathrm{CC}}$ optimizes also has a Bayesian interpretation. Specifically, consider a Bayesian who has a prior over the potential outcomes and wishes to minimize the expected conservativeness of the variance estimator, subject to the constraint that it is weakly upwards-biased for all possible potential outcomes. We show that $\widehat{V}_{\mathrm{CC}}$ corresponds to the optimal solution (among quadratic-form estimators) for a Bayesian whose prior suggests that cross-component heterogeneity in treatment effects is small and ``equally spread out'' across components in a sense that we make precise below. Since $\widehat{V}_{\mathrm{CC}}$ reduces to the Neyman variance estimator in a completely randomized experiment and the cluster-robust variance estimator in cluster-randomized experiments, this result also gives a novel Bayesian justification for these commonly-used variance estimators.
In \Cref{sec:comparison-in-experiments}, we compare the CC estimator to other design-based estimators that have been proposed for general assignment schemes, focusing primarily on the estimators proposed by \citet{aronow_estimating_2017} and \citet{harshaw_optimized_2026}. When specialized to a completely randomized experiment, the CC estimator reduces to the Neyman variance, with conservativeness proportional to the variance of treatment effects, $ \frac{1}{N}\operatorname{Var}_i(\tau_i)$. By contrast, we show that the \citet{aronow_estimating_2017} estimator has conservativeness proportional to the uncentered second moment of treatment effects, $\frac{1}{N} E_i(\tau_i^2)$. Hence, we expect the CC estimator to be much less conservative when treatment effects are large but relatively homogeneous, so that $\operatorname{Var}_i(\tau_i) \ll E_i(\tau_i^2)$.\footnote{We note that $\operatorname{Var}_i(\tau_i)$ is defined to be the finite-population variance, with denominator $N-1$ rather than $N$, so it is not the case that $\operatorname{Var}_i(\tau_i) \leq E_i(\tau_i^2)$ for all choices of potential outcomes. Hence, the conservativeness of the \citet{aronow_estimating_2017} estimator can be slightly smaller than that of the CC estimator, but only by the factor $(N-1)/N$.} \citet{harshaw_optimized_2026} propose a class of estimators based on minimizing different measures of the slackness of the upper bound, and in their simulations consider implementations using simple norms such as the trace and Frobenius norms. In a completely randomized experiment with equal-sized treatment arms, we show that the \citet{harshaw_optimized_2026} estimator using the Frobenius norm objective corresponds exactly to the AS estimator, and hence the same comparison applies to this estimator. Turning to more general designs, any variance estimator corresponding to a solution to the class of optimizations considered in \citet{harshaw_optimized_2026} is admissible in the class of quadratic-form variance estimators, as is the CC estimator, and thus neither dominates the other for all choices of potential outcomes. The CC estimator is exact when treatment effects are homogeneous across components, while this is not guaranteed for the \citet{harshaw_optimized_2026} estimator based on the trace or Frobenius norms, and thus we expect the CC estimator to be less conservative when effects are relatively homogeneous across components, while admissibility implies that the trace and Frobenius-norm based estimators must be less conservative under at least some forms of treatment effect heterogeneity. We further show in \Cref{subsec:complex-designs} that all valid quadratic-form variance estimators that are exact under the sharp null of no treatment effects have conservativeness that depends on the heterogeneity of treatment effects across connected components. Specifically, the conservativeness of any such estimator is proportional to a quadratic form involving the component-level sums of treatment effects and a correlation matrix. The CC estimator corresponds to a particular choice of the correlation matrix, which we show has the minimum operator norm among correlation matrices corresponding to variance estimators that are exact under component homogeneity.\footnote{Estimators based on other correlation matrices, which optimize different objectives, are also possible; see \Cref{rem:alt-estimators}.}
In \Cref{sec:consistency}, we consider finite-population asymptotics and provide conditions under which the CC estimator is consistent for an upper bound on the asymptotic variance of $\hat\tau$, employing similar arguments for the consistency of Horvitz-Thompson estimators as in \citet{delevoye_consistency_2020}. In \Cref{sec:non-ht-estimators}, we show that our approach to variance estimation for Horvitz-Thompson estimators is also asymptotically valid for other estimators that have an asymptotic linear representation. We provide high-level conditions for consistency for generic asymptotically linear estimators, and lower-level conditions for estimators based on linear regressions.
In \Cref{sec:arbitrary-clusters}, we consider the setting where the joint unobservability graph is connected ($C=1$), in which case the CC estimator is not well-defined. However, we show that if the units in the population are partitioned into clusters, then using similar ideas it is possible to construct a variance estimator whose conservativeness intuitively depends on (a) the heterogeneity of treatment effects across the clusters, and (b) the number of cross-cluster pairs $k,k'$ that have $\pi_{kk'}=0$. The resulting estimator will thus not be too conservative if the clusters are relatively homogeneous and have relatively few joint unobservability links between them. Owing to term (b), the resulting estimator will not be exact even under the sharp null of no treatment effects. We show, however, that some conservativeness under the sharp null appears to be necessary in the challenging case when $C=1$: the only quadratic-form variance estimator that is exact under the sharp null is the trivial estimator $\hat\tau^2$, which always produces a $t$-statistic with absolute value of 1.
We find that the CC estimator performs well in calibrated simulations and applications to two applied settings, as described in \Cref{sec:sims-and-apps}. We first revisit and extend \citet{harshaw_optimized_2026}'s simulation study calibrated to \citet{paluck_changing_2016}, who studied an anti-bullying intervention in schools which is expected to have social spillovers. We find that confidence intervals based on the CC estimator have good coverage and are shorter on average than the \citet{aronow_estimating_2017} estimator and two implementations of the \citet{harshaw_optimized_2026} approach based on simple norms in all of the simulation specifications taken from \citet{harshaw_optimized_2026}. In line with the intuition from a completely randomized experiment, the differences tend to be larger when the treatment effect is large but relatively homogeneous. Applied to the actual \citet{paluck_changing_2016} data, the CC estimator yields the shortest standard error in 7 of 8 specifications; in the most extreme case, the CC SE is less than half that of the other estimators. The CC estimator can also be computed in seconds in all specifications, compared with over 30 hours in some cases for the \citet{harshaw_optimized_2026} approach, which relies on nonlinear optimization. Our second simulation design is based on the spatial treatment setting studied in \citet{pollmann_causal_2026}, who investigated the impact of having a nearby grocery store on foot traffic to restaurants. We calibrate our simulations based on a public-use synthetic dataset provided by \citet{pollmann_causal_2026}, which was simulated to match features of the proprietary data used in his actual application. We find in our simulations that the CC variance estimator yields confidence intervals with close to nominal coverage, and that are approximately three times shorter on average than those obtained using the variance estimator proposed by \citet{pollmann_causal_2026}. We likewise obtain standard errors that are approximately three times smaller when applied to Pollmann's dataset. We also find in our simulations that the CC estimator is competitive with spatial HAC \citep{conley_gmm_1999} using a design-motivated choice of kernel and bandwidth \citep{wang_design-based_2025} under homogeneous effects and simple forms of treatment effect heterogeneity, while remaining valid under non-homophilic forms of heterogeneity that cause HAC to be downward-biased. We emphasize that all these comparisons across estimators may be specific to the simulation designs considered, as theory generally suggests that there is no uniform dominance relationship across the estimators considered.
\paragraph{Related Literature.} The most closely related work is \citet{aronow_estimating_2017} and \citet{harshaw_optimized_2026}, discussed above, who also propose variance estimation methods for causal settings with known but potentially arbitrary assignment probabilities and a finite number of known treatment exposure levels. \citet{chang_design-based_2025} also studies complex experimental designs; his focus is mainly on estimation, although he also extends the \citet{aronow_estimating_2017} approach to variance estimation to a broader class of estimators. Our contribution relative to this work is to propose a variance estimator that generalizes the well-known properties of the Neyman variance estimator to more complicated designs. Another closely related paper is \citet{pollmann_causal_2026}, who considers a design-based setting with spatial treatments that, as we describe in \Cref{sec:sims-and-apps} below, can be cast as a special case of the setting we consider.
Our work relates to the broader literatures on inference in settings with spatial treatments or network spillovers. One common approach is to use HAC-style variance estimators, such as \citet{conley_gmm_1999}. HAC-style estimators were originally developed in a model-based framework in which one restricts the dependence of errors over space; a challenge in implementing these estimators in practice is choosing the distance metric, bandwidth, and kernel, which together govern the allowed spatial correlation. \citet{xu_design-based_2026} study HAC-style estimators from the design-based perspective and provide general mixing conditions under which spatial HAC is asymptotically valid, but choosing a bandwidth and kernel for good finite-sample performance remains a practical challenge. Other recent work has provided a design-based justification for HAC-style variance estimators with particular choices of distance metric, bandwidth, and kernel in specific network \citep{leung_causal_2022} and spatial-treatment \citep{wang_design-based_2025} settings.\footnote{The appendices to \citet{harshaw_optimized_2026} and \citet{pollmann_causal_2026} respectively discuss how the frameworks of \citet{leung_causal_2022} and \citet{wang_design-based_2025} relate to our design-based framework. \citet{harshaw_optimized_2026}, for example, find that the \citet{leung_causal_2022} variance estimator can be anti-conservative in their simulation design, which we adopt in \Cref{sec:sims-and-apps} below. They also discuss how the \citet{leung_causal_2022} estimator may be preferable in other designs where the exposure mapping is mis-specified.} One advantage of these approaches relative to ours is that they typically can allow for some mis-specification of the exposure mapping (e.g. approximate network interference), whereas our framework relies on a correctly-specified exposure mapping; see \citet{savje_average_2021} for discussion of how inference can be distorted under mis-specified exposure mappings. On the other hand, these results typically rely on specific treatment assignment mechanisms (e.g. Bernoulli assignment of nodes or candidate sites) and require some restrictions on treatment effects, either in terms of how they decay with distance in the network, or that treatment effects satisfy a spatial homophily assumption. The CC estimator, by contrast, can accommodate complex treatment assignment mechanisms and does not require restrictions on treatment effect heterogeneity or the choice of a distance metric, bandwidth or kernel, although it does require that $C>1$ and may be conservative when treatment effects are very heterogeneous. We thus view the two approaches as complementary in spatial settings since they have non-nested use-cases; see \Cref{subsec:complex-designs} for further discussion. An alternative design-based approach in settings with spillovers is to use randomization tests, as in \citet{athey_exact_2018}, \citet{basse_randomization_2019}, \citet{zhong_unconditional_2024}, which yield finite-sample valid inference for (partially) sharp nulls in settings with interference; by contrast, we focus on developing standard errors for weak nulls on average effects. See \citet{bhadra_causal_2025} for a review of methods for inference in the presence of spillovers, including model-based approaches that make assumptions about the formation of the network and potential outcomes.
\section{Set-up\label{sec:setup}}
For our formal results, we will consider the problem of estimating the mean of an attribute in a finite population under a general sampling scheme; as we will show shortly, this framework nests a variety of causal settings if each item in the finite population corresponds to a potential outcome of interest. Formally, we suppose there exists a finite population of $K$ items, where the $k$th item has attribute $Y_k$. These attributes are viewed as non-stochastic (or conditioned on). We are interested in the population average $\tau := \tfrac{1}{K} \sum_k Y_k$. We denote by $W_k$ the random variable indicating whether the $k$th item is observed. We denote the first-order and second-order sampling probabilities by $\pi_k := P(W_k =1)$ and $\pi_{kk'} := P(W_k = 1, W_{k'}=1)$, respectively, and assume that these probabilities are known (or simulatable to arbitrary accuracy). We assume that $\pi_k>0$ for all $k$, so that all attributes are sampled with positive probability.\footnote{When $k$ indexes potential outcomes, the requirement that $\pi_k>0$ for all $k$ imposes the overlap restriction that each unit receives each treatment level with positive probability. Since the design is assumed to be known, the researcher can restrict attention to the sub-population for which this overlap condition holds. For ease of exposition, we assume that the researcher has already restricted attention to the population satisfying overlap.} We consider the Horvitz-Thompson estimator $\hat\tau = \frac{1}{K} \sum_k \frac{W_k}{\pi_k} Y_k$; see \Cref{sec:non-ht-estimators} below for extensions to other estimators.
We can nest a variety of causal settings in this framework by having $k$ index the relevant potential outcomes and setting $Y_k$ to be an appropriate transformation of the potential outcomes. Consider a setting where unit $i$ receives treatment exposure level $D_i \in \{0,\ldots, Q\}$, with realized outcome $Y_i = Y_i(D_i)$, where $Y_i(\cdot)$ is $i$'s potential outcome, which following the design-based paradigm we view as fixed (or conditioned on). If we let $k=(i,d)$ index the $K=N (Q+1)$ combinations of $i$ and $d$, and set $Y_k = \omega_{i,d}Y_i(d)$ for some weights $\omega_{i,d}$, then $\tau$ corresponds to the causal contrast $\frac{1}{K} \sum_{i} \sum_d \omega_{i,d} Y_i(d)$, and $\hat{\tau}$ corresponds to the Horvitz-Thompson estimator $\frac{1}{K} \sum_i \sum_d \frac{1\{ D_i =d \} }{ P(D_i =d)} \omega_{i,d} Y_i$. To be more concrete, we now outline several examples fitting this framework.
\begin{example}[Completely randomized experiment] \label{example:cre}
Consider a completely randomized experiment (CRE) in which $N_1$ out of $N$ units are randomly selected to receive a binary treatment ($D_i \in \{0,1\}$) with equal probability. (We will assume throughout that $2 \leq N_1 \leq N-2$, so there are at least two units in each treatment arm.) Let $k$ index the possible values of $(i,d)$ for $i=1,\ldots,N$ and $d = 0,1$, so that $K= 2N$. For $k=(i,d)$, define $Y_k = 2Y_i(1)$ if $d=1$ and $-2 Y_i(0)$ if $d=0$, and let $W_k = 1\{D_i =d \}$. Then $$\tau = \frac{1}{K} \sum_k Y_k = \frac{1}{N} \sum_i Y_i(1) - Y_i(0), $$
the finite-population average treatment effect. For $k=(i,d)$ and $k'=(i',d')$, we have $\pi_k = P(D_i =d)$ and $\pi_{kk'} = P(D_i=d, D_{i'} = d')$. For example, if $d=d' = 1$ and $i \neq i'$, then $\pi_k = P(D_i=1) = \frac{N_1}{N}$ and $\pi_{kk'} = P(D_i = 1, D_{i'}= 1) = \frac{N_1}{N} \frac{N_1 - 1}{N-1}$.
\end{example}
\begin{example}[Cluster randomized experiment] \label{example:cluster-experiment}
As in the previous example, suppose there is a binary treatment. However, we now suppose that the $N$ units are divided into $G$ different clusters (or groups), and denote by $g(i) \in \{1,\ldots,G\}$ the cluster of unit $i$. For simplicity, we assume each cluster $g$ has exactly $m=N/G$ units. We suppose that we have a cluster randomized experiment in which $G_1$ out of the $G$ clusters are selected for treatment, and unit $i$ is treated if and only if $D_{g(i)}=1$, where $D_g$ is an indicator for whether cluster $g$ is selected for treatment. (We will assume throughout that $2 \leq G_1 \leq G-2$, so there are at least 2 clusters in each treatment arm.) We can define $k$, $Y_k$, and $W_k$ as in the previous example, so that the target parameter is still the ATE. The only change is that the $\pi_{k}, \pi_{kk'}$ now correspond to the treatment probabilities under \emph{clustered} assignment. For example, if $k = (i,1)$ and $k' = (i',1)$ for $i$ and $i'$ in the same cluster ($g(i)= g(i')$), then $\pi_{kk'}= P(D_{g(i)} =1, D_{g(i')}=1) = \frac{G_1}{G}$. Likewise, if $k=(i,1)$ and $k'=(i',0)$ for $i$ and $i'$ again in the same cluster, then $\pi_{kk'}=0$, since units in the same cluster cannot have different treatment statuses.
\end{example}
\begin{example}[Spatial treatments] \label{example:spatial}
Similar to \citet{pollmann_causal_2026} and \citet{wang_design-based_2025}, suppose we have locations $i=1,\ldots,N$ that are located at various points on a map. We are interested in the effect of having a local amenity (or disamenity) nearby; for example, we may wish to know how the presence of a grocery store affects foot traffic, or how a polluting plant affects property values. We suppose that the sites of the amenities are randomly selected from one of $S$ possible sites, and let $T_s \in \{0,1\}$ be an indicator for whether site $s$ is selected. We may then say that location $i$ is treated if there is a closeby selected site, so that $D_i = 1\{ \exists s \text{ s.t. } T_s =1, d(s,i) \leq d_0 \}$, where $d(s,i)$ is the distance between $s$ and $i$ and $d_0$ is some threshold. (We could likewise define treatment to depend on the \emph{number} of closeby sites, or even the set of sites within distance $d_0$ of $i$.) We can again define $k$ and $Y_k$ as in the previous examples so that the target parameter is the ATE. The key change is that the $\pi_{k}$ and $\pi_{kk'}$ are now determined by the spatial nature of the treatment. For $k=(i,1)$ and $k'=(i',1)$, for example, $\pi_{kk'}$ corresponds to the probability that there is both a site close to $i$ that is treated and a site close to $i'$ that is treated.
\end{example}
\begin{example}[Interference] \label{example:interference}
In many settings, we may expect the treatment assignment of one unit to affect the outcome of another unit. For example, if $T_i$ denotes whether individual $i$ is randomized to receive a vaccine, and unit $j$ is a family member of unit $i$, then we may expect $T_i$ to affect the outcome for unit $j$. A popular framework for studying spillovers is the exposure mapping framework of \citet{aronow_estimating_2017}, in which we specify that unit $i$'s treatment depends on some function $h_i(T_i,T_{-i})$. In the vaccine context, we expect that unit $i$'s outcome depends on both whether they were vaccinated and whether any of their neighbors were vaccinated, in which case $h$ takes on four possible values (for each combination of own-treatment and neighbors' treatment). In this case, we can let $k$ index the relevant values for $(i,h)$ and let $Y_k = \omega_{i,h} Y_i(h)$, in which case $\tau$ estimates the contrast $\frac{1}{K} \sum_i \sum_h \omega_{i,h} Y_i(h)$. For example, if we're interested in the $(1,1)$ vs. $(0,0)$ contrast, then setting $\omega_{i,h}=2$ if $h=(1,1)$ and $-2$ if $h=(0,0)$ yields the average total effect, $\tau = \frac{1}{N} \sum_i Y_i(1,1) - Y_i(0,0)$. In this case, the $\pi_k$ and $\pi_{kk'}$ are determined both by the treatment assignment $T_i$ and by the network.
\end{example}
\paragraph{Other examples.} Although we focus on the four examples given above, we briefly sketch several other classes of applications in which the assignment mechanism is known but may exhibit dependence across units. First, a variety of quasi-random processes will induce treatments with spatial correlation, beyond the example with site-level treatments in \Cref{example:spatial}. For example, previous papers have argued that the timing of hurricanes is quasi-random \citep[e.g.][]{Deryugina2017}; the treatment assignments induced by randomly assigning dates to hurricane paths will exhibit complicated dependence owing to multiple places being hit by the same hurricane. Second, consider the literature that has exploited randomness in centralized assignment mechanisms that, say, assign students to schools or doctors to residencies \citep{abdulkadiroglu_research_2017,chen_nonparametric_2026}. It is straightforward to simulate from the assignment mechanism many times, and therefore to learn assignment probabilities to arbitrary accuracy. However, owing to capacity constraints, assignments are likely to exhibit complicated dependence across individuals. Third, consider experiments in two-sided markets---e.g. a randomized intervention to drivers on Uber or sellers on eBay \citep{johari_experimental_2021}. If we are interested in the impact of this intervention on the other side of the market---e.g. the impact on Uber riders or eBay buyers---then the exposures to treatment likely exhibit complicated dependence patterns since, for example, two riders who had the same driver will have correlated exposures. Fourth, in a variety of settings, researchers employ potentially complex treatment assignment protocols involving re-randomization, stratification, or sequential treatment assignment rules \citep[e.g.][]{li_asymptotic_2018,bai_optimality_2022}. Finally, we note that in a variety of settings, researchers run permutation tests that simulate counterfactual treatment exposures by permuting the observed treatment exposures along various dimensions: the assignment of treatments over this distribution is known, but may again exhibit complicated dependence across units.
With these examples in mind, we now describe the challenge of estimating the variance of $\hat\tau$. Since the $Y_k$ are fixed, as shown in \citet{horvitz_generalization_1952},
\begin{equation} \operatorname{Var}(\hat\tau) = \frac{1}{K^2} \sum_{k,k'} v_{kk'} Y_k Y_{k'}, \qquad \text{ where } \qquad v_{kk'} := \underbracket{\frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}}}_{= \operatorname{Cov}\big(\frac{W_k}{\pi_k}, \frac{W_{k'}}{\pi_{k'}} \big)}. \label{eqn:var-tauhat}
\end{equation}
We can decompose this expression into sums involving pairs $(k,k')$ with $\pi_{kk'}>0$ and those with $\pi_{kk'} = 0$,
$$ \operatorname{Var}(\hat\tau) = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} v_{kk'} Y_k Y_{k'} + \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} =0 } } -Y_k Y_{k'}.$$
The first sum is straightforward to estimate since it is a weighted sum of terms of the form $Y_k Y_{k'}$, where $Y_k$ and $Y_{k'}$ are simultaneously observable with positive probability ($\pi_{kk'}>0$). These terms can be unbiasedly estimated by Horvitz-Thompson weighting, $\operatorname{E}[ v_{kk'} \frac{W_k W_{k'} }{ \pi_{kk'} } Y_k Y_{k'} ]= v_{kk'} Y_k Y_{k'}$. The challenge arises from the second term, which involves products $-Y_k Y_{k'}$ which are never directly observed since $\pi_{kk'}=0$. The key challenge is to find an estimable upper bound on the sum of these unobservable terms.
\section{Proposed estimator and its properties\label{sec:estimator}}
To derive an upper bound on the jointly unobservable product terms that enter $\operatorname{Var}(\hat\tau)$, we first formalize the structure of which potential outcomes cannot be observed together using what we call the \emph{joint unobservability graph}. We first define this graph, then show how it can be used to obtain a variance estimator.
\paragraph{Joint Unobservability Graph.} We say that $k$ and $k'$ are jointly unobservable if $\pi_{kk'}=0$. We then consider the \emph{joint unobservability graph} in which the nodes $k$ and $k'$ are connected if and only if $k,k'$ are jointly unobservable. This graph has adjacency matrix $A$ such that $A_{kk'} = 1$ if and only if $\pi_{kk'} =0$. We can then partition the indices $\{1,...,K\}$ into the connected components of the graph: i.e. the finest partition $\mathcal{K}_1,...,\mathcal{K}_C$ of $\{1,...,K\}$ such that if $\pi_{kk'} =0$, then $k,k' \in \mathcal{K}_c$ for some $c$. We will assume for now that $C \geq 2$, so the joint unobservability graph is not connected. See \Cref{sec:arbitrary-clusters} below for an extension to the setting where the joint unobservability graph may be connected but has ``weakly connected'' clusters.
\paragraph{Examples (revisited).} Revisiting our examples from earlier, we note that in a completely randomized experiment, $k$ and $k'$ are jointly unobservable if and only if they correspond to the potential outcomes for the same unit ($Y_i(0),Y_i(1)$). There are thus $N$ connected components of the graph, where each connected component contains the two potential outcomes corresponding to a single unit. In a cluster-randomized experiment, $k$ and $k'$ are jointly unobservable if either (i) they correspond to the two potential outcomes for a single unit $i$, or (ii) they correspond to potential outcomes $Y_i(1)$ and $Y_j(0)$ for $i$ and $j$ in the same cluster. It follows that the number of connected components equals the number of clusters ($G$), with each component collecting all the potential outcomes associated with a particular cluster $g$. Next consider our spatial treatment example (\Cref{example:spatial}). As before, the two potential outcomes for the same unit are jointly unobservable. Moreover, if every possible site $s$ that is within distance $d_0$ of unit $i$ is also within distance $d_0$ of unit $j$, then it is impossible to jointly observe $Y_i(1)$ and $Y_j(0)$. The joint unobservability graph thus contains edges generated from these two types of pairs of potential outcomes. The number of connected components then depends on the spatial orientation of the sites and units. Finally, consider the setting with interference (\Cref{example:interference}). Suppose the exposure mapping $h$ captures own-treatment and whether any neighbors were treated, every unit has at least one neighbor, and we are interested in the $(1,1)$ vs. $(0,0)$ contrast. Suppose the individual-level treatments $T_i$ are assigned by Bernoulli randomization. In this case, it turns out that each connected component of the joint unobservability graph collects all the potential outcomes for units in a given connected component of the original network. To see why, observe that as in our other examples, the two potential outcomes for the same unit $i$ are jointly unobservable. Moreover, $Y_i(1,1)$ and $Y_j(0,0)$ are not jointly observable if $i$ and $j$ are neighbors, since it is not possible for $i$ to be treated while $j$ has no treated neighbors. This implies that all the potential outcomes for units who are connected in the original network must be connected in the joint unobservability graph. On the other hand, it is straightforward to verify that all potential outcomes for units that are not connected in the original graph are jointly observable under Bernoulli assignment. Thus, the connected components of the joint unobservability graph are one-to-one with those of the original network. We note that for other contrasts, the connected components of the joint unobservability graph need not align with the connected components of the original graph; see our application to \citet{paluck_changing_2016} in \Cref{sec:sims-and-apps} for an example.
\begin{rem}[Relationship to dependency graphs]
We note that the joint unobservability graph is distinct from the dependency graph considered in e.g. \citet{baldi_asymptotic_1989} and \citet{leung_treatment_2020}. In a dependency graph, $k$ and $k'$ are adjacent if $W_k$ and $W_{k'}$ are dependent.\footnote{The definition of a dependency graph typically requires that for disjoint sets $S_1,S_2$ with no links between them, $(W_k)_{k \in S_1} \perp \!\!\! \perp (W_k)_{k\in S_2}$, which implies that $k$ and $k'$ are adjacent if they are dependent but is actually a stronger condition.} By contrast, the joint unobservability graph defines $k$ and $k'$ to be adjacent only if $\pi_{kk'} = P(W_k=1,W_{k'}=1)=0$, which is a stricter criterion. To appreciate the difference, note that in a completely randomized experiment, the treatment assignment indicators for all units are correlated owing to the constraint that $\sum_i D_i = N_1$. It follows that the dependency graph is connected. By contrast, as described above, the joint unobservability graph has $N$ connected components, with each component containing only the two potential outcomes for a given unit. This distinction arises because while treatment indicators are correlated across units, any pair of treatment assignments for units $i \neq j$ occurs with positive probability (i.e. $P(D_i = d, D_j = d') > 0$ for any $i \neq j$ and $d,d'\in \{0,1\}$). Thus pairs of potential outcomes involving distinct units are adjacent in the dependency graph but not in the joint unobservability graph.
\end{rem}
\paragraph{The CC Estimator.}
We now derive a variance estimator that is unbiased for an upper bound on $\operatorname{Var}(\hat\tau)$. The conservativeness of the variance estimator is proportional to the heterogeneity across the connected components of the joint unobservability graph.
Let $s_{kk'} = 1$ if $k$ and $k'$ are in the same component ($k,k' \in \mathcal{K}_c$ for some $c$), and let $s_{kk'} = - 1/(C-1)$ otherwise. (We will assume $C \geq 2$ so that division by $C-1$ is well-defined.) Let $b_{kk'} = v_{kk'}+ s_{kk'}$. Define the feasible variance estimator
$$ \widehat{V}_{\mathrm{CC}} = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} b_{kk'} \frac{W_k W_{k'}}{ \pi_{kk'} } Y_k Y_{k'}, $$
which we refer to as the connected components (or CC) estimator.
\begin{prop} \label{prop:cc-valid}
Suppose $C > 1$. Then
\begin{equation}
\operatorname{E}[\widehat{V}_{\mathrm{CC}}] = \operatorname{Var}(\hat\tau) + \frac{C}{K^2} \operatorname{Var}_c(\tau_c) \geq \operatorname{Var}(\hat\tau) \label{eqn:conservativeness-vhatcc}
\end{equation}
where $$\operatorname{Var}_c(\tau_c) = \frac{1}{C-1} \sum_c (\tau_c - \bar{\tau})^2 $$
is the finite-population variance of the component sums $\tau_c:=\sum_{k\in\mathcal{K}_c}Y_k$ across components (with $\bar{\tau}:= \frac{1}{C} \sum_c \tau_c$).
\end{prop}
As we now show, $\widehat{V}_{\mathrm{CC}}$ recovers the standard variance estimators for completely randomized and cluster-randomized experiments, where the conservativeness of the usual estimator is proportional to the treatment effect heterogeneity across either individuals or clusters.
\paragraph{\Cref{example:cre} (revisited).} In a completely randomized experiment, it is straightforward to verify that $\widehat{V}_{\mathrm{CC}}$ corresponds to the usual Neyman variance estimator, i.e. $\widehat{V}_{\mathrm{CC}} = \tfrac{1}{N_1} \hat{s}_1^2 + \tfrac{1}{N_0} \hat{s}_0^2$ where $\hat{s}_d^2 = \tfrac{1}{N_d - 1} \sum_{i:D_i=d} (Y_i - \bar{Y}_d )^2$ is the sample variance for treatment arm $d$. In this case, the slackness term $\frac{C}{K^2} \operatorname{Var}_c(\tau_c)$ is equal to $\frac{1}{N} \operatorname{Var}_i(\tau_i)$, the variance of treatment effects $\tau_i = Y_i(1)- Y_i(0)$.\footnote{For any index $\ell=1,\ldots,L$ and fixed attributes $X_\ell$, we use $\operatorname{Var}_{\ell}(X_\ell)$ to denote the finite-population variance, $\frac{1}{L-1} \sum_\ell (X_\ell - \operatorname{E}_\ell[X_\ell])^2$, where $\operatorname{E}_\ell[X_\ell]:=\frac{1}{L}\sum_\ell X_\ell$ is the average of the $X_\ell$.} To see why this is the case, recall that the connected components $c$ correspond to the set of potential outcomes for a particular unit $i$. Thus, $\tau_c = \sum_{k \in \mathcal{K}_c} Y_k = 2 \tau_i$, so the variance of $\tau_c$ across components is proportional to the variance of $\tau_i$ across units. Equation \eqref{eqn:conservativeness-vhatcc} thus reduces to the classic result that the Neyman variance estimator's conservativeness is proportional to the variance of treatment effects.
\paragraph{\Cref{example:cluster-experiment} (revisited).} In a cluster-randomized experiment with equal cluster sizes, $\widehat{V}_{\mathrm{CC}}$ has two familiar representations. First, it is equal to the Neyman variance estimator applied to the cluster-averages $\bar{Y}_g = \frac{1}{m} \sum_{i: g(i) =g} Y_i$. Second, it is equal to the CR2 cluster-robust variance estimator applied to the representation of $\hat{\tau}$ as the OLS coefficient from the regression $Y_{i} = \alpha + \tau D_{g(i)} + \epsilon_i$.\footnote{\citet{samii_equivalencies_2012} show the equivalence between the Neyman variance estimator and the HC2 variance estimator in a completely randomized experiment. It is straightforward to verify a similar equivalence for CR2 and the Neyman variance estimator applied to cluster averages.} The conservativeness is proportional to the variance of cluster-averages of treatment effects, $\operatorname{Var}_g( \bar{\tau}_g )$. (Equivalently, it is proportional to the variance of cluster-level sums of treatment effects, $\operatorname{Var}_c(\tau_c) = 4 m^2 \operatorname{Var}_g(\bar\tau_g)$, since clusters are the same size $m$.) This is because each connected component $c$ corresponds to the set of potential outcomes for a particular cluster $g$, and thus $\tau_c$ is proportional to $\bar{\tau}_g$.
\paragraph{Conservativeness in general causal settings.} More generally, \Cref{prop:cc-valid} shows that the conservativeness of $\widehat{V}_{\mathrm{CC}}$ will generically be proportional to the heterogeneity of causal effects across components in causal settings. Recall that in causal settings, all the potential outcomes associated with a particular unit $i$ must be in the same connected component, since any two potential outcomes for the same unit are jointly unobservable. If the $Y_k$ are equal to $ \omega_{i,d} Y_i(d)$ such that $\tau_i = \sum_d \omega_{i,d} Y_i(d)$ is a causal contrast, then $\tau_c = \sum_{i \in \mathcal{I}_c} \tau_i$ is the sum of the individual contrasts over the set of individuals $i$ whose potential outcomes are in component $c$ (denoted $\mathcal{I}_c$). Thus, the conservativeness of $\widehat{V}_{\mathrm{CC}}$ in causal settings will generically be proportional to the heterogeneity of treatment effect sums across connected components.\footnote{Note that if the connected components are the same size, then the variance of treatment effect sums is proportional to the variance of treatment effect averages. This is why in \Cref{example:cluster-experiment}, which assumes equal cluster sizes for simplicity, the conservativeness is proportional to $\operatorname{Var}_g(\bar\tau_g)$ as described above.}
\subsection{Validity under sampling perspective}
Although we primarily focus on the design-based setting where the potential outcomes are viewed as fixed, we now show that the CC estimator, like the Neyman variance in a completely randomized experiment, is also valid under certain data-generating processes in which the potential outcomes are viewed as stochastic. In this subsection only, we consider the setting in which potential outcomes $Y=(Y_1,\ldots,Y_K)'$ and treatments $W=(W_1,\ldots,W_K)'$ are both viewed as stochastic. We assume that the treatment assignment mechanism does not depend on the sampled values of the potential outcomes $Y$, so that $Y \perp \!\!\! \perp W$. We will write $\operatorname{E}_{Y,W}[\cdot]$ and $\operatorname{Var}_{Y,W}(\cdot)$ to denote the expectation and variance over the joint distribution of $(Y,W)$. We note that the notation in previous sections implicitly conditions on $Y$ (since it is viewed as fixed), so expressions like $\operatorname{Var}(\hat\tau)$ in earlier sections can be viewed as $\operatorname{Var}_W(\hat\tau \mid Y)$ from the perspective in this subsection.
We now describe when the CC estimator is valid from this alternative sampling perspective. By the law of total variance, we have that
\begin{align}
\operatorname{Var}_{Y,W}(\hat\tau) &= \operatorname{E}_{Y}[ \overbracket{\operatorname{Var}_W(\hat\tau \mid Y)}^{\text{Design-based variance}} ] + \operatorname{Var}_{Y}\big( \operatorname{E}_W[\hat\tau \mid Y] \big) . \label{eqn:var-tauhat-sampling}
\end{align}
On the other hand, using \Cref{prop:cc-valid}, we have that
\begin{align}
\operatorname{E}_{Y,W}[\widehat{V}_{\mathrm{CC}}] = \operatorname{E}_{Y}[\operatorname{Var}_W(\hat\tau \mid Y)] + \frac{C}{K^2} \operatorname{E}_{Y}\big( \operatorname{Var}_c(\tau_c) \big) \label{eqn:e-vcc-sampling}.
\end{align}
It is immediate that the first terms on the right-hand side of \eqref{eqn:var-tauhat-sampling} and \eqref{eqn:e-vcc-sampling} are the same, and so $\widehat{V}_{\mathrm{CC}}$ is unbiased for $\operatorname{Var}_{Y,W}(\hat\tau)$ if the second terms are the same as well. Note, however, that since $\hat\tau$ is unbiased for $\tau$, we have that
\begin{equation}
\operatorname{Var}_{Y}\big( \operatorname{E}_W[\hat\tau \mid Y] \big) = \operatorname{Var}_{Y}(\tau) = \frac{C^2}{K^2} \operatorname{Var}_{Y}\bigg( \frac{1}{C} \sum_c \tau_c \bigg) ,
\label{eqn:e-hat-tau-given-y}
\end{equation}
where the second equality uses the fact that $\tau = \frac{1}{K} \sum_k Y_k = \frac{C}{K} \frac{1}{C} \sum_c \tau_c$. We thus see that $\widehat{V}_{\mathrm{CC}}$ is unbiased if
\[ \operatorname{E}_{Y}\big( \operatorname{Var}_c(\tau_c) / C \big) = \operatorname{Var}_{Y}\bigg(\frac{1}{C} \sum_c \tau_c \bigg) .\]
The term on the right-hand side is the variance of the sample mean of the $\tau_c$, while the term on the left-hand side is the sample variance of the $\tau_c$ (with denominator $C-1$) divided by $C$, i.e. the usual squared standard-error for the sample mean of $\tau_c$ under iid sampling. It follows that if the distribution of $Y$ is such that the component-level sums $\tau_c$ are iid, then $\widehat{V}_{\mathrm{CC}}$ is unbiased for $\operatorname{Var}_{Y,W}(\hat\tau)$.
The following proposition gives a somewhat more general statement: if the sampling scheme for $Y$ is such that $\tau_c$ and $\tau_{c'}$ are uncorrelated for $c \neq c'$, then $\widehat{V}_{\mathrm{CC}}$ estimates an upper bound on $\operatorname{Var}_{Y,W}(\hat\tau)$, with the conservativeness of the upper bound corresponding to the heterogeneity in $\mu_c = \operatorname{E}_{Y}[\tau_c]$. Hence, the CC estimator is valid from the sampling perspective if the potential outcome generation process is uncorrelated across connected components, and unbiased if there is no heterogeneity in the $\mu_{c}$. We note that if the potential outcomes are generated independently at the individual level, then this necessarily implies that they are uncorrelated across connected components (since all potential outcomes for the same unit are in the same connected component).
\begin{prop} \label{prop:super-pop-interp}
Suppose that the distribution of $Y$ is such that $\operatorname{Cov}_Y(\tau_c,\tau_{c'}) = 0$ for $c \neq c'$. Assume that $C>1$. Let $\mu_c = \operatorname{E}_{Y}[\tau_c]$. Then
\[ \operatorname{E}_{Y,W}[\widehat{V}_{\mathrm{CC}}] = \operatorname{Var}_{Y,W}(\hat\tau) + \frac{C}{K^2} \operatorname{Var}_c(\mu_c) \geq \operatorname{Var}_{Y,W}(\hat\tau). \]
\end{prop}
\subsection{Admissibility and Bayesian Interpretation}
\paragraph{Admissibility.} We now show that the proposed estimator $\widehat{V}_{\mathrm{CC}}$ is admissible among the class of quadratic-form estimators by the definition proposed in \citet{harshaw_optimized_2026}. As in \citet{harshaw_optimized_2026}, we will consider the class of variance estimators that use Horvitz-Thompson weighting over observable pairs to estimate a quadratic-form in the potential outcomes: specifically, letting $\mathcal{M}$ denote the set of symmetric $K \times K$ matrices with elements $m_{kk'}=0$ for all $k,k'$ with $\pi_{kk'}=0$, we will consider the class of estimators for $\operatorname{Var}(\hat\tau)$ of the form
\begin{equation}
\widehat{V}_M = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} m_{kk'} \frac{W_k W_{k'} }{ \pi_{kk'} } Y_k Y_{k'} , \qquad M=(m_{kk'})_{k,k'} \in \mathcal{M}, \label{eqn:hat-v-m}
\end{equation}
which we call the class of quadratic-form estimators for short. It is clear that $\widehat{V}_{\mathrm{CC}}$ above is a special case of this class in which $M=B$, for $B$ the matrix with elements $b_{kk'}= v_{kk'} + s_{kk'}$. Note that
$$\operatorname{E}[\widehat{V}_M] = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} m_{kk'} Y_k Y_{k'} = \frac{1}{K^2} Y' M Y, $$
where we use the fact that $m_{kk'}=0$ if $\pi_{kk'}=0$. Recall that
$$\operatorname{Var}(\hat\tau) = \frac{1}{K^2} \sum_{k,k'} v_{kk'} Y_k Y_{k'} = \frac{1}{K^2} Y' V Y,$$
where $V$ is the matrix with $kk'$th element $v_{kk'}$. It follows that $\operatorname{E}[\widehat{V}_M] \geq \operatorname{Var}(\hat\tau)$ for all possible values of $Y$ if and only if $M - V$ is positive semi-definite. (Henceforth, we will write that $M-V$ is PSD or $M-V \succeq 0$ to denote that the matrix $M-V$ is positive semi-definite.) \citet{harshaw_optimized_2026} define an estimator $\widehat{V}_M$ to be inadmissible if there exists another upwards-biased estimator $\widehat{V}_Q$ of the form \eqref{eqn:hat-v-m} such that $\operatorname{E}[\widehat{V}_Q] \leq \operatorname{E}[\widehat{V}_M]$ for all values of $Y$ with strict inequality for some value of $Y$. More precisely, $\widehat{V}_M$ is admissible (among the class of estimators of the form \eqref{eqn:hat-v-m}) if $M-V$ is PSD and there exists no matrix $Q \neq M$ such that, $Q-V$ is PSD, $Q_{kk'}=0$ if $\pi_{kk'}=0$, and $M - Q$ is PSD.
\citet{harshaw_optimized_2026} shows that $\widehat{V}_M$ is admissible if and only if it is the solution to an optimization of the form
\begin{equation}
\min_{M} g(M-V) \text{ s.t. } M_{kk'} = 0 \text{ if } \pi_{kk'}=0 \text{ and } M - V \text{ is PSD} , \label{eqn:g-optimization}
\end{equation}
for some function $g$ that is strictly monotone in the Loewner order in the sense that $g(M_1) \geq g(M_2)$ if $M_1 - M_2$ is PSD, with strict inequality if $M_1 \neq M_2$. \citet{harshaw_optimized_2026} propose to find admissible estimators by solving the optimization \eqref{eqn:g-optimization} for particular choices of $g$, such as the trace or Frobenius norm. This leads to a convex optimization over PSD matrices, that \citet{harshaw_optimized_2026} implement using a splitting conic solver, which may be computationally demanding.
We now show that the estimator $\widehat{V}_{\mathrm{CC}}$ corresponds to a closed-form solution to the optimization \eqref{eqn:g-optimization} for a particular choice of $g$, although not one of the objectives that \citet{harshaw_optimized_2026} consider in their applications. It follows immediately that $\widehat{V}_{\mathrm{CC}}$ is admissible.
\begin{prop} \label{prop:vhatcc-admissible}
Let $\langle M_1,M_2 \rangle = \operatorname{tr}(M_1 M_2)$ be the trace inner product. Suppose $C>1$ and each connected component $\mathcal{K}_c$ contains at least two items.\footnote{This is trivially satisfied in causal settings since all potential outcomes for the same unit must be in the same connected component.} Then there exists a positive definite matrix $\Sigma$ such that $B$ is the solution to \eqref{eqn:g-optimization} for $g(M-V) := \langle \Sigma, M-V \rangle$, which is strictly monotone in the Loewner order. Hence, $\widehat{V}_{\mathrm{CC}}$ is admissible (in the class of quadratic-form estimators).
\end{prop}
\paragraph{Bayesian interpretation.} Building on the discussion in \citet{harshaw_optimized_2026}, we can attribute a Bayesian interpretation to the optimization that yields $\widehat{V}_{\mathrm{CC}}$ as a solution. Consider a Bayesian who has a prior over $Y$ and wishes to choose a quadratic-form estimator $\widehat{V}_M$ that minimizes the expected conservativeness of $\widehat{V}_M$ over the prior, subject to the constraint that $\operatorname{E}[\widehat{V}_M] \geq \operatorname{Var}(\hat\tau)$ for all values of $Y$. Note that for any fixed value of $Y$, the conservativeness of $\widehat{V}_M$ is $$\operatorname{E}[\widehat{V}_M] - \operatorname{Var}(\hat\tau) = \frac{1}{K^2} Y'(M-V)Y = \frac{1}{K^2} \operatorname{tr}(YY' (M-V)) = \frac{1}{K^2} \langle YY', M-V \rangle.$$
It follows that given a prior $Y \sim \Pi$, the expected conservativeness of $\widehat{V}_M$ over $\Pi$ is $K^{-2} \langle \operatorname{E}_{\Pi}[YY'], M-V \rangle$. Hence, the Bayesian wishes to solve the problem \eqref{eqn:g-optimization} with $g(M-V)= \langle \operatorname{E}_{\Pi}[YY'], M-V \rangle$. \Cref{prop:vhatcc-admissible} shows that $\widehat{V}_{\mathrm{CC}}$ optimizes \eqref{eqn:g-optimization} for the objective $g(M-V) = \langle \Sigma, M-V \rangle$ for a particular positive-definite choice of $\Sigma$. Hence, $\widehat{V}_{\mathrm{CC}}$ minimizes expected conservativeness for a Bayesian with a prior $\Pi$ such that $\operatorname{E}_\Pi[Y]=0$ and $\operatorname{Var}_{\Pi}(Y) = \Sigma$. It is therefore instructive to consider the moments of this prior. The proof to \Cref{prop:vhatcc-admissible} considers the class of variance matrices of the form
$$\Sigma_{\alpha} = P_{S^\perp} + \alpha A, $$
where $\alpha$ is a positive constant and $P_{S^\perp}$ is the orthogonal projection into the null space of $S$ (for $S$ the matrix with elements $s_{kk'}$). The proof to \Cref{prop:vhatcc-admissible} then shows that $\widehat{V}_{\mathrm{CC}}$ corresponds to the solution to \eqref{eqn:g-optimization} with $g(M-V) = \langle \Sigma_{\alpha}, M-V \rangle$ for $\alpha$ sufficiently small. Intuitively, this corresponds to a class of priors that expect the heterogeneity in component sums $\tau_c$ to be relatively small and to be ``equally distributed'' across components in a particular sense. Specifically, for distinct components $c,c'$, if $\operatorname{Var}_\Pi(Y) = \Sigma_\alpha$, then
$$\operatorname{Var}_\Pi(\tau_c - \tau_{c'}) = \alpha (m_c + m_{c'}),$$
where $m_c = \sum_{k,k' \in \mathcal{K}_c} A_{kk'}$ is the number of pairs in component $c$ with $\pi_{kk'}=0$. Note that for $\alpha =0$, this imposes that $\operatorname{Var}_\Pi(\tau_c - \tau_{c'}) = 0$, in which case there is no cross-component heterogeneity in the component-sums of treatment effects $\tau_c$. Moreover, if $\alpha >0$ and the $m_c$ are constant across components (as occurs in a cluster-randomized experiment with equal cluster sizes, for example), then $\Pi$ expects the heterogeneity across components to be equally-distributed across pairs in the sense that $\operatorname{Var}_\Pi(\tau_c - \tau_{c'})$ is the same for all pairs. In sum, we see that $\widehat{V}_{\mathrm{CC}}$ corresponds to the optimal choice for a class of Bayesians who expect heterogeneity across clusters to be relatively small, and who distribute that heterogeneity equally across pairs of clusters that have the same number of zero-probability links.
\begin{rem}[Bayesian justification for the Neyman variance in experiments]
Recall that in a completely randomized experiment, the estimator $\widehat{V}_{\mathrm{CC}}$ corresponds to the usual Neyman variance estimator. In this setting, all the connected components are the same size, consisting of the two potential outcomes for each unit, and $\tau_c$ is proportional to $\tau_i$, the treatment effect for unit $i$. The discussion above thus implies that the Neyman variance is optimal (among the class of quadratic-form estimators) for a Bayesian who expects treatment effect heterogeneity to be sufficiently small, and whose prior is symmetric in the sense that $\operatorname{Var}_\Pi(\tau_i - \tau_j)$ is constant across distinct pairs $i \neq j$.
\end{rem}
\section{Comparison to other estimators\label{sec:comparison-in-experiments}}
We now provide a comparison of the CC estimator to other variance estimators that are applicable for complex designs. We begin by comparing the estimators in individual- and cluster-randomized experiments (\Cref{example:cre,example:cluster-experiment} above), which are simple enough that the estimators can be derived in closed form. We then provide a more general characterization of variance estimators that are exact under the sharp null of no treatment effects, and show how the CC estimator compares to other estimators in this class.
We focus our comparisons in simple experimental settings on three variance estimators: two implementations of the \citet{harshaw_optimized_2026} approach described above, using the trace and Frobenius norms as the objective, and the \citet{aronow_estimating_2017} estimator. Before proceeding, we briefly review the construction of the \citet{aronow_estimating_2017} estimator. \citet[][AS]{aronow_estimating_2017} make the insightful observation that by Young's inequality, $-Y_k Y_{k'} \leq \frac{1}{2}(Y_k^2 +Y_{k'}^2)$. We can therefore obtain an upper bound on $\operatorname{Var}(\hat\tau)$ by replacing the unobservable terms $-Y_kY_{k'}$ in the expression for $\operatorname{Var}(\hat\tau)$ with the upper bound $\frac{1}{2}(Y_k^2 +Y_{k'}^2)$. This leads to an estimator of the form $\widehat{V}_M$ given in \eqref{eqn:hat-v-m}, where $M= V + D + A$, where $D$ is the degree matrix of the joint unobservability graph (i.e. the diagonal matrix with $D_{kk} = \sum_{k'} A_{kk'}$). We thus see that the slack matrix of the AS estimator, $D+A \succeq 0$, is precisely the signless Laplacian of the joint unobservability graph. We refer to the resulting estimator as $\widehat{V}_{\mathrm{AS}}$.
We now turn to a comparison of the CC estimator and the AS and \citet{harshaw_optimized_2026} estimators in individual- and cluster-randomized experiments. We provide a summary of the comparisons of the estimators, with detailed derivations in \Cref{sec:cre-derivations}.
\subsection{Completely-randomized experiment.}
First, consider a completely randomized experiment, with $N_0,N_1 \geq 2$. The estimator $\hat\tau$ is the usual difference-in-means estimator.
\paragraph{The CC estimator.} As noted above, the CC estimator reduces to the usual Neyman variance estimator, $\frac{1}{N_1} \hat{s}_1^2 + \frac{1}{N_0} \hat{s}_0^2$. Its conservativeness is $\Delta_{\mathrm{CC}} = \frac{1}{N} \operatorname{Var}_i(\tau_i)$.
\paragraph{The AS estimator.} The AS estimator can be written as $\widehat{V}_{\mathrm{AS}} = \frac{N-1}{N} \widehat{V}_{\mathrm{Neyman}} + \frac{1}{N} \hat\tau^2$. Its conservativeness is $\Delta_{\mathrm{AS}} = \frac{1}{N} \operatorname{E}_i[\tau_i^2]$. We thus see that while the CC estimator's conservativeness is proportional to the finite-population variance of treatment effects, the AS estimator's conservativeness is proportional to the second moment of treatment effects, $\operatorname{E}_i[\tau_i^2]$. We thus expect the AS estimator to be substantially more conservative than the CC estimator when the magnitude of the average treatment effect ($\operatorname{E}_i[\tau_i]^2$) is large relative to the heterogeneity $\operatorname{Var}_i(\tau_i)$. We note that because $\operatorname{Var}_i(\tau_i)$ is defined as the finite-population variance with denominator $(N-1)$, it is not the case that the conservativeness of the CC estimator, $ \frac{1}{N} \operatorname{Var}_i(\tau_i) = \frac{1}{N} \frac{N}{N-1} (\operatorname{E}_i[\tau_i^2] - \operatorname{E}_i[\tau_i]^2) $ is necessarily weakly smaller than that of the AS estimator, $\frac{1}{N} \operatorname{E}_i[\tau_i^2]$. However, the conservativeness of the AS estimator can never be much smaller than that of the CC estimator: $\Delta_{\mathrm{AS}} \geq (1-\tfrac{1}{N}) \Delta_{\mathrm{CC}}$, so the proportional reduction in conservativeness from using AS over CC is never more than $1/N$. By contrast, if $\tau = \operatorname{E}_i[\tau_i]$ is non-zero but treatment effects are homogeneous, then the CC estimator is unbiased for the variance of $\hat\tau$ while the AS estimator is conservative by $\frac{1}{N} \tau^2$, which can be large if the treatment effect is large in magnitude. Thus, it is possible to have $\Delta_{\mathrm{AS}} \gg \Delta_{\mathrm{CC}}$ if treatment effects can be large but relatively homogeneous.
\paragraph{Frobenius Norm Minimization.} Consider the \citet{harshaw_optimized_2026} variance estimator obtained by minimizing the Frobenius norm of $M$.\footnote{As in \citet{harshaw_optimized_2026}, we consider minimizing the Frobenius norm of $||M||_F$, rather than of the slack $||M-V||_F$. This corresponds to setting $g(S)= ||S+V||_F$ in \eqref{eqn:g-optimization}.} \Cref{prop:collapsed-as-frobenius-optimal} in the appendix shows that the AS estimator is the unique solution to the Frobenius norm optimization when the experiment is balanced $(N_1=N_0 \geq 2)$. Thus, the comparison between CC and AS also applies to a comparison of CC and Frobenius-norm minimization in a balanced experiment.
\paragraph{Trace Norm Minimization.} Next consider the \citet{harshaw_optimized_2026} estimator using the trace as the objective. Applied to the CRE, it turns out that the trace optimization does not have a unique solution. Both the CC estimator and AS estimator can be shown to be trace-optimal. If one numerically solves the optimization \eqref{eqn:g-optimization} with the trace as the objective, which estimator is obtained will depend on how the solver chooses between multiple optima.
\subsection{Cluster-randomized experiment}
Next consider a cluster-randomized experiment with \(G\) clusters each of size \(m = \tfrac{N}{G}>1\), of which \(G_1\) are assigned to treatment and \(G_0\) are assigned to control.\footnote{The derivations in \Cref{sec:cre-derivations} allow for $m=1$ to nest the completely randomized experiment as a special case. Since we have already discussed the CRE, we focus on $m>1$ in this subsection for ease of exposition.} Let
\[
\bar Y_g(d) := \frac{1}{m} \sum_{i:g(i)=g}Y_i(d), \qquad \bar Y_g := \bar Y_g(D_g), \qquad \bar\tau_g = \frac{1}{m} \sum_{i:g(i)=g} (Y_i(1)-Y_i(0)) = \bar Y_g(1)- \bar Y_g(0).
\]
The estimator \(\hat\tau\) is the difference in means across treated and control units, $\hat\tau = \frac{1}{N_1} \sum_{i : D_{g(i)} =1} Y_i - \frac{1}{N_0} \sum_{i : D_{g(i)} =0} Y_i$, where $N_1=mG_1,N_0=mG_0$ denote the number of treated and untreated units. Equivalently, it is the difference-in-means across clusters of cluster-level averages, $\hat\tau =\frac{1}{G_1} \sum_{g: D_g =1} \bar Y_g - \frac{1}{G_0} \sum_{g: D_g =0} \bar Y_g$.
\paragraph{The CC estimator.} The CC estimator reduces to the usual Neyman variance estimator applied to cluster-level averages,
\[
\widehat{V}_{\mathrm{Neyman,cl}}
=
\frac{1}{G_1}\hat s_{1,\mathrm{cl}}^2
+
\frac{1}{G_0}\hat s_{0,\mathrm{cl}}^2,
\]
where \(\hat s_{d,\mathrm{cl}}^2\) is the sample variance of the observed cluster level averages \(\bar{Y}_g\) among clusters assigned to treatment arm \(d\). This is equivalent to applying the CR2 cluster-robust standard errors to the regression of $Y_i$ on $D_{g(i)}$. The conservativeness of the CC estimator is proportional to the variance of treatment effects across clusters, $\Delta_{\mathrm{CC}}=\frac{1}{G}\operatorname{Var}_g(\bar\tau_g)$.
\paragraph{The AS estimator.} There are two possible implementations of the AS estimator in the cluster-randomized setting. First, one can directly apply the AS estimator to the individual-level outcomes, in which case a pair $Y_i(1),Y_j(0)$ is jointly unobservable if and only if $i$ and $j$ are in the same cluster $g$. Alternatively, one can view the cluster-randomized trial as a completely-randomized experiment, viewing the cluster-level averages as the potential outcomes. Then the only jointly unobservable pairs are the two potential outcomes for the same cluster, $\bar{Y}_g(1),\bar{Y}_g(0)$. Denoting the direct and collapsed version of AS as $\widehat{V}_{\mathrm{AS}}$ and $\widehat{V}_{\mathrm{AS}}^{\mathrm{collapsed}}$, we have that
\[ \widehat{V}_{\mathrm{AS}}^{\mathrm{collapsed}} = \frac{G-1}{G} \widehat{V}_{\mathrm{Neyman,cl}} + \frac{1}{G} \hat\tau^2, \qquad \widehat{V}_{\mathrm{AS}} = \widehat{V}_{\mathrm{AS}}^{\mathrm{collapsed}} + \widehat{V}_{\mathrm{within}}, \]
\noindent where
\[
\widehat{V}_{\mathrm{within}}
=
\frac{m-1}{mG^2}
\left[
\frac{G}{G_1}\sum_{g:D_g=1}s_g^2
+
\frac{G}{G_0}\sum_{g:D_g=0}s_g^2
\right],
\qquad
s_g^2
=
\frac{1}{m-1}
\sum_{i:g(i)=g}
\left(Y_i-\bar Y_g\right)^2.
\]
The collapsed AS estimator is analogous to the AS estimator in a randomized experiment, replacing individual outcomes with cluster-level averages. The direct AS estimator contains the additional term, $\widehat{V}_{\mathrm{within}}$, which is a weighted sum of the sample variances of the outcome \emph{within} clusters, $s_g^2$. Since $\widehat{V}_{\mathrm{within}} \geq 0$, with strict inequality if there is within-cluster heterogeneity, $\widehat{V}_{\mathrm{AS}}^{\mathrm{collapsed}} \leq \widehat{V}_{\mathrm{AS}}$. Thus, the direct AS estimator is dominated by the collapsed one.
The comparison between the collapsed version of AS and the CC estimator is analogous to that in the individual-level CRE. The CC estimator's conservativeness is proportional to the cross-cluster variance of treatment effects, $\frac{1}{G} \operatorname{Var}_g(\bar\tau_g)$, while the collapsed AS estimator has conservativeness proportional to the uncentered second moment of cluster-level treatment effect averages, $\frac{1}{G} \operatorname{E}_g[\bar\tau_g^2]$. The conservativeness of the collapsed AS estimator can be smaller than that of the CC estimator only by the finite-sample factor \((G-1)/G\), attained when \(\tau=0\). By contrast, if cluster-level average treatment effects are homogeneous and nonzero, then the CC estimator is unbiased for the variance while collapsed AS is conservative by \(\tau^2/G\). The direct AS estimator is more conservative than the collapsed AS estimator, adding the within-cluster heterogeneity term \(\operatorname{E}[\widehat{V}_{\mathrm{within}}]\), which is strictly positive if there is within-cluster heterogeneity in potential outcomes. We thus see directly that the direct AS estimator is inadmissible, since it is dominated by the collapsed AS estimator.
\paragraph{Frobenius norm minimization.} Similar to the individual-level CRE, for balanced cluster assignment ($G_0=G_1$), the Frobenius norm optimization reduces to the collapsed AS estimator. The comparison of the CC estimator and collapsed AS estimator above thus also applies to a comparison of the CC estimator and Frobenius-norm optimizing estimator.
\paragraph{Trace norm minimization.} As in the individual-level CRE, trace minimization is not unique. Both the CC estimator and the collapsed AS estimator are among the optimal solutions. The non-collapsed AS estimator is not trace-optimal when \(m>1\), as it is dominated by the collapsed AS estimator. Thus, if one solves the trace optimization numerically, the selected estimator depends on how the solver chooses among multiple solutions.
\subsection{Comparisons in more complex designs\label{subsec:complex-designs}}
We next consider a comparison of the CC estimator versus other estimators in more complex designs.
We first note that, by definition, two admissible estimators cannot uniformly dominate each other for all sets of potential outcomes. As discussed above, the CC estimator is admissible (among quadratic-form estimators), as is any estimator that solves \eqref{eqn:g-optimization} for an objective $g$ that is strictly monotone in the Loewner order. Hence, neither the CC estimator nor the \citet{harshaw_optimized_2026} estimators based on the trace or Frobenius norms can uniformly dominate each other for all choices of potential outcomes. As shown above, the CC estimator is exact when the component-level treatment effect sums $\tau_c$ are homogeneous, and thus must perform at least as well as other quadratic-form estimators under homogeneity. Admissibility implies, however, that there must be orientations of the potential outcomes with sufficient heterogeneity in the $\tau_c$ such that the trace or Frobenius optimizers outperform the CC estimator.\footnote{This assumes that the Frobenius- and trace-optimal estimators are not equal to the CC estimator. We showed above that CC is one of the trace-optimal solutions in a CRE, for example.}
By contrast, the AS estimator need not be admissible. \citet{harshaw_optimized_2026} show, for example, that the AS estimator is inadmissible in an individual-level CRE with $N_1=N_0=1$. The following result provides a complete characterization of when the AS estimator is admissible: AS is admissible if and only if all connected components $\mathcal{K}_c$ have at most two elements.
\begin{prop} \label{prop:as-admissibility}
The Aronow-Samii estimator $\widehat{V}_{\mathrm{AS}}$ is admissible if and only if every connected component of the joint unobservability graph has 2 or fewer elements, $|\mathcal{K}_c| \leq 2$ for all $c$.
\end{prop}
Note that in causal settings, all of the potential outcomes for the same unit are jointly unobservable and hence in the same connected component. Thus, in a causal setting where we are interested in a contrast involving multiple potential outcomes for each unit, AS is admissible if and only if (a) the contrast involves two potential outcomes per unit, and (b) any pair of potential outcomes involving distinct units is jointly observable. This implies that AS is admissible in a completely randomized experiment with at least two units per treatment arm.\footnote{\citet{harshaw_optimized_2026} show that AS is inadmissible in a CRE involving only one unit per treatment arm. In this case, $Y_1(1),Y_2(1)$ are jointly unobservable, since units $i=1,2$ cannot simultaneously be treated.} It follows that the AS estimator will be inadmissible in complex experimental designs where there are jointly unobservable pairs involving multiple units, or the contrast involves more than two potential outcomes per unit, and hence there must be a quadratic-form estimator that dominates AS, although it need not necessarily be the CC estimator, or the Frobenius or trace-optimal solutions.
\paragraph{Comparison to other estimators that are exact under homogeneity.} A notable feature of the CC estimator is that it is exact under treatment effect homogeneity across connected components, which nests the special case of the sharp null of no treatment effects, since then $\tau_c =0$ for all $c$, and so $\operatorname{Var}_c(\tau_c)=0$. Our next result characterizes the class of quadratic-form variance estimators that are exact under the sharp null, as well as the subclass that are also exact under treatment effect homogeneity across components. This characterization implies a sense in which the CC estimator is optimal in the subclass of estimators that are exact under component homogeneity. For this comparison, we focus on the setting where the $Y_k$ correspond to a weighted average treatment effect. To be precise, we will say that the $Y_k$ correspond to a weighted average treatment effect if there exist potential outcomes $Y_i(q)$ for $i=1,\ldots,N$, $q=0,\ldots,Q$ such that $Y_k = \omega_{i,q} Y_i(q)$, where the weights $\omega_{i,q}$ are non-zero and sum to 0 within unit, $\sum_q \omega_{i,q} = 0$, and such that $W_k = 1\{D_i = q\}$.\footnote{One could allow for the number of potential outcomes to depend on the unit, so that unit $i$ has $Q_i$ potential outcomes, at the cost of additional notation. This would accommodate network settings where it may not be possible for certain units to realize certain exposure mappings.} We will then say that a variance estimator $\widehat{V}_M$ of the form in \eqref{eqn:hat-v-m} is \emph{sharp-null-exact} if for any $Y$ corresponding to a weighted average treatment effect for potential outcomes such that $Y_i(q)=Y_{i}(q')$ for all $i$ and all $q,q'$, it is the case that $\frac{1}{K^2}Y' M Y= \frac{1}{K^2}Y' V Y = \operatorname{Var}(\hat\tau)$. We will further say that $\widehat{V}_M$ is exact under component homogeneity if $\frac{1}{K^2}Y' M Y= \frac{1}{K^2}Y' V Y = \operatorname{Var}(\hat\tau)$ whenever $\tau_c = \sum_{k:k \in \mathcal{K}_c} Y_k$ is constant in $c$. We then have the following result, which characterizes the set of valid estimators $\widehat{V}_M$ that are sharp null exact and exact under component homogeneity.
\begin{prop} \label{prop:sharp-null-exact}
Suppose that the $Y_k$ correspond to a weighted average treatment effect, and consider an estimator $\widehat{V}_M$ that takes the form in \eqref{eqn:hat-v-m}.
\begin{enumerate}
\item \emph{(Characterization of sharp-null exact $\widehat{V}_M$).} $\widehat{V}_M$ is weakly upwards-biased, $\operatorname{E}[\widehat{V}_M] \geq \operatorname{Var}(\hat\tau)$, and sharp-null-exact if and only if there exists a $C \times C$ correlation matrix $\Lambda$ such that
\begin{equation} T(M) := M-V = H \Lambda H', \quad \Lambda \succeq 0, \quad \Lambda_{cc} =1 \text{ for all } c, \label{eqn:h-lambda-h}
\end{equation}
where $H$ is the $K \times C$ matrix indicating component-membership, i.e. such that $H_{kc} = 1\{k \in \mathcal{K}_c\}$. The conservativeness of any such $\widehat{V}_M$ depends only on the component-level sums, $\tau_c = \sum_{k: k \in \mathcal{K}_c} Y_k$. In particular,
\[ \operatorname{E}[\widehat{V}_M] - \operatorname{Var}(\hat\tau) = \frac{1}{K^2}\boldsymbol{\tau}' \Lambda \boldsymbol{\tau}, \]
where $\boldsymbol{\tau}$ is the vector with elements $\tau_c$. When $C>1$, the CC estimator corresponds to the choice $\Lambda_{CC}$ which has elements 1 on the diagonal and $-1/(C-1)$ on the off-diagonal.
\item
\emph{(CC estimator minimizes the operator norm among estimators exact under component homogeneity.)} Suppose that $\widehat{V}_M$ is weakly upwards-biased, $\operatorname{E}[\widehat{V}_M] \geq \operatorname{Var}(\hat\tau)$, and exact under component homogeneity. Then $T(M)$ takes the form in \eqref{eqn:h-lambda-h} for $\Lambda$ such that $\Lambda \mathbf{1}_C = 0$, for $\mathbf{1}_C$ the $C \times 1$ vector of 1s. Among such $\Lambda$, the choice $\Lambda_{CC}$ minimizes the operator norm $||\Lambda||_{op}$ (provided $C > 1$ so $\Lambda_{CC}$ is well-defined).
\end{enumerate}
\end{prop}
The first part of the result shows that all quadratic form variance estimators that are weakly upwards-biased and sharp-null exact have conservativeness equal to a quadratic-form of the component-level sums $\tau_c$ and a correlation matrix $\Lambda$. The CC estimator corresponds to the maximal negative equicorrelation choice with $-1/(C-1)$ on the off-diagonal. The second part of \Cref{prop:sharp-null-exact} shows that of all estimators $\widehat{V}_M$ that are unbiased when treatment effects are homogeneous across components, the choice $\Lambda_{CC}$ minimizes the operator norm.
\begin{rem}[Alternative variance estimators suggested by \Cref{prop:sharp-null-exact}]
The formulation in \eqref{eqn:h-lambda-h} also makes it straightforward to construct variance estimators that are exact under the sharp null and satisfy other criteria (e.g. minimize norms other than the operator norm). This formulation may help reduce the computational burden of implementing the numerical optimization suggested by \citet{harshaw_optimized_2026}; if we are content with a variance estimator that is sharp-null-exact, we can optimize only over the $C (C-1) / 2$ unique terms in $\Lambda$, rather than the $K(K+1)/2$ terms in $T$. \Cref{prop:sharp-null-exact} also suggests a way of constructing a variance estimator that is exact under \emph{individual}-level homogeneous effects rather than component-level homogeneous effects. Recall that the CC estimator is exact when the component-sums $\tau_c$ are homogenous. If components are the same size, then homogeneity of the $\tau_c$ is implied by homogeneous treatment effects at the individual level. However, if the components are different sizes, then homogeneity in treatment effects at the individual level will imply that $\tau_c \propto n_c$, for $n_c = |\mathcal{K}_c|$ the size of component $c$.\footnote{This assumes that $\omega_{i,d}\equiv \omega_d$, so all units get the same weight.} \Cref{prop:sharp-null-exact} implies that one could obtain a variance estimator that is exact under individual-level homogeneous effects by choosing a $\Lambda$ satisfying \eqref{eqn:h-lambda-h} and the constraint $\Lambda \mathbf{n} =0$, where $\mathbf{n}=(n_1,\ldots,n_C)'$, rather than the constraint $\Lambda \mathbf{1}_C = 0$ satisfied by the CC estimator.\footnote{It is possible to construct such a $\Lambda$ if and only if $\max_c n_c \leq \sum_c n_c /2$, so no component accounts for more than half of the total observations. This follows from the fact that $\Lambda$ is a correlation matrix and therefore $\Lambda_{cc'} \in [-1,1]$. The $c$th element of the constraint $\Lambda \mathbf{n}=0$ then yields $n_c = -\sum_{c': c' \neq c} \Lambda_{cc'} n_{c'}$, which together with the bound on $\Lambda$ implies that $n_c \leq \sum_{c': c' \neq c} n_{c'} = \sum_{c'} n_{c'} - n_c$, which yields the desired inequality. It is straightforward to verify that this condition is also sufficient by constructing the correlation matrix appropriately.}
\label{rem:alt-estimators}
\end{rem}
\paragraph{Comparison to spatial HAC estimators.} In settings with spatial treatments, similar to \Cref{example:spatial}, a common approach is to use spatial-HAC variance estimators, such as \citet{conley_gmm_1999}. Spatial-HAC estimators were originally developed from the model-based perspective, but recent work has also considered design-based justifications for spatial-HAC in specific settings. \citet{wang_design-based_2025} consider a setting similar to \Cref{example:spatial} and show that spatial HAC with a particular choice of bandwidth is valid for the variance of a H\'ajek-style estimator when (a) the site-level assignments $T_s$ are assigned independently across sites, and (b) treatment effects satisfy a spatial ``homophily'' assumption. In \Cref{sec:hac-comparison-details}, we explore the properties of spatial HAC variance estimators and compare them to the CC estimator. We show that when treatment assignments satisfy a local dependence assumption (which is implied by the Bernoulli assignment scheme in \citet{wang_design-based_2025}), spatial HAC is guaranteed to be upwards-biased for $\operatorname{Var}(\hat\tau)$ without restrictions on the potential outcomes if and only if the HAC-kernel $K_{\mathrm{HAC}}$ is PSD. When this is the case, we show that the conservativeness of the HAC estimator depends on the heterogeneity of treatment effects across a coarsening of the connected components used by the CC estimator, and thus the two estimators have similar properties. However, it need not necessarily be the case that $K_{\mathrm{HAC}}$ is PSD, in which case the HAC variance estimator requires restrictions on the potential outcomes, such as the homophily assumption in \citet{wang_design-based_2025}, to ensure its validity. We also show that under local dependence, the positive semidefinite projection of HAC is upwards-biased from the design-based perspective without restrictions on the potential outcomes, but may be conservative even under the sharp null. By contrast, the CC estimator does not require a local dependence assumption, is valid without restrictions on the potential outcomes, and is unbiased under the sharp null. On the other hand, the CC estimator requires multiple connected components ($C>1$), which is not required by spatial-HAC estimators. The two approaches therefore appear to be complementary in spatial settings.
\section{Consistency\label{sec:consistency}}
Our results so far give conditions under which $\widehat{V}_{\mathrm{CC}}$ is an unbiased estimator of an upper bound on $\operatorname{Var}(\hat\tau)$. We now formalize conditions under which it is asymptotically consistent for an upper bound on the asymptotic variance of $\hat\tau$. Recall that $\hat\tau$ is a Horvitz-Thompson estimator over items indexed by $k$, while $\widehat{V}_{\mathrm{CC}}$ is a Horvitz-Thompson estimator over items indexed by the pairs $(k,k')$ with $\pi_{kk'}>0$. Thus, general sufficient conditions for convergence of Horvitz-Thompson estimators, such as those in \citet{delevoye_consistency_2020}, apply immediately to both estimators. In what follows, we provide a simple set of sufficient conditions that yield consistency, rather than aiming for the most general set of conditions. We consider a triangular array of finite populations of growing size indexed by $K$, with outcomes $Y_{1,K},\ldots,Y_{K,K}$ and assignment probabilities $\pi_{kk',K}$ depending on the population size $K$. For ease of exposition, we leave the subscript $K$ implicit in the notation where not needed for clarity.
We begin with simple conditions under which $\hat\tau$ is consistent.
\begin{asm}[Assumptions for consistency of $\hat\tau$] \label{asm:first-order-ht-convergence}
\leavevmode
\begin{enumerate}[(i)]
\item
(Bounded operator norm of $V$) $s_K := ||V||_{op}/K \to 0$.
\item
(Bounded second moment of $Y$): $\frac{1}{K} \sum_k Y_k^2 = O(1)$.
\end{enumerate}
\end{asm}
\begin{lem} \label{lem:tauhat-consistency}
Under \Cref{asm:first-order-ht-convergence}, $\operatorname{Var}(\hat\tau) = O(s_K)$, and hence by Chebyshev's inequality, $\hat\tau - \tau = O_p(\sqrt{s_K})$.
\end{lem}
\noindent We note that $||V||_{op} \leq ||V||_{F}$, and so for (i) to hold it suffices that
$$\frac{1}{K} ||V||_F = \sqrt{\frac{1}{K^2} \sum_{k,k'} \operatorname{Cov}\bigg(\frac{W_k}{\pi_k}, \frac{W_{k'}}{\pi_{k'}} \bigg)^2 } \to 0 ,$$
so that the average squared covariance in normalized sampling indicators is small.
Having established conditions under which $\hat\tau$ is consistent, we now consider conditions under which the CC estimator is consistent for an upper bound on $\operatorname{Var}(\hat\tau)$ in the sense that $\widehat{V}_{\mathrm{CC}}/\operatorname{Var}(\hat\tau) \xrightarrow{p} r \geq 1$.
\begin{asm} \label{asm:second-order-ht-convergence}
\leavevmode
\begin{enumerate}[(i)]
\item (Limiting distribution of $\hat\tau$)
Suppose that $s_K^{-1/2} (\hat\tau - \tau) \xrightarrow{d} L$ for a non-degenerate distribution $L$ with $\operatorname{Var}(L) = \lim_{K \to \infty} s_K^{-1} \operatorname{Var}(\hat\tau) =: V^* \in (0,\infty)$.
\item (Limit of CC slack)
$s_K^{-1} (\operatorname{E}[\widehat{V}_{\mathrm{CC}}]- \operatorname{Var}(\hat\tau)) = s_K^{-1} \frac{C}{K^2} \operatorname{Var}_c(\tau_c) \to \delta \in [0,\infty)$.
\item (Bounded fourth moment of $Y$) $\frac{1}{K} \sum_k Y_k^4 = O(1)$.
\item (Limited pair-dependence in treatment indicators)
\begin{equation} \bigg(\frac{1}{K^4} \sum_{ \substack{(k,k'), (\ell,\ell') \\ \pi_{kk'}>0, \pi_{\ell \ell'}>0}} \left| \operatorname{Cov}\bigg(\frac{W_k W_{k'}}{\pi_{kk'}} b_{kk'}, \frac{W_{\ell}W_{\ell'} }{\pi_{\ell \ell'}} b_{\ell\ell'} \bigg)\right|^{4/3}\bigg)^{3/4} = o(s_{K}^2) .\label{eqn:pairwise-cov-bound}
\end{equation}
\end{enumerate}
\end{asm}
\begin{prop} \label{prop:vhat-consistency}
Suppose $C>1$ and \Cref{asm:first-order-ht-convergence,asm:second-order-ht-convergence} hold. Then $\widehat{V}_{\mathrm{CC}} = \operatorname{E}[\widehat{V}_{\mathrm{CC}}] + o_p(s_K)$ and hence $\frac{\widehat{V}_{\mathrm{CC}}}{\operatorname{Var}(\hat\tau)} \xrightarrow{p} 1+ \frac{\delta}{V^*} \geq 1$.
\end{prop}
Part (i) imposes regularity conditions on the asymptotic limit distribution of $\hat\tau$. Specifically, it imposes that after centering and normalizing $\hat\tau$ by the rate $s_K^{-1/2}$, we obtain a non-degenerate limit distribution $L$. If a central limit theorem were applicable to $s_K^{-1/2} (\hat\tau - \tau)$, then $L$ would correspond to a normal distribution with variance $V^*$. Establishing a central limit theorem typically requires stronger assumptions than establishing variance consistency, and so we simply require that $s_K^{-1/2} (\hat\tau - \tau)$ converges to \emph{some} limit distribution; we refer the reader to e.g. \citet{chandrasekhar_general_2024}, \citet{wang_design-based_2025}, \citet{xu_design-based_2026}, among others, for some sufficient conditions for central limit theorems under complicated assignment schemes. Part (i) also imposes that the variance of the limiting distribution is equal to the limit of the rescaled variance, $s_K^{-1} \operatorname{Var}(\hat\tau)$, as occurs if the sequence $s_K^{-1}(\hat\tau-\tau)^2$ is uniformly integrable. Part (ii) requires that the conservativeness of the CC estimator, which is proportional to the cross-component heterogeneity, is $O(s_K)$ and hence is the same (or smaller) order as $\operatorname{Var}(\hat\tau)$.\footnote{If we instead assume that $s_K^{-1} (\operatorname{E}[\widehat{V}_{\mathrm{CC}}]- \operatorname{Var}(\hat\tau)) \to \infty$, then under the remaining conditions of \Cref{prop:vhat-consistency}, $\widehat{V}_{\mathrm{CC}} / \operatorname{Var}(\hat\tau) \xrightarrow{p} \infty$, so the CC estimator is asymptotically valid but infinitely conservative.} Part (iii) bounds the fourth moment of $Y_k$.
Part (iv) is perhaps the most substantive assumption, which intuitively, imposes that knowing that one pair $(k,k')$ is observed cannot, on average, be very informative about whether another pair $(\ell,\ell')$ is observed. Recall that to obtain consistency of $\hat\tau$, which is a Horvitz-Thompson estimator over items indexed by $k$, we had to restrict the dependence between the events that $k$ is observed and $k'$ is observed. However, $\widehat{V}_{\mathrm{CC}}$ is a Horvitz-Thompson estimator over items indexed by observable pairs $(k,k')$. Thus, to obtain consistency, we have to restrict the dependence between pairs of pairs. More precisely, \eqref{eqn:pairwise-cov-bound} restricts the dependence of pairwise indicators multiplied by the CC-variance weights $b_{kk'}$. Note that when $v_{kk'} \approx 0$, $b_{kk'} \approx -1/(C-1)$ if $k,k'$ are in distinct components, and thus weighting by $b_{kk'}$ reduces the weight placed on covariances involving cross-component pairs whose treatment indicators are weakly correlated.
\section{Extensions}
\subsection{Extension to other estimators\label{sec:non-ht-estimators}}
For simplicity, we have focused our attention on estimating the variance of a Horvitz-Thompson estimator, and provided a variance estimator that is upwards-biased in finite samples. We now show that our approach can be applied to obtain asymptotically valid variance estimators for estimators that have an influence function representation under an appropriate asymptotic sequence of finite populations. We first provide a general sketch for how the approach can be applied to any estimator with an influence function representation under high-level conditions, and then provide lower-level conditions for the case of ordinary least squares estimators.
\paragraph{Application to estimators with influence function representation.} First, observe that a sum of the form $\frac{1}{K} \sum_{k} W_k Y_k$ can be written as a Horvitz-Thompson estimator $\frac{1}{K} \sum_k \frac{W_k}{\pi_k} \tilde{Y}_k$ for the transformed outcome $\tilde{Y}_k := \pi_k Y_k$. However, a variety of estimators can be well-approximated by a sum of influence functions,
$$ \hat\theta - \theta = \underbracket{ \frac{1}{K} \sum_k W_k \phi(Z_k; \eta)}_{=O_p(r_K)} + o_p(r_K) $$
where $\hat\theta$ is the estimator and $\theta$ its estimand, $r_K$ is the rate of convergence of $\hat\theta$, and $\phi(Z_k;\eta)$ is the influence function, which takes as input observable features $Z_k$ and a nuisance parameter $\eta$, and is such that $\operatorname{E}[\sum_k W_k \phi(Z_k;\eta)] = \sum_k \pi_k \phi(Z_k;\eta) =0$. The leading term in this expansion can be written as
$$ \hat\theta^{Lin}_K = \frac{1}{K} \sum_k W_k \phi(Z_k; \eta) = \frac{1}{K} \sum_k \frac{W_k}{\pi_k} \psi_k, \qquad \psi_k := \pi_k \phi(Z_k;\eta),$$
a Horvitz-Thompson estimator for the average of the transformed influence functions $\psi_k := \pi_k \phi(Z_k;\eta)$. Under some regularity conditions, it will typically be the case that the asymptotic variance of $r_K^{-1}(\hat\theta - \theta)$ is equal to the asymptotic variance of $r_K^{-1} \hat\theta^{Lin}_K$, so that $r_K^{-1} (\hat\theta - \theta) \xrightarrow{d} L$ where $\operatorname{Var}(L) = \lim_{K \to \infty} r_K^{-2} \operatorname{Var}(\hat\theta^{Lin}_K)$. If the $\psi_k$ were directly observed, then it would be straightforward to apply the CC estimator (or other estimators described above) to estimate $\operatorname{Var}(\hat\theta^{Lin}_K)$ by setting $Y_k= \psi_k$, i.e. using $\widehat{V}_{\psi} := \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} b_{kk'} \frac{W_{k} W_{k'}}{ \pi_{kk'}} \psi_k \psi_{k'}$. This would yield a consistent estimator for an upper bound on the asymptotic variance if the design satisfies the conditions in \Cref{prop:vhat-consistency} above with $Y_k = \psi_k$. In practice $\psi_k$ is typically not observed since it depends on the nuisance parameter $\eta$, and so we replace $\psi$ with $\hat\psi_k := \pi_k \phi(Z_k;\hat\eta) = \psi_k + \epsilon_k$, for $\epsilon_k:= \hat\psi_k - \psi_k$ the estimation error. This leads to the feasible variance estimator
$$ \widehat{V}_{\hat\psi} := \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} b_{kk'} \frac{W_{k} W_{k'}}{ \pi_{kk'}} \hat\psi_k \hat\psi_{k'} = \widehat{V}_{\psi} + \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} b_{kk'} \frac{W_{k} W_{k'}}{ \pi_{kk'}} (2 \psi_k \epsilon_{k'} + \epsilon_k \epsilon_{k'}) . $$
Under regularity conditions that control the estimation error $\epsilon_k$, the equation in the previous display will imply $\widehat{V}_{\hat\psi} = \widehat{V}_{\psi} + o_p(r_K^2)$, so that the error from plugging-in $\hat\psi$ is asymptotically negligible and hence $\widehat{V}_{\hat\psi}$ is also consistent for an upper bound. For example, suppose that $\epsilon_k = -x_k (\hat\eta - \eta)$ and $\hat\eta - \eta = O_p(r_K)$ (e.g. because $\hat\eta$ is also a Horvitz-Thompson estimator). In that case, we have that
$$\widehat{V}_{\hat\psi} = \widehat{V}_{\psi} - \underbracket{(\hat\eta - \eta)}_{=O_p(r_K)} \cdot \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} \frac{W_{k} W_{k'}}{ \pi_{kk'}} 2 b_{kk'} \psi_k x_{k'} + \underbracket{(\hat\eta-\eta)^2}_{=O_p(r_K^2)} \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} \frac{W_{k} W_{k'}}{ \pi_{kk'}} b_{kk'} x_k x_{k'},$$
and so it suffices to establish that the pair-level Horvitz-Thompson estimators for $b_{kk'} \psi_k x_{k'} $ and $b_{kk'} x_k x_{k'}$ are $o_p(r_K)$ and $o_p(1)$, respectively.
\paragraph{OLS Estimators.} We now focus on the case of ordinary least squares (OLS) estimators, and give explicit conditions that ensure consistency. Suppose each unit $k$ is associated with a scalar outcome $y_k$ and a vector of covariates $X_k$, and consider the OLS estimator $\hat\theta = \big( \frac{1}{K} \sum_k W_k X_k X_k' \big)^{-1} \frac{1}{K} \sum_k W_k X_k y_k$. For example, if $d$ is a binary treatment and $k=(i,d)$ indexes the potential outcomes $Y_i(d)$, then setting $y_k = Y_i(d)$ and $X_k = (1,d,x_i')'$, for unit-level covariates $x_i$, yields $\hat\theta$ equal to the OLS estimator from the regression of $Y_i= Y_i(D_i)$ on $(1,D_i,x_i)$ in the observed sample. Let $$\theta = \operatorname{E}\bigg[ \frac{1}{K} \sum_k W_k X_k X_k' \bigg]^{-1} \operatorname{E} \bigg[ \frac{1}{K} \sum_k W_k X_k y_k\bigg] = \underbracket{ \bigg( \frac{1}{K} \sum_k \pi_k X_k X_k' \bigg)^{-1} }_{=:\Gamma_K^{-1}} \frac{1}{K} \sum_k \pi_k X_k y_k$$
be the OLS estimand, and define the residual $\epsilon_k = y_k - X_k'\theta$, which from the least-squares first order conditions satisfies the orthogonality condition $\sum_k \pi_k X_k \epsilon_k = 0$. We can then write
\[\hat\theta - \theta =\underbracket{ \bigg( \frac{1}{K} \sum_k \frac{W_k}{\pi_k} \pi_k X_k X_{k}'\bigg)^{-1} }_{=: \hat\Gamma_K^{-1} } \underbracket{ \frac{1}{K }\sum_k \frac{W_k}{\pi_k} \pi_k X_k \epsilon_k}_{=:\hat{U}_K} . \]
Observe that $\hat\Gamma_K$ and $\hat{U}_K$ are Horvitz-Thompson estimators for averages of the outcomes $\pi_k X_k X_{k}'$ and $\pi_k X_k \epsilon_k$, and thus if the regularity conditions in \Cref{lem:tauhat-consistency} apply to these outcomes, we will have that $\hat\Gamma_K = \Gamma_K + O_p(\sqrt{s_K})$ and $\hat{U}_K = O_p(\sqrt{s_K})$. Suppose for simplicity that we are interested in a scalar contrast of the form $a'\theta$. If $\Gamma_K$ converges to a positive-definite limit $\Gamma$, then by Slutsky's lemma we expect that $\operatorname{Var}(a'\hat\theta - a'\theta) \approx a'\hat\Gamma_K^{-1} \operatorname{Var}(\hat{U}_K) \hat\Gamma_K^{-1} a$. This suggests a ``sandwich form'' estimator
$$\widehat{V}_{\epsilon} := a' \hat\Gamma_{K}^{-1} \hat\Omega_{\epsilon} \hat\Gamma_{K}^{-1} a, \quad \text{ where } \quad \hat\Omega_{\epsilon} = \frac{1}{K^2} \sum_{k,k': \pi_{kk'}>0} b_{kk'} \frac{W_k W_{k'}}{\pi_{kk'}} \pi_k X_k \epsilon_k \pi_{k'} X_{k'}' \epsilon_{k'}.$$
$\hat\Omega_{\epsilon}$ is the CC estimator for the variance of $\hat{U}_K$, i.e. a Horvitz-Thompson estimator with outcome $\pi_k X_k \epsilon_k$. The estimator $\widehat{V}_{\epsilon}$ is not feasible since it requires knowledge of the true residuals $\epsilon_k$. We therefore will consider the plug-in replacement $\widehat{V}_{\hat\epsilon}$ that replaces $\epsilon_k = y_k - X_k'\theta$ with $\hat\epsilon_k = y_k - X_k'\hat\theta$. Since $\hat\epsilon_k - \epsilon_k = -X_k'(\hat\theta - \theta)$, which is pointwise $o_p(1)$, we expect the error from plugging in $\hat\epsilon$ to be small asymptotically under appropriate regularity conditions. The following lemma formalizes sufficient conditions under which the resulting estimator $\widehat{V}_{\hat\epsilon}$ is consistent for an upper bound on the asymptotic variance.
\begin{lem}[Consistency for OLS] \label{lem:ols-consistency}
Consider an asymptotic sequence where $\dim(X)$ is fixed for all $K$. Suppose that $C>1$ for all $K$ and that
\begin{enumerate}[(i)]
\item
$s_K \to 0$.
\item
(Bounded moments for $\pi_k X_k X_k'$, $\pi_k X_k\epsilon_k$) $\frac{1}{K} \sum_k || \pi_k X_k X_k'||^4_F = O(1)$ and $\frac{1}{K} \sum_k || \pi_k X_k \epsilon_k||^4_F = O(1)$.
\item
(Limited pair dependence in treatment indicators)
Assignment probabilities satisfy \eqref{eqn:pairwise-cov-bound}.
\item
(Limiting distribution for $\hat{U}_K$)
$\Gamma_K \to \Gamma \succ 0$ and $s_K^{-1/2} \hat{U}_K \xrightarrow{d} L$ for some distribution $L$ with $\operatorname{Var}(L) = \lim_{K \to \infty} s_K^{-1} \operatorname{Var}(\hat{U}_K) =: \Omega$ for finite $\Omega \succ 0$.
\item
(Convergence of CC slack at rate $s_K$.) For $u_k = \pi_k X_k \epsilon_k$ and $z_k = \mathrm{vec}(\pi_k X_k X_k')$, define $\tau_c^{(u,z)} := \sum_{k \in \mathcal{K}_c} (u_k' , z_k')'$. We have that $s_K^{-1} \frac{C}{K^2} \operatorname{Var}_c( \tau_c^{(u,z)} ) \to S^*_{(u,z)}$ for a finite matrix $S^*_{(u,z)}$, so that the CC conservativeness for the HT estimator for $\frac{1}{K} \sum_k (u_k',z_k')'$ is $O(s_K)$.
\end{enumerate}
Then for any vector $a \neq 0$, we have that
$$a' s_K^{-1/2} (\hat\theta - \theta) \xrightarrow{d} L^* := a' \Gamma^{-1} L, \text{ where } \operatorname{Var}(L^*) = a' \Gamma^{-1} \Omega \Gamma^{-1} a$$
and $ s_K^{-1} \widehat{V}_{\hat\epsilon} \xrightarrow{p} V^{ub} := \operatorname{Var}(L^*) + a' \Gamma^{-1} S_u^* \Gamma^{-1} a \geq \operatorname{Var}(L^*)$, where $S^*_u$ is the upper left block of $S^*_{(u,z)}$.
\end{lem}
\subsection{Using arbitrary clustering instead of connected components\label{sec:arbitrary-clusters}}
The construction of the CC estimator requires that the joint unobservability graph has at least two connected components. We now show how a similar idea can be extended to the case where the joint unobservability graph is connected, but one has a partition of the units into clusters $\tilde{\mathcal{K}}_1,\ldots, \tilde{\mathcal{K}}_{\tilde{C}}$ for some $\tilde{C}>1$. We then show how one can construct a variance estimator whose conservativeness depends on (a) the variance of sums $\tau_{\tilde{c}}$ across the $\tilde{\mathcal{K}}_{\tilde{c}}$, and (b) the conservativeness of a separate upper bound (e.g. Young's inequality) applied to pairs of $(k,k')$ that cross clusters and have $\pi_{kk'}=0$. The estimator will therefore tend to be close to exact if cross-cluster heterogeneity is small \emph{and} there are relatively few links across clusters.
To construct the estimator, observe that for any partition $\tilde{\mathcal{K}}_1,\ldots, \tilde{\mathcal{K}}_{\tilde{C}}$ with $\tilde{C}>1$, we can write the variance of $\hat\tau$ as
\begin{align*}
\operatorname{Var}(\hat\tau) = \bigg( \overbracket{\frac{1}{K^2} \sum_{\substack{ k,k' \\ \pi_{kk'} >0 }} v_{kk'} Y_k Y_{k'}}^{\text{Jointly observable pairs}} + \overbracket{ \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} =0 \\ k,k' \in \tilde{\mathcal{K}}_{\tilde{c}}} } (-1) Y_k Y_{k'} }^{\text{Same-cluster unobservable pairs}}\bigg) + \overbracket{\frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} =0 \\ k \in \tilde{\mathcal{K}}_{\tilde{c}}, k' \in \tilde{\mathcal{K}}_{\tilde{c}'} \neq \tilde{\mathcal{K}}_{\tilde{c}} } } (-1) Y_k Y_{k'} }^{\text{Cross-cluster unobservable pairs}} . \label{eqn:vhat-within-across-cluster}
\end{align*}
We will denote the term in parentheses by $V_{\text{within}}$ and the term outside the parentheses $V_{\text{across}}$. Observe that if the $\tilde{\mathcal{K}}_{\tilde{c}}$ are the connected components, then the $V_{\text{across}}$ term vanishes since there are no cross-component pairs with $\pi_{kk'}=0$. We now show that an analog to the CC estimator that uses the $\tilde{\mathcal{K}}_{\tilde{c}}$ instead of the $\mathcal{K}_{c}$ provides an upwards-biased estimator for $V_{\text{within}} - \frac{1}{\tilde{C}-1} V_{\text{across}}$ whose conservativeness depends on the heterogeneity across the $\tilde{\mathcal{K}}_{\tilde{c}}$. Specifically, define $s^{\tilde{\mathcal{K}}}_{kk'}$ equal to 1 if $k,k' \in \tilde{\mathcal{K}}_{\tilde{c}}$ for some $\tilde{c}$ and equal to $-1/(\tilde{C}-1)$ otherwise, and define
$$\widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}} =\frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} >0 }} (v_{kk'}+s^{\tilde{\mathcal{K}}}_{kk'}) \frac{W_k W_{k'}}{ \pi_{kk'} } Y_k Y_{k'}. $$
Note that the estimator $\widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}}$ corresponds exactly to $\widehat{V}_{\mathrm{CC}}$ when the $\tilde{\mathcal{K}}_{\tilde{c}}$ are the connected components. The following lemma characterizes the expectation of $\widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}}$.
\begin{lem}[Expectation of $\widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}}$]
\label{lem:expectation-vwithin}
$$\operatorname{E}[\widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}}] = V_{\text{within}} - \frac{1}{ \tilde{C} -1 } V_{\text{across}} + \frac{\tilde{C}}{K^2} \operatorname{Var}_{\tilde{c}}(\tau_{\tilde{c}}) $$
where $\tau_{\tilde{c}} = \sum_{k \in \tilde{\mathcal{K}}_{\tilde{c}}} Y_k$ is the within-cluster sum and $\operatorname{Var}_{\tilde{c}}(\cdot)$ denotes the finite-population variance across the $\tilde{\mathcal{K}}_{\tilde{c}}$ (with denominator $\tilde{C}-1$).
\end{lem}
It follows that we can obtain an upwards-biased estimator for $\operatorname{Var}(\hat\tau)$ if we can obtain an upwards-biased estimator for $V_{\text{across}}$. One simple approach would be to use Young's inequality, as in AS, but applied only to the cross-cluster components. That is, we could use
$$\widehat{V}_{\text{across}} = \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} =0 \\ k \in \tilde{\mathcal{K}}_{\tilde{c}}, k' \in \tilde{\mathcal{K}}_{\tilde{c}'} \neq \tilde{\mathcal{K}}_{\tilde{c}} } } \frac{1}{2}\big( \frac{W_k}{\pi_k} Y_k^2 + \frac{W_{k'}}{\pi_{k'}} Y_{k'}^2 \big) $$
which has expectation
$$\operatorname{E}[\widehat{V}_{\text{across}} ]= V_{\text{across}} + \frac{1}{K^2} \sum_{ \substack{ k,k' \\ \pi_{kk'} =0 \\ k \in \tilde{\mathcal{K}}_{\tilde{c}}, k' \in \tilde{\mathcal{K}}_{\tilde{c}'} \neq \tilde{\mathcal{K}}_{\tilde{c}} } } \frac{1}{2}\big( Y_k +Y_{k'})^2 =: V_{\text{across}} + \Delta_{\text{Young},\text{across}}. $$
If a small share of the pairs $(k,k')$ are cross-cluster with $\pi_{kk'}=0$, then we will generally expect $ \Delta_{\text{Young},\text{across}}$ to be small, since it is normalized by $1/K^2$. The variance estimator $\widehat{V}^{\tilde{\mathcal{K}}} := \widehat{V}^{\tilde{\mathcal{K}}}_{\text{within}} + \frac{\tilde{C}}{\tilde{C}-1} \widehat{V}_{\text{across}}$ is then upwards-biased for $\operatorname{Var}(\hat\tau)$ with conservativeness $\frac{\tilde{C}}{K^2} \operatorname{Var}_{\tilde{c}}(\tau_{\tilde{c}}) + \frac{\tilde{C}}{\tilde{C}-1} \Delta_{\text{Young},\text{across}} $. We therefore expect the conservativeness of this estimator to be relatively small if it is true that both (a) the $\tau_{\tilde{c}}$ are close to homogeneous across clusters, so that $\operatorname{Var}_{\tilde{c}}(\tau_{\tilde{c}})$ is small, and (b) there are relatively few cross-cluster pairs $(k,k')$ with $\pi_{kk'}=0$, so that $ \Delta_{\text{Young},\text{across}} $ is small. Other approaches to obtain an upper bound on $V_{\text{across}}$ are possible, too: for example, one could solve an optimization similar to that in \citet{harshaw_optimized_2026} to find a feasible quadratic-form estimator for $V_{\text{across}}$ that minimizes a particular matrix norm.
While the CC estimator is exact when components are homogeneous, the estimator $\widehat{V}^{\tilde{\mathcal{K}}}$ just described will generally be conservative owing to the $\Delta_{\text{Young},\text{across}}$ term. It appears, however, that some conservativeness even under homogeneous effects is necessary when there is only one connected component, at least if we focus on quadratic-form estimators. Specifically, \Cref{prop:sharp-null-exact} above characterizes the set of quadratic-form variance estimators for causal contrasts that are exact under the sharp null of no treatment effects. A corollary of this result is that the only quadratic-form variance estimator that is sharp-null-exact when $C=1$ is the variance estimator $\widehat{V} = \hat\tau^2$.
\begin{cor} \label{cor:cequals1-trivial-sharp-null-exact}
Consider the setting in \Cref{prop:sharp-null-exact} and suppose $C=1$. Then $\widehat{V}_M$ is weakly upwards-biased, $\operatorname{E}[\widehat{V}_M] \geq \operatorname{Var}(\hat\tau)$, and sharp-null-exact if and only if $\widehat{V}_M = \hat\tau^2$.
\end{cor}
We note that since $\operatorname{Var}(\hat\tau) = \operatorname{E}[\hat\tau^2] - \tau^2$, the variance estimator $\hat\tau^2$ is trivially upwards-biased and exact when $\tau=0$. It is not a very useful variance estimator, however. For one, if $\tau \neq 0$, then it has conservativeness $\tau^2$, regardless of sample size, and therefore does not converge to zero even as $K \to \infty$. Moreover, the implied $t$-statistic always has absolute value equal to 1, and therefore it is not useful for $t$-based inference. \Cref{cor:cequals1-trivial-sharp-null-exact} shows that any quadratic-form estimator that does not have these undesirable properties will necessarily not be sharp-null exact when $C=1$.
\begin{rem}[Choice of clusters]
The above discussion allows an arbitrary partition of units into clusters $\tilde{\mathcal{K}}_1,\ldots, \tilde{\mathcal{K}}_{\tilde{C}}$. This raises the question of how to choose the clusters in practice. To minimize conservativeness, the researcher would ideally choose clusters that (a) have relatively homogeneous effects and (b) for which there are few cross-cluster pairs $k,k'$ with $\pi_{kk'}=0$. In some settings, there may be natural choices that approximately satisfy these criteria. For example, in spatial settings there may be geographic units (say counties) that we expect to have relatively homogeneous effects and for which there are relatively few cross-border pairs of jointly unobservable potential outcomes. In principle, one could also choose clusters algorithmically based on the joint unobservability graph, e.g. using balanced graph-partitioning methods \citep{kernighan_efficient_1970,buluc_recent_2016} to find clusters of approximately equal size with few cross-cluster edges. We leave a formal analysis of the optimal choice of algorithm to future work.
\end{rem}
\section{Simulations and Applications\label{sec:sims-and-apps}}
We now consider simulations and applications based on two applied settings, one involving interference and the other a spatial treatment, similar to \Cref{example:spatial,example:interference}.
\subsection{A setting with interference}
We first revisit the simulation designs in \citet{harshaw_optimized_2026}, which are based on the social network experiment in \citet{paluck_changing_2016}.
\paragraph{Setting.} As a brief review, \citet{paluck_changing_2016} conduct a randomized trial in schools intended to reduce bullying and conflict. They employ a two-stage stratified randomization procedure, in which 28 of 56 schools are randomly selected for the intervention; among the treated schools, a random subset of students are selected to participate in an anticonflict intervention. Prior to the intervention, the researchers collected information on the social network of students (based on questions about which students they spent time with). In this setting it is natural to expect that the intervention may have spillovers across students. \citet{harshaw_optimized_2026} specify an exposure mapping where student $i$'s outcome is assumed to depend on (i) whether they are in a treated school, (ii) whether they were directly selected for the intervention, and (iii) whether any of their friends were selected for the intervention.\footnote{This is a somewhat simplified version of the exposure mapping in \citet{paluck_changing_2016}, who further subdivide (iii) into levels for whether any peer was treated and for whether any popular peer was treated.}
\paragraph{Simulation designs.} \citet{harshaw_optimized_2026} focus on two treatment contrasts. First, they study the $(1,1,0)$ vs. $(0,0,0)$ contrast, which corresponds to the effect of being directly treated without any treated peers versus being in a control school (i.e. the ``direct effect'' of treatment). Second, they study the $(1,1,0)$ vs. $(1,0,1)$ contrast, which compares the effect of being directly treated (without any treated peers) to not being directly treated but having a peer treated (i.e. a comparison of the ``direct'' vs. ``indirect'' effects). \citet{harshaw_optimized_2026} approximate the $\pi_{kk'}$ by simulating 5 million times from the assignment mechanism. For each treatment contrast, they restrict attention to students who have at least a 0.005 probability of receiving both treatment levels in the contrast. They construct point estimates of the relevant contrasts using the Horvitz-Thompson estimator. They then implement five variance estimators, the \citet{aronow_estimating_2017} estimator and the variance estimator defined by the optimization in \eqref{eqn:g-optimization} for four different objectives. For simplicity of presentation, we focus on the AS estimator and those based on the trace and ``composite'' objectives, which of the four objectives considered in \citet{harshaw_optimized_2026} yield the most and least conservative variance estimators in their main results reported in Table 1. (The composite objective is a weighted combination of the trace and Frobenius norms.) \citet{harshaw_optimized_2026} consider four designs for the potential outcomes. In two of the designs, they impose the sharp null of no treatment effects by setting all potential outcomes equal to an observed outcome in \citet{paluck_changing_2016}. They do this for two observed outcomes: whether a student wore an anticonflict wristband, and whether there were any disciplinary actions against the student. In the third simulation design, they used ``mixed'' potential outcomes where one of the potential outcomes is set to the observed Wristband outcome and the other is set to the observed Disciplinary outcome. The final simulation design is based on drawing the potential outcomes as $(X' \beta_0 +\epsilon_0, X' \beta_1 + \epsilon_1)$ where $X$ are observed covariates and the $\epsilon$'s are independent Gaussian noise (see \citet{harshaw_optimized_2026} for details on the choice of $X$ and $\beta$). The first two designs have homogeneous (zero) treatment effects, whereas the third and fourth designs have heterogeneous effects. To explore how the estimators behave under homogeneous but non-zero effects, we also consider a modification to the sharp-null simulations where we add a constant effect $\tau$ to the $Y(1,1,0)$ potential outcome for all units. Below we focus on $\tau=0.5$, which is close to the estimated average direct effect for the wristband outcome in the \citet{paluck_changing_2016} data.
\paragraph{Simulation implementation.} We implement the AS and trace and composite objective estimators identically to the simulations in \citet{harshaw_optimized_2026} using their replication code.\footnote{\citet{harshaw_optimized_2026} impose a constraint such that the variance estimator does not put any weight on pairs with $\pi_{kk'}<0.002$, and so we impose this constraint on these estimators as well.} For comparison, we implement the CC variance estimator using the same design. For the $(1,1,0)$ vs. $(0,0,0)$ contrast, there are 56 connected components corresponding exactly to the 56 schools in the data. This is because it is not possible for one student in a school to be directly treated while another student in the same school is in a control school; hence, all the potential outcomes for students in the same school are connected in the joint unobservability graph. For the $(1,1,0)$ vs. $(1,0,1)$ contrast, there are 468 connected components. The connected components do not reduce to the school for this outcome because the two relevant potential outcomes are not necessarily jointly unobservable for students in the same school; it is possible that unit $i$ in one school is directly treated with no treated neighbors and unit $j$ in the same school is not directly treated but has treated neighbors.\footnote{This is possible unless all of $j$'s friends are also $i$'s friends. In that case, $j$ can only have a treated neighbor if $i$ also has a treated neighbor. Hence, there are some joint unobservability links across units in the same school.} We conduct 100,000 simulations using the simulation designs in \citet{harshaw_optimized_2026} described above. For each variance estimator, we report the same statistics as in \citet{harshaw_optimized_2026}: (i) the bias of the variance estimator as a fraction of the true variance, namely $(\hat{\operatorname{E}}[\widehat{V}] - \widehat{\operatorname{Var}}(\hat\tau))/\widehat{\operatorname{Var}}(\hat\tau)$ (where $\hat{\operatorname{E}}$, $\widehat{V}$ denote mean and variance over the simulations), (ii) the precision of the variance estimator $\widehat{\mathrm{SD}}(\widehat{V})/ \widehat{\operatorname{Var}}(\hat\tau)$, (iii) coverage of a nominal 95\% interval of the form $\hat\tau \pm 1.96 \sqrt{\widehat{V}}$, and (iv) the average width of the nominal 95\% CI relative to that for the AS estimator.\footnote{Horvitz-Thompson-based variance estimators are generally not guaranteed to be non-negative with probability 1. If $\hat{V}<0$, we set $\sqrt{\hat{V}}=0$ when calculating the confidence intervals and average length. This occurs for the CC estimator in only 2 simulation draws across all specifications (which each use 100,000 simulations), and does not occur for the other variance estimators.}
\paragraph{Simulation results.} The results for the two contrasts are shown in \Cref{tab:simulation-results-30,tab:simulation-results-32}, respectively. Recall from \Cref{prop:cc-valid} that the CC estimator should be unbiased under the sharp null. In line with this, the estimated proportional bias in Panels A and B, which impose the sharp null, is very small in magnitude, less than 1\% of $\operatorname{Var}(\hat\tau)$ (with differences from zero likely due to simulation error). As it turns out, the CC estimator is also close to unbiased under the synthetic outcomes shown in Panel D, with upwards bias of 1-2\%. Intuitively, this occurs because the noise generated by the $\epsilon$ draws is independent across units, and thus with the specified error process there does not wind up being a large variance in the total treatment effect across components. In the Mixed PO design shown in Panel C, there is relatively more heterogeneity, and so the CC estimator has larger upward bias (64 and 16 percent, respectively). Coverage for the CC estimator is at least 93.5\% for nominal 95\% intervals in all specifications. Comparing to the other estimators, we find that the CC estimator has the lowest bias and shortest resulting confidence interval widths in all simulation designs. The gains relative to the AS estimator can be large, with the width of the CC estimator CI ranging from 44 to 84\% of that for the AS estimator across designs. The gains relative to the \citet{harshaw_optimized_2026} estimators using the trace and composite objectives are much more modest in Panels A-D (which match the original \citet{harshaw_optimized_2026} design), with the CI width relative to AS typically differing only by single-digit percentage points. We find larger reductions in average length from using the CC estimator in Panels E and F, which consider DGPs with a constant non-zero treatment effect, with differences of at least 10\% of the AS estimator length in all cases and over 30\% in one case. To better understand these results, \Cref{fig:bias-as-function-of-tau} plots the normalized bias of each estimator as a function of the constant effect size $\tau$, which shows that the bias of the CC estimator increases much more slowly in $\tau$ as for the other estimators. To gain intuition, recall that the bias of the CC estimator is proportional to $\operatorname{Var}(\tau_c)$, the variance of the sum of treatment effects across connected components. Hence, if all components were the same size, the CC estimator would be unbiased for all values of $\tau$; in practice, there is modest variation in component size in this application, and so the bias of the CC estimator is not completely insensitive to the magnitude of $\tau$, but much less sensitive than the other estimators. In addition, we find that the precision of the CC estimator is generally similar to that of the two \citet{harshaw_optimized_2026} estimators, and generally smaller than that of AS; the reductions in length do not appear to come at a large cost in terms of estimator precision. Finally, we comment briefly on computation time. Solving the optimization \eqref{eqn:g-optimization} to compute the trace-optimal and composite-objective-optimal variance estimators for the direct effect took 31 and 2 hours, respectively, on a MacStudio; for the direct vs. indirect contrast, the analogous runtimes were 9 and 2 hours. In comparison, computation of the analogous matrices for the CC and AS approach took under 1 second for both contrasts.\footnote{These runtimes compare how long it takes to compute the matrix $M$ given the second-order treatment probabilities ($\pi_{kk'}$) as input. In practice we calculate these via 5 million simulations, which takes about 5 minutes when run in parallel using 20 cores.}
\paragraph{Application results.} \Cref{tab:paluck-application} applies the four estimators studied in the simulations to the actual data from \citet{paluck_changing_2016}. We report results for four outcomes, the two studied in the simulation (Wristband and Discipline) and the two additional outcomes reported in Figure 3 of \citet{paluck_changing_2016}. These are an indicator for whether a student's friends talked to them about conflict, and an index capturing whether students perceived social norms against conflict. The first column shows the Horvitz-Thompson point estimates while the remaining columns show the standard errors (SEs) from the four estimators. Overall, we find that the CC estimator SEs are shorter than the AS and trace-optimal SEs for all eight specifications, and smaller than those for the Composite objective in 7 of the 8. The differences relative to AS are large in many cases, with the CC SEs being less than half those of the AS ones in several specifications. The CC SEs are also substantively shorter than the trace-optimal and composite-objective optimal SEs for the two contrasts involving the Wristband outcome. The CC SE is less than half the magnitude of the other estimators for the direct effect, and about one-third shorter than the composite-objective SE for the comparison of the direct and indirect effects. For the latter comparison, the implied $t$-statistic is 2.4 using the CC SE, compared to 1.97 and 1.6 for the trace and composite-objective estimators. The differences for the other outcomes are much smaller in percentage terms. The direct effect for the Wristband outcome is the largest in magnitude (and has the largest $t$-statistic), so these results are consistent with the findings in \Cref{fig:bias-as-function-of-tau} in which the gap between the CC and other estimators grew as the treatment effect became larger.
\subsection{A setting with spatial treatments\label{subsec:sims-pollmann}}
We next consider simulations based on \citet{pollmann_causal_2026}, who studies settings in which treatments occur at specific points in space and affect nearby locations, similar to \Cref{example:spatial}.
\paragraph{Setting.} Pollmann supposes there is a finite set $\mathcal{S}$ of possible treatment sites, with $T_s = 1$ denoting that site $s$ is treated. Outcomes are measured for a separate set of units indexed by $i$, and each unit's outcome $Y_i$ may be affected by treatment at nearby sites. This is a bipartite experiment in the terminology of \citet{zigler_bipartite_2021}. In his application, Pollmann is interested in the causal effect of having a nearby grocery store on the foot traffic to restaurants during the COVID-19 pandemic; hence, $\mathcal{S}$ indicates possible locations for grocery stores, while $Y_i$ denotes foot traffic to restaurant $i$. Pollmann considers the setting where $\pi_s = P(T_s =1)$ is treated as known, and assumes that the $T_{s}$ are independent across $s$. For inference, Pollmann imposes that there is some distance $d_0$ such that the outcome for unit $i$ is only affected by sites within distance $d_0$ of $i$ (see his Assumption 3(a)) . We can thus write the potential outcomes as $Y_i(\mathbf{s}_i)$ where $\mathbf{s}_i$ denotes the set of sites within distance $d_0$ of $i$ that are treated---that is $\mathbf{s}_i \subseteq \mathcal{N}_i$, where $\mathcal{N}_i$ is the set of sites within distance $d_0$ of $i$.
Pollmann considers the H\'ajek-style estimator
\begin{equation}
\hat\tau(d)
=
\frac{
\sum_{s\in\mathcal S}T_s\sum_i w_i(s,d)Y_i
}{
\sum_{s\in\mathcal S}T_s\sum_i w_i(s,d)
}
-
\frac{
\sum_{s\in\mathcal S}(1-T_s)\frac{\pi_s}{1-\pi_s}
\sum_i w_i(s,d)Y_i
}{
\sum_{s\in\mathcal S}(1-T_s)\frac{\pi_s}{1-\pi_s}
\sum_i w_i(s,d)
}, \label{eqn:tauhat-d-pollmann}
\end{equation}
where for a reference distance $d$, $w_i(s,d) \geq 0$ denotes the weight placed on site $s$ for unit $i$. For example, we may set $w_i(s,d)$ to be positive if and only if site $s$ is between 0.075 and 0.1 miles of $i$. For simplicity of exposition, we will assume that the weights $w_i(s,d)$ are non-zero only if $s \in \mathcal{N}_i$, so we place no weight on sites that are assumed not to affect the outcome of unit $i$. Under appropriate regularity conditions, Pollmann shows that $\hat\tau(d)$ is consistent for the estimand
\begin{equation}
\tau_{\mathrm{marginal}}(d) = \frac{
\sum_{\mathbf{s} \subseteq \mathcal{S} } P(\mathbf{S} = \mathbf{s} ) \sum_{s\in \mathbf{s} } \sum_i w_i(s,d) (Y_i(\mathbf{s}_i) -Y_i( \mathbf{s}_i \setminus \{s\} ))
}{
\sum_{\mathbf{s} \subseteq \mathcal{S} } P(\mathbf{S} = \mathbf{s} ) \sum_{s\in \mathbf{s} } \sum_i w_i(s,d)
} \label{eqn:tau-marginal-def}
\end{equation}
where $\mathbf{S}$ is the random set denoting which sites are treated, and $\mathbf{s}_i = \mathbf{s} \cap \mathcal{N}_i$. In words, $\tau_{\mathrm{marginal}}(d)$ is a weighted average of the marginal effect of including site $s$ in the treated set of sites for nearby restaurants.
As we describe in detail in \Cref{sec:pollmann-details}, Pollmann also shows that under appropriate regularity conditions, the H\'ajek-style estimator $\hat\tau$ is asymptotically equivalent to a Horvitz-Thompson estimator $\tilde\tau$, where
$$\tilde\tau - \tau_{\mathrm{marginal}}(d) = \frac{1}{K} \sum_{k} \frac{W_k}{\pi_k} \psi_k, $$
where $k$ indexes pairs of the form $(i,s,\mathbf{s}_i)$, $K$ is the total number of such configurations, $W_k = 1\{\mathbf{S}_i = \mathbf{s}_i\}$, $\pi_k = \mathop{}\!\textnormal{P}(W_k=1)$, and the $\psi_k$ are fixed scores defined in \eqref{eqn:defn-psi-k} in the appendix. Since $\tilde\tau$ is a Horvitz-Thompson estimator, its variance is given by the expression in \eqref{eqn:var-tauhat}, with $Y_k = \psi_k$. As we describe in \Cref{sec:pollmann-details}, Pollmann obtains an upper bound on $\operatorname{Var}(\tilde\tau)$ by using a weighted version of Young's inequality to upper bound the terms involving $\psi_k \psi_{k'}$ for pairs with $\pi_{kk'} =0$. If the $\psi_k$ were known, it would then be straightforward to estimate this upper bound using a Horvitz-Thompson variance estimator for the pairs with $\pi_{kk'}>0$; one could likewise estimate $\operatorname{Var}(\tilde\tau)$ using the CC estimator. In practice, the scores $\psi_k$ are not known exactly since they depend on $\mu_t(d), \mu_c(d)$, which are the values that are obtained if $Y_i(\mathbf{s}_i) -Y_i( \mathbf{s}_i \setminus \{s\} )$ is replaced with $Y_i(\mathbf{s}_i) $ and $Y_i( \mathbf{s}_i \setminus \{s\} )$, respectively, in the definition of $\tau_{\mathrm{marginal}}(d)$. Pollmann thus considers a feasible version of the variance estimator that replaces $\psi_k$ with $\hat\psi_k$, where $\hat\psi_k$ is the estimated score that replaces $\mu_t(d), \mu_c(d)$ with plug-in estimates $\hat\mu_t(d), \hat\mu_c(d)$ corresponding to the first and second terms in the expression for $\hat\tau(d)$. A plug-in version of the CC estimator can be constructed analogously. We also consider a leave-out version of the Pollmann and CC variance estimators where the score $\hat\psi_{k,-c}$ for potential outcome $k$ in component $c$ is constructed only using observations $k$ that are not in component $c$.
\paragraph{Data.} Pollmann uses proprietary data on foot traffic to businesses from Safegraph. However, he posted on Github a dataset that is ``simulated to resemble the spatial distributions of businesses'' in his data. We base our simulations and application on this publicly available dataset. The data contains information on the foot traffic $Y_i$ to 1,575 restaurants. It also contains the locations of 481 potential sites for grocery stores, of which 167 actually have a grocery store. For each site $s$, the data contain $\pi_s$, the probability that site $s$ would have a grocery store.\footnote{In practice, Pollmann estimates the $\pi_s$ using a convolutional neural network trained on characteristics of the sites. Following Pollmann's ``hypothetical experiment perspective'', however, we treat the $\pi_s$ as known. Pollmann also provides conditions under which one can ignore estimation error in the $\pi_s$ if one uses an estimator based on a doubly-robust score; in his Table I, he finds similar results using the doubly-robust score as he does treating the $\pi_s$ as known. We therefore focus on the hypothetical experiment perspective for simplicity, although one could likewise construct the CC estimator based on the doubly robust score.}
\paragraph{Simulation design.} We consider simulations under the sharp null of no treatment effects based on Pollmann's posted data. That is, the potential outcomes $Y_i(\mathbf{s}_i)$ are set equal to the observed value $Y_i$ in the data for all $\mathbf{s}_i$. We then simulate different data realizations by re-drawing the treatment statuses $T_s \sim \mathrm{Bernoulli}\left(\pi_s\right)$, independently across $s$. We first calculate Pollmann's variance estimator and the CC estimator based on the true scores $\psi_k$ (which are known since we know $\mu_t(d),\mu_c(d)$) and compare these relative to the estimator $\tilde\tau$. We then calculate the feasible variance estimators for $\hat\tau$ that plug-in the estimated score and compare these to the estimator $\hat\tau$. Comparing these two approaches illustrates the quality of the first-order approximation to the distribution of $\hat\tau$ based on estimated scores. Following the example, on Pollmann's Github page, we impose that there is no treatment effect for sites more than 0.1 miles away ($d_0=.1$), and we consider the weights $w_i(s,d)$ corresponding to indicators for whether the distance between $i$ and $s$ falls in the following bins: $(0,0.025],(0.025,.05], (0.05,0.075], (0.075,.1]$.
We construct the connected components as follows. Recall that the potential outcome for restaurant $i$, $Y_i(\mathbf{s_i})$, depends on the set of grocery stores within distance $d_0=0.1$ miles. It follows that restaurants $i$ and $j$ have some pair of potential outcomes that are jointly unobservable if and only if they are both within distance $d_0$ of the same site, i.e. $\mathcal{N}_i \cap \mathcal{N}_j \neq \emptyset$. The potential outcomes for units $i$ and $j$ are thus in the same connected component if there is a path of restaurants from $i$ to $j$ that each share a neighboring site. When focusing on a particular distance bin, we can restrict attention to the subset of restaurants that fall within that distance of some grocery store, and hence get non-zero weight in the estimand. The number of connected components thus differs across distance bins, ranging from 89 to 121.
\paragraph{Simulation results.} \Cref{tab:pollmann-sims-1bin} shows the simulated results for the $(0.075,0.1]$ bin. Panel A shows the results for estimating the variance of the Horvitz-Thompson estimator $\tilde\tau$ using the true scores. Since the simulation design imposes the sharp null, we expect the CC estimator to be unbiased, and consistent with this, the simulated relative bias of the CC estimator is very close to zero. Coverage for CIs based on the CC estimator is also very close to the nominal level of 95\%. By contrast, the Pollmann variance estimator is upwards-biased by approximately 7 times the actual variance of $\tilde\tau$. This leads to coverage of 100\% in our simulations. The average length of the CI based on the CC estimator is approximately three times smaller than that based on Pollmann's variance estimator. Panel B shows the results using estimated scores as estimators for $\hat\tau$. Using the naive plug-in CC estimator, we find moderate downward-bias for the CC estimator, with coverage dropping to about 92 percent for nominal 95\% intervals. This is natural, since we expect terms of the form $(Y_k - \hat\mu)^2$ to be downward biased for $(Y_k - \mu)^2$ when the outcome $Y_k$ contributes to the estimate $\hat\mu$, owing to overfitting. Using the leave-out estimates for the scores, we find that the CC estimator is now moderately upwards-biased, and coverage is close to the nominal level (94\%). Both the plug-in and leave-out versions of Pollmann's estimator remain extremely conservative, with coverage of 100 percent and average length roughly 3 times that of the CC estimators. We find similar results for the other bins (\Cref{tab:pollmann-sims-all-bins-true,tab:pollmann-sims-all-bins-feasible}).
\paragraph{Additional simulations.} In \Cref{tab:pollmann-sims-hac} in the appendix, we present some extensions of the simulation design in which we consider designs with heterogeneous treatment effects and also report results for spatial HAC estimators (\citet{conley_gmm_1999}) using the choice of bandwidth and kernel recommended in \citet{wang_design-based_2025}. We find that spatial HAC and the CC estimator both have good coverage and yield CIs of a similar length under the sharp null and under a simple form of treatment effect heterogeneity, but that spatial HAC can be strongly downward biased under treatment effect heterogeneity that violates the homophily assumption in \citet{wang_design-based_2025}. The positive semi-definite projection of HAC has coverage above the nominal level in all our simulation designs, but tends to be somewhat more conservative than the CC estimator. In all our simulation designs, Pollmann's variance estimator has large upwards bias, by a factor of 5 or more relative to the true variance.
\paragraph{Application.} \Cref{tab:pollmann-application} shows the point estimates ($\hat\tau$) and SEs when applied to Pollmann's data (which recall is simulated, but calibrated based on the real Safegraph data). We focus on comparing the SEs using leave-out scores, although the results are similar with plug-in scores. For all of the distance bins, we find that the CC SEs are approximately 3 times smaller than those using Pollmann's approach. For the first three bins, the estimates are statistically significantly different from zero at the 5\% level when using the CC SEs, while none of the estimates are significant using Pollmann's SEs.
\begin{table}[htbp]
\centering
\caption{Simulation results based on \citet{paluck_changing_2016} for the direct effect ($(1,1,0)$ vs. $(0,0,0)$ contrast)}
\label{tab:simulation-results-30}
\begin{tabular}{lrrrr}
\toprule
Estimator & Bias & Precision & Coverage & Width \\
\midrule
\multicolumn{5}{l}{\textit{Panel A: Outcome: Wristband (Sharp null)}} \\
Aronow--Samii & 0.8411 & 0.285 & 0.992 & 1.000 \\
Harshaw--trace & 0.0460 & 0.201 & 0.960 & 0.752 \\
Harshaw--composite & 0.0182 & 0.181 & 0.951 & 0.743 \\
CC & -0.0090 & 0.206 & 0.940 & 0.732 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel B: Outcome: Disciplinary (Sharp null)}} \\
Aronow--Samii & 2.0277 & 0.628 & 1.000 & 1.000 \\
Harshaw--trace & 0.0879 & 0.334 & 0.970 & 0.596 \\
Harshaw--composite & 0.0432 & 0.339 & 0.955 & 0.583 \\
CC & -0.0061 & 0.376 & 0.939 & 0.567 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel C: Outcome: Mixed POs}} \\
Aronow--Samii & 2.7559 & 0.700 & 1.000 & 1.000 \\
Harshaw--trace & 0.7573 & 0.416 & 0.992 & 0.682 \\
Harshaw--composite & 0.7520 & 0.388 & 0.991 & 0.682 \\
CC & 0.6424 & 0.440 & 0.984 & 0.659 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel D: Outcome: Synthetic Linear}} \\
Aronow--Samii & 4.0343 & 0.768 & 1.000 & 1.000 \\
Harshaw--trace & 0.1114 & 0.416 & 0.968 & 0.465 \\
Harshaw--composite & 0.0277 & 0.347 & 0.954 & 0.448 \\
CC & 0.0162 & 0.384 & 0.945 & 0.444 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel E: Outcome: Wristband (Constant effect = 0.5)}} \\
Aronow--Samii & 1.4795 & 0.541 & 0.995 & 1.000 \\
Harshaw--trace & 0.7884 & 0.497 & 0.979 & 0.846 \\
Harshaw--composite & 0.4934 & 0.414 & 0.972 & 0.773 \\
CC & 0.0050 & 0.292 & 0.940 & 0.634 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel F: Outcome: Disciplinary (Constant effect = 0.5)}} \\
Aronow--Samii & 2.5925 & 0.862 & 0.999 & 1.000 \\
Harshaw--trace & 1.5259 & 0.678 & 0.993 & 0.837 \\
Harshaw--composite & 0.9360 & 0.603 & 0.985 & 0.731 \\
CC & 0.0265 & 0.435 & 0.941 & 0.528 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes:} Results are based on 100,000 simulation draws. Bias is the bias of the variance estimator relative to the true variance. Precision is the standard deviation of the variance estimator relative to the true variance. Coverage is the coverage rate of nominal 95\% confidence intervals. Width is average confidence-interval width relative to Aronow--Samii. All estimators use assignment probabilities calculated from 5 million design draws.
\end{minipage}
\end{table}
\begin{table}[htbp]
\centering
\caption{Simulation results based on \citet{paluck_changing_2016} for comparison of direct vs. indirect effect ($(1,1,0)$ vs. $(1,0,1)$ contrast)}
\label{tab:simulation-results-32}
\begin{tabular}{lrrrr}
\toprule
Estimator & Bias & Precision & Coverage & Width \\
\midrule
\multicolumn{5}{l}{\textit{Panel A: Outcome: Wristband (Sharp null)}} \\
Aronow--Samii & 1.3665 & 1.298 & 0.992 & 1.000 \\
Harshaw--trace & 0.0554 & 0.649 & 0.947 & 0.664 \\
Harshaw--composite & 0.3428 & 0.672 & 0.970 & 0.756 \\
CC & -0.0019 & 0.661 & 0.941 & 0.644 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel B: Outcome: Disciplinary (Sharp null)}} \\
Aronow--Samii & 1.3292 & 1.650 & 0.990 & 1.000 \\
Harshaw--trace & 0.0596 & 0.810 & 0.944 & 0.671 \\
Harshaw--composite & 0.2876 & 0.793 & 0.962 & 0.749 \\
CC & -0.0007 & 0.803 & 0.935 & 0.652 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel C: Outcome: Mixed POs}} \\
Aronow--Samii & 0.6537 & 0.787 & 0.983 & 1.000 \\
Harshaw--trace & 0.1996 & 0.441 & 0.970 & 0.858 \\
Harshaw--composite & 0.1873 & 0.448 & 0.963 & 0.853 \\
CC & 0.1634 & 0.496 & 0.958 & 0.842 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel D: Outcome: Synthetic Linear}} \\
Aronow--Samii & 1.3774 & 1.503 & 0.998 & 1.000 \\
Harshaw--trace & 0.0894 & 0.733 & 0.963 & 0.674 \\
Harshaw--composite & 0.0401 & 0.645 & 0.958 & 0.662 \\
CC & 0.0174 & 0.711 & 0.950 & 0.650 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel E: Outcome: Wristband (Constant effect = 0.5)}} \\
Aronow--Samii & 1.6601 & 1.200 & 0.994 & 1.000 \\
Harshaw--trace & 1.0938 & 0.715 & 0.988 & 0.895 \\
Harshaw--composite & 0.9633 & 0.707 & 0.987 & 0.865 \\
CC & 0.1041 & 0.626 & 0.944 & 0.639 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel F: Outcome: Disciplinary (Constant effect = 0.5)}} \\
Aronow--Samii & 1.6382 & 1.377 & 0.993 & 1.000 \\
Harshaw--trace & 1.2596 & 0.789 & 0.991 & 0.937 \\
Harshaw--composite & 0.9649 & 0.757 & 0.986 & 0.872 \\
CC & 0.1264 & 0.688 & 0.945 & 0.649 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes:} Results are based on 100,000 simulation draws. Bias is the bias of the variance estimator relative to the true variance. Precision is the standard deviation of the variance estimator relative to the true variance. Coverage is the coverage rate of nominal 95\% confidence intervals. Width is average confidence-interval width relative to Aronow--Samii. All estimators use assignment probabilities calculated from 5 million design draws.
\end{minipage}
\end{table}
\begin{figure}
\includegraphics[width = 0.9 \linewidth]{harshaw-output/cc-harshaw-current/figure-analytic-bias-wristband-constant-effects.pdf}
\caption{Bias of variance estimators as a function of constant effect $\tau$ using wristband outcome}
\label{fig:bias-as-function-of-tau}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes: } This figure reports the relative bias of the variance estimators as we vary the constant effect $\tau$ for the wristband outcome. $\tau=0$ corresponds to the sharp null (Panel A of \Cref{tab:simulation-results-30,tab:simulation-results-32}) while $\tau=0.5$ corresponds to the results in Panel E of \Cref{tab:simulation-results-30,tab:simulation-results-32}. We calculate the relative bias analytically (by dividing $1/K^2 Y' (M-V) Y$ by $1/K^2 Y' V Y$), although we have verified that the results are quite similar using simulation.
\end{minipage}
\end{figure}
\begin{table}
\centering
\caption{Comparison of standard errors using actual \citet{paluck_changing_2016} data}
\label{tab:paluck-application}
\begin{tabular}{lrrrrr}
\toprule
& & \multicolumn{4}{c}{SE} \\
\cmidrule(lr){3-6}
Outcome & Estimate & CC & Trace & Comp. & AS \\
\midrule
\multicolumn{6}{l}{\textit{Panel A: Direct Effect: $(1,1,0)$ vs. $(0,0,0)$}} \\
Wristband & 0.546 & 0.043 & 0.111 & 0.091 & 0.128 \\
Discipline & -0.047 & 0.042 & 0.046 & 0.045 & 0.082 \\
Friends talked & 0.236 & 0.111 & 0.112 & 0.110 & 0.240 \\
Perceived social norm & -0.143 & 0.317 & 0.349 & 0.332 & 0.720 \\
\addlinespace
\multicolumn{6}{l}{\textit{Panel B: Direct vs Indirect Effect: $(1,1,0)$ vs. $(1,0,1)$}} \\
Wristband & 0.144 & 0.060 & 0.073 & 0.089 & 0.111 \\
Discipline & -0.038 & 0.027 & 0.029 & 0.030 & 0.042 \\
Friends talked & 0.083 & 0.127 & 0.131 & 0.145 & 0.180 \\
Perceived social norm & 0.526 & 0.373 & 0.391 & 0.403 & 0.568 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes: } This table reports the Horvitz-Thompson point estimates and SEs for different methods using the actual data from \citet{paluck_changing_2016}. The implementation of the estimators is the same as in the simulations. The Wristband and Discipline outcomes are the same as in the simulations; the remaining two outcomes are taken from Figure 3 of \citet{paluck_changing_2016}: `Friends talked' is an indicator for whether a student talked with friends about conflict; `Perceived social norm' is an index with higher values indicating a student is more anti-conflict.
\end{minipage}
\end{table}
\begin{table}
\caption{Simulation results based on \citet{pollmann_causal_2026}}
\label{tab:pollmann-sims-1bin}
\begin{tabular}{lrrrr}
\toprule
Estimator & Bias & Precision & Coverage & Width \\
\midrule
\multicolumn{5}{l}{\textit{Panel A: Variance estimates for $\tilde\tau$ using true scores}} \\
Pollmann & 7.0222 & 2.372 & 1.000 & 1.000 \\
CC & 0.0024 & 1.124 & 0.949 & 0.330 \\
\addlinespace
\multicolumn{5}{l}{\textit{Panel B: Variance estimates for $\hat\tau$ using plug-in/leave-out scores}} \\
Pollmann (plug-in) & 7.4809 & 2.127 & 1.000 & 0.988 \\
Pollmann (leave-out) & 8.1547 & 2.822 & 1.000 & 1.023 \\
CC (plug-in) & -0.0625 & 0.845 & 0.918 & 0.312 \\
CC (leave-out) & 0.2241 & 1.383 & 0.939 & 0.348 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes:} Results are based on 100,000 simulation draws. Panel A shows results where the variance of the Horvitz-Thompson $\tilde\tau$ is estimated using the true scores $\psi$. Panel B shows results where the variance of the H\'ajek-style estimator $\hat\tau$ is estimated using feasible versions of the scores that replace $\mu_t,\mu_c$ with feasible estimates. We consider two feasible variance estimators, a direct plug-in and leave-one-out version where the estimated value for units in component $c$ excludes units in component $c$. All CI widths are normalized with respect to the CI width for the Pollmann estimator for $\tilde\tau$. Results are shown for the $(0.075,.1]$ mile bin.
\end{minipage}
\end{table}
\begin{table}
\begin{tabular}{lrrrrr}
\toprule
& & \multicolumn{4}{c}{SE} \\
\cmidrule(lr){3-6}
Distance bin & Estimate & CC-PI & CC-LO & P-PI & P-LO \\
\midrule
$(0.000,0.025]$ & 6.014 & 2.413 & 2.645 & 7.188 & 7.421 \\
$(0.025,0.050]$ & 5.210 & 1.609 & 1.788 & 5.360 & 5.514 \\
$(0.050,0.075]$ & 4.566 & 1.383 & 1.549 & 4.795 & 4.950 \\
$(0.075,0.100]$ & 3.052 & 1.559 & 1.656 & 5.385 & 5.486 \\
\bottomrule
\end{tabular}
\caption{Application to Pollmann's (simulated) data} \label{tab:pollmann-application}
\begin{minipage}{0.85\textwidth}
\vspace{0.5em}
\footnotesize
\textit{Notes:} This table shows the estimated causal effect of having a grocery store at the indicated distance bin on restaurant foot traffic. We use the public version of Pollmann's data, which is simulated to match characteristics of the proprietary data he uses in his application. CC-PI and CC-LO denote the SE based on the CC variance estimator using plug-in and leave-out estimates of the scores, respectively. P-PI and P-LO are defined analogously for Pollmann's variance estimator.
\end{minipage}
\end{table}
\clearpage
\section{Conclusion}
This paper develops a variance estimator that is applicable in general settings where the treatment assignment scheme is known but may exhibit complicated dependence across units. The proposed estimator reduces to the common Neyman variance estimator in a completely randomized experiment and a cluster-robust variance estimator in a cluster-randomized experiment. It moreover generalizes three well-known properties of the Neyman variance estimator to more general designs: it is valid from the design-based perspective, its conservativeness depends on treatment effect heterogeneity, and it is also valid from a super-population sampling perspective. We find the proposed estimator performs well in simulations and two applications.
Future work could extend the analysis in this paper in a variety of directions. First, we have focused on quadratic-form variance estimators. In the context of a completely randomized experiment, \citet{aronow_sharp_2014} and \citet{abadie_sampling-based_2020} respectively show that one can obtain tighter variance estimators by using distributional information or observable covariates to lower-bound treatment effect heterogeneity and hence the conservativeness of the Neyman variance estimator. Future work might explore whether analogous approaches could be used to lower-bound the conservativeness of the CC estimator in more general designs, which likewise depends on treatment effect heterogeneity. Second, we have focused on settings where the exposure to treatment is discrete. Future work might extend the framework to settings with continuous exposures (possibly under functional form assumptions), which would, for example, incorporate shift-share settings with known shock distributions \citep{borusyak_nonrandom_2023}. Third, while we have shown consistency of the CC estimator for an upper bound on the asymptotic variance of $\hat\tau$, formal justification for normal-based confidence intervals requires a central limit theorem (CLT); such CLTs are currently available for some complicated design-based settings, but not (to our knowledge) for all settings matching the general set-up here. Future work may develop more general CLTs for this setting. Finally, our framework assumes correctly-specified exposure mappings. Future work may explore the properties of the CC estimator when the exposure mapping may be mis-specified, as in \citet{savje_average_2021} or \citet{leung_causal_2022}.
\bibliography{main.bib}