# Design-based theory for cluster rerandomization

Xin Lu1, Tianle Liu2, Hanzhong Liu1, Peng Ding3

1 Center for Statistical Science, Department of Industrial Engineering,
Tsinghua University, Beijing, China
2 Department of Statistics, Harvard University, Cambridge, MA, U.S.A.
3 Department of Statistics, University of California, Berkeley, California,
U.S.A.

###### Abstract

Complete randomization balances covariates on average, but covariate imbalance
often exists in finite samples. Rerandomization can ensure covariate balance
in the realized experiment by discarding the undesired treatment assignments.
Many field experiments in public health and social sciences assign the
treatment at the cluster level due to logistical constraints or policy
considerations. Moreover, they are frequently combined with rerandomization in
the design stage. We refer to cluster rerandomization as a cluster-randomized
experiment compounded with rerandomization to balance covariates at the
individual or cluster level. Existing asymptotic theory can only deal with
rerandomization with treatments assigned at the individual level, leaving that
for cluster rerandomization an open problem. To fill the gap, we provide a
design-based theory for cluster rerandomization. Moreover, we compare two
cluster rerandomization schemes that use prior information on the importance
of the covariates: one based on the weighted Euclidean distance and the other
based on the Mahalanobis distance with tiers of covariates. We demonstrate
that the former dominates the latter with optimal weights and orthogonalized
covariates. Last but not least, we discuss the role of covariate adjustment in
the analysis stage and recommend covariate-adjusted procedures that can be
conveniently implemented by least squares with the associated robust standard
errors.

Key words: cluster randomization; cluster-robust standard error; constrained
randomization; covariate balance

##  1 Introduction

Cluster randomization has been widely used in public health (Donner and Klar,,
2000; Turner et al., 2017a, ; Turner et al., 2017b, ; Hayes and Moulton,,
2017) and social sciences (Raudenbush,, 1997; Schochet,, 2013; Athey and
Imbens,, 2017; Schochet,, 2020). It assigns the treatments at the cluster
level, with units within a cluster receiving the same treatment or control
condition, which helps to avoid interference within clusters and is applicable
when individual-level assignments are logistically infeasible.

Before experiments, researchers often collect covariates at the individual or
cluster level. For instance, in clinical trials, individual-level covariates
may include gender, race, age, and education of patients while cluster-level
covariates may include the capacity of the clinics and whether they are in a
metropolitan area (Li et al.,, 2016). The probability of having imbalance in
covariates after treatment assignments is non-negligible and would complicate
the interpretation of the experimental results, especially if some particular
covariates are predictive to the outcome.

A common approach to tackle covariate imbalance is to perform regression
adjustment in the analysis stage of an experiment (Lin,, 2013; Su and Ding,,
2021). In contrast, rerandomization actively reduces the imbalance in the
design stage. Rerandomization, also known as constrained randomization
(Moulton,, 2004; Morgan and Rubin,, 2012; Li et al.,, 2016, 2017, 2018), is a
procedure where one discards the undesired assignments that result in highly
unbalanced covariates and repeats randomization until a proper assignment
appears. The existing design-based theory for rerandomization assumes that the
treatments are assigned at the individual level (Morgan and Rubin,, 2012; Li
et al.,, 2018), leaving it unclear how rerandomization would theoretically
work in cluster experiments. The overarching goal of this paper is to fill
this gap.

We study cluster rerandomization with both individual- and cluster-level
covariates, and derive an asymptotic theory for estimators either with or
without regression adjustment. Our theory precisely quantifies the asymptotic
efficiency gain of performing cluster rerandomization at both levels. We
further show that under cluster rerandomization, the adjusted estimators based
on linear regressions recommended by Su and Ding, (2021) still have desired
theoretical guarantees. All these inferential procedures can be conveniently
realized by standard statistical packages. Moreover, previous asymptotic
theory mainly focuses on the rerandomization scheme based on the Mahalanobis
distance (Morgan and Rubin,, 2012; Li et al.,, 2018, 2020) while we also
discuss a distinct but widely-adopted practical scheme based on the weighted
Euclidean distance where one could utilize prior knowledge of the covariate
importance (Raab and Butcher,, 2001; Wight et al.,, 2002; Althabe et al.,,
2008; de Hoop et al.,, 2012; Li et al.,, 2016, 2017; Hayes and Moulton,, 2017;
Dempsey et al.,, 2018). We show that for orthogonalized covariates with
optimal weights, it dominates the cluster rerandomization with tiers of
covariates, another scheme accounting for the covariate importance, the non-
cluster version of which is proposed in Morgan and Rubin, (2015).

To facilitate the discussion, we introduce the following notation. Let
$\operatorname{det}(\cdot)$ denote the determinant of a matrix. For a real-
valued matrix $A$, we write $A>0$ if $A$ is strictly positive definite. For a
finite population $\\{a_{1},\ldots,a_{M}\\}$, let
$\operatorname{var}_{\textup{f}}(a)=(M-1)^{-1}\sum_{i=1}^{M}(a_{i}-\bar{a})^{2}$
denote its finite-population variance, where
$\bar{a}=M^{-1}\sum_{i=1}^{M}a_{i}$. Analogously, let
$\operatorname{cov}_{\textup{f}}$ denote the finite-population covariance. For
two random sequences $A$ and $B$, let $A\mathrel{\dot{\sim}}B$ indicate that
they have the same asymptotic distribution. Let
$\textnormal{pr}_{\textnormal{a}}$, $\textnormal{var}_{\textnormal{a}}$ and
$\textnormal{cov}_{\textnormal{a}}$ denote the asymptotic probability,
variance, and covariance, respectively. Besides, our asymptotic theory
requires mild moment conditions, but we relegate them to the appendix to avoid
clutter in the main text.

##  2 Notation and review for cluster-randomized experiments

Consider a cluster-randomized experiment with a finite population of $N$ units
grouped into $M$ clusters. The experimenter randomly assigns $M_{1}$ clusters
to the treatment arm and $M_{0}=M-M_{1}$ clusters to the control arm, with
units in the same cluster receiving the same treatment or control. For cluster
$i$, let $Z_{i}$ be the treatment indicator with $Z_{i}=1$ when it is assigned
to the treatment arm and $Z_{i}=0$ when it is assigned to the control arm. Let
$n_{i}$ be the size of cluster $i$, with $\sum_{i=1}^{M}n_{i}=N$. Let $(i,j)$
or simply $ij$ to index unit $j$ in cluster $i$. Let $Z_{ij}$ be the treatment
indicator for unit $(i,j)$ and by design, $Z_{ij}=Z_{i}$ for
$j=1,\ldots,n_{i}$. Let $N_{1}=\sum_{i=1}^{M}Z_{i}n_{i}$ and $N_{0}=N-N_{1}$
denote the numbers of units assigned to the treatment and control arms,
respectively. Importantly, $N_{1}$ and $N_{0}$ are random if the $n_{i}$’s
vary. Let $Y_{ij}(1)$ and $Y_{ij}(0)$ denote the potential outcomes of unit
$(i,j)$ under the treatment and control, respectively. We adopt the design-
based framework by conditioning on the potential outcomes with the $Z_{i}$’s
being the only source of randomness. The observed outcome
$Y_{ij}=Z_{i}Y_{ij}(1)+(1-Z_{i})Y_{ij}(0)$ is random in general due to the
random $Z_{i}$’s. We are interested in estimating the average treatment effect

$\tau=N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}\\{Y_{ij}(1)-Y_{ij}(0)\\}$

based on the observed data.

The difference-in-means estimator of $\tau$ is given by

$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}=N_{1}^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}Z_{ij}Y_{ij}-N_{0}^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}(1-Z_{ij})Y_{ij},$

with the subscript “haj” signifying that it is the Hajek estimator based on
the inverse of the treatment probabilities (Su and Ding,, 2021). As pointed
out by Middleton and Aronow, (2015), the random denominators in
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ can result in bias when the
$n_{i}$’s vary. The unbiased Horvitz–Thompson estimator replaces the random
denominators $N_{1}$ and $N_{0}$ with their expectations:

$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}=(NM_{1}/M)^{-1}\sum_{i=1}^{M}Z_{i}\sum_{j=1}^{n_{i}}Y_{ij}-(NM_{0}/M)^{-1}\sum_{i=1}^{M}(1-Z_{i})\sum_{j=1}^{n_{i}}Y_{ij}.$

Define the scaled cluster total potential outcome as
$\tilde{Y}_{i\cdot}(z)=\sum_{j=1}^{n_{i}}Y_{ij}(z)M/N$ so that
$\tau=M^{-1}\sum_{i=1}^{M}\\{\tilde{Y}_{i\cdot}(1)-\tilde{Y}_{i\cdot}(0)\\}.$
Based on the observed analog
$\tilde{Y}_{i\cdot}=Z_{i}\tilde{Y}_{i\cdot}(1)+(1-Z_{i})\tilde{Y}_{i\cdot}(0)$,
the Horvitz–Thompson estimator simplifies to

$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}=M_{1}^{-1}\sum_{i=1}^{M}Z_{i}\tilde{Y}_{i\cdot}-M_{0}^{-1}\sum_{i=1}^{M}(1-Z_{i})\tilde{Y}_{i\cdot}.$

This is unbiased for $\tau$ because it is the scaled difference in means of
the potential outcomes if each cluster is treated as a whole.

Both $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ are asymptotically normal
under some mild regularity conditions (Su and Ding,, 2021). Let
$e_{z}=M_{z}/M$ be the proportion of clusters assigned to the treatment arm
$z$ ($z=0,1$). Denote the normalized cluster size vector by
$\tilde{\omega}=(\tilde{\omega}_{1},\ldots,\tilde{\omega}_{M})^{{\mathrm{\scriptscriptstyle
T}}}$, with $\tilde{\omega}_{i}=n_{i}M/N$. Let
$\bar{Y}(z)=N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}Y_{ij}(z)$. Define
$\tilde{\varepsilon}(z)=(\tilde{\varepsilon}_{1\cdot}(z),\ldots,\tilde{\varepsilon}_{M\cdot}(z))^{{\mathrm{\scriptscriptstyle
T}}}$ and
$\tilde{Y}(z)=(\tilde{Y}_{1\cdot}(z),\ldots,\tilde{Y}_{M\cdot}(z))^{{\mathrm{\scriptscriptstyle
T}}}$, where

$\tilde{\varepsilon}_{i\cdot}(z)=\tilde{Y}_{i\cdot}(z)-\tilde{\omega}_{i}\bar{Y}(z)=\sum_{j=1}^{n_{i}}\varepsilon_{ij}(z)M/N,\quad\varepsilon_{ij}(z)=Y_{ij}(z)-\bar{Y}(z).$

Su and Ding, (2021) showed that $M^{1/2}(\hat{\tau}_{\star}-\tau)\mathrel{\
\dot{\sim}\ }\mathcal{N}(0,V_{\star,\tau\tau})$ for
$\star=\textnormal{ht},\textnormal{haj}$, where

$\displaystyle
V_{\textnormal{ht},\tau\tau}=e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}(1)\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}(0)\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}(1)-\tilde{Y}(0)\\},$
(1) $\displaystyle
V_{\textnormal{haj},\tau\tau}=e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{{\varepsilon}}(1)\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{{\varepsilon}}(0)\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{{\varepsilon}}(1)-\tilde{{\varepsilon}}(0)\\}.$
(2)

Therefore, $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ is unbiased in
finite samples and $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ is
unbiased asymptotically. By Middleton and Aronow, (2015) and Su and Ding,
(2021), $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ often practically
outperforms $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ in terms of the
mean squared error.

##  3 Cluster rerandomization using the Mahalanobis distance

###  3.1 Design

In a cluster-randomized experiment, suppose we observe baseline covariates at
the individual or cluster level. Let
$x_{ij}=(x_{ij1},\ldots,x_{ijK})^{{\mathrm{\scriptscriptstyle T}}}\\!\\!$
denote the individual-level covariates for unit $(i,j)$, and
$c_{i}=(c_{i1},\ldots,c_{iK})^{{\mathrm{\scriptscriptstyle T}}}$ be the
cluster-level covariates for cluster $i$. To simplify the presentation, we
center them to ensure that
$\bar{x}=N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}x_{ij}=0$ and
$M^{-1}\sum_{i=1}^{M}c_{i}=0$, and use the same $K$ to denote their
dimensions, although they may differ. Let
$\tilde{x}_{i\cdot}=\sum_{j=1}^{n_{i}}x_{ij}M/N$ be the scaled cluster totals,
concatenated as
$\tilde{X}=(\tilde{x}_{1\cdot},\ldots,\tilde{x}_{M\cdot})^{{\mathrm{\scriptscriptstyle
T}}}$. Similarly, let $C=(c_{1},\ldots,c_{M})^{{\mathrm{\scriptscriptstyle
T}}}$.

Two types of covariates suggest two schemes of cluster rerandomization. The
first one is based on the Mahalanobis distance of covariates at the cluster
level. Define

${\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}}=M_{1}^{-1}\sum_{i=1}^{M}Z_{i}c_{i}-M_{0}^{-1}\sum_{i=1}^{M}(1-Z_{i})c_{i}$

as the Horvitz–Thompson estimator applied to the cluster-level covariates. For
a pre-determined threshold $a>0$, cluster rerandomization based on cluster-
level covariates accepts the treatment assignment if and only if the following
event happens:

$\mathcal{M}_{c}=\\{{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}}^{{\mathrm{\scriptscriptstyle
T}}}\operatorname{cov}({\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}})^{-1}\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\leq a\\}=\\{e_{1}e_{0}M\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}^{{\mathrm{\scriptscriptstyle
T}}}\operatorname{cov}_{\textnormal{f}}({C})^{-1}\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\leq a\\}.$

The second one is based on the Mahalanobis distance of covariates at the
individual level. Define

${\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}=N_{1}^{-1}\sum_{i=1}^{M}Z_{i}\sum_{j=1}^{n_{i}}x_{ij}-N_{0}^{-1}\sum_{i=1}^{M}(1-Z_{i})\sum_{j=1}^{n_{i}}x_{ij}$

as the difference in means of the individual-level covariates. For a pre-
determined threshold $a>0$, cluster rerandomization based on individual-level
covariates accepts the treatment assignment if and only if the following event
happens:

$\mathcal{M}_{x}=\\{{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{{\mathrm{\scriptscriptstyle
T}}}\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})^{-1}\
{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\leq a\\}=\\{e_{1}e_{0}M\
{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{{\mathrm{\scriptscriptstyle
T}}}\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})^{-1}{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\leq
a\\}.$ (3)

We use the same threshold $a$ for two types of cluster rerandomization to
simplify the notation. Due to the complicated form of the exact covariance of
${\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}$, we have used its asymptotic
analog in (3), with more details given in Proposition 1 in the next section.

The definition of cluster-level or individual-level covariates can be subtle,
and different levels of covariates could be used interchangeably. On the one
hand, we could take some statistics of aggregated individual-level covariates
within a cluster to construct cluster-level covariates. For instance, we could
define ${c}_{i}=\tilde{{x}}_{i\cdot}$. On the other hand, it is common that
only cluster-level covariates $c_{i}$ are observed prior to rerandomization,
in which case we can define individual-level covariates as $x_{ij}=c_{i}$ for
all $j$ in cluster $i$. Hence, both cluster rerandomization schemes are
applicable.

###  3.2 Asymptotic distributions

We have introduced two estimators
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ for $\tau$, and two cluster
rerandomization schemes $\mathcal{M}_{c}$ and $\mathcal{M}_{x}$.
Theoretically, they generate four combinations. Nevertheless, it is more
natural to consider the design and analysis with the same type. Therefore, we
will focus on the asymptotic distributions of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ given $\mathcal{M}_{c}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ given $\mathcal{M}_{x}$. We
start with unconditional joint asymptotic distributions of
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}},{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})$
and
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}},{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}})$ as shown in Proposition 1 below.

######  Proposition 1.

Under regularity conditions,

$\displaystyle
M^{1/2}\left(\begin{array}[]{c}\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau\\\
{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\end{array}\right)\mathrel{\
\dot{\sim}\
}\mathcal{N}\left(0,\begin{bmatrix}V_{\textnormal{haj},\tau\tau}&V_{\textnormal{haj},\tau
x}\\\ V_{\textnormal{haj},x\tau}&V_{\textnormal{haj},xx}\end{bmatrix}\right),$
$\displaystyle
M^{1/2}\left(\begin{array}[]{c}\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau\\\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\end{array}\right)\mathrel{\
\dot{\sim}\
}\mathcal{N}\left(0,\begin{bmatrix}V_{\textnormal{ht},\tau\tau}&V_{\textnormal{ht},\tau
c}\\\ V_{\textnormal{ht},c\tau}&V_{\textnormal{ht},cc}\end{bmatrix}\right),$

where $V_{\textnormal{ht},\tau\tau}$ and $V_{\textnormal{haj},\tau\tau}$ are
defined in (1) and (2), respectively, and

$\displaystyle V_{\textnormal{haj},x\tau}$ $\displaystyle=$
$\displaystyle(V_{\textnormal{haj},\tau x})^{{\mathrm{\scriptscriptstyle
T}}}=e_{1}^{-1}\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}(1)\\}+e_{0}^{-1}\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}(0)\\},\quad
V_{\textnormal{haj},xx}=(e_{1}e_{0})^{-1}\operatorname{cov}_{\textnormal{f}}(\tilde{{X}}),$
$\displaystyle V_{\textnormal{ht},c\tau}$ $\displaystyle=$
$\displaystyle(V_{\textnormal{ht},\tau c})^{{\mathrm{\scriptscriptstyle
T}}}=e_{1}^{-1}\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}(1)\\}+e_{0}^{-1}\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}(0)\\},\quad
V_{\textnormal{ht},cc}=(e_{1}e_{0})^{-1}\operatorname{cov}_{\textnormal{f}}({C}).$

By Proposition 1, the Mahalanobis distances based on
${\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}$ and
${\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}$ converge in distribution to
$\chi_{K}^{2}$. This provides a guidance for choosing $a$. For example, we can
choose $a$ as the $\alpha$th quantile of $\chi_{K}^{2}$ to ensure an
asymptotic acceptance rate of $\alpha$. Morgan and Rubin, (2012) suggested
$\alpha=0.001$ when the cluster numbers are moderate or large. However,
practitioners need to be aware that with small $M$ and $a$, we do not have
enough randomizations and asymptotic approximation can be poor. We suggest the
Fisher randomization tests in such situations, where the threshold $a$ is
chosen to ensure non-trivial powers (Johansson et al.,, 2021).

Proposition 1 provides the basis for deriving the conditional asymptotic
distributions, $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ given
$\mathcal{M}_{c}$ and $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ given
$\mathcal{M}_{x}$, which are the asymptotic distributions of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ under cluster
rerandomization. Extending Li et al., (2018), we introduce two squared
multiple correlations:

$\displaystyle R_{c}^{2}$ $\displaystyle=$
$\displaystyle\textnormal{cov}_{\textnormal{a}}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}},{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}})\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}})^{-1}\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}},\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}})/\textnormal{var}_{\textnormal{a}}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}),$
$\displaystyle R_{x}^{2}$ $\displaystyle=$
$\displaystyle\textnormal{cov}_{\textnormal{a}}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}},{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})^{-1}\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}},\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}})/\textnormal{var}_{\textnormal{a}}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}).$

Moreover, let $L_{K,a}\sim D_{1}\mid{D}^{{\mathrm{\scriptscriptstyle
T}}}{D}\leq a$ where ${D}=(D_{1},\ldots,D_{K})^{{\mathrm{\scriptscriptstyle
T}}}$ is a $K$-dimensional standard normal random vector, and let $\epsilon$
be a standard normal random variable independent of $L_{K,a}$. These
quantities play important roles in Theorem 1 below.

######  Theorem 1.

Under regularity conditions,

$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{M}_{x}$
$\displaystyle\mathrel{\ \dot{\sim}\ }$
$\displaystyle(V_{\textnormal{haj},\tau\tau})^{1/2}\\{(1-R_{x}^{2})^{1/2}\epsilon+R_{x}L_{K,a}\\},$
$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau)\mid\mathcal{M}_{c}$
$\displaystyle\mathrel{\ \dot{\sim}\ }$
$\displaystyle(V_{\textnormal{ht},\tau\tau})^{1/2}\\{(1-R_{c}^{2})^{1/2}\epsilon+R_{c}L_{K,a}\\}.$

Theorem 1 extends Li et al., (2018, Theorem 1) for rerandomization with
treatments assigned at the individual level. We comment on the result on
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ in Theorem 1, and analogous
discussion also applies to $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$.
It demonstrates that the asymptotic distribution of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ is a linear combination of
two independent random variables: a standard normal $\epsilon$ and a truncated
normal $L_{k,a}$. The truncated normal random variable is due to the
projection of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ onto
${\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}$, which is affected by cluster
rerandomization. The untruncated normal random variable is due to the residual
from the projection, which is unaffected by cluster rerandomization. By
Theorem 1, rerandomization narrows the quantile ranges of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$, and with proper moment
estimators of the asymptotic distributions, we can construct narrower
confidence intervals for $\tau$.

Theorem 1 also allows for a direct comparison of the two cluster
rerandomization schemes. We construct cluster-level covariates from
individual-level ones by defining ${c}_{i}=\tilde{{x}}_{i\cdot}$ for
$\mathcal{M}_{c}$, while we use the individual-level covariates ${x}_{ij}$ for
$\mathcal{M}_{x}$ to enable a fair comparison. Corollary 1 below demonstrates
the superiority of the cluster rerandomization scheme $\mathcal{M}_{c}$.

######  Corollary 1.

Under regularity conditions, if
${c}_{i}=(n_{i},\tilde{{x}}_{i\cdot}^{{\mathrm{\scriptscriptstyle
T}}})^{{\mathrm{\scriptscriptstyle T}}}$, then
$V_{\textnormal{haj},\tau\tau}(1-R_{x}^{2})\geq
V_{\textnormal{ht},\tau\tau}(1-R_{c}^{2}).$

With a small threshold $a>0$, the variance of $L_{K,a}$ is small, and
therefore, $V_{\textnormal{haj},\tau\tau}(1-R_{x}^{2})$ and
$V_{\textnormal{ht},\tau\tau}(1-R_{c}^{2})$ are the leading terms of the
asymptotic variances of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$, respectively. Corollary 1
demonstrates that based on the asymptotic variances, $\mathcal{M}_{c}$ with
cluster-level covariates dominates $\mathcal{M}_{x}$ with individual-level
covariates if the cluster-level covariates include the cluster size and scaled
cluster totals of the individual-level covariates. A key requirement of
Corollary 1 is that the cluster size must be used as a cluster-level
covariate. This parallels the theory for regression adjustment in the analysis
of cluster-randomized experiments that the regression-adjusted estimator based
on scaled cluster totals dominates the regression-adjusted estimator based on
individual-level data with properly defined covariates (Su and Ding,, 2021).

##  4 Cluster rerandomization with prior knowledge on covariate importance

###  4.1 Weighted Euclidean distance criterion

The cluster rerandomization schemes in Section 3 view all covariates as
equally important. Although they have the advantage of being invariant to non-
degenerate linear transformation of the covariates, they are not ideal in
experiments with prior knowledge about the relative importance of the
covariates. In those cases, a better choice is cluster rerandomization with
the weighted Euclidean distance, which has been frequently used in practice
(Wight et al.,, 2002; Althabe et al.,, 2008; Li et al.,, 2016, 2017; Hayes and
Moulton,, 2017; Dempsey et al.,, 2018). We will study its design-based
properties in this subsection.

Consider cluster rerandomization schemes defined by general quadratic forms of
the measures of covariate imbalance:

$\mathcal{D}_{x}(A_{x})=\\{M{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}\ {\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\leq
a\\},\quad\mathcal{D}_{c}(A_{c})=\\{M{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}}^{{\mathrm{\scriptscriptstyle T}}}A_{c}\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\leq a\\},$

where $A_{x}>0$ and $A_{c}>0$. They reduce to cluster rerandomization schemes
in Section 3 if
$A_{x}=M^{-1}\textnormal{cov}_{\textnormal{a}}({\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})^{-1}$
and $A_{c}=M^{-1}\operatorname{cov}({\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}})^{-1}$, respectively. The weighted Euclidean distances correspond to
diagonal $A_{x}$ and $A_{c}$. Theorem 2 below provides the basis for inference
under $\mathcal{D}_{x}(A_{x})$ and $\mathcal{D}_{c}(A_{c})$.

######  Theorem 2.

Under regularity conditions,

$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{D}_{x}(A_{x})$
$\displaystyle\mathrel{\ \dot{\sim}\
}V_{\textnormal{haj},\tau\tau}^{1/2}\big{\\{}(1-R_{x}^{2})^{1/2}\epsilon+R_{x}\mu_{x}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{1/2}A_{x}V^{1/2}_{\textnormal{haj},xx}{\eta}\leq
a\big{\\}},$ $\displaystyle
M^{1/2}(\hat{\tau}_{\textnormal{ht}}-\tau)\mid\mathcal{D}_{c}(A_{c})$
$\displaystyle\mathrel{\ \dot{\sim}\
}V_{\textnormal{ht},\tau\tau}^{1/2}\big{\\{}(1-R_{c}^{2})^{1/2}\epsilon+R_{c}\mu_{c}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}{V_{\textnormal{ht},cc}^{1/2}A_{c}V^{1/2}_{\textnormal{ht},cc}}{\eta}\leq
a\big{\\}},$

where ${\eta}=(\eta_{1},\ldots,\eta_{K})^{{\mathrm{\scriptscriptstyle T}}}$,
$\epsilon,\eta_{1},\ldots,\eta_{K}$ are independent standard normal random
variables, and

${\mu}_{x}^{{\mathrm{\scriptscriptstyle T}}}=(V_{\textnormal{haj},\tau
x}V^{-1}_{\textnormal{haj},xx}V_{\textnormal{haj},x\tau})^{-1/2}V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2},\quad{\mu}_{c}^{{\mathrm{\scriptscriptstyle
T}}}=(V_{\textnormal{ht},\tau
c}V^{-1}_{\textnormal{ht},cc}V_{\textnormal{ht},c\tau})^{-1/2}V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1/2}.$

Theorem 2 is similar to Theorem 1 in that the asymptotic distributions have
two components: one is normal, and the other is truncated normal. It allows us
to derive properties of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ and
$\hat{\tau}_{\textnormal{ht}}$ under cluster rerandomization.

######  Proposition 2.

Under regularity conditions, (i) the asymptotic distributions in Theorem 2 are
symmetric around zero and unimodal, and (ii)
$\textnormal{pr}_{\textnormal{a}}\\{M^{1/2}|\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau|<\delta\mid\mathcal{D}_{x}(A_{x})\\}$
is a non-decreasing function of $R^{2}_{x}$ and
$\textnormal{pr}_{\textnormal{a}}\\{M^{1/2}|\hat{\tau}_{\textnormal{ht}}-\tau|<\delta\mid\mathcal{D}_{c}(A_{c})\\}$
is a non-decreasing function of ${R}_{c}^{2}$ for any fixed $\delta>0$.

Proposition 2 (i) ensures that the asymptotic distributions in Theorem 2 are
both bell-shaped, although they differ from normal distributions. Proposition
2 (ii) ensures that the asymptotic distributions in Theorem 2 are more
concentrated at zero than those under standard cluster randomization.

To compare the asymptotic efficiency of the two cluster rerandomization
schemes, we can compare their variance reductions given the same acceptance
rate. Let $\alpha$ denote the common asymptotic acceptance rate:

$\alpha=\textup{pr}_{\textup{a}}\\{\mathcal{D}_{x}(A_{x})\\}=\textup{pr}_{\textup{a}}(M\hat{\tau}_{\textnormal{haj},x}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}\hat{\tau}_{\textnormal{haj},x}\leq
a),\quad\alpha=\textup{pr}_{\textup{a}}\\{\mathcal{D}_{c}(A_{c})\\}=\textup{pr}_{\textup{a}}(M{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}}^{{\mathrm{\scriptscriptstyle
T}}}A_{c}{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\leq a).$

Let $\Gamma(\cdot)$ be the Gamma function and

$p_{K}=\frac{2\pi}{K+2}\left\\{\frac{2\pi^{K/2}}{K\Gamma(K/2)}\right\\}^{-2/K}.$

The expansions in Theorem 3 below provide the basis for comparing asymptotic
efficiency.

######  Theorem 3.

Under regularity conditions,

$\displaystyle\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{D}_{x}(A_{x})\\}={V}_{\textnormal{haj},\tau\tau}\\{(1-R^{2}_{x})+{R}^{2}_{x}p_{K}{\nu}_{x}(A_{x})\alpha^{2/K}+o(\alpha^{2/K})\\},$
$\displaystyle\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\textnormal{ht}}-\tau)\mid\mathcal{D}_{c}(A_{c})\\}={V}_{\textnormal{ht},\tau\tau}\\{(1-R^{2}_{c})+{R}^{2}_{c}p_{K}{\nu}_{c}(A_{c})\alpha^{2/K}+o(\alpha^{2/K})\\},$

for a small $\alpha$, where

$\displaystyle\nu_{x}(A_{x})$ $\displaystyle=$
$\displaystyle\frac{V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}A_{x}^{-1}V_{\textnormal{haj},xx}^{-1}V_{\textnormal{haj},x\tau}\det(A_{x})^{1/K}\det(V_{\textnormal{haj},xx})^{1/K}}{V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V_{\textnormal{haj},x\tau}},$
$\displaystyle\nu_{c}(A_{c})$ $\displaystyle=$
$\displaystyle\frac{V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}A_{c}^{-1}V_{\textnormal{ht},cc}^{-1}V_{\textnormal{ht},c\tau}\det(A_{c})^{1/K}\det(V_{\textnormal{ht},cc})^{1/K}}{V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}V_{\textnormal{ht},c\tau}}.$

As a sanity check, the cluster rerandomization schemes using the Mahalanobis
distances lead to ${\nu}_{x}(A_{x})=\nu_{x}(V_{\textnormal{haj},xx}^{-1})=1$
and ${\nu}_{c}(A_{c})=\nu_{c}(V_{\textnormal{ht},cc}^{-1})=1$. Given the same
$R^{2}_{x}$ or $R^{2}_{c}$ in Theorem 3, the asymptotic efficiency of using
$\mathcal{D}_{x}(A_{x})$ or $\mathcal{D}_{c}(A_{c})$ depends only on
$\nu_{x}(A_{x})$ or $\nu_{c}(A_{c})$, respectively. Cluster rerandomization
schemes using the weighted Euclidean distances have diagonal $A_{x}$ and
$A_{c}$, which allows us to derive the optimal weighting matrices as shown in
Theorem 4 below. Let $\xi_{k}$ denote a $K$-dimensional vector with $1$ at the
$k$th dimension and $0$ at other dimensions.

######  Theorem 4.

Under regularity conditions, if ${V}_{\textnormal{haj},\tau
x}{V}_{\textnormal{haj},xx}^{-1}\xi_{k}$ and ${V}_{\textnormal{ht},\tau
c}{V}_{\textnormal{ht},cc}^{-1}\xi_{k}$ are nonzero for all $k=1,\ldots,K$,
then $\nu_{x}\\{{\operatorname{diag}}(w_{1},\ldots,w_{K})\\}$ reaches minimum
if $w_{k}\propto(V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}\xi_{k})^{2}$ for $k=1,\ldots,K$, and
$\nu_{c}\\{{\operatorname{diag}}(w_{1},\ldots,w_{K})\\}$ reaches minimum if
$w_{k}\propto(V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}\xi_{k})^{2}$ for $k=1,\ldots,K$.

Based on Theorem 4, let $A^{\textnormal{opt}}_{x}$ and
$A^{\textnormal{opt}}_{c}$ denote the optimal weighting matrices at the
individual level and cluster level, respectively. Recall Proposition 1 to
obtain that

$\displaystyle V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}=\operatorname{cov}_{\textnormal{f}}\\{e_{1}\tilde{\varepsilon}(0)+e_{0}\tilde{\varepsilon}(1),\tilde{{X}}\\}\\{\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})\\}^{-1},$
$\displaystyle V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}=\operatorname{cov}_{\textnormal{f}}\\{e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1),{C}\\}\\{\operatorname{cov}_{\textnormal{f}}({C})\\}^{-1}.$

Therefore, the optimal weights in $A^{\textnormal{opt}}_{x}$ and
$A^{\textnormal{opt}}_{c}$ correspond to the squared coefficients obtained
from regressing $e_{1}\tilde{\varepsilon}(0)+e_{0}\tilde{\varepsilon}(1)$ on
$\tilde{{X}}$ and regressing $e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)$ on
$\tilde{{C}}$, respectively.

######  Remark 1.

Orthogonalizing the covariates ensures that $V_{\textnormal{haj},xx}^{-1}$ and
$V_{\textnormal{ht},cc}^{-1}$ are both diagonal. Then the optimality of
$A_{x}^{\textnormal{opt}}$ or $A_{c}^{\textnormal{opt}}$ implies
$\nu_{x}(A_{x}^{\textnormal{opt}})\leq\nu_{x}(V_{\textnormal{haj},xx}^{-1})=1$
or
$\nu_{c}(A_{c}^{\textnormal{opt}})\leq\nu_{c}(V_{\textnormal{ht},cc}^{-1})=1$,
respectively. Therefore, with orthogonalized covariates, the cluster
rerandomization schemes based on the optimal weighted Euclidean distances are
superior to those based on the Mahalanobis distances. However, this conclusion
does not hold without orthogonalizing the covariates. We give a counterexample
in the appendix.

###  4.2 Comparison with cluster rerandomization with tiers of covariates

With prior knowledge on the importance of the covariates, Morgan and Rubin,
(2015) proposed rerandomization with tiers of covariates as an alternative to
rerandomization with the weighted Euclidean distance. Morgan and Rubin, (2015)
and Li et al., (2018) derived its properties under the assumption that
treatments are assigned at the individual level. The theory for the cluster
analog is still missing. Moreover, no comparison has been made between these
two rerandomization schemes. We fill these gaps in this subsection.

Similar to Remark 1, the comparison does not yield definite answers with
general covariates. However, we can demonstrate that rerandomization with the
optimal weighted Euclidean distance is superior to rerandomization with tiers
of covariates when the covariates are orthogonalized. In this subsection, we
focus on cluster rerandomization with individual-level covariates; in the
appendix, we give the corresponding results for cluster rerandomization with
cluster-level covariates.

Orthogonalization is applied on the scaled cluster totals. The procedure
proceeds as follows: first, find the upper triangular matrix $U$ by applying
Gram–Schmidt orthogonalization to ${\tilde{x}}_{i\cdot}$ such that
$U{\tilde{x}}_{i\cdot}$ has orthogonal components, and then transform the
individual-level covariates to $Ux_{ij}$. To simplify the presentation, we
assume that $x_{ij}=(x_{ij1},\ldots,x_{ijK})^{\mathrm{\scriptscriptstyle T}}$
are already orthogonalized by the aforementioned procedure and are ordered in
decreasing importance.

We first re-stated the results for cluster rerandomization with orthogonalized
covariates. Let

$V_{\textnormal{haj},xx}=\operatorname{diag}(V_{\textnormal{haj},x_{1}x_{1}},\ldots,V_{\textnormal{haj},x_{K}x_{K}}),\quad
V_{\textnormal{haj},\tau x}=(V_{\textnormal{haj},\tau
x_{1}},\ldots,V_{\textnormal{haj},\tau x_{K}}).$

Define the squared multiple correlation for the $k$th covariate as

$R_{x_{k}}^{2}=V_{\textnormal{haj},\tau
x_{k}}^{2}/(V_{\textnormal{haj},\tau\tau}V_{\textnormal{haj},x_{k}x_{k}}),\quad
k=1,\ldots,K.$

Then by Theorem 4, the cluster rerandomization scheme based on the optimal
weighted Euclidean distance is
$\mathcal{D}_{x}(A_{x}^{\textnormal{opt}})=\\{M{\hat{\tau}_{\textnormal{haj},x}}^{\mathrm{\scriptscriptstyle
T}}A_{x}^{\textnormal{opt}}{\hat{\tau}_{\textnormal{haj},x}}\leq a\\}$ with
the optimal diagonal matrix

$A_{x}^{\textnormal{opt}}=\operatorname{diag}\\{(V_{\textnormal{haj},\tau
x_{1}}V_{\textnormal{haj},x_{1}x_{1}}^{-1})^{2},\ldots,(V_{\textnormal{haj},\tau
x_{k}}V_{\textnormal{haj},x_{K}x_{K}}^{-1})^{2}\\}.$

Corollary 2 below is a special case of Theorem 3.

######  Corollary 2.

Under regularity conditions with orthogonalized covariates and optimal
weighted Euclidean distance, if ${V}_{\textnormal{haj},\tau
x_{k}}{V}_{\textnormal{haj},x_{k}x_{k}}^{-1}$ are nonzero for all
$k=1,\ldots,K$, then

$\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{D}_{x}(A_{x}^{\textnormal{opt}})\\}=V_{\textnormal{haj},\tau\tau}\Big{\\{}(1-R_{x}^{2})+K\Big{(}\prod_{k=1}^{K}R_{x_{k}}^{2}\Big{)}^{1/K}p_{K}\alpha^{2/K}+o(\alpha^{2/K})\Big{\\}},$

for a small $\alpha$.

By Theorem 3, we can obtain that

${\textnormal{var}_{\textnormal{a}}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{M}_{x}\\}=V_{\textnormal{haj},\tau\tau}\\{(1-R_{x}^{2})+R_{x}^{2}p_{K}\alpha^{2/K}+o(\alpha^{2/K})\\}.$

Comparing it with Corollary 2 above, we can see that cluster rerandomization
with the optimal weighted Euclidean distance is indeed superior to cluster
rerandomization with the Mahalanobis distance due to
$K\big{(}\prod_{k=1}^{K}R_{x_{k}}^{2}\big{)}^{1/K}\leq R_{x}^{2}$ by the
inequality of arithmetic and geometric means, which echoes Remark 1.

In rerandomization with tiers of covariates, Morgan and Rubin, (2015)
partitioned the covariates into $L$ tiers by their importance with $K_{l}$
covariates in tier $l$, for $l=1,\ldots,L$ and $\sum_{l=1}^{L}K_{l}=K$. Let
${x}_{[l]}$ be the covariates in the $l$th tier. Then
${x}=({x}_{[1]},\ldots,{x}_{[L]})$. Morgan and Rubin, (2015) originally used
the block-wise Gram–Schmidt orthogonalization, but we use the above element-
wise orthogonalization because it is notationally simpler and leads to the
same asymptotic distribution. Within tier $l$, define
$\hat{\tau}_{\textnormal{haj},x_{[l]}}$ as the difference in means of
covariates, $V_{\textnormal{haj},x_{[l]}x_{[l]}}$ as its covariance,
$\mathcal{M}_{[l]}$ as the event that the corresponding Mahalanobis distance
is smaller than or equal to the threshold $a_{[l]}$, $R_{x_{[l]}}^{2}$ as the
corresponding $R^{2}$, and $\alpha_{[l]}$ as the asymptotic acceptance rate.
The overall asymptotic acceptance rate satisfies
$\alpha=\prod_{l=1}^{L}\alpha_{[l]}$ due to the orthogonalization of the
covariates.

######  Corollary 3.

Under regularity conditions,

$\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{M}_{[1]},\ldots,\mathcal{M}_{[L]}\\}=V_{\textnormal{haj},\tau\tau}\Big{\\{}(1-R_{x}^{2})+\sum_{l=1}^{L}R_{x_{[l]}}^{2}p_{K_{l}}\alpha_{[l]}^{2/K_{l}}+o(\alpha_{[l]}^{2/K_{l}})\Big{\\}},$

for small $\alpha_{[l]}$, $l=1,\ldots,L$.

Given $\alpha=\prod_{l=1}^{L}\alpha_{[l]}>0$, to minimize the second term in
the asymptotic variance in Corollary 3, we must choose

$\alpha_{[l]}=\left({c_{0}{R}_{x_{[l]}}^{2}p_{K_{l}}}/{K_{l}}\right)^{-K_{l}/2}\quad(l=1,\ldots,L),$

for some positive constant $c_{0}$. Theorem 5 below compares the second terms
in the asymptotic variances in Corollaries 2 and 3.

######  Theorem 5.

Under regularity conditions with orthogonalized covariates, the following
inequality holds for any $\alpha_{[1]},\ldots,\alpha_{[L]}$ satisfying
$\alpha=\prod_{l=1}^{L}\alpha_{[l]}$:

$\displaystyle\sum_{l=1}^{L}{R}_{x_{[l]}}^{2}p_{K_{l}}\alpha_{[l]}^{2/K_{l}}\geq
K\Big{(}\prod_{k=1}^{K}R_{x_{k}}^{2}\Big{)}^{1/K}p_{K}\alpha^{2/K}.$

Based on the comparison of the asymptotic variances, Theorem 5 quantifies the
superiority of cluster rerandomization with the optimal weighted Euclidean
distance when the covariates are orthogonalized.

##  5 Rerandomization and regression adjustment

Rerandomization uses covariates in the design stage (Morgan and Rubin,, 2012),
and regression adjustment uses covariates in the analysis stage (Lin,, 2013).
Li et al., (2018) pointed out that they are dual in the design and analysis of
randomized experiments with treatments assigned at the individual level.
Moreover, Li and Ding, (2020) showed that they could be used simultaneously.
In this section, we will show that analogous results hold under cluster
rerandomization and point out some subtle differences.

Under $\mathcal{D}_{x}(A_{x})$, we can use the coefficient of $Z_{ij}$ in the
least squares fit of $Y_{ij}$ on $(1,Z_{ij},x_{ij},Z_{ij}x_{ij})$ to estimate
$\tau$, and use the cluster-robust standard error to approximate the true
asymptotic standard error (Liang and Zeger,, 1986); under
$\mathcal{D}_{c}(A_{c})$, we can use the coefficient of $Z_{i}$ in the least
squares fit of $\tilde{Y}_{i\cdot}$ on $(1,Z_{i},c_{i},Z_{i}c_{i})$ to
estimate $\tau$, and use the heteroskedasticity-robust standard error to
approximate the true asymptotic standard error (Huber,, 1967; White,, 1980).
Denote the two sets of regression coefficient and variance estimator by
$({\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}},\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}})$
and
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}},\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}})$,
respectively. Let $z_{\varsigma}$ ($0<\varsigma<1$) be the $\varsigma$th
quantile of a standard normal distribution. Theorem 6 below summarizes the
important results about cluster rerandomization combined with regression
adjustment.

######  Theorem 6.

Assume regularity conditions hold.

  * (i)

Under $\mathcal{D}_{c}(A_{c})$, the estimator
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}$
is consistent for $\tau$ and asymptotically normal, the probability limit of
$M\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}}$ is larger
than or equal to the true asymptotic variance of
$M^{1/2}\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}$,
and the $1-\varsigma$ confidence interval

$\left[\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}})^{1/2}z_{\varsigma/2},\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}})^{1/2}z_{1-\varsigma/2}\right]$

has asymptotic coverage rate $\geq 1-\varsigma$;

  * (ii)

Under $\mathcal{D}_{x}(A_{x})$, the estimator
${\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}$
is consistent for $\tau$ and its asymptotic distribution is a convolution of
normal and truncated normal, the probability limit of
$M\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}}$ is larger
than or equal to the true asymptotic variance of
$M^{1/2}{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}$,
and the $1-\varsigma$ confidence interval

$\left[{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}})^{1/2}z_{\varsigma/2},{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}})^{1/2}z_{1-\varsigma/2}\right]$

has asymptotic coverage rate $\geq 1-\varsigma$;

  * (iii)

If ${c}_{i}=(n_{i},\tilde{{x}}_{i\cdot}^{{\mathrm{\scriptscriptstyle
T}}})^{{\mathrm{\scriptscriptstyle T}}}$, the asymptotic distribution of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}\mid\mathcal{D}_{c}(A_{c})$
is more concentrated at $\tau$ than
${\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}\mid\mathcal{D}_{x}(A_{x})$,
in the sense that for any $\delta>0$, we have

$\displaystyle\textnormal{pr}_{\textnormal{a}}\\{M^{1/2}|{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}-\tau|<\delta\mid\mathcal{D}_{x}(A_{x})\\}\leq\textnormal{pr}_{\textnormal{a}}\\{M^{1/2}|\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau|<\delta\mid\mathcal{D}_{c}(A_{c})\\}.$

To avoid notation complexity, we relegate the details of the asymptotic
distributions of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}$
given $\mathcal{D}_{c}(A_{c})$ and
${\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}$
given $\mathcal{D}_{x}(A_{x})$ to Theorem 7 in the appendix. Theorem 6
parallels the key results in Su and Ding, (2021) without rerandomization. It
justifies the standard Wald-type inference based on appropriate regressions
and standard errors, which can be conveniently obtained by standard
statistical software packages. It also extends the results in Li and Ding,
(2020) from classic rerandomization to cluster rerandomization. However, some
asymptotic results differ from those under rerandomization with treatments
assigned at the individual level. In particular, with cluster rerandomization
and regression adjustment using individual-level covariates, the asymptotic
distribution is not normal anymore; see Theorem 7 in the appendix for more
details. Moreover, Theorem 6 (iii) extends the theory in Su and Ding, (2021)
that regression adjustment using cluster size and scaled cluster totals
outperforms regression adjustment using individual-level covariates.

######  Remark 2.

Li and Ding, (2020) also discussed the scenarios that the designer and
analyzer do not communicate well, allowing the covariates used in the analysis
stage be different from those used in the design stage. One important special
case is that the set of covariates used in the analysis stage is empty. To be
specific, under $\mathcal{D}_{x}(A_{x})$, we can use the coefficient of
$Z_{ij}$ in the least squares fit of $Y_{ij}$ on $(1,Z_{ij})$ to estimate
$\tau$, and use the cluster-robust standard error to approximate the true
asymptotic standard error; under $\mathcal{D}_{c}(A_{c})$, we can use the
coefficient of $Z_{i}$ in the least squares fit of $\tilde{Y}_{i\cdot}$ on
$(1,Z_{i})$ to estimate $\tau$, and use the heteroskedasticity-robust standard
error to approximate the true asymptotic standard error. These point
estimators correspond to the Hajek and Horvitz–Thompson estimators discussed
in Section 2, with standard errors
$(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}})^{1/2}$ and
$(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}})^{1/2}$,
respectively. Analogous to Theorem 6, these point estimators are still
consistent for $\tau$ and converge in distribution to convolutions of normal
and truncated normal under cluster rerandomization. Moreover, the normal-based
confidence intervals are still asymptotically conservative. However, different
from the results in Su and Ding, (2021), the standard errors
$(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}})^{1/2}$ and
$(\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}})^{1/2}$ are
overly conservative because they ignore the truncated normal components in the
asymptotic distributions. Due to the technical complexities, we relegate the
theory and improved inference methods to the appendix.

######  Remark 3.

Our asymptotic theory requires a large number of clusters, which may be
unrealistic in many cluster randomized experiments. When $M$ is small, one
alternative approach is to use a mixed-effects model by imposing modelling
assumptions on the data generating process. Another alternative approach is to
use Fisher randomization tests that deliver finite-sample exact $p$-values
under the sharp null hypothesis. See Zhao and Ding, (2021, Section 4) for more
details. Importantly, the Fisher randomization tests must follow the same
treatment assignment rule of cluster rerandomization.

##  6 Numeric Examples

###  6.1 Simulation

In this section, we conduct simulation to assess the finite-sample
performances of all eight combinations with binary choices for three
orthogonal axes, individual-level (X) versus cluster-level (C), the
Mahalanobis distance (M) versus the optimal weighted Euclidean distance
without orthogonalization (W), and using regression adjustment (.adj) or not.
We also consider two baseline methods: Hajek (Haj) and Horvitz–Thompson (HT)
estimators without using cluster rerandomization for better comparison. The
names of the methods are written as ReMC, ReWC, ReMI, ReWI, ReMC.adj,
ReWC.adj, ReMI.adj, and ReWI.adj, where the first two letters refer to cluster
rerandomization and the rest correspond to the aforementioned three axes.

We generate potential outcomes from the following model:

$\displaystyle Y_{ij}(z)=g(n_{i})+{x}^{\mathrm{\scriptscriptstyle
T}}_{ij}\beta_{iz}+\varepsilon_{ij}(z)\quad(i=1,\ldots,M,\ j=1,\ldots,n_{i},\
z=0,1).$

Here $g(n_{i})$ captures the cluster effects on individuals and
$\varepsilon_{ij}(z)$’s are independent $\mathcal{N}(0,16)$. Pre-treatment
covariates $x_{ij}$ are generated from a $K$-dimensional normal distribution
with mean $0$ and covariance matrix $\Sigma$, where
$\Sigma_{ij}=(1-\rho)\delta_{ij}+\rho$, with $\delta_{ij}=1$ if $i=j$ and $0$
otherwise. We fix $M=100$, with $50$ assigned to the treatment arm and $50$ to
the control arm. We also conduct simulations with relatively small $M$, and
assess the validity of such asymptotic inference. The results are shown in the
appendix. The size of each cluster is sampled uniformly from
$\\{m\in\mathbb{N}\mid 4\leq m\leq 10\\}$. The coefficients are generated from
$\beta_{iz}=\beta_{z}+U(-0.1,0.1)$, where $U(\cdot,\cdot)$ denotes uniform
distribution. Here each component of $\beta_{1}$ is sampled from
$\\{0.5\gamma,\gamma,1.5\gamma\\}$ with equal probabilities, and
$\beta_{0}=2\gamma{1}_{K}-\beta_{1}$, where $1_{K}$ is a $K$-dimensional all-
one vector, and $\gamma$ is chosen to ensure the proportion of variance
contributed by the covariates to be $0.5$. These $K$ covariates are used to
generate the potential outcomes, conduct cluster rerandomization, and perform
regression adjustment. We use the scaled cluster totals of individual-level
covariates together with cluster size as cluster-level covariates. Each
cluster rerandomization is conducted with an acceptance rate $\alpha=0.1\%$.
We consider four scenarios with parameters given in Table 1.

Table 1: Parameters of four scenarios Scenario | $K$ | $\rho$ | $\gamma$ | $g(n_{i})$
---|---|---|---|---
1 | $7$ | $0.8$ | $1$ | $(n_{i}-7)/2$
2 | $7$ | $-0.15$ | $5$ | $(n_{i}-7)/2$
3 | $12$ | $0.4$ | $0.5$ | $6$
4 | $12$ | $-0.09$ | $12$ | $6$

Once generated, the potential outcomes and covariates are fixed throughout the
simulation. We repeat the simulation for each scenario, draw treatment
assignments $1000$ times, and report the biases, standard deviations (SD),
root mean squared errors (RMSE), and empirical coverage probabilities (CP) of
95% confidence intervals, and mean lengths of such intervals (CI Length). In
particular, coverage probability is the percentage of confidence intervals
that cover the true value of the average treatment effect from the $1000$
replications. We consider two methods for constructing confidence intervals:
the normal-based method given in Theorem 6 and Remark 2 and improved method
provided in the appendix.

Figure 1: Comparison results in simulation. Table 2: Comparison results in
simulation

Method | Bias | SD | RMSE | Normal-based | Improved
---|---|---|---|---|---
CP | CI Length | CP | CI Length
Scenario 1
ReMC | -0.01 | 0.29 | 0.29 | 1.00 | 2.53 | 0.97 | 1.24
ReWC | 0.00 | 0.23 | 0.23 | 1.00 | 2.53 | 0.98 | 1.07
ReMX | -0.01 | 0.30 | 0.30 | 1.00 | 2.53 | 0.97 | 1.28
ReWX | 0.01 | 0.27 | 0.27 | 1.00 | 2.54 | 0.96 | 1.20
Haj | -0.02 | 0.62 | 0.62 | 0.96 | 2.53 | – | –
HT | -0.03 | 0.62 | 0.62 | 0.95 | 2.53 | – | –
ReMC.adj | 0.00 | 0.22 | 0.22 | 0.98 | 1.04 | – | –
ReWC.adj | 0.00 | 0.22 | 0.22 | 0.97 | 1.05 | – | –
ReMX.adj | 0.00 | 0.27 | 0.27 | 0.98 | 1.29 | 0.96 | 1.19
ReWX.adj | 0.00 | 0.28 | 0.28 | 0.98 | 1.29 | 0.97 | 1.22
Scenario 2
ReMC | 0.00 | 0.24 | 0.24 | 1.00 | 2.18 | 0.98 | 1.14
ReWC | 0.01 | 0.30 | 0.30 | 1.00 | 2.17 | 0.96 | 1.23
ReMX | -0.01 | 0.27 | 0.27 | 1.00 | 2.17 | 0.97 | 1.21
ReWX | 0.00 | 0.32 | 0.32 | 1.00 | 2.18 | 0.95 | 1.31
Haj | -0.01 | 0.38 | 0.38 | 1.00 | 2.17 | – | –
HT | -0.01 | 0.41 | 0.41 | 0.99 | 2.16 | – | –
ReMC.adj | 0.00 | 0.22 | 0.22 | 0.98 | 1.04 | – | –
ReWC.adj | 0.02 | 0.23 | 0.23 | 0.97 | 1.05 | – | –
ReMX.adj | -0.01 | 0.28 | 0.28 | 0.98 | 1.29 | 0.96 | 1.20
ReWX.adj | 0.00 | 0.27 | 0.27 | 0.99 | 1.30 | 0.97 | 1.20
Scenario 3
ReMC | 0.00 | 0.27 | 0.27 | 1.00 | 2.12 | 0.97 | 1.19
ReWC | 0.00 | 0.20 | 0.20 | 1.00 | 2.12 | 0.99 | 1.03
ReMX | -0.01 | 0.24 | 0.24 | 1.00 | 1.58 | 0.97 | 1.07
ReWX | 0.01 | 0.20 | 0.20 | 1.00 | 1.58 | 0.99 | 1.02
Haj | -0.01 | 0.35 | 0.35 | 0.97 | 1.58 | – | –
HT | 0.01 | 0.50 | 0.50 | 0.95 | 2.10 | – | –
ReMC.adj | 0.00 | 0.21 | 0.21 | 0.98 | 0.99 | – | –
ReWC.adj | -0.01 | 0.21 | 0.21 | 0.98 | 1.01 | – | –
ReMX.adj | 0.00 | 0.20 | 0.20 | 1.00 | 1.16 | 0.99 | 1.00
ReWX.adj | -0.01 | 0.21 | 0.21 | 1.00 | 1.16 | 0.98 | 1.01
Scenario 4
ReMC | -0.01 | 0.27 | 0.27 | 1.00 | 4.50 | 0.98 | 1.34
ReWC | 0.01 | 0.41 | 0.41 | 1.00 | 4.49 | 0.97 | 1.76
ReMX | 0.01 | 0.23 | 0.23 | 1.00 | 4.29 | 0.99 | 1.20
ReWX | 0.01 | 0.31 | 0.31 | 1.00 | 4.30 | 0.97 | 1.49
Haj | -0.01 | 0.32 | 0.32 | 1.00 | 4.27 | – | –
HT | 0.01 | 0.49 | 0.49 | 1.00 | 4.48 | – | –
ReMC.adj | 0.00 | 0.20 | 0.20 | 0.98 | 0.99 | – | –
ReWC.adj | 0.00 | 0.21 | 0.21 | 0.98 | 0.99 | – | –
ReMX.adj | 0.00 | 0.20 | 0.20 | 1.00 | 1.16 | 0.99 | 1.01
ReWX.adj | 0.01 | 0.21 | 0.21 | 1.00 | 1.16 | 0.98 | 1.00

As shown in Table 2 and Figure 1, the combination of cluster rerandomization
and regression adjustment results in better performance for all four
scenarios. However, in some occasions, the combination has small improvement
over cluster rerandomization only. This is because using cluster
rerandomization only already has significant improvement, as the standard
deviations are 30%–50% smaller than the two baseline methods. Moreover, the
biases are negligible for all methods.

In Scenarios 1 and 3, cluster rerandomization schemes based on the optimal
weighted Euclidean distances outperform cluster rerandomization schemes based
on the Mahalanobis distances. This is due to the positively correlated
covariates. On the contrary, in Scenarios 2 and 4, cluster rerandomization
schemes based on the Mahalanobis distances are better because the covariates
are negatively correlated. We conclude that the correlation structure between
covariates can severely influence the efficiency of the balancing criterion.

Moreover, cluster rerandomization with cluster-level covariates outperforms
cluster rerandomization with individual-level covariates in Scenarios 1 and 2,
whereas cluster rerandomization with individual-level covariates is slightly
better in Scenarios 3 and 4. To illustrate this, if the potential outcomes are
not influenced by cluster size like in Scenarios 3 and 4, especially if the
potential outcomes have a small variation across units, cluster
rerandomization with individual-level covariates could be better. This is
because they are based on the Hajek estimator, which partially eliminates the
effect of cluster size by using $\tilde{{\varepsilon}}$ instead of $Y$ and
fits a more precise model for this particular situation. On the other hand, we
recommend cluster rerandomization with cluster-level covariates, especially
when the potential outcomes depend heavily on cluster size, because in such
settings, the effect of cluster size is not directly accounted for if we use
cluster rerandomization with individual-level covariates. Finally, by
Corollary 1, with a sufficient small threshold, cluster rerandomization with
cluster-level covariates is more efficient than cluster rerandomization with
individual-level covariates in terms of variance reduction.

###  6.2 Child growth monitoring data

Tembo et al., (2017) conducted a cluster-randomized experiment in developing
countries to study the effect of community-based monitoring on child growth
faltering. We use part of this dataset to illustrate our theoretical results.
We use two arms of the study, community-based monitoring as treatment and no
monitoring as control. In the experiment, $41$ communities were randomly
assigned to the treatment arm while $41$ to the control arm. The mean, median,
and maximum of cluster size are $3.085$, $3$, and $8$, respectively. We choose
$6$ individual-level covariates, including the child age, family background
information, and some baseline growth data. We use cluster size and scaled
cluster totals of individual-level covariates as cluster-level covariates. The
height of a child is used as the outcome. As the potential outcomes are not
completely available, we fit two simple linear regression models to impute the
missing potential outcomes, and perform simulation based on the imputed
dataset. The overall acceptance rate of all cluster rerandomization schemes
are approximately $0.1\%$, computed from the quantile of the empirical
distances. We use $4$ covariates at the design stage for cluster
rerandomization and $6$ covariates at the analysis stage for regression
adjustment.

Figure 2: Comparison results on the child monitoring data.

Table 3 and Figure 2 show that cluster rerandomization always improves the
performance, and cluster rerandomization schemes based on the weighted
Euclidean distances are better than cluster rerandomization schemes based on
the Mahalanobis distances. In the case when cluster-level covariates are used,
the standard deviations under cluster rerandomization schemes using the
weighted Euclidean distances are about half of those using the Mahalanobis
distances. The combinations of cluster rerandomization and regression
adjustment perform the best, and they have almost the same standard deviations
and mean squared errors.

Table 3: Comparison results on the child monitoring data Method | Bias | SD | RMSE | Normal-based | Improved
---|---|---|---|---|---
CP | CI Length | CP | CI Length
ReMC | -0.01 | 0.25 | 0.25 | 1.00 | 4.63 | 0.94 | 0.95
ReWC | -0.01 | 0.14 | 0.14 | 1.00 | 4.63 | 0.93 | 0.54
ReMX | -0.01 | 0.13 | 0.13 | 0.99 | 0.72 | 0.94 | 0.51
ReWX | 0.00 | 0.13 | 0.13 | 0.99 | 0.72 | 0.95 | 0.52
Haj | 0.00 | 0.17 | 0.17 | 0.96 | 0.72 | – | –
HT | -0.05 | 1.20 | 1.20 | 0.94 | 4.60 | – | –
ReMC.adj | 0.00 | 0.07 | 0.07 | 0.97 | 0.31 | – | –
ReWC.adj | 0.00 | 0.07 | 0.07 | 0.97 | 0.31 | – | –
ReMX.adj | 0.00 | 0.06 | 0.06 | 0.99 | 0.34 | 0.99 | 0.31
ReWX.adj | 0.00 | 0.06 | 0.06 | 0.99 | 0.34 | 0.98 | 0.32

For this particular dataset, the covariates and outcomes are all positive. All
cluster-level covariates and outcomes are strongly correlated with cluster
size and there exists strong collinearity between covariates, which helps to
explain why cluster rerandomization with the weighted Euclidean distance is
much better than cluster rerandomization with the Mahalanobis distance. For
cluster rerandomization with cluster-level covariates, the existence of
collinearity and the fact that the weighted Euclidean distance puts more
weights on more important covariates give rise to a dramatic increase of
efficiency. The outcomes in this dataset, children’s heights, have a small
variation across units, which helps to explain why the Hajek estimator
performs much better than the Horvitz–Thompson estimator.

## Acknowledgement

Dr. Hanzhong Liu was supported by the National Natural Science Foundation of
China, Grant No. 12171476. Dr. Peng Ding was partially supported by the U.S.
National Science Foundation.

## References

  * Althabe et al., (2008) Althabe, F., Buekens, P., Bergel, E., Belizán, J. M., Campbell, M. K., Moss, N., Hartwell, T., and Wright, L. L. (2008).  A behavioral intervention to improve obstetrical care.  N. Eng. J. Med., 358:1929–40.
  * Alzer, (1997) Alzer, H. (1997).  On some inequalities for the gamma and psi functions.  Math. Comput., 66:373–89.
  * Athey and Imbens, (2017) Athey, S. and Imbens, G. W. (2017).  The Econometrics of Randomized Experiments.  North-Holland, Amsterdam.
  * Das Gupta et al., (1972) Das Gupta, S., Eaton, M. L., Olkin, I., Perlman, M., Savage, L. J., and Sobel, M. (1972).  Inequalities on the probability content of convex regions for elliptically contoured distributions.  In Proc. Sixth Berkeley Symp. Math. Statist. Prob., volume 2. Berkeley, California: University of California Press, pp. 241–65.
  * de Hoop et al., (2012) de Hoop, E., Teerenstra, S., van Gaal, B. G., Moerbeek, M., and Borm, G. F. (2012).  The “best balance” allocation led to optimal balance in cluster-controlled trials.  J. Clin. Epidemiol., 65:132–37.
  * Dempsey et al., (2018) Dempsey, A. F., Pyrznawoski, J., Lockhart, S., Barnard, J., Campagna, E. J., Garrett, K., Fisher, A., Dickinson, L. M., and O’Leary, S. T. (2018).  Effect of a health care professional communication training intervention on adolescent human papillomavirus vaccination: a cluster randomized clinical trial.  JAMA Pediatr., 172:e180016–e180016.
  * Dharmadhikari and Joag-Dev, (1988) Dharmadhikari, S. and Joag-Dev, K. (1988).  Unimodality, Convexity, and Applications.  Elsevier, Amsterdam.
  * Donner and Klar, (2000) Donner, A. and Klar, N. (2000).  Design and Analysis of Cluster Randomization Trials in Health Research.  Arnold Publishers Limited, London.
  * Hayes and Moulton, (2017) Hayes, R. J. and Moulton, L. H. (2017).  Cluster Randomised Trials.  CRC Press, Florida.
  * Huber, (1967) Huber, P. J. (1967).  The behavior of maximum likelihood estimates under nonstandard conditions.  In Proc. Fifth Berkeley Symp. Math. Statist. Prob., volume 1. Berkeley, California: University of California Press, pp. 221–33.
  * Johansson et al., (2021) Johansson, P., Rubin, D. B., and Schultzberg, M. (2021).  On optimal rerandomization designs.  J. R. Statist. Soc. B, 83:395–403.
  * Li et al., (2016) Li, F., Lokhnygina, Y., Murray, D. M., Heagerty, P. J., and DeLong, E. R. (2016).  An evaluation of constrained randomization for the design and analysis of group-randomized trials.  Statist. Med., 35:1565–79.
  * Li et al., (2017) Li, F., Turner, E. L., Heagerty, P. J., Murray, D. M., Vollmer, W. M., and DeLong, E. R. (2017).  An evaluation of constrained randomization for the design and analysis of group-randomized trials with binary outcomes.  Statist. Med., 36:3791–806.
  * Li and Ding, (2017) Li, X. and Ding, P. (2017).  General forms of finite population central limit theorems with applications to causal inference.  J. Am. Statist. Assoc., 112:1759–69.
  * Li and Ding, (2020) Li, X. and Ding, P. (2020).  Rerandomization and regression adjustment.  J. R. Statist. Soc. B, 82:241–68.
  * Li et al., (2018) Li, X., Ding, P., and Rubin, D. B. (2018).  Asymptotic theory of rerandomization in treatment-control experiments.  Proc. Nat. Acad. Sci., 115:9157–62.
  * Li et al., (2020) Li, X., Ding, P., and Rubin, D. B. (2020).  Rerandomization in $2^{\textup{K}}$ factorial experiments.  Ann. Statist., 48:43–63.
  * Liang and Zeger, (1986) Liang, K.-Y. and Zeger, S. L. (1986).  Longitudinal data analysis using generalized linear models.  Biometrika, 73:13–22.
  * Lin, (2013) Lin, W. (2013).  Agnostic notes on regression adjustments to experimental data: Reexamining Freedman’s critique.  Ann. Appl. Statist., 7:295–318.
  * Middleton and Aronow, (2015) Middleton, J. A. and Aronow, P. M. (2015).  Unbiased estimation of the average treatment effect in cluster-randomized experiments.  Statist. Polit. Policy, 6:39–75.
  * Morgan and Rubin, (2012) Morgan, K. L. and Rubin, D. B. (2012).  Rerandomization to improve covariate balance in experiments.  Ann. Statist., 40:1263–82.
  * Morgan and Rubin, (2015) Morgan, K. L. and Rubin, D. B. (2015).  Rerandomization to balance tiers of covariates.  J. Am. Statist. Assoc., 110:1412–21.
  * Moulton, (2004) Moulton, L. H. (2004).  Covariate-based constrained randomization of group-randomized trials.  Clin. Trials, 1:297–305.
  * Raab and Butcher, (2001) Raab, G. M. and Butcher, I. (2001).  Balance in cluster randomized trials.  Statist. Med., 20:351–65.
  * Raudenbush, (1997) Raudenbush, S. W. (1997).  Statistical analysis and optimal design for cluster randomized trials.  Psychol. Methods, 2:173–85.
  * Schochet, (2013) Schochet, P. Z. (2013).  Estimators for clustered education RCTs using the Neyman model for causal inference.  J. Educ. Behav. Statist., 38:219–38.
  * Schochet, (2020) Schochet, P. Z. (2020).  Analyzing grouped administrative data for RCTs using design-based methods.  J. Educ. Behav. Statist., 45:32–57.
  * Su and Ding, (2021) Su, F. and Ding, P. (2021).  Model-assisted analyses of cluster-randomized experiments.  J. R. Statist. Soc. B, 83:994–1015.
  * Tembo et al., (2017) Tembo, Sarah, Levenson, Rachel, Rockers, Peter, C., Fink, and Gunther (2017).  Home- and community-based growth monitoring to reduce early life growth faltering: an open-label, cluster-randomized controlled trial.  Am. J. Clin. Nutr., 106:1070–77.
  * (30) Turner, E. L., Li, F., Gallis, J. A., Prague, M., and Murray, D. M. (2017a).  Review of recent methodological developments in group-randomized trials: part 1—design.  Am. J. Public Health, 107:907–15.
  * (31) Turner, E. L., Prague, M., Gallis, J. A., Li, F., and Murray, D. M. (2017b).  Review of recent methodological developments in group-randomized trials: part 2—analysis.  Am. J. Public Health, 107:1078–86.
  * White, (1980) White, H. (1980).  A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  Econometrica, 48:817–38.
  * Wight et al., (2002) Wight, D., Raab, G. M., Henderson, M., Abraham, C., Buston, K., Hart, G., and Scott, S. (2002).  Limits of teacher delivered sex education: interim behavioural outcomes from randomised trial.  Br. Med. J., 324:1430–35.
  * Zhao and Ding, (2021) Zhao, A. and Ding, P. (2021).  Covariate-adjusted fisher randomization tests for the average treatment effect.  J. Econom., page in press.

##  Appendix A Additional Simulation

Let $I_{K}$ denote a $K\times K$ identity matrix. We generate potential
outcomes from the following model:

$Y_{ij}(z)=f(2+\epsilon_{iz}+x_{ij}^{\mathrm{\scriptscriptstyle
T}}\beta_{iz})+\varepsilon_{ij}(z)\quad(i=1,\ldots,M,\ j=1,\ldots,n_{i},\
z=0,1),$

where $\beta_{iz}=\beta_{z}+\varsigma_{iz}$ with the components of $\beta_{z}$
$(z=0,1)$ being independently generated from the $t$-distribution with three
degrees of freedom and $\varsigma_{iz}$ $(i=1,\ldots,M,z=0,1)$ being
independently generated from $\mathcal{N}(0,I_{K})$. We independently generate
$x_{ij}$ $(i=1,\ldots,M,\ j=1,\ldots,n_{i})$ from $\mathcal{N}(0,I_{K})$ and
$\epsilon_{iz}$ $(i=1,\ldots,M,z=0,1)$ from $\mathcal{N}(0,1)$. The error
terms $\varepsilon_{ij}(z)$ $(i=1,\ldots,M,\ j=1,\ldots,n_{i},\ z=0,1)$ are
independently generated from $N(0,\sigma^{2}_{z})$, where $\sigma_{z}^{2}$
($z=0,1$) are chosen so that the
$f(2+\epsilon_{iz}+x_{ij}^{\mathrm{\scriptscriptstyle T}}\beta_{iz})$
$(i=1,\ldots,M,\ j=1,\ldots,n_{i})$ constitute half of the variation of
$Y_{ij}(z)$. We use $(\tilde{x}_{i\cdot},n_{i})\in\mathbb{R}^{K+1}$ in cluster
rerandomization using Mahalanobis distance and regression adjustment. We view
simulation as a full factorial experiment and generate data under all
combinations of the following five factors:

  * •

Number of clusters $M=\\{20+4k\mid k\in\mathbb{N},\ 0\leq k\leq 15\\}$;

  * •

Variance of cluster sizes
$\textnormal{vn}=\\{\textnormal{H},\textnormal{L}\\}$. If vn $=$ H, $n_{i}$’s
are generated uniformly from $\\{m\in\mathbb{N}\mid 2\leq m\leq 16\\}$; if vn
$=$ L, $n_{i}$’s are generated uniformly from $\\{m\in\mathbb{N}\mid 4\leq
m\leq 10\\}$;

  * •

Type of outcome-generating functions
$\textnormal{fn}=\\{\textnormal{linear},\textnormal{nonlinear}\\}$. $f(t)=t$
for fn $=$ linear; $f(t)=t^{3}$ for fn $=$ nonlinear;

  * •

Dimension of covariates $K=\\{1,5\\}$;

  * •

Acceptance rate $\alpha=\\{0.001,0.1\\}$.

Each of the $16\times 2\times 2\times 2\times 2=256$ scenarios is replicated
under $100$ random seeds. In each scenario under each random seed, we draw
assignments for $1000$ times with half of the clusters being assigned to the
treatment arm. We construct $95\%$ normal-based confidence intervals by the
method provided in Section C based on the regression residuals. We compute the
empirical coverage probabilities (CP) and percentages of reduction in root
mean squared errors (RMSE) of the combination of cluster rerandomization and
regression adjustment relative to two benchmark estimators without using
cluster rerandomization, i.e., the Horvitz–Thompson estimator
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ and Hajek estimator
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$ under classic cluster
randomization.

As shown in Figure 3, the asymptotic inference does not vary much between
vn$=$H and vn$=$L and between fn$=$linear and fn$=$nonlinear. The dimension of
covariates plays an important role in the validity of asymptotic inference.
When $K=1$, $M=20$ is enough for the validity of the confidence interval under
half of the random seeds and $M=24$ is enough for the validity of the
confidence interval under most of the random seeds. When $K=5$, $M$ must be
greater than $44$ to ensure the validity of asymptotic inference.

Figure 3: Coverage probabilities for different values of $M$, $K$, and
variance of cluster sizes. Figure 4: Percentages of reduction in RMSE relative
to the Horvitz–Thompson estimator without using cluster rerandomization for
different values of $M$, $K$, and variance of cluster sizes. Figure 5:
Percentages of reduction in RMSE relative to the Hajek estimator without using
cluster rerandomization for different values of $M$, $K$, and variance of
cluster sizes.

Figures 4–5 show the percentages of reduction in RMSE. We can see that the
combination of cluster rerandomization and regression adjustment is beneficial
in reducing RMSE relative to the two benchmark estimators under most of the
random seeds for all values of $M$. The benefit is more evident relative to
the Horvitz–Thompson estimator than to the Hajek estimator.

The asymptotic inference is more likely to be valid and the regression
adjustment provides more efficiency gains under more balanced assignments with
$\alpha=0.001$ than with $\alpha=0.1$. This demonstrates the benefit of
cluster rerandomization.

##  Appendix B More detailed results on cluster rerandomization and regression
adjustment

###  B.1 Asymptotic properties of regression-adjusted estimators under cluster
rerandomization

In general, the covariates used in the analysis stage can be different from
those used in the design stage. In the analysis stage, we use ${w}_{ij}$ for
individual-level covariates and ${v}_{i}$ for cluster-level covariates.
Suppose that all these covariates are centered such that
$\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}{w}_{ij}=0$ and $\sum_{i=1}^{M}{v}_{i}=0$.
Define $\tilde{{w}}_{i\cdot}=\sum_{j=1}^{n_{i}}{w}_{ij}M/N$. Let ${V}$ and
${\tilde{W}}$ be the concatenations of ${v}_{i}$ and $\tilde{w}_{i\cdot}$
respectively.

With a slight abuse of notation, let
${\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}$
denote the coefficient of $Z_{ij}$ obtained from the ordinary least squares
fit of $Y_{ij}$ on $(1,Z_{ij},{w}_{ij},Z_{ij}{w}_{ij})$ in the individual-
level regression adjustment, and
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}$
denote the coefficient of $Z_{i}$ obtained from the ordinary least squares fit
of $\tilde{Y}_{i\cdot}$ on $(1,Z_{i},{v}_{i},Z_{i}{v}_{i})$ in the cluster-
level regression adjustment.

Li and Ding, (2020) derived the asymptotic distributions of the regression-
adjusted estimators under rerandomization with treatments assigned at the
individual level. We extend their results to cluster rerandomization. Define,
for $z=0,1$,

$\displaystyle{\beta}_{w}(z)=$
$\displaystyle\arg\min_{{\beta}}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}\left\\{Y_{ij}(z)-\bar{Y}(z)-{w}_{ij}^{\mathrm{\scriptscriptstyle
T}}{\beta}\right\\}^{2},$ $\displaystyle{\beta}_{v}(z)=$
$\displaystyle\arg\min_{\beta}\sum_{i=1}^{M}\left\\{\tilde{Y}_{i\cdot}(z)-\bar{Y}(z)-{v}_{i}^{\mathrm{\scriptscriptstyle
T}}{\beta}\right\\}^{2}.$

Let $V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$,
$V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau x}$,
$V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},x\tau}$,
$V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$, and
$(R_{x}^{\textnormal{adj}})^{2}$ be the quantities defined similarly to those
in Proposition 1 except that we replace $Y_{ij}(z)$ by
$Y_{ij}(z)-{w}_{ij}^{{\mathrm{\scriptscriptstyle T}}}{\beta}_{w}(z)$, $z=0,1$.
Let $V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$,
$V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau c}$,
$V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},c\tau}$,
$V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ and
$(R_{c}^{\textnormal{adj}})^{2}$ be the quantities defined similarly to those
in Proposition 1 except that we replace $\tilde{Y}_{i\cdot}(z)$ by
$\tilde{Y}_{i\cdot}(z)-{v}_{i}^{{\mathrm{\scriptscriptstyle
T}}}{\beta}_{v}(z)$, $z=0,1$. Their limits exist under regularity conditions
given in Section F.1, and we use the same notation to denote their limits if
no confusion arises. We have the following result.

######  Theorem 7.

Under regularity conditions, cluster rerandomization with individual-level
covariates satisfies

$\displaystyle
M^{1/2}({\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}-\tau)\mid\mathcal{M}_{x}\mathrel{\
\dot{\sim}\
}(V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}})^{1/2}\big{[}\\{1-(R^{\textnormal{adj}}_{x})^{2}\\}^{1/2}\epsilon+R^{\textnormal{adj}}_{x}L_{K,a}\big{]},$
$\displaystyle
M^{1/2}({\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}-\tau)\mid\mathcal{D}_{x}(A_{x})\mathrel{\
\dot{\sim}\
}(V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}})^{1/2}\big{[}\\{1-(R^{\textnormal{adj}}_{x})^{2}\\}^{1/2}\epsilon+R^{\textnormal{adj}}_{x}{\mu}_{x}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2}{\eta}\leq
a\big{]},$

and cluster rerandomization with cluster-level covariates satisfies

$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)\mid\mathcal{M}_{c}\mathrel{\
\dot{\sim}\
}(V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}})^{1/2}\big{[}\\{1-(R^{\textnormal{adj}}_{c})^{2}\\}^{1/2}\epsilon+R^{\textnormal{adj}}_{c}L_{K,a}\big{]},$
$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)\mid\mathcal{D}_{c}(A_{c})\mathrel{\
\dot{\sim}\
}(V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}})^{1/2}\big{[}\\{1-(R^{\textnormal{adj}}_{c})^{2}\\}^{1/2}\epsilon+R^{\textnormal{adj}}_{c}{\mu}_{c}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{ht},cc}^{1/2}A_{c}V_{\textnormal{ht},cc}^{1/2}{\eta}\leq
a\big{]},$

where $\eta=(\eta_{1},\ldots,\eta_{K})^{\mathrm{\scriptscriptstyle T}}$,
$\epsilon$, $\eta_{1},\ldots,\eta_{K}$ are independent standard normal random
variables, and

$\displaystyle{\mu}_{x}^{\mathrm{\scriptscriptstyle
T}}=V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2}/\\{V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},x\tau}\\}^{1/2},\quad{\mu}_{c}^{\mathrm{\scriptscriptstyle
T}}=V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1/2}/\\{V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},c\tau}\\}^{1/2}.$

###  B.2 The gains from cluster rerandomization and regression adjustment

To compare the efficiency, Li and Ding, (2020) introduced the notion of
S-optimality to refer to the coefficient vector ${\beta}$ that could lead to
the shortest asymptotic quantile range among all adjusted estimators. For
cluster rerandomization with cluster-level covariates, the coefficient vector
with respect to S-optimal regression-adjusted average treatment effect
estimator coincides with ${\beta}_{v}(z)$ when ${v}_{i}=B{c}_{i}$ or
${c}_{i}=B{v}_{i}$ for certain matrix $B$. For cluster rerandomization with
individual-level covariates, it may not be ${\beta}_{w}(z)$, which is the
coefficient vector of $w_{ij}$ from the ordinary least squares fit of
$Y_{ij}(z)$ on $(1,{w}_{ij})$. Instead, it equals the coefficient vector of
${\tilde{w}}_{i\cdot}$ from the ordinary least squares fit of
$\tilde{Y}_{i\cdot}(z)$ on $(1,{\tilde{w}}_{i\cdot})$. Su and Ding, (2021)
concluded that the individual-level regression adjustment is suboptimal, so we
will not discuss the efficiency gains from the individual-level regression
adjustment.

Li et al., (2020) proved the efficiency gain for the combination of
rerandomization and regression adjustment with treatments assigned at the
individual level if we use more or less information in the regression
adjustment. We reexamine their results under cluster rerandomization using
cluster-level covariates. We consider two special cases of cluster-level
regression adjustment: ${v}_{i}=B{c}_{i}$ or ${c}_{i}=B{v}_{i}$. The former
indicates that there is less information at the analysis stage, whereas the
latter indicates that there is more information at the analysis stage. For the
latter case, Corollary 4 below shows that regression adjustment can further
improve the efficiency.

######  Corollary 4.

Under regularity conditions, if ${c}_{i}=B{v}_{i}$,

$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)\mid\mathcal{D}_{c}(A_{c})\mathrel{\
\dot{\sim}\ }(V_{\textnormal{ht},\tau\tau})^{1/2}(1-R^{2}_{v})^{1/2}\epsilon,$

where $\epsilon$ is a standard normal random variable, and $R^{2}_{v}$ is the
squared multiple correlation defined similarly to $R_{c}^{2}$ with $c$ being
replaced by $v$.

Corollary 4 can be interpreted as follows. If we use cluster-level covariates
and have more information in the analysis stage, then using the combination of
cluster rerandomization and regression adjustment is asymptotically equivalent
to using regression adjustment only, which is no worse than using cluster
rerandomization only.

For the case we have less information in the analysis stage, we have the
following corollaries.

######  Corollary 5.

Under regularity conditions, if ${v}_{i}=B{c}_{i}$,

$\displaystyle
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)\mid\mathcal{D}_{c}(A_{c})\mathrel{\
\dot{\sim}\
}(V_{\textnormal{ht},\tau\tau})^{1/2}\big{\\{}(1-R^{2}_{c})^{1/2}\epsilon+R_{c}{\mu}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{ht},cc}^{1/2}A_{c}V_{\textnormal{ht},cc}^{1/2}{\eta}\leq
a\big{\\}},$

where $\eta=(\eta_{1},\ldots,\eta_{K})^{\mathrm{\scriptscriptstyle T}}$,
$\epsilon$, $\eta_{1},\ldots,\eta_{K}$ are independent standard normal random
variables,

${\mu}^{\mathrm{\scriptscriptstyle
T}}=(\mu_{1},\ldots,\mu_{K})=\\{V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}V_{\textnormal{ht},c\tau}\\}^{-1/2}V_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1/2}H=\mu_{c}^{\mathrm{\scriptscriptstyle T}}H,$

and

$H=I_{K}-V_{\textnormal{ht},cc}^{1/2}B^{{\mathrm{\scriptscriptstyle
T}}}(BV_{\textnormal{ht},cc}B^{{\mathrm{\scriptscriptstyle
T}}})^{-1}BV_{\textnormal{ht},cc}^{1/2}.$

######  Corollary 6.

Suppose that the covariates are orthogonalized. If we use the first $J<K$
covariates in the regression adjustment, then under regularity conditions,

$\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)\mid\mathcal{D}_{c}(A_{c})\\}\leq\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau)\mid\mathcal{D}_{c}(A_{c})\\}.$

By Corollary 5, the asymptotic distribution of
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)$
under $\mathcal{D}_{c}(A_{c})$ is determined by the direction and length of
the vector ${\mu}$. For an ellipsoidal acceptance region of the cluster
rerandomization criterion, some directions lead to variance reduction while
others do not. Regression adjustment is equivalent to a projection of
${\mu_{c}}$ onto a direction determined by $H$. However, such projection might
not be in a favorable direction and there is no guarantee of efficiency gain
in general. For some special cases as stated in Corollary 6, the asymptotic
variance can be reduced. Moreover, if we use cluster rerandomization based on
the Mahalanobis distance, the acceptance region is a ball and regression
adjustment as a projection decreases the length of the vector ${\mu_{c}}$,
leading to variance reduction. This echoes the result in Li et al., (2020)
that regression adjustment improves the precision when we have less
information in the analysis stage.

The gain from cluster rerandomization can be seen directly from Theorem 7 and
Proposition 2 where $R^{\textnormal{adj}}_{c}$ and $R^{\textnormal{adj}}_{x}$
are equal to $0$ without cluster rerandomization and greater than or equal to
$0$ with cluster rerandomization.

###  B.3 Estimating the sampling distributions of regression-adjusted
estimators

In the main text, we discuss the validity of
$\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}}$ and
$\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}}$ when
$c_{i}=v_{i}$ and $x_{ij}=w_{ij}$. Theorem 8 below states general results
without assuming $c_{i}=v_{i}$ and $x_{ij}=w_{ij}$.

######  Theorem 8.

Under regularity conditions,

$\displaystyle
M\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}}-V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\mid\mathcal{D}_{x}(A_{x})\mathrel{\geq_{\textup{p}}}0,\quad
M\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}}-V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\mid\mathcal{D}_{c}(A_{c})\mathrel{\geq_{\textup{p}}}0,$

where a random variable $\geq_{\textup{p}}0$ means that its probability limit
is greater than or equal to $0$.

By Theorem 8, the two types of robust standard errors are asymptotically
conservative under cluster rerandomization.

##  Appendix C Improvement for variance estimation and confidence intervals

Although $\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{HW}$}}$ and
$\hat{V}^{\textnormal{adj}}_{\scalebox{0.5}{$\textnormal{LZ}$}}$ are valid
under cluster rerandomization, they can be overly conservative partly because
they ignore that the asymptotic distributions have truncated normal components
due to cluster rerandomization. We can derive better variance and confidence
interval estimators by better estimating the sampling distributions.

In Theorem 7, we have derived the asymptotic distributions of
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}-\tau)$
and
$M^{1/2}({\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}-\tau)$
under cluster rerandomization. First, we derive better estimators of
$V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ and
$V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$. Define
the adjusted potential outcomes as
$\tilde{Y}^{\textnormal{adj}}_{i\cdot}(z)=\tilde{Y}_{i\cdot}(z)-{v}_{i}^{{\mathrm{\scriptscriptstyle
T}}}{\beta}_{v}(z)$ and
$\varepsilon_{ij}^{\textnormal{adj}}(z)={\varepsilon}_{ij}(z)-{w}_{ij}^{{\mathrm{\scriptscriptstyle
T}}}{\beta}_{w}(z)$, for $z=0,1$. Let
$\tilde{\varepsilon}^{\textnormal{adj}}_{i\cdot}(z)$ be the scaled cluster
total of ${\varepsilon}_{ij}(z)$. Denote $\tilde{Y}^{\textnormal{adj}}(z)$ and
$\tilde{\varepsilon}^{\textnormal{adj}}(z)$ as the vectors of
$\tilde{Y}^{\textnormal{adj}}_{i\cdot}(z)$ and
$\tilde{\varepsilon}^{\textnormal{adj}}_{i\cdot}(z)$, respectively. Recall
that $V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ and
$V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ are
defined similarly as $V_{\textnormal{haj},\tau\tau}$ and
$V_{\textnormal{ht},\tau\tau}$ except that we replace $\varepsilon_{ij}(z)$ by
$\varepsilon_{ij}^{\textnormal{adj}}(z)$ and $\tilde{Y}_{i\cdot}(z)$ by
$\tilde{Y}^{\textnormal{adj}}_{i\cdot}(z)$, for $z=0,1$.

Define the variance of the linear projection of
$\tilde{Y}^{\textnormal{adj}}(z)$ on ${C}$ as

$\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(z)\mid{C}\\}=\operatorname{cov}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(z),{C}\\}\\{\operatorname{cov}_{\textnormal{f}}({C})\\}^{-1}\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}^{\textnormal{adj}}(z)\\},~{}~{}z=0,1.$

Similarly, define

$\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(z)\mid\tilde{{X}}\\}=\operatorname{cov}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(z),\tilde{{X}}\\}\\{\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})\\}^{-1}\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}^{\textnormal{adj}}(z)\\},~{}~{}z=0,1.$

Let $D$ be the vector of residuals of the cluster-level regression adjustment.
Let ${u}_{ij}$ be the residual of the individual-level regression adjustment,
$\tilde{u}_{i\cdot}$ be its scaled cluster total, and $\tilde{U}$ be the
vector of $\tilde{u}_{i\cdot}$. Denote by
$\hat{\operatorname{var}}_{\textnormal{f,z}}(\cdot)$ and
$\hat{\operatorname{cov}}_{\textnormal{f,z}}(\cdot)$ the sample variance and
covariance under treatment arm $z$. Let $\tilde{G}$ be the union of
$\tilde{X}$ and $\tilde{W}$. To simplify the presentation, we assume that
$\operatorname{cov}_{\textnormal{f}}(\tilde{{G}})$ is invertible.

Following the results of Li et al., (2018), less conservative estimators of
$V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ and
$V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ are

$\displaystyle\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$
$\displaystyle={e_{1}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(0)\\}-\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{G}\\}},$
$\displaystyle\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$
$\displaystyle=e_{1}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(0)\\}-\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid
C\\},$

where
$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(z)\\}=\hat{\operatorname{var}}_{\textnormal{f,z}}(\tilde{U}),$
$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(z)\\}=\hat{\operatorname{var}}_{\textnormal{f,z}}(D),$
for $z=0,1,$ and

$\displaystyle\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}$
$\displaystyle=\big{\\{}\hat{\operatorname{cov}}_{\textnormal{f,1}}(D,{C})-\hat{\operatorname{cov}}_{\textnormal{f,0}}(D,{C})\big{\\}}\\{\operatorname{cov}_{\textnormal{f}}({C})\\}^{-1}\big{\\{}\hat{\operatorname{cov}}_{\textnormal{f,1}}({C},D)-\hat{\operatorname{cov}}_{\textnormal{f,0}}({C},D)\big{\\}},$
$\displaystyle\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{G}}\\}$
$\displaystyle=\\{\hat{\operatorname{cov}}_{\textnormal{f,1}}(\tilde{U},\tilde{{G}})-\hat{\operatorname{cov}}_{\textnormal{f,0}}(\tilde{U},\tilde{{G}})\\}\\{\operatorname{cov}_{\textnormal{f}}(\tilde{{G}})\\}^{-1}\\{\hat{\operatorname{cov}}_{\textnormal{f,1}}(\tilde{{G}},\tilde{U})-\hat{\operatorname{cov}}_{\textnormal{f,0}}(\tilde{{G}},\tilde{U})\\}.$

Second, we estimate $R^{\textnormal{adj}}_{x}$ and $R^{\textnormal{adj}}_{c}$.
Recall that $R^{\textnormal{adj}}_{x}$ and $R^{\textnormal{adj}}_{c}$ are
defined similarly as $R_{x}$ and $R_{c}$ except that we replace $Y_{ij}(z)$ by
$Y_{ij}(z)-{w}_{ij}^{{\mathrm{\scriptscriptstyle T}}}{\beta}_{w}(z)$ and
$\tilde{Y}_{i\cdot}(z)$ by
$\tilde{Y}^{\textnormal{adj}}_{i\cdot}(z)=\tilde{Y}_{i\cdot}(z)-{v}_{i}^{{\mathrm{\scriptscriptstyle
T}}}{\beta}_{v}(z)$, for $z=0,1$. Applying Proposition 1 of Li et al., (2018),
we obtain the following results for $R^{\textnormal{adj}}_{x}$ and
$R^{\textnormal{adj}}_{c}$.

######  Proposition 3.

Under regularity conditions, $(R^{\textnormal{adj}}_{x})^{2}$ or
$(R^{\textnormal{adj}}_{c})^{2}$ can be expressed in terms of the variances of
the adjusted potential outcomes and of their projections on $\tilde{{X}}$ or
$C$:

$(R^{\textnormal{adj}}_{x})^{2}=\frac{e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)\mid\tilde{{X}}\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}}{e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(0)\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\\}},$
(4)

$(R^{\textnormal{adj}}_{c})^{2}=\frac{e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)\mid{C}\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}}{e_{1}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(0)\\}-\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\\}}.$
(5)

The denominator of (4) is equal to
$V_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$ and the
denominator of (5) is equal to
$V_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}$. Following
the results of Li et al., (2018), we can estimate
$\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(z)\mid{C}\\}$
by

$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(z)\mid{C}\\}=\hat{\operatorname{cov}}_{\textnormal{f,z}}(D,{C})\\{\hat{\operatorname{cov}}_{\textnormal{f,z}}({C})\\}^{-1}\hat{\operatorname{cov}}_{\textnormal{f,z}}({C},D).$

Similarly,
$\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}$
can be estimated by
$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}$.
Then we get an estimator for the numerator of (5). Thus, we estimate
$(R^{\textnormal{adj}}_{c})^{2}$ by

$(\hat{R}^{\textnormal{adj}}_{c})^{2}=\big{[}e_{1}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)\mid{C}\\}+e_{0}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}-\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{Y}^{\textnormal{adj}}(1)-\tilde{Y}^{\textnormal{adj}}(0)\mid{C}\\}\big{]}\Big{/}{\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}}.$

For cluster rerandomization using individual-level covariates,
$\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(z)\mid\tilde{{X}}\\}$
can be estimated by

$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(z)\mid\tilde{{X}}\\}=\hat{\operatorname{cov}}_{\textnormal{f,z}}(\tilde{U},\tilde{{X}})\\{\hat{\operatorname{cov}}_{\textnormal{f,z}}(\tilde{{X}})\\}^{-1}\hat{\operatorname{cov}}_{\textnormal{f,z}}(\tilde{{X}},{\tilde{U}}).$

Similarly,
$\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}$
can be estimated by

$\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}=\\{\hat{\operatorname{cov}}_{\textnormal{f,1}}(\tilde{U},\tilde{{X}})-\hat{\operatorname{cov}}_{\textnormal{f,0}}(\tilde{U},\tilde{{X}})\\}\\{\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})\\}^{-1}\\{\hat{\operatorname{cov}}_{\textnormal{f,1}}(\tilde{{X}},\tilde{U})-\hat{\operatorname{cov}}_{\textnormal{f,0}}(\tilde{{X}},\tilde{U})\\}.$

Again we derive a consistent estimator for the numerator of (4). Thus, we
estimate $(R^{\textnormal{adj}}_{x})^{2}$ by

$(\hat{R}^{\textnormal{adj}}_{x})^{2}=\big{[}{e_{1}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)\mid\tilde{{X}}\\}+e_{0}^{-1}\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}-\hat{\operatorname{var}}_{\textnormal{f}}\\{\tilde{\varepsilon}^{\textnormal{adj}}(1)-\tilde{\varepsilon}^{\textnormal{adj}}(0)\mid\tilde{{X}}\\}}\big{]}\Big{/}{\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}}.$

Finally, we estimate $\mu_{x}$ and $\mu_{c}$ defined in Theorem 7. Recall that

$\displaystyle{\mu}_{x}^{\mathrm{\scriptscriptstyle T}}$ $\displaystyle=$
$\displaystyle V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2}/\\{V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},x\tau}\\}^{1/2},$
(6) $\displaystyle{\mu}_{c}^{\mathrm{\scriptscriptstyle T}}$ $\displaystyle=$
$\displaystyle V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1/2}/\\{V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau
c}V_{\textnormal{ht},cc}^{-1}V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},c\tau}\\}^{1/2}.$
(7)

The denominators of (6) and (7) equal the square roots of the numerators of
(4) and (5), respectively. We have derived their consistent estimators. Recall
that

$\displaystyle
V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},x\tau}$
$\displaystyle=(V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{haj},\tau
x})^{{\mathrm{\scriptscriptstyle
T}}}=e_{1}^{-1}\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}^{\textnormal{adj}}(0)\\},$
$\displaystyle V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},c\tau}$
$\displaystyle=(V^{\scalebox{0.6}{$\textnormal{adj}$}}_{\textnormal{ht},\tau
c})^{{\mathrm{\scriptscriptstyle
T}}}=e_{1}^{-1}\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}^{\textnormal{adj}}(1)\\}+e_{0}^{-1}\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}^{\textnormal{adj}}(0)\\}.$

We estimate
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}^{\textnormal{adj}}(z)\\}$
by $\hat{\operatorname{cov}}_{\textnormal{f,z}}(\tilde{X},\tilde{U})$ and
$\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}^{\textnormal{adj}}(z)\\}$
by $\hat{\operatorname{cov}}_{\textnormal{f,z}}(C,D)$. Plugging in these
terms, we derive consistent estimators of $\mu_{x}$ and $\mu_{c}$, denoted as
$\hat{\mu}_{x}$ and $\hat{\mu}_{c}$, respectively.

Recall that $\xi_{k}$ is a $K$-dimensional vector with $1$ at the $k$th
dimension and $0$ at other dimensions and $I_{K}$ is a $K\times K$ identity
matrix. Let $q_{\zeta}(R^{2},\mu,A,a)$ be the $\zeta$th quantile of
$(1-R^{2})\epsilon+R^{2}\mu^{\mathrm{\scriptscriptstyle
T}}\eta\mid\eta^{\mathrm{\scriptscriptstyle T}}A\eta\leq a$. Theorem 9 below
provides asymptotically conservative confidence intervals.

######  Theorem 9.

Under regularity conditions,

$\displaystyle\Bigl{[}{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+\bigl{(}\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{x})^{2},\xi_{1},I_{K},a\bigr{\\}},\
\ \
{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+\bigl{(}\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{1-\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{x})^{2},\xi_{1},I_{K},a\bigr{\\}}\Bigr{]},$
$\displaystyle\Bigl{[}\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+\bigl{(}\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{c})^{2},\xi_{1},I_{K},a\bigr{\\}},\
\ \
\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+\bigl{(}\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{1-\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{c})^{2},\xi_{1},I_{K},,a\bigr{\\}}\Bigr{]}$

are asymptotically conservative $1-\varsigma$ confidence intervals of $\tau$
under $\mathcal{M}_{x}$ and $\mathcal{M}_{c}$, respectively;

$\displaystyle\Bigl{[}{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+\bigl{(}\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{x})^{2},\hat{\mu}_{x},V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2},a\bigr{\\}},\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad{\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}+\bigl{(}\hat{V}_{\textnormal{haj},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{1-\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{x})^{2},\hat{\mu}_{x},V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2},a\bigr{\\}}\Bigr{]},$
$\displaystyle\Bigl{[}\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+\bigl{(}\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{c}\\}^{2},\hat{\mu}_{c},V_{\textnormal{ht},cc}^{1/2}A_{c}V_{\textnormal{ht},cc}^{1/2},a\bigr{\\}},\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}+\bigl{(}\hat{V}_{\textnormal{ht},\tau\tau}^{\scalebox{0.6}{$\textnormal{adj}$}}\bigm{/}M\bigr{)}^{1/2}q_{1-\varsigma/2}\bigl{\\{}(\hat{R}^{\textnormal{adj}}_{c})^{2},\hat{\mu}_{c},V_{\textnormal{ht},cc}^{1/2}A_{c}V_{\textnormal{ht},cc}^{1/2},a\bigr{\\}}\Bigr{]}$

are asymptotically conservative $1-\varsigma$ confidence intervals of $\tau$
under $\mathcal{D}_{x}(A_{x})$ and $\mathcal{D}_{c}(A_{c})$, respectively.

##  Appendix D Results for cluster rerandomization with cluster-level
covariates that are omitted in Section 4.2

We assume that the $c_{i}$’s are orthogonalized by Gram-Schmidt
orthogonalization. Let

$V_{\textnormal{ht},cc}=\operatorname{diag}(V_{\textnormal{ht},c_{1}c_{1}},\ldots,V_{\textnormal{ht},c_{K}c_{K}}),\quad
V_{\textnormal{ht},\tau c}=(V_{\textnormal{ht},\tau
c_{1}},\ldots,V_{\textnormal{ht},\tau c_{K}}).$

Define the squared multiple correlation for the $k$th covariate as
$R_{c_{k}}^{2}=V_{\textnormal{ht},\tau
c_{k}}^{2}/(V_{\textnormal{ht},\tau\tau}V_{\textnormal{ht},c_{k}c_{k}}).$ The
optimal cluster rerandomization scheme is
$\mathcal{D}_{c}(A_{c}^{\textnormal{opt}})=\\{M{\hat{\tau}_{\textnormal{ht},c}}^{\mathrm{\scriptscriptstyle
T}}A_{c}^{\textnormal{opt}}{\hat{\tau}_{\textnormal{ht},c}}<a\\}$ with the
optimal diagonal matrix

$A_{c}^{\textnormal{opt}}=\operatorname{diag}\\{(V_{\textnormal{ht},\tau
c_{1}}V_{\textnormal{ht},c_{1}c_{1}}^{-1})^{2},\ldots,(V_{\textnormal{haj},\tau
c_{k}}V_{\textnormal{haj},c_{K}c_{K}}^{-1})^{2}\\}.$

Corollary 7 below is the counterpart of Corollary 2 in the main text. Recall
that

$p_{K}=\frac{2\pi}{K+2}\left\\{\frac{2\pi^{K/2}}{K\Gamma(K/2)}\right\\}^{-2/K}.$

######  Corollary 7.

Under regularity conditions with orthogonalized covariates and optimal
weighted Euclidean distance, if ${V}_{\textnormal{haj},\tau
c_{k}}{V}_{\textnormal{ht},c_{k}c_{k}}^{-1}$’s are nonzero for all
$k=1,\ldots,K$, then

$\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau)\mid\mathcal{D}_{c}(A_{c}^{\textnormal{opt}})\\}=V_{\textnormal{ht},\tau\tau}\Big{\\{}(1-R_{c}^{2})+K\Big{(}\prod_{k=1}^{K}R_{c_{k}}^{2}\Big{)}^{1/K}p_{K}\alpha^{2/K}+o(\alpha^{2/K})\Big{\\}},$

for a small $\alpha$.

Similar to Section 4.2, we partition the covariates into $L$ tiers with
$K_{l}$ covariates in tier $l$, for $l=1,\ldots,L$ and
$\sum_{l=1}^{L}K_{l}=K$. Let ${c}_{[l]}$ be the covariates in the $l$th tier.
Then ${c}=({c}_{[1]},\ldots,{c}_{[L]})$. Within tier $l$, define
$\hat{\tau}_{\textnormal{ht},c_{[l]}}$ as the Horvitz–Thompson estimator for
the covariates, $V_{\textnormal{ht},c_{[l]}c_{[l]}}$ as its covariance,
$\mathcal{M}_{[l]}$ as the event that the corresponding Mahalanobis distance
is smaller than or equal to the threshold $a_{[l]}$, $R_{c_{[l]}}^{2}$ as the
corresponding $R^{2}$, and $\alpha_{[l]}$ as the asymptotic acceptance rate.
The overall asymptotic acceptance rate satisfies
$\alpha=\prod_{l=1}^{L}\alpha_{[l]}$ due to the orthogonalization of the
covariates. Corollary 8 below is the counterpart of Corollary 3 for cluster
rerandomization with cluster-level covariates.

######  Corollary 8.

Under regularity conditions,

$\textnormal{var}_{\textnormal{a}}\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau)\mid\mathcal{M}_{[1]},\ldots,\mathcal{M}_{[L]}\\}=V_{\textnormal{ht},\tau\tau}\Big{\\{}(1-R_{c}^{2})+\sum_{l=1}^{L}R_{c_{[l]}}^{2}p_{K_{l}}\alpha_{[l]}^{2/K_{l}}+o(\alpha_{[l]}^{2/K_{l}})\Big{\\}},$

for small $\alpha_{[l]}$, $l=1,\ldots,L$.

Given $\alpha=\prod_{l=1}^{L}\alpha_{[l]}>0$, to minimize the second term in
the asymptotic variance in Corollary 8, we must choose

$\alpha_{[l]}=\left({c_{0}{R}_{c_{[l]}}^{2}p_{K_{l}}}/{K_{l}}\right)^{-K_{l}/2}\quad(l=1,\ldots,L),$

for some positive constant $c_{0}$. Theorem 10 below is the counterpart of
Theorem 5 for cluster rerandomization with cluster-level covariates.

######  Theorem 10.

Under regularity conditions with orthogonalized covariates, the following
inequality holds for any $\alpha_{[1]},\ldots,\alpha_{[L]}$ satisfying
$\alpha=\prod_{l=1}^{L}\alpha_{[l]}$:

$\displaystyle\sum_{l=1}^{L}{R}_{c_{[l]}}^{2}p_{K_{l}}\alpha_{[l]}^{2/K_{l}}\geq
K\Big{(}\prod_{k=1}^{K}R_{c_{k}}^{2}\Big{)}^{1/K}p_{K}\alpha^{2/K}.$

Theorem 10 quantifies the superiority of cluster rerandomization scheme based
on the optimal weighted Euclidean distance when the cluster-level covariates
are orthogonalized.

##  Appendix E A counterexample for Remark 1

Naturally, we want to compare the optimal $A_{x}$ or $A_{c}$ with the one used
in the Mahalanobis distance. Unfortunately, there could not be a universal
result since $\nu_{x}(A_{x}^{\textnormal{opt}})$ or
$\nu_{c}(A_{c}^{\textnormal{opt}})$ could be smaller than, equal to, or larger
than $\nu_{x}(V_{\textnormal{haj},xx}^{-1})=1$ or
$\nu_{c}(V^{-1}_{\textnormal{ht},cc})=1$, respectively, as shown in the
example below.

######  Example 1.

Suppose that

$V_{\textnormal{haj},xx}=\begin{bmatrix}4&\ \ \delta\\\ \delta&\ \
4\end{bmatrix},\quad V_{\textnormal{haj},\tau x}=(1,1).$

Then

$\nu_{x}(A_{x}^{\textnormal{opt}})/\nu_{x}(V_{\textnormal{haj},xx}^{-1})=\left(\frac{4-\delta}{4+\delta}\right)^{1/2}.$

When $\delta>0$, i.e., there is a positive correlation between covariates,
$\nu_{x}(A_{x}^{\textnormal{opt}})$ is smaller than
$\nu_{x}(V_{\textnormal{haj},xx}^{-1})$. When $\delta<0$, the opposite result
holds.

##  Appendix F Proofs

###  F.1 Regularity conditions

This section provides the regularity conditions for the theoretical results.
These conditions are required for using the finite population central limit
theorem. Let $||\cdot||_{\infty}$ denote the $\ell_{\infty}$ norm of a vector.

The regularity conditions below are used for deriving the asymptotic
properties of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$ under classic
cluster randomization.

######  Condition 1.

As $M\rightarrow\infty$, $e_{1}$ has a limit in $(0,1)$.

######  Condition 2.

As $M\rightarrow\infty$,

  1. (i)

the finite population variance
$\operatorname{var}_{\textnormal{f}}\\{\tilde{Y}(z)\\}$ $(z=0,1)$ and
covariance
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{Y}(0),\tilde{Y}(1)\\}$ have
finite limiting values, and the limit of $V_{\textnormal{ht},\tau\tau}$ is
positive;

  2. (ii)

$\max_{z=0,1}\max_{1\leq i\leq
M}\\{\tilde{Y}_{i\cdot}(z)-\bar{Y}(z)\\}^{2}/M\rightarrow 0$.

Condition 1 indicates that we have enough clusters assigned to both the
treatment and control arms. Condition 2 restricts the moments of the potential
outcomes.

The asymptotic properties of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$
under classic cluster randomization require Condition 1 and Condition 3 below.

######  Condition 3.

As $M\rightarrow\infty$,

  1. (i)

the finite population variance,
$\operatorname{var}_{\textnormal{f}}\\{\tilde{\varepsilon}(z)\\}$ $(z=0,1)$,
and covariances,
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{\varepsilon}(0),\tilde{\varepsilon}(1)\\}$
and
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{\varepsilon}(z),\tilde{\omega}\\}$,
have finite limiting values, and the limit of $V_{\textnormal{haj},\tau\tau}$
is positive;

  2. (ii)

$M^{-1}\max_{1\leq i\leq M}(\tilde{\omega}_{i}-1)^{2}\to 0$ and
$M^{-1}\max_{z=0,1}\max_{1\leq i\leq
M}\big{\\{}\tilde{\varepsilon}_{i\cdot}(z)\big{\\}}^{2}\to 0$.

The asymptotic properties of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}$
in Theorems 1–4 and 10, Corollaries 1 and 7–8, and Propositions 1–2 require
Conditions 1–2 and Condition 4 below.

######  Condition 4.

As $M\rightarrow\infty$,

  1. (i)

the finite population covariances, $\operatorname{cov}_{\textnormal{f}}({C})$
and $\operatorname{cov}_{\textnormal{f}}\\{{C},\tilde{Y}(z)\\}$ $(z=0,1)$, and
$A_{c}$ have finite limiting values, and the limits of
$\operatorname{cov}_{\textnormal{f}}({C})$ and $A_{c}$ are nonsingular;

  2. (ii)

$M^{-1}\max_{1\leq i\leq M}\|{{c}}_{i}\|_{\infty}^{2}\rightarrow 0$.

The asymptotic properties of $\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}$
in Theorems 1–5, Corollaries 1–3, and Propositions 1–2 require Conditions 1,
3, and Condition 5 below.

######  Condition 5.

As $M\rightarrow\infty$,

  1. (i)

the finite population covariances,
$\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})$,
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{{X}},\tilde{\varepsilon}(z)\\}$
$(z=0,1)$, and
$\operatorname{cov}_{\textnormal{f}}(\tilde{{X}},\tilde{\omega})$, and $A_{x}$
have finite limiting values, and the limits of
$\operatorname{cov}_{\textnormal{f}}(\tilde{{X}})$ and $A_{x}$ are
nonsingular;

  2. (ii)

$M^{-1}\max_{1\leq i\leq M}\|\tilde{{x}}_{i\cdot}\|_{\infty}^{2}\rightarrow
0$.

The asymptotic properties of
$\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}$
in Theorems 6 and 7–9, Corollaries 4–6, and Proposition 3 require Conditions
1, 2, 4, and Condition 6 below.

######  Condition 6.

As $M\rightarrow\infty$,

  1. (i)

$M^{-1}\max_{1\leq i\leq M}\left\|{v}_{i}\right\|^{2}_{\infty}=o(1)$;

  2. (ii)

the finite population covariances, $\operatorname{cov}_{\textnormal{f}}({V})$,
$\operatorname{cov}_{\textnormal{f}}\\{{V},\tilde{Y}(z)\\}$ ($z=0,1$), and
$\operatorname{cov}_{\textnormal{f}}({V},{C})$, have finite limiting values,
and the limit of $\operatorname{cov}_{\textnormal{f}}({V})$ is nonsingular;

The asymptotic properties of
${\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}^{\scalebox{0.6}{$\textnormal{adj}$}}}$
in Theorems 6 and 7–9 and Proposition 3 require Conditions 1, 3, 5, and
Condition 7 below.

######  Condition 7.

As $M\rightarrow\infty$,

  1. (i)

$N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}{w}_{ij}\varepsilon_{ij}(z)$ ($z=0,1$)
converges to a finite vector, and
$N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}{w}_{ij}{w}_{ij}^{\mathrm{T}}$
converges to a finite and invertible matrix;

  2. (ii)

$\max_{1\leq i\leq M}\tilde{\omega}_{i}=o(M^{1/3})$,
$N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}\left\|{w}_{ij}\right\|_{\infty}^{4}=O(1)$,
$M^{-1}\max_{1\leq i\leq
M}\left\|{\tilde{{w}}}_{i\cdot}\right\|^{2}_{\infty}=o(1)$,
$N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{n_{i}}\\{Y_{ij}(z)\\}^{4}=O(1)\ (z=0,1)$;

  3. (iii)

the finite population covariances,
$\operatorname{cov}_{\textnormal{f}}(\tilde{{W}})$,
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{{W}},\tilde{{\varepsilon}}(z)\\}$
($z=0,1$),
$\operatorname{cov}_{\textnormal{f}}\\{\tilde{{W}},\tilde{\omega}\\}$, and
$\operatorname{cov}_{\textnormal{f}}(\tilde{{W}},\tilde{{X}})$, have finite
limiting values, and the limit of
$\operatorname{cov}_{\textnormal{f}}(\tilde{{W}})$ is nonsingular.

###  F.2 Some preliminary results on weak convergence

Before proving the main results, we introduce some preliminary results on weak
convergence derived by Li et al., (2018). Consider the combination
$(\hat{\tau},\hat{\tau}_{\star})$ as the estimator of the average treatment
effects for the outcomes and covariates, for example,
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}},{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}})$ and
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}},{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}})$.
Consider a balance criterion $\phi(M^{1/2}\hat{\tau}_{\star},A)$ which is a
binary function, and $A$ is the matrix used to conduct cluster
rerandomization. Denote the event that a treatment assignment is acceptable
under the balance criterion $\phi(M^{1/2}\hat{\tau}_{\star},A)$ as
$\phi(M^{1/2}\hat{\tau}_{\star},A)=1$, and let
$\mathfrak{G}=\left\\{{\mu}:\phi\left({\mu},A\right)=1\right\\}$ be the
acceptance region for $M^{1/2}\hat{\tau}_{\star}$. We write $A>0$ if $A$ is
strictly positive definite. Suppose that Condition 8 below on
$\phi(\cdot,\cdot)$ holds.

######  Condition 8.

  1. (i)

$\phi(\cdot,\cdot)$ is almost surely continuous;

  2. (ii)

if ${T}_{\star}\sim\mathcal{N}\left(0,V_{\star\star}\right)$, then
$\textnormal{pr}\left\\{\phi\left({T}_{\star},A\right)=1\right\\}>0$ for any
$A>0$, and
$\operatorname{var}\left\\{{T}_{\star}\mid\phi\left({T}_{\star},A\right)=1\right\\}$
is a continuous function of $V_{\star\star};$

  3. (iii)

$\phi\left({\mu},A\right)=\phi\left(-{\mu},A\right)$, for all ${\mu}$ and
$A>0.$

In Condition 8, (i) and (ii) impose smoothness constraints on $\phi$ and
prevent the acceptance region from being a set of measure zero; (iii) is a
symmetry consideration to ensure that relabel the treatment and control units
does not change the balance criterion. Cluster rerandomization schemes based
on the Mahalanobis distances and weighted Euclidean distances satisfy
Condition 8. We denote the variance of
$M^{1/2}(\hat{\tau},\hat{\tau}_{\star})$ as $V$, and we use the same notation
to represent its limit:

${V}=\begin{bmatrix}V_{\tau\tau}&V_{\tau\star}\\\
V_{\star\tau}&V_{\star\star}\end{bmatrix}.$

Assume that $V$ is strictly positive definite. Let $A_{\infty}$ be the limit
of $A$, and
$\mathfrak{G}_{\infty}=\left\\{{\mu}:\phi\left({\mu},{A}_{\infty}\right)=1\right\\}$
be the limit of $\mathfrak{G}$. Proposition 4 below is a direct result of
Proposition A1 in Li et al., (2018).

######  Proposition 4.

Under Condition 8, as $M\rightarrow\infty$, if

$M^{1/2}\left(\begin{array}[]{c}\hat{\tau}-\tau\\\
\hat{\tau}_{\star}\end{array}\right)\mathrel{\ \dot{\sim}\ }\mathcal{N}(0,V),$

then

$\displaystyle M^{1/2}\left(\begin{array}[]{c}\hat{\tau}-\tau\\\
\hat{\tau}_{\star}\end{array}\right)\biggm{|}M^{1/2}\hat{\tau}_{\star}\in\mathfrak{G}\mathrel{\
\dot{\sim}\ }\left(\begin{array}[]{c}{T}\\\
{{T}}_{\star}\end{array}\right)\biggm{|}{{T}}_{\star}\in\mathfrak{G},$

where $({T},{{T}_{\star}^{\mathrm{\scriptscriptstyle
T}}})^{{\mathrm{\scriptscriptstyle T}}}\sim\mathcal{N}({0},{V})$.

###  F.3 Proof of Proposition 1

We need the following Lemma from Su and Ding, (2021).

######  Lemma 1.

Under Conditions 1 and 3,

$\displaystyle
M^{1/2}\Bigl{[}\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\Big{\\{}M_{1}^{-1}\sum_{i=1}^{N}Z_{i}\tilde{\varepsilon}_{i\cdot}-M_{0}^{-1}\sum_{i=1}^{N}(1-Z_{i})\tilde{\varepsilon}_{i\cdot}\Big{\\}}\Bigr{]}=o_{\textup{p}}(1).$

Su and Ding, (2021) proved Lemma 1 under slightly stronger conditions. After
carefully examining their proof, we found that it can be easily generalized to
the case of our weaker conditions, so we omit the proof of Lemma 1.

###### of Proposition 1.

Applying the vector-form finite population central limit theorem developed by
Li and Ding, (2017) to the scaled cluster-total potential outcomes
$\tilde{Y}_{i\cdot}(z)$, $z=0,1$, and cluster-level covariates ${c}_{i}$, we
have, under Conditions 1, 2, and 4,

$M^{1/2}\left(\begin{array}[]{c}\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau\\\
{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$ }}}\end{array}\right)\mathrel{\
\dot{\sim}\
}\mathcal{N}\left(0,\begin{bmatrix}V_{\textnormal{ht},\tau\tau}&V_{\textnormal{ht},\tau
c}\\\ V_{\textnormal{ht},c\tau}&V_{\textnormal{ht},cc}\end{bmatrix}\right).$

We could get a similar result for the joint asymptotic distribution of
$(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}},{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{\mathrm{\scriptscriptstyle
T}})^{\mathrm{\scriptscriptstyle T}}$ after replacing
$\tilde{Y}_{i\cdot}(z)-\bar{Y}(z)$ by
$\tilde{\varepsilon}_{i\cdot}(z)=\tilde{Y}_{i\cdot}(z)-\tilde{\omega}_{i}\bar{Y}(z)$
and applying Lemma 1. Specifically, applying Lemma 1 to each element of
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau,{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{\mathrm{\scriptscriptstyle
T}})^{{\mathrm{\scriptscriptstyle T}}}$, we have, under Conditions 1 and 3,

$M^{1/2}\left(\begin{array}[]{c}\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau\\\
{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\end{array}\right)=M^{1/2}\Bigg{(}\begin{gathered}\scalebox{1.0}{$M_{1}^{-1}\sum\nolimits_{i=1}^{M}Z_{i}\tilde{\varepsilon}_{i\cdot}(1)-M_{0}^{-1}\sum\nolimits_{i=1}^{M}(1-Z_{i})\tilde{\varepsilon}_{i\cdot}(0)$}\\\
\scalebox{1.0}{$M_{1}^{-1}\sum\nolimits_{i=1}^{M}Z_{i}{\tilde{x}}_{i\cdot}-M_{0}^{-1}\sum\nolimits_{i=1}^{M}(1-Z_{i}){\tilde{x}}_{i\cdot}$}\end{gathered}\Bigg{)}+o_{\textup{p}}(1).$

Under Conditions 1, 3, and 5, applying the vector-form finite population
central limit theorem to $\tilde{\varepsilon}_{i\cdot}(z)$, $z=0,1$, and
${\tilde{x}}_{i\cdot}$, we have

$M^{1/2}\Bigg{(}\begin{gathered}\scalebox{1.0}{$M_{1}^{-1}\sum\nolimits_{i=1}^{M}Z_{i}\tilde{\varepsilon}_{i\cdot}(1)-M_{0}^{-1}\sum_{i=1}^{M}(1-Z_{i})\tilde{\varepsilon}_{i\cdot}(0)$}\\\
\scalebox{1.0}{$M_{1}^{-1}\sum_{i=1}^{M}Z_{i}{\tilde{x}}_{i\cdot}-M_{0}^{-1}\sum_{i=1}^{M}(1-Z_{i}){\tilde{x}}_{i\cdot}$}\end{gathered}\Bigg{)}\mathrel{\
\dot{\sim}\ }\mathcal{N}\left(0,\Bigg{[}\begin{aligned}
V_{\textnormal{haj},\tau\tau}\ &V_{\textnormal{haj},\tau x}\\\
V_{\textnormal{haj},x\tau}\
&V_{\textnormal{haj},xx}\end{aligned}\Bigg{]}\right).$

By Slutsky’s theorem, the conclusion holds. ∎

###  F.4 Proof of Theorem 1

Proposition 5 below follows from Propositions 1 and 4.

######  Proposition 5.

Under Conditions 1, 3, and 5,

$\displaystyle
M^{1/2}\begin{pmatrix}\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau\\\
\hat{\tau}_{\textnormal{haj},x}\end{pmatrix}\biggm{|}M^{1/2}\hat{\tau}_{\textnormal{haj},x}\in\mathfrak{G}_{x}\mathrel{\
\dot{\sim}\ }\begin{pmatrix}{T}\\\
T_{x}\end{pmatrix}\biggm{|}T_{x}\in\mathfrak{G}_{x},$

and under Conditions 1, 2, and 4,

$\displaystyle
M^{1/2}\left(\begin{array}[]{c}\hat{\tau}_{\textnormal{ht}}-\tau\\\
\hat{\tau}_{\textnormal{ht},c}\end{array}\right)\biggm{|}M^{1/2}\hat{\tau}_{\textnormal{ht},c}\in\mathfrak{G}_{c}\mathrel{\
\dot{\sim}\ }\left(\begin{array}[]{c}{T}\\\
{{T}}_{c}\end{array}\right)\biggm{|}{{T}}_{c}\in\mathfrak{G}_{c},$

where $({T},{{T}_{x}^{\mathrm{\scriptscriptstyle
T}}})^{{\mathrm{\scriptscriptstyle
T}}}\sim\mathcal{N}({0},{V}_{\textnormal{haj}})$ and
$({T},{{T}_{c}^{\mathrm{\scriptscriptstyle T}}})^{{\mathrm{\scriptscriptstyle
T}}}\sim\mathcal{N}({0},{V}_{\textnormal{ht}})$ with

${V}_{\textnormal{haj}}=\Bigg{[}\begin{aligned} V_{\textnormal{haj},\tau\tau}\
&V_{\textnormal{haj},\tau x}\\\ V_{\textnormal{haj},x\tau}\
&V_{\textnormal{haj},xx}\end{aligned}\Bigg{]},\quad{V}_{\textnormal{ht}}=\begin{bmatrix}V_{\textnormal{ht},\tau\tau}&V_{\textnormal{ht},\tau
c}\\\ V_{\textnormal{ht},c\tau}&V_{\textnormal{ht},cc}\end{bmatrix}.$

###### of Theorem 1.

The asymptotic distribution of
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau)\mid\mathcal{M}_{c}$
is a direct result of applying Theorem 1 of Li et al., (2018) to
$\tilde{Y}_{i\cdot}(z)$, $z=0,1$, and ${c}_{i}$.

To derive the asymptotic distribution of
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{M}_{x}$,
similar to the proof of Proposition 1, we have
$\\{M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau),M^{1/2}{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}},A_{x}\\}$
has the same asymptotic distribution as $\\{T,{T}_{x},A_{x}\\}$ with
$A_{x}=V_{\textnormal{haj},xx}^{-1}$. Therefore,
$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid
M{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}\leq a$ has the same
asymptotic distribution as $T\mid{T}_{x}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}{T}_{x}\leq a$. Applying Theorem 1 of Li et al., (2018) to
$\tilde{\varepsilon}_{i\cdot}(z)$, $z=0,1$, and ${\tilde{x}}_{i\cdot}$, we
have

$T\mid{T}_{x}^{{\mathrm{\scriptscriptstyle T}}}A_{x}{T}_{x}\leq a\mathrel{\
\dot{\sim}\
}V_{\textnormal{ht},\tau\tau}^{1/2}\\{(1-R_{x}^{2})^{1/2}\epsilon+R_{x}L_{K,a}\\}.$

The conclusion follows immediately.

∎

###  F.5 Proof of Corollary 1

###### of Corollary 1.

For the cluster-level covariates ${c}_{i}=(n_{i},{\tilde{x}}_{i\cdot})$, we
have

$\displaystyle
V_{\textnormal{ht},\tau\tau}=(e_{1}e_{0})^{-1}\operatorname{var}_{\textnormal{f}}\big{\\{}e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)\big{\\}},$
$\displaystyle
V_{\textnormal{ht},\tau\tau}R^{2}_{c}=(e_{1}e_{0})^{-1}\operatorname{cov}_{\textnormal{f}}\big{\\{}e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1),{C}\big{\\}}\operatorname{cov}_{\textnormal{f}}({C})^{-1}\operatorname{cov}_{\textnormal{f}}\big{\\{}{C},e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)\big{\\}}.$

Therefore, $e_{1}e_{0}V_{\textnormal{ht},\tau\tau}(1-R_{c}^{2})$ is the
variance of the residuals when we regress
$e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)$ on $C$ with intercept, i.e.,

$e_{1}e_{0}V_{\textnormal{ht},\tau\tau}\left(1-R_{c}^{2}\right)=\min_{\beta\in\mathbb{R}^{K+1},\
\beta_{0}\in\mathbb{R}}\operatorname{var}_{\textnormal{f}}\left\\{e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)-{1}_{M}\beta_{0}-(N/M\tilde{\omega},\tilde{{X}})\beta\right\\},$

where ${1}_{M}$ is an $M$-dimensional all-one vector. Similarly, we have

$\displaystyle
e_{1}e_{0}V_{\textnormal{haj},\tau\tau}\left(1-R_{x}^{2}\right)$
$\displaystyle=$ $\displaystyle\min_{\beta\in\mathbb{R}^{K},\
\beta_{0}\in\mathbb{R}}\operatorname{var}_{\textnormal{f}}\left\\{e_{1}\tilde{\varepsilon}(0)+e_{0}\tilde{\varepsilon}(1)-{1}_{M}\beta_{0}-\tilde{{X}}\beta\right\\}$
$\displaystyle=$ $\displaystyle\min_{\beta\in\mathbb{R}^{K},\
\beta_{0}\in\mathbb{R}}\operatorname{var}_{\textnormal{f}}\left[e_{1}\tilde{Y}(0)+e_{0}\tilde{Y}(1)-\\{e_{1}\bar{Y}(0)+e_{0}\bar{Y}(1)\\}\tilde{\omega}-{1}_{M}\beta_{0}-\tilde{{X}}\beta\right],$

where the second equality follows from
$\tilde{\varepsilon}_{i\cdot}(z)=\tilde{Y}_{i\cdot}(z)-\tilde{\omega}_{i}\bar{Y}(z).$
Therefore,

$V_{\textnormal{ht},\tau\tau}\left(1-R_{c}^{2}\right)\leq
V_{\textnormal{haj},\tau\tau}\left(1-R_{x}^{2}\right),$

which completes the proof. ∎

###  F.6 Proof of Theorem 2

###### of Theorem 2.

By Proposition 1,

$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{ht}$}}-\tau,{\hat{\tau}_{\scalebox{0.6}{{ht},${c}$
}}}^{\mathrm{\scriptscriptstyle T}})\mathrel{\ \dot{\sim}\
}\left(T,{T}_{c}^{\mathrm{\scriptscriptstyle T}}\right),\quad
M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau,{\hat{\tau}_{\scalebox{0.6}{{haj,${x}$}}}}^{\mathrm{\scriptscriptstyle
T}})\mathrel{\ \dot{\sim}\ }\left(T,{T}_{x}^{\mathrm{\scriptscriptstyle
T}}\right),$

where $(T,{T}_{c}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle
T}}\sim\mathcal{N}(0,V_{\textnormal{ht}})$ and
$(T,{T}_{x}^{\mathrm{\scriptscriptstyle T}})^{\mathrm{\scriptscriptstyle
T}}\sim\mathcal{N}(0,V_{\textnormal{haj}})$. We only give the proof for
cluster rerandomization with individual-level covariates, as the proof for
cluster rerandomization with cluster-level covariates is similar. By
Proposition 5, we have

$\displaystyle
M^{1/2}(\hat{\tau}_{\textnormal{haj},x}-\tau)\mid\mathcal{D}_{x}(A_{x})$
$\displaystyle=$ $\displaystyle
M^{1/2}(\hat{\tau}_{\textnormal{haj},x}-\tau)\mid(M^{1/2}\hat{\tau}_{\textnormal{haj},x})^{{\mathrm{\scriptscriptstyle
T}}}A_{x}(M^{1/2}\hat{\tau}_{\textnormal{haj},x})\leq a$
$\displaystyle\mathrel{\ \dot{\sim}\ }$ $\displaystyle
T\mid{T}_{x}^{{\mathrm{\scriptscriptstyle T}}}A_{x}{T}_{x}\leq a.$

We decompose $T$ as

$T=T-V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x}+V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x}.$

Note that $T-V_{\textnormal{haj},\tau x}V_{\textnormal{haj},xx}^{-1}{T}_{x}$
is independent of ${T}_{x}$, and thus, independent of
$V_{\textnormal{haj},\tau x}V_{\textnormal{haj},xx}^{-1}{T}_{x}$ and
${T}_{x}^{{\mathrm{\scriptscriptstyle T}}}A_{x}{T}_{x}$. Therefore,

$T-V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x}\mid{T}_{x}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}{T}_{x}\leq a\sim T-V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x}.$

By the definition of $R_{x}^{2}$, we have $R_{x}^{2}=V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V_{\textnormal{haj},x\tau}/V_{\textnormal{haj},\tau\tau}$,
and therefore

$\operatorname{var}(T-V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x})=V_{\textnormal{haj},\tau\tau}-V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V_{\textnormal{haj},x\tau}=(1-R_{x}^{2})V_{\textnormal{haj},\tau\tau}.$

Thus,

$T-V_{\textnormal{haj},\tau x}V_{\textnormal{haj},xx}^{-1}{T}_{x}\sim
V_{\textnormal{haj},\tau\tau}^{1/2}(1-R_{x}^{2})^{1/2}\epsilon.$

Since $V_{\textnormal{haj},xx}^{-1/2}{T}_{x}$ is a $K$-dimensional standard
normal random vector,

$\displaystyle V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}{T}_{x}\mid{T}_{x}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}{T}_{x}\leq a$ $\displaystyle\sim$ $\displaystyle
V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2}(V_{\textnormal{haj},xx}^{-1/2}{T}_{x})\mid({T}_{x}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{-1/2})(V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2})(V_{\textnormal{haj},xx}^{-1/2}{T}_{x})\leq
a$ $\displaystyle\sim$ $\displaystyle V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2}{\eta}\leq
a$ $\displaystyle\sim$ $\displaystyle
V_{\textnormal{haj},\tau\tau}^{1/2}(R_{x}{\mu}_{x}^{{\mathrm{\scriptscriptstyle
T}}}{\eta})\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2}{\eta}\leq
a,$

where ${\eta}$ is a $K$-dimensional standard normal random vector and is
independent of $\epsilon$ and

${\mu}_{x}^{{\mathrm{\scriptscriptstyle T}}}=(V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1}V_{\textnormal{haj},x\tau})^{-1/2}V_{\textnormal{haj},\tau
x}V_{\textnormal{haj},xx}^{-1/2}.$

Therefore,

$V_{\textnormal{haj},\tau\tau}^{1/2}\big{\\{}(1-R_{x}^{2})^{1/2}\epsilon+R_{x}{\mu}_{x}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}{V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2}}_{x}{\eta}\leq
a\big{\\}}\mathrel{\ \dot{\sim}\ }T\mid{T}_{x}^{{\mathrm{\scriptscriptstyle
T}}}A_{x}{T}_{x}\leq a,$

which implies that

$M^{1/2}(\hat{\tau}_{\scalebox{0.6}{$\textnormal{haj}$}}-\tau)\mid\mathcal{D}_{x}(A_{x})\mathrel{\
\dot{\sim}\
}V_{\textnormal{haj},\tau\tau}^{1/2}\big{\\{}(1-R_{x}^{2})^{1/2}\epsilon+R_{x}{\mu}_{x}^{{\mathrm{\scriptscriptstyle
T}}}{\eta}\mid{\eta}^{{\mathrm{\scriptscriptstyle
T}}}V_{\textnormal{haj},xx}^{1/2}A_{x}V_{\textnormal{haj},xx}^{1/2}{\eta}\leq
a\big{\\}}.$

∎

###  F.7 Proof of Proposition 2

First, we introduce some useful lemmas. Lemma 2 below is from Li et al.,
(2018, Lemma A8).

######  Lemma 2.

If $A$ and $B$ are independent, and are both symmetric around zero and
unimodal, then $A+B$ is symmetric around zero and unimodal.

Recall that the asymptotic distributions in Theorem 2 are composed of two
independent components: a normal random variable and a truncated normal random
variable. To simplify the presentation, we say that a random variable is
symmetric unimodal if it is symmetric around zero and unimodal. Clearly, the
normal random variable is symmetric unimodal. To prove the symmetric
unimodality of the truncated normal component, we follow the proof of Li et
al., (2020, Proposition 2) using properties of symmetric unimodality of random
vectors. The definition of symmetric unimodal random vector follows from
Dharmadhikari and Joag-Dev, (1988) as an extension of univariate case.

######  Definition 1.

For a set $\mathcal{B}$ of distributions on $\mathbb{R}^{K}$, we say that
$\mathcal{B}$ is closed convex if it satisfies two properties: (i) for any
distributions $\nu_{1},\nu_{2}\in\mathcal{B}$ and for any $\lambda\in(0,1)$,
the distribution $(1-\lambda)\nu_{1}+\lambda\nu_{2}\in\mathcal{B}$, and (ii) a
distribution $\nu\in\mathcal{B}$ if there exist a sequence of distributions in
$\mathcal{B}$ converging weakly to $\nu$.

For any set $\mathcal{C}$ of distributions, the closed convex hull of
$\mathcal{C}$ is the smallest closed convex set containing $\mathcal{C}$. A
compact convex set in Euclidean space $\mathbb{R}^{K}$ is called a convex body
if it has a nonempty interior. A set $\mathcal{K}\subset\mathbb{R}^{K}$ is
symmetric if $\mathcal{K}=\\{-a:a\in\mathcal{K}\\}$.

######  Definition 2.

A distribution on $\mathbb{R}^{K}$ is symmetric unimodal if it is in the
closed convex hull of $\mathcal{U}$, where $\mathcal{U}$ is the set of all
uniform distributions on symmetric convex bodies in $\mathbb{R}^{K}$.

Lemmas 3 and 4 below are from Li et al., (2020, Lemma A6) and Li et al.,
(2020, Lemma A8) respectively.

######  Lemma 3.

If $\psi\in\mathbb{R}^{K}$ is a symmetric unimodal random vector, then for any
non-random vector $b\in\mathbb{R}^{K}$, $b^{\mathrm{\scriptscriptstyle
T}}\psi$ is a symmetric unimodal random variable.

######  Lemma 4.

If a random vector in $\mathbb{R}^{K}$ has a log-concave density, then it is
symmetric unimodal.

Lemma 5 below is new which shows the symmetric unimodality of a truncated
normal random vector.

######  Lemma 5.

The truncated normal random variable $\rho b^{\mathrm{\scriptscriptstyle
T}}\eta\mid\eta^{\mathrm{\scriptscriptstyle T}}A\eta\leq a$ is symmetric