Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
114,488 characters · 7 sections · 72 citation commands
Randomization Inference of Heterogeneous Treatment Effects under Network Interference
\onehalfspacing
{3ex}
\setcounter{page}{1} \doublespacing
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Introduction} The no-interference assumption is a fundamental premise in causal inference, particularly in experiments where individuals are randomly assigned to treatments. It posits that an individual's treatment assignment does not affect the outcomes of others cox1958planning. However, this assumption is often unrealistic in settings where individuals interact, such as social networks, markets, or community-based interventions. When treatment effects propagate through such connections, interference arises, complicating causal inference. For instance, in a public health intervention, the treatment received by one individual may influence the health behaviors or outcomes of their peers. Accounting for interference is, therefore, essential to accurately characterize the data-generating process (DGP) and obtain credible causal estimates.
This paper adds to the growing literature on causal inference under network interference.\footnote{See, for example, halloran1995causal,hudgens2008toward, tchetgen2012causal, manski2013identification, aronow2017estimating, leung2020treatment, vazquez2023identification.} Our distinct focus is on developing statistical methods for testing heterogeneous treatment effects (HTEs) and zero treatment effects in experimental network data. Specifically, we propose a randomization testing procedure that applies to a broad class of non-sharp null hypotheses of constant treatment effects (CTEs) in the presence of interference.
We adopt a randomization-based inference approach for two reasons. First, the interconnected nature of social networks violates the assumption of cross-sectional independence, making classical asymptotic approximation methods challenging to apply. Second, randomization tests are nonparametric and yield exact p-values for sharp null hypotheses, where all potential outcomes can be imputed.
However, unlike the classical Fisher null hypothesis of zero treatment effect under no interference fisher1925statistical, the null hypotheses in this paper are not sharp for two fundamental reasons. First, they may depend on unknown parameters, making it possible to impute only a subset of potential outcomes under the null. Second, network interference leads to many potential outcomes, rendering the complete imputation of potential outcomes under the null hypothesis infeasible. As a result, the classical Fisher randomization testing procedure is not directly applicable.
Several papers in statistics have studied conditional randomization procedures, where inference is based on a subset of treatment assignment vectors and units, to handle non-sharp null hypotheses aronow2012general, athey2018exact, basse2019randomization, puelz2022graph.\footnote{See Section 8 of the Supplementary Material for an overview of the literature on conditional randomization tests.} The key challenge in designing such procedures lies in selecting the subset of treatment assignment vectors and units, often referred to as focal assignments and units. If focal units depend on the functions of the observed treatment assignment, they tend to vary across focal assignments, complicating randomization-based inference. Consequently, athey2018exact cautions against selecting focal units based on observed treatment assignment. However, many null hypotheses, such as those concerning constant treatment effects in the presence of interference, inherently depend on treatment assignment functions. This necessitates the development of randomization testing procedures that condition on functions of observed treatment assignment while maintaining desirable finite-sample properties. This paper makes three main contributions.
First, we introduce a broad class of null hypotheses that formalize various notions of constant treatment effects (CTEs) under network interference. These hypotheses encompass distinct forms of homogeneity and allow researchers to test for different dimensions of treatment effect constancy. When tested individually, they enable inference on specific homogeneity conditions; when tested jointly, they facilitate the identification of sources of heterogeneity. To the best of our knowledge, this is the first formal analysis of such a class of null hypotheses. Related work by owusu2024nonparametric proposes asymptotic tests for covariate-defined heterogeneous treatment effects (HTEs) in clustered networks.
Second, we develop a new randomization testing procedure for non-sharp null hypotheses. In contrast to existing approaches, our procedure selects focal units and assignments as functions of the observed treatment assignment. We establish three main theoretical results: (i) for a class of non-negative test statistics satisfying a pairwise stochastic dominance condition, we prove that the proposed test is unconditionally asymptotically valid for any nominal level $\alpha\in (0,1),$ (ii) for a broader class of test statistics, we establish asymptotic validity at some nominal levels $\alpha\in (0,1),$ and, (iii) we characterize the asymptotic size distortion of the proposed testing procedure in terms of observables.
Third, we show that our randomization procedure can be integrated with the confidence interval method of berger1994p to account for nuisance parameters under the null. We prove the asymptotic validity of the resulting procedure, further broadening the applicability of our approach.
The remainder of the paper is organized as follows. Section (ref) presents the framework and formulates the hypothesis testing problem. Section (ref) introduces the proposed testing procedure and establishes the main results. Section (ref) \!\! provides implementation guidelines. Section (ref) \!\! reports results from a Monte Carlo study. In Section (ref), we illustrate the practical application of the proposed method using the experimental data of cai2015social, testing for heterogeneous treatment effects. Section (ref) \!\! concludes. Proofs of the main results are presented in the Appendix. Extensions and proofs of auxiliary results are provided in the Supplementary Material.
\paragraph*{Notation} For any positive integer $N,$ let $[N]=\{1,2, \dots N\}.$ For any random variables/vectors $W$ and $V,$ $\Pr_{_V}$ denotes the probability with respect to $V,$ $\Pr_{_{V|W}}$ denotes the conditional probability with respect to $V$ given $W.$ We use an analogous notation for the expectation of random variables/vectors. For any set $\mathcal{B},$ $|\mathcal{B}|$ denotes the cardinality of the set. Finally, $\mathbbm{I}(\cdot)$ represents the indicator function.
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Framework}
Consistent with the framework in athey2018exact, we consider the following setting. Suppose we have a population of $N$ units (with $i$ indexing the units) connected through a single acyclic exogenous network denoted by a $N\times N$ adjacency matrix $\mathbf{A} \in \mathbbm{A},$ where $\mathbbm{A}$ is the space of possible adjacency matrices. The $(i,j)$-th element of the adjacency matrix, $A_{ij}$, equals one if units $i$ and $j$ interact, and zero otherwise. Hereafter, we refer to units $i$ and $j$ as neighbors if $A_{ij} = 1$.
Consider an experimental setting in which each unit \( i \in [N] \) is randomly assigned to one of two treatments: \( T_i = 0 \) (control) or \( T_i = 1 \) (treatment). Let \( N_0 \) and \( N_1 \) denote, respectively, the number of control and treated units, so that \( N_0 + N_1 = N \). The vector of treatment indicators denoted by \( \mathbf{T} \in \{0,1\}^N \) follows an assignment mechanism \( p : \{0,1\}^N \to [0,1] \), where \( p(\mathbf{t}) \) represents the probability that \( \mathbf{T} = \mathbf{t} \). We define the support of \( \mathbf{T} \) as $\mathcal{T}_0 = \left\{ \mathbf{t} \in \{0,1\}^N : p(\mathbf{t}) > 0 \right\}.$ For any assignment \( \mathbf{T} \in \mathcal{T}_0 \), let \( \Lambda(\mathbf{T}) \) denote a permutation of \( \mathbf{T} \) such that \( p(\Lambda(\mathbf{T})) > 0 \). Similarly, \( \Lambda(T_i) \) refers to the treatment status of unit \( i \) under the permuted assignment \( \Lambda(\mathbf{T}) \).
We adopt the Rubin--Neyman potential outcomes framework to formalize the causal setting. Specifically, there exists a mapping of potential outcomes \( \mathbf{Y} : \mathcal{T}_0 \to \mathbbm{Y}^N \subset \mathbbm{R}^N \), where \( \mathbf{Y}(\mathbf{t}) \) denotes the vector of potential outcomes corresponding to assignment \( \mathbf{t} \in \mathcal{T}_0 \). The \( i^{\text{th}} \) element of \( \mathbf{Y}(\mathbf{t}) \) is written as \( Y_i(\mathbf{t}) \), representing the potential outcome of unit \( i \) under assignment \( \mathbf{T} = \mathbf{t} \). The notation \( Y_i(\mathbf{t}) \) allows for the possibility that the potential outcome for unit $i$ may depend on the treatment assignment of other units in the population, thereby violating the classical Stable Unit Treatment Value Assumption (SUTVA) as formulated by cox1958planning.
The observed treatment assignment vector denoted by \( \mathbf{T}^{\mathrm{obs}} \) is randomly drawn from \( p(\mathbf{T}) \). We let \( \mathbf{Y}^{\mathrm{obs}} =\mathbf{Y}(\mathbf{T}^{\mathrm{obs}}) \) represent the observed outcome vector. Analogously, for each unit \( i\in [N] \), \( Y_i = Y_i(\mathbf{T}^{\mathrm{obs}}) \) denotes the observed outcome.
In addition to observed outcomes and treatment variables, we assume that for each unit \( i \in [N] \), the researcher observes a vector of $L$ pretreatment covariates \( X_i \in \mathbbm{X} \subset \mathbbm{R}^L \), where $\mathbbm{X}$ is a finite set. Let \( \mathbf{X} \) denote the \( N \times L \) matrix collecting the covariates of all units. Thus, the observed data consist the quadruple \( (\mathbf{Y}^{\mathrm{obs}}, \mathbf{T}^{\mathrm{obs}}, \mathbf{A}, \mathbf{X}) \). We adopt a design-based approach where treatment assignment is the sole source of randomness.
A defining feature of our setting is the presence of interference, whereby the potential outcome of a given unit may depend not only on its treatment status but also on the treatment assignments of other units. Interference fundamentally exacerbates the missing potential outcomes problem: without any restrictions on the nature of the interference, each unit has \( 2^N \) potential outcomes, yet only one is observed. Consequently, studying causal effects under interference typically requires additional structure to render the analysis tractable.
A common approach is to impose an exposure mapping, as proposed by aronow2017estimating. Under this framework, the dependence of each unit's outcome on the treatment of others is summarized by a low-dimensional exposure variable. For instance, leung2020treatment considered settings where the fraction of treated neighbors is a sufficient statistic for the influence of network peers' treatments on a unit's outcome.
In this paper, we adopt the exposure mapping approach and formally define the network exposure mapping as
where $ \mathbf{\Pi}$ is an arbitrary finite set of exposure values. We assume that the exposure mapping is the same for all units. However, since the network structure \( \mathbf{A} \) is fixed (exogenous), we simplify notation and write \( \pi(i, \mathbf{T}, \mathbf{A}) = \pi_i(\mathbf{T}) \), unless otherwise necessary. Moreover, we let $\Pi_i \in \mathbf{\Pi}$ denote the resulting random variable, with $\Pi_i=\pi_i(\mathbf{T}).$ Hereon, we refer to \( \Pi_i \) as the network exposure variable. Borrowing the terminology of manski2013identification, the pair \( (T_i, \Pi_i) \in \{0,1\} \times \mathbf{\Pi} \) is called the effective treatment of unit \( i \). The effective treatment jointly captures a unit's treatment status and exposure to others' treatment assignments, as determined by the network structure.
The following assumptions formally characterize the network structure and the role of exposure mapping in the framework.
Assumption (ref) restricts interference to first-order spillovers. Specifically, it allows the treatment status of a unit’s immediate neighbors to affect its potential outcome but assumes that treatments assigned to neighbors-of-neighbors have no effect. This restriction is convenient and empirically testable, as discussed in athey2018exact. Moreover, it implies that the observed adjacency indicator, $A_{ij}$, coincides with the interference dependence variable defined in savje2021average. While we adopt this formulation for concreteness, all theoretical results in this paper extend to alternative definitions of the interference dependence variable.
Assumption (ref) requires that the researcher correctly specifies the exposure mapping given the network structure. This assumption is standard in the literature on causality under interference; see, for example, aronow2017estimating and leung2020treatment. Although the assumption is strong, recent methodological contributions---such as hoshino2023randomization---develop procedures for assessing or selecting plausible exposure mappings in applied settings. We examine the consequences of exposure mapping misspecification in Section 5 of the Supplementary Material.
Given an arbitrary exposure mapping, we consider a general class of non-sharp null hypotheses defined by
These hypotheses assert that differences in potential outcomes across effective treatment realizations are governed by a function of the effective treatments and unit-level covariates. They possess several distinct features. First, the function \(\tau(\cdot)\) is typically unknown and unobserved in practice. Second, the restrictions imposed by the null do not generally permit the imputation of all potential outcomes in the population. Third, because exposure values depend explicitly on the assigned treatment vector, the set of treatment assignments consistent with the null varies across units. We elaborate on these features in Section (ref).
The hypothesis \( H^{G}_{0} \) encompasses several null hypotheses of substantive interest, which may be broadly classified into three categories: (i) those asserting no/zero treatment effects under network interference; (ii) those characterizing various forms of constant direct treatment effects under network interference; and (iii) those asserting constant indirect (spillover) effects under network interference. For brevity, we defer detailed discussion of each class to Section 7 of the Supplementary Material.
Throughout the paper, we focus on a representative null hypothesis, \( H_0 \), that illustrates the core challenges associated with the broader class \( H^{G}_{0} \). The testing procedures we develop for \( H_0 \) extend directly to other hypotheses encompassed by \( H^{G}_{0} \). Specifically, we consider $$H_0: Y_i(1, \pi_k) - Y_i(0, \pi_k) = \tau(\pi_k) \quad \text{for some function } \tau(\cdot), \text{ and for some } \pi_k \in \boldsymbol{\Pi}, \forall\, i\, \in [N]. \vspace{-0.1cm} $$
In other words, $H_0$ posits that individual-level direct effects depend solely on exposure values. Testing this hypothesis is of direct relevance for designing treatment assignment policies aimed at maximizing direct welfare; see Section 7 of the Supplementary Material for more discussion.
In contrast to the sharp null of no treatment effects considered by fisher1925statistical, the hypothesis $H_0 $ and its variants do not permit the imputation of all potential outcomes under alternative treatment assignments for any given experimental design. That is, $H_0$ is not sharp. The following example illustrates the nature of this non-sharpness.
The non-sharpness exhibited in Example (ref) is a general feature of the null hypotheses encompassed by $H_0^G.$ These hypotheses belong to a broader class referred to as partial nulls by zhang2023randomization. Because potential outcomes cannot be fully imputed under such nulls, the classical unconditional Fisher randomization test is not applicable. Moreover, as we demonstrate in the next section, existing conditional Fisher randomization procedures may exhibit low statistical power, even when applicable. These limitations motivate the development of specialized testing procedures tailored to the null hypotheses considered in this paper.
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{The Testing Procedure: Randomization Inference} A fundamental component of any testing procedure is the choice of test statistic. Throughout, we assume the existence of a predetermined two-sample-type statistic $z(\mathbf{Y}, \mathbf{T}, \mathbf{X}, \mathbf{A})$ appropriate for the null hypothesis under consideration. Under Assumptions (ref) and (ref), this statistic simplifies to \( z(\mathbf{Y}, \mathbf{T}, \boldsymbol{\pi}(\mathbf{T}), \mathbf{X}) \). The main theoretical results in this paper are derived under general conditions that allow for a broad class of test statistics. Nevertheless, the choice of statistics plays a central role in determining the power of the test.
Another essential component of randomization-based inference is the treatment assignment mechanism. We assume this mechanism is known and may follow an arbitrary but pre-specified design. The results in this paper accommodate a broad class of designs, including complete randomization, Bernoulli, paired, and stratified designs; see imbens2015causal.\footnote{Under complete randomization, a sparsity condition on the network is required to ensure that dependencies among individual treatments induced by design do not obscure the underlying network structure savje2021average.} We impose the following condition.
Assumption (ref) strengthens the standard overlap condition frequently imposed in causal inference; see imbens2015causal. It requires that, for each assignment vector generated by the experimental design, a strictly positive fraction of units within each effective treatment and covariate-defined subgroup is available to compute the test statistic. In the case of the null hypothesis $H_0$, which concerns a single exposure value $\pi_k,$ the condition is only required to hold for that value.
As shown in the main results, Assumption (ref) is sufficient to guarantee the unconditional validity of the proposed testing procedures. However, in contrast to the standard overlap condition, Assumption (ref) may fail under common experimental designs, particularly when the exposure mapping is multi-valued. One approach to address this issue is to restrict attention to treatment assignments that satisfy conditions (ref) and (ref). Consequently, the proposed procedures remain conditionally valid, even in traditional designs where Assumption (ref) does not hold globally.
Example (ref) shows that the null hypothesis \( H_0 \) is not sharp for two distinct reasons. For clarity of exposition, the next subsection focuses on testing procedures for null hypotheses that do not involve unknown parameters. Subsection (ref) then extends the analysis to nulls that include unknown parameters.
In this subsection, we test the null hypotheses in which all parameters are known a priori. Specifically, we consider a variant of $H_0:$
The null $H^{\tau}_{0}$ is similar to the CTE null hypothesis in settings without interference; see ding2016randomization and chung2021permutation. A formal comparison between the testing problem considered in this paper and those in the aforementioned references is provided in Section 6 of the Supplementary Material.
Revisiting Example (ref), observe that under \( H^{\tau}_{0} \), all missing potential outcomes in Table (ref) that were previously unobserved due to the unknown functional form of \( \tau(\cdot) \) become known. In particular, the entries marked \textcolor{blue}{!} in Table (ref) are equal to \( \tau \) under \( H^{\tau}_{0} \). Consequently, the potential outcomes of some units are observed under the randomized treatment assignment \( \Lambda(\mathbf{T}) \). This feature---referred to as “partial imputability”---motivates the conditional randomization inference (CRI) procedures commonly employed in the literature, where inference is conducted using focal units and focal assignments. Applying most existing CRI procedures presents two key technical challenges.
First, the subset of treatment assignment vectors that guarantees non-empty groups within each effective treatment arm---required to compute a two-sample test statistic---depends on the observed assignment through the realized exposure values. These exposure values determine which randomized assignments yield non-empty treatment and control groups within the relevant subpopulation. For example, in Table (ref), the observed exposure values imply that only units 1, 3, 7, 8, and 9 are relevant for testing the null \( H_0^1 \) with \( \tau(1) = \tau \). A randomized assignment qualifies as a focal assignment if and only if at least two of these units are assigned different treatments while maintaining their exposure status. Under this criterion, the assignment vector \( \tilde{\mathbf{t}} \) in Table (ref) does not qualify as a focal assignment for \( H_0^1 \) under the observed treatment, although it may qualify under a different observed assignment. Consequently, it is challenging to construct focal assignment selection rules that are both universal and independent of the observed assignment. Moreover, when focal assignment sets are defined based on exposure values, they may overlap across different observed assignments, violating the non-overlapping partition condition required for unconditional validity in CRI procedures.
Second, the subset of units for which potential outcomes can be imputed under \( H^{\tau}_0 \) may vary across both observed and randomized treatment assignments. This is because exposure values, which determine the potential outcomes that can be imputed under the null, are functions of the treatment assignment vector. As a result, the set of focal units is inherently assignment-dependent, complicating their selection in traditional CRI procedures. For example, Table (ref), which builds on Example (ref), displays two randomized treatment vectors. Under the first, only units 1 and 7 have observable outcomes under the null, while under the second, units 1, 3, 7, 8, and 9 do. Thus, the definition of focal units varies across randomized assignments. However, athey2018exact recommends against selecting focal units based on functions of the assignment vector, creating a tension between practical feasibility and theoretical validity.
Three existing CRI procedures may, in principle, be applicable to testing \( H_0^\tau \).
The first, which we refer to as the naive method, defines the focal units as those with observed exposure equal to \( \pi_k \). The focal assignments are then restricted to treatment vectors under which these units retain their observed exposure value \( \pi_k \). While straightforward to implement, this approach may yield an empty set of focal assignments, particularly in settings with complex network structures and multi-valued exposure mappings.
The second feasible procedure, known as the intersection method zhang2023randomization, partitions the space of possible treatment assignments into non-overlapping subsets, with the subset containing the observed assignment designated as the focal assignment set. The focal units are defined as the intersection of units that, across all focal assignments, maintain a fixed exposure value \( \pi_k \). This approach suffers from two main limitations. First, constructing a reasonable partition of the treatment space is non-trivial and inherently hypothesis-specific. Second, by design, the resulting set of focal units is often “near empty” or empty, leading to a test with little or no power zhang2023randomization.
The third feasible procedure is the graph-theoretic approach proposed by puelz2022graph, commonly referred to as the biclique method. This is the state-of-the-art method that constructs a bipartite graph---termed the null exposure graph---by linking units and treatment assignments that are jointly consistent with the null hypothesis. The selection of focal units and focal assignments is then reduced to identifying large bicliques within this graph. While conceptually appealing, the method suffers from a critical limitation: in complex networks with multi-valued exposure mappings, the largest bicliques are often small, resulting in tests with trivial statistical power.
In this paper, we propose a procedure that addresses the challenges outlined above while overcoming the limitations of existing feasible methods. The approach is grounded in the idea of conditional randomization but is designed to deliver unconditional validity without requiring a partition of the treatment assignment space. We begin by presenting an abstract formulation of the procedure and introducing key definitions to facilitate comparison with conventional CRI methods. The exposition then specializes to the case of \( H_0^\tau \), allowing for clearer interpretation and implementation.
Following the notation of basse2019randomization, our procedure begins by defining an event space \( \mathbbm{C} = \{\mathcal{C} = (\mathcal{S}, \mathcal{U}, \mathcal{T}) : \mathcal{S} \in \mathbbm{U}, \mathcal{U} \in \mathbbm{U}, \mathcal{T} \in \mathbbm{T} \} \), where \( \mathcal{S} \subseteq [N] \) denotes the set of super-focal units---the pool from which focal units are drawn---\( \mathcal{U} \subseteq \mathcal{S} \) denotes the set of focal units,\footnote{Focal units are referred to as imputable units in zhang2023randomization.} and \( \mathcal{T} \subseteq \mathcal{T}_0 \) denotes the set of focal assignments. Next, we specify a conditioning mechanism \( m(\mathcal{C} \mid \mathbf{T}^{\text{obs}}) \) over the event space. The mechanism is chosen such that, conditional on \( \mathcal{C} \), the predetermined test statistic \( z(\cdot) \) is well-defined and computable under the null hypothesis. This framework provides a principled basis for selecting focal units and assignments.
In its general form, the proposed conditioning mechanism decomposes as
where $\Tilde{f}$ and $\Tilde{h}$ are distributions over $\mathbbm{U}, $ and $\tilde{g}$ is a distribution over $\mathbbm{T}$. This decomposition implies that the proposed conditioning process involves the following steps:
We now sequentially apply the three-step conditioning procedure described above to the null hypothesis \(H_0^\tau\), thereby establishing the technical details of our approach.
\paragraph*{Step 1: Selection of Super-Focal Units} As \( H_0^{\tau} \) pertains only to the potential outcomes of units with exposure value \( \pi_k \), the units relevant for testing the null are those whose observed exposure equals \( \pi_k \). Formally, given an arbitrary observed assignment \( \mathbf{t}^{\text{obs}} \), we define the set of super-focal units as
Assumption (ref) ensures that $|\mathcal{S}(\mathbf{t}^{obs})|$ is sufficiently large, thereby enabling the computation of the observed test statistic for any $\mathbf{t}^{obs} \in \mathcal{T}_0$. Note that for unit $i$ in the population, the probability of being a super-focal unit is $\Pr{}_{\mathbf{T}^{obs}}(\pi_i(\mathbf{T}^{obs})=\pi_k)= \sum_{\mathbf{t} \in \mathcal{T}_0}\mathbbm{I}\{\pi_i(\mathbf{t})=\pi_k\}\cdot \Pr{}_{\mathbf{T} }(\mathbf{T}=\mathbf{t}).$ Based on Assumption (ref), these inclusion probabilities are strictly positive for all units.
\paragraph*{Step 2: Selection of Focal Assignments} Exposure values vary across treatment assignments; as a result, some randomized assignments may alter the exposure status of all units in \( \mathcal{S}(\mathbf{t}^{\text{obs}}) \), while others may leave these exposures unchanged. Assignment vectors that preserve the exposure status of all super-focal units constitute the focal assignment set in the naive method. However, as previously noted, such sets tend to have limited cardinality in complex networks with rich exposure mappings, resulting in low statistical power and increased computational burden.
Between these extremes, certain randomized assignments may change the exposure of only a subset of units in \( \mathcal{S}(\mathbf{t}^{\text{obs}}) \), while others retain exposure \( \pi_k \). Our proposed conditioning procedure selects subsets of assignments \( \mathcal{T} \subseteq \mathbbm{T} \) for which the number of super-focal units retaining their exposure status exceeds a pre-specified threshold determined by a tuning parameter \( \epsilon \).
Formally, for an arbitrary observed assignment, $\mathbf{t}^{obs},$ we define the focal assignment set for $H_0^\tau$ as:
where $$ \hat{R}(t,\mathbf{t}';\mathcal{S}(\mathbf{t}^{obs}))\coloneqq{\sum_{i\in \mathcal{S}(\mathbf{t}^{obs})} \mathbbm{I}\{t'_i=t,\pi_i(\mathbf{t}')=\pi_k \}},\,\, \text{for}\,\, t=0, 1, \vspace{-0.1cm} $$ defines the count of super-focal units assigned to treatment $t$ and exposure value $\pi_k$ under a randomized assignment $\mathbf{t}'.$ For $t\in \{0,1\},$ $\mathcal{I}_{t,\epsilon}$ denotes closed intervals whose widths are governed by the tuning parameter $\epsilon\geq 0.$ These intervals ensure the computability of the test statistics at the focal assignments. For instance, $\mathcal{I}_{t,\epsilon}$ may represent the exact central interval of the distribution of $\hat{R}(t,\mathbf{t}';\mathcal{S}(\mathbf{t}^{obs}))$ excluding extreme values $0$ and $ |\mathcal{S}(\mathbf{t}^{obs})|.$ Thus, unlike the naive method, this selection rule does not limit focal assignments to those where focal units are fixed and equal to the super-focal units.
The tuning parameter $\epsilon$ controls the balance of treated and untreated units for each focal assignment, analogous to the role of the minimum biclique size tuning parameter in the biclique method. For any $\mathbf{t} \in \mathcal{T}_0$, $\hat R(1,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs}))+ \hat R(0,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs}))\leq|\mathcal{S}(\mathbf{t}^{obs})|.$ Hence, extreme high values of $\hat R(0,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs}))$ can lead to small complementary values of $\hat R(1,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs})),$ and vice versa. To mitigate such imbalanced sample sizes---which can render test statistics non-computable or substantially reduce their precision---we select assignments where $\hat R(0,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs}))$ and $R(1,\mathbf{t};\mathcal{S}(\mathbf{t}^{obs}))$ lie within meaningful intervals, ensuring adequate sample size balance.
Larger values of $\epsilon$---corresponding to narrower intervals $\mathcal{I}_{0,\epsilon}$ and $\mathcal{I}_{1,\epsilon}$---help mitigate sample size imbalances between treatment arms. This improves the precision of the test statistic and reduces the sensitivity of the inference to variability in effective sample sizes. However, narrower intervals also restrict the focal assignment set, thereby reducing statistical power and increasing computational burden. The choice of \( \epsilon \) thus entails a fundamental trade-off between estimation precision and statistical power. Practical guidance for selecting \( \mathcal{I}_{t, \epsilon} \) is provided in Section 2 of the Supplementary Material.
It is worth noting that for any \( \mathbf{t} \in \mathcal{T}_0 \) and \( t \in \{0,1\} \), the probability that the realized sample sizes \( \hat{R}(0, \mathbf{t}; \mathcal{S}(\mathbf{T}^{\text{obs}})) \) and \( \hat{R}(1, \mathbf{t}; \mathcal{S}(\mathbf{T}^{\text{obs}})) \) fall within the respective intervals \( I_{0,\epsilon} \) and \( I_{1,\epsilon} \) is defined by
We assume that all assignments in the design have a positive probability of selection as focal assignments; that is, \( 0 < \phi_{\mathbf{t}}^\epsilon < 1 \) for all \( \mathbf{t} \in \mathcal{T}_0 \). If \( \phi_{\mathbf{t}}^\epsilon = 0 \) for some \( \mathbf{t} \), the focal assignment set must exclude those assignments, or the intervals \( I_{t,\epsilon} \) must be redefined accordingly.
\paragraph*{Step 3: Selection of Focal Units} Given the focal assignment set defined in (ref), the subset of units in \( \mathcal{S}(\mathbf{t}^{\text{obs}}) \) for which \( \pi_i(\mathbf{t}) = \pi_k \) will, in general, vary across assignments \( \mathbf{t} \in \mathcal{T}_\epsilon(\mathbf{t}^{\text{obs}}) \). As a result, the proposed focal unit selection rule implies that focal units are inherently assignment-dependent. This assignment dependence facilitates more flexible and efficient data use relative to traditional CRI methods, potentially enhancing statistical power.
Formally, for a given observed assignment \( \mathbf{t}^{\text{obs}} \), the set of focal units under any \( \mathbf{t} \in \mathcal{T}_\epsilon(\mathbf{t}^{\text{obs}}) \) is defined by
Accordingly, the probability that unit \( i \in [N] \) is selected as a focal unit can be expressed as $\Pr(\pi_i(\mathbf{T}^{obs})=\pi_k, \pi_i(\mathbf{T})=\pi_k)= \Pr{}_{\mathbf{T}}(\pi_i(\mathbf{T}^{obs})=\pi_k)\cdot Pr{}_{\mathbf{T}|\mathcal{T}_\epsilon(\mathbf{t}^{obs}) }(\pi_i(\mathbf{T})=\pi_k),$ where $\Pr{}_{\mathbf{T}|\mathcal{T}_\epsilon(\mathbf{t}^{obs}) }(\pi_i(\mathbf{T})=\pi_k)=\sum_{\mathbf{t} \in \mathcal{T}_\epsilon(\mathbf{t}^{obs})}\mathbbm{I}\{\pi_i(\mathbf{t})=\pi_k\}\cdot \Pr{}_{\mathbf{T}|\mathcal{T}_\epsilon(\mathbf{t}^{obs}) }(\mathbf{T}=\mathbf{t}).$ Under Assumption (ref), these inclusion probabilities are strictly positive for all \( i \in [N] \).
Applying the selection rules outlined above, we introduce a novel randomization testing procedure for \( H_0^\tau \), summarized in Procedure (ref). This procedure serves to illustrate the unique challenges associated with the proposed conditioning framework and provides a foundation for the methodological developments that follow.
\RestyleAlgo{ruled} \SetKwComment{Comment}{/* }{ */}
Several features of Procedure (ref) depart deliberately from conventional CRI methods, giving rise to nontrivial technical challenges that must be addressed to ensure control of the Type I error rate. First, the set of focal units varies across focal assignments. Thus, even under a true null hypothesis, the distribution of the observed test statistic computed with the super-focal units may differ from that of the imputed test statistics. Second, under Assumption (ref), the observed test statistic is well-defined under all assignments in the design. By contrast, the imputed test statistics are computable only over the subset of focal assignments, which are not independent and identical (i.i.d) draws from the full assignment space. This may introduce a systematic discrepancy between the observed and randomized reference distributions. Third, the proposed procedure does not induce a partition of \( \mathcal{T}_0 \) into mutually exclusive focal assignment sets. In particular, a single assignment may belong to multiple focal assignment sets. As shown by hennessy2016conditional and zhang2023randomization, such a partitioning property is central to establishing the unconditional validity of CRI procedures via conditional arguments.
We address the first two challenges by employing inverse probability weighted (IPW) estimators for both test statistics and p-values horvitz1952generalization, hajek1971comment. First, observe that super-focal and focal units are sampled without replacement from the finite population, with unequal inclusion probabilities given by $\Pr{}_{\mathbf{T}}(\pi_i(\mathbf{T}^{obs})=\pi_k)$ and $\Pr{}_{\mathbf{T}}(\pi_i(\mathbf{T}^{obs})=\pi_k)\cdot \Pr{}_{\mathbf{T}|\mathcal{T}_\epsilon(\mathbf{t}^{obs}) }(\pi_i(\mathbf{T})=\pi_k)$ respectively, for each $i \in [N].$ Incorporating these probabilities yields consistent (and unbiased) estimators of population parameters. Similarly, focal assignments are drawn without replacement from the finite population of assignments induced by the experimental design, again with unequal inclusion probabilities \( \phi_{\mathbf{t}}^\epsilon \), defined in (ref), for all \( \mathbf{t} \in \mathcal{T}_0 \). Hence, IPW estimators can be used to construct consistent (and unbiased) estimators of the p-value.
To address the third challenge, we propose test statistics that ensure unconditional validity directly without relying on conditional validity as an intermediate step. Proposition (ref) provides a general sufficient condition for finite-sample unconditional validity. This condition takes the form of a pairwise pointwise dominance requirement and does not require the standard exchangeability between observed and imputed test statistics typically invoked to justify the validity of randomization tests.\footnote{Following the release of an earlier version of this paper on arXiv, zhong2024unconditional pursued this direction by proposing test statistics specifically constructed to satisfy the exchangeability condition required for finite-sample validity. However, to attain a test of size $\alpha\in (0,1)$, his procedure requires a target nominal level of $\alpha/2$. This not only complicates interpretation but also undermines the conventional theoretical justification for level-$\alpha$ tests, as the procedure no longer guarantees control at the nominal level specified by the researcher. }
The Appendix contains proofs of all main results; all remaining proofs are given in the Supplementary Material.
Proposition (ref) establishes that Procedure (ref), when implemented with unbiased estimators of test statistics and p-values, achieves finite-sample unconditional validity under the null hypothesis. The key requirement is that the conditional null distribution of the imputed test statistic must not exceed the observed distribution at any realized value of the observed test statistic, as expressed in inequality (ref). This condition is strictly weaker than pairwise first-order stochastic dominance: while first-order dominance implies (ref), the converse does not hold. Indeed, (ref) may be satisfied even when the null and observed distributions cross, allowing for some forms of second-order dominance. Thus, Proposition (ref) provides a flexible criterion that broadens the class of valid testing procedures beyond those justified by stronger distributional assumptions.
Identifying test statistics that satisfy condition (ref) in finite samples is generally difficult, as the relevant null conditional probabilities typically lack closed-form representations. To circumvent this challenge, we focus on constructing test statistics for which condition (ref) holds asymptotically. Thus, we aim to establish validity in the limit as $N \to \infty$, in the spirit of wu2021randomization. Our asymptotic framework follows the classical sequence-of-populations approach initiated by brewer1979class and extended by isaki1982survey. Specifically, we consider a sequence of nested finite populations indexed by size \( N \), where the treatment assignment mechanism is independently reapplied to each population in the sequence. The following theorem shows that for a broad class of test statistics, condition (ref) is satisfied in the limit as \( N \to \infty \).
Theorem (ref) asserts that any consistent non-negative test statistic that converges in probability under the null to the minimum value in its support satisfies the pairwise dominance condition in the limit as \( N \to \infty \). Thus, for this class of test statistics, the randomization test described in Procedure (ref) is unconditionally valid at any significant level $\alpha\in(0,1)$ as $N\to\infty.$
Next, we introduce a set of primitive statistics that serve as building blocks for the construction of test statistics satisfying the sufficient conditions of Theorem (ref). Let \( \bar{y}_{t}(\pi_k) \) denote the Horvitz--Thompson (HT) estimator of the mean potential outcome under effective treatment \( (t, \pi_k) \), defined as $\bar{y}_{t}(\pi_k):=N^{-1}\sum_{i=1}^NY_i\cdot D_i/\Pr{}_{\mathbf{T}}(D_i=1)$ where $D_i= \mathbbm{I}\{T_i=t,\Pi_i=\pi_k\}.$ Also, let $\hat{\sigma}^2_{t}(\pi_k)$ denote the Yates-Grundy consistent and unbiased estimator of the population variance of the potential outcome under effective treatment $(t, \pi_k),$ i.e., $\hat{\sigma}^2_{t}(\pi_k):=(N(N-1))^{-1}\sum_{i=1}^N\sum_{j>i}D_iD_j (Y_i-Y_j)^2/\Pr{}_{\mathbf{T}}(D_i=1, D_j=1)$ yates1953selection. Finally, let $\hat{F}_{t,\pi_k}(\cdot)$ denote the HT estimator of the empirical distribution of the potential outcome under effective treatment $(t, \pi_k),$ i.e., $\hat{F}_{t,\pi_k}(y):=N^{-1}\sum_{i=1}^N D_i\mathbbm{I}\{Y_i\leq y\}/\Pr{}_{\mathbf{T}}(D_i=1).$
We consider three feasible statistics to test for $H_0^\tau.$ Specifically, we let $$z_{_{AR}}(\mathbf{Y}^{obs}, \mathbf{T}^{obs} , \boldsymbol{\pi}(\mathbf{T}^{obs});\mathcal{S}(\mathbf{T}^{obs}) ):=\left|\frac{\hat{\sigma}^2_1(\pi_k)- \hat{\sigma}^2_0(\pi_k) }{ \hat{\sigma}^2_0(\pi_k)}\right|$$ denote the “absolute ratio of variance”(AR) statistic, $$z_{_{MR}}(\mathbf{Y}^{obs}, \mathbf{T}^{obs} , \boldsymbol{\pi}(\mathbf{T}^{obs});\mathcal{S}(\mathbf{T}^{obs}) ):=\max\Bigg\{\frac{\hat{\sigma}^2_1(\pi_k)}{ \hat{\sigma}^2_0(\pi_k)}, \frac{\hat{\sigma}^2_0(\pi_k)}{ \hat{\sigma}^2_1(\pi_k)} \Bigg\}$$ denote the “maximum ratio”(MR) statistic, and $$z_{_{SK}}(\mathbf{Y}^{obs}, \mathbf{T}^{obs} , \boldsymbol{\pi}(\mathbf{T}^{obs});\mathcal{S}(\mathbf{T}^{obs}) ):=\max_{y}|\hat{F}_{0,\pi_k}(y)-\hat{F}_{1,\pi_k}(y+\tau)|$$ denote the shifted Kolmogorov--Smirnov statistic.
To analytically justify that the foregoing statistics satisfy the sufficient conditions of Theorem (ref)---particularly the convergence in probability condition---we introduce the following regularity conditions.
Assumption (ref) imposes standard boundedness conditions on potential outcomes and exposure probabilities. Assumption (ref) restricts the extent of dependence induced by the design and the exposure mapping across quadruples of units. In particular, it ensures that, although individual pairs of units may exhibit nontrivial clustering in their exposures, the joint dependence between distinct pairs becomes asymptotically negligible as the population size grows. This condition is less restrictive than the local dependence condition required for standard asymptotic-based inference under network interference; see aronow2017estimating, liu2014large, leung2020treatment and leung2022causal.
Theorem (ref) identifies three test statistics\footnote{The list is not exhaustive; additional statistics satisfy the sufficient conditions of Theorem (ref) under \( H_0^\tau \) and related null hypotheses. For instance, under no treatment effect nulls, the absolute difference-in-means satisfies these conditions.} for $H_0^\tau$ that guarantees asymptotic validity of Procedure (ref) at all nominal levels. We numerically assess the finite sample performance of these test statistics in Section (ref).
The class of test statistics that satisfy the sufficient conditions in Theorem (ref) may nonetheless be restrictive. In the next theoretical result, we establish the asymptotic validity of Procedure (ref) for a broader class of test statistics. Specifically, we show that validity holds asymptotically over a restricted set of significance levels, which may depend on the choice of statistic. While this result does not guarantee uniform validity over the entire unit interval, it remains practically relevant since standard nominal levels used in empirical applications typically lie in the interval $(0, 0.1]$.
Corollary (ref) states that if the test statistic is consistent and its probability limit lies strictly below the supremum of its support, then Procedure (ref) is asymptotically valid for nominal levels approximately between zero and the limiting survival function evaluated at the point of convergence. For example, if the asymptotic distribution is symmetric with a mean value equal to the limit \( z^* \), then the procedure is valid for all significance levels \( \alpha \) in the interval \( (0, 0.5] \). Corollary (ref) generalizes Theorem (ref): when \( z^* = \inf I \) (the infimum of the support), we recover the earlier result; when \( z^* = \sup I \), the procedure fails to control size for any \( \alpha \in (0,1) \).
Most appropriately scaled test statistics satisfy the sufficient conditions stated in Corollary (ref). In contrast to Theorem (ref), however, verifying these conditions analytically---particularly the range of significance levels for which asymptotic validity holds---requires deriving the limiting distribution of the test statistic. This task typically demands stronger assumptions on the network structure, such as the widely used local dependence condition. Similar to the result in Theorem (ref), under local dependence and additional regularity conditions, IPW versions of classical statistics for testing equality of variances---such as the variance ratio, Pitman--Morgan pitman1939note, morgan1939test, and Levene’s test levene1960robust---can be shown to satisfy the conditions in Corollary (ref).
It is worth emphasizing that, unlike the classical Fisher randomization test---where the p-value distribution stochastically dominates the uniform distribution solely due to the discreteness of the assignment mechanism---the test statistics that satisfies the sufficient conditions of in Theorem (ref) and Corollary (ref) may yield even greater divergence from uniformity under the null. As a result, Procedure (ref) is valid in a conservative sense: the Type I error rate may fall strictly below the nominal level $\alpha$. This is akin to the validity results of weak nulls in wu2021randomization.
In this subsection, we characterize the asymptotic size distortion of Procedure (ref). Our objective is to derive informative bounds on the limiting size distortion associated with various test statistics. These bounds, expressed in terms of observable data features, serve to guide practitioners in identifying settings where the proposed procedure performs optimally or exhibits diminished control over Type I errors. We begin by introducing the following definition.
We introduce additional notation to formalize the size distortion of Procedure (ref). Let $ F_{Z|\mathcal{T}_\epsilon(\mathbf{t}^{obs})}(z^{obs}) = \Pr{}_{\mathbf{T}|\mathcal{T}_\epsilon(\mathbf{t}^{obs})}( z\big( \mathbf{Y}(\mathbf{T}), \mathbf{T}, \boldsymbol{\pi}(\mathbf{T}); \mathcal{U}(\mathbf{t}^{obs}, \mathbf{T}) \big) \leq z^{obs} \,\big|\, \mathcal{C}(\mathbf{t}^{obs}), H_0^{\tau} ) $ denote the conditional distribution of the imputed test statistic evaluated at an arbitrary observed value \( z^{obs} \), under the null hypothesis \( H_0^{\tau} \). Similarly, define $ F_{Z^{obs}}(z^{obs}) = \Pr{}_{\mathbf{T}^{obs}}( z( \mathbf{Y}^{obs},\mathbf{T}^{obs}, \boldsymbol{\pi}(\mathbf{T}^{obs}); \mathcal{S}(\mathbf{T}^{obs})) \leq z^{obs} \,\big|\, H_0^{\tau} ) $ as the marginal distribution of the observed test statistic under \( H_0^{\tau} \), taken over the assignment mechanism. Finally, define the mixture null distribution by averaging the conditional distributions of imputed test statistics across all possible treatment assignments: $ F_{\mathrm{mix}}(z^{obs}) = \sum_{\mathbf{t}^{obs} \in \mathcal{T}_0} F_{Z|\mathcal{T}_\epsilon(\mathbf{t}^{obs})}(z^{obs})\cdot \Pr{}_{\mathbf{T}^{obs}}(\mathbf{T}^{obs}=\mathbf{t}^{obs}).$
In the presence of multiple imputed distributions, a key challenge is to identify an appropriate reference distribution against which the distribution of the observed test statistic can be meaningfully compared. For test statistics that satisfy the pairwise dominance condition, the mixture distribution \( F_{\text{mix}}(\cdot) \) serves as a suitable reference distribution for this purpose. See Section 4.1 of the Supplementary Material for numerical justification. The degree of size distortion of the proposed procedure can then be measured by the discrepancy between \( F_{\text{mix}}(\cdot) \) and \( F_{Z^{\text{obs}}}(\cdot) \); the smaller this discrepancy, the better the size control.
The next theorem provides a general bound on the size distortion of the proposed procedure, using the mixture randomization distribution, \( F_{\mathrm{mix}}(\cdot), \) as the benchmark distribution for the imputed test statistics.
Theorem (ref) establishes a general uniform bound on the asymptotic size distortion of the proposed testing procedure in Procedure (ref). The bound consists of two components. The first term in (ref) reflects the conventional size distortion of randomization tests arising from the discreteness of the treatment assignment. The second term captures the contribution of the pairwise dominance condition satisfied by the test statistics under the sufficient conditions of Theorem (ref). This term is determined by the limiting distribution of the test statistic and the Wasserstein distance between the mixture randomization distribution and the limiting distribution of the observed test statistic as \( N \to \infty \).
If the observed test statistic converges weakly to a non-degenerate distribution, such as the normal, gamma, or beta distribution, then Stein's method stein1972bound can be employed to derive an informative upper bound on the second term, $\lim_{N \to \infty} \{2\kappa \cdot d_W(F_{\mathrm{mix}}, F_{Z^{obs}})\}^{1/2}$. In the following corollary, we apply this technique to the case where the test statistic is asymptotically standard normal. The resulting bound depends on observable features of the data, thereby offering practical insight into the conditions under which the proposed testing procedure achieves better size control.
Corollary (ref) states that for test statistics expressible as sums of random variables and converging under the null to the standard normal distribution, the asymptotic size distortion depends on the distribution of focal units through $ W_i$. In particular, the distortion decreases as the number of focal units increases. Furthermore, the size distortion worsens as network density increases, as captured by $A_{\mathrm{max}},$ the maximum degree in the network.
In many applications, the functional form of $\tau(t, t', \pi, \pi', X_i)$ in $H_0^G$ is unknown a priori. A natural approach to testing CTEs is to posit $\tau(\cdot)$ as the conditional average treatment effect (CATE) function. However, this introduces a nuisance parameter into the testing problem, as the CATE function is typically unknown and must be estimated from the data.
A naive strategy involves replacing the CATE function with its sample analog and proceeding with the proposed randomization inference procedure. However, as noted by ding2016randomization, such plug-in approaches typically lack theoretical guarantees for valid inference, particularly regarding the control of Type I error in both finite samples and asymptotic regimes.
We propose a tractable and easily implementable procedure that jointly addresses the two sources of non-sharpness: the multiplicity of potential outcomes and the presence of nuisance parameters. The procedure combines the conditioning strategy outlined in Procedure (ref) with the confidence interval (CI) method of berger1994p, as adapted to the randomization testing framework by ding2016randomization.
The core idea of the CI method is to compute the p-value corresponding to the least favorable value of the nuisance parameter within a prespecified confidence region, thereby ensuring uniformly valid inference over all values in the region.
Before presenting the validity results associated with applying the CI method to the proposed randomization testing procedure---referred to hereafter as the CRI-CI procedure---we formally introduce the method in the context of $H_0^{\tau}$, where $\tau = \text{ATE}(\pi_k)$ denotes the average direct treatment effect at exposure value $\pi_k$.
Let $\mathrm{CI}_\gamma$ denote the $(1-\gamma)$ confidence interval for the nuisance parameter $\tau,$ where $\gamma \in (0,1).$ The p-value under the CRI-CI procedure is then defined as:
where $pval_k( \cdot , \tau')$ denotes the p-value when $\tau =\tau'.$\footnote{In practice, computing p-values for every possible value of \( \tau \) within the confidence interval is computationally infeasible. Following ding2016randomization, we implement the resulting procedure by evaluating p-values over a finite uniform grid within the estimated confidence interval. This approximation preserves the theoretical validity of the method. } Procedure (ref) summarizes the CRI-CI procedure for testing $H_0^{\tau}$ in the presence of the nuisance parameter $\tau.$
The resulting p-values from the CRI-CI procedure represent worst-case values over the confidence region and are therefore conservative by construction ding2016randomization. Developing methods that mitigate the associated power loss remains an open question and is beyond the scope of this paper. Nonetheless, Theorem (ref) establishes that the CRI-CI procedure is asymptotically valid.
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Implementation Guidelines}
In this section, we provide practical guidelines for implementing the testing procedures developed in the preceding sections using Monte Carlo methods. In particular, we focus on estimating test statistics and p-values via the Monte Carlo approach introduced by dwass1957modified. For a comprehensive review of Monte Carlo p-value estimation and its applications in econometrics, see dufour2001monte. To save space, we relegate other implementation issues to Section 2 of the Supplementary Material.
Computing exact p-values is often infeasible due to the large number of possible focal assignments, even in relatively small populations. This challenge is exacerbated by the fact that unbiased and consistent estimation of the test statistics and p-values in Procedures (ref) and (ref) requires knowledge of the inclusion probabilities associated with focal units, super-focal units, and focal assignments.
We consider two estimation strategies: (i) IPW estimators, which are consistent when the relevant inclusion probabilities are known, and (ii) uniformly weighted estimators, which are computationally more tractable but generally biased and inconsistent. We analytically characterize the magnitude and direction of this bias and justify the use of uniformly weighted estimators in settings where the computational burden of IPW estimation is prohibitive.
\paragraph*{Inverse Probability Weighted Estimators} Based on procedures (ref) and (ref), recall that super-focal units, focal units and assignments are sampled without replacement from their respective finite populations, typically with unequal probabilities. As discussed in Section (ref), IPW estimators remain consistent and unbiased for both the observed and imputed test statistics, as well as for the resulting p-values, under this sampling design. However, computing the exact inclusion probabilities is only feasible when the sample size is small. We propose the following Monte-Carlo estimation procedure that involves:
The estimators of the inclusion probabilities mentioned above are consistent; see fattorini2006applying and aronow2017estimating for formal proofs.
Using the estimated inclusion probabilities of units in steps 2 and 4, we can compute the unbiased IPW estimators of observed and imputed test statistics. In addition, using the inclusion probability of the randomly drawn assignments in step 5, we can compute the unbiased Horvitz--Thompson estimator of the p-value defined as
or the more efficient sarndal2003model, biased but consistent H\'ajek ratio estimator
where $ \hat{S}^{HT}=\sum_{\mathbf{t} \in \widehat{\mathcal{T}}_\epsilon(\mathbf{t}^{obs})}1/\hat{\phi}_\mathbf{t}^\epsilon$ is the Horvitz--Thompson estimator of the $ |\widehat{\mathcal{T}}_\epsilon(\mathbf{t}^{obs})| .$
If the focal assignment set, $\widehat{\mathcal{T}}_\epsilon(\mathbf{t}^{obs}),$ is large, a more computationally efficient alternative for computing the p-value is to use a smaller random sample of $B$ i.i.d. draws---with $B < |\widehat{\mathcal{T}}_\epsilon(\mathbf{t}^{obs})|$---from the focal assignment set to approximate the null distribution. We defer the discussion of the choice of $B$ to Section 2 of the Supplementary Material.
\paragraph*{Uniformly Weighted Estimators} The computational demands of the IPW estimators can be prohibitive. This motivates the consideration of uniformly weighted estimators of test statistics and p-values, which treat the super-focal units, focal units, and focal assignments as if they were sampled with replacement with equal probability. These estimators are attractive for their computational simplicity, as they do not require knowledge of inclusion probabilities.
It is well known in the survey sampling literature (see, e.g., sarndal2003model) that uniformly weighted estimators are generally biased and inconsistent when sampling is conducted without replacement and with unequal inclusion probabilities. Nevertheless, in some settings and for certain statistics, the bias introduced by this approximation may be sufficiently small to be inconsequential in practice. To clarify when this holds in the context of variance-based statistics, we formally characterize the bias of the uniformly weighted sample variance estimator in the following proposition.
Proposition (ref) asserts that the bias of the uniformly weighted sample variance estimator when computed from a non-i.i.d. sample is governed by two principal factors. First, the bias decreases as the variability in inclusion probabilities \( p_i \) declines---that is, as the sampling design approaches uniformity---since the deviations \( p_i - n_{tk}/N \) and \( p_{ij} - n_{tk}(n_{tk}-1)/N(N-1) \) become smaller. Second, because these deviations are mean-zero (zero-sum) over the population, the bias is further mitigated when inclusion probabilities are uncorrelated with the squared outcomes and cross-products of outcomes. In particular, the bias vanishes when the inclusion probabilities are independent of \( Y_i(t, \pi_k) \), holding the sampling fraction fixed.
These observations imply that, in network settings where the degree distribution is approximately constant across units, the bias of the uniformly weighted sample variance estimator is negligible. Moreover, in randomized experiments, treatment assignments are independent of potential outcomes by design, further mitigating the bias introduced by uniform weighting. Consequently, under such conditions, sample variance-based test statistics and their associated p-values can be computed using uniform weights at minimal accuracy cost, providing a computationally efficient alternative to the more intensive IPW estimators.
Similar to the IPW estimators, uniformly weighted p-values are also computationally infeasible if the focal assignment set is large. In practice, however, one can use the Monte Carlo method to estimate the p-values in Procedures (ref)--(ref) by generating i.i.d. draws from the focal assignment set using rejection or importance sampling methods as described in branson2019randomization. Specifically, for uniformly weighted p-values, we recommend using a random sample of focal assignments of size $B < |\mathcal{T}_\epsilon(\mathbf{t}^{obs})|$---that meet the eligibility requirements of the focal assignment set to compute the p-values. The resulting Monte-Carlo p-value is of the form specified in lehmann2022general:
Under the assumption that focal units and assignments are approximately i.i.d draws, the approximated p-value is “approximately” consistent.\footnote{hennessy2016conditional and hoshino2023randomization show that when focal units and assignments are i.i.d draws, then the approximated p-value is a consistent estimator of the true p-value.}
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Simulation} In this section, we present the design and results of Monte Carlo experiments that evaluate the finite-sample performance of the proposed testing procedures for $H_0^\tau$. Specifically, we compare the performance of several feasible test statistics in terms of both empirical size and power. In addition, we benchmark the proposed method against the biclique method.
Our simulation design builds on that in ding2016randomization, with key modifications to incorporate network interference. For each unit $i\in [N],$ the potential outcomes and treatment effects are defined as:
where $\psi_0,$ $\psi_1$ and $\sigma_\tau$ parameterize different forms of treatment effect heterogeneity. Throughout, we set $\psi_0=\psi_1=0$ and $N=200.$
Treatment is assigned in two stages. First, we generate assignments under complete randomization with $N_0=N_1=100.$ Second, we restrict attention to assignments satisfying the overlap condition in Assumption (ref). To estimate the inclusion probabilities, we randomly draw 5000 treatment assignments and follow the steps outlined in Section (ref).
Throughout, we use a binary-valued exposure mapping defined as:\footnote{Results for a multi-valued exposure mapping defined as the number of treated neighbors are provided in Section 4.2 of the Supplementary Material.} $$\pi_i(\mathbf{T}):= \mathbbm{I}(\sum_{j=1}^NT_jA_{ij}/\sum_{j=1}^NA_{ij} >0.5).$$ Based on the exposure mapping and the population size, we use the intervals $\mathcal{I}_{t,\epsilon}=\mathcal{I}_{t,20} =[20,\,\, |\mathcal{S}(\mathbf{t}^{obs})|-20]$ for $t=0,1$.\footnote{In Section 2.1 of the Supplementary Material, we numerically examine how variations in $\mathcal{I}_{t,\epsilon}$ influence the size and power of the test.} Finally, we draw the errors $\{\varepsilon_i\}_{i\in[N]}$ from a multivariate normal distribution whose covariance structure reflects both individual variability and network-based correlation.
All rejection rates are computed as the proportion of rejections at the 5% significance level based on 1,000 Monte Carlo replications, yielding a simulation standard error of approximately 0.00689 under the null.
In this subsection, we assess the finite sample performance of some feasible statistics for testing $H_0^\tau$ using Procedure (ref). Specifically, we compute the empirical rejection probabilities using the absolute ratio of variance, maximum ratio, shifted Kolmogorov--Smirnov, Levene, and Pitman--Morgan statistics. Recall that the absolute ratio of variance, maximum ratio, and shifted Kolmogorov--Smirnov statistics satisfy the sufficient conditions of Theorem (ref). On the other hand, the Levene and Pitman--Morgan statistics satisfy the sufficient conditions of Corollary (ref). See Section 3 of the Supplementary Material for a formal discussion of the Pitman--Morgan statistic, a description of how it is adapted to the present setting, and its limitations.
For the results reported in this subsection, we utilize an adjacency matrix where each unit is connected to at least one other unit. In addition, $\Psi_i(\pi)'$s are i.i.d draws from a right-skewed Beta distribution with shape parameters $0.5$ and $2+\pi$ for all $\pi \in \boldsymbol{\Pi}.$
In Figure (ref)(a), we compare the performance of the absolute ratio of variance (ratio) and maximum ratio test statistics. Specifically, we report the empirical rejection probabilities by varying $\sigma_\tau$ over the range $ \{0.00, 0.01, 0.02, \dots, 1.00 \}.$ As expected, the result shows that the statistics tend to under-reject when the null is true. The statistics based on the H\'ajek estimator tend to exhibit lower rejection rates compared to their uniformly weighted counterparts, although the differences are negligible. The difference is negligible since the exposure mapping is binary, and the variability in inclusion probabilities is small. This corroborates the theoretical findings in Section (ref).
In Figure (ref)(b), we compare the rejection rates of all the statistics computed using uniformly weighted estimators. Our results show that the Pitman--Morgan statistics tend to have the lowest size distortion and also outperform the rest in terms of power. This is unsurprising since the Pitman--Morgan statistics is related to the likelihood ratio criterion (see morgan1939test); as such, it satisfies the Neyman--Pearson Lemma neyman1933ix. In both figures, p-values are estimated using the H\'ajek estimator.
In this subsection, we compare the finite sample performance of four procedures: Biclique-Oracle (B-Oracle), super-focal-Oracle (SF-Oracle), Plug-in, and Confidence Interval (CRI-CI).\footnote{For CRI-CI, we set $\gamma=0.0001$ and confidence interval grid size of 151.} The first two procedures rely on knowledge of the nuisance parameter and are infeasible in many practical settings. B-Oracle uses the biclique-based conditioning procedure of puelz2022graph, while SF-Oracle is based on Procedure (ref). On the other hand, the Plug-in method estimates the ATE and substitutes it directly into the null hypothesis, while the CRI-CI method is based on Procedure (ref).
Throughout this subsection, we employ an adjacency matrix where each unit is connected to at most five other units to mimic the observed network structure of our empirical application data in the next section. In addition, we consider two distributions of $\Psi_i(\pi)'$s: (i) i.i.d draws from a “near-symmetric” Beta distribution with shape parameters $10$ and $10+\pi$ for all $\pi \in \boldsymbol{\Pi},$ and (ii) i.i.d draws from a right-skewed Beta distribution with shape parameters $0.5$ and $2+\pi$ for all $\pi \in \boldsymbol{\Pi}.$ Finally, we use the H\'ajek absolute ratio test statistic.
Due to the high computational time requirement of the B-Oracle procedure, we only compare its empirical size to that of SF-Oracle. We set the number of randomizations in both procedures to $ 50$.
The results in Table (ref) show that the Biclique-Oracle procedure overrejects under the null. We conjecture that this is primarily attributable to the small size of the bicliques generated by the algorithm, which limits the amount of information used for inference and thus reduces power. Across 1000 replications, we record an average biclique size (the number of focal units and assignments) of 6.944. In contrast, the SF-Oracle procedure uses 50 focal assignments to estimate the reference null distribution in every iteration. Thus, as expected, the proposed SF-Oracle procedure is valid. In Section 8.1 of the Supplementary Material, we present a formal analytical comparison of the power properties of the B-Oracle and the SF-Oracle.
Next, we extend our simulation exercise by computing the empirical rejection probabilities of the SF-Oracle, Plug-in, and CRI-CI procedures by varying $\sigma_\tau$ over the range $ \{0.00, 0.01, 0.02, \dots, 1.00 \}.$ For all the results reported below, we set $B=199.$
Using symmetrically distributed untreated potential outcomes, the power curves in Figure (ref)(a) show that the Plugin, SF-Oracle, and CRI-CI procedures are all valid, with the CRI-CI procedure uniformly exhibiting lower power as expected. Notably, the Plugin and the \textit{SF-Oracle} procedures share an identical power function. This equivalence aligns with findings in ding2016randomization for symmetric distributions. The \textit{CRI-CI} procedure has low statistical power for parameter values close to the null parameter value $(\sigma_\tau=0)$, corroborating our theory. In Figure (ref)(b), we display the power curves under the \textit{asymmetrically distributed untreated potential outcomes}. In general, the results are similar to those in Figure (ref)(a).
Thus, among procedures that guarantee theoretical control of Type I error, the proposed SF-Oracle procedure is most suitable when all parameters are known, whereas the CRI-CI procedure is appropriate in the presence of nuisance parameters.
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Empirical Application} In this section, we illustrate the application of the proposed testing procedures using real data from the field experiment conducted by cai2015social. The study was designed to assess whether farmers’ understanding of a weather insurance policy influences their decision to purchase the product. Specifically, the authors evaluate the effects of two types of information sessions on insurance adoption among 4,902 households residing in 173 small rice-producing villages---nested within 47 administrative villages---across three regions in Jiangxi Province, China. Their findings indicate that the format of the information session not only has a direct impact on participants’ adoption behavior but also exerts significant peer effects on the adoption decisions of their named friends.
The dataset includes detailed social network information at the household level, as well as a range of pre-treatment covariates, including age, gender, rice cultivation area, risk aversion score, and the fraction of household income derived from rice production. The outcome of interest is binary, indicating whether the household purchased the insurance policy or not. The resulting social network is sparse because households could nominate up to five friends.
In each village, the experimental design consisted of two rounds of information sessions introducing the insurance product. During each round, two sessions were conducted simultaneously: one providing basic information (the simple session) and the other offering more detailed information (the intensive session). The second round of sessions was administered three days after the first. This delay was sufficiently long to allow participants to share information with their direct friends but not enough to fully diffuse shared information throughout the broader network via friends-of-friends; see cai2015social. Consequently, Assumption (ref) is plausibly satisfied.
To illustrate the proposed CRI-CI testing procedure, we restrict attention to three of the largest villages in the study---Dukou, Yazhou, and Yongfeng---comprising 502 households. Because the overall social network is clustered, we can analyze a subset of villages and their corresponding sub-networks in isolation.
A specification test of the interference structure by hoshino2023randomization suggests that the correct network exposure mapping for this dataset is a threshold function of the number of neighbors who attended the first-round intensive sessions. This is formally defined as $ \pi_i(\mathbf{T}):= \mathbbm{I}(\sum_{j=1}^NT_jP_jA_{ij}/ \sum_{j=1}^NP_jA_{ij}\geq 0), $ where for $j\in [N],$ $T_j$ and $P_j$ are the treatment indicators: $T_j = 1$ if household $j$ attended the intensive session (and zero otherwise), and $P_j = 1$ if household $j$ attended the first-round sessions (and zero otherwise).\footnote{We do not account for selection in the test, as post-selection inference is beyond the scope of this paper.}
We test three null hypotheses: (i) $H_0,$ (ii) $H^{\Pi}_{0}: Y_i(1, \pi) - Y_i(0, \pi)= \tau (\pi) \,\,\text{for some function}\,\, $\\$ \tau(\cdot), \,\, \forall\,\, \pi \in \mathbf{\Pi},\,\, and\,\,\forall \,\,i\in [N],$ and (ii) $H^{X,\Pi}_{0}: Y_{i}(1, \pi) - Y_{i}(0, \pi)= \tau(\pi, X_i) \, for some \, \tau(\cdot, \cdot),\, $\\$\forall\pi \in \mathbf{\Pi},\,\forall\, X_i \in \mathbbm{X}\,and \,\forall\,i\in [N];$ see Section 7 of the Supplementary Material for discussions of these nulls. For $H^{X,\Pi}_{0}$, we use the binary covariate insurance_repay, which equals one if a household previously received a payout from an insurance policy and zero otherwise.
As described in cai2015social, the treatment pair $(T, P)$ was assigned to households using a stratified randomization design using household size and area of rice production per capita. We treat $T_i$ as the sole source of randomness, conditioning on the observed values of $P_i,$ which are held fixed and included as covariates for all $i\in [N].$ In particular, we define the null hypotheses with respect to units for which $P_i=0$, i.e., we focus on the subset of households that attended the second-round sessions.
We use the focal assignment set defined as $\mathcal{T}_\epsilon(\mathbf{t}^{obs}):= \{\mathbf{t}' \in \mathcal{T}_0,\,\, \epsilon \leq \hat{R}(0,\mathbf{t}';\mathcal{S}(\mathbf{t}^{obs}))\leq|\mathcal{S}(\mathbf{t}^{obs})|-\epsilon,\,\, \text{and}\,\, \epsilon \leq \hat{R}(0,\mathbf{t}';\mathcal{S}(\mathbf{t}^{obs}))\leq|\mathcal{S}(\mathbf{t}^{obs})|-\epsilon \},$ where $\epsilon$ is varied between 20 to 50, for robustness check.
Tables (ref)–(ref) report p-values for tests of the null hypotheses evaluated at each exposure value, using the H\'ajek absolute ratio of variance statistic and the CRI-CI procedure with $\gamma=0.0001, B=399$ and confidence interval grid size of 151. To conduct joint inference across exposure values, one may apply the standard multiple-testing procedures (MTPs) that control either the family-wise error rate or the false discovery rate; see romano2010multiple for an overview of MTPs.
Based on the p-values in Table (ref), we fail to reject the null hypothesis of constant treatment effect across the population. As a result, there may be no heterogeneous effect on the decision to purchase weather insurance among participants in the second round. Based on the decision from the test of $H^{\tau}_0,$ we must also fail to reject $H^{\Pi}_0$ and $H^{X,\Pi}_0.$ The p-values in Tables (ref)--(ref) corroborate this assertion.
\ifnum0=1 \addtocounter{alphasect}{1} \fi \oldsection{Conclusion} This paper develops randomization-based inference procedures for testing heterogeneous treatment effects in the presence of network interference. Within the exposure mapping framework, we formulate a general class of non-sharp null hypotheses that encompass various notions of constant treatment effects in networked populations. These hypotheses depend on unknown functions that may act as nuisance parameters, and the restrictions they impose do not permit full imputation of unobserved potential outcomes. Existing conditional randomization procedures either lack statistical power, are invalid, or are inapplicable in this setting.
We propose a novel randomization testing procedure that constructs a data-dependent focal assignment set tailored to the observed treatment assignment and exposure configuration. In contrast to conventional approaches that fix focal units and assignment sets ex-ante, our method allows the set of focal units to vary across focal assignments in accordance with the restrictions imposed by the null hypothesis. These adaptive features introduce technical complications that render the use of uniformly weighted estimators of test statistics and p-values invalid. To overcome this, we propose consistent and unbiased inverse probability-weighted estimators. Under general conditions on the test statistic, we establish the asymptotic validity of the procedure and characterize the limiting size distortion in terms of observable quantities.
We illustrate the proposed procedure using data from the field experiment of cai2015social, which investigates how information influences the adoption of weather insurance among rice farmers in rural China. Finally, we present results from an extensive Monte Carlo study that corroborate the theoretical findings.
The randomization testing procedure developed in this paper is broadly applicable to a wide class of partial null hypotheses, including those arising in settings without interference. Its flexibility in accommodating non-sharp nulls while preserving validity makes it a promising foundation for addressing more challenging testing problems, such as Neyman's weak null hypothesis neyman1923application. An important direction for future research is the development of refined procedures that mitigate the finite-sample size distortion observed under the proposed approach, potentially through improved conditioning schemes or alternative test statistics.