\section{Introduction}
The split-plot design has been widely used in agricultural sciences \citep{fisher1925statistical,Yates1937} and industrial experiments \citep{yates1935complex, jones}, and is gaining increasing popularity in social and biomedical sciences \citep{olken2007monitoring, meon, breza}.
The experimental units, also known as the {\it subplots}, are nested within groups known as the {\it whole plots}.
The split-plot design assigns the whole-plot intervention at the whole-plot level via a cluster randomization, and the subplot intervention at the subplot level via a stratified randomization.
By design, subplots within the same whole plot receive the same level of the whole-plot intervention. This provides a convenient way to accommodate hard-to-change factors and avoid interference within whole plots.

\cite{kempthorne1952design} initiated the discussion on design-based inference
of split-plot designs under the assumption of additive treatment effects.
\cite{zhao2018randomization} loosened the requirement on additivity, and established the theory for finite-sample exact inference in {\it uniform} split-plot designs, i.e., the whole-plot sizes and proportions of treated units for the subplot intervention within each whole plot are constant across whole plots.
\cite{mukerjee2019causal} extended the discussion to possibly nonuniform split-plot designs, and established  the finite-sample exact theory for the Horvitz--Thompson estimator.
\cite{zhao2021reconciling} extended the theory to the Hajek estimator and established the consistency and asymptotic normality of the Horvitz--Thompson and Hajek estimators in possibly nonuniform split-plot designs.

In split-plot designs, experimenters often collect baseline covariates  at both the whole-plot and subplot levels.
For example, in a split-plot design with students as subplots nested within whole plots of classes, class characteristics such as class size and teacher experience are whole-plot covariates, whereas student characteristics like race and gender are subplot covariates.
These baseline covariates are measured prior to the physical implementation of treatment assignments and hence not affected by the treatment. Randomization ensures that covariates are balanced across treatment levels on average.
However, covariate imbalance often exists in a particular treatment allocation, and can complicate the interpretation of the experimental results   \citep{rubin2008comment,Morgan2012,krieger2019bka}.
Rerandomization arose in such context and enforces covariate balance
in the design stage of randomized experiments \citep{Morgan2012}.
It has drawn much attention in the field of experimental design recently
and is shown to ensure efficiency gains in various settings \citep[see, e.g.,][]{moulton2004covariate,Morgan2015,Li2018,Li2020factorial,Wang2021,Zhu2021,zdfactx,zdrep,Lu2022}.

The existing literature of rerandomization focuses on treatments assigned at either the unit or the group level, but not both, leaving the corresponding theory for rerandomization in split-plot designs an open problem.
To fill this gap, we define {\it split-plot rerandomization} as a split-plot design compounded with rerandomization to balance covariates, and propose two split-plot rerandomization schemes based on the Mahalanobis distances of the Horvitz--Thompson and Hajek estimators of contrasts of covariate means to the origin, respectively.
We derive the asymptotic distributions of the Horvitz--Thompson and Hajek estimators for the average treatment effects under split-plot rerandomization and demonstrate the efficiency gains relative to split-plot randomization.

Regression adjustment
 is another approach to dealing with covariate imbalance, taking place in the analysis stage.
The existing literature sees efficiency gains by regression adjustment in various randomized experiments, including completely randomized experiments  \citep{lin2013,bloniarz2015lasso,Lei2020,zdfactx},
stratified randomized experiments \citep{Liu2019,Zhu2021block,Liu2022,Ma2022}, cluster randomized experiments \citep{Su2021,Lu2022},
completely or stratified randomized factorial experiments \citep{lu2016covariate,lu2016randomization,Liu2021,zdfact}, and split-plot experiments \citep{zhao2021reconciling}.
In particular, \cite{zhao2021reconciling} studied several specifications for regression adjustment in split-plot designs
and recommended an aggregate specification with full treatment-covariate interactions to ensure
efficiency gains when only whole-plot covariates are used.
Recent work by \cite{Li2020}, \cite{Wang2021}, and \cite{zdfactx, zdrep} further recommended combining rerandomization and regression adjustment in randomized experiments with treatments assigned at the unit level.

In this paper, we propose a novel alternative to regression adjustment for covariate adjustment in the analysis stage, and provide a design-based theory for the combination of rerandomization and covariate adjustment in split-plot designs.
We consider two  strategies for covariate adjustment in the analysis stage
and derive their asymptotic distributions under split-plot rerandomization.
The first strategy
follows the regression formulation by \cite{zhao2021reconciling}, and ensures efficiency gains when only whole-plot covariates are used.
Different from rerandomization with treatments assigned at the unit level, the asymptotic distributions of the regression-adjusted estimators under split-plot rerandomization are not normal, but convolutions of a normal component and a truncated normal component.
Moreover, the regression adjustment
may degrade efficiency if heterogeneous sub-plot covariates are used.
The second strategy is new, and approaches covariate adjustment from a projection or conditional inference perspective.
It adjusts an estimator for its asymptotic conditional bias given contrasts of covariate means, and yields estimators that are consistent and asymptotically normal under split-plot rerandomization with guaranteed efficiency gains.

We use the following notation.
Let $\mathcal{I}(\cdot)$ be the indicator function. Let $\chi_n^2$ denote the chi-squared distribution with $n$ degrees of freedom.
Let $0_m$ and $0_{m \times n} $ be the $m \times 1$ vector and $m \times n$ matrix of zeros, respectively. Let $1_m$ and $1_{m \times n} $ be the $m \times 1$ vector and $m \times n$ matrix of ones, respectively. Let $I_m$ be the $m \times m$ identity matrix. We suppress the dimensions when no confusion arises.
Let $\otimes$ and $\circ$ denote the Kronecker and Hadamard products of matrices, respectively. For two matrices $D_1$ and $D_2$, write $D_1 \geq D_2$ if $D_1 - D_2 $ is positive semi-definite.
Let  $||\cdot||_\infty$ denote the $\ell_\infty$ norm.
Let $\rightsquigarrow$ denote convergence in distribution. For a sequence of random variables $(U_n)_{n=1}^\infty$, write  $U_n \rightsquigarrow  U$ if as $n$ goes to $\infty$, the asymptotic distribution of $U_n$ equals the distribution of $U$.
Let $\textnormal{pr}_\textup{a}$ and $\textnormal{cov}_\textup{a}$ denote the asymptotic probability and covariance, respectively.

\section{Review of the split-plot design}

We follow the framework and notation in \cite{zhao2021reconciling}.
Consider a $2^2$ split-plot design with two binary factors of interest, indexed by A, B $\in \{0,1\}$.
This defines four treatment combinations, $\mathcal{T} = \{z = (a,b) :a,b=0,1\}$, where $a$ and $b$ index the levels of factors A and B, respectively.
We abbreviate $(a,b)$ as $(ab)$ when no confusion would arise.
Assume a study population of $N$ units nested in $W$ groups of possibly different sizes $M_w$ ($w = 1, \ldots,W$; $\sum_{w=1}^{W} M_w = N$).
We refer to each group as a whole plot and each unit as a subplot.
Index by $ws$ the $s$th subplot in whole plot $w$, and let $\mathcal{S} = \{ws: w=1,\ldots,W;\ s=1,\ldots,M_w\}$ denote the entire population.
The $2^2$ split-plot design assigns the units to different treatment combinations in two stages:
\begin{itemize}
    \item[(I)] the first stage assigns factor A  at the whole-plot level by a cluster randomization; that is,  it randomly assigns $W_a$ whole plots to receive level $a\in\{0,1\}$ of factor A for prespecified $W_a$'s with $W_0 + W_1 = W$;
    \item[(II)] the second stage assigns factor B at the subplot level by a stratified randomization; that is, it randomly assigns $M_{wb}$ units in whole plot $w$ to receive level $b\in\{0,1\}$ of factor B  for prespecified $M_{wb}$'s with $M_{w0} + M_{w1} = M_w, w=1,\ldots,M$, and the assignments across different whole plots are independent.
\end{itemize}
The final treatment of subplot $ws$, denoted by $Z_{ws} \in \mathcal T$, is then a combination of the level of factor A received by whole plot $w$ in stage (I) and the level of factor B received by itself in stage (II).
Refer to factor A and factor B as the whole-plot and subplot factors, respectively.
The probability of a whole plot assigned to level $a$ of factor A is $p_a=W_a/W$ for $a=0, 1$.
The probability of a subplot in whole plot $w$ assigned to level $b$ of factor B is $q_{wb}=M_{wb}/M_w$ for $w=1,\ldots,W$ and $b=0, 1$.
Assume that the cluster and stratified randomizations are independent throughout. The probability of subplot $ws$ assigned to treatment $z=(ab)$ is $p_{ws}(z)=p_a q_{wb}$.

Let $\bar M = N/W$ denote the average size of the whole plots, and let $\alphaw = M_w/\bar M$ denote the whole-plot size factor with  $W^{-1} \sum_{w=1}^{W}\alphaw = 1$. We call a split-plot design {\it uniform} if $M_w$ and $M_{wb}$ are constants across $w = 1, \ldots, W$. A uniform design has $\alphaw = 1$ for all $w$.

We define treatment effects using the potential outcomes framework \citep{Neyman:1923,Rubin:1974}. Denote by $Y_{ws}(z)$ the potential outcome of subplot $ws$ if assigned to treatment $z \in \mathcal {T}$, and let $\bar Y(z) = N^{-1} \sum_{ws \in \mathcal{S}} Y_{ws} (z)$ be the finite population average.
The main effects and interaction under the $2^2$ split-plot design are
\begin{align}
    \tau_{\textsc{a}} & = 2^{-1}\left\{\bar Y(10) + \bar Y(11) \right\}
    - 2^{-1}\left\{\bar Y(00) + \bar Y(01) \right\}, \notag \\
    \tau_{\textsc{b}} & = 2^{-1}\left\{\bar Y(01) + \bar Y(11) \right\}
    - 2^{-1}\left\{\bar Y(00) + \bar Y(10) \right\}, \notag \\
    \tau_{\textsc{ab}} & = \left\{\bar Y(00) + \bar Y(11) \right\}
    - \left\{\bar Y(01) + \bar Y(10) \right\} \notag
\label{Def tau}
\end{align}
\citep{mukerjee2019causal,zhao2021reconciling}.
Let $\bar{Y} = (\bar{Y}(00), \bar{Y}(01), \bar{Y}(10),\bar{Y}(11))^\text{T}$ vectorize the $\bar Y(z)$'s in lexicographical order of $z$.
We write the three effects in vector form as
\begin{eqnarray*}
\tau =(\tau_{\textsc{a}}, \tau_{\textsc{b}}, \tau_{\textsc{ab}})^\text{T}= G \bar Y
\end{eqnarray*} with  $G = (g_{\textsc{a}}, g_{\textsc{b}}, g_{\textsc{ab}})^\text{T}$ and $g_{\textsc{a}} = 2^{-1}(-1, -1, 1, 1)^\text{T}$, $g_{\textsc{b}} = 2^{-1}(-1, 1, -1, 1)^\text{T}$, $g_{\textsc{ab}} = (1, -1, -1, 1)^\text{T}$.
There are other effects of interest, $\tau_g = g^\text{T} \bar Y$, where $g$ is a $4\times 1$ contrast vector with $g^\text{T} 1_4 = 0$ \citep{Cuesta2019,zdfact}. Such a $g$ can be represented by a linear combination of $g_{\textsc{a}}$, $g_{\textsc{b}}$, and $g_{\textsc{ab}}$ such that $\tau_g$ is a linear transformation of $\tau$. To simplify the presentation, we focus on $\tau$ in this paper.

The observed outcome for subplot $ws$ is $Y_{ws} = \sum_{z \in \mathcal{T} } \mathcal{I}(Z_{ws}=z) Y_{ws}(z)$.
Let $\mathcal{S}(z) = \{ws: Z_{ws} = z, \ ws \in \mathcal{S} \}$ denote the set of subplots assigned to treatment $z \in \mathcal{T}$.
The Horvitz--Thompson estimator for $\bar{Y}(z)$ is
\begin{eqnarray*}
\hat Y_{\textnormal{ht}} (z) = N^{-1}\sum_{ws \in \mathcal{S}(z)}{p^{-1}_{ws}}(z){Y_{ws}} = N^{-1}\sum_{ws \in \mathcal{S}}\frac{\mathcal{I}(Z_{ws}=z)}{p_{ws}(z)}Y_{ws}(z),
\end{eqnarray*}
{and is unbiased under the $2^2$ split-plot randomization.}
Let $\hat Y_{\textnormal{ht}}$ be the vectorization of $\{ \hat Y_{\textnormal{ht}} (z) \}_{z \in \mathcal{T}}$ in lexicographical order of $z$.
We call $\hat \tau_{\HT} = G \hat Y_{\textnormal{ht}}$ the Horvitz--Thompson estimator of $\tau$, which is unbiased under the split-plot randomization.
A major drawback of the Horvitz--Thompson estimator is that it is not invariant to location shifts \citep{Fuller2009}. To address this issue, another widely used estimator, the Hajek estimator, is defined as
$$
 \hat Y_{\textnormal{haj}} (z)=\frac{\hat Y_{\textnormal{ht}}(z)}{\hat 1_{\textnormal{ht}}(z)},
$$
where $\hat 1_{\textnormal{ht}}(z) = N^{-1}\sum_{ws \in \mathcal{S}(z)}{p^{-1}_{ws}}(z)$ is the Horvitz--Thompson estimator of constant 1.
As pointed out by \cite{zhao2021reconciling}, the Hajek estimator is a ratio estimator for $\bar{Y}(z) = \bar{Y}(z) / 1$ with the numerator and denominator estimated by their Horvitz--Thompson estimators, respectively. Let $\hat Y_{\textnormal{haj}}$ be the vectorization of $\{ \hat Y_{\textnormal{haj}} (z) \}_{z \in \mathcal{T}}$ in lexicographical order of $z$. We call $\hat \tau_{\haj} = G \hat Y_{\textnormal{haj}}$ the Hajek estimator of $\tau$.

We adopt the design-based framework, which conditions on the potential outcomes and evaluates the sampling properties of $\hat \tau_{\HT}$ and $\hat \tau_{\haj}$ over the joint distribution of $Z_{ws}$'s.
Let $\bar{Y}_w(z) = M_w^{-1}\sum_{s=1}^{M_w} Y_{ws}(z)$ be the average potential outcome in whole plot $w$. The covariances of $\hat \tau_{\HT}$ and $\hat \tau_{\haj}$ under split-plot randomization depend on the {\it scaled between-} and {\it within-whole-plot covariances} of $\{Y_{ws}(z): ws\in\mathcal S; \, z\in\mathcal {T}\}$ defined as follows: $S_{\textnormal{ht}} = (S_{\textnormal{ht}}(z,z'))_{z,z'\in\mathcal {T}}$, $S_{\textnormal{haj}}= (S_{\textnormal{haj}}(z,z'))_{z,z'\in\mathcal {T}}$, and $S_w = (S_w(z,z'))_{z,z'\in\mathcal {T}}$ for $w = 1, \ldots, W$, where
\begin{eqnarray*}
S_{\textnormal{ht}}(z,z^\prime) &=& (W-1)^{-1}\sum_{w=1}^W\left\{\alpha_w \bar Y_w(z) - \bar Y(z) \right\} \left\{\alpha_w \bar Y_w(z^\prime) - \bar Y(z^\prime)\right\}, \\
S_{\textnormal{haj}}(z,z^\prime) &=& (W-1)^{-1}\sum_{w=1}^W\alpha_w^2\left\{ \bar Y_w(z) - \bar Y(z) \right\} \left\{\bar Y_w(z^\prime) - \bar Y(z^\prime)\right\},\\
S_w(z,z^\prime) &=& (M_w-1)^{-1}\sum_{s=1}^{M_w}\alpha_w^2\left\{  Y_{ws}(z) - \bar Y_w(z) \right\} \left\{ Y_{ws}(z^\prime) - \bar Y_w(z^\prime)\right\}
\end{eqnarray*}
for $z, z^\prime \in \mathcal{T}$ \citep{mukerjee2019causal, zhao2021reconciling}.

\begin{comment}
Let $H=\text{diag}(p_0^{-1},p_1^{-1})\otimes1_{2\times 2} - 1_{4\times 4}$, $H_w=\text{diag}(p_0^{-1},p_1^{-1})\otimes \{\text{diag}(q_{w0}^{-1},q_{w1}^{-1}) - 1_{2\times 2} \}$, and $\Psi =W^{-1}\sum_{w=1}^WM_w^{-1}(H_w\circ S_w)$. \cite{zhao2021reconciling} showed that {\color{blue} under mild moment conditions  [---do we need ``split-plot randomization" here---]}, $\sqrt{W} ( \hat \tau_{*} - \tau ) \rightsquigarrow  \mathcal{N}(0,\Sigma_{*, \tau \tau}) $ for $* = \textnormal{ht}, \textnormal{haj}$, where
$
\Sigma_{*, \tau \tau}  = G ( H\circ S_{*}+\Psi) G^\text{T}.
$
{\color{blue}  Here and in what follows, we use the same notation to denote the limit of the covariance in the asymptotic distribution, such as $\Sigma_{*, \tau \tau}$, if no confusion would arise.}

 {\tiny\dotfill}

\begin{itemize}\color{blue}
\item I know 100\% what you mean, but this sentence is still quite confusing. Very hard to introduce the subtlety without mentioning finite population asymptotic scheme.
\item[$\cdot$] Alternatively, split the current Condition 1 into two parts. Give the first half from Zhao and Ding here, and introduce the result under split-plot randomization formally as a lemma.
\item[$\cdot$] With the first half of Condition 1, we can say that Condition xxx ensures that $\Sigma_{*, \tau\tau}$ has a finite limit under the finite population asymptotic scheme. We use the same notation to also denote its limiting value when no confusion would arise. See below for an attempt.
\item[$\cdot$] If this helps, then we introduce the rest of the original Condition 2, and use Conditions 1--2 to replace the original Condition 1.
\end{itemize}

 {\tiny\dotfill}

\end{comment}

Let $H=\text{diag}(p_0^{-1},p_1^{-1})\otimes1_{2\times 2} - 1_{4\times 4}$, $H_w=\text{diag}(p_0^{-1},p_1^{-1})\otimes \{\text{diag}(q_{w0}^{-1},q_{w1}^{-1}) - 1_{2\times 2} \}$, and $\Psi =W^{-1}\sum_{w=1}^WM_w^{-1}(H_w\circ S_w)$.
Let $\overline{\alpha^k} = W^{-1} \sum_{w=1}^{W} \alphaw^k$ be the $k$th moment of $(\alphaw)_{w=1}^{W}$ for $k=1,2,4$, and let
$\overline{Y_{w.}^4(z)} = M_w^{-1} \sum_{s=1}^{M_w} Y_{ws}^4(z)$.
Condition \ref{cond::1} below  was proposed by \citet{zhao2021reconciling} and gives the regularity conditions for {\it finite population asymptotics} under split-plot randomization \citep{DingCLT}.

\begin{condition}\label{cond::1}
As $W$ goes to infinity, for $a,b=0,1$ and $z \in \mathcal{T}$,
\begin{itemize}
    \item[(i)] $\overline{\alpha^2}=O(1); \ \overline{\alpha^4}=o(W)$;
    \item[(ii)] $p_a$ has a limit in $(0,1)$; for all $w=1,\ldots,W$, $q_{wb}\in [c,1-c]$ for a constant $c \in (0,1/2]$ independent of $W$;
    \item[(iii)] for $*= \textnormal{ht},\textnormal{haj}$, $S_{*}$, $\bar Y$, and  $\Psi$  have finite limits;

    \item[(iv)]  $ W^{-1} \max_{w=1, \ldots, W} |\alpha_w \bar{Y}_{w}(z)-\bar{Y}(z)|^2=o(1)$;
       \item[(v)] $W^{-1} \sum_{w=1}^W\alpha_w^2\overline{Y_{w.}^4(z)}=O(1)$; $W^{-2} \sum_{w=1}^W\alpha_w^4\overline{Y_{w.}^4(z)}=o(1)$.
\end{itemize}
\end{condition}

Condition \ref{cond::1}(ii)--(iii) ensure that $\Sigma_{*, \tau\tau} = G ( H\circ S_{*}+\Psi) G^\text{T}$ has a finite
limit for $* = \textnormal{ht}, \textnormal{haj}$. We will use the same notation to also denote their respective limiting values when no confusion would arise.
Lemma \ref{lem:1} below follows from \cite{zhao2021reconciling} and ensures the consistency and asymptotic normality of $\hat \tau_{*} \ (* = \textnormal{ht}, \textnormal{haj})$ for estimating $\tau$.
\begin{lemma}\label{lem:1}
Under Condition \ref{cond::1}, $\sqrt{W} ( \hat \tau_{*} - \tau ) \rightsquigarrow  \mathcal{N}(0,\Sigma_{*, \tau \tau}) $ for $* = \textnormal{ht}, \textnormal{haj}$.
\end{lemma}

\section{Rerandomization in split-plot designs}

\subsection{Rerandomization schemes}
In split-plot designs, we often collect baseline covariates before the experiments, denoted by $x_{ws} = (x_{ws,1},\ldots,x_{ws,L})^\text{T} \in \mathbb{R}^{L}$.
The cluster randomization in stage (I) and stratified randomization in stage (II) ensure that the covariates are balanced on average at both the whole-plot and subplot levels. However, covariate imbalance often exists in finite samples and subjects subsequent inference to possibly large variability and conditional bias. Rerandomization provides a way to balance covariates in the design stage \citep[see, e.g.,][]{Morgan2012,Morgan2015,Li2018,Li2020factorial,Wang2021,zdrep}.
\citet{Morgan2012} suggested a rerandomization scheme using the Mahalanobis distance of the covariate means under different treatment arms to measure the covariate imbalance in a completely randomized treatment-control experiment. This motivates two rerandomization schemes under split-plot randomization.

Specifically, define
$$
 \hat{x}_{\textnormal{ht}} (z) = N^{-1}\sum_{ws \in \mathcal{S}(z)}{p^{-1}_{ws}}(z)x_{ws},  \quad  \hat {x}_{\textnormal{haj}} (z) = \frac{\hat {x}_{\textnormal{ht}} (z)}{\hat {1}_{\textnormal{ht}} (z)},
$$
as the Horvitz--Thompson and Hajek estimators of $\bar x = N^{-1} \sum_{ws \in \mathcal{S}} x_{ws}$ based on units under treatment $z$. Let $\hat x_{*} = ( \hat x_{*}(00), \hat x_{*}(01), \hat x_{*}(10), \hat x_{*}(11) )^\text{T} \in \mathbb{R}^{4 \times L}$ for $* = \textnormal{ht}, \textnormal{haj}$.
For a $4 \times 1$ contrast vector $g = (g_{00}, g_{01}, g_{10}, g_{11})^\text{T}$, the contrast of $\hat x^\text{T}_*(z)$'s,
$$
    g_{00}\hat x^\text{T}_*(00) + g_{01}\hat x^\text{T}_*(01) + g_{10}\hat x^\text{T}_*(10) + g_{11}\hat x^\text{T}_*(11) = g^\text{T} \hat x_* \in \mathbb{R}^{1 \times L},
$$
provides an intuitive measure of covariate balance under split-plot design. A balanced allocation intuitively has homogeneous $\hat x^\text{T}_*(z)$'s such that $g^\text{T} \hat x_*$ is close to $0_{1 \times L}$.
The contrasts that correspond to $g_{\textsc{a}}$, $g_{\textsc{b}}$, and $g_{\textsc{ab}}$ are
\begin{eqnarray*}
    g^\text{T}_{\textsc{a}} \hat x_* &=& 2^{-1}\{ \hat x_{*}^\text{T}(10) + \hat x_{*}^\text{T}(11) \} - 2^{-1}\{ \hat x_{*}^\text{T}(00) + \hat x_{*}^\text{T}(01) \} ,\\
   g^\text{T}_{\textsc{b}} \hat x_* &=& 2^{-1}\{ \hat x_{*}^\text{T}(01) + \hat x_{*}^\text{T}(11) \} - 2^{-1}\{ \hat x_{*}^\text{T}(00) + \hat x_{*}^\text{T}(10) \},\\
    g^\text{T}_{\textsc{ab}} \hat x_* &=& \{ \hat x_{*}^\text{T}(00) + \hat x_{*}^\text{T}(11) \} - \{ \hat x_{*}^\text{T}(01) + \hat x_{*}^\text{T}(10) \},
\end{eqnarray*}
respectively. Let $$ \hat \tau_{*,x} = ( g^\text{T}_{\textsc{a}} \hat x_*, g^\text{T}_{\textsc{b}} \hat x_*, g^\text{T}_{\textsc{ab}} \hat x_*)^\text{T} \in \mathbb{R}^{3L}$$
be their concatenation for $* = \textnormal{ht}, \textnormal{haj}$, which is intuitively close to $0_{3L}$ if the allocation is balanced.
We consider two rerandomization schemes based on the Mahalanobis distance between $\hat\tau_{*,x}$ and $0_{3L}$
under split-plot randomization.

The first scheme is based on the Mahalanobis distance between $\hat \tau_{\HT,x}$ and $0_{3L}$
under split-plot randomization: $\textsc{M}_{\textnormal{ht}}=\hat \tau_{\HT,x}^\text{T} \textnormal{cov}( \hat \tau_{\HT,x} )^{-1}  \hat \tau_{\HT,x}$. For a predetermined threshold $d > 0$,
rerandomization  accepts the treatment assignment if and only if {the following event happens:
$$
\mathcal{M}_{\textnormal{ht}} = \{ \textsc{M}_{\textnormal{ht}} \leq d \}.
$$
The second scheme is based on the Mahalanobis distance  between $\hat \tau_{\haj,x}$ and $0_{3L}$
under split-plot randomization: $\textsc{M}_{\textnormal{haj}}=\hat \tau_{\haj,x}^\text{T} \textnormal{cov}_\textup{a}( \hat \tau_{\haj,x} )^{-1}  \hat \tau_{\haj,x}$, and accepts the treatment assignment if and only if the following event happens:
$$
\mathcal{M}_{\textnormal{haj}} = \{ \textsc{M}_{\textnormal{haj}}\leq d \}.
$$
We define $\textsc{M}_{\textnormal{haj}}$ using the asymptotic covariance $\textnormal{cov}_\textup{a}( \hat \tau_{\haj,x} )$ due to the complicated form of the exact covariance $\textnormal{cov}(\hat \tau_{\haj,x})$; see Theorem~\ref{thm::clt} in Section~\ref{sec::asy} for more details.

Two treatment effect estimators, $\hat \tau_{\HT}$ and $\hat \tau_{\haj}$, and two rerandomization schemes, $\mathcal{M}_{\textnormal{ht}}$ and $\mathcal{M}_{\textnormal{haj}}$, give rise to four inferential strategies as their combinations. Nevertheless, it is more natural to consider design and analysis of the same type. Therefore,
{we will consider $\mathcal{M}_{\textnormal{ht}}$ for rerandomization  if using $\hat \tau_{\HT}$ for treatment effect estimation, and consider $\mathcal{M}_{\textnormal{haj}}$ if using $\hat \tau_{\haj}$.} To avoid confusion, we will henceforth use {\it classic split-plot randomization} to refer to the standard split-plot randomization without rerandomization.

\subsection{Asymptotic distribution}\label{sec::asy}
For $*=\textnormal{ht},\textnormal{haj}$, the asymptotic distribution of $\hat \tau_* $ under rerandomization scheme $\mathcal{M}_{*} $ is essentially the conditional asymptotic distribution of $\hat \tau_* $ under classic split-plot randomization given $\mathcal{M}_{*} $, denoted by $\hat \tau_* \mid \mathcal{M}_{*} $ \citep{Li2018}.
To study them, we start with the unconditional joint asymptotic distributions of  $(\hat\tau^\text{T}_*, \hat \tau^\text{T}_{*,x})^\text{T}$ under classic split-plot randomization.

Let $S_{\textnormal{ht},xx}$, $S_{\textnormal{haj},xx}$, $S_{w,xx}$, $S_{\textnormal{ht},xY(z)}$, $S_{\textnormal{haj},xY(z)}$, and $S_{w,xY(z)}$ be the scaled between and within whole-plot covariances of $(x_{ws})_{ws \in \mathcal{S}}$ with itself and with $\{Y_{ws}(z)\}_{ws \in \mathcal{S}}$, respectively, analogous to $S_{\textnormal{ht}}(z,z')$, $S_\textnormal{haj}(z,z')$, and $S_{w}(z,z')$.
To avoid too many formulas in the main paper, we relegate their explicit forms
to the supplementary materials.
Define
$$
    \Psi_{xx}=W^{-1}\sum_{w=1}^{W} M_w^{-1} (H_w\otimes S_{w,xx}), \quad
    \Psi_{ xY}=W^{-1}\sum_{w=1}^{W} M_w^{-1} (H_w\otimes 1_L)\circ(1_4\otimes S_{w,xY}),
$$
where $S_{w,xY} = ( S_{w,xY(00)}, S_{w,xY(01)}, S_{w,xY(10)}, S_{w,xY(11)}) \in \mathbb{R}^{L \times 4}$. For $* = \textnormal{ht}, \textnormal{haj}$, let
\begin{eqnarray*}
S_{*, xY} &=& ( S_{*,xY(00)}, S_{*,xY(01)}, S_{*,xY(10)}, S_{*,xY(11)}) \in \mathbb{R}^{L \times 4}, \\
\Sigma_{*, xx} & =  &  ( G \otimes I_L ) (H\otimes S_{*,xx}+\Psi_{xx})  ( G \otimes I_L)^\text{T},\\
\Sigma_{*, x\tau} &=& \Sigma_{*, \tau x}^\text{T} = ( G \otimes I_L) \{ (H\otimes 1_L)\circ(1_4\otimes S_{*,xY}) + \Psi_{xY} \} G^\text{T}.
\end{eqnarray*}
We require Condition~\ref{cond::1} and Condition~\ref{cond::2} below for deriving the joint asymptotic distribution of  $(\hat\tau^\text{T}_*, \hat \tau^\text{T}_{*,x})^\text{T}$ for $* = \textnormal{ht}, \textnormal{haj}$.
Let $\bar{x}_w = M_w^{-1}\sum_{s=1}^{M_w} x_{ws}$ and $\overline{ ||  x_{w.} ||_\infty^4 } = M_w^{-1} \sum_{s=1}^{M_w} || x_{ws} ||_\infty^4$.

\begin{condition}\label{cond::2}
As $W$ goes to infinity,
\begin{itemize}
    \item[(i)] for $*= \textnormal{ht},\textnormal{haj}$,  $S_{*,xx}, S_{*,xY}$,  $\Psi_{xx}$, and $\Psi_{xY}$  have finite limits; the limits of $S_{*,xx}$ and $\Sigma_{*, xx}$ are invertible;

    \item[(ii)] $ W^{-1} \max_{w=1, \ldots, W} || \alpha_w \bar{x}_{w}-\bar{x}||_\infty^2=o(1)$;
    \item[(iii)] $W^{-1} \sum_{w=1}^W\alpha_w^2\overline{|| x_{w.} ||_\infty^4}=O(1)$; $W^{-2} \sum_{w=1}^W\alpha_w^4\overline{|| x_{w.} ||_\infty^4}=o(1)$.
\end{itemize}
\end{condition}

Condition~\ref{cond::2} gives the analog of Condition~\ref{cond::1} for the covariates $x_{ws}$'s.
Condition~\ref{cond::1}(ii) and Condition~\ref{cond::2}(i) together ensure that $\Sigma_{*, xx}$, $\Sigma_{*, x \tau}$, and $\Sigma_{*, \tau x}$ all have finite limits for $* = \textnormal{ht}, \textnormal{haj}$. Again, we will use the same notation to also denote their respective limiting values when no confusion would arise.

\begin{theorem}\label{thm::clt}
Under Conditions~\ref{cond::1} and \ref{cond::2}, for $* = \textnormal{ht}, \textnormal{haj}$,
\begin{eqnarray}
   \sqrt{W} \left( \begin{array}{cc}
        \hat\tau_* - \tau  \\
        \hat \tau_{*,x}
   \end{array} \right) \rightsquigarrow
   \mathcal{N} \left( 0, \Sigma_* \right), \quad \Sigma_* = \left(\begin{array}{cc}
        \Sigma_{*, \tau\tau} & \Sigma_{*, \tau x}\\
        \Sigma_{*, x\tau} & \Sigma_{*, xx}\\
   \end{array} \right). \nonumber
\end{eqnarray}
\end{theorem}

Theorem~\ref{thm::clt} ensures the asymptotic joint normality of  $ \hat \tau_{*}$ and $\hat \tau_{*,x}$, and provides the basis for deriving the conditional asymptotic distribution of $\hat\tau_*$ given $\mathcal{M}_{*}$.
By Theorem~\ref{thm::clt}, the Mahalanobis distance
$
\textsc{M}_{*}= (\sqrt{W} \hat \tau_{*,x})^\text{T} \Sigma_{*, xx}^{-1}   (\sqrt{W} \hat \tau_{*,x})
$
converges in distribution to $\chi_{3L}^2$ for both $* = \textnormal{ht}, \textnormal{haj}$. Thus, we can choose the threshold $d$
as the $\alpha$th quantile of $\chi_{3L}^2$ to ensure an asymptotic acceptance rate of $\alpha$ for the rerandomization.

By Theorem \ref{thm::clt}, the linear projection of $\sqrt{W} \hat \tau_{*}$ onto $ \hat \tau_{*,x}$ equals
$
\textup{proj}(\sqrt{W} \htau_* \mid  \hat\tau_{*, x}) = \sqrt{W} \tau + \sqrt{W} \Sigma_{*, \tau x} \Sigma^{-1}_{*, xx}\hat\tau_{*, x}$ asymptotically.
Let $\Sigma_{*, \tau\tau}^{||} = \Sigma_{*, \tau x}\Sigma_{*,xx}^{-1}\Sigma_{*,x\tau}$ denote the asymptotic covariance of $ \textup{proj}( \sqrt{W}  \htau_* \mid \hat\tau_{*, x})$, and let $\Sigma_{*, \tau\tau}^{\bot}=\Sigma_{*, \tau\tau} - \Sigma_{*, \tau\tau}^{||}$ denote that of the residual
$\textup{res}(\sqrt{W}  \htau_* \mid \hat\tau_{*, x}) = \sqrt{W} \htau_* - \textup{proj}(\sqrt{W} \htau_* \mid \hat\tau_{*, x})$.

\begin{theorem}\label{thm::rerand}
Under Conditions~\ref{cond::1} and \ref{cond::2},  for $* = \textnormal{ht}, \textnormal{haj}$,
$$
\sqrt{W}(\hat\tau_{*}- \tau) \mid \mathcal{M}_{*} \rightsquigarrow  (\Sigma_{*,\tau\tau}^{\bot})^{1/2}\epsilon
            +\Sigma_{*,\tau x}\Sigma_{*,xx}^{-1/2}\zeta_{3L, d},
$$
where $\epsilon \sim \mathcal{N}(0, I_{3})$ is a 3-dimensional standard normal random vector, $\zeta_{3L, d}\sim D \mid D^\text{T} D \leq d$ is a $3L$-dimensional truncated normal random vector with $D \sim \mathcal{N}(0,I_{3L})$,
and $\epsilon$ and $\zeta_{3L, d}$ are independent.
\end{theorem}

Theorem~\ref{thm::rerand} indicates that the asymptotic distribution of $\hat \tau_{*}$ under rerandomization scheme $\mathcal{M}_{*}$ is the convolution of a normal component and a truncated normal component.
Observe that
\begin{eqnarray*}
\sqrt{W} (\htau_* - \tau)
&=& \textup{proj}(\sqrt{W} \htau_* \mid \hat\tau_{*, x}) + \textup{res}( \sqrt{W} \htau_* \mid \hat\tau_{*, x}) - \sqrt{W} \tau \\
&=& \sqrt{W} \Sigma_{*, \tau x} \Sigma^{-1}_{*, xx}\hat\tau_{*, x} + \textup{res}(\sqrt{W}\htau_* \mid \hat\tau_{*, x}).
\end{eqnarray*}
The term $\textup{res}(\sqrt{W}\htau_* \mid \hat\tau_{*, x})$ is asymptotically independent of $\hat\tau_{*, x}$ under split-plot randomization, and corresponds to the normal vector $(\Sigma_{*,\tau\tau}^{\bot})^{1/2}\epsilon$ unaffected by the rerandomization.
The term $\sqrt{W} \Sigma_{*, \tau x} \Sigma^{-1}_{*,xx}\hat\tau_{*, x}$ is affected by the rerandomization and corresponds to the truncated normal vector $\Sigma_{*,\tau x}\Sigma_{*,xx}^{-1/2}\zeta_{3L, d}$.
It extends the asymptotic theory of rerandomization with treatments assigned at only the unit level \citep{Li2018,Li2020factorial,Wang2021} or group level   \citep{Lu2022} to the split-plot designs.
Moreover, the asymptotic distributions in Theorem~\ref{thm::rerand} are central convex unimodal \citep[][Definition 2 and Proposition 2]{Li2020factorial}.

We use the following notion of peakedness \citep{sherman} to quantify the relative efficiency between different estimators \citep{Li2020factorial, zdrep}.

\begin{definition}
For two symmetric $m$-dimensional random vectors $U_1$ and $U_2$, we say that $U_1$ is more peaked than $U_2$ if  $\textnormal{pr} (U_1 \in \mathcal{K}) \geq \textnormal{pr} (U_2 \in \mathcal{K} )$ for every symmetric convex set $\mathcal{K} \subset \mathbb{R}^m$.
\end{definition}

Peakedness implies not only smaller covariance, but also  narrower central quantile regions.
It hence provides a more refined measure than covariance for comparing relative efficiency between estimators with nonnormal asymptotic distributions.
For $* = \textnormal{ht}, \textnormal{haj}$, we say that rerandomization improves the asymptotic efficiency of $\htau_* $ if the asymptotic distribution of $\htau_*  - \tau$ under rerandomization, namely $\htau_* -\tau \mid \mathcal M_*$, is more peaked than that of  $\htau_*  - \tau$ without rerandomization. Corollary~\ref{cor::varcom} below shows the improvement of asymptotic efficiency of $\htau_* $ by split-plot rerandomization.

\begin{corollary}\label{cor::varcom}
Under Conditions~\ref{cond::1} and \ref{cond::2}, for $*=\textnormal{ht},\textnormal{haj}$,
rerandomization by $\mathcal{M}_{*}$ improves the asymptotic efficiency of $\htau_*$
with
\begin{eqnarray*}
    W \left[\textnormal{cov}_\textup{a}( \hat \tau_{*}) - \textnormal{cov}_\textup{a}( \hat \tau_{*} \mid \mathcal{M}_{*}) \right] = ( 1 - r_{3L, d}) \Sigma_{*, \tau\tau}^{||} \geq 0,
\end{eqnarray*}
where $r_{3L, d} = \textnormal{pr}( \chi^2_{3L + 2} \leq d )/ \textnormal{pr}( \chi^2_{3L} \leq d ) \leq 1$.
\end{corollary}

\subsection{Estimation of the asymptotic distribution}

By Theorem~\ref{thm::rerand}, to infer $\tau$ based on $\hat \tau_{*}$ under rerandomization scheme $\mathcal{M}_{*}$, we need to estimate $\Sigma_{*, \tau\tau}^{\bot}$ and $\Sigma_{*, x\tau}$
for $*=\textnormal{ht},\textnormal{haj}$. By definition, it suffices to estimate $\Sigma_{*,\tau \tau}$ and $\Sigma_{*, x \tau}$.

Let $\hat Y_w(z) = M_{wb}^{-1} \sum_{s:Z_{ws} = z} Y_{ws}$ be the whole-plot sample mean under treatment $z=(ab)$, and let $A_{w}$ be the level of factor A received by whole plot $w$.
Define
\begin{eqnarray*}
\hat S_{\textnormal{ht}}(z,z^\prime) & = & (W_a-1)^{-1}\sum_{w: A_w=a}\left\{\alpha_w \hat Y_w(z) - \hat Y_{\textnormal{ht}}(z) \right\} \left\{\alpha_w \hat Y_w(z^\prime) - \hat Y_{\textnormal{ht}}(z^\prime)\right\},   \\
\hat S_{\textnormal{haj}}(z,z^\prime) & = & (W_a-1)^{-1}\sum_{w: A_w=a} \alphaw^2 \left\{ \hat Y_w(z) - \hat Y_{\textnormal{haj}}(z) \right\} \left\{ \hat Y_w(z^\prime) - \hat Y_{\textnormal{haj}}(z^\prime)\right\}
\end{eqnarray*}
as the sample analogs of $S_{\textnormal{ht}}(z,z^\prime)$ and $S_{\textnormal{haj}}(z,z^\prime)$ for $z=(ab)$ and $z^\prime=(ab^\prime)$ with the same level of factor A.
For $* = \textnormal{ht}, \textnormal{haj}$,
\citet[][Theorem 4.2]{zhao2021reconciling} ensures that
$$
\hat\Sigma_{*, \tau\tau}=G\left(
            \begin{array}{cc}
            p_0^{-1}\left(
            \begin{array}{cc}
            \hat S_{*}(00,00) &  \hat S_{*}(00,01) \\
            \hat S_{*}(00,01) &  \hat S_{*}(01,01) \\
            \end{array}\right)
            & 0_{2\times 2}\\
            0_{2\times 2} & p_1^{-1}\left(
            \begin{array}{cc}
            \hat S_{*}(10,10) &  \hat S_{*}(10,11) \\
            \hat  S_{*}(10,11) & \hat  S_{*}(11,11) \\
            \end{array}\right)\\
            \end{array}\right)G^\text{T}
$$
gives an asymptotically conservative estimator of $\Sigma_{*, \tau\tau}$ under classic split-plot randomization.

Let $\hat Y_{\textnormal{ht},ws}(z) =  \mathcal{I}(Z_{ws}=z) p_{ws}(z)^{-1}  Y_{ws}$ and $\hat Y_{\textnormal{ht},w}(z) = M_{w}^{-1} \sum_{s=1}^{M_w} \mathcal{I}(Z_{ws}=z) p_{ws}(z)^{-1}  Y_{ws}$
be the Horvitz--Thompson estimators of $Y_{ws}(z)$ and $\bar Y_w(z)$, respectively.
Let $\hat S_{\textnormal{ht},xY}$, $\hat S_{\textnormal{haj},xY}$, and $\hat S_{w,xY}$ be the sample analogs of $S_{\textnormal{ht},xY}$, $S_{\textnormal{haj},xY}$, and $S_{w,xY}$, respectively, with $Y_{ws}(z)$, $\bar Y_w(z)$, and $\bar Y(z)$ estimated by $\hat Y_{\textnormal{ht},ws}(z)$, $\hat Y_{\textnormal{ht},w}(z)$, and $\hat Y_{\textnormal{ht}}(z)$, respectively.
We can then estimate $\Sigma_{*,  x\tau}$ by
\begin{eqnarray*}
\hat\Sigma_{*,  x\tau}  = (G\otimes I_L)\left\{(H\otimes 1_L)\circ(1_4\otimes \hat S_{*, xY}) + \hat\Psi_{ xY}\right\}G^\text{T},
\end{eqnarray*}
where $\hat\Psi_{xY} = W^{-1}\sum_{w=1}^{W} M_w^{-1}(H_w\circ \hat S_{w, xY})$.
This yields
$$
 \hat\Sigma_{*} = \left(
            \begin{array}{cc}
            \hat\Sigma_{*, \tau\tau} & \hat\Sigma_{*,  \tau x}\\
            \hat\Sigma_{*,  x \tau} & \Sigma_{*, xx}\\
            \end{array}\right),
$$
where $\hat\Sigma_{*,  \tau x}  = \hat\Sigma^\text{T}_{*,  x \tau}$, as a plug-in estimator of $\Sigma_{*} \ (*  = \textnormal{ht}, \textnormal{haj})$.

\begin{theorem}\label{thm::varest}
Under Conditions~\ref{cond::1} and \ref{cond::2},
for $* = \textnormal{ht}, \textnormal{haj}$,
    \[
        \begin{split}
            ( \hat \Sigma_{*} - \Sigma_{*} )  \mid \mathcal{M}_{*} = \left(
            \begin{array}{cc}
            GS_{*}G^\text{T} & 0_{3 \times 3L}\\
            0_{3L \times 3} & 0_{3L \times 3L}\\
            \end{array}\right)
            + o_{\mathbb P}(1).
        \end{split}
    \]
\end{theorem}

As $GS_{*}G^\text{T}$ is positive semi-definite, Theorem~\ref{thm::varest} shows that $\hat \Sigma_{*, \tau\tau}$ is an asymptotically conservative estimator of $\Sigma_{*, \tau\tau}$ and $\hat\Sigma_{*, x\tau}$ is a consistent estimator of $\Sigma_{*, x\tau}$ under split-plot rerandomization. Thus, $\hat \Sigma_{*, \tau\tau}^{||} = {\hat \Sigma_{*,\tau x}} \Sigma_{*,xx}^{-1} \hat \Sigma_{*,x\tau} $ is a consistent estimator of $\Sigma_{*, \tau\tau}^{||} $ and $\hat \Sigma_{*, \tau\tau}^{\bot} = \hat \Sigma_{*, \tau\tau} - \hat \Sigma_{*, \tau\tau}^{||}$ is a conservative estimator of $\Sigma_{*, \tau\tau}^{\bot}$. Therefore, the asymptotic distribution of $\sqrt{W}( \hat \tau_{*} - \tau ) \mid \mathcal{M}_{*}$ can be conservatively estimated by
$
 \phi_* = (\hat \Sigma_{*,\tau\tau}^{\bot})^{1/2}\epsilon
            + \hat \Sigma_{*,\tau x}\Sigma_{*,xx}^{-1/2}\zeta_{3L, d}.
$
Suppose that the limit of $\Sigma_{*, \tau\tau}^{\bot}$ is invertable, then $\hat \Sigma_{*,\tau\tau}^{\bot}$ is invertable with probability tending to one. Let $\hat c_{*, 1 - \xi}$ and $\chi^2_{3, 1 - \xi}$  ($0 < \xi < 1$) be the $1 - \xi$ quantiles of $ \phi_*^\text{T} (\hat \Sigma_{*,\tau\tau}^{\bot})^{-1}  \phi_*$ and $\chi^2_{3}$,
respectively. Corollary~\ref{cor::CI} below provides asymptotically conservative confidence regions for $\tau$ and demonstrates that rerandomization generally improves the inference efficiency.

\begin{corollary}\label{cor::CI}
Suppose that the limit of $\Sigma_{*, \tau\tau}^{\bot}$ is invertable. Under Conditions~\ref{cond::1} and \ref{cond::2}, for $* = \textnormal{ht}, \textnormal{haj}$, the  Wald-type confidence region
$
\{ \tau:  W (\hat \tau_{*} - \tau)^\text{T} (\hat \Sigma_{*,\tau\tau}^{\bot})^{-1} (\hat \tau_{*} - \tau) \leq \hat c_{*, 1 - \xi}  \}
$
has asymptotic coverage rate greater than or equal to $1 - \xi$ under the corresponding  split-plot rerandomization scheme. Moreover, the area of the above confidence region is smaller than or equal to that of the confidence region
$
\{ \tau:  W (\hat \tau_{*} - \tau)^\text{T} \hat \Sigma_{*,\tau\tau}^{-1} (\hat \tau_{*} - \tau) \leq \chi^2_{3, 1 - \xi}  \}
$
under the classic split-plot randomization.
\end{corollary}

\section{Covariate Adjustment under Rerandomization}

The discussion so far concerned rerandomization that enforces covariate balance in the design stage.
Alternatively, we can adjust for covariate imbalance in the analysis stage.
\cite{Li2020} and \cite{Wang2021} showed the duality of rerandomization and regression adjustment for improving efficiency in completely randomized and stratified treatment-control experiments, respectively.
In this section, we extend the discussion to the method and design-based theory of the combination of rerandomization and covariate adjustment in $2^2$ split-plot designs.
We consider two strategies for covariate adjustment
 for each of the Horvitz--Thompson and Hajek estimators, and derive their design-based properties under split-plot rerandomization.
The first strategy follows the regression formulation by \cite{zhao2021reconciling}.
The second strategy is new and approaches covariate adjustment from a projection or conditional inference perspective.

Let $v_{ws} \in \mathbb{R}^{J}$ denote the covariates used in the analysis stage.
We allow the analysis stage to use more covariates than the design stage
in the sense that $ x_{ws} = C v_{ws} $ for some matrix $C \in \mathbb R^{L \times J}$ ($J \geq L$).
Let $\bar v = N^{-1} \sum_{ws \in \mathcal S} v_{ws}$, $\bar v_{w} = M_w^{-1} \sum_{s=1}^{M_w} v_{ws}$, and $\hat v_w (z) = M_{wb}^{-1} \sum_{s: Z_{ws}=z} v_{ws}$
for $w=1,\ldots,W$ and $z = (ab) \in \mathcal {T}$.
For $*=\textnormal{ht},\textnormal{haj}$, define  $S_{*,vv}$, $S_{w,vv}$, $\Psi_{vv}$, $S_{*,vY(z)}$, $S_{w,vY(z)}$, $S_{*,vY}$, $S_{w,vY}$, $\Psi_{ vY}$, $\Sigma_{*, vv}$, $\Sigma_{*,v \tau}$, $\Sigma_{*, \tau v}$, $\hat {v}_{*}(z)$, $\hat {v}_{*}$, and $\hat \tau_{*, v}$ similarly to $S_{*,xx}$, $S_{w,xx}$, $\Psi_{xx}$, $S_{*,xY(z)}$, $S_{w,xY(z)}$, $S_{*,xY}$, $S_{w,xY}$, $\Psi_{xY}$, $\Sigma_{*,xx}$, $\Sigma_{*,x \tau}$, $\Sigma_{*, \tau x}$, $\hat {x}_{*}(z)$, $\hat {x}_{*}$, and $\hat \tau_{*, x}$, with $x_{ws}$ replaced by $v_{ws}$.

\subsection{Regression with treatment-covariate interactions}

Regression adjustment provides a convenient
way to adjust for covariate imbalance in the analysis stage.
For observed data $\{(y_i, u_i):  i \in \mathcal J, \ y_i \in \mathbb R, \ u_i \in \mathbb R^m\}$, where $\mathcal J$ denotes the index set, denote by $y_i \sim u_i$  the linear regression of $y_i $ on $u_i$ over $i \in \mathcal J$.
\cite{zhao2021reconciling} showed that the Horvitz--Thompson and Hajek estimators $\hat\tau_* \ (* = \textnormal{ht}, \textnormal{haj})$ can be recovered from the ordinary least squares (ols) fit of the aggregate regression
\begin{eqnarray}\label{eq:ag}
        \alphaw \hat Y_w(A_w b) \sim \mathcal{I}(A_w b =00)+ \mathcal{I}(A_w b=01)+ \mathcal{I}(A_w b=10)+ \mathcal{I}(A_w b=11)
\end{eqnarray}
over $\{(w,b): w = 1,\ldots, W; \ b = 0,1\}$ and the weighted least squares (wls) fit of
\begin{eqnarray}\label{eq:wls}
        Y_{ws} \sim \mathcal{I}(Z_{ws}=00)+ \mathcal{I}(Z_{ws}=01)+ \mathcal{I}(Z_{ws}=10)+ \mathcal{I}(Z_{ws}=11)
\end{eqnarray}
over $ws \in \mathcal S$, respectively, and recommended including full interactions between the treatment indicators and centered covariates for regression adjustment.
In particular, let $\hat \beta_\ag $ and $\hat V_\ag$
be the ols coefficient vector and associated cluster-robust covariance from \eqref{eq:ag}, where we use
the subscript ``ag" to signify the use of whole-plot aggregate outcomes and covariates in forming the regression \citep{Abadie2008,Basse2018, imai2020, Su2021}.
Let $\hat \beta_\wls$ and $\hat V_\wls $ be the wls coefficient vector and associated cluster-robust covariance from \eqref{eq:wls}, where we weight subplot $ws$ by the inverse of its realized inclusion probability $p_{ws} (Z_{ws})$.
\cite{zhao2021reconciling} showed that
$\hat \beta_\ag  = \hat{Y}_{\textnormal{ht}}$ and $\hat \beta_\wls = \hat{Y}_{\textnormal{haj}}$, with $\hat V_\ag $ and $\hat V_\wls $ being asymptotically conservative for estimating the true sampling covariances. This justifies the large-sample Wald-type inference of $\tau$ based on $(G\hat \beta_\ag, G\hat V_\ag G^\text{T})$ and $(G\hat \beta_\wls, G\hat V_\wls G^\text{T})$.
Further let
\begin{eqnarray}\label{eq:ag_v}
 \alphaw \hat Y_w(A_w b) & \sim & \sum_{z \in \mathcal{T} } \mathcal{I}(A_w b =z) + \sum_{z \in \mathcal{T} } \mathcal{I}(A_w b =z)\alphaw \{\hat v_w(A_w b) - \bar v\}, \\
 Y_{ws} & \sim &  \sum_{z \in \mathcal{T} } \mathcal{I}(Z_{ws} =z) + \sum_{z \in \mathcal{T} } \mathcal{I}(Z_{ws} =z) ( v_{ws} - \bar v) \label{eq:wls_v}
\end{eqnarray}
be the fully interacted variants of \eqref{eq:ag} and \eqref{eq:wls}.
Let $\hat \beta_{\ag, \L} $ and $\hat \beta_{\wls, \L} $ denote the ols and wls coefficient vectors of $\{\mathcal{I}(A_w b =z)\}_{z \in \mathcal{T}}$ and $\{\mathcal{I}(Z_{ws} =z)\}_{z \in \mathcal{T}} $ from \eqref{eq:ag_v} and \eqref{eq:wls_v}, respectively, with $\hat V_{\ag, \L} $ and $\hat V_{\wls, \L} $ as the associated cluster-robust covariances.
They form the regression-adjusted counterparts of $(\hat \beta_\ag , \hat V_\ag )$ and $(\hat \beta_\wls , \hat V_\wls )$.
We use the subscript ``$\textsc{l}$" to signify \cite{lin2013}, who proposed the fully interacted adjustment under completely randomized experiments.

Let $\hat \tau_{\HT, \L} = G \hat \beta_{\ag, \L}  $ and $\hat \tau_{\haj, \L} = G \hat \beta_{\wls, \L} $ be the corresponding {\it regression-adjusted Horvitz--Thompson and Hajek estimators}  of $\tau$, with ${ \hat \Sigma}_{\textnormal{ht},\textsc{l}, \tau \tau} = W G \hat V_{\ag, \L}  G^\text{T} $ and ${ \hat \Sigma}_{\textnormal{haj},\textsc{l}, \tau \tau} = W G \hat V_{\wls, \L}  G^\text{T} $ as the associated cluster-robust covariance estimators up to a factor of $W$.
\citet[][Theorem 6.2]{zhao2021reconciling} ensured the asymptotic validity of $(\hat\tau_{*,\textsc{l}},   \hat \Sigma_{*,\textsc{l}, \tau \tau})$  for inferring $\tau$ under the classic split-plot randomization.
Theorem~\ref{thm::rerandreg} below extends their results and presents the asymptotic properties of $(\hat\tau_{*,\textsc{l}},   \hat \Sigma_{*,\textsc{l}, \tau \tau})$ under split-plot rerandomization.

Let $\hat \gamma_{\textnormal{ag},z}$ and $\hat \gamma_{\textnormal{wls},z}$ be the coefficient vectors of $\mathcal{I}(A_w b =z)\alphaw  \{ \hat v_w(A_w b) - \bar v\}$ and $\mathcal{I}(Z_{ws} =z) ( v_{ws} - \bar v)$ from the ols and wls fits of \eqref{eq:ag_v} and \eqref{eq:wls_v}, respectively.
Under Condition~\ref{cond::3} below, $\hat \gamma_{\textnormal{ag},z}$ and $\hat \gamma_{\textnormal{wls},z}$ have finite probability limits, denoted by $\gamma_{\textnormal{ag},z}$ and $\gamma_{\textnormal{wls},z}$ respectively, under split-plot rerandomization.
We give the exact formulas of $\gamma_{\textnormal{ag},z}$ and $\gamma_{\textnormal{wls},z}$ in the supplementary materials.
Define covariate-adjusted potential outcomes $Y_{ws}(z; \gamma_{\dagger,z}) = Y_{ws}(z) - (v_{ws} - \bar v)^\text{T} \gamma_{\dagger,z}$ for $\dagger = \textnormal{ag}, \textnormal{wls}$ and $z \in \mathcal{T}$.
Define
\begin{eqnarray*}
\Sigma_{*,\textsc{l}} = \left(\begin{array}{cc}
        \Sigma_{*, \textsc{l}, \tau\tau} & \Sigma_{*, \textsc{l}, \tau x}\\
        \Sigma_{*, \textsc{l}, x \tau} &  \Sigma_{*, xx} \\
   \end{array} \right) \quad (* = \textnormal{ht}, \textnormal{haj})
\end{eqnarray*}
similarly to $\Sigma_{*}$ with $Y_{ws}(z)$ replaced by $Y_{ws}(z; \gamma_{\textnormal{ag},z})$ and $Y_{ws}(z; \gamma_{\textnormal{wls},z})$, respectively, for $* = \textnormal{ht}$ and $* = \textnormal{haj}$.
Applying Theorem \ref{thm::clt} to the covariate-adjusted potential outcomes ensures that $\Sigma_{*, \textsc{l}}$ gives the asymptotic covariance matrix of $\sqrt W(\hat \tau_{*,\textsc{l}}^\text{T}, \hat \tau_{*,x}^\text{T})^\text{T}$.
The $\hat \Sigma_{*,\textsc{l},\tau\tau}$ from regression thus gives a convenient estimator of  $\Sigma_{*, \textsc{l}, \tau\tau} = W\textnormal{cov}_\textup{a}(\hat \tau_{*,\textsc{l}})$.

Let $Q_{vv} = (N-1)^{-1} \sum_{ws \in \mathcal{S}} (v_{ws}-\bar v)(v_{ws}-\bar v)^\text{T}$ and $Q_{vY(z)} = (N-1)^{-1} \sum_{ws \in \mathcal{S}}  (v_{ws} - \bar v)Y_{ws}(z)$ be the finite population covariances of $(v_{ws})_{ws \in \mathcal{S}}$ with itself and $\{Y_{ws}(z)\}_{ws \in \mathcal{S}}$, respectively.

\begin{condition}\label{cond::3} (i) Condition~\ref{cond::2} holds with $x_{ws}$ replaced by $v_{ws}$; (ii) as $M$ goes to infinity, $Q_{vv}$ and $Q_{vY(z)}$
have finite limits, and the limit of $Q_{vv}$ is invertible.
\end{condition}

Conditions~\ref{cond::1}--\ref{cond::3} ensure that $\Sigma_{*, \textsc{l}, \tau\tau}$, $ \Sigma_{*, \textsc{l}, \tau x}$,  $\Sigma_{*, \textsc{l}, x \tau}$, and $ \Sigma_{*,  \textsc{l}, xx}$ all have finite limits  for $* = \textnormal{ht}, \textnormal{haj}$. We will use the same notation to also denote their respective limiting values when no confusion would arise. Recall that $\hat \Sigma_{*,\textsc{l},\tau\tau}$ gives a convenient estimator of  $\Sigma_{*, \textsc{l}, \tau\tau} = W\textnormal{cov}_\textup{a}(\hat \tau_{*,\textsc{l}})$.
Let $\Sigma_{*, \textsc{l}, \tau\tau}^{||} = \Sigma_{*, \textsc{l},\tau x}\Sigma_{*,xx}^{-1}\Sigma_{*,\textsc{l},x\tau}$ denote the covariance of the linear projection of $\sqrt{W} \hat \tau_{*,\textsc{l}}$ onto $\hat \tau_{*,x}$ analogous to $\Sigma_{*, \tau\tau}^{||}$, and let $\Sigma_{*, \textsc{l},\tau\tau}^{\bot}=\Sigma_{*, \textsc{l},\tau\tau} - \Sigma_{*, \textsc{l}, \tau\tau}^{||}$ denote the corresponding covariance of the residual.
Let $\hat \Sigma_{*, \textsc{l}, \tau x} = \hat \Sigma_{*, \textsc{l},  x\tau}^\text{T}$  be the plug-in estimators of $ \Sigma_{*, \textsc{l}, \tau x}=  \Sigma_{*, \textsc{l},  x\tau}^\text{T}$, which are defined similarly to $\hat \Sigma_{*, \tau x} = \hat \Sigma_{*,  x\tau}^\text{T}$ with  $Y_{ws}(z)$ replaced by $Y_{ws}(z; \hat \gamma_{\textnormal{ag},z})$ and $Y_{ws}(z; \hat \gamma_{\textnormal{wls},z})$, respectively, for $* = \textnormal{ht}$ and $* = \textnormal{haj}$.
Let $\hat \Sigma_{*, \textsc{l},\tau\tau}^{\bot} = \hat\Sigma_{*, \textsc{l},\tau\tau} - \hat\Sigma_{*, \textsc{l},\tau x}\Sigma_{*,xx}^{-1}\hat\Sigma_{*,\textsc{l},x\tau}$ be the corresponding estimator of $\Sigma_{*, \textsc{l},\tau\tau}^{\bot}$.

\begin{theorem}\label{thm::rerandreg}
Under Conditions~\ref{cond::1}--\ref{cond::3},
for $* = \textnormal{ht}, \textnormal{haj}$,
\begin{eqnarray*}
 \sqrt{W}(\hat\tau_{*, \textsc{l}}- \tau) \mid \mathcal{M}_{*} &\rightsquigarrow & (\Sigma_{*, \textsc{l}, \tau\tau}^{\bot})^{1/2}\epsilon+\Sigma_{*, \textsc{l}, \tau x}\Sigma_{*,  xx}^{-1/2}\zeta_{3L, d},
\end{eqnarray*}
where $\epsilon \sim \mathcal{N}(0, I_{3})$ is a 3-dimensional standard normal random vector, $\zeta_{3L, d}\sim D \mid D^\text{T} D \leq d$ is a $3L$-dimensional truncated normal random vector with $D \sim \mathcal{N}(0,I_{3L})$,
and $\epsilon$ and $\zeta_{3L, d}$ are independent. Moreover,
  \[
        \begin{split}
         ( { \hat \Sigma}_{*,\textsc{l}, \tau \tau} -  \Sigma_{*,\textsc{l}, \tau \tau} )  \mid \mathcal{M}_{*} &=  GS_{*,\textsc{l}}G^T +  o_{\mathbb P}(1), \quad ( \hat \Sigma_{*, \textsc{l}, \tau x} -  \Sigma_{*, \textsc{l}, \tau x} )  \mid \mathcal{M}_{*} =  o_{\mathbb P}(1), \\
        \hat \Sigma_{*, \textsc{l}, \tau\tau}^{\bot} - \Sigma_{*, \textsc{l}, \tau\tau}^{\bot} \mid \mathcal{M}_{*} &=  GS_{*,\textsc{l}}G^T +  o_{\mathbb P}(1),
        \end{split}
    \]
where $S_{*,\textsc{l}}$ is a positive semi-definite matrix.
\end{theorem}

\begin{comment}
{\color{blue}
 {\tiny\dotfill}

For Theorem 4:
\begin{itemize}
\item helpful to also add results on  $\hat \Sigma_{*, \textsc{l}, \tau\tau}^{\bot} - \Sigma_{*, \textsc{l}, \tau\tau}^{\bot}$?
\item helpful to split the theorem into two? --- asymptotic distribution + covariance estimation for inference?
\end{itemize}
 {\tiny\dotfill}
}
\end{comment}

Theorem~\ref{thm::rerandreg} implies that
the cluster-robust covariance estimator $\hat \Sigma_{*,\textsc{l}, \tau \tau}$
is asymptotically conservative for $\Sigma_{*,\textsc{l}, \tau \tau}$.
As the truncated normal distribution is more peaked than the normal distribution and ${ \hat \Sigma}_{*,\textsc{l}, \tau \tau}$ $ \geq \Sigma_{*,\textsc{l}, \tau \tau} \geq \textnormal{cov}_\textup{a}\{ \sqrt{W}(\hat\tau_{*, \textsc{l}}- \tau) \mid \mathcal{M}_{*}\}$ holds in probability,
we can still use the normal approximation with the cluster-robust covariance to construct Wald-type confidence regions as
$
\{ \tau:  W (\hat \tau_{*,\textsc{l}} - \tau)^\text{T} { \hat \Sigma}_{*,\textsc{l}, \tau\tau}^{-1} (\hat \tau_{*,\textsc{l}} - \tau) \leq \chi^2_{3, 1 - \xi}  \}$.
Such confidence regions, whereas asymptotically valid, are overconservative.
A less conservative confidence region is
$
\{ \tau:  W (\hat \tau_{*,\textsc{l}} - \tau)^\text{T} (\hat \Sigma_{*,\textsc{l},\tau\tau}^{\bot})^{-1} (\hat \tau_{*,\textsc{l}} - \tau) \leq \hat c_{*, \textsc{l}, 1 - \xi}  \},
$
where $\hat c_{*, \textsc{l}, 1 - \xi} $ is defined similarly to $\hat c_{*, 1 - \xi}$, i.e.,  the $1 - \xi$ quantile of $ \phi_{*, \textsc{l}}^\text{T} (\hat \Sigma_{*, \textsc{l}, \tau\tau}^{\bot})^{-1}  \phi_{*, \textsc{l}}$ with $\phi_{*, \textsc{l}} = (\hat \Sigma_{*, \textsc{l},  \tau\tau}^{\bot})^{1/2}\epsilon
            + \hat \Sigma_{*, \textsc{l}, \tau x}\Sigma_{*,xx}^{-1/2}\zeta_{3L, d}$.

Theorem~\ref{thm::rerandreg} extends \citet[][Theorem 1]{Li2018} and \citet[][Theorem 2]{Li2020factorial}
to rerandomization under split-plot designs.
Distinct from these previous results, the asymptotic distributions of the regression-adjusted estimators under split-plot rerandomization are generally not normal, but convolutions of a normal component and a truncated normal component.
The reason is as follows: as shown in the supplementary materials, for $* = \textnormal{ht}, \textnormal{haj}$, the regression adjustments are equivalent to linearly projecting
$\hat {Y}_{*}(z)$ onto $ \hat {v}_{*}(z)$ for $z \in  \mathcal{T}$ separately; however, the separate projection differs from the joint projection of $\hat Y_{*}$ onto $ \hat {v}_{*}$ due to the dependence structure of $\{ \hat {v}_{*}(z)\}_{z \in \mathcal{T}}$, such that $\Sigma_{*,\textsc{l},  \tau x} \neq 0$ in general.
Moreover, the regression-adjusted estimators cannot guarantee efficiency gains over the unadjusted counterparts.

In some special cases, for example, when only whole-plot covariates
 are used with $v_{ws} = \bar v_{w}$ or more generally, $\Psi_{vv} = o(1)$,

 the truncated normal component can disappear and the regression-adjusted Horvitz--Thompson estimator $\hat\tau_{\textnormal{ht}, \textsc{l}}$ is asymptotically more efficient than its unadjusted counterpart under split-plot rerandomization.  Corollary~\ref{cor::varl} below shows the asymptotic distribution of $\hat\tau_{\textnormal{ht}, \textsc{l}} \mid \mathcal{M}_{\textnormal{ht}}$ when $\Psi_{vv} = o(1)$, and ensures its efficiency gain over the unadjusted counterpart

\begin{corollary}\label{cor::varl}
Under Condition~\ref{cond::1}--\ref{cond::3},
if $\Psi_{vv} = o(1)$, then $\Sigma_{\textnormal{ht},\textsc{l}, \tau x} = o(1)$, $\Sigma_{\textnormal{ht}, \textsc{l}, \tau\tau}^{\bot} = \Sigma_{\textnormal{ht}, \textsc{l}, \tau\tau} + o(1)$,
\begin{eqnarray*}
 \sqrt{W}(\hat\tau_{\textnormal{ht}, \textsc{l}}- \tau) \mid \mathcal{M}_{\textnormal{ht}}   \rightsquigarrow    (\Sigma_{\textnormal{ht}, \textsc{l}, \tau\tau}^{\bot})^{1/2}\epsilon, \quad ( { \hat \Sigma}_{\textnormal{ht},\textsc{l}, \tau \tau} -  \Sigma^{\bot}_{\textnormal{ht},\textsc{l}, \tau \tau} )  \mid \mathcal{M}_{\textnormal{ht}} = GS_{\textnormal{ht},\textsc{l}}G^\text{T} +   o_{\mathbb P}(1).
\end{eqnarray*}
Moreover, $\Sigma_{\textnormal{ht}, \tau\tau}^{\bot} \geq \Sigma_{\textnormal{ht}, \textsc{l}, \tau\tau}^{\bot}$ and
\begin{eqnarray*}
    W\left[\textnormal{cov}_\textup{a}(\hat\tau_{\textnormal{ht}} \mid \mathcal{M}_{\textnormal{ht}}) - \textnormal{cov}_\textup{a}(\hat\tau_{\textnormal{ht}, \textsc{l}} \mid \mathcal{M}_{\textnormal{ht}}) \right] &=& \Sigma_{\textnormal{ht}, \tau\tau}^{\bot} - \Sigma_{\textnormal{ht}, \textsc{l}, \tau\tau}^{\bot} + r_{3L, d}\Sigma_{\textnormal{ht}, \tau\tau}^{||} \geq 0.
\end{eqnarray*}
\end{corollary}

Two sufficient conditions for $\Psi_{vv} = o(1)$ are (i) $v_{ws} = \bar v_w$ and (ii) $(S_{w,vv})_{w=1}^W$ are uniformly bounded while $M_w$ goes to infinity for all $w$.
 Corollary~\ref{cor::varl} implies that under either of these two conditions,
 we can ensure efficiency gain of the Horvitz--Thompson estimator by regression with treatment-covariate interactions under split-plot rerandomization.

We cannot guarantee efficiency improvement for the regression-adjusted Hajek estimator under the condition $\Psi_{vv} = o(1)$.
However, when the whole-plot total potential outcomes are more heterogeneous than the whole-plot average potential outcomes,
$\hat \tau_{\textnormal{haj},\textsc{l}}$ can be more efficient than $\hat \tau_{\textnormal{ht},\textsc{l}}$ under their corresponding rerandomization schemes.

\subsection{Covariate adjustment by removing the conditional bias}

By Theorem~\ref{thm::rerandreg}, the regression-adjusted estimators cannot guarantee efficiency gains when heterogeneous unit-level covariates are used in the analysis stage under split-plot randomization or rerandomization. To address this issue, we propose a new covariate-adjusted estimator based on a projection or conditional inference perspective.

Applying Theorem~\ref{thm::clt} to $Y_{ws}(z)$ and $v_{ws}$, $\sqrt{W}( (\hat \tau_{*} - \tau)^\text{T}, \hat \tau^\text{T}_{*,v})^\text{T}$ is asymptotically jointly normal. Then conditional on $\hat \tau_{*,v}$, $\sqrt{W}  ( \hat \tau_{*} - \tau )$ is asymptotically normal with mean $ \sqrt{W} \Sigma_{*,\tau v} \Sigma_{*, vv}^{-1}  \hat \tau_{*,v}$ and covariance $\Sigma_{*,\textsc{p},\tau \tau}^\bot = \Sigma_{*, \tau \tau} - \Sigma_{*, \tau v} \Sigma_{*, vv}^{-1} \Sigma_{*, v \tau} \leq \Sigma_{*, \tau \tau}$.
Let $\hat \Sigma_{*,\tau v} = \hat \Sigma^\text{T}_{*,v \tau} $ be a consistent estimator of $\Sigma_{*,\tau v} = \Sigma^\text{T}_{*,v \tau}$, defined similarly to $\hat \Sigma_{*,\tau x} = \hat \Sigma^\text{T}_{*,x \tau}$ with $x_{ws}$ replaced by $v_{ws}$.
We define
\begin{eqnarray*}
\hat \tau_{*,\textsc{p}} = \hat \tau_{*} - \hat \Sigma_{*,\tau v}  \Sigma_{*, vv}^{-1}  \hat \tau_{*,v}
\end{eqnarray*}
as a conditionally consistent estimator of $\tau$.
Since $W \textnormal{cov}_\textup{a}( \hat \tau_{*,\textsc{p}} )= \Sigma_{*,\textsc{p}, \tau \tau}^\bot = W \min_{\Gamma} \textnormal{cov}_\textup{a}( \hat \tau_{*} - \Gamma  \hat \tau_{*,v} )$, $\hat \tau_{*,\textsc{p}}$
is asymptotically equivalent to the linear projection of $ \hat \tau_{*} $ onto $ \hat \tau_{*, v} $, referred to as the {\it projection estimator} of $\tau$. Let $\hat \Sigma_{*,\textsc{p}, \tau \tau}^\bot = \hat \Sigma_{*, \tau \tau} - \hat \Sigma_{*,\tau v} \Sigma_{*, vv}^{-1} \hat \Sigma_{*, v \tau}$.

\begin{theorem}\label{thm::project}
Under Conditions~\ref{cond::1}--\ref{cond::3}, for $* = \textnormal{ht}, \textnormal{haj}$,
\begin{eqnarray*}
 \sqrt{W}(\hat\tau_{*, \textsc{p}}- \tau) \mid \mathcal{M}_{*} &\rightsquigarrow &
 (\Sigma_{*,\textsc{p}, \tau \tau}^\bot)^{1/2}\epsilon, \quad ( \hat \Sigma_{*,\textsc{p}, \tau \tau}^\bot -  \Sigma_{*,\textsc{p}, \tau \tau}^\bot )  \mid \mathcal{M}_{*} = GS_{*}G^\text{T} +   o_{\mathbb P}(1) .
\end{eqnarray*
Moreover,
\begin{eqnarray*}
    W\left[\textnormal{cov}_\textup{a}(\hat\tau_{*} \mid \mathcal{M}_{*}) - \textnormal{cov}_\textup{a}(\hat\tau_{*,\textsc{p}} \mid \mathcal{M}_{*}) \right] =
    \Sigma_{*,\tau \tau}^{\bot} - \Sigma_{*,\textsc{p},\tau\tau}^{\bot} + r_{3L, d} \Sigma_{*,\tau \tau}^{||} \geq 0.
\end{eqnarray*}
\end{theorem}

Theorem~\ref{thm::project} implies that, under the rerandomization scheme $\mathcal{M}_{*}$, the treatment effect estimator $\hat\tau_{*, \textsc{p}}$ is consistent and asymptotically normal, and the covariance estimator $\hat \Sigma_{*,\textsc{p}, \tau \tau}^\bot$ is asymptotically conservative. Moreover, $\hat\tau_{*, \textsc{p}}$ improves the efficiency of $\hat \tau_{*}$ without requiring $\Psi_{vv} = o(1)$.
Based on this theorem, an asymptotically conservative Wald-type confidence region for $\tau$ is
$
\{ \tau:  W (\hat \tau_{*,\textsc{p}} - \tau)^\text{T} (\hat \Sigma_{*,\textsc{p}, \tau \tau}^\bot)^{-1} (\hat \tau_{*,\textsc{p}} - \tau) \leq \chi^2_{3, 1 - \xi}  \}.
$

\subsection{Relative efficiency of different rerandomization and estimation schemes}

We have introduced the regression-adjusted and projection-based variants for both the Horvitz--Thompson and Hajek estimators of the average treatment effects.
Corollary~\ref{cor::comrr} below gives the relative efficiency between the Horvitz--Thompson and Hajek estimators either with or without covariate adjustment under their respective rerandomization schemes.

Let
$Q_{\textnormal{in},vv} = (N-1)^{-1} \sum_{ws \in \mathcal{S}} (v_{ws}-\bar v_w)(v_{ws}-\bar v_w)^\text{T}$ be a variant of $Q_{vv}$ with $v_{ws}$ centered by the whole-plot average $\bar v_w$ instead of $\bar v$.
It is then a weighted average of the $S_{w,vv}$'s with $Q_{\textnormal{in},vv} = (N-1)^{-1}\sum_{w=1}^{W} (M_w-1) \alpha_w^{-2} S_{w,vv}$.
Similarly define $Q_{\textnormal{in},vY(z)}=(N-1)^{-1} \sum_{ws \in \mathcal{S}}  (v_{ws}-\bar v_w)\left\{ Y_{ws}(z)-\bar Y_w (z)\right\}$ and $Q_{\textnormal{in}} (z,z^\prime) = (N-1)^{-1} \sum_{ws \in \mathcal{S}}  \left\{ Y_{ws}(z)-\bar Y_w (z)\right\} \left\{ Y_{ws}(z^\prime)-\bar Y_w (z^\prime)\right\}$
for $z, z^\prime \in \mathcal{T}$. We use the subscript ``in" to signify within whole-plot covariances.

\begin{condition}\label{cond::Qcond}
As $W$ goes to infinity, $Q_{\textnormal{in},vv} = o(1)$ and $Q_{\textnormal{in}} (z,z)  = O(1)$ for all $z \in \mathcal{T}$.
\end{condition}

\begin{remark}
If only whole-plot covariates are used, then $Q_{\textnormal{in},vv} = 0$ and $\Psi_{vv} =0$. Both $Q_{\textnormal{in},vv}$ and $\Psi_{vv}$ measure the variability of covariates within whole plots, but $Q_{\textnormal{in},vv}=o(1)$ is a stricter condition than $\Psi_{vv} = o(1)$. See the supplementary materials for details.

\end{remark}

\begin{corollary}\label{cor::comrr}
Under Conditions~\ref{cond::1}--\ref{cond::3},
\begin{itemize}
    \item[\textup{(i)}]
    $\textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj}} \mid
    \mathcal{M}_{\textnormal{haj}}) = \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht}}  \mid
    \mathcal{M}_{\textnormal{ht}})$ and $\textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj},\textsc{p}}\mid \mathcal{M}_{\textnormal{haj}}) = \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht},\textsc{p}}\mid \mathcal{M}_{\textnormal{ht}})$
    if $\bar x = 0$, $\bar Y(z)=0$ for all $z$ or $\alphaw=1$ for all $w$;\\
    $\textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj},\textsc{l}}\mid \mathcal{M}_{\textnormal{haj}}) = \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht},\textsc{l}}\mid \mathcal{M}_{\textnormal{ht}}) $ if the design is uniform

    and Condition~\ref{cond::Qcond} holds;
    \item[\textup{(ii)}] Further assume that $\Psi_{vv} = o(1)$, then
\begin{eqnarray*}
\textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj}} \mid
    \mathcal{M}_{\textnormal{haj}}) \leq \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht}}  \mid
    \mathcal{M}_{\textnormal{ht}}), \quad \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj},\diamond}\mid \mathcal{M}_{\textnormal{haj}}) \leq \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht},\diamond}\mid \mathcal{M}_{\textnormal{ht}})\quad (\diamond= \textsc{l}, \textsc{p})
\end{eqnarray*} if $\bar Y_w(z)$ are constant over all $w$, and
 \begin{eqnarray*}
 \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj}} \mid
    \mathcal{M}_{\textnormal{haj}}) \geq \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht}}  \mid
    \mathcal{M}_{\textnormal{ht}}), \quad \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{haj},\diamond}\mid \mathcal{M}_{\textnormal{haj}}) \geq \textnormal{cov}_\textup{a}(\hat \tau_{\textnormal{ht},\diamond}\mid \mathcal{M}_{\textnormal{ht}})\quad (\diamond = \textsc{l}, \textsc{p})
    \end{eqnarray*}
     if $\alphaw \bar Y_w(z)$ are constant over all $w$.
\end{itemize}
\end{corollary}

Corollary~\ref{cor::comrr}{(i)} implies that $\hat \tau_{\textnormal{haj},\textsc{p}}$ and $\hat \tau_{\textnormal{ht},\textsc{p}}$ are asymptotically equally efficient if the whole plots are of equal sizes, and $\hat \tau_{\textnormal{haj},\textsc{l}}$ and $\hat \tau_{\textnormal{ht},\textsc{l}}$ are asymptotically equally efficient under uniform design and Condition~\ref{cond::Qcond}.
Suppose that the within whole-plot covariance of covariates is neglectable, i.e., $\Psi_{vv} = o(1)$. Corollary~\ref{cor::comrr}(ii) implies that, under split-plot rerandomization,
$\hat \tau_{\textnormal{haj}}$, $\hat \tau_{\textnormal{haj},\textsc{p}}$, and $\hat \tau_{\textnormal{haj},\textsc{l}}$ are asymptotically more efficient than $\hat \tau_{\textnormal{ht}}$, $\hat \tau_{\textnormal{ht},\textsc{p}}$ and $\hat \tau_{\textnormal{ht},\textsc{l}}$, respectively, if the whole plots have similar average potential outcomes, and
vice versa if the whole plots have similar total potential outcomes.
As the whole-plot totals are often more heterogeneous than the whole-plot averages in practice, we prefer the Hajek estimators and the associated rerandomization scheme over the Horvitz--Thompson estimators and the associated rerandomization scheme in general.

Next, we study the relative efficiency of the regression-adjusted estimators versus the projection estimators.
Let $\hat \tau_{\textnormal{ht},\textsc{l},\alpha}$ denote the analog of $\hat \tau_{\textnormal{ht},\textsc{l}}$ that further includes the centered whole-plot size factor $\alphaw - 1$  as an additional covariate in the regression formula \eqref{eq:ag_v}.

\begin{corollary}\label{cor::comrrALL}
Under Conditions~\ref{cond::1}--\ref{cond::3},
if $\Psi_{vv} = o(1)$, then $\hat \tau_{\textnormal{ht},\textsc{l},\alpha} \mid \mathcal{M}_{\textnormal{ht}}$ is the most peaked around $\tau$ among the set of estimators:
\begin{eqnarray*}
  \{(\hat \tau_{\textnormal{ht},\textsc{l},
  \alpha} \mid \mathcal{M}_{\textnormal{ht}}), \ \hat \tau_{*} , \ (\hat \tau_{*} \mid \mathcal{M}_{*}), \ \hat \tau_{*,\textsc{l}},\ (\hat \tau_{*,\textsc{l}} \mid \mathcal{M}_{*}), \ \hat \tau_{*,\textsc{p}},\ (\hat \tau_{*,\textsc{p}} \mid \mathcal{M}_{*}): * = \textnormal{ht}, \textnormal{haj} \}.
\end{eqnarray*}
\end{corollary}

Corollary~\ref{cor::comrrALL} establishes the optimality of $\hat \tau_{\textnormal{ht},\textsc{l},\alpha} \mid \mathcal{M}_{\textnormal{ht}}$  among all considered estimators when $\Psi_{vv} = o(1)$, and highlights the utility of including $\alphaw -1$ as an additional covariate in the aggregate regression for ensuring additional efficiency. Intuitively, the unadjusted Hajek estimator $\hat\tau_{\textnormal{haj}}$ implicitly adjusts for the whole-plot sizes, and is hence in general better than the unadjusted Horvitz--Thompson estimator $\hat\tau_{\textnormal{ht}}$; see the comments after Corollary \ref{cor::comrr}.
The $\hat\tau_{\textnormal{ht}, \textsc{l},\alpha}$, on the other hand, gives a more efficient way of adjusting for the whole-plot sizes than the Hajek estimator when $\Phi_{vv} = o(1)$.
  We thus recommend the split-plot rerandomization scheme $\mathcal{M}_{\textnormal{ht}}$ and the associated regression-adjusted estimator $\hat \tau_{\textnormal{ht},\textsc{l},\alpha}$  when the covariates are relatively homogeneous within whole plots or when only whole-plot covariates are used. When the covariates vary greatly within whole plots such that $\Phi_{vv} = o(1)$ does not hold, the projection estimators $\hat \tau_{*,\textsc{p}} \ (* = \textnormal{haj},\textnormal{ht})$ always improve the efficiency under rerandomization, whereas the regression-adjusted estimators $\hat \tau_{\textnormal{ht},\textsc{l},\alpha}$ and $\hat \tau_{*,\textsc{l}}$ may degrade efficiency compared to the unadjusted counterparts.
This gives an advantage of projection adjustment over regression adjustment. We illustrate this by simulation.

\section{Numerical Examples}

\subsection{Simulation}
In this section, we conduct simulation to assess the finite-sample performance of the unadjusted and covariate-adjusted estimators under split-plot rerandomization.
We set $W=600$, $(W_1, W_0) = (0.3W, 0.7W)$,  and generate $(M_{w0}, M_{w1}, M_{w})_{w=1}^W$ as $M_{w0} = \max(2,\zeta_{w0})$, $M_{w1} =\max(2,\zeta_{w1})$, and $M_w = M_{w0} + M_{w1}$, where $\zeta_{w0}$'s are independent Poisson(5) and $\zeta_{w1}$'s are independent Poisson(3).
For $w=1, \ldots, W$, we draw $v_w = (v_{w1}, v_{w2})^\text{T}$ independently from $\mathcal N((0.6,0.6)^\text{T}, 0.8I_2)$, and use the following two methods to construct subplot covariates {$v_{ws} = (v_{ws,1}, v_{ws,2})^\text{T}$}: (i) $v_{ws}=v_w$ for $s = 1,\ldots,M_w$, which
corresponds to the case where only whole-plot covariates are used and ensures $\Psi_{vv}= o(1)$;
(ii) $v_{ws}=v_w + \delta_{ws}$ for $ws \in \mathcal{S}$, where $\delta_{ws}$'s are independent $\mathcal{N} ( 0_2, 0.5I_2 )$, so that the covariates vary within each whole plot.
{We use $v_{ws}$ for covariate adjustment in the analysis stage, and set $x_{ws} = v_{ws,1}$ for rerandomization in the design stage.}
The potential outcomes are then generated as
    \[
        \begin{split}
            Y_{ws}(00) &= \theta_w + 0.5 + 2v_{ws,1}^2 + 2v_{ws,2}^2 +\epsilon_{ws},\\
            Y_{ws}(01) &= -0.5\theta_w + 1 + v_{ws,1}^2 + v_{ws,2}^2 +\epsilon_{ws},\\
            Y_{ws}(10) &= 0.5\theta_w + 1 - v_{ws,1}^2 - v_{ws,2}^2 +\epsilon_{ws},\\
            Y_{ws}(11) &= \theta_w + 2 + 2v_{ws,1}^2 + 2v_{ws,1}^2 +\epsilon_{ws}
        \end{split}
    \]
for $ws \in \mathcal{S}$, where $\theta_w$'s are independent $\mathcal{N} ( 2\text{max}(M_w)/M_w, 0.2 )$ and $\epsilon_{ws}$'s are independent Uniform$(-1, 1)$.
The covariates and potential outcomes are generated once and then kept fixed.
We perform $2^2$ split-plot randomization and two types of split-plot rerandomization 2,000 times, respectively, and summarize the operating characteristics of $\htau_*$, $\hat\tau_{*, \textsc{l}}$,  $\hat\tau_{*, \textsc{p}}$, and $\hat\tau_{\textnormal{ht},\textsc{l},\alpha}$ for $* = \textnormal{ht}, \textnormal{haj}$.
For rerandomization criteria, we set $d$ to be the 1st percentile
of $\chi_3^2$, implying an asymptotic acceptance rate of $1\%$.

Figure~\ref{Simulation1} shows the comparison between estimators under  split-plot randomization and rerandomization when only whole-plot covariates are used.
The first row illustrates the biases of the covariate-adjusted estimators in finite samples.
These estimators are asymptotically unbiased, but can have small finite-sample biases \citep{lin2013}.
The second row shows the standard deviations, illustrating the efficiency gain by rerandomization and covariate adjustment.
Among them, $\hat\tau_{\textnormal{ht},\textsc{l},\alpha}$ under rerandomization is the most efficient, which is coherent with the result of Corollary~\ref{cor::comrrALL}.
The third row shows the positive empirical biases of standard deviation estimators, implying the conservativeness of distribution estimation.
The fourth row shows the coverage rates of the constructed $95\%$ confidence intervals, and suggests the validity of all estimators under rerandomization.
The fifth row shows the average confidence interval lengths, which illustrates the efficiency gain by conducting inference with both rerandomization and  covariate adjustment.

Figure \ref{Simulation2} shows the analogous results when covariates vary within each whole plot.
We can see that $\hat\tau_{\textnormal{ht},\textsc{l},\alpha}$ is no longer the most efficient, as $\Psi_{vv} = o(1)$ is not satisfied.
In this case, the projection estimators $\hat \tau_{*,\textsc{p}} \ (* = \textnormal{haj},\textnormal{ht})$ always improve the efficiency, but the regression-adjusted estimators $\hat \tau_{*,\textsc{l}}$ may degrade efficiency compared to the unadjusted estimator under  rerandomization. We present an example in the supplementary materials.

\begin{figure}
\includegraphics[width=\linewidth]{RRSP_sim1.pdf}
\caption{
Comparison of the estimators under $2^2$ split-plot randomization and rerandomization with $v_{ws} = v_{w}$.
The row ``bias" summarizes the average deviations of the point estimators from the true values. The row ``sd" summarizes the standard deviations of the point estimators.
The row ``ese" summarizes the average errors of the standard deviation estimators. The row ``coverage" summarizes the coverage rates of the 95\% confidence intervals. The row ``length" summarizes the average confidence interval lengths of 95\% confidence intervals. The column ``ht.rnd" stands for $\hat\tau_\textnormal{ht}$ under classic split-plot randomization, ``ht" stands for $\hat\tau_{\textnormal{ht}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.P" stands for $\hat\tau_{\textnormal{ht},\textsc{p}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L" stands for $\hat\tau_{\textnormal{ht},\textsc{l}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L.a" stands for $\hat\tau_{\textnormal{ht},\textsc{l},\alpha} \mid \mathcal{M}_{\textnormal{ht}}$, ``haj.rnd" stands for $\hat\tau_{\textnormal{haj}}$ under classic split-plot randomization, ``haj" stands for $\hat\tau_{\textnormal{haj}} \mid \mathcal{M}_{\textnormal{haj}}$, ``haj.P" stands for $\hat\tau_{\textnormal{haj},\textsc{p}} \mid \mathcal{M}_{\textnormal{haj}}$, and ``haj.L" stands for $\hat\tau_{\textnormal{haj},\textsc{l}} \mid \mathcal{M}_{\textnormal{haj}}$.
\label{Simulation1}}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{RRSP_sim2.pdf}
\caption{Comparison of the estimators under $2^2$ split-plot randomization and rerandomization with varying $v_{ws}$ within each whole plot. The row ``bias" summarizes the average deviations of the point estimators from the true values. The row ``sd" summarizes the standard deviations of the point estimators. The row ``ese" summarizes the average errors of the standard deviation estimators. The row ``coverage" summarizes the coverage rates of the 95\% confidence intervals. The row ``length" summarizes the average confidence interval lengths of 95\% confidence intervals. The column ``ht.rnd" stands for $\hat\tau_\textnormal{ht}$ under classic split-plot randomization, ``ht" stands for $\hat\tau_{\textnormal{ht}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.P" stands for $\hat\tau_{\textnormal{ht},\textsc{p}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L" stands for $\hat\tau_{\textnormal{ht},\textsc{l}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L.a" stands for $\hat\tau_{\textnormal{ht},\textsc{l},\alpha} \mid \mathcal{M}_{\textnormal{ht}}$, ``haj.rnd" stands for $\hat\tau_{\textnormal{haj}}$ under classic split-plot randomization, ``haj" stands for $\hat\tau_{\textnormal{haj}} \mid \mathcal{M}_{\textnormal{haj}}$, ``haj.P" stands for $\hat\tau_{\textnormal{haj},\textsc{p}} \mid \mathcal{M}_{\textnormal{haj}}$, and ``haj.L" stands for $\hat\tau_{\textnormal{haj},\textsc{l}} \mid \mathcal{M}_{\textnormal{haj}}$.
\label{Simulation2}}
\end{figure}

\subsection{Real data illustration}
In this section, we analyze a real data set to assess the performance of different estimators under split-plot randomization and rerandomization. \cite{olken2007monitoring} conducted a randomized field experiment on reducing corruption in 608 Indonesian village road projects. We consider two interventions of the study: increasing the probability of external government audits (``audits") and distributing invitations to accountability meetings (``invitations").

The villages are nested in subdistricts, and there was a concern of the spillover effect of audits.
Therefore, the randomization of audits was clustered by subdistrict. On the other hand, the randomization of invitations was done village by village. This defines a  nonuniform split-plot experiment with the audits and invitations constituting the whole-plot and subplot factors, respectively.

Before the experiment, \cite{olken2007monitoring} conducted a survey to collect
ten village characteristics, including village population, village head education, village head salary, and total budget.
To measure the corruption level as the primary outcome of interest, \cite{olken2007monitoring} constructed an independent estimate of the amount each project actually cost to build and then compared it with what the village reported it spent on the project. The percent missing, defined as the difference between the log of the reported amount and the log of the actual amount, is the main measure of corruption level used in the experiment.

We fill in the missing potential outcomes before the analysis.
In the data set, there are subdistricts containing only one village, for which we can not calculate covariances within whole plot such as $S_w$ and $S_{w,xY}$. We leave out those subdistricts, and there are 136 subdistricts and 550 villages left. The missing potential outcomes are filled by linear regression based on treatments and ten covariates.

In our analysis, we include village population and village head salary as covariates used for both rerandomization ($x$)
and
covariate adjustment ($v$), and focus on the missing percent for materials in road project as the outcome.
We then perform $2^2$ split-plot randomization and two types of split-plot rerandomization 1,000 times, respectively. For rerandomization criteria, we set $d$ to be the  1st percentile
of $\chi_6^2$.

Figure~\ref{Olken1} shows the results. Here we use relative standard deviation and average confidence interval length compared to the Horvitz--Thompson estimator under classic split-plot randomization to display the results more clearly. From the second and fifth rows we can see that rerandomization gains estimation and inference efficiency for both the Horvitz--Thompson and Hajek estimators. For example, the standard deviation and average confidence interval length of the Horvitz--Thompson estimator are reduced by rerandomization by approximately $8\%$ for the main effect of audits. Overall, the Hajek estimator performs better than the Horvitz--Thompson estimator. This may be because the
subdistricts, as the whole plots in our example, have similar average potential outcomes. The covariate-adjusted estimators $\hat \tau_{\textnormal{ht},\textsc{l}, \alpha}$, $\hat \tau_{\textnormal{haj},\textsc{l}}$, and $\hat \tau_{\textnormal{haj},\textsc{p}}$ perform similarly and are the best methods.

\begin{figure}
\includegraphics[width=\linewidth]{RRSP_Olken.pdf}
\caption{ Comparison of different estimators using \cite{olken2007monitoring}'s data set. Factor A and factor B represent audits and invitations, respectively.
The row ``bias" summarizes the average deviations of the point estimators from the true values. The row ``ese" summarizes the average errors of the standard deviation estimators. The row ``coverage" summarizes the coverage rates of the 95\% confidence intervals. The row ``sd (relative)" and ``length (relative)" summarizes the standard deviations and average confidence interval lengths of 95\% confidence intervals divided by that of Horvitz--Thompson estimator under classic $2^2$ split-plot randomization (``ht.rnd").  The numbers in parentheses are the absolute values for standard deviations and interval lengths. The column ``ht.rnd" stands for $\hat\tau_\textnormal{ht}$ under classic split-plot randomization, ``ht" stands for $\hat\tau_{\textnormal{ht}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.P" stands for $\hat\tau_{\textnormal{ht},\textsc{p}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L" stands for $\hat\tau_{\textnormal{ht},\textsc{l}} \mid \mathcal{M}_{\textnormal{ht}}$, ``ht.L.a" stands for $\hat\tau_{\textnormal{ht},\textsc{l},\alpha} \mid \mathcal{M}_{\textnormal{ht}}$, ``haj.rnd" stands for $\hat\tau_{\textnormal{haj}}$ under classic split-plot randomization, ``haj" stands for $\hat\tau_{\textnormal{haj}} \mid \mathcal{M}_{\textnormal{haj}}$, ``haj.P" stands for $\hat\tau_{\textnormal{haj},\textsc{p}} \mid \mathcal{M}_{\textnormal{haj}}$, ``haj.L" stands for $\hat\tau_{\textnormal{haj},\textsc{l}} \mid \mathcal{M}_{\textnormal{haj}}$.
\label{Olken1}}
\end{figure}

\section{Discussion}
We investigated the asymptotic properties of rerandomization and covariate adjustment under split-plot designs.
Based on the asymptotic results, we recommend the use of  rerandomization scheme based on the Horvitz--Thompson estimator if the whole plots have similar total potential outcomes, and  rerandomization scheme based on the Hajek estimator if the whole plots have similar average potential outcomes.
In the analysis stage, we recommend the fully interacted aggregate regression after adjusting for the whole-plot sizes if only whole-plot covariates are used or more generally,
$\Psi_{vv} = o(1)$,
and the projection estimator otherwise. The resulting inference is model-free, and remains valid regardless of how well the regression specifications represent the true data generating process of the outcome, treatments, and covariates.

\bigskip
\begin{center}
{\large\bf Supplementary Materials}
\end{center}

\begin{description}
\item The supplementary materials provide additional simulation results and proofs.
\end{description}

\bibliographystyle{agsm}