# Multiple conditional randomization tests††thanks: We thank Jos Twisk for
sharing two clinical trial datasets with us. We thank Stephen Bates, Jing Lei,
and Paul Rosenbaum for reading and commenting on a previous version of this
manuscript.

Yao Zhang  Department of Applied Mathematics and Theoretical Physics,
University of Cambridge  Qingyuan Zhao  Statistical Laboratory, Department of
Pure Mathematics and Mathematical Statistics, University of Cambridge

We establish a general sufficient condition on constructing multiple “nearly
independent” conditional randomization tests, in the sense that the joint
distribution of their p-values is almost uniform under the global null. This
property implies that the tests are jointly valid and can be combined using
standard methods. Our theory generalizes existing techniques in the literature
that use independent treatments, sequential treatments, or post-randomization,
to construct multiple randomization tests. In particular, it places no
condition on the experimental design, allowing for arbitrary treatment
variables, assignment mechanisms and unit interference. The flexibility of
this framework is illustrated through developing conditional randomization
tests for lagged treatment effects in stepped-wedge randomized controlled
trials. A weighted $Z$-score test is further proposed to maximize the power
when the tests are combined. We compare the efficiency and robustness of the
commonly used mixed-effect models and the proposed conditional randomization
tests using simulated experiments and real trial data.

##  1 Introduction

Randomization tests, as conceived by Fisher (1935, Section 21), have an
enduring impact on statistics and its applications. By basing the inference on
the physical act of randomization, randomization tests can control type I
errors in finite samples without any distributional or modeling assumptions.
Early literature on experimental design focused on exploring the asymptotic
equivalence (or non-equivalence) of randomization tests and tests based the
normal linear model (Pitman, 1937; Welch, 1937; Kempthorne, 1955). With the
advent of computers, randomization test became an attractive alternative to
model-based methods that are often sensitive to model misspecification. A
review of this literature can be found in Section 2 and an accompanying
article (Zhang and Zhao, 2022).

More recently, there is a surge of interest in developing conditional
randomization tests (CRT) to address challenges in analyzing modern, complex
experiments. Some notable advances include conditional randomization tests for
meta-analysis (Zheng and Zelen, 2008), covariate adjustment (Hennessy et al.,
2016), and unit interference (Athey et al., 2018; Basse et al., 2019; Puelz et
al., 2019). These recent developments have opened up many promising
applications of conditional randomization tests. In these problems, exploring
the role of conditioning is the key to testing some targeted causal hypotheses
that do not fully specify the potential outcomes.

Following this line of work, an important but underexplored problem is the
construction of multiple CRTs. Generally speaking, complex experiments offer
scattered evidence of causation. In model-based inference, this is often
tackled by specifying a shared parameter for the treatment effect. In
randomization inference, combining the scatter pieces of evidence is not
straightforward and requires a careful consideration of the experimental
design. A good example of this is the construction of evidence factors for
different statistical hypotheses derived from the same causal theory
(Rosenbaum, 2010, 2017; Karmakar and Small, 2020). There, independence is
crucial to interpret them as corroborating pieces of evidence. Another example
is the crossover trial design in which experimental units receive a sequence
of different treatments, so the evidence of causation spreads across time.

Our main theoretical contribution in this article is a general sufficient
condition in Section 3 under which the p-values of the CRTs are jointly valid
and “nearly independent”, meaning that their joint distribution is nearly
uniform under the global null. In such cases, the tests can be combined by
standard methods such as Fisher’s combined probability test and can be easily
extended to simultaneous testing using the closed testing procedure (Marcus et
al., 1976). Our theory places the condition only on the construction of
randomization tests but not on the underlying design. Thus, it is potentially
applicable to any design with arbitrary treatment variable, assignment
mechanism and unit interference.

The flexibility of our theoretical result is demonstrated by some
methodological development for testing lagged treatment effect in stepped-
wedge randomized controlled trials. In stepped-wedge trials, patients
crossover from control to treatment but the treatment initating times are
randomized along a time grid (Hussey and Hughes, 2007). This design offers a
good opportunity to test lagged effect: “Of particular interest … might be
whether the intervention has a delayed effect” (Hemming et al., 2018).
However, the most straightforward permutation tests that compare untreated
patients and treated patients with the specified time lag are not independent.
In Section 4, we develop an algorithm to construct independent CRTs for lagged
effect and propose an efficient way to combine the p-values. The robustness
and improved power of our method is validated by a set of numerical
experiments and real trial data in Section 5.

Finally, Section 6 concludes the article with some discussion on connections
with evidence factors and potential avenues for future research.

Notation. The set of natural integers is denoted by $\mathbb{N}$ and the set
of real numbers by $\mathbb{R}$. We use calligraphic letters for other sets,
boldface letters for vectors, upper-case letters for random quantities, and
lower-case letters for fixed quantities. We use $\bm{1}_{l}$ to denote a
length $l$ vector of ones and $\bm{0}_{l}$ a length $l$ vector of zeros. The
subscript $l$ is often omitted if the vector’s dimension is clear from the
context. We use a single integer in a pair of square brackets as a shorthand
notation for the indexing set from $1$: $[N]=\\{1,\dotsc,N\\}$. We use set-
valued subscript to denote a sub-vector; for example,
$\bm{Y}_{\\{1,3,5\\}}=(Y_{1},Y_{3},Y_{5})$.

##  2 Review: A single conditional randomization test

This section provides some background on conditional randomization tests; A
more thorough review can be found in an accompanying article (Zhang and Zhao,
2022). Consider an experiment on $N$ units. Let $\bm{Z}\in\mathcal{Z}$ denote
the realized treatment assignment of all the units. A common setup for causal
inference is to assume that $\bm{Z}=(Z_{1},\dotsc,Z_{N})$ is an $N$-vector and
there is _no interference_ among the units, but this is not essential for the
theory below. Each unit $i$ has a vector of real-valued _potential outcomes_
$(Y_{i}(\bm{z})\mid\bm{z}\in\mathcal{Z})$ and an _observed outcome_
$Y_{i}=Y_{i}(\bm{Z})$. This is often referred to as the _consistency
assumption_ in the causal inference literature. The well-known _stable unit
treatment value assumption_ (SUTVA) assumes both no interference and
consistency (Rubin, 1980). Denote
$\bm{Y}(\bm{z})=(Y_{1}(\bm{z}),\dotsc,Y_{N}(\bm{z}))\in\mathcal{Y}\subseteq\mathbb{R}^{N}$
and $\bm{Y}=(Y_{1},\dotsc,Y_{N})\in\mathcal{Y}$. The _potential outcomes
schedule_ for all units and treatment assignments is denoted by
$\bm{W}=(\bm{Y}(\bm{z}):\bm{z}\in\mathcal{Z})\in\mathcal{W}$.

We make the following minimal assumption the _treatment assignment mechanism_.

######  Assumption 1 (Randomized experiment).

$\bm{Z}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{W}$ and the density function $\pi(\cdot)$
of $\bm{Z}$ (with respect to some reference measure on $\mathcal{Z}$) is known
and positive everywhere.

Independence between $\bm{Z}$ and $\bm{W}$ is justified by the randomization.
We may view $\bm{Z}$ and $\bm{W}$ as defined on two separate probability
spaces as $\bm{W}$ represents randomness in nature and consider the product
probability space. If $\bm{Z}$ is randomized according to some $\bm{X}$
(possibly some observed attributes of the units), we may assume that
$\bm{Z}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{W}\mid\bm{X}$. To simplify the exposition,
we will simply treat $\bm{X}$ as fixed, so $\bm{Z}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{W}$ is true in the conditional probability
space given $\bm{X}$.

###  2.1 Partially sharp null hypotheseses

Treatment effects are defined as differences between potential outcomes, but
only one potential outcome $Y_{i}=Y_{i}(\bm{Z})$ is observed for every unit
$i$ under the consistency assumption. A causal hypothesis defines a number of
relationships between the potential outcomes. Each relationship allows us to
impute some of the potential outcomes if another potential outcome is
observed. We can summarize these relationships using a set-valued mapping.

######  Definition 1.

A (partially) sharp null hypothesis $H$ defines an _imputability mapping_

$\displaystyle\mathcal{H}:\leavevmode\nobreak\ $
$\displaystyle\mathcal{Z}\times\mathcal{Z}\to 2^{[N]},$
$\displaystyle(\bm{z},\bm{z}^{*})\mapsto\mathcal{H}(\bm{z},\bm{z}^{*}),$

where $\mathcal{H}(\bm{z},\bm{z}^{*})$ is the largest subset of $[N]$ such
that $\bm{Y}_{\mathcal{H}(\bm{z},\bm{z}^{*})}(\bm{z}^{*})$ is imputable from
$\bm{Y}(\bm{z})$ under $H$.

Given the observed assignment $\bm{Z}$ and another assignment
$\bm{z}^{*}\in\mathcal{Z}$, $\mathcal{H}(\bm{Z},\bm{z}^{*})$ informs us the
subset of units whose outcome
$\bm{Y}_{\mathcal{H}(\bm{Z},\bm{z}^{*})}(\bm{z}^{*})$ is imputable from the
observed outcome $\bm{Y}=\bm{Y}(\bm{Z})$. A randomization test for $H$ is
implemented by comparing the test statistics under the realized $\bm{Z}$ and
all other $\bm{z}^{*}\in\mathcal{Z}$. Thus, the mapping
$\mathcal{H}(\bm{Z},\bm{z}^{*})$ essentially tells us the largest subset of
units we can use in a randomization test that considers comparing $\bm{Z}$
with $\bm{z}^{*}$. Following Definition 1, we call a hypothesis $H$ _sharp_ if
$\mathcal{H}(\bm{z},\bm{z}^{*})=[N]$ for any
$\bm{z},\bm{z}^{*}\in\mathcal{Z}$. Otherwise, we call $H$ _partially sharp_.
Some examples are provided to illustrate the definitions.

######  Example 1 (Simple experiment with two treatment arms).

Suppose the treatment for each unit is binary, i.e., $Z_{i}=0$ (control) or
$1$ (treated). Consider the constant treatment effect hypothesis
$H:Y_{i}(1)=Y_{i}(0)+\tau,\forall i\in[N]$. With no interference and
consistency, we can impute any potential outcome by
$Y_{i}(z_{i}^{*})=Y_{i}+(z_{i}^{*}-Z_{i})\tau$. Thus,
$\mathcal{H}(\bm{z},\bm{z}^{*})=[N]$ for all
$\bm{z},\bm{z}^{*}\in\mathcal{Z}=\\{0,1\\}^{N}$ and the hypothesis $H$ is
sharp.

######  Example 2 (Interference).

Consider a randomized experiment in which units interfere with others in the
same cluster (Hudgens and Halloran, 2008). For every unit $i$, $Z_{i}$ can
take four different values $0$, $1$, $2$ and $3$, denoting control, neighbour-
treated, itself-treated, neighbour- and itself-treated respectively. Consider
a hypothesis $H:Y_{i}(1)=Y_{i}(0)+\tau,\forall i\in[N].$ If $Z_{i}=2$ or $3$,
the potential outcomes $Y_{i}(0)$ and $Y_{i}(1)$ are not imputable from
$Y_{i}$ under $H$. For any $i$ such that $z_{i}\not\in\\{0,1\\}$ and
$z_{i}^{*}\in\\{0,1\\}$, we have $i\not\in\mathcal{H}(\bm{z},\bm{z}^{*})$.
Thus, $\mathcal{H}(\bm{z},\bm{z}^{*})$ is a strict subset of $[N]$ in general,
implying that $H$ is partially sharp. As a comparison, a fully sharp
hypothesis is $H:Y_{i}(3)=Y_{i}(2)=Y_{i}(1)=Y_{i}(0)+\tau,\forall i\in[N]$,
which assumes that the treatment effect is the same for any treatment statuses
of a unit and its neighbour.

######  Example 3 (Stepped-wedge trial (Brown and Lilford, 2006) or staggered
adoption (Abraham and Sun, 2021; Athey and Imbens, 2022)).

Consider a sequential experiment that randomizes each unit’s treatment
starting time $Z_{i}\in[T]\cup\\{\infty\\}$; $Z_{i}=\infty$ denotes never-
treated. Suppose we only observe an outcome at the end of the experiment,
$Y_{iT},\forall i\in[N]$. Consider the hypothesis
$H:Y_{iT}(1)=Y_{iT}(\infty)+\tau,\forall i\in[N]$, i.e., treating at time $1$
has a constant effect $\tau$. This hypothesis is partially sharp because
$Y_{i}(z_{i}^{*})$ is not imputable from $Y(z_{i})$ if
$z_{i}\in\\{2,3,\dotsc,T\\}$ and $z_{i}^{*}\in\\{1,\infty\\}$. As a
comparison, a sharp hypothesis is given by
$H:Y_{iT}(1)=\dotsc=Y_{iT}(T)=Y_{iT}(\infty)+\tau,\forall i\in[N]$, which
further assumes that the treatment starting time does not alter the treatment
effect.

###  2.2 Conditional randomization tests

A common feature of the complex designs above in Examples 2 and 3 is that the
units have more than two potential outcomes, but the causal hypotheses do not
involve all of them. In this case, $\mathcal{H}(\bm{z},\bm{z}^{*})$ depends on
$\bm{z}$ and $\bm{z}^{*}$ in a non-trivial way. Conditioning, as described
below, offers a principled way to address nonimputable potential outcomes
arising from testing partially sharp hypotheses.

A conditional randomization test involves conditioning on a subset of
treatment assignments in a countable partition
$\mathcal{R}=\\{\mathcal{S}_{m}\\}_{m=1}^{\infty}$ of $\mathcal{Z}$ where
$\mathcal{S}_{1},\mathcal{S}_{2},\dotsc$ are (measurable) subsets. We require
that $\int_{\mathcal{S}_{m}}\pi(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}>0$ for all
$m$ to avoid conditioning on zero probability events, although this can be
relaxed if suitable conditional density functions can be defined. For each
$m$, a conditional randomization distribution can be obtained by evaluating a
test statistic $T_{m}$ under all the assignments
$\bm{z}^{*}\in\mathcal{S}_{m}$. If the observed assignment
$\bm{z}\in\mathcal{S}_{m}$, the conditional randomization test will compare
observed test statistic with the corresponding conditional randomization
distribution. For any $\bm{z}\in\mathcal{S}_{m}$, we write
$\mathcal{S}_{\bm{z}}=\mathcal{S}_{m}$ and
$T_{\bm{z}}(\cdot,\cdot)=T_{m}(\cdot,\cdot)$. This allows us to define the
p-value below by conditioning on $\bm{Z}^{*}\in\mathcal{S}_{\bm{Z}}$.

######  Definition 2.

A _conditional randomization test_ (CRT) for a treatment $\bm{Z}$ is defined
by a countable _partition_ $\mathcal{R}=\\{\mathcal{S}_{m}\\}_{m=1}^{\infty}$
of $\mathcal{Z}$ such that $\mathcal{S}_{m},m=1,2,\dotsc,$ of $\mathcal{Z}$
(disjoint subsets such that
$\bigcup_{m=1}^{\infty}\mathcal{S}_{m}=\mathcal{Z}$) and a collection of _test
statistics_ $(T_{m}(\cdot,\cdot))_{m=1}^{\infty}$, where
$T_{m}:\mathcal{Z}\times\mathcal{W}\to\mathbb{R}$ is a (measurable)
function.The _p-value_ of the CRT is given by

$P(\bm{Z},\bm{W})=\mathbb{P}^{*}\\{T_{\bm{Z}}(\bm{Z}^{*},\bm{W})\leq
T_{\bm{Z}}(\bm{Z},\bm{W})\mid\bm{Z}^{*}\in\mathcal{S}_{\bm{Z}},\bm{W}\\},$ (1)

where $\bm{Z}^{*}$ is an independent copy of $\bm{Z}$ conditional on $\bm{W}$.

In this definition, we use the notation $\mathbb{P}^{*}$ to emphasize that the
probability is taken over the randomness of $\bm{Z}^{*}$. When
$\mathcal{S}_{\bm{z}}=\mathcal{Z}$ for all $\bm{z}\in\mathcal{Z}$ so the
partition is trivial, the CRT reduces to an unconditional randomization test.

Denote the $\sigma$-algebra generated by the conditioning events in the
p-value (1) by

$\mathcal{G}=\sigma\left(\\{\bm{Z}\in\mathcal{S}_{m}\\}_{m=1}^{\infty}\right)$

Because $\\{\mathcal{S}_{m}\\}_{m=1}^{\infty}$ is a partition of
$\mathcal{Z}$, $\mathcal{G}$ consists of all countable unions of
$\\{\mathcal{S}_{m}\\}_{m=1}^{\infty}$. This allows us the formalize the
validity of a CRT as follows (a proof can be found in Zhang and Zhao (2022)).

######  Theorem 1.

Consider a CRT defined as above. The p-value $P(\bm{Z},\bm{W})$ stochastically
dominates the uniform distribution on $[0,1]$ in the following sense:

$\operatorname{\mathbb{P}}\left(P(\bm{Z},\bm{W})\leq\alpha\mid\mathcal{G},\bm{W}\right)\leq\alpha,\quad\forall\alpha\in[0,1].$
(2)

In consequence, given Assumption 1 and a partially sharp null hypothesis $H$,
if $P(\bm{Z},\bm{W})$ is computable, then

$\mathbb{P}\left\\{P(\bm{Z},\bm{Y})\leq\alpha\right\\}\leq\alpha,\quad\forall\alpha\in[0,1].$
(3)

The conditional probability on the left-hand side of (2) is a well defined
random variable (by the Radon-Nikodym theorem as $P(\bm{Z},\bm{W})$ is
measurable). Because $\mathcal{G}$ is generated by a countable partition, the
left-hand side can be written as

$\operatorname{\mathbb{P}}\left(P(\bm{Z},\bm{W})\leq\alpha\mid\mathcal{G},\bm{W}\right)=\sum_{m=1}^{\infty}1_{\\{\bm{Z}\in\mathcal{S}_{m}\\}}\operatorname{\mathbb{P}}\left(P(\bm{Z},\bm{W})\leq\alpha\mid\bm{Z}\in\mathcal{S}_{m},\bm{W}\right).$

Then (2) is equivalent to

$\mathbb{P}\left\\{P(\bm{Z},\bm{W})\leq\alpha\mid\bm{Z}\in\mathcal{S}_{\bm{z}},\bm{W}\right\\}\leq\alpha,\quad\forall\alpha\in[0,1],\bm{z}\in\mathcal{Z}.$
(4)

In (3), we denote the p-value $P(\bm{Z},\bm{W})$ by $P(\bm{Z},\bm{Y})$ with an
abuse of notation since it can be computed from the observed outcomes
$\bm{Y}$.

A level-$\alpha$ CRT rejects the null hypothesis if $P(\bm{Z},Y)\leq\alpha$.
By further randomizing the rejection when $P(\bm{Z},\bm{Y})$ is close to
$\alpha$, (3) holds with equality for any $\alpha\in[0,1]$; see Lehmann and
Romano (2006, Equation 5.51). To construct a computable and valid CRT, we
generally need to restrict the dependence of $T_{\bm{Z}}(\bm{Z}^{*},\bm{W})$
on $\bm{W}$. For instance, to test the partially sharp null hypothesis in
Example 3, we can compute the test statistics only with the units that are
treated at time $1$ or never treated. To improve the power of a CRT, the
partition $\mathcal{R}$ need to be carefully designed. For more discussion and
examples, see Zhang and Zhao (2022, Sections 3 and 4) and references therein.

##  3 Multiple conditional randomization tests

We next establish a sufficient condition under which multiple CRTs are jointly
valid and nearly independent. Consider $K$ conditional randomization tests for
$K$ possibly different and related hypotheses $H^{(k)}$ that are defined by
partitions
$\mathcal{R}^{(k)}=\left\\{\mathcal{S}^{(k)}_{m}\right\\}_{m=1}^{\infty}$ and
test statistics $(T^{(k)}_{m}(\cdot,\cdot))_{m=1}^{\infty}$, $k=1,\dotsc,K$.
We denote the corresponding p-values in (1) as
$P^{(1)}(\bm{Z},\bm{W}),\dotsc,P^{(K)}(\bm{Z},\bm{W})$.

###  3.1 Main theorem

Some additional notations are needed before stating the main theorem. For any
subset of tests $\mathcal{J}\subseteq[K]$, we define the _union_ ,
_refinement_ and _coarsening_ of the conditioning sets as

$\mathcal{R}^{\mathcal{J}}=\bigcup_{k\in\mathcal{J}}\mathcal{R}^{(k)},\quad\underline{\mathcal{R}}^{\mathcal{J}}=\left\\{\bigcap_{j\in\mathcal{J}}\mathcal{S}_{\bm{z}}^{(j)}:\bm{z}\in\mathcal{Z}\right\\},\quad\text{and}\quad\overline{\mathcal{R}}^{\mathcal{J}}=\left\\{\bigcup_{j\in\mathcal{J}}\mathcal{S}^{(j)}_{\bm{z}}:\bm{z}\in\mathcal{Z}\right\\}.$

As before, let
$\mathcal{G}^{(k)}=\sigma(\\{\bm{Z}\in\mathcal{S}\\}:\mathcal{S}\in\mathcal{R}^{(k)})$
denote the $\sigma$-algebra generated by the conditioning events for test $k$.
Let $\mathcal{G}^{\mathcal{J}}$ be the $\sigma$-algebra generated by the sets
in $\mathcal{R}^{\mathcal{J}}$, so
$\mathcal{G}^{\mathcal{J}}=\sigma\big{(}\\{\bm{Z}\in\mathcal{S}\\}:\mathcal{S}\in\mathcal{R}^{\mathcal{J}}\big{)}$.
Similarly, the refinement and coarsening $\sigma$-algebras are defined as
$\underline{\mathcal{G}}^{\mathcal{J}}=\sigma\big{(}\\{\bm{Z}\in\mathcal{S}\\}:\mathcal{S}\in\underline{\mathcal{R}}^{\mathcal{J}}\big{)}$
and
$\overline{\mathcal{G}}^{\mathcal{J}}=\sigma\big{(}\\{\bm{Z}\in\mathcal{S}\\}:\mathcal{S}\in\overline{\mathcal{R}}^{\mathcal{J}}\big{)}$.

######  Theorem 2.

Suppose the following two conditions

$\underline{\mathcal{R}}^{\\{j,k\\}}\subseteq\mathcal{R}^{\\{j,k\\}},$ (5)
$T^{(j)}_{\bm{Z}}(\bm{Z},\bm{W})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}T_{\bm{Z}}^{(k)}(\bm{Z},\bm{W})\mid\underline{\mathcal{G}}^{\\{j,k\\}},\bm{W},$
(6)

are satisfied for all $j,k\in[K]$, $j\neq k$. Then we have

$\begin{split}\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{W})\leq\alpha^{(1)},\dotsc,P^{(K)}(\bm{Z},\bm{W})\leq\alpha^{(K)}\mid\overline{\mathcal{G}}^{[K]},\bm{W}\right\\}\leq\prod_{k=1}^{K}\alpha^{(k)},\\\
\forall\alpha^{(1)},\dotsc,\alpha^{(K)}\in[0,1].\end{split}$ (7)

In consequence, given Assumption 1 and that the null hypotheses
$H^{(1)},\dotsc,H^{(K)}$ are true, if the CRTs are computable, then

$\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{Y})\leq\alpha^{(1)},\dotsc,P^{(K)}(\bm{Z},\bm{Y})\leq\alpha^{(K)}\right\\}\leq\prod_{k=1}^{K}\alpha^{(k)},\quad\forall\alpha^{(1)},\dotsc,\alpha^{(K)}\in[0,1].$
(8)

Theorem 2 provides a sufficient condition for jointly valid CRTs: the
partitions $\mathcal{R}^{(k)},k\in[K]$ are nested as in (5), and the test
statistics are conditional independent as in (6). We illustrate these
conditions by considering $K=2$. In this case, condition (5) amounts to
assuming a nested structure between the two partitions:
$\mathcal{S}^{(1)}_{\bm{z}}\subseteq\mathcal{S}^{(2)}_{\bm{z}}$ or
$\mathcal{S}^{(1)}_{\bm{z}}\supseteq\mathcal{S}^{(2)}_{\bm{z}}$ for all
$\bm{z}\in\mathcal{Z}$. In other words,
$\mathcal{S}^{(1)}_{\bm{z}}\cap\mathcal{S}^{(2)}_{\bm{z}}=\mathcal{S}^{(1)}_{\bm{z}}$
or $\mathcal{S}^{(2)}_{\bm{z}}$ for all $\bm{z}\in\mathcal{Z}$. This allows
$\mathcal{S}^{(1)}_{\bm{z}}\subseteq\mathcal{S}^{(2)}_{\bm{z}}$ for some
$\bm{z}$ but
$\mathcal{S}^{(1)}_{\bm{z}^{*}}\supseteq\mathcal{S}^{(2)}_{\bm{z}^{*}}$ for
another $\bm{z}^{*}\neq\bm{z}$. Furthermore, when $K=2$, condition (6) is
equivalent to assuming

$T_{\bm{z}}^{(1)}(\bm{Z},\bm{W})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}T_{\bm{z}}^{(2)}(\bm{Z},\bm{W})\mid\bm{Z}\in\mathcal{S}^{(1)}_{\bm{z}}\cap\mathcal{S}^{(2)}_{\bm{z}},\bm{W},\quad\forall\bm{z}\in\mathcal{Z}.$
(9)

Finally, the main conclusion (7) is equivalent to

$\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{W})\leq\alpha^{(1)},P^{(2)}(\bm{Z},\bm{W})\leq\alpha_{2}\mid\bm{Z}\in\mathcal{S}^{(1)}_{\bm{z}}\cup\mathcal{S}^{(2)}_{\bm{z}},\bm{W}\right\\}\leq\alpha^{(1)}\alpha^{(2)}$

for all $\bm{z}\in\mathcal{Z}$ and $\alpha^{(1)},\alpha^{(2)}\in[0,1]$.

###  3.2 Special cases

To better understand this general result, we first provide some special cases
and simple extensions. In particular, see Section 3.2.2 for a simple proof of
the $K=2$ case under stronger conditions. To simplify the exposition, we
assume that the test statistics satisfy
$T_{m}^{(k)}(\cdot,\cdot)=T^{(k)}(\cdot,\cdot)$ for all $k$ and $m$.

####  3.2.1 Independent treatment variables

We begin by noting that equation (5) is immediately satisfied if all the
randomization tests are unconditional (i.e.,
$\mathcal{S}_{\bm{z}}^{(k)}=\mathcal{Z}$ for all $\bm{z}$ and $k$).
Furthermore, the conditional independence (6) may be satisfied if the problem
structure allows the treatment $\bm{Z}$ to be decomposed into independent
components. This strategy is summarized in the proposition below.

######  Proposition 1.

The conditions (5) and (6) are satisfied for all $j,k\in[K],j\neq k$ if

  1. (i)

$\mathcal{S}_{\bm{z}}^{(k)}=\mathcal{Z}$ for all $k$ and $\bm{z}$; and

  2. (ii)

$T^{(k)}(\bm{Z},\bm{W})$ only depends on $\bm{Z}$ through
$\bm{Z}^{(k)}=h^{(k)}(\bm{Z})$ for all $k$ and
$\bm{Z}^{(j)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{Z}^{(k)}$ for all $j\neq k$.

This result can be easily extended to the case if the same partition is used
for all the conditional randomization tests, i.e.,
$\mathcal{R}^{(1)}=\dotsb=\mathcal{R}^{(K)}$.

One way to meet condition (ii) in Proposition 1 is through post-treatment
randomization Bates et al. (2020). They consider the natural “experiment” of
Mendelian inheritance, in which an unobserved variable $U$ satisfies

$\bm{Z}^{(1)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\dotsb\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{Z}^{(K)}\mid U.$

We can not directly apply Proposition 1 because $U$ is unobserved. However,
the joint distribution of $(U,\bm{Z})$ is given (or well approximated) by our
knowledge of Mendelian inheritance. The key idea in Bates et al. (2020) is to
construct a post-treatment random variable $G=g(\bm{Z},V)$ ($V$ is randomized
by the analyst), where the function $g(\cdot,\cdot)$ and the distribution of
$V$ are chosen such that $G$ has the same conditional distribution as $U$
given $\bm{Z}$. Then we immediately have

$\bm{Z}^{(1)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\dotsb\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\bm{Z}^{(K)}\mid G.$ (10)

Thus, we can apply Proposition 1 by treating $G$ as fixed. We refer readers to
Zhang and Zhao (2022, Section 3.5) for more details of how $G$ is sampled
using the known distribution of $(U,\bm{Z})$.

The nontrivial idea here is that when conditional independence between
treatment variables requires some unobserved variable $U$, we can instead
generate another variable $G$ that follow the same distribution of $U$ given
$Z$ and treat it as given. This proposal might seem astonishing at first, but
notice that the power of the test depends crucially on $G$ and on how well it
resembles $U$. If $G$ is similar to $U$, the randomized CRTs may have very
little power.

####  3.2.2 Sequential CRTs

Our next example is a sequence of CRTs, where each CRT conditions on the
randomness utilized by the preceding CRTs. Mathematically, we assume that the
CRTs are constructed in a nested way such that
$\mathcal{G}^{(1)}\subseteq\dotsb\subseteq\mathcal{G}^{(K)}$. By definition,
this is equivalent to using a nested sequence of partitions so that
$\mathcal{S}_{\bm{z}}^{(1)}\supseteq\dotsb\supseteq\mathcal{S}_{\bm{z}}^{(K)}$
for all $\bm{z}$. Under this construction, it is easy to see that the
refinement of two partitions is given by
$\underline{\mathcal{R}}^{\\{j,k\\}}=\mathcal{R}^{(\max(j,k))}$ and (5) is
immediately satisfied. Moreover, (6) is automatically satisfied if
$T^{(j)}(\bm{z},\bm{W})$ does not depend on $\bm{z}$ when
$\bm{z}\in\mathcal{S}_{m}^{(k)}$ for any $m$ and $k>j$. In other words,
$T^{(j)}(\bm{Z},\bm{W})$ is $\mathcal{G}^{(k)}$-measurable given $\bm{W}$.
Given this, the condition independence

$T^{(j)}(\bm{Z},\bm{W})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}T^{(k)}(\bm{Z},\bm{W})\mid\bm{Z}\in\mathcal{S}_{m}^{(k)},\bm{W},$

is immediately satisfied because $T^{(j)}(\bm{Z},\bm{W})$ is just a constant
given $\bm{Z}\in\mathcal{S}_{m}^{(k)}$.

To summarize, we have the following result.

######  Proposition 2.

The conditions (5) and (6) are satisfied for all $j,k\in[K],j\neq k$ if

  1. (i)

$\mathcal{S}_{\bm{z}}^{(1)}\supseteq\dotsb\supseteq\mathcal{S}_{\bm{z}}^{(K)}$
for all $\bm{z}\in\mathcal{Z}$; and

  2. (ii)

$T^{(j)}(\bm{z},\bm{W})$ does not depend on $\bm{z}$ when
$\bm{z}\in\mathcal{S}_{m}^{(k)}$ for all $m$ and $k>j$.

Under the stronger conditions in Proposition 2, the proof of Theorem 2 can be
greatly simplified. To illustrate this, consider the $K=2$ case. Let
$\psi^{(k)}(\bm{Z},\bm{W})=1_{\\{P^{(k)}(\bm{Z},\bm{W})\leq\alpha^{(k)}\\}}$
for $k=1,2$ denote the test functions. By Proposition 2(i),
$\mathcal{G}^{(1)}\subseteq\mathcal{G}^{(2)}$. By Proposition 2(ii) and the
definition of the p-value (1), $P^{(1)}(\bm{Z},\bm{w})$ and thus
$\psi^{(1)}(\bm{Z},\bm{w})$ are $\mathcal{G}^{(2)}$-measurable for any fixed
$\bm{w}$. Then by the law of iterated expectation, for any
$\bm{w}\in\mathcal{W}$,

$\displaystyle\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{w})\leq\alpha^{(1)},P^{(2)}(\bm{Z},\bm{w})\leq\alpha^{(2)}\mid\mathcal{G}^{(1)}\right\\}$
$\displaystyle=\ $
$\displaystyle\mathbb{E}\left\\{\psi^{(1)}(\bm{Z},\bm{w})\psi^{(2)}(\bm{Z},\bm{w})\mid\mathcal{G}^{(1)}\right\\}$
$\displaystyle=\ $
$\displaystyle\mathbb{E}\left\\{\mathbb{E}\left[\psi^{(1)}(\bm{Z},\bm{w})\psi^{(2)}(\bm{Z},\bm{w})\mid\mathcal{G}^{(2)}\right]\mid\mathcal{G}^{(1)}\right\\}$
$\displaystyle=\ $
$\displaystyle\mathbb{E}\left\\{\psi^{(1)}(\bm{Z},\bm{w})\mathbb{E}\left[\psi^{(2)}(\bm{Z},\bm{w})\mid\mathcal{G}^{(2)}\right]\mid\mathcal{G}^{(1)}\right\\}$
$\displaystyle\leq\ $
$\displaystyle\alpha^{(2)}\mathbb{E}\left\\{\psi^{(1)}(\bm{Z},\bm{w})\mid\mathcal{G}^{(1)}\right\\}$
$\displaystyle\leq\ $ $\displaystyle\alpha^{(1)}\alpha^{(2)}.$

Proposition 2 generalizes the knit product structure described by Rosenbaum
(2017), which essentially requires that permutation tests are constructed in a
sequential manner. Theorem 2 further extends this argument by considering the
conditioning structure in (5).

###  3.3 Proof outline

We now outline a proof of Theorem 2. We start with the following observation.
(All proofs can be found in the Appendix.)

######  Lemma 1.

Suppose (5) is satisfied. Then for any $\mathcal{J}\subseteq[K]$ and
$\mathcal{S},\mathcal{S}^{\prime}\in\mathcal{R}^{\mathcal{J}}$, the sets
$\mathcal{S}$ and $\mathcal{S}^{\prime}$ are either disjoint or nested, that
is,

$\mathcal{S}\cap\mathcal{S}^{\prime}\in\\{\varnothing,\mathcal{S},\mathcal{S}^{\prime}\\}.$

Furthermore, we have $\underline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$
and $\overline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$.

The sets in $\mathcal{R}^{[K]}$ can be partially ordered by set inclusion.
This induces a graphical structure on $\mathcal{R}^{[K]}$:

######  Definition 3.

The _Hasse diagram_ for
$\mathcal{R}^{[K]}=\\{\mathcal{S}_{\bm{z}}^{(k)}:\bm{z}\in\mathcal{Z},k\in[K]\\}$
is a graph where each node in the graph is a set in $\mathcal{R}^{[K]}$ and a
directed edge $\mathcal{S}\rightarrow\mathcal{S}^{\prime}$ exists between two
distinct nodes $\mathcal{S},\mathcal{S}^{\prime}\in\mathcal{R}^{[K]}$ if
$\mathcal{S}\supset\mathcal{S}^{\prime}$ and there is no
$\mathcal{S}^{\prime\prime}\in\mathcal{R}^{[K]}$ such that
$\mathcal{S}\supset\mathcal{S}^{\prime\prime}\supset\mathcal{S}^{\prime}$.

It is straightforward to show that all edges in the Hasse diagram for
$\mathcal{R}^{[K]}$ are directed and this graph has no cycles. Thus, the Hasse
diagram is a directed acyclic graph. For any node
$\mathcal{S}\in\mathcal{R}^{[K]}$ in this graph, we can further define its
parent set as
$\text{pa}(\mathcal{S})=\\{\mathcal{S}^{{}^{\prime}}\in\mathcal{R}^{[K]}:\mathcal{S}^{\prime}\to\mathcal{S}\\}$,
child set as
$\text{ch}(\mathcal{S})=\\{\mathcal{S}^{{}^{\prime}}\in\mathcal{R}^{[K]}:\mathcal{S}\to\mathcal{S}^{\prime}\\}$,
ancestor set as
$\text{an}(\mathcal{S})=\\{\mathcal{S}^{\prime}\in\mathcal{R}^{[K]}:\mathcal{S}^{\prime}\supset\mathcal{S}\\}$,
and descendant set as
$\text{de}(\mathcal{S})=\\{\mathcal{S}^{\prime}\in\mathcal{R}^{[K]}:\mathcal{S}\supset\mathcal{S}^{\prime}\\}$.

Notice that one conditioning set $\mathcal{S}\in\mathcal{R}^{[K]}$ can be used
in multiple CRTs. To fully characterize this structure, we introduce an
additional notation.

######  Definition 4.

For any $\mathcal{S}\in\mathcal{R}^{[K]}$, let
$\mathcal{K}(\mathcal{S})=\\{k\in[K]:\mathcal{S}\in\mathcal{R}^{(k)}\\}$ be
the collection of indices such that $\mathcal{S}$ is a conditioning set in the
corresponding test. Furthermore, for any collection of conditioning sets
$\mathcal{R}\subseteq\mathcal{R}^{[K]}$, denote
$\mathcal{K}(\mathcal{R})=\cup_{\mathcal{S}\in\mathcal{R}}\mathcal{K}(\mathcal{S})$.

######  Lemma 2.

Suppose (5) is satisfied. Then for any $\mathcal{S}\in\mathcal{R}^{[K]}$, we
have

  1. (i)

If $\text{ch}(\mathcal{S})\neq\varnothing$, then $\text{ch}(\mathcal{S})$ is a
partition of $\mathcal{S}$;

  2. (ii)

$\\{\mathcal{K}(\text{an}(\mathcal{S})),\mathcal{K}(\mathcal{S}),\mathcal{K}(\text{de}(\mathcal{S}))\\}$
forms a partition of $[K]$.

  3. (iii)

For any $\mathcal{S}^{\prime}\in\text{ch}(\mathcal{S})$,
$\mathcal{K}(\text{an}(\mathcal{S}^{\prime}))=\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\})$
and
$\mathcal{K}(\\{\mathcal{S}^{\prime}\\}\cup\text{de}(\mathcal{S}^{\prime}))=\mathcal{K}(\text{de}(\mathcal{S}))$.

Using Lemma 2, we can prove the following key lemma that establishes the
conditional independence between two p-values.

######  Lemma 3.

Suppose (5) and (6) are satisfied. Then for any
$\mathcal{S}\in\mathcal{R}^{[K]}$, $j\in\mathcal{K}(\mathcal{S})$, and
$k\in\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\})\setminus\\{j\\}$,
we have

$P^{(j)}(\bm{Z},\bm{W})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}P^{(k)}(\bm{Z},\bm{W})\mid\bm{Z}\in\mathcal{S},\bm{W}.$

Finally, we state a general result based on the above Hasse diagram, from
which Theorem 2 almost immediately follows.

######  Lemma 4.

Given conditions (5) and (6), we have, for any
$\mathcal{S}\in\mathcal{R}^{[K]}$,

$\begin{split}&\leavevmode\nobreak\
\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{W})\leq\alpha^{(k)},\dotsc,P^{(K)}(\bm{Z},\bm{W})\leq\alpha^{(K)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\\\
\leq&\leavevmode\nobreak\
\mathbb{P}\left\\{P^{(k)}(\bm{Z},\bm{W})\leq\alpha^{(k)}\text{ for
}k\in\mathcal{K}(\text{an}(\mathcal{S}))\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\\{\mathcal{S}\\}\cup\text{de}(\mathcal{S}))}\alpha_{j}.\end{split}$
(11)

##  4 Testing lagged treatment effect in stepped-wedge randomized trials

We next demonstrate the versatility of our theory using the stepped-wedge
randomized controlled trials, a widely used experimental design in medical and
policy research. Mixed-effects models are commonly used to analyze stepped-
wedge randomized trials (Hussey and Hughes, 2007; Hemming et al., 2018; Li et
al., 2021). However, it is well recognized that the statistical inference
tends to be biased (such as inflated type I error and poor confidence
coverage) when the model is misspecified (Thompson et al., 2017; Ji et al.,
2017). In light of this, unconditional randomization tests have been proposed
to test the null hypothese of no treatment effect whatsoever by several
authors (Ji et al., 2017; Wang and De Gruttola, 2017; Thompson et al., 2018;
Hughes et al., 2020). However, this approach cannot be used to test fine-
grained hypotheses such as those involving lagged treatment effect. Following
our theory, we develop a method (Algorithm 1) to construct nearly independent
CRTs for lagged treatment effects in Section 4.2. We also consider methods to
combine the tests in Section 4.3.

###  4.1 Hypotheses for lagged treatment effect

The stepped-wedge randomized controlled trial is a monotonic cross-over design
in which all units start out in the control and then cross over to the
treatment at staggered times (Gambia Hepatitis Study Group et al., 1987;
Hussey and Hughes, 2007). Figure 1 shows the treatment schedule for a typical
stepped-wedge trial; the name “stepped-wedge” refers to the wedge shape given
by the treated groups over subsequent crossover points. This design is also
known as _staggered adoption_ in econometrics (Abraham and Sun, 2021; Athey
and Imbens, 2022). Many stepped-wedge trials are cluster-randomized, i.e., the
units in the same cluster always have the same treatment status. For
simplicity, here we assume that the cross-over times are randomized at the
unit level. The same method below applies to cluster-randomized trials design
by considering aggregated cluster-level outcomes (Middleton and Aronow, 2015;
Thompson et al., 2018).

Figure 1: Stepped-wedge randomized trials. A group of $9$ units or clusters
switch from control to treatment at each time point and then remain exposed to
the treatment. This figure shows the treatment status of the units right
before the treatment assignment at each time point.

Consider a stepped-wedge trial on $N$ units over some evenly spaced time grid
$[T]=\\{1,\cdots,T\\}$. Let $N_{1},\dotsc,N_{T}$ be positive integers that
satisfy $\sum_{t=1}^{T}N_{t}=N$. The cross-overs can be represented by a
binary matrix $\bm{Z}\in\\{0,1\\}^{N\times T}$, where $Z_{it}=1$ indicates
that unit $i$ crosses over from control to treatment at time $t$. We assume
that the treatment assignment mechanism is given by the uniform distribution
over

$\mathcal{Z}=\left\\{\bm{z}\in\\{0,1\\}^{N\times
T}:\bm{z}\bm{1}=\bm{1},\bm{1}^{\top}\bm{z}=(N_{1},\cdots,N_{T})\right\\}.$
(12)

The set $\mathcal{Z}$ contains all binary matrices with exactly one element of
$1$ in each row and $N_{t}$ elements of $1$ in the $t$-th column, so
$|\mathcal{Z}|=N!/(N_{1}!\dotsb N_{T}!)$.

To construct multiple CRTs, a key observation is that this stepped-wedge trial
can be viewed as a sequentially randomized experiment. Let $\bm{Z}_{t}$ be the
$t$-th column of $\bm{Z}$. It is not difficult to show that the uniform
distribution over $\bm{Z}$ over $\mathcal{Z}$ is equivalent to randomly
assigning treatment to $N_{t}$ control units at time $t=1,\dotsc,T$. Denote
$\bm{z}_{[t]}=(\bm{z}_{1},\dotsc,\bm{z}_{t})\in\\{0,1\\}^{N\times t}$ and
$N_{[t]}=\sum_{s=1}^{t}N_{s}$. More precisely, we may decompose the uniform
distribution $\pi(\bm{z})$ over $\bm{z}\in\mathcal{Z}$ as

$\frac{N_{1}!\dotsb
N_{T}!}{N!}=\pi(\bm{z})=\pi_{1}(\bm{z}_{1})\pi_{2}(\bm{z}_{2}\mid\bm{z}_{1})\dotsb\pi_{T}(\bm{z}_{T}\mid\bm{z}_{[T-1]}),$

where

$\pi_{t}(\bm{z}_{t}\mid\bm{z}_{[t-1]})=\frac{N_{t}!(N-N_{[t]})!}{(N-N_{[t-1]})!}\leavevmode\nobreak\
\text{for all $\bm{z}_{t}\in\\{0,1\\}^{N}$ such that
$(\bm{1}_{N}-\bm{z}_{[t-1]}\bm{1}_{t-1})^{T}\bm{z}_{t}=N_{t}$}.$

After the cross-overs at time $t$, we assume that the experimenter immediately
measures the outcomes of all units, which are denote by
$\bm{Y}_{t}\in\mathbb{R}^{N}$. The corresponding vector of potential outcomes
under the intervention $\bm{Z}=\bm{z}$ is denoted by
$\bm{Y}_{t}(\bm{z})=(Y_{1t}(\bm{z}),\dotsc,Y_{nt}(\bm{z}))$. Let
$\bm{Y}=(\bm{Y}_{1},\dotsc,\bm{Y}_{T})$ be all the realized outcomes and
$\bm{Y}(\bm{z})$ be the collection of potential outcomes under treatment
assignment $\bm{z}\in\mathcal{Z}$. As before, we maintain the consistency
assumption that $\bm{Y}=\bm{Y}(\bm{Z})$.

Recall that a null hypothesis $H$ is fully sharp if its imputability mapping
$\mathcal{H}(\bm{z},\bm{z}^{*})=[N]$ for any $\bm{z},\bm{z}^{*}\in\mathcal{Z}$
(Section 2.1). An example of the fully sharp null hypothesis in this setting
is that the treatment has no effect whatsoever, i.e.,

$H:\bm{Y}(\bm{z})=\bm{Y}(\bm{z}^{*}),\quad\forall\bm{z},\bm{z}^{*}\in\mathcal{Z},i\in[N],\text{
and }t\in[T].$ (13)

This hypothesis can be tested by an unconditional randomization test in a
straightforward manner. Because $\bm{Z}$ is completely randomized, this is a
permutation test and has been previously studied in the literature (Ji et al.,
2017; Wang and De Gruttola, 2017; Hughes et al., 2020). The test statistic is
usually obtained by fitting some linear mixed-effects model. However, this
sharp null hypothesis is rather restrictive. In particular, the staggered
treatment assignment offers an opportunity to investigate how soon the
treatment takes effect, which is ignored in (13). This is what we will
consider next.

To make the problem more tractable, we assume no interference and no
anticipation effect (Athey and Imbens, 2022) in the following sense:

######  Assumption 2.

For all $i\in[N]$ and $t\in[T]$, $Y_{it}(\bm{z})$ only depends on $\bm{z}$
through $\bm{z}_{i,[t]}$, the treatment history of unit $i$.

The hypothesis we consider assumes that the lag $l$ treatment effect is equal
to a constant $\tau_{l}$ ($l$ is a fixed integer between $0$ and $T-1$). More
precisely, consider the following sequence of partially sharp null hypotheses,

$H^{(t)}:Y_{i,t+l}((\bm{0}_{t-1},1,\bm{0}_{l}))-Y_{i,t+l}(\bm{0}_{t+l})=\tau_{l},\
\forall i\in[N],$ (14)

for $t=1,\dotsc,T-l$. The hypothesis $H^{(t)}$ states that the treatment that
cross-overs at time $t$ has a constant effect $\tau_{l}$ on the outcome at
time $t+l$ for all the units. We are interested in testing the intersection of
$H^{(1)},\dotsc,H^{(T-l)}$. Note that this intersection is not the same as the
fully sharp hypothesis in (13) as it only concerns the lag-$l$ effect.

###  4.2 Multiple conditional randomization tests for testing lagged treatment
effect

We first note that $H^{(t)}$ is partially sharp. In fact, if we denote the
imputability mapping of $H^{(t)}$ as $\mathcal{H}^{(t)}(\bm{z},\bm{z}^{*})$
(see Definition 1), then unit has imputable potential outcome (i.e.,
$i\in\mathcal{H}^{(t)}(\bm{z},\bm{z}^{*})$) if and only if $i$ crosses over
either at time $t$ or after $t+l$ under both assignments $\bm{z}$ and
$\bm{z}^{*}$. Therefore,

$\mathcal{H}^{(t)}(\bm{z},\bm{z}^{*})=\left\\{i\in[N]:\bm{z}_{i,[t+l]},\bm{z}^{*}_{i,[t+l]}\in\left\\{(\bm{0}_{t-1},1,\bm{0}_{l}),\bm{0}_{t+l}\right\\}\right\\}.$

This motivates us to condition on the units that have imputable potential
outcomes. We let the CRT for $H^{(t)}$ use the conditioning set

$\mathcal{S}_{\bm{Z}}^{(t)}=\\{\bm{z}^{*}\in\mathcal{Z}:g^{(t)}(\bm{z}^{*})=g^{(t)}(\bm{Z})\\}$

where $g^{(t)}(\bm{Z})$ is a subset of units that cross over at time $t$ or
after $t+l$,

$g^{(t)}(\bm{Z})=\left\\{i\in[N]:\bm{Z}_{i,[t+l]}=(\bm{0}_{t-1},1,\bm{0}_{l})\text{
or }\bm{0}_{t+l}\right\\}.$ (15)

The conditioning sets $\mathcal{S}_{\bm{z}}^{(t)},\bm{z}\in\mathcal{Z}$, form
a partition
$\mathcal{R}^{(k)}=\\{\mathcal{S}_{\bm{z}}^{(t)}:\bm{z}\in\mathcal{Z}\\}$ of
$\mathcal{Z}$ based on the value of $g^{(t)}$. The function $g$ is often known
as an _exposure mapping_ in the causal inference literature (see Zhang and
Zhao, 2022, Section 4 and references therein). We may then compute the p-value
(1) since both outcomes in (14) are either observed or imputable. Since the
treatment assignment is completely randomized, the CRT is essentially a
permutation test with $\mathcal{S}_{\bm{Z}}$ consisting of all the possible
permutations of units that cross over at time $t$ and after $t+l$.

However, a careful examination reveals that the conditioning sets
$\mathcal{S}_{\bm{Z}}^{(1)},\dotsc,\mathcal{S}_{\bm{Z}}^{(T-l)}$ fail to
satisfy the nested condition (5) in Theorem 2 when lag $l\geq 1$. Consider the
case of $l=1$. The $t$-th CRT (i.e. CRT $t$) is a permutation test that
compares the time $t+1$ outcome of units that cross over at $t$ with the units
that cross over after $t+1$: (i) The first CRT compares units that cross over
at time $1$ and units that cross over at time $3,4,\dotsc$; (ii) The second
CRT compares units that cross over at time $2$ and units that cross over at
time $4,5\dotsc$; and so on (see Figure 2(a) for an illustration). It is easy
to see that the corresponding conditioning sets are not nested. For example,
the conditioning set of CRT 1 includes cross over (i.e. starting the
treatment) at time 3 but not time 2, while the reverse is true for the
conditioning set of CRT 2. Independence of the CRTs thus cannot be
established.

The discussion above also suggests that no single CRT can be constructed to
test the intersection of the hypotheses $H^{(t)},t=1,\dotsc,T$ at multiple
times. For instance, suppose we would like to test $H^{(1)}\cap H^{(2)}$ with
lag $l=1$ by using potential outcomes at $t=2$ and $3$. Notice that the
potential outcome at time $3$, $Y_{i3}(\bm{Z}^{*})$ is not imputable when
$Z_{i}=1$ and $Z^{*}_{i}\neq Z_{i}$ because the lag is $l=1$ (see CRT2 in
Figure 2(a)), so units treated at time $1$ must be excluded from the test (by
conditioning on $\\{i:Z_{i}=1\\}$). The dilemma, however, is that those are
precisely the treated units for outcomes measured at time $2$.

(a) Non-nested CRTs.

(b) Nested CRTs.

Figure 2: Illustration of the CRTs for the lagged treatment effect in a
stepped-wedge randomized trial with $T=8$ time points (lag $l=1$).

To resolve this problem, we can gain some ideas by considering the $l=0$ case.
The conditioning variable $g^{(t)}(\bm{Z})$ in (15) actually induces nested
CRTs for $l=0$, as follows: (i) The first CRT compares units that cross over
at time $1$ and units that cross over at time $2,3,\dotsc$; (ii) The second
CRT compares units that cross over at time $2$ and units that cross over at
time $3,4\dotsc$; and so on. Intuitively, the first CRT utilizes the
randomness in $\bm{Z}_{1}$, the second CRT utilizes the randomness in
$\bm{Z}_{2}$ given $\bm{Z}_{1}$, and so on. It is easy to see that the
conditions in Proposition 2 are satisfied if the $t$-th CRT only depends on
$\bm{Z}$ through $\bm{Z}_{t}$ given $\bm{Z}_{1},\dotsc,\bm{Z}_{t-1}$.

Algorithm 1 Multiple conditional randomization tests (MCRTs) for testing
lag-$l$ treatment effect in stepped-wedge randomized controlled trials

1:Input: Number of units $N$, Number of time steps $T$, Time lag $l$, Outcomes
$\bm{Y}=(Y_{it}:i\in[N],t\in[T])$, Treatment assignments
$\bm{Z}=(Z_{it}:i\in[N],t\in[T])$, Test statistics $T$.

2:Initialization: $J\leftarrow\min(l+1,T-l-1)$ and
$\mathcal{I}_{t}\in\\{i\in[N]:Z_{it}=1\\},\forall t\in[T]$

3:for $j\in[J]$ do $\triangleright$ Divide the time steps $[T]$ into $J$
subsets

4: $t\leftarrow j$, $\mathcal{C}_{j}\leftarrow\\{t\\}$

5: while $t+l+1\leq T$ do

6: $t\leftarrow t+l+1$, $\mathcal{C}_{j}\leftarrow\mathcal{C}_{j}\cup\\{t\\}$

7: end while

8:end for

9:for $j\in[J]$ do

10: for $k\in\mathcal{C}_{j}$ do $\triangleright$ Define the $k$-th
permutation test

11: $\mathcal{T}^{(k)}\leftarrow\big{\\{}t\in\mathcal{C}_{j}:t\geq k\big{\\}}$

12: $\mathcal{I}_{1}^{(k)}\leftarrow\mathcal{I}_{k}$,
$\bm{Y}_{\text{treated}}^{(k)}\leftarrow\big{\\{}Y_{i,k+l}:i\in\mathcal{I}_{1}^{(k)}\big{\\}}$
$\triangleright$ Create a treated group

13:
$\mathcal{I}_{0}^{(k)}\leftarrow(\bigcup_{t\in\mathcal{T}^{(k)}}\mathcal{I}_{t}\big{)}\setminus\mathcal{I}_{k}$,
$\bm{Y}_{\text{control}}^{(k)}\leftarrow\big{\\{}Y_{i,k+l}:i\in\mathcal{I}_{0}^{(k)}\big{\\}}$
$\triangleright$ Create a control group

14: $P^{(k)}(\bm{Z},\bm{Y})\leftarrow\text{Permutation
Test}\big{(}\bm{Y}_{\text{treated}}^{(k)},\bm{Y}_{\text{control}}^{(k)};T\big{)}$

15: end for

16:end for

17:Output: P-values $\big{\\{}P^{(k)}:=P^{(k)}(\bm{Z},\bm{Y})\big{\\}}$

This inspires us to modify the tests for the $l=1$ case as follows: (i) The
first CRT compares units that cross over at time $1$ and units that cross over
at time $3,5,\dotsc$; (ii) The second CRT compares units that cross over at
time $2$ and units that cross over at time $4,6,\dotsc$; (iii) The third CRT
compares units that cross over at time $3$ and units that cross over at time
$5,7,\dotsc$; and so on (see Figure 2(b) for an illustration). By considering
smaller conditioning sets, the odd tests become a sequence of nested CRTs and
so do the even tests. Thus, the nesting condition (5) in Theorem 2 holds.
Moreover, because all the CRTs furhter condition on

$g(\bm{Z})=\left\\{i\in[N]:Z_{i,t}=1\text{ for some odd $t$}\right\\},$

it is not difficult to show that

$(\bm{Z}_{1},\bm{Z}_{3},\dotsc)\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}(\bm{Z}_{2},\bm{Z}_{4},\dotsc)\mid g(\bm{Z}).$
(16)

For every $t$, the $t$-th CRT only depends on $\bm{Z}$ through $\bm{Z}_{t}$
and uses the randomization distribution of $\bm{Z}_{t}$ given $\bm{Z}_{[t-1]}$
and $g(\bm{Z})$. Then it is straightforward to verify that every pair of odd
tests are conditionally independent and thus satisfy the condition (6) in
Theorem 2. Similarly, (6) also holds for the even tests. Finally, any pair of
odd and even tests satisfy (6) due to (16). Therefore, the p-values of the
modified CRTs satisfy (8) and can be combined using the standard global
testing methods such as Fisher’s combined probability test (Fisher, 1925).

The main idea in the discussion above is that we can further restrict the
conditioning sets when the most obvious conditioning sets are not nested. This
argument can be easily extended to a longer time lag $l$. When $l=1$, we have
divided $[T]$ into two subsets of cross-over times (odd and even). To test the
lag $l$ effect for $l>1$, we can simply increase the gap and divide $[T]$ into
disjoint subsets:
$\mathcal{C}_{1}=\\{1,l+2,2l+3,\dotsc\\},\mathcal{C}_{2}=\\{2,l+3,2l+4,\dotsc\\},\dotsc$.
A formal algorithm for general $l$ is given in Algorithm 1 and will be
referred to as multiple conditional randomization tests (MCRTs) in what
follows.

###  4.3 Method of combining p-values

A remaining question is how to combine the permutation tests obtained from
Algorithm 1 (MCRTs) efficiently. Heard and Rubin-Delanchy (2018) compared
several p-value combination methods in the literature. By recasting them as
likelihood ratio tests, they demonstrated that the power of a combiner
crucially depends on distribution of the p-values under the alternative
hypotheses. In large samples, the behaviour of permutation tests is well
studied in the literature. Lehmann and Romano (2006, Theorem 15.2.3) showed
that if the test statistics in a permutation test converges in distribution,
the permutation distribution will converge to the same limiting distribution
in probability. These results form the basis of our investigation of combining
CRTs.

Suppose that $K$ permutation p-values $P^{(1)},\dotsc,P^{(K)}$ are obtained
from from Algorithm 1. Suppose the $k$th test statistic, when scaled by
$\sqrt{N}$, is asymptotically normal so that it converges in distribution to
$T_{\infty}^{(k)}\sim\mathcal{N}(\tau_{l},V_{\infty}^{(k)})$ under the null
hypothesis ($\tau_{l}=0$) and some local alternative hypothesis
($\tau_{l}=h/\sqrt{N}$) indexed by $h$. The limiting distributions have the
same mean but different variances. Suppose that we define the p-value
likelihood ratio as the product of the alternative p-value distributions for
all $k$ divided by the product of the null p-value distributions for all $k$.
Since the permutation p-values are standard uniform, it is easy to verify that
the logarithm of this p-value likelihood ratio is proportion to a weighted sum
of $Z$-scores, $T_{\infty}=\sum_{k=1}^{K}w_{\infty}^{(k)}\Phi^{-1}(P^{(k)})$.
This motivates a weighted Stouffer’s method (Stouffer et al., 1949) for
combing the p-values. The weights are non-negative and sum to one, thus
$T_{\infty}\sim\mathcal{N}(0,1)$ if the p-values are independent. Then we can
reject $\bigcap_{k\in[K]}H_{0}^{(k)}$ if $\Phi(T_{\infty})\leq\alpha$. This
test has been shown to be uniformly most powerful for all $h$ in the normal
location problem (Heard and Rubin-Delanchy, 2018).

To formalize the discussion above, next we consider the difference-in-means
statistics as an example. Following the notation in Algorithm 1 (see lines 6
and 11), the treated group $\mathcal{I}_{1}^{(k)}$ for the $k$th CRT crosses
over at time $k$ and the control group $\mathcal{I}_{0}^{(k)}$ crosses over
time $t\in\mathcal{T}^{(k)}\setminus\\{k\\}$, where
$\mathcal{T}^{(k)}\setminus\\{k\\}$ is a subset of $\mathcal{C}_{j}$ that
collects some time points after $k+l$. Let
$N_{1}^{(k)}=|\mathcal{I}_{1}^{(k)}|$, $N_{0}^{(k)}=|\mathcal{I}_{0}^{(k)}|$
and $N^{(k)}=N_{0}^{(k)}+N_{1}^{(k)}$ Let $A_{i}^{(k)}=1\text{ or }0$ denote
unit $i\in\mathcal{I}_{0}^{(k)}\cup\mathcal{I}_{1}^{(k)}$ starting the
treatment at time $k$ or after $k+l$, i.e.,
$\bm{Z}_{i,k+l}=(\bm{0}_{k-1},1,\bm{0}_{l})\text{ or }\bm{0}_{k+l}$. To
simplify the exposition, we use the abbreviations $Y_{i}^{(k)}=Y_{i,k+l}$,
$Y_{i}^{(k)}(1)=Y_{i,k+l}(\bm{0}_{k-1},1,\bm{0}_{l})$ and
$Y_{i}^{(k)}(0)=Y_{i,k+l}(\bm{0}_{k+l})$ in the results below.

We further assume that the $N$ units are i.i.d draws from a super-population
model, and every unit’s potential outcome $Y_{i}^{(k)}(a)$ is a random copy of
some generic $Y^{(k)}(a)$ for $a\in\\{0,1\\}$ and $k\in{K}$.

######  Proposition 3.

Suppose that the $K$ permutation tests in Algorithm 1 use the difference-in-
means statistics. We assume that for every $k\in[K]$ and $a\in\\{0,1\\}$,

  1. (i)

$Y^{(k)}(a)$ has a finite and nonzero variance;

  2. (ii)

$N_{1}^{(k)}/N$ and $N_{0}^{(k)}/N$ are fixed as $N\rightarrow\infty$;

The weighted $Z$-score combiner $T_{\infty}$ is then given with weights
$w_{\infty}^{(k)}=\sqrt{\frac{\Lambda^{(k)}}{\sum_{j=1}^{K}\Lambda^{(j)}}},$
where

$\Lambda^{(k)}=\left(\frac{N}{N_{0}^{(k)}}\operatorname{Var}\big{[}Y^{(k)}(1)\big{]}+\frac{N}{N_{1}^{(k)}}\operatorname{Var}\big{[}Y^{(k)}(0)\big{]}\right)^{-1}$
(17)

is the inverse of the asymptotic variance $V_{\infty}^{(k)}$ of the statistics
in the $k$-th test.

The weights in (17) are obtained by inverting the asymptotic variance of the
difference-in-means statistics; see Imbens and Rubin (2015, Section 6.4) and
Li and Ding (2017). We can estimate
$\operatorname{Var}\big{[}Y^{(k)}(a)\big{]}$ in $\Lambda^{(k)}$ consistently
using the sample variance of $Y_{i}^{(k)},\forall i\in\mathcal{I}_{a}^{(k)}$.
We denote the estimator of $\Lambda^{(k)}$ by $\hat{\Lambda}^{(k)}$. The
empirical version of $T_{\infty}$ is given by

$\hat{T}=\sum_{k=1}^{K}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})\quad\text{where}\quad\hat{w}^{(k)}=\sqrt{\frac{\hat{\Lambda}^{(k)}}{\sum_{j=1}^{K}\hat{\Lambda}^{(j)}}}.$
(18)

######  Proposition 4.

Suppose that the p-values $P^{(1)},\dotsc,P^{(K)}$ from Algorithm 1 are valid,
the combined p-value $\hat{P}(\bm{Z},\bm{Y})=\Phi(\hat{T})$ is also valid such
that under the null hypotheses,

$\mathbb{P}\\{\hat{P}(\bm{Z},\bm{Y})\leq\alpha\\}\leq\alpha.$

In general, maximizing the statistical power in combining multiple CRTs is a
independent task on top of choosing the most efficient test statistics for the
permutation tests. Different test statistics may have different asymptotical
distributions and the optimal p-value combiners (if exist) may be different.
Nonetheless, a key insight from the discussion above is that we should weight
the CRTs appropriately, often according to their sample size.

##  5 Numerical experiments

Next, we examine the empirical performance of methods proposed in Section 4 by
comparing the power of them with Bonferroni’s Method (Simulation I), checking
their validity when there are unit-by-time interactions (Simulation II), and
investigating several real trial datasets. We let MCRTs+F and MCRTs+Z denote
the combination of MCRTs with Fisher’s combination method (Fisher, 1925) and
the weighted Z-score method introduced above, respectively. In both MCRTs+F
and MCRTs+Z, the test statistic is the simple difference-in-means.

###  5.1 Simulation I: Size and power

As illustrated in Figure 2, the non-nested CRTs have larger control groups
than our nested CRTs created in MCRTs. However, the non-nested CRTs are not
independent and there are limited ways to combine them. The most widely used
method in this case is Bonferroni’s method that rejects the intersection of
$K$ hypotheses if the smallest p-value is less than $\alpha/K$. One remaining
question is whether the reduced sample size in MCRTs can indeed be compensated
by a better p-value combiner. Simulation I is designed to answer this question
empirically by investigating varying the sample size $N$, trial length $T$,
time lag $l$, and effect sizes $\tau_{l}$, respectively.

In this simulation, we fixed $N_{1}=\dotsc=N_{T}=N/T$ and
$N_{T}=N-\sum_{t=1}^{T-1}N_{t}$. The treatment assignment $\bm{Z}$ was
randomly drawn from $\mathcal{Z}$ given in (12). Let $A_{i}$ to denote the
treatment starting time of unit $i$, i.e., the index of the non-zero entry of
$Z_{i}$. Outcomes were generated by a linear mixed-effects model,

$Y_{it}=\mu_{i}+0.5(X_{i}+t)+\sum_{l=0}^{T-1}1_{\\{A_{i}-t=l\\}}\tau_{l}+\epsilon_{it},\leavevmode\nobreak\
i=1,\dotsc,N,\leavevmode\nobreak\ t=0,\dotsc,T,$

where $\mu_{i}\sim\mathcal{N}(0,0.25)$, $X_{i}\sim\mathcal{N}(0,0.25)$ and
$\epsilon_{it}\sim\mathcal{N}(0,0.1)$. This assumes that the baseline outcome
$Y_{i0}$ is measured, which is not uncommon in real clinical trials. Basic
parameters were varied in the following ranges: (i) number of units
$N\in\\{100,200,300,400,500\\}$; (ii) number of time steps
$T\in\\{4,6,8,10,12\\}$; (iii) time lag $l\in\\{0,1,2,3,4\\}$; effect size
$\tau_{l}\in\\{0.01,0.02,0.03,0.04,0.05\\}$. Empirical performance of the
methods was examined when one of $N$, $T$, and $l$ is changed while the other
two are fixed at the median of their ranges, respectively. For example, we
increased $N$ from 100 to 500 while keeping $T=8$ and $l=2$. In these
simulations, $\tau_{l}$ was set to $0$ and $0.03$ to investigate type I error
and power of the methods, respectively. One more simulation was created to
study the power as the effect size $\tau_{l}$ varies, in which we keep the
first three basic parameters at their median values $(N=300,T=8,l=2)$ and
increase $\tau_{l}$ from 0.01 to 0.05. For all permutation tests, difference-
in-means was used as the test statistic and the number of permutations was
fixed at 1000.

Figure 3: Performance of MCRTs+F, MCRTs+Z and Bonferroni’s method: type I
error rates and powers in testing lagged effects at five different numbers of
units, numbers of time steps, time lags and effect sizes. The results were
averaged over 1000 independent runs.

The upper panels of Figure 3 show that all the methods control the type I
errors at any number of units, time steps and time lags. The lower panel shows
that our methods are more powerful than Bonferroni’s method in all the
simulations. MCRTs+Z is slightly more powerful than MCRTs+F in all
experiments. This shows that the weighted z-score is more effective than
Fisher’s combination method for MCRTs. Panels (d) and (g) show that our
methods outperform Bonferroni’s method by a wider margin as the sample size or
the effect size increases. Panel (e) establishes the same observation for
trial length, which can be explained by the fixed sample size (so fewer units
start treatment at each step as trial length increases). This is not an ideal
scenario for applying Bonferroni’s method, as the individual tests have
diminishing power. Finally, panel (f) shows that all the methods have smaller
power as the lag size increases, and our methods are particularly powerful
when the time lag is small. These results are due to facts. First, there are
fewer permutation tests available for larger time lags. Moreover, only
including the units that are treated after a large time lag, the control
groups in the permutation tests are small.

Overall, these simulations supports the conclusion that that the sample
splitting in MCRTs is a worthy sacrifice as a more powerful p-value combiner
can be applied. Note that every outcome relevant to lag $l$ effect is still
used in at least one of our permutation tests in MCRTs. So the reduced sample
size in MCRTs may not be as damaging as it might first appear.

###  5.2 Simulation II: Overcome misspecification of mixed-effects models

Simulation II is designed to investigate the finite-sample properties of
confidence intervals obtained by inverting conditional randomization tests
(see Appendix B for more details). In particular, we are interested in
comparing its efficiency and robustness with mixed-effects models for
estimating lagged treatment effect.

Figure 4: Performance of MCRTs+F, MCRTs+Z and the mixed-effects model:
coverage rates and lengths of confidence intervals (CIs) under the outcome
generating processes with no interaction effect and with three different types
of covariate-and-time interactions. The results are averaged over 1000
independent runs.

In Simulation II, the treatment generating process was kept the same as
Simulation I. Simulation parameters are set as $N=200$, $T=8$ and lagged
effects $(\tau_{0},\dotsc,\tau_{7})=(0.1,0.3,0.6,0.4,0.2,0,0,0)$. The
treatment effect is gradually realized and then decayed to 0 over time. For
every $i\in[N]$ and $t=0,1,\dotsc,T$, we generate the outcome $Y_{it}$ using a
mixed-effects model,

$Y_{it}=\mu_{i}+0.5(1-0.1\cdot 1_{\\{m\neq
0\\}})(X_{i}+t)+0.1f_{m}(X_{i}+t)+\sum_{l=0}^{7}1_{\\{A_{i}-t=l\\}}\tau_{l}+\epsilon_{it},$

where $\mu_{i}\sim\mathcal{N}(0,0.25)$, $X_{i}\sim\mathcal{N}(0,0.25)$,
$\epsilon_{it}\sim\mathcal{N}(0,0.1)$ and the unit-by-time interaction is
given by

$f_{m}(X_{i}+t)=\begin{cases}0,&\quad\text{if }m=0,\\\
(X_{i}+t)^{2},&\quad\text{if }m=1,\\\ 2\exp[(X_{i}+t)/2],&\quad\text{if
}m=2,\\\ 5\tanh(X_{i}+t),&\quad\text{if }m=3,\\\ \end{cases}$

which correspond to no, quadratic, exponential and hyperbolic tangent
interactions.

We tasked our methods (MCRTs+F and MCRTs+Z) and a mixed-effects model
described below to construct valid 90%-confidence interval (CIs) for the
lagged effects $\tau_{0},\dotsc\tau_{4}.$ Our methods were implemented in
exactly the same way as in Simulation I except that the permutation tests were
inverted to obtain interval estimators. The mixed-effects model takes the form

$Y_{it}=\beta_{0i}+\beta_{1}x_{i}+\beta_{2}t+\sum_{l=0}^{7}\xi_{l}1_{\\{A_{i}-t=l\\}}+\epsilon_{it},$
(19)

where $\beta_{0i}$ is a random unit effect, $\beta_{1},\beta_{2}$ are fixed
effects and $\xi_{0},\dotsc,\xi_{7}$ are lagged effects, and were fitted using
the R-package lme4 (Bates et al., 2015). Recently, Kenny et al. (2021)
proposed mixed-effects models that can leverage the shapes of time-varying
treatment effects. Their models are given by specifying some parametric effect
curves with the help of basis functions (e.g. cubic spline). Besides the unit-
by-time interaction, the model (19) is already correctly specified for
modelling the treatment effects and other parts of the outcome generating
process above. Excluding some of the lagged effect parameters
$\xi_{5},\dotsc,\xi_{7}$ from the model may change the effect estimates but
not the validity of its CIs.

As the unit-by-time interactions are unknown and fully specifying them would
render the model unidentifiable, they are typically not considered in mixed-
effect models. If the stepp-wedge trial is randomized at the cluster levels,
one can introduce cluster-by-time interaction effects in a mixed-effects
model; see (Ji et al., 2017, Equation 3.1) as an example. However, the
cluster-by-time interactions are not exactly the same as the interaction
between time and some covariates of the units. It depends on if the
interaction varies within each cluster and the clusters are defined by the
covariates which interact with time.

The results of Simulation II are reported in Figure 4. The upper panels show
that the CIs of the mixed-effects model only achieve the target coverage rate
of $90\%$ for all the lagged effects when there is no unit-by-time interaction
(panel a). By contrast, the CIs of our methods fulfil the target coverage rate
of $90\%$ under all scenarios. The lower panels show that the valid CIs given
by our methods have reasonable lengths between 0.05 and 0.10 for covering the
true lagged effects $(\tau_{0},\dotsc,\tau_{4})=(0.1,0.3,0.6,0.4,0.2)$. We
also note from panel (e) that when the linear mixed-effects model is specified
correctly, it gives valid CIs that are narrower than our nonparametric tests.

###  5.3 Real data applications

We next demonstrate an application of MCRTs+Z in comparison to mixed-effects
models (with and without time effect parameters) on real data collected from
two stepped-wedge randomized trials. The results are reported as
90%-confidence intervals (CIs) of lagged effects in Figure 5.

Figure 5: Effect estimates from MCRTs+Z and mixed-effects models with and
without time effect parameters: 90%-confidence intervals (CIs) of lagged
effects on real data collected from four different stepped-wedge randomized
trials.

Trial I was conducted to examine if a chronic care model was more effective
than usual care in the Netherlands between 2010 and 2012 (Muntinga et al.,
2012). The study consisteds of 1,147 frail older adults in 35 primary care
practices. The primary outcome in the dataset was quality of life measured by
a mental health component score (MCS) in a 12-item Short Form questionnaire
(SF-12). The outcome for each adult was measured at a baseline time point and
every 6 months over the study. The trial randomly assigned some practices to
start the chronic care model every 6 months. We treated each practice as an
experimental unit. We used MCRTs+Z and two mixed-effects model to estimate lag
$0,1$ and $2$ effects on the average outcome of each practice,
$Y_{it},i=1,\dotsc,35$ and $t=0,1,\dotsc,4$. In panel (a) of Figure 5, the
mixed-effects model without using a time effect parameter would conclude that
the chronic care model improved the quality of life significantly. However, as
noted by Twisk (2021, Section 6.3.1), this model fails to take into account
the increase in the quality of life over time, irrespective of the
intervention. In contrast, adding a time effect parameter to the mixed-effects
model and applying MCRTs+Z produced similar CIs, both showing that the effect
of the chronic care model on quality of life is not statistically significant.

Trial II was a stepped-wedge trial performed over the pain clinics of 17
hospitals to estimate the effectiveness of an intervention in reducing chronic
pain for patients (Twisk, 2021, Chapter 6.4). The outcome was a pain score
from 1 (least pain) and 5 (most pain) and six measurements were taken over
time: one at the baseline and five more equally spaced over a period of twenty
weeks. After each measurement, intervention was started in a few randomly
selected untreated hospitals. We considered each hospital as an experimental
unit and estimated the lag $0,1,2$ and $3$ effects using hospital-level
average outcomes. Panel (b) of Figure 5 shows that the CIs of the mixed-effect
models with and without a time effect do not overlap. The CIs of MCRTs+Z lie
roughly between those, showing that the effectiveness of the intervention
appear to increase over time.

##  6 Discussion

In this article, we have developed an overarching theory for understanding
multiple conditional randomization tests and applied it to testing lagged
treatment effects in stepped-wedge randomized trials. A particularly
attractive feature of randomization tests is that they do not require
specifying any models and are thus robust to model misspecification, which is
confirmed in our numerical experiments.

The main theoretical result (Theorem 2) places no restriction on the
assignment mechanism, whether there is interference, or the nature of the
intervention. It provides a general principle for constructing nearly
independent CRTs in complex experiments. More user-friendly corollaries were
provided in Section 3.2. Sample splitting is another useful technique when the
most obvious CRTs do not satisfy these conditions, which is employed by our
method for the stepped-wedge design in Section 4.2.

Our theory assumes a randomised experiment, and it is interesting to study if
the same theory also applies to randomization inference for observational
studies. Matching by observed confounders is often used to mimic an randomized
experiment using observational study. The complication here is that the
matched sets are often a random subset of units that depend on the treatment
assignment, so our theory cannot directly apply because it assumes
$\mathcal{Z}$ is fixed. This subtle issue is usually ignored in the
randomization inference literature, which typically treats the matched
observational sample as a block-randomized experiment (Rosenbaum, 2002).

Rosenbaum (2010, 2017) developed a concept called “evidence factors” for
observational studies, which are essentially different pieces of evidence
regarding different aspects of an abstract causal hypothesis. Roughly
speaking, each evidence factor is characterized by an p-value upper bound
obtained from a sensitivity analyses for a different (partially) sharp null
hypothesis. Rosenbaum’s key observation is that in certain problems, these
evidence factors are nearly independent in the sense that the p-value upper
bounds satisfy Proposition 2. The “knit-product” structure identified in
Rosenbaum (2017) motivated the consideration of sequential CRTs in Section
3.2.2 and is essentially a special case of Proposition 2 when the treatment
assignment mechanism $\pi(\cdot)$ is uniform over $\mathcal{Z}$ (so the CRT is
a permutation test). Based on our more general theory, new and more powerful
evidence factors can potentially be developed.

In a parallel line of research, it has been shown that a randomization test
wit a carefully consutrcted test statistic is asymptotically valid for testing
average treatment effects (Cohen and Fogarty, 2022; Fogarty, 2021; Zhao and
Ding, 2021), “bounded” null hypotheses (Caughey et al., 2021), or quantiles of
individual treatment effects (Caughey et al., 2021). These techniques can be
used in conjunction with our theory to construct nearly independent CRTs
beyond (partially) sharp null hypotheses.

## References

  * Abraham and Sun (2021) Abraham, S. and L. Sun (2021).  Estimating dynamic treatment effects in event studies with heterogeneous treatment effects.  Journal of Econometrics 225(2), 175–199.
  * Athey et al. (2018) Athey, S., D. Eckles, and G. W. Imbens (2018).  Exact p-values for network interference.  Journal of the American Statistical Association 113(521), 230–240.
  * Athey and Imbens (2022) Athey, S. and G. W. Imbens (2022).  Design-based analysis in difference-in-differences settings with staggered adoption.  Journal of Econometrics 226(1), 62–79.
  * Basse et al. (2019) Basse, G., A. Feller, and P. Toulis (2019).  Randomization tests of causal effects under interference.  Biometrika 106(2), 487–494.
  * Bates et al. (2015) Bates, D., M. Mächler, B. Bolker, and S. Walker (2015).  Fitting linear mixed-effects models using lme4.  Journal of Statistical Software 67(1), 1–48.
  * Bates et al. (2020) Bates, S., M. Sesia, C. Sabatti, and E. Candès (2020).  Causal inference in genetic trio studies.  Proceedings of the National Academy of Sciences 117(39), 24117–24126.
  * Brown and Lilford (2006) Brown, C. A. and R. J. Lilford (2006).  The stepped wedge trial design: a systematic review.  BMC Medical Research Methodology 6(1), 1–9.
  * Caughey et al. (2021) Caughey, D., A. Dafoe, X. Li, and L. Miratrix (2021).  Randomization inference beyond the sharp null: Bounded null hypotheses and quantiles of individual treatment effects.  arXiv preprint arXiv:2101.09195.
  * Cohen and Fogarty (2022) Cohen, P. L. and C. B. Fogarty (2022).  Gaussian prepivoting for finite population causal inference.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), to appear.
  * Ernst (2004) Ernst, M. D. (2004).  Permutation methods: a basis for exact inference.  Statistical Science 19(4), 676–685.
  * Fisher (1925) Fisher, R. A. (1925).  Statistical methods for research workers. Oliver and Boyd, Edinburgh and London.
  * Fisher (1935) Fisher, R. A. (1935).  Design of Experiments.  Edinburgh: Oliver and Boyd.
  * Fogarty (2021) Fogarty, C. B. (2021).  Prepivoted permutation tests.  arXiv preprint arXiv:2102.04423.
  * Gambia Hepatitis Study Group et al. (1987) Gambia Hepatitis Study Group et al. (1987).  The gambia hepatitis intervention study.  Cancer Research 47(21), 5782–5787.
  * Garthwaite (1996) Garthwaite, P. H. (1996).  Confidence intervals from randomization tests.  Biometrics, 1387–1393.
  * Haines et al. (2017) Haines, T. P., K.-A. Bowles, D. Mitchell, L. O’Brien, D. Markham, S. Plumb, K. May, K. Philip, R. Haas, M. N. Sarkies, et al. (2017).  Impact of disinvestment from weekend allied health services across acute medical and surgical wards: 2 stepped-wedge cluster randomised controlled trials.  PLoS medicine 14(10), e1002412.
  * Heard and Rubin-Delanchy (2018) Heard, N. A. and P. Rubin-Delanchy (2018).  Choosing between methods of combining-values.  Biometrika 105(1), 239–246.
  * Hemming et al. (2018) Hemming, K., M. Taljaard, J. E. McKenzie, R. Hooper, A. Copas, J. A. Thompson, M. Dixon-Woods, A. Aldcroft, A. Doussau, M. Grayling, C. Kristunas, C. E. Goldstein, M. K. Campbell, A. Girling, S. Eldridge, M. J. Campbell, R. J. Lilford, C. Weijer, A. B. Forbes, and J. M. Grimshaw (2018).  Reporting of stepped wedge cluster randomised trials: Extension of the consort 2010 statement with explanation and elaboration.  BMJ 363, k1614.
  * Hennessy et al. (2016) Hennessy, J., T. Dasgupta, L. Miratrix, C. Pattanayak, and P. Sarkar (2016).  A conditional randomization test to account for covariate imbalance in randomized experiments.  Journal of Causal Inference 4(1), 61–80.
  * Hudgens and Halloran (2008) Hudgens, M. G. and M. E. Halloran (2008).  Toward causal inference with interference.  Journal of the American Statistical Association 103(482), 832–842.
  * Hughes et al. (2020) Hughes, J. P., P. J. Heagerty, F. Xia, and Y. Ren (2020).  Robust inference for the stepped wedge design.  Biometrics 76(1), 119–130.
  * Hussey and Hughes (2007) Hussey, M. A. and J. P. Hughes (2007).  Design and analysis of stepped wedge cluster randomized trials.  Contemporary Clinical Trials 28(2), 182–191.
  * Imbens and Rubin (2015) Imbens, G. W. and D. B. Rubin (2015).  Causal inference in statistics, social, and biomedical sciences.  Cambridge University Press.
  * Ji et al. (2017) Ji, X., G. Fink, P. J. Robyn, D. S. Small, et al. (2017).  Randomization inference for stepped-wedge cluster-randomized trials: an application to community-based health insurance.  The Annals of Applied Statistics 11(1), 1–20.
  * Karmakar and Small (2020) Karmakar, B. and D. S. Small (2020).  Assessment of the extent of corroboration of an elaborate theory of a causal hypothesis using partial conjunctions of evidence factors.  The Annals of Statistics 48(6), nil.
  * Kempthorne (1955) Kempthorne, O. (1955).  The randomization theory of experimental inference.  Journal of the American Statistical Association 50(271), 946–967.
  * Kenny et al. (2021) Kenny, A., E. Voldal, F. Xia, P. J. Heagerty, and J. P. Hughes (2021).  Analysis of stepped wedge cluster randomized trials in the presence of a time-varying treatment effect.  arXiv preprint arXiv:2111.07190.
  * Lehmann and Romano (2006) Lehmann, E. L. and J. P. Romano (2006).  Testing Statistical Hypotheses.  Springer Science & Business Media.
  * Li et al. (2021) Li, F., J. P. Hughes, K. Hemming, M. Taljaard, E. R. Melnick, and P. J. Heagerty (2021).  Mixed-effects models for the design and analysis of stepped wedge cluster randomized trials: An overview.  Statistical Methods in Medical Research 30(2), 612–639.
  * Li and Ding (2017) Li, X. and P. Ding (2017).  General forms of finite population central limit theorems with applications to causal inference.  Journal of the American Statistical Association 112(520), 1759–1769.
  * Marcus et al. (1976) Marcus, R., E. Peritz, and K. R. Gabriel (1976).  On closed testing procedures with special reference to ordered analysis of variance.  Biometrika 63(3), 655–660.
  * Middleton and Aronow (2015) Middleton, J. A. and P. M. Aronow (2015).  Unbiased estimation of the average treatment effect in cluster-randomized experiments.  Statistics, Politics and Policy 6(1-2), 39–75.
  * Muntinga et al. (2012) Muntinga, M. E., E. O. Hoogendijk, K. M. Van Leeuwen, H. P. Van Hout, J. W. Twisk, H. E. Van Der Horst, G. Nijpels, and A. P. Jansen (2012).  Implementing the chronic care model for frail older adults in the netherlands: study protocol of act (frail older adults: care in transition).  BMC geriatrics 12(1), 1–10.
  * Pitman (1937) Pitman, E. J. G. (1937).  Significance tests which may be applied to samples from any populations. ii. the correlation coefficient test.  Supplement to the Journal of the Royal Statistical Society 4(2), 225–232.
  * Puelz et al. (2019) Puelz, D., G. Basse, A. Feller, and P. Toulis (2019).  A graph-theoretic approach to randomization tests of causal effects under general interference.  Journal of the Royal Statistical Society: Series B (Statistical Methodology).
  * Rosenbaum (2002) Rosenbaum, P. R. (2002).  Observational studies.  Springer.
  * Rosenbaum (2010) Rosenbaum, P. R. (2010).  Evidence factors in observational studies.  Biometrika 97(2), 333–345.
  * Rosenbaum (2017) Rosenbaum, P. R. (2017).  The general structure of evidence factors in observational studies.  Statistical Science 32(4), 514–530.
  * Rubin (1980) Rubin, D. B. (1980).  Comment on "randomization analysis of experimental data: The fisher randomization test".  Journal of the American Statistical Association 75(371), 591.
  * Stouffer et al. (1949) Stouffer, S. A., E. A. Suchman, L. C. DeVinney, S. A. Star, and R. M. Williams Jr (1949).  The american soldier: Adjustment during army life.(studies in social psychology in world war ii), vol. 1.
  * Thompson et al. (2018) Thompson, J., C. Davey, K. Fielding, J. Hargreaves, and R. Hayes (2018).  Robust analysis of stepped wedge trials using cluster-level summaries within periods.  Statistics in medicine 37(16), 2487–2500.
  * Thompson et al. (2017) Thompson, J. A., K. L. Fielding, C. Davey, A. M. Aiken, J. R. Hargreaves, and R. J. Hayes (2017).  Bias and inference from misspecified mixed-effect models in stepped wedge trial analysis.  Statistics in Medicine 36(23), 3670–3682.
  * Twisk (2021) Twisk, J. W. (2021).  Analysis of Data from Randomized Controlled Trials: A Practical Guide.  Springer Nature.
  * Wang and De Gruttola (2017) Wang, R. and V. De Gruttola (2017).  The use of permutation tests for the analysis of parallel and stepped-wedge cluster-randomized trials.  Statistics in medicine 36(18), 2831–2843.
  * Welch (1937) Welch, B. L. (1937).  On the z-test in randomized blocks and latin squares.  Biometrika 29(1/2), 21–52.
  * Zhang and Zhao (2022) Zhang, Y. and Q. Zhao (2022).  What is a randomization test?  arXiv preprint arXiv:2203.10980.
  * Zhao and Ding (2021) Zhao, A. and P. Ding (2021).  Covariate-adjusted fisher randomization tests for the average treatment effect.  Journal of Econometrics.
  * Zheng and Zelen (2008) Zheng, L. and M. Zelen (2008).  Multi-center clinical trials: Randomization and ancillary statistics.  The Annals of Applied Statistics 2(2), 582.

##  Appendix A Technical Proofs

###  A.1 Proof of Lemma 1

###### Proof.

Consider $\mathcal{S}\in\mathcal{R}^{(j)}$ and
$\mathcal{S}^{\prime}\in\mathcal{R}^{(k)}$ for some $j,k\in[K]$ and
$\mathcal{S}\cap\mathcal{S}^{\prime}\neq\emptyset$. By the definition of
refinement and (5),
$\mathcal{S}\cap\mathcal{S}^{\prime}\in\underline{R}^{\\{j,k\\}}\subseteq\mathcal{R}^{(j)}\cup\mathcal{R}^{(k)}$,
so there exists integers $m$ such that
$\mathcal{S}\cap\mathcal{S}^{\prime}=\mathcal{S}^{(j)}_{m}$ or
$\mathcal{S}^{(k)}_{m}$. Because $\mathcal{R}^{(j)}$ and $\mathcal{R}^{(k)}$
are partitions, this means that $\mathcal{S}=\mathcal{S}^{(j)}_{m}$ or
$\mathcal{S}^{\prime}=\mathcal{S}^{(k)}_{m}$. In either case,
$\mathcal{S}\cap\mathcal{S}^{\prime}=\mathcal{S}$ or $\mathcal{S}^{\prime}$.

Now consider any $\bm{z}\in\mathcal{Z}$ and $j,k\in[K]$. Because
$\mathcal{S}_{\bm{z}}^{(j)},\mathcal{S}_{\bm{z}}^{(k)}\in\mathcal{R}^{[K]}$
and $\mathcal{S}_{\bm{z}}^{(j)}\cap\mathcal{S}_{\bm{z}}^{(k)}\neq\emptyset$
(because they contain at least $\bm{z}$), either
$\mathcal{S}_{\bm{z}}^{(j)}\supseteq\mathcal{S}_{\bm{z}}^{(k)}\text{or
}\mathcal{S}_{\bm{z}}^{(k)}\supseteq\mathcal{S}_{\bm{z}}^{(j)}$ must be true.
This means that we can order the conditioning events
$\mathcal{S}_{\bm{z}}^{(k)},k\in[K]$ according to the relation $\supseteq$.
Without loss of generality, we assume that, at $\bm{z}$,

$\mathcal{S}_{\bm{z}}^{(1)}\supseteq\mathcal{S}_{\bm{z}}^{(2)}\supseteq\cdots\supseteq\mathcal{S}_{\bm{z}}^{(K-1)}\supseteq\mathcal{S}_{\bm{z}}^{(K)}.$

Then $\bigcap_{k=1}^{K}\mathcal{S}_{\bm{z}}^{(k)}=\mathcal{S}_{\bm{z}}^{(K)}$
and $\bigcup_{k=1}^{K}\mathcal{S}_{\bm{z}}^{(k)}=\mathcal{S}_{\bm{z}}^{(1)}$.
Thus the intersection and the union of
$\\{\mathcal{S}_{\bm{z}}^{(k)}\\}_{k=1}^{K}$ are contained in
$\mathcal{R}^{[K]}$ which collects all $\mathcal{S}_{\bm{z}}^{(k)}$ by the
definition. As this is true for all $\bm{z}\in\mathcal{Z}$, we have
$\underline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$ and
$\overline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$. ∎

###  A.2 Proof of Lemma 2

###### Proof.

(i) Suppose $\mathcal{S}^{\prime}$, $\mathcal{S}^{\prime\prime}$ are two
distinct nodes in $\text{ch}(\mathcal{S})$. By Lemma 1, they are either
disjoint or nested. If they are nested, without loss of generality, suppose
$\mathcal{S}^{\prime\prime}\supset\mathcal{S}^{\prime}$. However, this
contradicts with the definition of the edge
$\mathcal{S}\rightarrow\mathcal{S}^{\prime}$, as by Definition 3 there should
be no $\mathcal{S}^{\prime\prime}$ satisfying
$\mathcal{S}\supset\mathcal{S}^{\prime\prime}\supset\mathcal{S}^{\prime}$.
This shows that any two nodes in $\text{ch}(\mathcal{S})$ are disjoint.

Next we show that the union of the sets in $\text{ch}(\mathcal{S})$ is
$\mathcal{S}$. Suppose there exists $\bm{z}\in\mathcal{S}$ such that
$\bm{z}\not\in\mathcal{S}^{\prime}$ for all
$\mathcal{S}^{\prime}\in\text{ch}(\mathcal{S})$. In consequence,
$\bm{z}\not\in\mathcal{S}^{\prime}$ for all
$\mathcal{S}^{\prime}\in\text{de}(\mathcal{S})$. Similar to the proof of Lemma
1, we can order $\mathcal{S}_{\bm{z}}^{(k)},k\in[K]$ according to set
inclusion. Without loss of generality, suppose

$\mathcal{S}_{\bm{z}}^{(1)}\supseteq\mathcal{S}_{\bm{z}}^{(2)}\supseteq\cdots\supseteq\mathcal{S}_{\bm{z}}^{(K-1)}\supseteq\mathcal{S}_{\bm{z}}^{(K)}.$

This shows that $\mathcal{S}=\mathcal{S}_{\bm{z}}^{(K)}$, so
$\text{ch}(\mathcal{S})=\text{de}(\mathcal{S})=\varnothing$. This contradicts
the assumption.

(ii) Consider any $\mathcal{S}\in\mathcal{R}^{[K]},$
$\mathcal{S}^{\prime}\in\text{an}(\mathcal{S})$ and
$\mathcal{S}^{\prime\prime}\in\text{de}(\mathcal{S})$. By definition,
$\mathcal{S}^{\prime}\supset\mathcal{S}\supset\mathcal{S}^{\prime\prime}$.
Because the sets in any partition $\mathcal{R}^{(k)}$ are disjoint, this shows
that no pairs of $\mathcal{S},\mathcal{S}^{\prime},\mathcal{S}^{\prime\prime}$
can belong to the same partition $\mathcal{R}^{(k)}$. Thus,
$\mathcal{K}(\mathcal{S}^{\prime}),\mathcal{K}(\mathcal{S}),\mathcal{K}(\mathcal{S}^{\prime\prime})$
are disjoint. Because this is true for any
$\mathcal{S}^{\prime}\in\text{an}(\mathcal{S})$ and
$\mathcal{S}^{\prime\prime}\in\text{de}(\mathcal{S})$, this shows
$\mathcal{K}(\text{an}(\mathcal{S}))$, $\mathcal{K}(\mathcal{S})$, and
$\mathcal{K}(\text{de}(\mathcal{S}))$ are disjoint.

Now consider any $\bm{z}\in\mathcal{S}$. By the proof of (i), $\mathcal{S}$ is
in a nested sequence consisting of
$\mathcal{S}_{\bm{z}}^{(1)},\dotsc,\mathcal{S}_{\bm{z}}^{(K)}$. Hence,
$\mathcal{K}(\text{an}(\mathcal{S}))\cup\mathcal{K}(\mathcal{S})\cup\mathcal{K}(\text{de}(\mathcal{S}))=\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\}\cup\text{de}(\mathcal{S}))=[K]$.

(iii) The first result follows immediately from the fact that
$\text{an}(\mathcal{S}^{\prime})=\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\}$.
The second result is true because both
$\mathcal{K}(\\{\mathcal{S}^{\prime}\\}\cup\text{de}(\mathcal{S}^{\prime}))$
and $\mathcal{K}(\text{de}(\mathcal{S}))$ are equal to
$[K]\setminus\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\})$. ∎

###  A.3 Proof of Lemma 3

###### Proof.

First, we claim that for any $j,k\in[K]$ and $\bm{z}\in\mathcal{Z}$, given
that $\bm{Z}\in\mathcal{S}^{(j)}_{\bm{z}}\cap\mathcal{S}^{(k)}_{\bm{z}}$,
$P^{(k)}(\bm{Z},\bm{W})$ only depends on $\bm{Z}$ through
$T_{\bm{z}}^{(k)}(\bm{Z},\bm{W})$. This is true because of the definition of
the p-value (1) and the invariance of the conditioning sets and test
statistics (so $\mathcal{S}_{\bm{Z}}^{(k)}=\mathcal{S}^{(k)}_{\bm{z}}$ and
$T_{\bm{Z}}(\cdot,\cdot)=T_{\bm{z}}(\cdot,\cdot)$). By using this claim, the
conditional independence (6) implies that

$P^{(j)}(\bm{Z},\bm{W})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}P^{(k)}(\bm{Z},\bm{W})\mid\mathcal{S}_{\bm{z}}^{(j)}\cap\mathcal{S}_{\bm{z}}^{(k)},\bm{W},\quad\forall
j,k\in[K],\bm{z}\in\mathcal{Z}.$

Now fix $\mathcal{S}\in\mathcal{R}^{[K]}$ and consider any
$j\in\mathcal{K}(\mathcal{S})$ and
$k\in\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\})\setminus\\{j\\}$.
By Lemma 1,
$\mathcal{S}_{\bm{z}}^{(j)}\cap\mathcal{S}_{\bm{z}}^{(k)}=\mathcal{S}_{\bm{z}}^{(j)}=\mathcal{S}$
for any $\bm{z}\in\mathcal{S}$. Then the claimed independence follows
immediately. ∎

###  A.4 Proof of Lemma 4

###### Proof.

Let
$\psi^{(k)}=\psi^{(k)}(\bm{Z},\bm{W})=1_{\\{P^{(k)}(\bm{Z},\bm{W})\leq\alpha^{(k)}\\}}$
for $k\in[K]$, so the left-hand side of (11) can be written as

$\mathbb{P}\left\\{P^{(1)}(\bm{Z},\bm{W})\leq\alpha^{(1)},\dotsc,P^{(K)}(\bm{Z},\bm{W})\leq\alpha^{(K)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}=\mathbb{E}\Bigg{\\{}\prod_{k=1}^{K}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\Bigg{\\}}.$

We prove Lemma 4 by induction. First, consider any leaf node in the Hasse
diagram, that is, any $\mathcal{S}\in\mathcal{R}^{[K]}$ such that
$\text{ch}(\mathcal{S})=\varnothing$. By Lemma 2,
$\mathcal{K}(\mathcal{S})\cup\mathcal{K}(\text{an}(\mathcal{S}))=[K]$. Then by
Lemma 3, $\psi^{(j)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\psi^{(k)}$ for any
$j\in\mathcal{K}(\mathcal{S})$ and $k\neq j$. By using the validity of each
CRT (Theorem 1), we obtain

$\displaystyle\mathbb{E}\Bigg{\\{}\prod_{k=1}^{K}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\Bigg{\\}}$
$\displaystyle=\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S}))}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\mathcal{S})}\mathbb{E}\left\\{\psi^{(j)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}$
$\displaystyle\leq\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S}))}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\mathcal{S})}\alpha^{(j)}.$

This is exactly (11) for a leaf node. Now consider a non-leaf node
$\mathcal{S}\in\mathcal{R}^{[K]}$ (so $\text{ch}(\mathcal{S})\neq\emptyset$)
and suppose (11) holds for any descendant of $\mathcal{S}$. We have

$\displaystyle\
\mathbb{E}\Bigg{\\{}\prod_{k=1}^{K}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\Bigg{\\}}$
$\displaystyle=$ $\displaystyle\
\mathbb{E}\left\\{\sum_{\mathcal{S}^{\prime}\in\text{ch}(\mathcal{S})}1_{\\{\bm{Z}\in\mathcal{S}^{\prime}\\}}\mathbb{E}\left\\{\prod_{k=1}^{K}\psi^{(k)}\mid\bm{Z}\in\mathcal{S}^{\prime},\bm{W}\right\\}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\qquad\quad\text{(By
\lx@cref{creftypecap~refnum}{lem:hasse-structure}\ref{equ:hasse1})}$
$\displaystyle\leq$ $\displaystyle\
\mathbb{E}\left\\{\sum_{\mathcal{S}^{\prime}\in\text{ch}(\mathcal{S})}1_{\\{\bm{Z}\in\mathcal{S}^{\prime}\\}}\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S}^{\prime}))}\psi^{(k)}\mid\bm{Z}\in\mathcal{S}^{\prime},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\\{\mathcal{S}^{\prime}\\}\cup\text{de}(\mathcal{S}^{\prime}))}\alpha^{(j)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}$
(By the induction hypothesis) $\displaystyle=$ $\displaystyle\
\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S})\cup\\{\mathcal{S}\\})}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\text{de}(\mathcal{S}))}\alpha^{(j)}\qquad\qquad\qquad\quad\text{(By
\lx@cref{creftypecap~refnum}{lem:hasse-structure}\ref{equ:hasse3})}$
$\displaystyle=$ $\displaystyle\
\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S}))}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\mathcal{S})}\mathbb{E}\left\\{\psi^{(j)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\text{de}(\mathcal{S}))}\alpha^{(j)}$
(By Lemma 2(ii) and Lemma 3) $\displaystyle\leq$ $\displaystyle\
\mathbb{E}\left\\{\prod_{k\in\mathcal{K}(\text{an}(\mathcal{S}))}\psi^{(k)}\mid\bm{Z}\in\mathcal{S},\bm{W}\right\\}\prod_{j\in\mathcal{K}(\\{\mathcal{S}\\}\cup\text{de}(\mathcal{S}))}\alpha^{(j)}.\qquad\qquad\qquad\qquad\text{(By
\lx@cref{creftypecap~refnum}{thm:valid})}$

By induction, this shows that (11) holds for all
$\mathcal{S}\in\mathcal{R}^{[K]}$. ∎

###  A.5 Proof of Theorem 2

###### Proof.

Lemma 1 shows that $\overline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$,
thus (11) holds for every
$\mathcal{S}\in\overline{\mathcal{R}}^{[K]}\subseteq\mathcal{R}^{[K]}$.
Moreover, any $\mathcal{S}\in\overline{\mathcal{R}}^{[K]}$ has no superset in
$\mathcal{R}^{[K]}$ and thus has no ancestors in the Hasse diagram. This means
that the right-hand side of (11) is simply $\prod_{k=1}^{K}\alpha^{(k)}$.
Because $\overline{\mathcal{G}}^{[K]}$ is the $\sigma$-algebra generated by
the events $\\{\bm{Z}\in\mathcal{S}\\}$ for
$\mathcal{S}\in\overline{\mathcal{R}}^{[K]}$, equation (7) holds. Finally,
equation (8) holds trivially by taking the expectation of (7) over
$\overline{\mathcal{G}}^{[K]}$. ∎

###  A.6 Proof of Proposition 3

We denote the units used in the $k$-th test by
$\mathcal{I}^{(k)}=\mathcal{I}_{0}^{(k)}\cup\mathcal{I}_{1}^{(k)}$. We denote
the treatment variables and outcomes used in the $k$-th test by
$\bm{A}^{(k)}=(A_{i}:i\in\mathcal{I}^{(k)})$ and
$\bm{Y}^{(k)}=(Y_{i}^{(k)}:i\in\mathcal{I}^{(k)})$, respectively. We define
the difference-in-means statistics in the $k$-th test as

$T(\bm{A}^{(k)},\bm{Y}^{(k)})=\sqrt{N}(\bar{Y}_{1}^{(k)}-\bar{Y}_{0}^{(k)}),$
(20)

where the average outcomes are given by

$\bar{Y}_{1}^{(k)}=(N_{1}^{(k)})^{-1}\sum_{i\in\mathcal{I}^{(k)}}A_{i}^{(k)}Y_{ik}\quad\text{and}\quad\bar{Y}_{0}^{(k)}=(N_{0}^{(k)})^{-1}\sum_{i\in\mathcal{I}^{(k)}}(1-A_{i}^{(k)})Y_{ik}.$

Let
$\bm{W}^{(k)}=\big{(}Y_{i}^{(k)}(0),Y_{i}^{(k)}(1):i\in\mathcal{I}^{(k)}\big{)}$.
The randomization distribution in the $k$-th test is

$\hat{G}^{(k)}(b^{(k)})=\mathbb{P}^{*}\left\\{T\big{(}\bm{A}_{*}^{(k)},\bm{Y}^{(k)}(\bm{A}_{*}^{(k)})\big{)}\leq
b^{(k)}\mid\bm{W}^{(k)}\right\\}$

where $\bm{A}_{*}^{(k)}$ is a permutation of $\bm{A}^{(k)}$, i.e.,
$\bm{A}_{*}^{(k)}$ is drawn from the same uniform assignment distribution as
$\bm{A}^{(k)}$. Under assumption (i), using the bivariate Central Limit
Theorem,

$\Big{(}T\big{(}\bm{A}^{(k)},\bm{Y}^{(k)}\big{)},T\big{(}\bm{A}_{*}^{(k)},\bm{Y}^{(k)}(\bm{A}_{*}^{(k)})\big{)}\Big{)}\xrightarrow{d}\big{(}T_{\infty}^{(k)},T_{\infty,*}^{(k)}\big{)},$

where $T_{\infty}^{(k)}$ and $T_{\infty,*}^{(k)}$ are independent, each with a
common c.d.f $G^{(k)}(\cdot).$ By Lehmann and Romano (2006, Theorem 15.2.3),

$\hat{G}^{(k)}(b^{(k)})\xrightarrow{P}G^{(k)}(b^{(k)})$

for every $b^{(k)}$ which is a continuity point of $G^{(k)}(\cdot)$.

Under assumptions (i), (ii) and the zero-effect null hypothesis $H_{0}^{(k)}$,
following Lehmann and Romano (2006, Theorem 15.2.5), we have

$T_{\infty}^{(k)}\sim g_{0}^{(k)}=\mathcal{N}(0,V_{\infty}^{(k)}),$ (21)

where the variance $V_{\infty}^{(k)}$ takes the form

$V_{\infty}^{(k)}=1/\Lambda^{(k)}=\frac{N}{N_{0}^{(k)}}\operatorname{Var}\big{[}Y^{(k)}(1)\big{]}+\frac{N}{N_{1}^{(k)}}\operatorname{Var}\big{[}Y^{(k)}(0)\big{]}.$

This expression of the asymptotic variance can be found in Imbens and Rubin
(2015, Section 6.4). For completeness, an alternative derivation is provided
in Section A.6.1. Under assumptions (i), (ii) and the constant effect
alternative hypothesis $H_{1}^{(k)}$ (with $\tau=h/\sqrt{N}$), we have

$T_{\infty}^{(k)}\sim g_{1}^{(k)}=\mathcal{N}(h,V_{\infty}^{(k)}).$ (22)

The c.d.f of $T_{\infty}^{(k)}$ under $H_{0}^{(k)}$ and $H_{1}^{(k)}$ are
given by

$G_{0}^{(k)}(b^{(k)})=\Phi\left(b^{(k)}/\sqrt{V_{\infty}^{(k)}}\right)\quad\text{and}\quad
G_{1}^{(k)}(b^{(k)})=\Phi\left(\Big{[}b^{(k)}-h\Big{]}/\sqrt{V_{\infty}^{(k)}}\right).$

Let $\tilde{b}^{(k)}=\big{[}G_{0}^{(k)}\big{]}^{-1}(p^{(k)})$. The p-value
density function under $H_{1}^{(k)}$ can be rewritten as

$f_{1}^{(k)}(p^{(k)})=g_{1}^{(k)}\big{(}\tilde{b}^{(k)}\big{)}\left|\frac{d\tilde{b}^{(k)}}{dp^{(k)}}\right|=g_{1}^{(k)}(\tilde{b}^{(k)})\left|\bigg{(}\big{[}G_{0}^{(k)}\big{]}^{\prime}(\tilde{b}^{(k)})\bigg{)}^{-1}\right|=g_{1}^{(k)}(\tilde{b}^{(k)})/g_{0}^{(k)}(\tilde{b}^{(k)}).$

Since the p-values from the permutation tests follow a standard uniform
distribution under the null, the log-likelihood ratio of the p-values
$P^{(1)}\dotsc,P^{(K)}$ is given by

$\sum_{k=1}^{K}\log\left[f_{1}^{(k)}\big{(}p^{(k)}\big{)}\right]=\sum_{k=1}^{K}\log\left[\frac{g_{1}^{(k)}(\tilde{b}^{(k)})}{g_{0}^{(k)}(\tilde{b}^{(k)})}\right]\propto\sum_{k=1}^{K}\sqrt{\Lambda^{(k)}}\Phi^{-1}(p^{(k)}).$

Using the p-value weights $\Lambda^{(k)},k\in[K],$ from the log-likelihood
ratio, we obtain

$T_{\infty}=\sum_{k=1}^{K}w_{\infty}^{(k)}\Phi^{-1}(P^{(k)})\quad\text{where}\quad
w_{\infty}^{(k)}=\sqrt{\frac{\Lambda^{(k)}}{\sum_{j=1}^{K}\Lambda^{(j)}}}.$
(23)

####  A.6.1 Additional proof for Equation 21

Lehmann and Romano (2006, Theorem 15.2.3) uses the fact that under assumptions
(i) and (ii),
$V_{\infty}^{(k)}=\operatorname{Var}[T\big{(}\bm{A}^{(k)},\bm{Y}^{(k)}\big{)}]$,
but leaves the calculation of
$\operatorname{Var}[T\big{(}\bm{A}^{(k)},\bm{Y}^{(k)}\big{)}]$ (their Equation
15.15) as an exercise for readers. Here we provide the calculation to complete
our proof for (21). To simplify the exposition, we let $m=N_{1}^{(k)}$,
$n=N_{0}^{(k)}$. Suppose that $\mathcal{I}^{(k)}=[m+n]$,
$\mathcal{I}_{1}^{(k)}=[m]$ and $\mathcal{I}_{0}^{(k)}=\\{m+1,\dotsc,m+n\\},$
and that

$\bm{Y}^{(k)}=(\underbrace{Y_{1}^{(k)},\cdots,Y_{m}^{(k)}}_{\text{treated
outcomes}},\underbrace{Y_{m+1}^{(k)},\dotsc,Y_{m+n}^{(k)}}_{\text{control
outcomes}}).$

Let $\Pi=\big{[}\Pi(1),\dotsc,\Pi(m+n)\big{]}$ be an independent random
permutation of $1,\dotsc,m+n.$ We rewrite the difference-in-means statistics
(20) as

$T:=T\big{(}\bm{A}^{(k)},\bm{Y}^{(k)}\big{)}=\frac{\sqrt{N}}{m}\sum_{i=1}^{m+n}E_{i}Y_{i}^{(k)},\text{
where }E_{i}=\begin{cases}1&\text{if }\Pi(i)\leq m\\\
-m/n&\text{otherwise}\end{cases},\forall i\in[m+n].$

Let $D$ be the number of $i\leq m$ such that $\Pi(i)\leq m$. The value of $T$
is determined by $m-D$, the number of units swapped between the treated group
($i=1,\dotsc m$) and control group ($i=m+1,\dotsc,m+n$). Since all the treated
(control) units have the same outcome variance, it does not matter which units
are swapped in computing the variance of $T$.

We first consider the case that $m\leq n$. Using the expectation of the
hypergeometric distribution, we know

$\mathbb{E}[D]=\sum_{d=0}^{m}\frac{{m\choose d}{n\choose m-d}}{{m+n\choose
m}}d=\frac{m^{2}}{m+n}.$

Let $\sigma_{1}^{2}=\operatorname{Var}[Y^{(k)}(1)]$ and
$\sigma_{0}^{2}=\operatorname{Var}[Y^{(k)}(0)].$ The variance of $T$ is given
by

$\displaystyle\operatorname{Var}(T)$
$\displaystyle=\frac{N}{m^{2}}\sum_{d=0}^{m}\frac{{m\choose d}{n\choose
m-d}}{{m+n\choose
m}}\left[d\sigma_{1}^{2}+(m-d)\sigma_{0}^{2}+(m-d)\frac{m^{2}}{n^{2}}\sigma_{1}^{2}+(n-m+d)\frac{m^{2}}{n^{2}}\sigma_{0}^{2}\right]$
$\displaystyle=\frac{N}{m^{2}}\left[\frac{m^{2}}{m+n}\sigma_{1}^{2}+\frac{mn}{m+n}\sigma_{0}^{2}+\frac{mn}{m+n}\frac{m^{2}}{n^{2}}\sigma_{1}^{2}+\frac{m^{2}}{m+n}\sigma_{0}^{2}\right]$
$\displaystyle=\frac{N}{m^{2}}\left[\frac{m^{2}}{n}\sigma_{1}^{2}+m\sigma_{0}^{2}\right]$
$\displaystyle=\frac{N}{n}\sigma_{1}^{2}+\frac{N}{m}\sigma_{0}^{2}.$

If $m>n$, we let $D$ denote the number of $i\in\\{m+1,\dotsc,m+n\\}$ such that
$\Pi(i)\in\\{m+1,\dotsc,m+n\\}$. Then, $\mathbb{E}[D]=\frac{n^{2}}{m+n}$,  and

$\displaystyle\operatorname{Var}(T)$
$\displaystyle=\frac{N}{m^{2}}\sum_{d=0}^{n}\frac{{m\choose d}{n\choose
m-d}}{{m+n\choose
n}}\left[d\frac{m^{2}}{n^{2}}\sigma_{0}^{2}+(n-d)\frac{m^{2}}{n^{2}}\sigma_{1}^{2}+(n-d)\sigma_{0}^{2}+(m-n+d)\sigma_{1}^{2}\right]$
$\displaystyle=\frac{N}{m^{2}}\left[\frac{m^{2}}{m+n}\sigma_{0}^{2}+\frac{mn}{m+n}\frac{m^{2}}{n^{2}}\sigma_{1}^{2}+\frac{mn}{m+n}\sigma_{0}^{2}+\frac{m^{2}}{m+n}\sigma_{1}^{2}\right]$
$\displaystyle=\frac{N}{m^{2}}\left[\frac{m^{2}}{n}\sigma_{1}^{2}+m\sigma_{0}^{2}\right]$
$\displaystyle=\frac{N}{n}\sigma_{1}^{2}+\frac{N}{m}\sigma_{0}^{2}.$

The variance is the same for $m\leq n$ and $m>n$. This proves our claim and
Equation 15.15 in Lehmann and Romano (2006):
$\operatorname{Var}(m^{-1/2}\sum_{i=1}^{m+n}E_{i}Y_{i}^{(k)})=\operatorname{Var}(\sqrt{\frac{m}{N}}T)=\frac{m}{n}\sigma_{1}^{2}+\sigma_{0}^{2}.$

###  A.7 Proof of Proposition 4

In Algorithm 1, the time steps in $\mathcal{C}_{j}$ define a sequence of
nested permutation tests. For example, $\mathcal{C}_{1}=\\{1,3,5,7\\}$ defines
CRTs 1,3,5 and $\mathcal{C}_{2}=\\{2,4,6,8\\}$ defines CRTs 2,4,6 in Figure 2.
Suppose that we only have one subset $\mathcal{C}=\\{c_{1},\dotsc,c_{K}\\}$
consists of $K$ time points. We next prove the tests defined on $\mathcal{C}$
are valid to combine. It is suffices to show that the test statistics
$\hat{T}=\sum_{k=1}^{K}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})$ stochastically
dominates the random variable

$\tilde{T}:=\sum_{k=1}^{K}\hat{w}^{(k)}\Phi^{-1}(U^{(k)})\sim\mathcal{N}\left(0,\sum_{k=1}^{K}\big{[}\hat{w}^{(k)}\big{]}^{2}\right)=\mathcal{N}(0,1).$

where each $U^{(k)}$ is a standard uniform random variable. The second
equality is achieved by the definition of $\hat{w}^{(k)},k\in[K].$

By conditioning on the potential outcomes $(\bm{W}^{(k)},k\in[K])$, the
weights $\hat{w}^{(k)},k\in[K],$ are fixed. The derivation follows the same
steps as Proposition 2 (the conditioning on $\bm{W}$’s is suppressed). By
construction, $c_{1}<\dotsc<c_{K}$ and the conditioning sets
$\mathcal{S}^{(1)}\supseteq\cdots\supseteq\mathcal{S}^{(K)}$ for all
$\bm{z}\in\mathcal{Z}$. This implies that
$\mathcal{G}^{(1)}\subseteq\cdots\subseteq\mathcal{G}^{(K)}$. Conditioning on
$\mathcal{G}^{(k^{\prime})}$, the term
$\sum_{k=1}^{k^{\prime}-1}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})$ is fixed. By the
law of iterated expectation,

$\displaystyle\mathbb{E}\left[1\left\\{\hat{T}\leq b\right\\}\right]=\ $
$\displaystyle\mathbb{E}\left[1\left\\{\hat{w}^{(K)}\Phi^{-1}(P^{(K)})\leq
b-\sum_{k=1}^{K-1}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})\right\\}\right]$
$\displaystyle=\ $
$\displaystyle\mathbb{E}\left(\mathbb{E}\left[1{\left\\{\hat{w}^{(K)}\Phi^{-1}(P^{(K)})\leq
b-\sum_{k=1}^{K-1}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})\right\\}}\ \Big{|}\
\mathcal{G}^{(K)}\right]\right)$ $\displaystyle\leq\ $
$\displaystyle\mathbb{E}_{U^{(K)}}\left(\mathbb{E}\left[1\left\\{\hat{w}^{(K)}\Phi^{-1}(U^{(K)})\leq
b-\sum_{k=1}^{K-1}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})\right\\}\ \Big{|}\
U^{(K)}\right]\right)$ $\displaystyle=\ $
$\displaystyle\mathbb{E}_{U^{(K)}}\bigg{(}\mathbb{E}\bigg{[}1\bigg{\\{}\hat{w}^{(K-1)}\Phi^{-1}(P^{(K-1)})\leq
b-\sum_{k=1}^{K-2}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})$ $\displaystyle\hskip
150.0pt-\hat{w}^{(K)}\Phi^{-1}(U^{(K)})\bigg{\\}}\ \Big{|}\
\mathcal{G}^{(K-1)},U^{(K)}\bigg{]}\bigg{)}$ $\displaystyle\leq\ $
$\displaystyle\mathbb{E}_{U^{(K-1)},U^{(K)}}\bigg{(}\mathbb{E}\bigg{[}1\bigg{\\{}\hat{w}^{(K-1)}\Phi^{-1}(U^{(K-1)})\leq
b-\sum_{k=1}^{K-2}\hat{w}^{(k)}\Phi^{-1}(P^{(k)})$ $\displaystyle\hskip
150.0pt-\hat{w}^{(K)}\Phi^{-1}(U^{(K)})\bigg{\\}}\ \Big{|}\
U^{(K-1)},U^{(K)}\bigg{]}\bigg{)}$ $\displaystyle\vdots$ $\displaystyle\leq\ $
$\displaystyle\mathbb{E}_{U^{(1)},\dotsc,U^{(K)}}\bigg{(}\mathbb{E}\bigg{[}1\bigg{\\{}\sum_{k=1}^{K}\hat{w}^{(k)}\Phi^{-1}(U^{(k)})\leq
b\bigg{\\}}\ \Big{|}\ U^{(1)},\dotsc,U^{(K)}\bigg{]}\bigg{)}$ $\displaystyle=\
$ $\displaystyle\mathbb{E}\left[1{\left\\{\tilde{T}\leq b\right\\}}\right].$

The inequalities are attained by the validity of the p-values
$P^{(1)},\dotsc,P^{(K)}$, i.e., each $P^{(k)}$ stochastically dominates the
standard uniform variable $U^{(k)}$. Since $\hat{T}$ stochastically dominates
the standard normal random variable $\tilde{T}$,

$\mathbb{P}\\{\hat{P}\leq\alpha\\}=\mathbb{P}\\{\Phi(\hat{T})\leq\alpha\\}\leq\mathbb{P}\\{\Phi(\tilde{T})\leq\alpha\\}=\alpha$
(24)

The last equality is achieved by the fact that $\Phi(\tilde{T})$ is a standard
uniform random variable.

Suppose that Algorithm 1 creates multiple subsets $\mathcal{C}_{j},j\in[J].$
Applying the same proof to the tests defined on each $\mathcal{C}_{j}$, we
will have

$\mathbb{E}\left[1\left\\{\hat{T}_{j}\leq
b\right\\}\right]\leq\mathbb{E}\left[1\left\\{\tilde{T}_{j}\leq
b\right\\}\right]$

where $\hat{T}_{j}=\sum_{c\in\mathcal{C}_{j}}\hat{w}^{(c)}\Phi(P^{(c)})$ and
$\tilde{T}_{j}=\sum_{c\in\mathcal{C}_{j}}\hat{w}^{(c)}\Phi(U^{(c)})\sim\mathcal{N}\left(0,\sum_{c\in\mathcal{C}_{j}}\big{[}\hat{w}^{(c)}\big{]}^{2}\right)$.
Because of sample splitting, we can combine the p-values from the tests based
on different $\mathcal{C}_{j}$. The test statistics
$\hat{T}=\sum_{j\in\mathcal{J}}\hat{T}_{j}$ stochastically dominates the
random variable

$\tilde{T}=\sum_{j\in\mathcal{J}}\tilde{T}_{j}\sim\mathcal{N}\bigg{(}0,\sum_{j\in\mathcal{J}}\sum_{c\in\mathcal{C}_{j}}\big{[}\hat{w}^{(c)}\big{]}^{2}\bigg{)}=\mathcal{N}(0,1),$

which implies that (24) still holds for the combined p-value
$\hat{P}=\Phi(\hat{T})$.

##  Appendix B Confidence intervals from randomization tests

Here we describe how to invert (a combination of) permutation tests and the
complexity involves; see also Ernst (2004, Section 3.4) for an introduction.
Consider a completely randomized experiment with treated outcomes
$(Y_{1},\dotsc,Y_{m})$ and control outcomes $(Y_{m+1},\dotsc,Y_{m+n})$.
Consider the constant effect null hypothesis

$H_{0}:Y_{i}(1)=Y_{i}(0)+\Delta,\forall i\in[m+n],$

We can implement a permutation test for $H_{0}$ by testing the zero-effect
hypothesis on the shifted outcomes
$\bm{Y}_{\Delta}=(Y_{1}-\Delta,\dotsc,Y_{m}-\Delta,Y_{m+1},\dotsc,Y_{m+n})$.

To construct a confidence interval for the true constant treatment effect
$\tau$, we can consider all values of $\Delta$ for which we do not reject
$H_{0}$. We define the left and right tails of the randomization distribution
of $T(\bm{Z}^{*},\bm{Y}_{\Delta})$ by

$P_{1}(\Delta)=\mathbb{P}^{*}\\{T(\bm{Z}^{*},\bm{Y}_{\Delta})\leq
T(\bm{Z},\bm{Y}_{\Delta})\\}\text{ and
}P_{2}(\Delta)=\mathbb{P}^{*}\\{T(\bm{Z}^{*},\bm{Y}_{\Delta})\geq
T(\bm{Z},\bm{Y}_{\Delta})\\},$

where $\bm{Z}=[\bm{1}_{m}/m,-\bm{1}_{n}/n]$ is observed assignment,
$\bm{Z}^{*}$ is a permutation of $\bm{Z}$, $T$ is the test statistics, e.g.,
difference-in-means, $T(\bm{Z},\bm{Y}_{\Delta})=\bm{Z}^{\top}Y_{\Delta}$. The
complexity of computing $P_{1}(\Delta)$ and $P_{2}(\Delta)$ with $b$ different
permutations is $O(mb+nb)$. The two-sided $(1-\alpha)$-confidence interval for
$\Delta$ is given by

$[\Delta_{L},\Delta_{U}]:=\left[\min_{P_{2}(\Delta)>\alpha/2}\Delta,\max_{P_{1}(\Delta)>\alpha/2}\Delta\right].$

We use a grid search to approximately find the minimum and maximum $\Delta$ in
$[\Delta_{L},\Delta_{U}]$. Since $P_{1}(\Delta)$ and $P_{2}(\Delta)$ are
monotone functions of $\Delta$, we can also consider obtaining the optimal
$\Delta$’s via a root-finding numerical method (see e.g. Garthwaite, 1996).

Inverting a combination of permutation tests can be done similarly. For
example, by searching the same $\Delta$’s for every permutation test, the
lower confidence bound $\Delta_{L}$ is given by the minimum $\Delta$ under
which the p-value of the combined test (e.g. $\hat{P}$ in Proposition 4) is
larger than $\alpha/2$. The complexity of inverting a combined test scales
linearly in terms of the number of tests. The total complexity is manageable
as long as the numbers of units and permutations are not very large at the
same time. If we have covariates information in the dataset, we can consider
fitting a linear regression model (with basis functions for nonlinearity). We
can compute the matrix inversion in the least-squares solution once, then
update the solution by shifting the outcome vector with different $\Delta$’s.
Inverting the test returns an interval with coverage probability approximately
equal to $1-\alpha$. When the probability is slightly below $1-\alpha$, we can
decrease $\alpha$ by a small value (e.g. 0.0025) gradually and reconstruct the
interval based on the previously searched $\Delta$’s. We stop if we obtain an
interval with coverage probability above $1-\alpha$.