Imputation-based randomization tests for randomized experiments with interference

Consider a finite population ${\mathbb{U}}=\{1,\ldots,N\}$ with $N$ units, to each of which $K$ different treatments in the treatment set ${\mathcal{T}}=\{t_{1},\ldots,t_{K}\}$ can be potentially assigned. For simplicity, we mainly consider the simple case with $K=2$ in this study, though our approach can be easily extended to more complicated scenarios where $K>2$ (as demonstrated in the real data application in Section 5). Let ${\mathcal{T}}=\{0,1\}$ be the binary treatment set with $0$ as control treatment and 1 as active treatment, and $Z_{i}\in{\mathcal{T}}$ be the indicator of treatment assignment for the $i$-th unit in ${\mathbb{U}}$. The collective treatment assignment for all $N$ units is represented by the binary vector ${\boldsymbol{Z}}=(Z_{1},\ldots,Z_{N})\in{\mathcal{Z}}\subseteq\{0,1\}^{N}$. Let ${\boldsymbol{Z}}^{{\rm{obs}}}=(Z_{1}^{\rm{obs}},\ldots,Z_{N}^{\rm{obs}})\in{\mathcal{Z}}$ be the treatment assignment actually implemented in the randomized experiment. Throughout this study, we always assume that ${\boldsymbol{Z}}^{\rm{obs}}$ is a random sample from a known randomized design ${\boldsymbol{F}}_{\mathcal{Z}}$ defined over ${\mathcal{Z}}$.

In the presence of interference, each unit $i$ has a collection of potential outcomes, i.e., $\{Y_{i}({\boldsymbol{z}})\}_{{\boldsymbol{z}}\in{\mathcal{Z}}}$, in general, because there could be a total of $2^{N}$ distinct exposures for each unit in the most general case where interference occurs between any two units. In practice, however, interference typically occurs only among certain units. In this case, the number of distinct exposures for each unit would be much smaller than $2^{N}$. In practice, researchers often assume that only a small exposure set ${\mathcal{E}}$ is available, with the actual exposure received by unit $i$ under treatment assignment vector ${\boldsymbol{Z}}$ determined by an exposure mapping $m_{i}:{\mathcal{Z}}\rightarrow{\mathcal{E}}$. This leads to the assumption of correct exposure mapping as outlined below (Manski,, 2013; Aronow and Samii,, 2017; Hoshino and Yanagi,, 2023). Apparently, the potential outcomes of unit $i$ reduce to a much smaller set $\{Y_{i}(a)\}_{a\in{\mathcal{E}}}$ under such an assumption.

Different exposure mappings define different interference mechanisms. A general framework for establishing an interference mechanism is to introduce an interference network ${\mathcal{N}}$ over units of interest where units $i$ and $j$ are connected (denoted as $i\sim j$) if and only if interference may occur between them, and define

Under such an interference mechanism, we have three types of treatment exposures for each unit $i$, i.e., ${\mathcal{E}}=\{0,1,2\}$, where exposure $2$ implies that unit $i$ itself is assigned to treatment, exposure $1$ implies that units $i$ is assigned to control while one of its direct neighbors is under treatment, and exposure $0$ implies that unit $i$ and all its direct neighbors are all under control. Many popular interference mechanisms can be viewed as special cases of this framework. Below are three concrete examples.

In case of no interference, a classic problem in a randomized case-control study is to test the sharp null hypothesis of no treatment effect for all units. Following the notation in Example 3,

where $Y_{i}(1)$ and $Y_{i}(0)$ denote the potential outcomes of unit $i$ under treatment and control, respectively. Under $H_{0}^{\mathcal{F}}$, define $Y_{i}(1)=Y_{i}(0)=\theta_{i}$ for all $i\in{\mathbb{U}}$ and ${\boldsymbol{\theta}}=(\theta_{1},\cdots,\theta_{N})$. With the observed outcomes ${\boldsymbol{Y}}^{{\rm{obs}}}=\{Y_{i}^{{\rm{obs}}}\}_{i\in{\mathbb{U}}}$ where $Y_{i}^{{\rm{obs}}}=Z_{i}^{\rm{obs}}Y_{i}(1)+(1-Z_{i}^{\rm{obs}})Y_{i}(0)=\theta_{i}$, ${\boldsymbol{\theta}}$ is fully observed. Consider a finite-population framework where potential outcomes are fixed quantities, and the only source of randomness is the treatment assignment. The classic difference-in-means test statistic under a treatment assignment ${\boldsymbol{z}}$ is

Because ${\boldsymbol{\theta}}$ is a fixed vector after the randomization experiment is conducted, the null distribution of $T({\boldsymbol{Z}},{\boldsymbol{\theta}})$ is determined by the randomization distribution ${\boldsymbol{F}}_{\mathcal{Z}}$. Based on these facts, Fisher randomization test (FRT) can be conducted by calculating the p-value below:

where ${\boldsymbol{Z}}^{{\rm{obs}}}$ stands for the realized assignment, and ${\mathbb{P}}_{\boldsymbol{Z}}$ is the probability defined by ${\boldsymbol{Z}}\sim{\boldsymbol{F}}_{\mathcal{Z}}$. According to Fisher, (1935), FRT is finite-sample exact in the sense that under $H_{0}^{\mathcal{F}}$,

Hypothesis $H_{0}^{\mathcal{F}}$ is commonly referred to as the sharp null hypothesis, since the reference distribution of the test statistic $T({\boldsymbol{Z}},{\boldsymbol{\theta}})$ is fully determined under $H_{0}^{\mathcal{F}}$.

In presence of interference, however, the situation becomes more complicated. Given the treatment exposure mapping in Eq. (1), we aim to test the null hypothesis on the contrast between any two exposures $a$ and $b$ in ${\mathcal{E}}$, i.e.,

For example, we are often interested in testing the null hypothesis of no spillover effect, i.e., $H_{0}^{\{0,1\}}$. Under $H_{0}^{\{a,b\}}$, define $Y_{i}(a)=Y_{i}(b)=\theta_{i}$ for all $i\in{\mathbb{U}}$ and ${\boldsymbol{\theta}}=(\theta_{1},\cdots,\theta_{N})$, then the classic difference-in-means test statistic under treatment assignment ${\boldsymbol{z}}$ becomes

where ${\mathbb{U}}_{b}({\boldsymbol{z}})=\{i\in{\mathbb{U}}:m_{i}({\boldsymbol{z}})=b\}$, $N_{b}({\boldsymbol{z}})=|{\mathbb{U}}_{b}({\boldsymbol{z}})|$, ${\mathbb{U}}_{a}({\boldsymbol{z}})=\{i\in{\mathbb{U}}:m_{i}({\boldsymbol{z}})=a\}$, and $N_{a}({\boldsymbol{z}})=|{\mathbb{U}}_{a}({\boldsymbol{z}})|$. $T_{a,b}({\boldsymbol{z}},{\boldsymbol{\theta}})$ is well defined as long as $N_{a}({\boldsymbol{z}}),N_{b}({\boldsymbol{z}})>0$ for all ${\boldsymbol{z}}\in{\mathcal{Z}}$.

A critical challenge in testing $H_{0}^{\{a,b\}}$ in a randomized experiment with interference, however, lies in the fact that even under $H_{0}^{\{a,b\}}$, not all elements of ${\boldsymbol{\theta}}$ are observable after the experiment is conducted with the realized assignment ${\boldsymbol{Z}}^{{\rm{obs}}}$. To be concrete, define ${\mathbb{U}}^{{\rm{obs}}}_{a,b}={\mathbb{U}}_{a}({\boldsymbol{Z}}^{{\rm{obs}}})\cup{\mathbb{U}}_{b}({\boldsymbol{Z}}^{{\rm{obs}}})$, the partially observed version of ${\boldsymbol{\theta}}$ can be expressed as

where the question mark ‘$?$’ stands for the missing values. Define $N^{{\rm{obs}}}_{a,b}=|{\mathbb{U}}^{{\rm{obs}}}_{a,b}|$. Because only $N^{{\rm{obs}}}_{a,b}$ elements of ${\boldsymbol{\theta}}$ are observable, the null hypothesis $H_{0}^{\{a,b\}}$ is no longer sharp, resulting in an unknown reference distribution of the test statistic $T_{a,b}({\boldsymbol{Z}},{\boldsymbol{\theta}})$ under $H_{0}^{\{a,b\}}$.

To deal with this critical issue, Aronow, (2012) and Athey et al., (2018) proposed conditional randomization tests (CRTs) to substitute the classic FRT that is not applicable in this setting. The main idea is to avoid nuisance unknowns by focusing only on a subset of carefully selected units, referred to as “focal units”, and conducting randomized experiments on them with a randomization distribution restricted to a subset of treatment assignments, referred to as “focal assignments”. This approach ensures that the null hypothesis $H_{0}^{\{a,b\}}$ becomes sharp in the conditional inference. Because the reference distribution of the test statistic with focal units is precisely known conditional on the focal assignments and under the null hypothesis, an FRT-like testing procedure can be conducted. Basse et al., (2019) has shown that CRTs are exact both conditionally and unconditionally under mild conditions. Nonetheless, because CRTs employ only a restrictive subset of treatment assignments on a small subset of units of interest, they often suffer from severe power loss. To mitigate power loss, it is advisable to utilize as large a set of focal units associated with a set of focal assignments as possible, because a larger focal unit set often ensures a higher power, whereas a larger focal assignment set offers a higher resolution in obtaining the p-value. To achieve this, Puelz et al., (2022) proposed the concept of null exposure graph (NEG), which is a bipartite graph between units and assignments with an edge between unit $i\in{\mathbb{U}}$ and assignment ${\boldsymbol{z}}\in{\mathcal{Z}}$ if and only if unit $i$ is exposed to either $a$ or $b$ under ${\boldsymbol{z}}$. Encoding the impact of treatment assignments on units in terms of the observability of $\theta_{i}$’s under a given treatment assignment, NEG offers a way to better select focal units and focal assignments. For example, every biclique in NEG defines a set of appropriate focal units and focal assignments. Such an insight leads to a series of CRTs based on biclique decomposition of NEG, whose efficiency is further improved by Yanagi and Sei, (2022) via power analysis. More recently, Zhong, (2024) expanded the core concept of CRTs into a more general framework named partial null randomization tests (PNRTs), which reports an ensemble p-value by averaging multiple p-values obtained from a series of CRTs. Each CRT utilizes a small focal assignment set comprising only two focal assignments, yet incorporates an enlarged set of focal units including all focal units compatible with the specific pair of focal assignments. In PNRTs, the larger focal unit set of each individual CRT enhances testing power, while the ensemble p-value, obtained through averaging, addresses the issue of extremely low resolution in the individual p-values.

Despite these efforts, randomization tests based on the idea of conditional inference still suffer from diminishing power. For instance, consider a scenario where a high proportion of treatment $z_{i}=1$ arises from ethical concerns (e.g., the treatment significantly benefits participants). This typically leads to a high proportion of exposure $m_{i}({\boldsymbol{z}})=2$. If the goal is to test for no spillover effect $H_{0}^{{0,1}}$, the proportion of units with observable potential outcomes under this null hypothesis, $N_{0,1}^{\rm{obs}}/N$, tends to be small for many assignments ${\boldsymbol{z}}\in\mathcal{Z}$. This results in a sparse NEG, posing significant challenges in identifying sufficiently large bicliques within the NEG to conduct a randomization test with realistic power.

For a parametric statistical model $f_{\boldsymbol{\theta}}(x)$ with ${\boldsymbol{\theta}}\in\Theta$ as the vector of parameters, consider testing the following complex null hypothesis $H_{0}:{\boldsymbol{\theta}}\in\Theta_{0},$ where $\Theta_{0}$ is a non-empty subset of the parameter space $\Theta$. Given a collection of independent and identically distributed (i.i.d.) samples ${\boldsymbol{X}}^{{\rm{obs}}}=(X_{1},\cdots,X_{n})$ from model $f_{\boldsymbol{\theta}}(x)$, researchers often test $H_{0}$ based on a parameter-dependent test statistics $D({\boldsymbol{X}}^{{\rm{obs}}},{\boldsymbol{\theta}})$. In case that $\Theta_{0}$ contains only one element and thus ${\boldsymbol{\theta}}$ is fully specified under $H_{0}$, the p-value for testing $H_{0}$ can be naturally established as below:

where the probability ${\mathbb{P}}_{{\boldsymbol{X}}}$ is derived by random sample ${\boldsymbol{X}}\sim f_{\boldsymbol{\theta}}(x)$. In case that $\Theta_{0}$ contains more than one element and thus ${\boldsymbol{\theta}}$ is not fully specified under $H_{0}$, the most classic way to establish the p-value for $H_{0}$ is to find a proper estimate of ${\boldsymbol{\theta}}$ (e.g., the MLE) first, and then replace the unspecified ${\boldsymbol{\theta}}$ in (8) by its estimate $\hat{\boldsymbol{\theta}}$, resulting in the classic plug-in p-value (PIP) $p({\boldsymbol{X}}^{{\rm{obs}}};\hat{\boldsymbol{\theta}})$.

From the Bayesian perspective, $p({\boldsymbol{X}}^{{\rm{obs}}};\hat{\boldsymbol{\theta}})$ imputes the unspecified parameter ${\boldsymbol{\theta}}$ with a point estimation based on the observed data, which is obviously suboptimal. A natural alternative is to impute ${\boldsymbol{\theta}}$ based on its posterior distribution under $H_{0}$. Based on this insight, Meng, (1994) and Gelman et al., (1996) proposed to test $H_{0}$ by the posterior predictive p-value (PPP) defined as:

where $\pi({\boldsymbol{\theta}}\mid{\boldsymbol{X}}^{{\rm{obs}}},\Theta_{0})\propto f({\boldsymbol{X}}^{{\rm{obs}}}\mid{\boldsymbol{\theta}})\pi({\boldsymbol{\theta}}\mid\Theta_{0})$ is the posterior distribution of ${\boldsymbol{\theta}}$ given ${\boldsymbol{X}}^{obs}$ and a prior distribution $\pi({\boldsymbol{\theta}}\mid\Theta_{0})$ defined over $\Theta_{0}$. Notably, because typically $\pi({\boldsymbol{\theta}}\mid{\boldsymbol{X}}^{{\rm{obs}}},\Theta_{0})$ and $\hat{{\boldsymbol{\theta}}}$ both converge to the true parameter ${\boldsymbol{\theta}}$ as the sample size approaches infinity, the PPP ${p}_{B}({\boldsymbol{X}}^{{\rm{obs}}};\Theta_{0})$ and the PIP ${p}({\boldsymbol{X}}^{{\rm{obs}}};\hat{\theta})$ often shares the same asymptotical behavior. A key advantage of PPP over PIP is its ability to address uncertainty in specifying ${\boldsymbol{\theta}}$ via the posterior distribution $\pi({\boldsymbol{\theta}}\mid{\boldsymbol{X}}^{{\rm{obs}}},\Theta_{0})$, which makes the inference procedure more robust.

Many causal inference problems with complicated structures conform to the framework of PPP. For example, Rubin, (1998) utilized PPP to improve the power of FRT in randomized experiments with one-sided noncompliance. Under the null hypothesis of no treatment effect for compliers, the unknown compliance types of units serve as the nuisance parameter ${\boldsymbol{\theta}}$, whose value can be imputed based on its posterior predictive distribution under a parametric model with an appropriate prior distribution. Averaging the pristine p-values across the imputed datasets according to Eq. (9), PPP is obtained to summarize the observed evidence against the null hypothesis. Earlier efforts of utilizing multiple imputation based on posterior predictive distribution for Bayesian model checking were summarized in Guttman, (1967) and Rubin, (1981, 1984). In this paper, we aim to utilize the imputation idea to address the unobservable potential outcomes in randomization tests with interference.

Following the spirit of PPP, we propose the following imputation-based randomization test. First, we impute the missing elements of the partially observed ${\boldsymbol{\theta}}$ according to an imputation distribution $\pi$ to get a complete version of ${\boldsymbol{\theta}}$ referred to as ${\boldsymbol{\theta}}^{{\rm{imp}}}$ hereafter, and conduct the classic FRT based on the single imputation for the p-value below:

with the partially observed ${\boldsymbol{\theta}}$ in Eq. (4) substituted by the fully imputed ${\boldsymbol{\theta}}^{{\rm{imp}}}$. Next, we average $p_{{\rm{SI}}}$ according to the imputation distribution $\pi(\cdot\mid{\boldsymbol{\theta}}^{\rm{obs}})$, which is dependent on the realized assignment and observed outcomes, to get the following ensemble p-value:

and reject $H_{0}^{\{a,b\}}$ if $p_{{\rm{MI}}}$ is smaller than a pre-specified significance level $\alpha$. We term the above testing procedure as imputation-based randomization test (IRT). Considering that it’s difficult to calculate $p_{{\rm{MI}}}({\boldsymbol{Z}}^{{\rm{obs}}};{\boldsymbol{\theta}}^{\rm{obs}})$ exactly because the involved FRT p-values typically have no analytical form, we can approximate $p_{{\rm{MI}}}({\boldsymbol{Z}}^{{\rm{obs}}};{\boldsymbol{\theta}}^{\rm{obs}})$ via Monte-Carlo simulation, as suggested by Rubin, (1998), after conducting the randomized experiment according to the realized assignment ${\boldsymbol{Z}}^{{\rm{obs}}}$. The algorithm below summarizes the procedure of IRT. Define ${\boldsymbol{I}}(\cdot)$ as the indicator function.

Apparently, the imputation distribution $\pi$ plays an important role in the construction of IRT. Assuming that units in ${\mathbb{U}}$ are random samples from a super population ${\mathcal{P}}$ with $g(\cdot)$ as the marginal density of $\theta_{i}$ for $i\in{\mathbb{U}}$, the structure of the hypothesis testing problem defined in Eq. (5) naturally suggests the following strategy to impute ${\boldsymbol{\theta}}=(\theta_{1},\cdots,\theta_{N})$, given its observed version ${\boldsymbol{\theta}}^{\rm{obs}}$: we set $\theta_{i}^{{\rm{imp}}}=\theta_{i}^{{\rm{obs}}}$ for $i\in{\mathbb{U}}^{{\rm{obs}}}_{a,b}$, and we impute $\theta_{i}^{{\rm{imp}}}$ with a random sample from $g(\cdot)$ for $i\notin{\mathbb{U}}^{{\rm{obs}}}_{a,b}$, leading to the oracle imputation distribution

where $\delta(\cdot)$ is the Dirac delta function. We refer to the IRT based on $\pi_{o}({\boldsymbol{\theta}}^{{\rm{imp}}}\mid{\boldsymbol{\theta}}^{\rm{obs}})$ as IRT_o.

In practice, however, it’s often the case that distribution $g(\cdot)$ is unknown and needs to be estimated based on the observed subset in ${\boldsymbol{\theta}}^{{\rm{obs}}}$, namely $\tilde{{\boldsymbol{\theta}}}^{\rm{obs}}=\{\theta_{i}^{\rm{obs}}\}_{i\in{\mathbb{U}}^{\rm{obs}}_{a,b}}$ according to Eq. (7). Let $\hat{g}(\cdot\mid\tilde{\boldsymbol{\theta}}^{\rm{obs}})$ be the estimated marginal distribution of $\theta_{i}$, we get the following imputation distribution in action:

When distribution $g(\cdot)$ has a parametric form with ${\boldsymbol{\gamma}}$ as the model parameters, we can learn $g(\cdot)$ from $\tilde{\boldsymbol{\theta}}^{{\rm{obs}}}$ by inferring its parameter ${\boldsymbol{\gamma}}$ in either frequentist or Bayesian manner. We term this version of IRT as IRT_p.

If no appropriate parametric model is available, we can estimate $g(\cdot)$ via non-parametric approaches instead. For example, we can approximate the unknown $g(\cdot)$ by the empirical distribution of the observed potential outcomes $\tilde{\boldsymbol{\theta}}^{{\rm{obs}}}$, i.e.,

Alternatively, when $g(\cdot)$ is known to be a continuous distribution, we can rely on the smoothed version of the empirical distribution $\hat{g}$ instead as suggested by van der Vaart, (1994), i.e.,

where $K(\cdot)$ is a symmetric kernel function with bandwidth $h_{N}$ satisfying $\int K(u)du=1$, $\int u^{2}K(u)du<\infty$, and $K(u)=K(-u)$. We refer to the IRTs based on $\hat{g}_{e}$ and $\hat{g}_{s}$ as IRT_e and as IRT_s, respectively.

In this section, we present the theoretical results for IRT under the null hypothesis in Eq. (5), when the oracle imputation distribution is precisely known or needs to be estimated based on the observables. Some of the results are closely related to the results for PPP in Meng, (1994). A critical difference between IRT and PPP, however, is that the imputation distribution in IRT would not converge to a point mass, as in the classic PPP, with the increase of sample size. Such a fact leads to some unique issues in the theoretical analysis of IRT, making it non-trivial to establish the validity of IRT. First, we start with the simple case where IRT is conducted with the oracle imputation distribution, i.e., IRT_o.

Theorem 1 shows that IRT_o is more centered around $1/2$ than $U\sim{\rm{Unif}}(0,1)$ when the null hypothesis holds. The proof is detailed in the Supplementary Material. The theorem immediately leads to the following corollary based on Lemma 1 in Meng, (1994).

Corollary 1 ensures that IRT_o can control the type I error rate at most twice the significance level, i.e., $\alpha$, a result that is well-known for PPP. Although such a result cannot be further improved to more positive forms from the theoretical perspective, e.g., ${\mathbb{P}}\left(p_{{\rm{MI}}}({\boldsymbol{Z}}^{{\rm{obs}}};{\boldsymbol{\theta}}^{{\rm{obs}}})\leq\alpha\right)\leq\alpha$ for all $\alpha>0$, according to the counterexample in Dahl, (2006), our empirical studies in Sections 4 and 5 show that the type I error rate of IRT_o and its extensions are often well controlled at $\alpha$. Interestingly, such a phenomenon also occurs in other model-free, finite-sample exact inference methods, such as cross-conformal prediction, jackknife+, and CV+ (Vovk et al.,, 2018; Vovk and Wang,, 2020; Barber et al.,, 2021; Guan,, 2024). For instance, Barber et al., (2021) recommended the jackknife+ prediction interval for its typical empirical coverage of $1-\alpha$, although theoretically, it can only guarantee $1-2\alpha$ coverage due to certain pathological cases. These practices provide extra support for the application of IRT.

Next, we focus on the more challenging cases where IRT is conducted with the estimated imputation distribution as shown in Eq. (13). Define

as the ensemble p-value difference between the two versions of IRT with the estimated and oracle imputation distributions. Intuitively, if the distribution of $\Delta({\boldsymbol{Z}}^{\rm{obs}};{\boldsymbol{\theta}}^{\rm{obs}})$, which is determined by the joint distribution of $({\boldsymbol{Z}}^{\rm{obs}},{\boldsymbol{\theta}})$, concentrates at 0, the IRT under the estimated imputation distribution will closely resemble IRT_o with similar properties. The following lemma provides a positive answer to the above intuition, ensuring that if $\Delta({\boldsymbol{Z}}^{{\rm{obs}}};{\boldsymbol{\theta}}^{{\rm{obs}}})$ converges to 0 in probability, the asymptotic type I error rate of IRT with the corresponding estimated imputation distribution is also bounded by twice the significance level.

Denote the sampling distribution of $N_{a,b}^{\rm{obs}}$ derived by ${\boldsymbol{Z}}^{\rm{obs}}\sim{\boldsymbol{F}}_{\mathcal{Z}}$ as ${\boldsymbol{F}}_{N}$, and the corresponding probability mass function (p.m.f.) as

where ${\mathbb{P}}_{\boldsymbol{Z}}$ is defined over ${\boldsymbol{Z}}\sim{\boldsymbol{F}}_{\mathcal{Z}}$. The following theorem provides sufficient conditions under which $\Delta({\boldsymbol{Z}}^{{\rm{obs}}};{\boldsymbol{\theta}}^{{\rm{obs}}})=o_{p}(1)$, and thus Eq. (16) holds.

The above sufficient conditions imply that when the likelihood ratio between the estimated distribution $\hat{g}$ and the true distribution $g$ on the unobservables converges to 1 in $L^{1}$, IRT with the estimated imputation distribution will exhibit similar performance as IRT_o, whose behavior has been established previously. The corollary below shows that in a simple case where $g(\cdot)$ represents a normal distribution with known variance, the sufficient conditions can be theoretically validated under specific conditions.

For more complex cases, theoretical verification of the sufficient conditions in Theorem 2 becomes challenging, but can be approximated by numerical verification when $r_{N}^{\rm{mis}}$ converges to 0 at a sufficiently fast rate, e.g., $O(N^{-1/2})$. Further details are provided in the Supplementary Material.

In this section, we examine the proposed method under various interference mechanisms via simulation. First, we consider the clustered interference as defined in Example 1. For a group of $N=300$ units equally divided into $K\in\{30,75,150\}$ clusters, we are interested in testing the null hypothesis of no spillover effect $H_{0}^{\{0,1\}}$. Since the potential outcome $Y_{i}(2)$ is irrelevant to our interest, we only generate potential outcomes $Y_{i}(0)$ and $Y_{i}(1)$ by

where $\tau\in\{0,0.1,\cdots,1\}$. To identify the spillover effect under clustered interference, we consider the two-stage randomized experiment suggested by Hudgens and Halloran, (2008), where we randomly assign $\lfloor K/2\rfloor$ clusters to the treatment group and the remaining clusters to the control group in the first stage, and then randomly select one unit from each cluster in the treatment group to receive the treatment in the second stage with all the other units untreated. For each of the $3\times 11$ configurations of $(K,\tau)$, we generate 10 independent datasets according to the data simulation model in Eq. (20). For each simulated dataset with a specific set of potential outcomes, we conduct 1,000 parallel randomized experiments with the treatment assignment ${\boldsymbol{Z}}$ randomly specified based on the randomized design ${\boldsymbol{F}}_{\mathcal{Z}}$ according to the two-stage experiment.

Using the difference-in-means in Eq. (6) as the test statistic with the significance level to reject the null hypothesis set to $\alpha=0.05$, we compare the proposed IRT methods, including IRT_o, IRT_p, IRT_e and IRT_s, to a group of CRT methods with focal units and assignments selected by various algorithms (including CRT_o, CRT_b, CRT_db and CRT_pb) and two PNRT methods (namely PNRT_p and PNRT_m). Detailed configurations of the 10 involved competing methods can be found in Table 1. For each of these competing randomization tests, we measure its type I error rate via simulated datasets with $\tau=0$, and power via simulated datasets with $\tau\neq 0$.

We present the type I error rate across 10 datasets and the average power of the 10 competing methods in Fig. 1. For cases where $K=150$, however, the results of CRT_b and CRT_pb are absent because these two methods become infeasible in these cases. When $K=150$, the average missing rate of $\theta_{i}$ reaches up to 25%, leading to a null exposure graph with a density lower than 0.75. Consequently, it is impractical to find a biclique decomposition that meets the minimum size requirements established by Puelz et al., (2022) and Yanagi and Sei, (2022), which specify at least 20 focal assignments and 20 focal units, rendering CRT_b and CRT_pb infeasible. From the results, we make the following observations. First, while IRT and CRT methods successfully control the type I error rate around 0.05 across all 10 datasets, the two PNRT methods fail to do so and tend to yield type I error rates that are much smaller than 0.05 in all cases. Although our theory guarantees only a $2\alpha$ bound of the type I error rate for IRT methods asymptotically, these simulation results suggest that IRT can often successfully control the type I error rate below $\alpha$ in practice. Second, the IRT methods demonstrate significantly higher power than CRT and PNRT methods across all cases, especially for cases with a large number of clusters $K$. Although CRT_o and CRT_db also yield competitively high power in these cases, they are specifically designed for clustered interference only and can not be extended to other settings. Third, the IRT methods based on estimated distribution – namely IRT_e, IRT_p and IRT_s – exhibit performance similar to IRT_o, even though the missing rate of nuisance parameters is far from 0. These results confirm that the proposed IRTs are powerful tools for casual inference in presence of cluster interference.

Additional experiments where $Y_{i}(0)$ values are generated from heavier tail distributions, such as $\chi^{2}$ or $t$ distributions, but still modeled as a normal distribution, are reported in the Supplementary Materials, showing that IRT_p is robust to model mis-specification.

Next, we consider the spatial interference as defined in Example 2. For a group of $N=1,000$ units with spatial interference, the potential outcomes are simulated according to Eq. (20). We adopt the similar setting as in Yanagi and Sei, (2022) to specify the interference mechanism: each unit $i$ is associated with a 2-dimensional coordinate $\boldsymbol{x}_{i}\in{\mathbb{R}}^{2}$ according to

where $I_{2}$ denotes the $2\times 2$ identity matrix, and allow interference to occur between two unit $i$ and $j$ if and only if their Euclidean distance $d(i,j)=||\boldsymbol{x}_{i}-\boldsymbol{x}_{j}||_{2}\leq r$, with the distance of interference $r\in\{0.005,0.01,0.02\}$. For each configuration of $(r,p,\tau)$, we generate 10 independent datasets of sample size $N$ according to Eq. (20). For each simulated dataset, the treatment assignment ${\boldsymbol{Z}}=(Z_{1},\ldots,Z_{n})$ is sampled randomly for 1,000 times according to an independent Bernoulli trial where $Z_{i}\sim\text{Bernoulli}(p)$ for $i=1,\ldots,N$ with $p\in\{0.2,0.5,0.8\}$.

Using the difference-in-means in Eq. (6) as the test statistic, we aim to test the null hypothesis of no spillover effect, i.e., $H_{0}^{\{0,1\}}$, with a significance level of $\alpha=0.05$. The type I error rate and average power across 10 datasets for the methods described in Section 4.1 are measured for comparison. The results are presented in Fig. 2, where CRT_o and CRT_db are not included due to their inapplicability under spatial interference. In this figure, subfigures A, B, and C depict the type I error rate and power for $p=0.2$, $0.5$, and $0.8$, respectively. For $p=0.5$ and $0.8$, the results of CRT_b and CRT_pb are absent due to the high average missing rate of nuisance parameters, which can reach up to 50% and 80%.

These results show that the proposed IRTs outperform existing methods significantly under spatial interference as well as under clustered interference. First, the proposed IRT methods effectively control the type I error rate and exhibit significantly higher power compared to CRTs and PNRTs, especially in the case of $p=0.8$, where CRTs are not applicable and PNRTs have almost no power at all. Second, similar to the results under clustered interference, the other IRT methods exhibit similar performance to IRT_o even with a large missing rate of nuisance parameters. Additional results for $Y_{i}(0)$ generated from $\chi^{2}(4)$ and $t(4)$ are detailed in the Supplementary Material, which show that IRT methods remain robust across different distributions in the case of spatial interference as well.

Because all simulation results show that IRTs based on different estimated imputation distributions perform similarly to IRT_o, we recommend using IRT_e in practice for its ease of implementation and independence from additional assumptions for the marginal distribution of $\theta_{i}$.

We re-analyze the insurance purchase data collected from a two-round randomized experiment with multiple treatments conducted by Cai et al., (2015) to study the influence of social network among farmers on their decisions to purchase a new weather insurance product. In this experiment, $N=4,902$ households of rice farmers were randomly assigned to the first or the second rounds of insurance introduction sessions. The households assigned to both rounds were randomized to receive either a simple or an intensive introduction to the insurance product. The households assigned to the second round, however, were further divided into three random groups and informed additional information about the purchase decisions of the first-round in three different ways: no additional information at all (NoInfo), summary information in terms of the overall purchase rate (Overall), or detailed information on concrete purchase decisions of all households in the first round (Indiv). At the end, all involved households were randomized into $2+2\times 3=8$ treatment groups, as summarized in Fig 3. A stratified randomized experiment was conducted to randomly assign the 8 treatments to households.

Households were asked to make a purchase decision immediately after the insurance introduction sessions, making sure that there was no interference between households in the same round. However, there was a 3-day gap between the two rounds to allow households from the first round to personally discuss the insurance with their friends in the second round, implying a one-way interference from the first round to the second round. The social network of the 4,902 households was known in the experiment. Let ${\boldsymbol{W}}=(W_{ij})_{N\times N}$ be the adjacent matrix of the social network. We have $W_{ij}=1$ if two households $i$ and $j$ are friends, and $W_{ij}=0$ otherwise. The purchase decisions of all households, i.e., $Y_{i}$’s, were recorded as the output of the experiment, where $Y_{i}=1$ if household $i$ made purchase and $Y_{i}=0$ otherwise.

Containing rich ingredients, the insurance purchase data can be analyzed in different ways to answer different scientific questions. Here, we are interested in testing the spillover effect from the first-round households to the second-round households. To be concrete, let ${\mathbb{U}}=\{1,\dots,N\}$ be the $N=4,902$ households of interest, ${\mathcal{T}}=\{1,\ldots,8\}$ be the treatment set showed in Fig. 3, $Z^{\rm{obs}}_{i}\in{\mathcal{T}}$ be the observed treatment assignment of household $i\in{\mathbb{U}}$ and ${\boldsymbol{Z}}^{\rm{obs}}=(Z^{\rm{obs}}_{1},\ldots,Z^{\rm{obs}}_{N})$. Due to the possible interference between households through their social network, a household assigned to treatment $t\in{\mathcal{T}}$ may receive distinct exposures depending on how much prior information its friend households can bring to it before treatment: if some of its friend households get access to the intensive introduction of the insurance product a priori in the first stage, treatment $t$ would be conducted with strong prior information, leading to exposure $t^{(s)}$; if some of its friend households get access to the simple introduction a priori only, treatment $t$ would be conducted with weak prior information, leading to exposure $t^{(w)}$; if none of its friend households get access to the simple or intensive introduction a priori, treatment $t$ would be conducted with no prior information at all, leading to exposure $t^{(n)}$. Although we allow insurance content to spill over through social networks, for simplicity, we assume that households cannot obtain information about the first-round purchase decisions from social networks; they can only obtain this information from the experimenter as intended. Based on these insights, we work on the following exposure mapping from ${\mathcal{Z}}$ to exposure set ${\mathcal{E}}=\cup_{t\in{\mathcal{T}}}\{t^{(n)},t^{(w)},t^{(s)}\}$:

where indicator $a_{i}$ takes value in set $\{n,w,s\}$ highlighting the level of prior information household $i$ receives before treatment. Note that because $a_{i}$’s are fully determined by the treatment assignment ${\boldsymbol{z}}$ and the social network of the households, $m_{i}({\boldsymbol{z}})$ in Eq. (21) is a well defined exposure mapping.

For treatment $t\in{\mathcal{T}}$, we consider testing the following null hypothesis of no spillover effect, which is the contrast between $t^{(w)}$ and $t^{(s)}$ within treatment group $t$:

For this purpose, we use the difference-in-means between the $t^{(s)}$ exposure group and the $t^{(w)}$ exposure group as the estimator of the spillover effect and the test statistic for testing $H_{0}^{t}$ for treatment $t$, and establish two-sided p-value of the test via various randomization tests. For IRTs, we consider IRT_e and IRT_p with the following Beta-Binomial model for $\theta_{i}$’s: $\theta_{i}\sim\text{Binomial}(1,q)$, $q\sim\text{Beta}(1,1).$ For PNRTs, we try both PNRT_p and PNRT_m under the default setting. For CRTs, however, we find that all members in this method family are not applicable to this real data example: CRT_o and CRT_db cannot be applied to two-round experiments with multiple treatments and one-way interference, and CRT_b and CRT_pb are infeasible because the low density (e.g., 0.107 when $t=3$) of the null exposure graph makes the biclique decomposition infeasible. Therefore, we only compare the performance of IRT_e and IRT_p versus PNRT_p and PNRT_m here.

Table 2 summarizes the spillover effect for $t\in\{1,\cdots,8\}$ and the corresponding p-values obtained by different methods. From the table, we can see the following facts. First, IRTs, as well as PNRTs, suggest that spillover effects are not significant for treatment $1$ and $2$, which is consistent with the domain knowledge of no interference within the same round. Second, both IRT methods reveal significant spillover effects for treatment 3 and treatment 7, although only the spillover effect for treatment 3 is significant after considering the Bonferroni adjustment for multiple testing. Since households in the treatment 3 group received only a simple introduction with no additional information on the purchase decision in the first round, it’s reasonable that leaked information from the social network would significantly influence their purchase decision. However, both PNRT methods miss these signals, with all p-values from PNRT_m far from 0.05.

To further assess the performance of the proposed methods in practice, we train a predictive regression model for the potential outcomes based on the real data and generate simulated data from the fitted model for performance evaluation and comparison. To achieve this goal, we select 9 relevant covariates of the households according to Cai et al., (2015) as predictors to build a logistic regression model for the binary potential outcomes, resulting in a $4902\times 9$ covariate matrix ${\boldsymbol{X}}$. Because 6 of the 9 columns in ${\boldsymbol{X}}$ have missing values, we would like to impute these missing values for convenient model training and data simulation. Here, we follow the strategy suggested by Zhao and Ding, (2024) to impute all missing values in ${\boldsymbol{X}}$ with 0, and include the binary missing data indicators (1 for missing elements and 0 otherwise) of the 6 incomplete covariates as additional covariates, resulting in an extra $4902\times 6$ covariate matrix ${\boldsymbol{M}}$. Moreover, dummy variables recording actual exposures received by households in the experiment, i.e., ${\boldsymbol{E}}=(E_{i,e})_{i\in{\mathbb{U}},e\in{\mathcal{E}}}$, are also included in the model as covariates, where $E_{i,e}=1$ if household $i$ received exposure $e$ in the experiment. With these covariates, we fit the logistic model for the observed outcomes ${\boldsymbol{Y}}^{{\rm{obs}}}$ as our base model for data generation: $\text{logistic}({\boldsymbol{Y}}^{{\rm{obs}}}\sim 1+{\boldsymbol{E}}+{\boldsymbol{X}}+{\boldsymbol{M}}).$ Simulation datasets can be generated based on the above base model to evaluate the performance of various methods for testing spillover effects. For illustrative purpose, we focus on testing $H_{0}^{t}$ with $t=3$, which examines the spillover effect between exposures $3^{(s)}$ and $3^{(w)}$. To be specific, let coef$(s)$ and coef$(w)$ be the coefficients of dummy variables $E_{:,3^{(s)}}$ and $E_{:,3^{(w)}}$ obtained in the fitted base model, and $\Delta(s,w)$ be their difference. We modify the base model by resetting coef$(s)=$ coef$(w)+\Delta_{s,w}\cdot r$ for some $r\in(0,1)$, and generate potential outcomes ${\boldsymbol{Y}}(3^{(s)})$ and ${\boldsymbol{Y}}(3^{(w)})$ of the same sample size as the original experiment according to the modified model. If $r=0$, coef$(s)=$ coef$(w)$ in the modified model, and $H_{0}$ holds; if $r>0$, $H_{0}$ fails. Here, we adjust $r$ to align the differences in means of the generated ${\boldsymbol{Y}}(3^{(s)})$ and ${\boldsymbol{Y}}(3^{(w)})$ (denoted as $\tau$) to eight different values, ranging from 0 to 0.27, to simulate scenarios with varying signal strengths.

We apply PNRT_p, PNRT_m, IRT_e, and IRT_p to this simulated dataset, and visualize their type I error rate (for cases with $\tau=0$) and power (for cases with $\tau>0$) in Fig. 4. The results indicate that IRTs can effectively control the type I error rate and has high power to detect spillover effects, while PNRTs have almost no power at all.