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.
68,000 characters · 22 sections · 35 citation commands
Randomization Inference For the Always-Reporter Average Treatment Effect
Sample attrition, in which some units’ outcomes are unobserved after randomization, is common in field experiments and can introduce selection bias if it is systematically related to potential outcomes. A popular framework for addressing the attrition problem, since lee2009training, imposes a monotonicity assumption under which treatment affects reporting behavior in only one direction.\footnote{See also zhang2003estimation.} Under this assumption, the observed outcome distributions imply bounds on the average treatment effect of always-reporters (AR-ATE hereafter) —units who would report outcomes regardless of treatment assignment.\footnote{Other common approaches to addressing the sample attrition problem include worst-case bounds horowitz2000nonparametric and inverse-probability weighting under a conditional ignorability assumption robins1994estimation.}
Existing methods for conducting inference for AR-ATE under monotonicity rely on asymptotic approximations lee2009training, imbens2004confidence, stoye2009more. This paper proposes an alternative randomization-based testing procedure for inference on the AR-ATE. The primary advantage of our approach is that it delivers strong finite-sample guarantees under the sharp null hypothesis (e.g., no treatment effects for all always-reporters), while maintaining asymptotic validity when treatment effects may be heterogeneous. Specifically, our procedure is finite-sample valid for testing the sharp-null hypothesis and remains asymptotically valid under the weak null hypothesis that the average treatment effect for always-reporters is zero.
The principle that randomization inference based on properly chosen (e.g. studentizied) statistics delivers the dual guarantees as described above has long been known (see the related literature section below). This paper applies this principle to randomization inference in the presence of sample attrition. However, the structure of our problem differs substantially from previously studied settings. In particular, the target subpopulation (always-reporters) is unobservable, and the outcome distributions of always-reporters cannot be recovered exactly even under the sharp-null hypothesis.
To address this challenge, we adopt a worst-case p-value approach and consider a worst-case randomization test that maximizes the randomization p-value over all always-reporter configurations consistent with the observed data. Importantly, our model implies two balance conditions that can be tested under the sharp-null hypothesis: (i) balance of outcomes, as in standard randomization-inference settings, and (ii) balance in the number of always-reporters between the treatment and control groups. The second balance condition can be used to construct a two-step testing procedure or combined with the first condition to yield a chi-square–type test statistic. This stands in contrast to the existing literature, where the sharp null hypothesis typically implies only a balance condition on the outcomes. The second balance follows from the observation that, while always-reporters cannot be individually identified, random assignment implies that their expected counts are equal across the treatment and control groups.
The worst-case p-value approach delivers finite-sample valid inference under the sharp null hypothesis via a simple argument. Establishing asymptotic validity under mild conditions, however, requires more refined analysis than in the existing literature. In particular, we show that our procedure is uniformly asymptotically valid over a parameter region in which the differential attrition rate between treatment and control groups may be negligible. This empirically relevant region is generally ruled out by the analysis of existing procedures lee2009training.\footnote{A notable exception is semenova2025generalized, which proposes a pretesting-based approach.} Incorporating this region requires a refined conditional analysis, as it may induce nonstandard asymptotic distributions for the component of the statistic that tests balance in the number of always-reporters, depending on the particular sequence of parameters under consideration.
In addition to the statistical results, we provide computationally feasible algorithms for implementing the proposed statistical procedure. For outcomes with finite support (e.g., binary or categorical), we exploit the symmetry of randomization distribution to compute p-values via exhaustive enumeration. For continuous outcomes, we reformulate the p-value computation as a collection of smaller optimization problems, each solved using integer linear programming.
The remainder of the paper is organized as follows. Section (ref) discusses related literature, Section (ref) describes the setup and standard assumptions, Section (ref) introduces our statistical procedures, Section (ref) states additional assumptions and establishes the statistical guarantees, Section (ref) discusses computational implementation and Section(ref) presents simulation results.
Asymptotic inference methods for (trimming) bounds have been studied in imbens2004confidence,lee2009training,stoye2009more,samii2023generalizing,semenova2025generalized. Randomization inference and design-based inference with sample attrition or missing outcomes has been examined in lin2017placement,ivanova2022randomization,heussen2024randomization,heng2025design,li2025randomization. To our knowledge, none of the existing papers establish the dual guarantee in the monotone attrition framework as provided by our procedure. kline2025finite also consider randomization-based procedures for weak null hypotheses and their approaches requires a user-specified apriori bounds on the outcome support.
In general, heterogeneity-robust permutation/randomization inference is studied in neuhaus1993conditional, janssen1997studentized, janssen1999testing, janssen2003bootstrap, chung2013exact, canay2017randomization, wu2021randomization, zhao2021covariate, cohen2022gaussian, tuvaandorj2024robust and aronow2024randomization. The primary difference of our paper with respect to this literature has been discussed in the introduction.
We consider an experiment with $n$ units, where each unit is randomized to the treatment group or the control group. Let $y_i(d)$ denote the potential outcome of unit $i$ under treatment status $d\in \{0,1\}$. Let $r_i(d)$ denote the reporting status of the $i$th unit when under treatment status $d$. We write $r_i(1) = 1$ if the $i$th unit is assigned to treatment and is present in the following-up survey, and we write $r_i(1)=0$ if it is not present. The reporting status $r_i(0)$ is defined analogously.
As in lee2009training, we make the following assumption on the reporting status.
The monotonicity assumption rules out the existence of units who would report if untreated but not report if treated. We note that the direction of the monotonicity is not critical here. With slight modification, the same identification argument and statistical procedure work if one assumes that $r_i(1)\leq r_i(0)$ for all $i\in[n]$. The key requirement is that treatment affects reporting status in the same direction for all units.
Given Assumption (ref), an unit $i$ belongs to one of the three principal strata based on its potential reporting status. Specifically, it can be:
We refer to the sets of always-reporters, if-reporters, and never-reporters as the principal reporting strata frangakis2002principal. They are the groups defined with respect to the potential reporting status.
Let $\mathcal{A}$ denote the set of always-reporters, defined as:
The parameter of interest is the average treatment effect (ATE) of the always-reporters, defined as:
where we implicitly assume that $|\mathcal{A}|\geq 1$ so that the quantity is well-defined. We refer to this parameter as the always-reporter average treatment effect (AR-ATE).\footnote{For if-reporters, we never observe their control outcomes; for never-reporters, we observe neither their treated nor control outcomes. For the ATE parameters of these groups, it should be clear that we can do no better than worst-case bounds under a bounded outcome assumption. }
Under Assumption (ref) and using the argument in lee2009training, one can establish that AR-ATE is partially identified and a sharp bound on AR-ATE can be derived.\footnote{For example, see Section 7.4 of gerber2012field. } In particular, the logic of lee2009training implies
For statistical inference on the identified set, lee2009training uses asymptotic methods. In this paper, we retain the identification framework of lee2009training but replace asymptotics methods with a randomization-based inferential approach.
We assume that the experiments under consideration are completely randomized: among $n$ units, researchers assign $n_1$ units to treatment, chosen uniformly at random. For each unit $i$, we let $D_i$ be a random variable such that $D_i=1$ if unit $i$ is assigned to treatment, and $D_i=0$ otherwise. We note that $\textnormal{P}(D_i =1 ) = n_1/n$. We denote a random assignment vector following a complete randomization of $n_1$ out of $n$ units as $D=\left(D_i\right)_{i=1}^n\sim \textrm{CR}(n,n_1)$.
A sample of $n$ units is associated with an unknown potential-outcome-potential-reporting-status table $\left(\left(y_i(1),y_i(0),r_i(1),r_i(0)\right)\right)_{i=1}^n$. The units are randomly assigned to the treatment according to a realization of the assignment vector $D^{\textnormal{obs}}=(D^{\textnormal{obs}}_i)_{i=1}^n \sim \textrm{CR}(n,n_1)$. We observe reporting status for all units given by $R_i=r_i\left(D_i\right)$. For reported units with $R_i=1$, we observe their outcomes $Y_i=y_i\left(D_i\right)\in \mathbb{R}$. For units with $R_i=0$, we do not observe their outcomes and we denote their outcomes as $\mathbf{NA}$. The observed dataset hence consists of a set of outcome--assignment--reporting-status triples
We will also write $\mathcal{D}=\left(Y,D^\textnormal{obs},R\right)$ where $Y=(Y_i)_{i=1}^n$, $D^\textnormal{obs}=(D^\textnormal{obs}_i)_{i=1}^n$ and $R=(R_i)_{i=1}^n$ if needed.
We adopt the finite-population (design-based) framework, treating the potential outcomes and potential reporting status as fixed parameters imbens2015causal. The sole source of randomness in our model is the vector of random treatment assignments, and statistical uncertainties are evaluated exclusively with respect to this randomness.
With the setup above, we now state our problem. Given a sample of $n$ units and the set of always-reporters $\mathcal{A}$ as defined in ((ref)), we are interested in designing statistical inferential procedures that are valid as a test (but with different guarantees) for both the sharp-null hypothesis:
and the weak-null hypothesis:
We shall propose below procedures that are finite-sample valid for testing the sharp-null hypothesis and asymptotically valid for the weak-null hypothesis.
We first give a heuristic explanation for our inferential procedure. We denote $A_i=1$ if the unit $i$ is an always-reporter and $A_i=0$ if otherwise. Denote the binary vector of always-report indicators as $A=\left(A_i\right)_{i=1}^n\in\{0,1\}^n$. We note that $A$ is unknown to the researchers and is only partially revealed by the realized assignments (e.g. see discussion in Section (ref)).
If we know the set of always-reporters, randomization inference under the sharp-null hypothesis is straightforward. For example, we can use the (possibly studentized) absolute value of the difference-in-means statistic
as our test statistic, where $\tilde{D}=(\tilde{D}_i)_{i=1}^n \in \{0,1\}^n$ is an arbitrary treatment assignment, $Y=\left(Y_i\right)_{i=1}^n$ is the vector of (possibly missing) observed outcomes, and $A$ is the binary vector of always-report indicators defined above. We can obtain the randomization distribution of the statistics under the sharp-null hypothesis and reject the sharp-null hypothesis ((ref)) if the p-value $p(A)$ is less than a pre-specified level $\alpha$. A standard argument imbens2015causal shows that this test is a finite-sample valid level-$\alpha$ test for testing the sharp-null hypothesis. However, this procedure is infeasible because it relies on the knowledge of the set of always-reporters, which is unknown in practice.
To address this difficulty, we take a worst-case approach: we calculate worst-case p-value as the maximum across all randomization p-values based on sets of always-reporters consistent with the observed data. We shall explain below how to construct such sets of always-reporters in Section (ref). For now, let $\mathbb{A}(D^\textnormal{obs},R)$ be the set of sets of always-reporters that are consistent with the observed data. Formally, the worse-case p-value is defined as
For statistical testing, we will reject the sharp-null hypothesis if $ p^{\textrm{worst}}\leq \alpha$.
The statistical inferential algorithm is detailed in Algorithm (ref). Given a sample size $n$, we denote an arbitrary test statistic as:
which is a function mapping (possibly missing) observed outcomes, treatment assignments and always-reporters indicators to the extended real line. We discuss the choice of test statistics in Section (ref).
\begin{remark'} In practice, one does not need to iterate over all possible reporting tables. One can terminate the algorithm and output $0$ as soon as a single p-value is above $\alpha-\beta$. \end{remark'}
In what follows, Section (ref) describes the construction of all reporting tables compatible with the observed assignments and reporting statuses. Section (ref) introduces several test statistics.
To construct the worst-case p-value, we need to first construct the set of sets of always-reporters that are consistent with the observed data. Recall that $R_i$ is the reporting status of the $i$th unit and $D_i$ is its treatment assignment.
Conditioning on treatment assignments, it is possible to assign some subjects to principal reporting strata based on realized assignments and reporting statuses:
The toy example in Table (ref) illustrates the attribution procedure. Note that the always-reporters in the control group can be identified exactly. The ambiguity comes from the treated units with $D_i=1$ and $R_i=1$, which are a mixture of always-reporters and if-reporters. Hence the set of sets of always reporters that are consistent with the data can be described as:
We shall call each $A\in\mathbb{A}(D,R)$ a reporting table, which is a vector of always-reporter indicators.
Based on the observed missing pattern, we can pretest the number of always-reporters in the data. Because the treatment assignment is randomized independent of the reporting status, the fractions of always-reporters should on average be the same. This argument allows one to prune probabilistically improbable tables in $\mathbb{A}(D,R)$ that contain too many or too few always-reporters. This is an application of thr Bergers-Boo procedure in our setting berger1994p. For example, given a table $A\in\mathbb{A}\left(D,R\right)$, define the test-statistic:
and we reject a table if the p-value based on its randomization distribution is smaller than $\beta$. A detailed procedure is included in Algorithm (ref).\footnote{Note that the randomization distribution of $\widehat{\textrm{DIM}}(D,A)$ depends only on the number of always-reporters in $A$ because of the symmetry of complete randomization. Hence, all tables with the same number of always-reporters yield identical rejection decisions. This observation can substantially reduce the computational burden.} Other pre-testing procedures, such as the one based on tail bounds or test inversions rigdon2015randomization can also be employed.
Given a reporting table $A$ that is consistent with the data, define the Hajek estimator $\widehat{\tau}^{\textnormal{hj}}_n(Y,D,A)$:
and its variance estimator:
where,
We define the absolute value of the studentized-Hajek statistic as:
Apart from the balance of outcomes among always-reporters, we can include the balance of the number of always-reporters between the treated and control groups, leading to variations of Wald statistics. The variance of the statistic $ \widehat{\tau}^a_n$, $\operatorname{V}_n(A)$, is defined as
We now define two test statistics:
and
where we define $\lfloor x\rfloor_-=\max\{0,-x\}$ and $0/0=0$.
Let $\Pi_n$ be a distribution of random assignment variables $\{D_i\}_{i=1}^n$ implementing a completely randomized design. We first define a general parameter space encoding Assumption (ref).
Note that each element $\theta_n\in \Theta_n$ is a combination of potential outcomes, potential reporting status, and an experimental design. Each element $\theta_n$ completely determines distributions of our test statistics. For the asymptotic guarantee of testing the weak-null hypothesis, we need the following assumptions on our parameter space:
\begin{remark'} Condition (i) assumes that the always-reporters take an non-negligible share of the sample. Condition (ii) rules out extreme correlations between the potential treatment and control outcomes of the compliers, a common assumption in the finite-population setting. Condition (iii) is a standard technical condition that arises when applying the central limit theorem for triangular arrays, and may be thought of as excluding heavy-tailed or sparse data. \end{remark'}
We denote the parameter space of potential outcomes, potential reporting statuses, and experimental designs satisfying Assumptions 1, 2 and 3 with $n$ units and constants $\delta$, $r$, and $B$ as $\Theta_n\left(\delta,s,r,B\right)\subset \Theta_n $. For the weak-null hypothesis, we further define the parameter space:
For the sharp-null hypothesis, we define the parameter space:
Note that the parameter space for the sharp-null hypothesis needs not satisfy Assumption (ref) and Assumption (ref).
The following theorem states that using test statistics $\mathcal{T}_0$, $\mathcal{T}_1$ and $\mathcal{T}_2$ introduced in ((ref)), ((ref)), and ((ref)), the statistical procedure as in Algorithm (ref) provides finite-sample valid level-$\alpha$ test against the sharp-null hypothesis and asymptotically-valid level-$\alpha$ test against the weak-null hypothesis. \begin{theorem'} Consider the statistical procedure in Algorithm (ref) with test statistics $\mathcal{T}_n^0$, $\mathcal{T}_n^1$ and $\mathcal{T}_n^2$ defined in ((ref)), ((ref)), and ((ref)). Fix a significance level $\alpha\in(0,0.25]$ and a pre-testing level $\beta\in [0,\alpha)$. Let $\Theta^s_n$ be the parameter space of the sharp-null hypothesis defined in ((ref)). We have
Fix $\delta <1$, $s\in(0,1]$, $r \in (0,1/2]$, and $B > 0$. Let $\Theta^w_n(\delta,s,r,B)$ be the parameter space of the weak-null hypothesis defined in ((ref)). We have
\end{theorem'}
\begin{remark'} Our asymptotic uniform guarantee includes cases where the sample consists almost entirely of always-reporters and few (potentially zero) if-reporters. Other asymptotic methods generally rule out this portion of the parameter space lee2009training or rely on problem-specific tuning parameters SEMENOVA2025106055. Including this part of the parameter space necessitates a more involved mathematical analysis. \end{remark'}
We note that, as an intermediate step in the proof of Theorem (ref), we derived the asymptotic distributions (up to the slackness introduced by the variance bounds) of the statistics $\mathcal{T}_n^0$, $\mathcal{T}_n^1$, and $\mathcal{T}_n^2$ under the true but unknown always-reporter vector $A$. This result permits inference based on asymptotic critical values. Procedures based on critical values obtained from asymptotic approximations do not provide the finite-sample guarantee as in (ref); they retain only the asymptotic guarantee as in (ref). However, as we discuss in Section (ref), procedures based on asymptotic critical values are typically much easier to compute in practice.
The inferential algorithm based on asymptotic critical values is presented in Algorithm (ref). We first define three functions of \(n_A^1=\sum_{i=1}^n D_iA_i\), given \(n_1\), \(n_0\), and \(n_A\):
where $\textrm{V}_n(A)$ is defined in (ref) and we adopt the convention that $0/0=0$. These three functions correspond to the three statistics for balance in the number of always-reporters in test statistics (ref), (ref), and (ref).
\begin{theorem'} Consider the statistical procedure in Algorithm (ref). Fix a significance level $\alpha\in(0,0.25]$ and a pre-testing level $\beta\in [0,\alpha)$. Fix $\delta <1$, $s\in(0,1]$, $r \in (0,1/2]$, and $B > 0$. Let $\Theta^w_n(\delta,s,r,B)$ be the parameter space of the weak-null hypothesis defined in ((ref)). Define the event
Then we have,
\end{theorem'} \begin{remark'} We remark that the limiting distribution of $g_1(n_A^1)$ can depend on the sequence of always-reporter vectors. We give two simple examples here and other sequences are also possible and may yield different limits. If a sequence $\left(A^n\right)_{n=1}^{\infty}$ satisfies, for a constant $c<1$, \[ \frac{1}{n}\sum_{i=1}^n A_i^n \leq c, \quad \forall n, \] then $g_1(n_A^1)$ has a limiting chi-squared distribution with one degree of freedom. In contrast, if $\sum_{i=1}^n A_i^n=1$ for all $n$, then $g_1(n_A^1)$ has a point mass at $0$ as its limiting distribution. Both sequences are permitted by our parameter space. A similar comment applies to $g_2(n_A^1)$. Both the randomization-based procedure in Algorithm (ref) and the asymptotic procedure in Algorithm (ref) are agnostic to the particular sequence of always-reporter vectors, and both deliver uniformly asymptotically valid inference. \end{remark'}
In randomization inference, computing exact p-values by enumerating all possible treatment assignments is typically infeasible. Instead, p-values are usually approximated via Monte Carlo simulation, using a sufficiently large number of independent random assignments to approximate the randomization distribution.
Let $n_{mc}$ denote the number of Monte-Carlo simulations. For each simulation $s\in \{1,...,n_{mc}\}$, denote the $s$th simulated assignment variables as $D^s=\left(D^s_i\right)_{i=1}^n$. The computation of the worst-case p-value problem based on Monte-Carlo draws can be represented as:
with $\mathcal{T}_n^i, i\in\{0,1,2\}$ defined in (ref), (ref) and (ref) respectively, and $\mathbb{A}(D^\textnormal{obs},R)$ defined in (ref).
The optimization problem (ref) can be challenging to solve practically. If the number of reporting treated units is \(r_1=\sum_{i=1}^n D_i^{\textnormal{obs}}R_i\), then the set \(\mathbb{A}(D^\textnormal{obs},R)\) will generically contain \(2^{r_1}\) candidate always-reporter tables. Pre-testing the number of always-reporters in the treated group can shrink the search space, but typically not enough to make brute-force enumeration practical for moderately sized datasets---for example, even when \(r_1=50\).
We discuss two computation approaches in this section, designed for different outcome data variable types. Section (ref) considers the case where the observed outcome variables are discrete and have a small support. For example, this is the case when the outcome variables are binary, categorical or count-valued. In such a setting, exhaustive enumeration is feasible even with a large sample size, after exploiting the symmetry of complete randomization.
Section (ref) considers a second approach which can be used with continuous outcomes. It decomposes (ref) into smaller problems, each solvable using (linear) integer programming techniques. The resulting p-value, $p^{\textrm{worst,IP,mc}}$, has the guarantee that $p^{\textrm{worst,IP,mc}}\geq p^{\textrm{worst,mc}}$, hence enabling valid yet possibly conservative inference. In our simulations, we did not observe a major loss of power due to such relaxations, provided each subproblem is tight enough. See Lemma (ref).
We use the following notations through out the section. Given a dataset with a sample size $n$ and a vector of always-reporter indicators $A=\left(A_i\right)_{i=1}^n$, let $n_A=\sum_{i=1}^n A_i$ denote the number of always-reporters. Denote the random assignment variables $\left(D_i\right)_{i=1}^n\sim \textrm{CR}(n,n_1)$. We define the random variables $n_{A,1}=\sum_{i=1}^n D_iA_i$ and $n_{A,0}=\sum_{i=1}^n \left(1-D_i\right)A_i$.
Datasets with discrete outcomes (e.g., binary, categorical or count-valued) are common. With such dataset, exhaustive enumeration is feasible even with a large sample size, after exploiting the symmetry of complete randomization.
Suppose that the support of the observed outcomes, $\mathcal{S}$, has cardinality $K$. We enumerate the support as $\mathcal{S}=\{v_1,...,v_K\}$. For the randomization distribution under the sharp-null hypothesis, we define the number of always-reporters with outcome $v_k$ as $n_A^k=\sum_{i=1}^n A_i\mathds{1}\{Y_i=v_k\}$, $k\in [K]$, and the number of always-reporters assigned to treatment and control groups with outcome $k$ as $n_{A}^{1,k}=\sum_{i=1}^n D_iA_i\mathds{1}\{Y_i=v_k\}$ and $n_{A}^{0,k}=\sum_{i=1}^n \left(1-D_i\right)A_i\mathds{1}\{Y_i=v_k\}$, where $D=\left(D_i\right)_{i=1}^n$ is a generic assignment vector.
Given $n_A$ and $\{n_{A}^k\}_{k\in [K]}$, the distribution of the test statistics (ref), (ref) and (ref) are completely determined by the distribution of $\{n_{A}^{1,k}\}_{k\in [K]}$. For example, the squared studentized-Hajek statistics (ref) can be written as a function of $n_A$, $\{n_{A}^k\}_{k\in [K]}$, $n_{A,1}$ and $\{n_{A}^{1,k}\}_{k\in [K]}$:
where, for $a\in \{0,1\}$,
with $n_{A,1}=\sum_{k=1}^Kn_{A}^{1,k}$, $n_{A,0}=n_A-n_{A,1}$, and $n_A^{0,k}=n_A^k-n_A^{1,k}$. The probability mass function of $\{n_{A}^{1,k}\}_{k\in [K]}$ is
which depends on $\{n,n_1,\{n_{A}^k\}_{k=1}^K\}$. Hence, the randomization distributions of test statistics (ref), (ref) and (ref) under the sharp-null hypothesis are completely determined by the parameters $\{n,n_1,\{n_{A}^k\}_{k=1}^K\}$.
Given an always-reporters vector $A=\left(A_i\right)_{i=1}^n\in \mathbb{A}(D^\textnormal{obs},R)$, we denote the vector of outcome counts as $c\left(A\right)=\left(n_{A}^k\right)_{k=1}^K$. By the argument above, the randomization-based $p$-value associated with a given always-reporter table depends on $A$ only through the outcome counts vector $c(A)$, together with the sample size $n$ and the number of treated units $n_1$:
Let $\mathbb{C}(D^\textnormal{obs},R)$ denote all possible counts vector that is compatible with the data:
We have the reduction,
If the support of observed outcomes has cardinality $K$, the worst-case cardinality of $\mathbb{C}(D^\textnormal{obs},R)$ is upper-bounded by $n^K$. Provided that we consider a regime where $K$ does not increase with $n$, an exhaustive search algorithm has a time complexity that is polynomial in $n$. When $K$ is small, for example, $K=2$ for binary outcomes, a simple implementation would lead to a practical algorithm. A pseudo algorithm is included as Algorithm (ref) in Appendix Section (ref).
When the observed outcome data have large support, the reduction in Section (ref) does not yield practically efficient algorithms. To optimize over a large search space, we employ integer programming techniques.
Because $\mathcal{T}_n^0$, $\mathcal{T}_n^1$, and $\mathcal{T}_n^2$ are nonlinear functions of the always-reporter indicators $A$, a direct integer-programming formulation is challenging. In particular, after factorization, $\mathcal{T}_n^0$ and $\mathcal{T}_n^1$ can be expressed as sums of ratios of polynomials in $A$ of degree up to 10, and $\mathcal{T}_n^2$ is further complicated by the floor operator $\lfloor \cdot \rfloor_{-}$. Standard linearization techniques in IP are not easily applicable even with small instances.
To address the computational challenge, we decompose (ref) into a collection of smaller subproblems that are substantially easier to solve in our simulations. Moreover, the subproblems can be solved independently, making this step readily distributable across computing resources. The decomposition proceeds in two steps. First, we rewrite (ref) as a set of subproblems with a simpler algebraic structure (for example, involving lower-degree polynomials). Second, we further partition these subproblems by the number of always-reporters \(n_A\). The two-step reductions are useful because they decompose the original problem into integer programming subproblems with only quadratic (second-order) polynomials, eliminating non-polynomial features such as floor functions. For details, see the discussion after Theorem (ref).
Let $\mathcal{T}_n$ denote, generically, any one of the test statistics $\mathcal{T}_n^0$, $\mathcal{T}_n^1$, and $\mathcal{T}_n^2$. Let $L,U\in[0,\infty]$ be two (possibly infinity-valued) scalars such that
Relatively tight bounds $L$ and $U$ are easy to obtain, either analytically or via numerical procedures. See Algorithm (ref), Algorithm (ref) and Section (ref).
Consider a partition of the interval $[L,U]$ with increasingly ordered endpoints $\{t_i\}_{i=0}^I$, where $t_0=L$ and $t_I=U$. For each subinterval $[t_{i-1},t_i]$, consider the following value associated with each always-reporter vector $A$:
where $D\sim \textrm{CR}(n,n_1)$. Also consider the optimization problem:
Lemma (ref) suggests that the collection of optimal values of (ref) across different pairs of consecutive endpoints is informative about $p^{\textrm{worst}}$.
\begin{lemma'} Given endpoints $\{t_i\}_{i=0}^I$ with $t_0=L$ and $t_I=U$ satisfying (ref), let $A^{i*}$ denote an optimal always-reporter vector for problem (ref) on the interval $[t_{i-1},t_i]$. We have the following inequalities:
In particular, let $i^*$ be the optimal index of the problem $ \max_{i\in [I]}v_i$. We have
\end{lemma'}
From a practical perspective, this lemma has three implications. First, if $\max_{i\in[I]} v_i \le \alpha$, then it certifies that $p^{\mathrm{worst}} \le \alpha$. Second, provided that the randomization distributions induced by the optimizing always-reporter vectors are not overly concentrated on the interval $[t_{i^*-1},\, t_{i^*})$, the quantity $\max_{i\in[I]} v_i$ should be close to the true $p^{\mathrm{worst}}$.\footnote{Ideally, we would like to show that, as the maximum interval length tends to zero, the additional term on the right-hand side of (ref) also tends to zero. If the underlying distribution were continuous, this would follow from an application of the dominated convergence theorem. In our setting, however, the randomization distributions are discrete with point masses, so such arguments are not directly applicable. We are not aware of a simple refinement that yields a useful bound for this term.} Third, it transforms the task of comparing two ratios of polynomials into comparing a ratio of polynomials with a scalar. This reduction simplifies the computation problem. After Step 2 below, we can clear the denominators (variance estimators) in the test statistics associated with outcome balance by multiplying them on both sides, yielding inequalities in which both sides are polynomials of degree at most $2$.
We note that the lemma and the accompanying discussion remain valid, after replacing $p^{\mathrm{worst}}$ with $p^{\mathrm{worst,mc}}$, when the expectation (average) is taken over simulated assignments, as opposed to the exact expectation under the complete randomization distribution.
To make the discussion simpler, we introduce some additional notations. For a given always-reporter vector $A$, we index the always-reporters by \(a\in[n_A]\) via a bijection \(\pi_A:[n_A]\to [n]\). Without loss of generality, we label indices for both \(i\) and \(a\) so that the first $r_0=\sum_{i=1}^n (1-D_i^{\mathrm{obs}})R_i$ elements correspond to always-reporters assigned to the control group. We write \(\pi_A(a)=i\) when unit \(i\) is labeled as the \(a\)-th always-reporter. Our indexing convention implies $\pi_A(i)=i$ for $i\leq r_0$ and all $A\in\mathbb{A}(D^\textnormal{obs},R)$.
The original assignment variables \(D=(D_i)_{i\in[n]}\) and outcome variables \(Y=(Y_i)_{i\in[n]}\) indexed by $i$ then induce assignment variables and outcome variables indexed by \(a\), \[ \widetilde{D}^A_a := D_{\pi_A(a)},\quad \widetilde{Y}_a := Y_{\pi_A(a)},\quad a\in[n_A]. \] We also introduce matching variables $\{x^A_{ai}\}_{a\in [n_A],i\in[n]}$, where $x^A_{ai}=1$ if $\pi_A(i)=a$ and $x^A_{ai}=0$ otherwise. By our indexing convention, we have $x^A_{ai}=1$ if $i=a$ and $i\leq r_0$ for every always-reporter table $A\in \mathbb{A}(D^{\textnormal{obs}},R)$. With the matching variables, the outcome for the $a$th always-reporter can be expressed as $Y_a=\sum_{i=1}^n x^A_{ai}Y_i$.
Each always-reporter table $A$ induces the assignment variables $\tilde{D}^A=\left(\tilde{D}^A_a\right)_{a\in [n_A]}$ and matching variables $\{x^A_{ai}\}_{a\in [n_A],i\in [n]}$. It is important to note that due to the symmetry of complete randomization, $\tilde{D}^A$ has the same distribution for all always-reporter tables with the same number of always-reporters. Motivated by this fact, we write $\tilde{D}^A$ as $\tilde{D}^{n_A}$. We denote the distribution of $\tilde{D}^{n_A}$ as $\mathcal{L}(n,n_1,n_A)$.
We now note that \(\mathcal{T}_n^0\), \(\mathcal{T}_n^1\), and \(\mathcal{T}_n^2\) depend on \((Y,D,A)\) only through \(Y\), the induced assignment vector \(\tilde{D}^{A}\), and the matching variables \(x^{A}=\{x^{A}_{ai}\}_{a,i}\). Define \[ n_{A,1}=\sum_{a=1}^{n_A} D_a \qquad\text{and}\qquad n_{A,0}=\sum_{a=1}^{n_A} (1-D_a). \] For example, \(\mathcal{T}_n^0\) in (ref) can be written as
where \[ \widehat{\mu}_n(Y,D^A,x^A) =\widehat{\mu}^1_n(Y,D^A,x^A)-\widehat{\mu}^0_n(Y,D^A,x^A), \] with \[ \widehat{\mu}^1_n(Y,D^A,x^A)=\frac{1}{n_{A,1}}\sum_{a,i} D_a x_{ai} Y_i, \quad \widehat{\mu}^0_n(Y,D^A,x^A)=\frac{1}{n_{A,0}}\sum_{a,i} (1-D_a) x_{ai} Y_i, \] and \[ \widehat{\sigma}^{2,\textnormal{hj}}_n\!\left(Y,D^A,x^A\right) =\widehat{v}^1_n(Y,D^A,x^A)+\widehat{v}^0_n(Y,D^A,x^A), \] where \[ \widehat{v}^1_n(Y,D^A,x^A) =\frac{1}{n_{A,1}^2}\sum_{a,i} D_a x_{ai} Y_i^2 -\frac{1}{n_{A,1}}\Big(\widehat{\mu}^1_n(Y,D^A,x^A)\Big)^2, \] \[ \widehat{v}^0_n(Y,D^A,x^A) =\frac{1}{n_{A,0}^2}\sum_{a,i} (1-D_a) x_{ai} Y_i^2 -\frac{1}{n_{A,0}}\Big(\widehat{\mu}^0_n(Y,D^A,x^A)\Big)^2. \]
Define the space of matching variables associated with all always-reporter tables with $k$ always-reporters as
We note that the test statistics $\mathcal{T}_n^1$ and $\mathcal{T}_n^2$ can be re-expressed, in a similar fashion, as
The functions $g_1$ and $g_2$ are defined in (ref) and (ref) in Appendix (ref). We omit their explicit forms for brevity. The only point we use is that, for a fixed always-reporter size $n_A$, both $g_1$ and $g_2$ depend on the data only through $n_{A,1}$.
Lemma (ref) summarizes the discussion above. \begin{lemma'} Let $\widetilde{\mathcal{T}}_n$ denote one of the test statistics $\widetilde{\mathcal{T}}^0_n$, $\widetilde{\mathcal{T}}^1_n$, or $\widetilde{\mathcal{T}}^2_n$ defined in (ref) and (ref). Given an always-reporter table $A=(A_i)_{i=1}^n$, observed outcomes $Y$, and observed assignments $D^{\textnormal{obs}}$, recall the p-value $p(A)$ defined in (ref) and the quantity $v(A,t_{i-1},t_i)$ defined in (ref) with endpoints $t_{i-1}$ and $t_i$.
Define the p-value associated with the matching variables $x=\{x_{ai}\}_{a,i}$ by
and, with endpoints $t_{i-1}$ and $t_i$, define
where $\widetilde{D}^{n_A}\sim \mathcal{L}(n,n_1,n_A)$ and $\widetilde{D}^{n_A}_\textnormal{obs}$ is the assignment vector indexed by $a$ that is induced by the observed assignments.
Let $x^A=\{x^A_{ai}\}_{a,i}$ denote the matching variables induced by $A$. Then \[ \widetilde{p}\!\left(x^A\right)=p(A) \qquad\text{and}\qquad \widetilde{v}\!\left(x^A,t_{i-1},t_i\right)=v\!\left(A,t_{i-1},t_i\right). \] \end{lemma'} Lemma (ref) and Lemma (ref) imply the following theorem. \begin{theorem'} Let $L,U\in[0,\infty]$ be two scalars satisfying (ref), and consider a partition of the interval $[L,U]$ with increasingly ordered endpoints $\{t_i\}_{i=0}^I$, where $t_0=L$ and $t_I=U$. Recall the definition of $v_i$ from (ref). We have,
\end{theorem'} The theorem implies that the optimization problem on the left-hand side can be decomposed into smaller subproblems indexed by $k$, the number of always-reporters, and by $i$, the interval index.
In addition, it converts the optimization over always-reporter indicators \(A=\{A_i\}_{i=1}^n\) into an optimization over the matching variables \(\{x_{ai}\}_{a,i}\). Importantly, the component of the test statistics \(\mathcal{T}^1_n\) and \(\mathcal{T}^2_n\) that assesses balance in the always-reporter indicators does not depend on \(\{x_{ai}\}_{a,i}\). To illustrate this point, consider the following IP formulation of the subproblem \(\max_{x\in x(k,D^\textnormal{obs},R)} v(x,t_{i-1},t_i)\) for some $k$ and $i$. Let \(\{\tilde{D}^k_s\}_{s=1}^{n_{mc}}\) be \(n_{mc}\) Monte Carlo draws of the assignment vector (indexed by \(a\)) from \(\mathcal{L}(n,n_1,k)\). The subproblem can be written as
subject to
where (i) the functions \(g_i(\cdot)\), $i\in \{0,1,2\}$, correspond to different test statistics $\mathcal{T}_n^0, \mathcal{T}_n^1$ and $\mathcal{T}_n^2$, and are defined in (ref) and (ref) and; (ii) each \(L_s\) is any constant satisfying \[ L_s \ \leq\ \min_{x\in x(k,D^\textnormal{obs},R)} \Bigl\{ \widehat{\mu}^2_n\!\left(Y,D^A_s,x\right) +\bigl(g_i(n_{A,s}^1)-t_{i-1}\bigr)\widehat{\sigma}^{2,\textnormal{hj}}_n\!\left(Y,D^A_s,x\right) \Bigr\}, \] and the collection \(\{L_s\}_{s=1}^{n_{mc}}\) can be obtained either analytically or computationally; (iii) the indicator \(I_s\) must be set to zero whenever
and is otherwise unconstrained.
We note that \(g_i(n_{A,s}^1)\) is independent of the decision variables and can therefore be treated as a known scalar for each \(s\). Moreover, both \(\widehat{\mu}^2_n(Y,D^A_s,x)\) and \(\widehat{\sigma}^{2,\textnormal{hj}}_n(Y,D^A_s,x)\) are quadratic functions of the matching variables \(x\). By contrast, obtaining a formulation with comparable structure in the original decision-variable space \(A=(A_i)_{i=1}^n\) appears challenging.
We conclude by noting that the reduction in Section (ref), together with a technique analogous to that in Section (ref), can be combined with a bisection search to solve Algorithm (ref), which relies on critical values from the asymptotic distribution for inference. In practice, we find the resulting computations to be highly efficient. This approach is useful because it provides a fast heuristic that delivers a reasonable lower bound for the worst-case Monte Carlo $p$-value \(p^{\textrm{worst,mc}}\). The resulting solution can be used to warm-start the integer-programming method to accelerate computation, or to certify non-rejection at the chosen significance level.
Readers can find a complete pseudo algorithm incorporating the discussions above in Algorithm (ref).
This section evaluates the finite-sample performance and computational cost of the proposed worst-case randomization tests. We report (i) empirical rejection rates under the sharp-null scenario (test size), (ii) empirical rejection rates under an alternative with a positive AR-ATE (power), and (iii) runtime of the implementation.
Each simulated dataset contains $n=100$ units, with complete randomization assigning $n_1=50$ to treatment and $n_0=50$ to control. Units belong to one of three principal reporting strata under Assumption (ref). We set the population shares to \[ \pi_{\mathrm{AR}}=0.9,\qquad \pi_{\mathrm{IR}}=0.05,\qquad \pi_{\mathrm{NR}}=0.05. \]
We generate potential outcomes for all units according to the following model: \[ y_i(0) \sim \textrm{N}(0,1), \qquad y_i(1)=y_i(0)+\tau, \] where $\textrm{N}(0,1)$ denote a standard normal random variable and $\tau$ is the treatemtent effect. We consider two scenarios. $\tau=0$ and $\tau=1$.
We implement the worst-case randomization test in Algorithm (ref) using the following the one-sided chi-square statistics $\mathcal T^2_n$ in (ref). We choose the pruning step at $\beta=0.05$ and the nominal size at $\alpha=0.05$. The number of simulated randomization draws is 1000. All computations were performed on a Linux compute node equipped with an Intel Xeon Platinum 8268 CPU (2.90 GHz), providing 8 logical cores.
In the null scenario (ATE =0), the rejection rate is 0.002. In most simulations, the algorithm terminates at the heuristic stage once a reporting configuration yields a randomization p-value above the rejection threshold, resulting in a median runtime of 35.26 seconds and a 90th percentile runtime of 72.11 seconds. There are rare cases where the null is rejected and the algorithm must exhaustively verify all admissible configurations; in two such instances this required roughly 10 hours.
In the alternative scenario (ATE = 1), the rejection rate is 0.9218. This scenario sees a median runtime of 308 seconds and a 90th percentile runtime of 7,080 seconds. Relative to the null case, a larger share of simulations require additional verification beyond the heuristic stage, resulting in longer runtimes and a substantially higher maximum computation time of approximately 75,000 seconds.
In general, we find that when the lower bound of the test statistic produced by the heuristic stage lies in the range of approximately 6–9, the algorithm may require substantial computation time to verify the rejection decision. By contrast, when the null hypothesis is not rejected or when the hypothesis is rejected with a sufficiently large test statistic (e.g., greater than or equal to 9), the algorithm typically terminates quickly.