# Bayesian criterion for Re-randomization

Zhaoyang Liu∗ Center for Statistical Science & Department of Industrial
Engineering, Tsinghua University  Tingxuan Han111These authors have
contributed equally to this work. Center for Statistical Science & Department
of Industrial Engineering, Tsinghua University  Donald B. Rubin  Yau
Mathematical Sciences Center, Tsinghua University  Ke Deng222Corresponding
author<EMAIL_ADDRESS>Center for Statistical Science & Department of
Industrial Engineering, Tsinghua University

###### Abstract

Re-randomization has gained popularity as a tool for experiment-based causal
inference due to its superior covariate balance and statistical efficiency
compared to classic randomized experiments. However, the basic re-
randomization method, known as ReM, and many of its extensions have been
deemed sub-optimal as they fail to prioritize covariates that are more
strongly associated with potential outcomes. To address this limitation and
design more efficient re-randomization procedures, a more precise
quantification of covariate heterogeneity and its impact on the causal effect
estimator is in a great appeal. This work fills in this gap with a Bayesian
criterion for re-randomization and a series of novel re-randomization
procedures derived under such a criterion. Both theoretical analyses and
numerical studies show that the proposed re-randomization procedures under the
Bayesian criterion outperform existing ReM-based procedures significantly in
effectively balancing covariates and precisely estimating the unknown causal
effect.

##  1 Introduction

Covariate imbalance is to be avoided in randomized experiments, even though
complete randomization balances all potential confounding factors on average
and therefore becomes the “gold standard” in causal inference (Holschuh,,
1980; Wu,, 1981; Urbach,, 1985; Imai et al.,, 2008; Cox,, 2009; Morgan and
Rubin,, 2012). A natural way to avoid treatment assignment allocations with
unbalanced covariates is to reject a “bad” allocation, and redo the
randomization until a “good” allocation with properly balanced covariates is
obtained, before the experiment is conducted. This is called re-randomization,
which was first suggested by Fisher, (1992). The general framework including
the sufficient condition to guarantee unbiased estimation and related benefits
from doing re-randomization when using Mahalanobis distance as the metric to
measure imbalance (referred to as ReM) is provided in Morgan and Rubin,
(2012).

Subsequent to the initial work, relevant literatures have shown extended
applications and theoretical properties of ReM. Unequal weights on covariates
are considered using _ReM in tiers of covariates_ (ReMT) in Morgan and Rubin,
(2015), where covariates are partitioned into tiers according to their
relative importance, with a tighter threshold of ReM for covariates in a more
important tier. After re-randomization in one tier, residuals orthogonal to
all previous tiers of covariates can be used to calculate Mahalanobis distance
in the next tier in ReMT, leading to reduced computational complexity and
refined statistical efficiency of ReM when a large number of covariates are
involved. Applications of ReM in factorial and sequential designs are analyzed
in Branson et al., (2016) and Zhou et al., (2018). Asymptotic theory for the
standard difference-in-means estimator of treatment effect under ReM is
derived by an orthogonal decomposition of the final sampling distribution as a
linear combination of a Gaussian random variable and a truncated Gaussian
random variable (Li et al.,, 2018, 2020). The asymptotic performance of ReM in
stratified randomization and survey experiments are studied in Wang et al.,
(2021) and Yang et al., (2021). Zhang et al., 2021b  recommended applying ReM
for the top-$k$ principal components of covariates only, leading to a variate
of ReM named PCA-ReM, which reduces computational complexity of ReM
significantly while retaining much of its benefits. Ridge-ReM proposed by
Branson and Shao, (2021) utilizes modified Mahalanobis distance to deal with
collinearities among covariates in situations with high-dimensional or highly-
correlated covariates. Considering that a small acceptance probability will
lead to heavy computational burden, Zhu and Liu, (2022) proposed a pair-
switching strategy to reduce the computational cost of ReM.

Beyond the ReM-based approaches, alternative strategies for re-randomization
based on other criteria have been proposed as well. Kallus, (2018) proposed an
optimal allocation based on a minimax criterion, which remains controversial,
as discussed in Johansson et al., (2021) and Kallus, (2021). To tackle with
this problem, Wang and Li, (2022) proposed a procedure allowing the covariate
imbalance to diminish as sample size increases. On the other hand, Li et al.,
(2021) suggested to estimate the density function of covariates in the
treatment group and the control group via kernel density estimation, and
measure the degree of covariate imbalance by comparing the estimated density
functions.

However, all these existing re-randomization methods are sub-optimal in sense
that they fail to prioritize covariates that are more strongly associated with
potential outcomes, while there is a clear insight that such covariates should
enjoy higher priority to be balanced. Although ReMT considers covariate
heterogeneity by partitioning the covariates into tiers and assigning higher
priority to balancing covariates in the leading tiers, such a solution is
incomplete because it treats covariates in the same tier equally. A more
precise quantification of covariate heterogeneity and its impact on the causal
effect estimator would lead to a better solution.

This study fills in this gap by establishing a Bayesian criterion for re-
randomization, which formulates our knowledge on the relative importance of
different covariates to the potential outcomes via a prior distribution, and
derives a novel re-randomization procedure referred to as ReB. We found that
ReB takes many existing ReM-based re-randomization procedures, e.g., ReM, PCA-
ReM and Ridge-ReM, as special cases, forming a unified framework for studying
various re-randomization approaches. Theoretical analyses show that the
proposed ReB outperforms the ReM-based re-randomization procedures in terms of
getting a more accurate causal effect estimator, as long as the prior
distribution is informative to highlight the relative importance of different
covariates.

In case that no such a nice informative prior distribution is available, we
propose a two-stage strategy to implement ReB, where a small pilot experiment
is conducted in the first stage to establish the informative prior
distribution we need, and ReB is conducted for the remaining samples in the
second stage. The final causal effect estimator can be obtained by properly
integrating the two estimators obtained in both stages. Theoretical analyses
and simulation studies show that the two-stage ReB procedure also achieves
superior covariate balance and thereby more precise estimation of causal
effects than existing methods. In the literature, the idea of using pre-
experimental data in re-randomaization has been studied. For example,
Johansson and Schultzberg, (2020) proposed a rank-based balance measure for
re-randomization and weight for each covariates was estimated using pre-
experimental data; Zhang et al., 2021a  introduced the response-adaptive
design for re-randomization and took ethics into consideration. Compared to
the proposed two-stage ReB procedure, however, these two methods do not make
full use of the information obtained from the pre-experiment and suffer from
the lack of theoretical guarantees.

The following part of this paper is organized as follows. Section 2 defines
the basic notations and briefly reviews the repeated sampling properties of
re-randomization criterion based on the Mahalanobis distance in the
literature. Section 3 proposes the Bayesian criterion for re-randomization,
establishes the corresponding theoretical properties and discusses its
connections to other re-randomization criteria in the literature. Section 4
introduces the framework of two-stage Bayesian re-randomization procedures and
gives their asymptotic properties. Section 5 and 6 support the proposed
methodology with simulation studies and a real data application. Finally, we
conclude with a discussion in Section 7. Detailed proofs and additional
simulation results are provided in the Supplementary Material.

##  2 Notation and preliminaries

Following the classic setting for causal inference in Rubin, (1974) and in
Imbens and Rubin, (2015), we consider statistical inference for the _average
causal effect_ (ACE) of an _active treatment_ $t$ with respect to a _control
treatment_ $c$ over a finite population of units
${\mathcal{P}}=\\{1,\cdots,N\\}$, i.e.,

$\tau=\frac{1}{N}\sum_{i=1}^{N}\tau_{i}=\frac{1}{N}\sum_{i=1}^{N}\big{(}Y_{i}(1)-Y_{i}(0)\big{)},$

where $Y_{i}(1)$ and $Y_{i}(0)$ are the potential outcomes of unit
$i\in{\mathcal{P}}$ under treatment $t$ and control $c$ respectively, and
$\tau_{i}=Y_{i}(1)-Y_{i}(0)$ is the unit-level causal effect.

Let $W_{i}$ be the _treatment assignment indicator_ of unit $i$ ($W_{i}=1$ if
unit $i$ is assigned to receive the active treatment $t$, and $W_{i}=0$ if the
control treatment $c$ is assigned), and
${\boldsymbol{W}}=(W_{1},\cdots,W_{N})^{T}$ be the _treatment assignment
vector_ for the experiment. In a randomized experiment with $N_{t}$ and
$N_{c}$ being the group size of the treatment group and control group
respectively (apparently $N_{t}+N_{c}=N$), the treatment assignment vector is
a random sample from the space:

$\begin{split}{\mathscr{W}}=\left\\{{\boldsymbol{W}}\in\\{0,1\\}^{N}:\sum_{i=1}^{N}W_{i}=N_{t}\right\\},\end{split}$

leading to the following unbiased estimator of the causal estimand $\tau$
under the Stable Unit Treatment Value Assumption (SUTVA) introduced by Rubin,
(1980):

$\hat{\tau}({\boldsymbol{W}})=\frac{1}{N_{t}}\sum_{i=1}^{N}Y_{i}(1)\cdot
W_{i}-\frac{1}{N_{c}}\sum_{i=1}^{N}Y_{i}(0)\cdot(1-W_{i}).$

###  2.1 Framework of re-randomization

Further assume that unit $i$ is associated with a group of covariates
${\boldsymbol{X}}_{i}=(X_{i1},\cdots,X_{ip})^{T}$, which may have prediction
power for potential outcomes $(Y_{i}(1),Y_{i}(0))$, and thus the causal effect
$\tau_{i}$. Let ${\boldsymbol{X}}=(X_{ij})_{1\leq i\leq N,1\leq j\leq p}$ be
the covariate matrix of the $N$ units in the population of interest. A re-
randomization procedure, which can be traced back to R.A. Fisher as pointed
out by Morgan and Rubin, (2012), defines a deterministic _acceptance rule_

$\phi(\boldsymbol{X},\boldsymbol{W})=\begin{cases}1,&\mbox{if assignment
${\boldsymbol{W}}$ is acceptable,}\\\ 0,&\mbox{otherwise,}\end{cases}$

to reject “bad” assignments (i.e., ones making
$\phi({\boldsymbol{X}},{\boldsymbol{W}})=0$), until a plausible assignment
${\boldsymbol{W}}$ satisfying $\phi({\boldsymbol{X}},{\boldsymbol{W}})=1$ is
obtained. Acceptance rule $\phi({\boldsymbol{X}},{\boldsymbol{W}})$ defines a
subset of acceptable assignments ${\mathscr{W}}_{\phi}\subseteq{\mathscr{W}}$,
which is referred to as the _acceptance region_ of rule $\phi$.

Morgan and Rubin, (2012) suggested specifying the re-randomization procedure
$\phi$ by controlling the Mahalanobis distance between

$\bar{{\boldsymbol{X}}}_{t}=\frac{1}{N_{t}}\sum_{i=1}^{N}{\boldsymbol{X}}_{i}W_{i}\quad\mbox{and}\quad\bar{{\boldsymbol{X}}}_{c}=\frac{1}{N_{c}}\sum_{i=1}^{N}{\boldsymbol{X}}_{i}(1-W_{i}),$

which is defined as

$d_{M}\triangleq{\boldsymbol{D}}^{T}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{-1}{\boldsymbol{D}}$

with ${\boldsymbol{D}}=\bar{{\boldsymbol{X}}}_{t}-\bar{{\boldsymbol{X}}}_{c}$
and
${\boldsymbol{\Sigma}}_{\boldsymbol{D}}=\operatorname{cov}\left({\boldsymbol{D}}\right)$
being the covariance matrix of ${\boldsymbol{D}}$, leading to the following
_re-randomization using the Mahalanobis distance_ (ReM):

$\phi_{M}({\boldsymbol{X}},{\boldsymbol{Z}})=I\big{(}{\boldsymbol{D}}^{T}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{-1}{\boldsymbol{D}}\leq
a\big{)},$

where $I(\cdot)$ is the 0-1 indicator function, and the threshold $a$ can be
determined based on the equation below for each specific acceptance rate
$\alpha\in(0,1)$:

${\mathbb{P}}(\phi_{M}=1)={\mathbb{P}}\big{(}{\boldsymbol{D}}^{T}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{-1}{\boldsymbol{D}}\leq
a\big{)}=\alpha.$

Denote $(\hat{\tau}\mid\phi_{M}=1)$ as the difference-in-mean estimator of
$\tau$ under the re-randomization mechanism $\phi_{M}$. Morgan and Rubin,
(2012) showed that $(\hat{\tau}\mid\phi_{M}=1)$ is an unbiased estimator of
$\tau$ satisfying ${\mathbb{E}}(\hat{\tau}\mid\phi_{M}=1)=\tau$ as long as
$N_{t}=N_{c}$; and, $(\hat{\tau}\mid\phi_{M}=1)$ is a more efficient estimator
of $\tau$ that the classic $\hat{\tau}$ under a completely randomized
experiment with the following _percent reduction in variance_ (PRIV):

$PRIV=1-\frac{{\text{var}}(\hat{\tau}\mid\phi_{M}=1)}{{\text{var}}(\hat{\tau})}=100\times(1-v_{\alpha,p})R^{2},$

where $R^{2}$ represents the squared multiple correlation between the observed
potential outcomes and the covariates within each treatment group, and

$v_{\alpha,p}\triangleq\frac{{\mathbb{E}}(\chi^{2}_{p}\mid\chi^{2}_{p}\leq\xi_{\alpha,p})}{p}=\frac{{\mathbb{P}}(\chi^{2}_{p+2}\leq\xi_{\alpha,p})}{{\mathbb{P}}(\chi^{2}_{p}\leq\xi_{\alpha,p})}=\frac{{\mathbb{P}}(\chi^{2}_{p+2}\leq\xi_{\alpha,p})}{\alpha},$

with $\chi^{2}_{p}$ being a chi-square random variable with degree of freedom
$p$ and $\xi_{\alpha,p}$ being the $\alpha$-quantile of $\chi^{2}_{p}$.

###  2.2 Repeated sampling inference under re-randomization

Hereinafter, we use notation
“${\boldsymbol{X}}_{n}\overset{\cdot}{\sim}{\boldsymbol{Y}}_{n}$” for the fact
that two sequences of random vectors ${\boldsymbol{X}}_{n}$ and
${\boldsymbol{Y}}_{n}$ converge weakly to the same distribution with sample
size $n$ going to infinite, and notation
“${\boldsymbol{X}}\sim{\boldsymbol{Y}}$” for the fact that two random vectors
${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$ follow the same distribution.

Under complete randomization, Theorem 3 in Li and Ding, (2017) implies that
$\sqrt{N}\big{(}\hat{\tau}-\tau,\bm{D}^{T}\big{)}$ has the following
covariance matrix:

$\bm{V}_{N}=\left(\begin{array}[]{cc}V_{\tau\tau,N}&\bm{V}_{\tau x,N}\\\
\bm{V}_{x\tau,N}&\bm{V}_{xx,N}\end{array}\right)=\left(\begin{array}[]{cc}r_{1}^{-1}S^{2}_{Y(1)}+r_{0}^{-1}S^{2}_{Y(0)}-S^{2}_{\tau}&r_{1}^{-1}\bm{S}^{2}_{Y(1),\bm{X}}+r_{0}^{-1}\bm{S}^{2}_{Y(0),\bm{X}}\\\
r_{1}^{-1}\bm{S}^{2}_{\bm{X},Y(1)}+r_{0}^{-1}\bm{S}^{2}_{\bm{X},Y(0)}&\bm{(}r_{0}r_{1})^{-1}\bm{S}^{2}_{\bm{X}}\end{array}\right),$
(1)

where $r_{0}=N_{c}/N$ and $r_{1}=N_{t}/N$, finite-population variance and
covariance matrices
$S^{2}_{Y(z)}=(N-1)^{-1}\sum_{i=1}^{N}\big{(}Y_{i}(z)-\bar{Y}(z)\big{)}^{2},\quad
S^{2}_{\tau}=(N-1)^{-1}\sum_{i=1}^{N}(\tau_{i}-\tau)^{2},\quad\bm{S}^{2}_{\bm{X}}=(N-1)^{-1}\sum_{i=1}^{N}\big{(}\bm{X}_{i}-\bar{\bm{X}}\big{)}\big{(}\bm{X}_{i}-\bar{\bm{X}}\big{)}^{T},\quad\bm{S}^{2}_{Y(z),\bm{X}}=\bm{S}^{2}_{\bm{X},Y(z)}=(N-1)^{-1}\sum_{i=1}^{N}\big{(}Y_{i}(z)-\bar{Y}(z)\big{)}\big{(}\bm{X}_{i}-\bar{\bm{X}}\big{)}^{T}$,
$z=0$ or $1$. Moreover, according to finite population central limit theorem,
i.e., Theorem 5 in Li and Ding, (2017), when Condition 1 below is satisfied,
$\sqrt{N}\big{(}\hat{\tau}-\tau,\bm{D}^{T}\big{)}$ will converge to a Gaussian
distribution weakly as the sample size $N$ approaches infinity, i.e.,

$\sqrt{N}\big{(}\hat{\tau}-\tau,\bm{D}^{T}\big{)}\overset{\cdot}{\sim}(A,{\boldsymbol{B}}^{T}),\text{
where }(A,{\boldsymbol{B}}^{T})\sim\mathcal{N}(\bm{0},\bm{V}_{N}).$ (2)

For simplicity, we omit $N$ in finite-population covariance matrices
$\bm{V}_{N},V_{\tau\tau,N},\bm{V}_{\tau x,N},\bm{V}_{x\tau,N}$ and
$\bm{V}_{xx,N}$, and denote their limits as
$\bm{V}_{\infty},V_{\tau\tau,\infty},\bm{V}_{\tau
x,\infty},\bm{V}_{x\tau,\infty}$ and $\bm{V}_{xx,\infty}$ hereinafter.

######  Condition 1.

As $N\rightarrow\infty$, for $z=0,1$, (1) $r_{z}$ has positive limits; (2) the
finite population variances and covariances $S^{2}_{Y(z)}$, $S^{2}_{\tau}$,
$\bm{S}^{2}_{\bm{X}}$ and $\bm{S}^{2}_{\bm{X},Y(z)}$ have finite limiting
values, and the limit of $\bm{S}^{2}_{\bm{X}}$ is nonsingular; and (3)
$\max_{1\leq i\leq N}|Y_{i}(z)-\bar{Y}(z)|^{2}/N\rightarrow 0$ and
$\max_{1\leq i\leq N}||\bm{X}_{i}-\bar{\bm{X}}||^{2}/N\rightarrow 0$.

Moreover, Li et al., (2018) has established more detailed asymptotic results
for a special family of re-randomization mechanisms whose acceptance rule
depends on $\sqrt{N}\bm{D}$ and $\bm{V}_{xx}$ only. Let $\Phi$ be the family
of all re-randomization mechanisms that is of the form
$\phi\big{(}\sqrt{N}\bm{D},\bm{V}_{xx}\big{)}$ and satisfies the Condition 2
below. Let $\mathcal{B}_{\phi}=\\{\bm{\mu}:\phi(\bm{\mu},\bm{V}_{xx})=1\\}$ be
the acceptance region for $\sqrt{N}\bm{D}$ induced by $\phi$, and
$\mathcal{B}_{\phi,\infty}=\\{\bm{\mu}:\phi(\bm{\mu},\bm{V}_{xx,\infty})=1\\}$
be the limit of $\mathcal{B}_{\phi}$. Li et al., (2018) proved that

$\left.\left(\begin{array}[]{c}\sqrt{N}(\hat{\tau}-\tau)\\\
\sqrt{N}\bm{D}\end{array}\right)\right|\sqrt{N}\bm{D}\in\mathcal{B}_{\phi}\overset{\cdot}{\sim}\left.\left(\begin{array}[]{c}A\\\
\bm{B}\end{array}\right)\right|\bm{B}\in\mathcal{B}_{\phi},\
\forall\phi\in\Phi,$ (3)

where $(A,\bm{B}^{T})\sim\mathcal{N}(\bm{0},\bm{V})$ with its limit
$(A_{\infty},\bm{B}_{\infty}^{T})\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{V}}_{\infty})$.

######  Condition 2.

The covariate balance criterion $\phi(\cdot,\cdot)$ satisfies: (1)
$\phi(\cdot,\cdot)$ is almost surely continuous; (2) for $\bm{B}\sim
N(0,\bm{U})$, ${\mathbb{P}}(\phi(\bm{B},\bm{U})=1)>0$ for any $\bm{U}>0$, and
${\text{var}}(\bm{B}|\phi(\bm{B},\bm{U})=1)$ is a continuous function of
$\bm{U}$; (3) $\phi(\bm{\mu},\bm{U})=\phi(-\bm{\mu},\bm{U})$ for all
$\bm{\mu}$ and $\bm{U}>0$.

The symmetry condition for $\phi$ in Condition 2 guarantees that $\hat{\tau}$
is an unbiased estimator of $\tau$ under any re-randomization mechanism
$\phi\in\Phi$. Moreover, based on (3), we have

$\begin{split}\big{.}\sqrt{N}(\hat{\tau}-\tau)\big{|}\sqrt{N}\bm{D}\in\mathcal{B}_{\phi}&\overset{\cdot}{\sim}\left.(A\right|\bm{B}\in\mathcal{B}_{\phi})\sim\epsilon+(\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{B}|\bm{B}\in\mathcal{B}_{\phi}),\end{split}$

where $\epsilon\sim\mathcal{N}(0,V_{\tau\tau}-\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau})$ is independent of $\bm{B}$. Thus, the
asymptotic sampling variance of $\sqrt{N}(\hat{\tau}-\tau)$ under a re-
randomization mechanism $\phi\in\Phi$ is

${\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}=(1-R^{2}_{\infty})V_{\tau\tau,\infty}+{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right),$
(4)

where $R^{2}_{\infty}$ and ${\boldsymbol{\beta}}_{\infty}$ are the limits of
the constant scalar and vector defined below:

$\displaystyle R^{2}$ $\displaystyle\triangleq$ $\displaystyle\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau}/V_{\tau\tau},$ (5) $\displaystyle\bm{\beta}$
$\displaystyle\triangleq$ $\displaystyle\bm{V}_{xx}^{-1}\bm{V}_{x\tau}.$ (6)

Clearly, $\bm{\beta}$ and $R^{2}$ are the regression coefficient for
individual causal effect and covariates, and the squared multiple correlation,
respectively.

In practice, as suggested by Li et al., (2018), we can measure the efficiency
improvement by a re-randomization procedure $\phi$ relative to complete
randomization via the _Percent Reduction in Asymptotic Sampling Variance_
(PRIASV) of $\sqrt{N}(\hat{\tau}-\tau)$ defined below:

$\text{PRIASV}_{\phi}\triangleq
100\times\left[1-\frac{{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\mid\phi=1\big{)}}{{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\big{)}}\right]=100\times\left(1-r_{\phi}\right)\times
R^{2}_{\infty},$ (7)

where

$r_{\phi}=\frac{{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right)}{{\boldsymbol{\beta}}^{T}_{\infty}\bm{V}_{xx,\infty}{\boldsymbol{\beta}}_{\infty}}=\frac{{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right)}{{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\right)}.$
(8)

##  3 Bayesian perspective for re-randomization

In this study, we aim to propose novel re-randomization mechanisms that are
more efficient than the classic ReM-based ones. The key idea roots in the
observation that all ReM-based re-randomization mechanisms, including ReM,
ReMT, PCA-ReM and Ridge-ReM, fail to take full consideration of the
heterogeneity of different covariates with respect to the causal effect, while
there is a clear insight that covariates associated with causal effect more
closely should enjoy higher priority to be balanced. Although ReMT considers
covariate heterogeneity by partitioning the covariates into tiers and giving
higher priority to balancing covariates in the leading tiers, such a strategy
is sub-optimal because it treats covariates in the same tier equally. We have
better choices if we can quantify the heterogeneity of covariates with respect
to causal effect more precisely.

###  3.1 Re-randomization under the oracle

To overcome the limitation of the ReM-based re-randomization mechanisms, we
seek a principled way to design more efficient re-randomization mechanisms.
For any $\phi\in\Phi$, we can define its asymptotic acceptance rate in the
sequence of finite populations as

$\gamma_{\phi}=\lim_{N\rightarrow\infty}{\mathbb{P}}\big{(}\sqrt{N}\bm{D}\in\mathcal{B}_{\phi}\big{)}={\mathbb{P}}\big{(}\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\big{)}.$

For any $\alpha\in(0,1)$, let

$\Phi_{\alpha}=\\{\phi\in\Phi:\gamma_{\phi}=\alpha\\}$

be the subset of $\Phi$ that covers all re-randomization mechanisms with
asymptotic acceptance rate $\alpha$. In principle, it is ideal to find the
most efficient re-randomization mechanism in $\Phi_{\alpha}$ that minimizes
the asymptotic conditional sampling variance
${\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}$, i.e,

$\phi^{*}=\arg\min_{\phi\in\Phi_{\alpha}}{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}.$

The following theorem gives the solution to this optimization problem.

######  Theorem 1.

Suppose Condition 1 holds and $\bm{\beta}$ is known. The following
$\bm{\beta}$-specific re-randomization procedure

$\phi_{\bm{\beta}}\big{(}\sqrt{N}\bm{D},\bm{V}_{xx}\big{)}=I\big{(}d_{\bm{\beta}}\leq\sigma^{2}_{\bm{\beta}}\cdot\xi_{\alpha,1}\big{)}$

satisfies $\phi_{\bm{\beta}}\in\Phi_{\alpha}$ and

$\begin{split}{\text{var}}_{a}\big{(}\big{.}\sqrt{N}(\hat{\tau}-\tau)\big{|}\phi_{\bm{\beta}}=1\big{)}\leq{\text{var}}_{a}\big{(}\big{.}\sqrt{N}(\hat{\tau}-\tau)\big{|}\phi=1\big{)},\
\forall\phi\in\Phi_{\alpha},\end{split}$

where $\xi_{\alpha,1}$ is the $\alpha$-quantile of the $\chi^{2}_{1}$
distribution,

$\displaystyle d_{{\boldsymbol{\beta}}}$ $\displaystyle=$
$\displaystyle\big{(}\sqrt{N}\bm{D}\big{)}^{T}\bm{\beta}\bm{\beta}^{T}\big{(}\sqrt{N}{\boldsymbol{D}}\big{)}=N\bm{D}^{T}\bm{\beta}\bm{\beta}^{T}\bm{D},$
$\displaystyle\text{and}\ \ \sigma^{2}_{\bm{\beta}}$ $\displaystyle=$
$\displaystyle\bm{\beta}^{T}\bm{V}_{xx}\bm{\beta}={\boldsymbol{V}}_{\tau
x}{\boldsymbol{V}}_{xx}^{-1}{\boldsymbol{V}}_{x\tau}.$

Theorem 1 shows that $\phi_{\bm{\beta}}$ is the optimal re-randomization
mechanism in $\Phi_{\alpha}$ when $\bm{\beta}$ is known, and
$\phi_{\bm{\beta}}$ leads to a ribbon-shaped acceptance region in the space of
$\bm{D}$, i.e.,

$\mathcal{B}_{\phi_{\bm{\beta}}}=\left\\{\sqrt{N}\bm{D}:\big{(}\sqrt{N}\bm{\beta}^{T}\bm{D}\big{)}^{2}\leq\bm{\beta}^{T}\bm{V}_{xx}\bm{\beta}\cdot\xi_{\alpha,1}\right\\},$

whose orientation is orthogonal to the direction of $\bm{\beta}$. In addition,
we note that

$\begin{split}\bm{\beta}^{T}\bm{D}=\frac{1}{N_{t}}\sum_{i=1}^{N}\bm{Z}_{i}(\bm{\beta})W_{i}-\frac{1}{N_{c}}\sum_{i=1}^{N}\bm{Z}_{i}(\bm{\beta})(1-W_{i})=\bar{Z}_{t}(\bm{\beta})-\bar{Z}_{c}(\bm{\beta}),\end{split}$

where $Z_{i}(\bm{\beta})=\bm{X}_{i}^{T}\bm{\beta}$ is the index variable
obtained by taking the inner product of $\bm{X}_{i}$ and $\bm{\beta}$.
Therefore, controlling
$d_{\bm{\beta}}=N\bm{D}^{T}\bm{\beta}\bm{\beta}^{T}\bm{D}$ as described in
$\phi_{\bm{\beta}}$ is equivalent to controlling the imbalance of the index
variable $Z_{i}(\bm{\beta})=\bm{X}_{i}^{T}\bm{\beta}$ in the Euclidean space.
Thus, when $\bm{\beta}$ is known, we can do dimension reduction first by
projecting the original high-dimensional covariates into an one-dimensional
sub-space spanned by $\bm{\beta}$, and then achieve more efficient covariate
balance after the dimension reduction. Hereinafter, we refer to
$\phi_{\bm{\beta}}$ as the _re-randomization under the oracle_ (ReO), as it is
the most efficient re-randomization mechanism achievable in the asymptotic
perspective when $\bm{\beta}$ is precisely known.

The theorem below gives the asymptotic distribution and PRIASV of
$\sqrt{N}(\hat{\tau}-\tau)$ under $\phi_{\bm{\beta}}$ relative to complete
randomization.

######  Theorem 2.

Under the oracle $\phi_{\bm{\beta}}$ and Condition 1,

$\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{{\boldsymbol{\beta}}}=1\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}}\big{(}\sqrt{1-R^{2}}\cdot\varepsilon_{0}+\sqrt{R^{2}}\cdot
L_{1,\xi_{\alpha,1}}\big{)},$ (9)

where $\varepsilon_{0}\sim\mathcal{N}(0,1)$ is independent of
$L_{1,\xi_{\alpha,1}}\sim(D|D^{2}\leq\xi_{\alpha,1})$ with
$D\sim\mathcal{N}(0,1)$. Moreover,

$\text{PRIASV}_{ReO}=100\times(1-v_{\alpha,1})R^{2}_{\infty},$ (10)

where
$v_{\alpha,1}={\mathbb{E}}\left(\chi^{2}_{1}|\chi^{2}_{1}\leq\xi_{\alpha,1}\right)$.
Compared with ReM, we have

$\frac{\text{PRIASV}_{ReM}}{\text{PRIASV}_{ReO}}=\frac{1-v_{\alpha,p}}{1-v_{\alpha,1}}=O\left(\frac{1}{\sqrt{p}}\right).$
(11)

Figure 1: Ratio of $\text{PRIASV}_{ReO}$ and $\text{PRIASV}_{ReM}$ for varying
covariate dimension $p$ and logarithm acceptance probability $\alpha$.

Figure 1 illustrates the superiority of ReO over ReM in terms of PRIASV by a
heat map of $\text{PRIASV}_{ReO}/\text{PRIASV}_{ReM}$ as a function of
covariate dimension $p$ and acceptance probability $\alpha$ based on Eq. (11).
From the figure, we observe an increasing performance advantage for ReO over
ReM with the increase of covariate dimension $p$.

###  3.2 Bayesian Criterion for Re-randomization

A critical limitation of ReO, however, is that we must know
${\boldsymbol{\beta}}$ for its implementation, which is often not realistic.
In this subsection, we retreat to a more realistic case where prior
information about ${\boldsymbol{\beta}}$, instead of ${\boldsymbol{\beta}}$
itself, can be obtained, and try to establish a Bayesian criterion for re-
randomization.

Suppose ${\boldsymbol{\beta}}_{\infty}$, which is the limit of
${\boldsymbol{\beta}}$, is associated with a prior distribution $\pi$. Under
the Bayesian framework, we quantify the performance of a re-randomization
mechanism $\phi$ under prior $\pi$ by averaging its
${\boldsymbol{\beta}}_{\infty}$-specific asymptotic sampling variance with
respect to $\pi$. Formally, define the prior-integrated asymptotic sampling
variance

$\begin{split}{\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\mid\phi=1\big{)}&={\mathbb{E}}_{\pi}\left[{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}\right]\\\
&\vspace{1cm}={\mathbb{E}}_{\pi}\Big{[}(1-R^{2}_{\infty})V_{\tau\tau,\infty}+{\mathbb{E}}\left(\bm{B}_{\infty}^{T}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right)\Big{]}\end{split}$

as the objective function for evaluating a re-randomization mechanism
$\phi\in\Phi$ under the Bayesian setting. We aim to search $\Phi_{\alpha}$ for
the solution of the optimization problem below:

$\phi^{*}=\arg\min_{\phi\in\Phi_{\alpha}}{\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}.$

The following theorem gives the solution to this optimization problem.

######  Theorem 3.

Suppose Condition 1 is satisfied and ${\boldsymbol{\beta}}_{\infty}$ follows
prior distribution $\pi$. Let $\bm{\mu}_{\pi}$ and
${\boldsymbol{\Sigma}}_{\pi}$ be the mean vector and covariance matrix of
$\pi$ respectively, and
${\boldsymbol{\Lambda}}_{\pi}=\bm{\mu}_{\pi}\bm{\mu}_{\pi}^{T}+{\boldsymbol{\Sigma}}_{\pi}$
be the corresponding characteristic matrix. Define the following prior-induced
distance

$d_{\pi}=\big{(}\sqrt{N}{\boldsymbol{D}}\big{)}^{T}\bm{\Lambda}_{\pi}\big{(}\sqrt{N}{\boldsymbol{D}}\big{)}$
(12)

as the measure of covariate imbalance. The Bayesian re-randomization mechanism

$\phi_{\pi}\big{(}\sqrt{N}\bm{D},\bm{V}_{xx}\big{)}=I\big{(}d_{\pi}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}$
(13)

satisfies $\phi_{\pi}\in\Phi_{\alpha}$, and

$\begin{split}{\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1\big{)}\leq{\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\mid\phi=1\big{)},\
\forall\phi\in\Phi_{\alpha},\end{split}$

where $\xi_{\alpha,{\boldsymbol{\lambda}}}$ is $\alpha$-quantile of
$\chi^{2}_{{\boldsymbol{\lambda}}}=\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2}$ with
${\boldsymbol{\lambda}}=(\lambda_{1},\ldots,\lambda_{p})$ being vector of
eigenvalues of matrix
$\bm{P}=\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}\bm{\Lambda}_{\pi}\bm{V}_{xx}^{1/2}$,
$Z_{1},...,Z_{p}$ being i.i.d. random variables from standard Gaussian
distribution, and $\bm{V}_{xx}^{1/2}$ being the Cholesky square root of
$\bm{V}_{xx}$ satisfying
$\bm{V}_{xx}^{1/2}\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}=\bm{V}_{xx}$.

The theorem tells us that $\phi_{\pi}$, the optimal re-randomization procedure
in $\Phi_{\alpha}$ under the Bayesian criterion, achieves the smallest prior-
integrated asymptotic sampling variance and leads to an elliptical acceptance
region, whose shape is determined by $\bm{\Lambda}_{\pi}$, i.e., the
characteristic matrix of the prior distribution $\pi$ (see Figure 2 for the
graphical illustration). Analogous to
$d_{\bm{\beta}}=N\bm{D}^{T}\bm{\beta}\bm{\beta}^{T}\bm{D}$ in ReO,
$d_{\pi}=N{\boldsymbol{D}}^{T}\bm{\Lambda}_{\pi}{\boldsymbol{D}}$ plays a
similar role in achieving covariance balance by projecting the high
dimensional $\bm{D}$ into a lower dimensional space. Hereinafter, we refer to
$\phi_{\pi}$ as the _re-randomization under the Bayesian criterion_ (ReB).

Figure 2: Acceptance region of ${\boldsymbol{D}}=(D_{1},D_{2})$ under
different criteria for re-randomization with only 2 covariates.
${\boldsymbol{\beta}}=c(1,0)^{T}$ in ReO and
$\bm{\Lambda}_{\pi}={\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}+\mbox{diag}\\{0.25,0.25\\}$
in ReB when acceptance probability $\alpha=0.5$.

Theorem 3 shows that the Bayesian criterion for re-randomization leads to the
acceptance rule
$\phi_{\pi}=I(N{\boldsymbol{D}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{D}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}})$,
which is completely determined by
${\boldsymbol{\Lambda}}_{\pi}={\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}+{\boldsymbol{\Sigma}}_{\pi}$,
the characteristic matrix of prior distribution $\pi$. By specifying
${\boldsymbol{\mu}}_{\pi}$ and ${\boldsymbol{\Sigma}}_{\pi}$ to different
values, we find that Bayesian re-randomization takes many other re-
randomization criteria as special cases. According to Table 1, ReB degenerates
to ReO, ReM, Ridge-ReM and PCA-ReM under different specifications of
${\boldsymbol{\mu}}_{\pi}$ and ${\boldsymbol{\Sigma}}_{\pi}$. Thus, ReB
provides us a unified framework to understand various re-randomization
criteria in the literature, and highlight their differences from the Bayesian
perspectives.

Table 1: ReB leads to other re-randomization criteria under specific prior specifications ${\boldsymbol{\mu}}_{\pi}$ | ${\boldsymbol{\Sigma}}_{\pi}$ | ${\boldsymbol{\Lambda}}_{\pi}$ | $d_{\pi}$ | Criteria
---|---|---|---|---
${\boldsymbol{\beta}}$ | $\bm{0}_{p\times p}$ | ${\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}$ | $N{\boldsymbol{D}}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}{\boldsymbol{D}}$ | ReO
$\bm{0}_{p}$ | ${\boldsymbol{V}}_{xx}^{-1}=(N{\boldsymbol{\Sigma}}_{\boldsymbol{D}})^{-1}$ | $(N{\boldsymbol{\Sigma}}_{\boldsymbol{D}})^{-1}$ | ${\boldsymbol{D}}^{T}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{-1}{\boldsymbol{D}}$ | ReM
$\bm{0}_{p}$ | $N^{-1}\big{(}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}+\lambda{\boldsymbol{I}}_{p}\big{)}^{-1}$ | $N^{-1}\big{(}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}+\lambda{\boldsymbol{I}}_{p}\big{)}^{-1}$ | ${\boldsymbol{D}}^{T}\big{(}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}+\lambda{\boldsymbol{I}}_{p}\big{)}^{-1}{\boldsymbol{D}}$ | Ridge-ReM
$\bm{0}_{p}$ | $\big{(}N{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{(k)}\big{)}^{-1}$ | $\big{(}N{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{(k)}\big{)}^{-1}$ | ${\boldsymbol{D}}^{T}\big{(}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{(k)}\big{)}^{-1}{\boldsymbol{D}}$ | PCA-ReM

An interesting property of ReB is that its acceptance region shows invariance
to any scale transformation of prior distribution $\pi$. Intuitively, ReB
enjoys such a property because only the direction, not the magnitude, of
${\boldsymbol{\beta}}$ matters in determining the relative importance of
covariates, and thus the acceptance region. The following corollary based on
Thm. 3 summarizes this phenomenon explicitly. The detailed proof can be found
in the Supplementary Material.

######  Corollary 1.

Suppose that $\phi_{\pi_{1}}$ and $\phi_{\pi_{2}}$ are two ReB procedures in
$\Phi_{\alpha}$ with different prior specifications $\pi_{1}$ and $\pi_{2}$,
respectively. We have

$\begin{split}\phi_{\pi_{2}}=\phi_{\pi_{1}}\ \text{as long as
}\bm{\Lambda}_{\pi_{2}}=r\bm{\Lambda}_{\pi_{1}}\ \text{for some
}r>0.\end{split}$

###  3.3 Frequentist properties of ReB

Although ReB is established from the Bayesian perspective, we can study its
performance with comparison to other re-randomization mechanisms, e.g., ReM,
from the frequentist perspective. According to (4), the relative efficiency of
$\phi_{\pi}$ with respect to complete randomization can be quantified in term
of PRIASV of $\sqrt{N}(\hat{\tau}-\tau)$:

$\text{PRIASV}_{ReB}=100\times(1-r_{\phi_{\pi}})R^{2}_{\infty},$ (14)

where

$r_{\phi_{\pi}}=v_{\alpha,\pi}\triangleq\frac{{\mathbb{E}}\big{(}\bm{B}_{\infty}^{T}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\mid\bm{B}_{\infty}\in\mathcal{B}_{\phi_{\pi},\infty}\big{)}}{{\mathbb{E}}\big{(}\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\big{)}}.$
(15)

When ${\boldsymbol{\Lambda}}_{\pi}$ is of full rank, the following theorem
gives the asymptotic distribution as well as the explicit form of
$v_{\alpha,\pi}$ for $\sqrt{N}(\hat{\tau}-\tau)$ under ReB.

######  Theorem 4.

Under the Bayesian setting with $\pi$ as the prior distribution of
${\boldsymbol{\beta}}_{\infty}$ and Condition 1, when
${\boldsymbol{\Lambda}}_{\pi}$ is symmetric and positive definite with
Cholesky decomposition
${\boldsymbol{\Lambda}}_{\pi}={\boldsymbol{M}}{\boldsymbol{M}}^{T}$, then
$\phi_{\pi}$ leads to

$\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}(1-R^{2})}\cdot\varepsilon_{0}+{\boldsymbol{\beta}}^{T}({\boldsymbol{M}}^{-1})^{T}\bm{Q}\mid\bm{Q}^{T}\bm{Q}\leq\xi_{\alpha,{\boldsymbol{\lambda}}},$
(16)

where
$\bm{Q}\sim\mathcal{N}(\bm{0},{\boldsymbol{M}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{M}})$,
and
${\boldsymbol{\beta}}^{T}({\boldsymbol{M}}^{-1})^{T}\bm{Q}\sim\mathcal{N}(0,{\boldsymbol{V}}_{\tau
x}{\boldsymbol{V}}_{xx}{\boldsymbol{V}}_{x\tau})$. Moreover,

$v_{\alpha,\pi}=\frac{{\boldsymbol{\beta}}^{T}_{\infty}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}_{\infty}\bm{C}_{\infty}\bm{\Gamma}_{\infty}\bm{M}^{-1}{\boldsymbol{\beta}}_{\infty}}{{\boldsymbol{V}}_{\tau
x,\infty}{\boldsymbol{V}}_{xx,\infty}{\boldsymbol{V}}_{x\tau,\infty}}=\frac{{\boldsymbol{\beta}}^{T}_{\infty}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}_{\infty}\bm{C}_{\infty}\bm{\Gamma}_{\infty}\bm{M}^{-1}{\boldsymbol{\beta}}_{\infty}}{{\boldsymbol{\beta}}^{T}_{\infty}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}_{\infty}\bm{\Lambda}_{\infty}\bm{\Gamma}_{\infty}\bm{M}^{-1}{\boldsymbol{\beta}}_{\infty}},$
(17)

where $\bm{\Lambda}=diag\\{\lambda_{1},\ldots,\lambda_{p}\\}$,
$\bm{C}=diag\\{c_{1},\ldots,c_{p}\\}$ with
$c_{i}={\mathbb{E}}\left(\lambda_{i}Z_{i}^{2}\mid\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)$
for $i=1,\ldots,p$, ${\boldsymbol{\lambda}}=(\lambda_{1},\ldots,\lambda_{p})$,
and $Z_{1},...,Z_{p}$ being i.i.d. standard normal distributed random
variables. $\bm{\Gamma}$ is an orthogonal matrix such that
$\bm{M}^{T}\bm{V}_{xx}\bm{M}=\bm{\Gamma}^{T}\bm{\Lambda}\bm{\Gamma}$.
${\boldsymbol{\Gamma}}_{\infty},\bm{C}_{\infty}$ and
${\boldsymbol{\Lambda}}_{\infty}$ are limits of ${\boldsymbol{\Gamma}},\bm{C}$
and ${\boldsymbol{\Lambda}}$ in the sequence of finite populations.

Theorem 4 gives the explicit expression of $v_{\alpha,\pi}$ when
${{\boldsymbol{\Lambda}}}_{\pi}$ is of full rank, which is easy to hold in
practice. The value of $v_{\alpha,\pi}$ depends on the eigenvalues
${\boldsymbol{\lambda}}$ of
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{{\boldsymbol{\Lambda}}}_{\pi}{\boldsymbol{V}}_{xx}^{1/2}$
via
${\mathbb{E}}\big{(}Z_{i}^{2}\mid\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)},i=1,\ldots,p.$
This expectation will degenerate to
$v_{\alpha,p}={\mathbb{E}}\big{(}Z_{i}^{2}\mid\sum_{j=1}^{p}Z_{j}^{2}\leq\xi_{\alpha,p}\big{)}$
in ReM, where ${{\boldsymbol{\Lambda}}}_{\pi}=\bm{V}_{xx}^{-1}$ and thus
$\lambda_{1}=\cdots=\lambda_{p}=1$.

To build a geometry intuition on the effect of the precision of the prior on
the performance of ReB, we further consider the case where the prior
distribution $\pi$ degenerates to a point mass. The theorem below shows a
simpler form of the asymptotic distribution and PRIASV of
$\sqrt{N}(\hat{\tau}-\tau)$ in this scenario.

######  Theorem 5.

Under the Bayesian setting where the prior distribution $\pi$ degenerates to a
point mass at ${\boldsymbol{\mu}}_{\pi}$, which means that
$\bm{\Sigma}_{\pi}=\bm{0}$ and
${{\boldsymbol{\Lambda}}}_{\pi}={\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}$,
we have

$\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}}\left(\sqrt{1-R^{2}_{\pi}}\cdot\varepsilon_{0}+\sqrt{R^{2}_{\pi}}\cdot
L_{1,\xi_{\alpha,1}}\right),$ (18)

where

$R^{2}_{\pi}=\frac{({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})^{2}}{({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})}\cdot
R^{2}=(\cos\theta)^{2}\cdot R^{2}$ (19)

with $\theta$ being the angle between
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\mu}}_{\pi}$ and
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\beta}}$. Denote
the limit of $\theta$ as $\theta_{\infty}$, then

$v_{\alpha,\pi}=1-(1-v_{\alpha,1})(\cos\theta_{\infty})^{2}.$ (20)

Intuitively, when
${{\boldsymbol{\Lambda}}}_{\pi}={\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}$,
${\boldsymbol{\mu}}_{\pi}$ becomes our prior guess for ${\boldsymbol{\beta}}$.
Thus, Thm. 5 tells us that the asymptotic performance of ReB only depends on
the limiting angle between ${\boldsymbol{\mu}}_{\pi}$ and
${\boldsymbol{\beta}}$ after an affine transformation. Apparently, a smaller
$\theta_{\infty}$ would lead to a re-randomization procedure that behaves more
like ReO, and thus enjoys a better chance to outperform ReM. The corollary
below gives the sufficient and necessary condition for ReB (under the
degenerated prior) to outperform ReM in terms of the admissible value range of
$\theta_{\infty}$.

######  Corollary 2.

Under the same setting as in Thm. 5 with $\theta_{\infty}$ standing for the
limit of the angle between
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\mu}}_{\pi}$ and
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\beta}}$. ReB
$\phi_{\pi}$ outperforms ReM $\phi_{M}$ if and only if

$(\cos\theta_{\infty})^{2}>\frac{1-v_{\alpha,p}}{1-v_{\alpha,1}}.$

Corollary 2 implies that once the probability mass of the prior distribution
$\pi$ concentrates in a cone centered in the specific direction along
${\boldsymbol{\beta}}$ (after the transformation according to
${\boldsymbol{V}}_{xx}^{1/2}$), ReB would benefit from the prior information
and outperform ReM. Figure 3 illustrates how the ratio between the opening
angle of the cone (when restricted on $[0,\pi/2]$) and $\pi/2$ changes with
$\alpha$ and $p$ when ${\boldsymbol{V}}_{xx}$ is an identity matrix. From the
figure, we can see that as the dimensionality $p$ gets larger, the admissible
limiting value range of $\theta$ enlarges to a wider region spanning from 0 to
nearly $\pi/2$, indicating that the superiority of ReB over ReM in practice
becomes less sensitive to the mis-specification of ${\boldsymbol{\beta}}$ as
$p$ gets larger. Such a result gives us the confidence that ReB would be a
useful re-randomization procedure in practice as long as our prior knowledge
about ${\boldsymbol{\beta}}$ is not that misleading.

Figure 3: Largest limiting angle between prior vector
${\boldsymbol{\mu}}_{\pi}$ and ${\boldsymbol{\beta}}$ after transformation
according to ${\boldsymbol{V}}_{xx}^{1/2}$ that guarantees out-performance of
ReB over ReM for varying covariate dimension $p$ and logarithm acceptance
probability $\alpha$.

##  4 Implementing ReB via a two-stage experiment

The theoretical analysis in the previous section shows that ReB is a more
efficient re-randomization procedure than ReM when the prior distribution
$\pi$ is informative to highlight ${\boldsymbol{\beta}}$. In case that such an
informative prior is not available in advance, we still can implement ReB if a
two-stage experiment is feasible, where the first stage involves a small pilot
randomized experiment with $N_{1}$ randomly selected units to learn an
informative prior for the unknown ${\boldsymbol{\beta}}$, while the second
stage implements a larger-scale re-randomized experiment with other $N_{2}$
randomly selected units via ReB guided by the learned prior in pursuit of a
more efficient causal estimator.

In the pilot experiment in the first stage, we implement a balanced completely
randomized experiment to get a preliminary causal effect estimator
$\hat{\tau}_{1}=\bar{Y}_{T}^{(1)}-\bar{Y}_{C}^{(1)}$ satisfying

${\mathbb{E}}\left(\hat{\tau}_{1}\right)=\tau_{1},\
{\text{var}}_{a}\big{(}\sqrt{N_{1}}(\hat{\tau}_{1}-\tau_{1})\big{)}={\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\big{)},$

where $\bar{Y}_{T}^{(1)}$ and $\bar{Y}_{C}^{(1)}$ are average responses in the
two arms of the first-stage experiment, $\tau_{1}$ and $\tau$ are average
treatment effects in the first stage and the complete study, respectively, and
$\hat{\tau}$ is the difference-in-means estimator of $\tau$ when we conduct a
complete randomization on the total $N=N_{1}+N_{2}$ units.

According to Eq’s. (1) and (6),

${\boldsymbol{\beta}}=r_{0}\cdot{\boldsymbol{S}}_{{\boldsymbol{X}}}^{-2}{\boldsymbol{S}}^{2}_{{\boldsymbol{X}},Y(1)}+r_{1}\cdot{\boldsymbol{S}}_{{\boldsymbol{X}}}^{-2}{\boldsymbol{S}}^{2}_{{\boldsymbol{X}},Y(0)}=r_{0}\cdot{\boldsymbol{\beta}}_{1}+r_{1}\cdot{\boldsymbol{\beta}}_{0},$

where ${\boldsymbol{\beta}}_{1}$ and ${\boldsymbol{\beta}}_{0}$ are projection
coefficient vector of ${\boldsymbol{Y}}(1)$ and ${\boldsymbol{Y}}(0)$ on
${\boldsymbol{X}}$, respectively. We can estimate ${\boldsymbol{\beta}}_{1}$
and ${\boldsymbol{\beta}}_{0}$ with observed outcomes and covariates in the
first stage experiment, denoted as $\hat{{\boldsymbol{\beta}}}_{1}$ and
$\hat{{\boldsymbol{\beta}}}_{0}$. Assisted with linear model, we obtain a
prior variance of $\hat{{\boldsymbol{\beta}}}_{1}$ and
$\hat{{\boldsymbol{\beta}}}_{0}$, denoted as ${\boldsymbol{V}}_{1}$ and
${\boldsymbol{V}}_{0}$. Then an informative prior distribution of
${\boldsymbol{\beta}}$ is

$\pi_{N_{1}}=\mathcal{N}\big{(}r_{0}\cdot\hat{{\boldsymbol{\beta}}}_{1}+r_{1}\cdot\hat{{\boldsymbol{\beta}}}_{0},r_{0}^{2}\cdot{\boldsymbol{V}}_{1}+r_{1}^{2}\cdot{\boldsymbol{V}}_{0}\big{)}.$
(21)

In the second stage, we implement ReB on the other $N_{2}$ units with the
$\pi_{N_{1}}$ derived in the first stage as the prior distribution of
${\boldsymbol{\beta}}$, leading to a Bayesian re-randomization procedure
referred to as $\phi_{\pi_{N_{1}}}$ with the balance criterion below:

$\phi_{\pi_{N_{1}}}\big{(}\sqrt{N_{2}}{\boldsymbol{D}}^{(2)},{\boldsymbol{V}}_{xx}^{(2)}\big{)}=I\big{(}N_{2}\cdot({\boldsymbol{D}}^{(2)})^{T}{\boldsymbol{\Lambda}}_{\pi_{N_{1}}}{\boldsymbol{D}}^{(2)}\leq\xi_{\alpha,{\boldsymbol{\lambda}}_{\pi_{N_{1}}}}\big{)},$

where $\bm{D}^{(2)}$ is the difference-in-means of covariates in the second-
stage experiment, $\bm{V}_{xx}^{(2)}$ is the sampling covariance of
$\sqrt{N_{2}}\cdot\bm{X}^{(2)}$, $\bm{X}^{(2)}$ is the covariate matrix in the
second stage and ${\boldsymbol{\lambda}}_{\pi_{N_{1}}}$ is the vector of
eigenvalues of matrix
$\bm{P}_{N_{1},{N_{2}}}=\Big{(}\big{(}{\boldsymbol{V}}_{xx}^{(2)}\big{)}^{1/2}\Big{)}^{T}{\boldsymbol{\Lambda}}_{\pi_{N_{1}}}\big{(}{\boldsymbol{V}}_{xx}^{(2)}\big{)}^{1/2}$.
This procedure leads to another causal effect estimator
$\hat{\tau}_{2}=\bar{Y}_{T}^{(2)}-\bar{Y}_{C}^{(2)}$ satisfying

${\mathbb{E}}\left(\hat{\tau}_{2}\right)=\tau_{2},\
{\text{var}}_{a}\big{(}\sqrt{N_{2}}(\hat{\tau}_{2}-\tau_{2})\big{)}=\left[1-(1-v_{\alpha,\pi_{N_{1}}})R^{2}_{\infty}\right]{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\big{)},$

where $\bar{Y}_{T}^{(2)}$ and $\bar{Y}_{C}^{(2)}$ are average responses in the
two arms of the second-stage experiment, and $\tau_{2}$ is the average
treatment effects in the second stage.

Denote the proportion of sample size in first stage as $\varrho$. Since
$\tau=\varrho\cdot\tau_{1}+(1-\varrho)\cdot\tau_{2}$,

$\hat{\tau}=\varrho\cdot\hat{\tau}_{1}+(1-\varrho)\cdot\hat{\tau}_{2}$

is an unbiased estimator of $\tau$. Denote the above two-stage re-
randomization procedure as $\phi_{BCRD\text{-}ReB}$. Similarly, BCRD in the
first stage can be replaced by ReM, resulting in an alternative two-stage
procedure referred to as $\phi_{ReM\text{-}ReB}$. The following theoren gives
the PRIASV of $\hat{\tau}$ under BCRD-ReB and ReM-ReB.

######  Theorem 6.

Suppose Condition 1 is satisfied, acceptance rate in the second stage of BCRD-
ReB is $\alpha_{1}$ and the acceptance rates in the two stages of ReM-ReB are
$\alpha_{2}$ and $\alpha_{3}$, respectively. When $N_{1}$ and $N_{2}$ both go
to infinity, $\varrho=N_{1}/N$ converges to a constant $\varrho^{*}$, the
PRIASV of $\hat{\tau}$ under BCRD-ReB is

$\text{PRIASV}_{BCRD\text{-}ReB}=100\times(1-\varrho^{*})(1-v_{\alpha_{1},1})R^{2}_{\infty}.$
(22)

Similarly, the PRIASV of $\tilde{\tau}_{N_{1},N_{2}}$ under ReM-ReB is

$\text{PRIASV}_{ReM\text{-}ReB}=100\times\left[\varrho^{*}(1-v_{\alpha_{2},p})+(1-\varrho^{*})(1-v_{\alpha_{3},1})\right]R^{2}_{\infty}.$
(23)

It is easy to see that both two-stage ReB procedures are more efficient than
complete randomization. ReM-ReB is more efficient than BCRD-ReB when the
acceptance rates in the second stage of both methods are the same, i.e.,
$\alpha_{1}=\alpha_{3}$. The following corollary is a simple application of
the above analysis and it shows that if the proportion of sample size in first
stage is small, BCRD-ReB and ReM-ReB with $\alpha_{1}=\alpha_{3}=\alpha$ are
asymptotically equivalent to ReO with acceptance rate $\alpha$.

######  Corollary 3.

Under the same settings in Theorem 6, if the acceptance rates in ReO, namely
$\alpha$, is the same as acceptance rate in the second stage of BCRD-ReB and
ReM-ReB, i.e., $\alpha=\alpha_{1}=\alpha_{3}$, then

$\text{PRIASV}_{BCRD\text{-}ReB}=\text{PRIASV}_{ReM\text{-}ReB}=\text{PRIASV}_{ReO},$

as $\varrho=N_{1}/N\rightarrow\varrho^{*}=0$.

Additionally, Eq. (21) allows defining a prior for ${\boldsymbol{\beta}}$ in
the two-stage procedure using a linear model. A simpler approach is to
degenerate $\pi_{N_{1}}$ into a point mass on
$r_{0}\cdot\hat{{\boldsymbol{\beta}}}_{1}+r_{1}\cdot\hat{{\boldsymbol{\beta}}}_{0}$,
using only the first moment. We refer to the procedure employing this
degenerated prior as BCRD-ReO and ReM-ReO, as their second-stage criterion
closely resembles that of ReO.

##  5 Simulation studies

In this section, we validate the theoretical analyses in previous sections and
compare various re-randomization procedures, including ReM, ReO, ReB and two-
stage ReB, via simulation. We generate simulated datasets with a linear model
for potential outcomes:

$\begin{split}Y_{i}(1)&=a+{\boldsymbol{\beta}}_{1}^{T}{\boldsymbol{X}}_{i}+\epsilon_{1,i},\
\ i=1,\ldots,N,\\\
Y_{i}(0)&={\boldsymbol{\beta}}_{0}^{T}{\boldsymbol{X}}_{i}+\epsilon_{0,i},\ \
i=1,\ldots,N.\end{split}$

Here, $a=5$, covariates
${\boldsymbol{X}}\sim\mathcal{N}(\bm{0}_{p},(1-\rho){\boldsymbol{I}}_{p}+\rho\bm{1}_{p}\bm{1}_{p}^{T})$,
$\epsilon_{1},\epsilon_{0}\sim\mathcal{N}(0,\sigma^{2}_{{\boldsymbol{\epsilon}}})$,
and regression coefficients
${\boldsymbol{\beta}}_{1}\sim\mathcal{N}(2\cdot\bm{1}_{p},2\sigma^{2}_{{\boldsymbol{\beta}}}\cdot\bm{I}_{p})$,
${\boldsymbol{\beta}}_{0}\sim\mathcal{N}(\bm{1}_{p},2\sigma^{2}_{{\boldsymbol{\beta}}}\cdot\bm{I}_{p})$.
We assume balanced treatment assignment with $r_{1}=r_{0}=0.5$, resulting in
the projection coefficient of treatment effect on difference-in-means of
covariates
${\boldsymbol{\beta}}=r_{1}\cdot{\boldsymbol{\beta}}_{0}+r_{0}\cdot{\boldsymbol{\beta}}_{1}\sim\mathcal{N}(1.5\cdot\bm{1}_{p},\sigma^{2}_{{\boldsymbol{\beta}}}\cdot\bm{I}_{p})$.
Denote the distribution of ${\boldsymbol{\beta}}$ as
$\pi({\boldsymbol{\beta}})$. We specify $p\in\\{2,5,10,20\\}$,
$\rho\in\\{0,0.2,0.5\\}$, $\sigma^{2}_{\boldsymbol{\beta}}\in\\{0,0.1,1\\}$.
For each combination of $(p,\rho,\sigma^{2}_{\boldsymbol{\beta}})$, we
porperly set $\sigma^{2}_{\boldsymbol{\epsilon}}$ to tune the average value of
$R^{2}$ in the datasets referred to as $\bar{R}^{2}\in\\{0.2,0.5,0.8\\}$, and
choose sample size $N\in\\{100,200,600,1200\\}$, leading to a configuration
space ${\mathcal{S}}$ containing $4\times 3\times 3\times 3\times 4=432$
combinations of the 5 controlling factors
$(p,\rho,\sigma^{2}_{\boldsymbol{\beta}},\bar{R}^{2},N)$. Table 2 summarizes
these factors, including a scheme factor for different re-randomization
procedures. Each configuration in ${\mathcal{S}}$ generates 100 independent
datasets, totaling 43,200 simulated datasets.

Table 2: Six major factors in the simulation study. Factor |  Levels |  Description
---|---|---
$p$ |  {2, 5, 10, 20} |  Dimensionality of ${\boldsymbol{X}}$
$\rho$ |  {0, 0.2, 0.5} |  Correlation coefficient of ${\boldsymbol{X}}$
$\sigma^{2}_{\boldsymbol{\beta}}$ |  $\\{0,0.1,1\\}$ |  Variance of $\pi({\boldsymbol{\beta}})$
$\bar{R}^{2}$ |  {0.2, 0.5, 0.8} |  Noise level of regression model
$N$ |  {100, 200, 600, 1200} |  Total sample size
Scheme |  {ReM, ReO, ReB} |  Re-randomization procedures

###  5.1 ReO and ReB versus ReM

First, we study the performance of ReO and ReB versus ReM in reducing the
estimation variance of $\hat{\tau}$ compared to Balanced Complete Randomized
Design (BCRD). These schemes are applied to our simulated datasets for
quantitative evaluation. We maintain a fixed acceptance probability of
$\alpha=0.05$, while assuming known covariate dimensionality $p$, sample size
$N$, and correlation $\rho$. In ReO, ${\boldsymbol{\beta}}$ is considered
known, while in ReB, the investigator has information about the sampling
distribution $\pi({\boldsymbol{\beta}})$ but not the specific
${\boldsymbol{\beta}}$ for each dataset, and ReM solely relies on simulated
data without any prior information about ${\boldsymbol{\beta}}$ or its
distribution.

For the $j$-th simulated dataset ${\mathcal{D}}_{ij}$ under the $i$-th
simulation configuration in ${\mathcal{S}}$, the estimation variance
${\text{var}}(\hat{\tau}\mid\phi=1)$ of a re-randomization procedure $\phi$ is
estimated based on 1000 accepted allocations from $\phi$ (referred to as
$\widehat{\text{var}}_{ij}(\hat{\tau}\mid\phi=1)$). Based on the estimated
$\widehat{\text{var}}_{ij}(\hat{\tau}\mid\phi=1)$, PRIV achieved by $\phi$
with respect to BCRD can be obtained for ${\mathcal{D}}_{ij}$ (denoted as
PRIV${}_{ij}(\phi)$), and the average performance of a re-randomization
procedure $\phi$ in simulation setting $i\in{\mathcal{S}}$ can be summarized
by

$\mbox{PRIV}_{i}(\phi)=\frac{1}{100}\sum_{j=1}^{100}\mbox{PRIV}_{ij}(\phi),$
(24)

which is equivalent to averaging the PRIV by prior of ${\boldsymbol{\beta}}$.

To assess the impact of the six factors in Table 2 on PRIV of $\hat{\tau}$, we
conducted an ANOVA for the estimated $\\{$PRIV${}_{ij}(\phi)\\}$ considering
these factors and their interactions. Table 3 presents the ANOVA outcomes with
the 5 most influential effects in explaining the variance of
PRIV${}_{ij}(\phi)$. The term “residual” encompasses all remaining effects
beyond these top 5 influential factors. Notably, the top 5 factors account for
nearly all PRIV${}_{ij}(\phi)$ variance, with the top 2 factors explaining
over 95% of the total variance. Specifically, $\bar{R}^{2}$ has the most
significant impact on PRIV variance, followed by scheme, covariate dimension
$p$, and their interactions, highlighting the re-randomization procedure’s
critical role. Moreover, sample size $N$ does not significantly affect PRIV.
These findings are consistent with the theoretical results presented in
Section 3.

Table 3: The 5 most influential effects for the ANOVA simulation results. Sources |  DF |  MS |  CP |  $F$-value
---|---|---|---|---
$\bar{R}^{2}$ |  2 |  3415.91 |  90.96 |  506659.51
Scheme |  2 |  206.21 |  96.45 |  30585.59
$p$ |  3 |  59.88 |  98.04 |  8881.31
$\mbox{Scheme}\times p$ |  6 |  32.32 |  98.91 |  4793.10
$\mbox{Scheme}\times\bar{R}^{2}$ |  4 |  21.20 |  99.47 |  3143.95
Residual |  129582 |  19.91 |  100 |
DF: “degrees of freedom”, MS: “mean square”, CP: “cumulative percentage of MS”
Table 4: Average PRIV achieved by different re-randomization procedures when $N=200$ and $\sigma^{2}_{\boldsymbol{\beta}}=1$. | $\bar{R}^{2}=0.2$ | $\bar{R}^{2}=0.5$ | $\bar{R}^{2}=0.8$
---|---|---|---
Scheme | $\rho=0$ | 0.2 | 0.5 | 0 | 0.2 | 0.5 | 0 | 0.2 | 0.5
| $p=2$
ReO | 19.36 | 19.61 | 20.15 | 49.50 | 49.71 | 49.75 | 79.74 | 79.75 | 79.78
ReB | 19.27 | 19.32 | 19.26 | 48.65 | 49.04 | 48.53 | 78.40 | 78.71 | 78.23
ReM | 19.03 | 18.98 | 19.05 | 48.07 | 48.24 | 48.03 | 77.57 | 77.60 | 77.56
| $p=5$
ReO | 20.12 | 20.45 | 19.90 | 49.71 | 49.69 | 49.91 | 79.69 | 79.83 | 79.80
ReB | 17.18 | 17.95 | 16.44 | 46.07 | 47.29 | 44.40 | 74.56 | 76.85 | 72.17
ReM | 15.53 | 15.53 | 15.50 | 41.18 | 40.94 | 41.16 | 66.75 | 66.81 | 66.86
| $p=10$
ReO | 19.54 | 19.85 | 19.04 | 49.14 | 49.11 | 49.15 | 79.59 | 79.60 | 79.66
ReB | 14.29 | 14.92 | 12.89 | 43.74 | 45.71 | 39.86 | 73.34 | 76.45 | 67.18
ReM | 10.21 | 10.21 | 10.28 | 31.98 | 32.03 | 32.03 | 54.03 | 53.99 | 54.17
| $p=20$
ReO | 19.35 | 19.88 | 20.36 | 50.30 | 49.70 | 49.99 | 79.91 | 80.04 | 79.84
ReB | 10.84 | 11.23 | 8.80 | 42.36 | 43.61 | 35.46 | 73.43 | 76.16 | 62.13
ReM | 6.66 | 6.68 | 6.66 | 24.61 | 24.45 | 24.46 | 42.02 | 42.02 | 41.95

Table 4 provides a detailed comparison of the average PRIV achieved by the
three competing re-randomization procedures under various simulation
complications. Given the insignificance of sample size $N$ in the ANOVA and
the greater challenge posed by larger $\sigma^{2}_{\boldsymbol{\beta}}$ for
Bayesian re-randomization procedures, we focus on reporting results for
settings where $N=200$ and $\sigma^{2}_{\boldsymbol{\beta}}=1$. Notable
findings include ReO consistently outperforming others and ReB surpassing
classic ReM, with its advantage increasing as $p$ grows.

Figure 4: The PRIV ratio of ReB over ReM when sample size $N$ is fixed at 600
and $\bar{R}^{2}$ fixed at 0.5 for different $p$ and $\rho$ with (A)
$\sigma^{2}_{{\boldsymbol{\beta}}}=0.1$ and (B)
$\sigma^{2}_{{\boldsymbol{\beta}}}=1$.

In Figure 4, we examine how ReB’s advantage over ReM varies with different
factors using the PRIV ratio of ReB to ReM. We maintain a fixed sample size of
600 and a constant $\bar{R}^{2}$ of 0.5 while varying $p$ and $\rho$. Two
scenarios are considered for $\sigma^{2}_{{\boldsymbol{\beta}}}$: 0.1 and 1,
and each box-plot summarizes results from 100 repeated experiments. The figure
shows that the majority of PRIV ratios are above 1 in all scenarios, even with
larger $\sigma^{2}_{{\boldsymbol{\beta}}}$, and exceptions become less common
as $\rho$ and $p$ increase.

Figure B.1 in the Supplement shows average PRIV for ReO and ReM across various
covariate dimensions ($p$) with different $\bar{R}^{2}$ values at $\rho=0$,
$\sigma^{2}_{\boldsymbol{\beta}}=0.1$ and $N=600$. Both methods show
substantial PRIV improvements, and ReM’s average PRIV decreases with higher
$p$, while ReO remains stable. These findings align with Theorem 2’s
predictions. Moreover, additional studies on acceptance probability precision
are detailed in Section B.2 of the Supplementary Material. As sample size $N$
grows, acceptance probabilities for the three procedures converge to
theoretical values, with ReO and ReB outperforming ReM.

###  5.2 Two-stage ReB versus ReM and ReO

In this subsection, we evaluate the performance of two-stage ReB procedures,
including BCRD-ReB, ReM-ReB, BCRD-ReO, and ReM-ReO, in comparison to ReO and
ReM. Our analysis focuses on simulated datasets with $\rho=0$,
$\sigma^{2}_{\boldsymbol{\beta}}=0$, $\bar{R}^{2}=0.5$, a first-stage sample
proportion $r\in{0.2,0.3,0.4}$, and a minimum sample size of $N\geq 200$,
excluding cases with $N=100$ due to potential two-stage ReB method limitations
in such scenarios. For a fair comparison, ReM’s acceptance probability is set
at $\alpha=0.05$, and other re-randomization methods’ acceptance probabilities
are adjusted numerically to ensure comparable computational costs, i.e., the
running time of sampling 1000 assignments. Each simulated dataset provides
PRIV estimates based on 1000 accepted allocations, and we calculate the
average PRIV for each simulation configuration. Detailed acceptance rates and
running times are documented in Section B.3 of the Supplementary Material.

Table 5 shows average PRIV for different methods in various simulation
scenarios, highlighting the best RPIV among ReM and two-stage procedures in
bold. Our observations are as follows: In most cases, ReM-ReB outperforms
other methods, except when the first-stage sample size is insufficient
relative to the covariate dimension ($p$) or when $p$ is small. Two-stage
methods exhibit a more pronounced advantage over ReM when dealing with high-
dimensional covariates. However, in scenarios with small sample sizes, BCRD-
ReO and ReM-ReO perform notably worse than BCRD-ReB and ReM-ReB. This suggests
that introducing variance through ReB improves causal effect estimation,
helping mitigate point estimator bias when sample sizes are limited.

Table 5: Average PRIV achieved by different re-randomization procedures when $\bar{R}^{2}=0.5$. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 48.08 | 48.08 | 48.08 | 48.90 | 48.90 | 48.90 | 48.55 | 48.55 | 48.55
ReM-ReB | 48.01 | 48.15 | 47.87 | 49.21 | 49.43 | 49.33 | 49.35 | 49.17 | 48.94
ReM-ReO | 47.08 | 46.94 | 47.59 | 49.03 | 48.82 | 48.95 | 49.11 | 49.07 | 49.24
BCRD-ReB | 38.94 | 34.39 | 29.57 | 40.04 | 34.20 | 30.98 | 39.63 | 34.84 | 29.64
BCRD-ReO | 36.95 | 32.78 | 28.52 | 39.80 | 34.20 | 30.79 | 39.55 | 34.40 | 29.72
| $p=5$
ReM | 41.27 | 41.27 | 41.27 | 41.94 | 41.94 | 41.94 | 42.56 | 42.56 | 42.56
ReM-ReB | 42.65 | 43.66 | 43.11 | 46.01 | 46.08 | 45.32 | 47.00 | 46.97 | 46.45
ReM-ReO | 37.92 | 41.16 | 41.73 | 45.22 | 45.22 | 44.49 | 46.81 | 46.76 | 46.07
BCRD-ReB | 37.36 | 33.12 | 28.74 | 38.68 | 34.16 | 29.63 | 39.78 | 35.38 | 30.93
BCRD-ReO | 31.44 | 29.37 | 26.46 | 36.37 | 33.62 | 29.16 | 39.32 | 34.97 | 30.11
| $p=10$
ReM | 32.74 | 32.74 | 32.74 | 33.83 | 33.83 | 33.83 | 34.48 | 34.48 | 34.48
ReM-ReB | 31.75 | 34.60 | 35.78 | 41.86 | 41.24 | 40.50 | 44.52 | 43.56 | 43.07
ReM-ReO | 21.18 | 28.35 | 32.29 | 40.02 | 40.45 | 39.59 | 43.67 | 43.10 | 42.75
BCRD-ReB | 28.67 | 28.13 | 25.66 | 35.99 | 31.77 | 28.25 | 38.58 | 34.07 | 29.64
BCRD-ReO | 17.20 | 21.62 | 20.90 | 33.38 | 30.16 | 27.08 | 37.26 | 33.66 | 28.69
| $p=20$
ReM | 23.59 | 23.59 | 23.59 | 26.41 | 26.41 | 26.41 | 26.26 | 26.26 | 26.26
ReM-ReB | $--$ | 22.32 | 25.05 | 34.50 | 35.39 | 34.43 | 40.11 | 38.69 | 37.42
ReM-ReO | $--$ | 14.10 | 20.45 | 29.13 | 32.63 | 33.44 | 38.45 | 38.10 | 37.54
BCRD-ReB | $--$ | 18.42 | 18.64 | 29.71 | 28.73 | 26.14 | 34.91 | 31.56 | 27.01
BCRD-ReO | $--$ | 9.82 | 12.69 | 24.80 | 25.57 | 23.49 | 32.44 | 30.38 | 26.74
Figure 5: Performance of BCRD-ReO, BCRD-ReB, ReM-ReO and ReM-ReB compared with
ReM with $r=0.2$. Sub-figure (A) explores the impact of factor $N$ with a
fixed $p=10$, while sub-figure (B) investigates factor $p$ with a fixed
$N=1200$.

To compare the performance of two-stage ReB procedures to ReM, Figure 5 shows
how different factors affect the PRIV ratio of BCRD-ReO, BCRD-ReB, ReM-ReO,
and ReM-ReB relative to ReM. With $r=0.2$ fixed, sub-figure (A) explores the
impact of factor $N$ with a fixed $p=10$, while sub-figure (B) investigates
factor $p$ with a fixed $N=1200$. The figure highlights that most PRIV ratios
of ReM-ReB over ReM are greater than 1, with exceptions diminishing as $N$ and
$p$ increase. In contrast, the other two-stage methods perform less
effectively, particularly when sample size or covariate dimensionality is
small.

In Section B.4 of the Supplement, we examine how varying sample proportions in
the first stage affect two-stage ReB procedures. We find that ReM-ReB with
around 20% of samples in the first stage consistently demonstrates the best
performance and robustness. Similar results also hold in more challenging
nonlinear scenarios, as shown in Section B.5 of the Supplementary Material. In
summary, simulations in both linear and nonlinear settings support the
superiority of the two-stage ReB approach over ReM, particularly when the
overall sample size is reasonably large (as indicated by Theorem 6). Among the
different two-stage ReB methods, ReM-ReB is generally the preferred choice in
most cases.

##  6 Real Data Application

In this section, we evaluate the performance of Bayesian re-randomization
procedures in a real-world application using a strategy based on Hill, (2011).
We employ covariates from the Infant Health and Development Program (IHDP), a
study investigating the impact of intensive child care and home visits on low-
birth-weight and premature infants. The IHDP dataset consists of 747 units and
25 pretreatment variables, including 6 continuous standardized covariates and
19 binary covariates. Here, we modify response surfaces based on Hill, (2011)
by adding more noise to response surface A. Specifically, we use the following
form for potential outcomes:

$\bm{Y}(0)\sim
N(\bm{X}{\boldsymbol{\beta}},3)\quad\text{and}\quad\bm{Y}(1)\sim
N(\bm{X}{\boldsymbol{\beta}}+4,3),$

where $\bm{X}$ is an expanded covariate matrix (747$\times$26) with first
column being a vector of ones, and ${\boldsymbol{\beta}}$ is a 26-dimensional
random coefficient vector with elements sampled from a discrete support
$\\{0,1,2,3,4\\}$ based on probability $(0.5,0.2,0.15,0.1,0.05)$.

To keep a balanced design, we exclude the last unit from the dataset, and
allocate the remaining 746 units into two groups of equal sample size (i.e.,
$N_{0}=N_{1}=373)$. Given the popularity of PCA-ReM and Ridge-ReM in high-
dimensional scenarios, we assess the performance of two-stage ReB procedures
(ReM-ReB, ReM-ReO, BCRD-ReB, and BCRD-ReO) alongside four different
randomization schemes: BCRD, ReM, PCA-ReM, and Ridge-ReM. We maintain ReM’s
acceptance probability at $\alpha=0.05$ while numerically tuning the
acceptance probabilities of other re-randomization procedures to match ReM’s
computational cost during the design phase. In PCA-ReM, we select principal
components capturing 95% of the variance. The running time ratio of these re-
randomization procedures relative to ReM is detailed in Section B.6 of the
Supplement. For the two-stage procedures, we vary the proportion of samples in
the first stage from 10% to 80% to assess their robustness. We calculate PRIV
for rerandomization procedures over BCRD based on 1000 accepted allocations
for each method over 100 iterations. The average PRIV for these methods is
presented in Figure 6.

Based on the figure, we can see the following facts. ReM-ReB and ReM-ReO
outperform ReM, PCA-ReM, and Ridge-ReM, even when allocating a substantial
proportion of samples to the first stage. Notably, Ridge-ReM fares poorly
under equivalent computational constraints. Additionally, ReM-ReB and ReM-ReO
exhibit greater robustness to changes in first-stage sample proportions
compared to BCRD-ReB and BCRD-ReO. Moreover, BCRD-ReO and ReM-ReO may struggle
with very small first-stage sample sizes. As a practical guideline, we
recommend using ReM-ReB with a 20% sample proportion in the first stage for
scenarios similar to this setting.

Figure 6: PRIV of ReM, BCRD-ReB, BCRD-ReO, ReM-ReB and ReM-ReO with different
proportion of sample in first stage $r$ in real application.

##  7 Conclusion and Discussions

In this study, we re-consider re-randomization from the Bayesian perspective.
Assuming that the importance vector of covariates ${\boldsymbol{\beta}}$ is
uncertain with respect to a prior distribution $\pi$, a Bayesian criterion is
established to derive a more efficient re-randomization procedure, referred to
as ReB, by minimizing the prior-integrated conditional asymptotic variance of
the casual effect estimator. We find that ReB only depends on the second
moment of prior distribution $\pi$, and enjoys a better efficiency than the
ReM-based re-randomization procedures in balancing high-dimensional covariates
and reducing the estimation variance of the difference-in-mean estimator for
causal effect, as long as the prior distribution is informative to highlight
the true importance vector. We also note that ReB is invariant to any scale
transformation of the prior distribution of ${\boldsymbol{\beta}}$, indicating
that re-scaling either $\pi$ or ${\boldsymbol{X}}$ would not change the re-
randomization procedure under the Bayesian criterion. Taking ReM, PCA-ReM,
Ridge-ReM as its special cases under different specifications of the prior
distribution $\pi$, ReB provides us a unified framework to understand and
interpret various re-randomization procedures. More importantly, ReB extends
the idea behind ReMT, and establishes a principled way to utilize prior
information about ${\boldsymbol{\beta}}$ in re-randomization.

Additionally, the Bayesian nature of ReB also provides an additional avenue
for designing more efficient experiments in sequential trials. The success of
the two-stage ReB gives us an important message: once we can observe the
responses of some units before the arrival of other experimental units, we
would be able to design a more efficient randomized experiment by following
the spirit of ReB. Now, suppose a typical scenario where experiment units
arrive in a sequential fashion with their potential outcomes observed along
the experiment process so that we can update our knowledge about
${\boldsymbol{\beta}}$ gradually. In such cases, we would be able to adjust
the measurement for covariate imbalance simultaneously, and thus improve the
overall balance of covariates more efficiently by carefully allocating
treatment assignments to upcoming units based on the adjusted imbalance
measurement induced by the Bayesian principle for re-randomization. Compared
to the classic approaches for sequential design, which typically do not take
into consideration the responses of units received along the experiment
process, the above alternative approach is apparently more attractive. Due to
the limited space of this paper, however, we reserve this exciting topic for
future research.

##  8 Funding

This work was supported by National Science Foundation of China under Grant
No. 11931001.

## References

  * Branson et al., (2016) Branson, Z., Dasgupta, T., and Rubin, D. B. (2016).  Improving covariate balance in 2k factorial designs via rerandomization.  The Annals of Applied Statistics, 10(4):1958–1976.
  * Branson and Shao, (2021) Branson, Z. and Shao, S. (2021).  Ridge rerandomization: an experimental design strategy in the presence of collinearity.  Journal of Statistical Planning and Inference, 211:287–314.
  * Cox, (2009) Cox, D. R. (2009).  Randomization in the design of experiments.  International Statistical Review, 77(3):415–429.
  * Fisher, (1992) Fisher, R. A. (1992).  The arrangement of field experiments.  Breakthroughs in statistics: Methodology and distribution, pages 82–91.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Holschuh, (1980) Holschuh, N. (1980).  Randomization and design: I.  In R.A. Fisher: An Appreciation, pages 35–45, New York, NY. Springer New York.
  * Imai et al., (2008) Imai, K., King, G., and Stuart, E. A. (2008).  Misunderstanding between experimentalists and observationalists about causal inference.  Journal of the Royal Statistical Society Series A, 171(2):481–502.
  * Imbens and Rubin, (2015) Imbens, G. W. and Rubin, D. B. (2015).  Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction.  Cambridge University Press.
  * Johansson et al., (2021) Johansson, P., Rubin, D. B., and Schultzberg, M. (2021).  On optimal rerandomization designs.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(2):395–403.
  * Johansson and Schultzberg, (2020) Johansson, P. and Schultzberg, M. (2020).  Rerandomization strategies for balancing covariates using pre-experimental longitudinal data.  Journal of Computational and Graphical Statistics, 29(4):798–813.
  * Kallus, (2018) Kallus, N. (2018).  Optimal a priori balance in the design of controlled experiments.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):85–112.
  * Kallus, (2021) Kallus, N. (2021).  On the optimality of randomization in experimental design: How to randomize for minimax variance and design-based inference.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(2):404–409.
  * Li and Ding, (2017) Li, X. and Ding, P. (2017).  General forms of finite population central limit theorems with applications to causal inference.  Journal of the American Statistical Association, 112(520):1759–1769.
  * Li et al., (2018) Li, X., Ding, P., and Rubin, D. B. (2018).  Asymptotic theory of rerandomization in treatment-control experiments.  Proceedings of the National Academy of Sciences, 115(37):9157–9162.
  * Li et al., (2020) Li, X., Ding, P., and Rubin, D. B. (2020).  Rerandomization in 2K factorial experiments.  The Annals of Statistics, 48(1):43–63.
  * Li et al., (2021) Li, Y., Kang, L., and Huang, X. (2021).  Covariate balancing based on kernel density estimates for controlled experiments.  Statistical Theory and Related Fields, 5(2):102–113.
  * Morgan and Rubin, (2012) Morgan, K. L. and Rubin, D. B. (2012).  Rerandomization to improve covariate balance in experiments.  The Annals of Statistics, 40(2):1263–1282.
  * Morgan and Rubin, (2015) Morgan, K. L. and Rubin, D. B. (2015).  Rerandomization to balance tiers of covariates.  Journal of the American Statistical Association, 110(512):1412–1421.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688–701.
  * Rubin, (1980) Rubin, D. B. (1980).  Randomization analysis of experimental data: The fisher randomization test comment.  Journal of the American Statistical Association, 75(371):591–593.
  * Urbach, (1985) Urbach, P. (1985).  Randomization and the design of experiments.  Philosophy of Science, 52(2):256–273.
  * Wang et al., (2021) Wang, X., Wang, T., and Liu, H. (2021).  Rerandomization in stratified randomized experiments.  Journal of the American Statistical Association, pages 1–10.
  * Wang and Li, (2022) Wang, Y. and Li, X. (2022).  Rerandomization with diminishing covariate imbalance and diverging number of covariates.  The Annals of Statistics, 50(6):3439–3465.
  * Wu, (1981) Wu, C.-F. (1981).  On the robustness and efficiency of some randomized designs.  The Annals of Statistics, pages 1168–1177.
  * Yang et al., (2021) Yang, Z., Qu, T., and Li, X. (2021).  Rejective sampling, rerandomization, and regression adjustment in survey experiments.  Journal of the American Statistical Association, pages 1–15.
  * (26) Zhang, H., Yin, G., et al. (2021a).  Response-adaptive rerandomization.  Journal of the Royal Statistical Society Series C, 70(5):1281–1298.
  * (27) Zhang, H., Yin, G., and Rubin, D. B. (2021b).  PCA rerandomization.  https://arxiv.org/abs/2102.12262.
  * Zhou et al., (2018) Zhou, Q., Ernst, P. A., Morgan, K. L., Rubin, D. B., and Zhang, A. (2018).  Sequential rerandomization.  Biometrika, 105(3):745–752.
  * Zhu and Liu, (2022) Zhu, K. and Liu, H. (2022).  Pair-switching rerandomization.  Biometrics.

Supplementary Material

In this document, Appendix A contains technical proofs, and Appendix B shows
additional simulation results.

##  Appendix A Technical Proofs

###  A.1 Useful Lemmas for Theorem 1

To prove Theorem 1, we propose the following two lemmas.

######  Lemma 1.

Let $\bm{Z}$ be a random vector defined on a probability space
$(\Omega,\mathcal{F},\mathcal{P})$, $g(\bm{Z})$ be a measurable function of
$\bm{Z}$ satisfying ${\mathbb{E}}|g(\bm{Z})|<\infty$, and

$\Psi_{\alpha}=\left\\{\phi(\bm{Z})=0\text{ or
}1:{\mathbb{P}}(\phi(\bm{Z})=1)=\alpha\right\\}$

be the set of re-randomization mechanisms for $\bm{Z}$ with acceptance rate
$\alpha$. Define

$\begin{split}\phi_{g,\alpha}(\bm{Z})\triangleq I\left(g(\bm{Z})\leq
a\right)\text{ with
}{\mathbb{P}}(\phi_{g,\alpha}(\bm{Z})=1)={\mathbb{P}}(g(\bm{Z})\leq
a)=\alpha\end{split}$

as re-randomization mechanism that defines the acceptance region based on the
contour line of $g(\cdot)$. We have:

$\begin{split}\phi_{g,\alpha}=\arg\min_{\phi\in{\Psi_{\alpha}}}{\mathbb{E}}\big{(}g(\bm{Z})|\phi(\bm{Z})=1\big{)}.\end{split}$

######  Lemma 2.

Suppose Condition 1 holds and $\bm{\beta}$ is known. Under complete
randomization, we have

$\begin{split}d_{\bm{\beta}}/\sigma^{2}_{\bm{\beta}}\xrightarrow{d}\chi^{2}_{1}\
\text{ as }N\rightarrow\infty,\end{split}$

where $d_{{\boldsymbol{\beta}}}$ and $\sigma^{2}_{\bm{\beta}}$ are defined in
Theorem 1.

###  A.2 Proof of Lemma 1

Assume that there exists another $\phi^{*}(\bm{Z})\in\Psi_{\alpha}$ that
satisfies $\phi^{*}\neq\phi_{g,\alpha}$ and

${\mathbb{E}}(g({\boldsymbol{Z}})\mid\phi^{*}(\bm{Z})=1)<{\mathbb{E}}(g({\boldsymbol{Z}})\mid\phi_{g,\alpha}(\bm{Z})=1).$

Let
$A\triangleq\left\\{\omega:\phi_{g,\alpha}(\bm{Z}(w))=1\right\\}\in\mathcal{F}$,
$A^{*}\triangleq\left\\{\omega:\phi^{*}(\bm{Z}(w))=1\right\\}\in\mathcal{F}$,
then ${\boldsymbol{Z}}$ is uniformly distributed on $A^{*}$ and $A$ under
$\phi^{*}$ and $\phi$, and we have

$\int_{A^{*}}g(\bm{Z})dF_{{\boldsymbol{Z}}}<\int_{A}g(\bm{Z})dF_{{\boldsymbol{Z}}}.$

Considering that

$\displaystyle\int_{A^{*}}g(\bm{Z})dF_{{\boldsymbol{Z}}}$ $\displaystyle=$
$\displaystyle\int_{A^{*}\cap A}g(\bm{Z})dF_{{\boldsymbol{Z}}}+\int_{A^{*}\cap
A^{c}}g(\bm{Z})dF_{{\boldsymbol{Z}}},$
$\displaystyle\int_{A}g(\bm{Z})dF_{{\boldsymbol{Z}}}$ $\displaystyle=$
$\displaystyle\int_{A^{*}\cap A}g(\bm{Z})dF_{{\boldsymbol{Z}}}+\int_{A\cap
A^{*c}}g(\bm{Z})dF_{{\boldsymbol{Z}}},$

and ${\mathbb{E}}|g(\bm{Z})|<\infty$, we have

$\int_{A^{*}\cap A^{c}}g(\bm{Z})dF_{{\boldsymbol{Z}}}<\int_{A\cap
A^{*c}}g(\bm{Z})dF_{{\boldsymbol{Z}}}.$

On the other hand, because $A=\left\\{\omega:g(\bm{Z}(w))\leq a\right\\}$, we
have

$\displaystyle\int_{A^{*}\cap A^{c}}g(\bm{Z})dF_{{\boldsymbol{Z}}}$
$\displaystyle>$ $\displaystyle a\cdot\mathbb{P}(A^{*}\cap A^{c}),$
$\displaystyle\int_{A\cap A^{*c}}g(\bm{Z})dF_{{\boldsymbol{Z}}}$
$\displaystyle\leq$ $\displaystyle a\cdot\mathbb{P}(A\cap A^{*c}).$

Considering that the acceptance rates of $\phi_{g}$ and $\phi^{*}$ are the
same, we have

$\mathbb{P}(A^{*}\cap A^{c})=\mathbb{P}(A\cap A^{*c}).$

Therefore,

$\displaystyle\int_{A^{*}\cap A^{c}}g(\bm{Z})dF_{{\boldsymbol{Z}}}>\int_{A\cap
A^{*c}}g(\bm{Z})dF_{{\boldsymbol{Z}}},$

and we reach a contradiction. $\blacksquare$

###  A.3 Proof of Lemma 2

According to the finite population central limit theorem, under regularity
conditions in Condition 1, the large sampling distribution, over all
randomizations, of $\sqrt{N}{\boldsymbol{D}}$ is asymptotically Gaussian:

$\begin{split}\sqrt{N}{\boldsymbol{D}}\overset{\cdot}{\sim}\bm{B},\
\bm{B}\sim\mathcal{N}\left(\bm{0},{\boldsymbol{V}}_{xx}\right),\end{split}$

where $\overset{\cdot}{\sim}$ means two random variables converge to the same
distribution weakly. Hence,

$\begin{split}\sqrt{N}{\boldsymbol{D}}^{T}{\boldsymbol{\beta}}\overset{\cdot}{\sim}\bm{B}^{T}{\boldsymbol{\beta}},\text{
where
}\bm{B}^{T}{\boldsymbol{\beta}}\sim\mathcal{N}\left(\bm{0},{\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}\right).\end{split}$

Continuous mapping theorem and Slutsky’s theorem imply that

$\begin{split}\frac{N{\boldsymbol{D}}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}{\boldsymbol{D}}}{{\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}}\overset{\cdot}{\sim}\frac{\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}}{{\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}},\text{
where
}\frac{\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}}{{\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}}\sim\chi^{2}_{1}.\end{split}$

Therefore,

$\begin{split}d_{\bm{\beta}}/\sigma^{2}_{\bm{\beta}}\xrightarrow{d}\chi^{2}_{1},\end{split}$

where $\sigma^{2}_{\bm{\beta}}=\bm{\beta}^{T}{\boldsymbol{V}}_{xx}\bm{\beta}$.
$\blacksquare$

###  A.4 Proof of Theorem 1

First of all, Lemma 2 implies that the asymptotic acceptance rate of
$\phi_{\bm{\beta}}$ is

$\begin{split}\gamma_{\phi_{\bm{\beta}}}=\lim_{N\rightarrow\infty}{\mathbb{P}}\big{(}d_{\bm{\beta}}/\sigma^{2}_{\bm{\beta}}\leq\xi_{\alpha,1}\big{)}=\alpha.\end{split}$

Hence, $\phi_{\bm{\beta}}\in\Phi_{\alpha}$. Since the first term in (4), i.e.,
$(1-R^{2}_{\infty})V_{\tau\tau,\infty}$, is a constant that does not depend on
the specification of the re-randomization mechanism $\phi$, and the asymptotic
sampling variance in (4) is determined by the interaction of two factors
$\bm{\beta}_{\infty}$ and ${\boldsymbol{B}}_{\infty}$ in the second term.
Therefore, minimizing
${\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}$ is equivalent
to minimizing

$\begin{split}{\mathbb{E}}\left(\bm{B}^{T}_{\infty}\bm{\beta}_{\infty}\bm{\beta}_{\infty}^{T}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right)\end{split}$

within $\phi\in\Phi_{\alpha}$. Direct application of Lemma 1 implies the
asymptotic acceptance region of the optimal criterion should take the form of

$\mathcal{B}_{\phi,\infty}=\big{\\{}\bm{\mu}:\bm{\mu}^{T}\bm{\beta}_{\infty}\bm{\beta}_{\infty}^{T}\bm{\mu}\leq
a\big{\\}},$ (A.1)

where the threshold $a$ is determined by solving

$\begin{split}{\mathbb{P}}\left(\bm{B}_{\infty}^{T}\bm{\beta}_{\infty}\bm{\beta}_{\infty}^{T}\bm{B}_{\infty}\leq
a\right)=\alpha.\end{split}$

Since
${\boldsymbol{B}}_{\infty}\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{V}}_{xx,\infty})$,
we have
$a=\bm{\beta}_{\infty}^{T}\bm{V}_{xx,\infty}\bm{\beta}_{\infty}\cdot\xi_{\alpha,1}$.
As the acceptance region of $\phi_{\bm{\beta}}$ is

$\begin{split}\mathcal{B}_{\phi_{\bm{\beta}}}=\left\\{\bm{\mu}:\bm{\mu}^{T}\bm{\beta}\bm{\beta}^{T}\bm{\mu}\leq\bm{\beta}^{T}\bm{V}_{xx}\bm{\beta}\cdot\xi_{\alpha,1}\right\\},\end{split}$

whose limit is consistent with the form of (A.1). Therefore,

$\phi_{{\boldsymbol{\beta}}}={\operatorname{arg\,min}}_{\phi\in\Phi_{\alpha}}{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}.$

$\blacksquare$

###  A.5 Proof of Theorem 2

Under Condition 1,

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{{\boldsymbol{\beta}}}=1&\overset{\cdot}{\sim}\epsilon+\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{{\boldsymbol{\beta}}}}\\\
&\sim\epsilon+{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}\mid({\boldsymbol{\beta}}^{T}{\boldsymbol{B}})^{2}\leq\sigma^{2}_{{\boldsymbol{\beta}}}\cdot\xi_{\alpha,1},\end{split}$

where $\epsilon\sim\mathcal{N}(0,V_{\tau\tau}-\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau})$ is independent of
$\bm{B}\sim\mathcal{N}(0,{\boldsymbol{V}}_{xx})$, and
$\sigma^{2}_{{\boldsymbol{\beta}}}={\text{var}}({\boldsymbol{\beta}}^{T}{\boldsymbol{B}})=\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau}$. Let
$\tilde{{\boldsymbol{B}}}=\sigma^{-1}_{{\boldsymbol{\beta}}}{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}$,
we have $\tilde{{\boldsymbol{B}}}\sim\mathcal{N}(0,1)$, and

$\begin{split}{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}\mid({\boldsymbol{\beta}}^{T}{\boldsymbol{B}})^{2}\leq\sigma^{2}_{{\boldsymbol{\beta}}}\cdot\xi_{\alpha,1}&\sim\sigma_{{\boldsymbol{\beta}}}\cdot\tilde{{\boldsymbol{B}}}\mid\tilde{{\boldsymbol{B}}}^{2}\leq\xi_{\alpha,1}.\end{split}$

Since $\sigma_{{\boldsymbol{\beta}}}^{2}=V_{\tau\tau}\cdot R^{2}$, we have

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{{\boldsymbol{\beta}}}=1&\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}(1-R^{2})}\cdot\varepsilon_{0}+\sqrt{V_{\tau\tau}\cdot
R^{2}}\cdot(D\mid D^{2}\leq\xi_{\alpha,1}),\end{split}$

where $\varepsilon_{0},D\sim\mathcal{N}(0,1)$ are independent. Moreover, by
Eq. (8),

$\begin{split}r_{\phi_{{\boldsymbol{\beta}}}}=\frac{{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi_{{\boldsymbol{\beta}}},\infty}\right)}{{\mathbb{E}}\left(\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\right)}=\frac{{\mathbb{E}}\left(\big{(}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}\big{)}^{2}|\big{(}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}\big{)}^{2}\leq{\boldsymbol{\beta}}^{T}_{\infty}\bm{V}_{xx,\infty}{\boldsymbol{\beta}}_{\infty}\cdot\xi_{\alpha,1}\right)}{{\mathbb{E}}\left(\big{(}{\boldsymbol{\beta}}^{T}_{\infty}\bm{B}_{\infty}\big{)}^{2}\right)},\end{split}$

where
$\left({\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\right)^{2}\sim{\boldsymbol{\beta}}^{T}_{\infty}\bm{V}_{xx,\infty}{\boldsymbol{\beta}}_{\infty}\cdot\chi_{1}^{2}$.
Hence, $r_{\phi_{{\boldsymbol{\beta}}}}=v_{\alpha,1}$ and

$\text{PRIASV}_{ReO}=100\times(1-v_{\alpha,1})R^{2}_{\infty}.$

Recall that ReM achieves

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{{\boldsymbol{\beta}}}=1\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}}\big{(}\sqrt{1-R^{2}}\cdot\varepsilon_{0}+\sqrt{R^{2}}\cdot
L_{1,\xi_{\alpha,p}}\big{)},\end{split}$

and

$\begin{split}\text{PRIASV}_{ReM}=100\times(1-v_{\alpha,p})R^{2}_{\infty},\end{split}$

where
$v_{\alpha,p}={\mathbb{E}}\left(\chi^{2}_{p}|\chi^{2}_{p}\leq\xi_{\alpha,p}\right)/p$.
Therefore,

$\frac{\text{PRIASV}_{ReM}}{\text{PRIASV}_{ReO}}=\frac{1-v_{\alpha,p}}{1-v_{\alpha,1}}.$

From the definition of $v_{\alpha,p}$, we have

$v_{\alpha,p}=\frac{{\mathbb{P}}(\chi^{2}_{p+2}\leq\xi_{\alpha,p})}{{\mathbb{P}}(\chi^{2}_{p}\leq\xi_{\alpha,p})}=\frac{1}{\alpha}\cdot\frac{\int_{0}^{\xi_{\alpha,p}}{x^{\frac{p}{2}}e^{-\frac{x}{2}}}dx}{2^{\frac{p+2}{2}}\cdot\Gamma(\frac{p+2}{2})}.$

Based on the integration by parts,

$\displaystyle\int_{0}^{\xi_{\alpha,p}}x^{\frac{p}{2}}e^{-\frac{x}{2}}dx$
$\displaystyle=$
$\displaystyle-2\left[x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}-\frac{p}{2}\int_{0}^{\xi_{\alpha,p}}x^{\frac{p}{2}-1}e^{-\frac{x}{2}}dx\right]$
$\displaystyle=$ $\displaystyle-2\cdot
x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}+p\cdot\frac{\int_{0}^{\xi_{\alpha,p}}x^{\frac{p}{2}-1}e^{-\frac{x}{2}}dx}{2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2})}\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2})$ $\displaystyle=$
$\displaystyle-2\cdot
x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}+p\cdot{\mathbb{P}}(\xi^{2}_{p}\leq\xi_{\alpha,p})\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2})$ $\displaystyle=$
$\displaystyle-2\cdot
x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}+p\cdot\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2}),$

indicating that

$\displaystyle v_{\alpha,p}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}\cdot\frac{-2\cdot
x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}+p\cdot\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2})}{2^{\frac{p+2}{2}}\cdot\Gamma(\frac{p+2}{2})}$
$\displaystyle=$
$\displaystyle-\frac{x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}}{\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p+2}{2})}+\frac{p\cdot\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p}{2})}{\alpha\cdot
2^{\frac{p+2}{2}}\cdot\Gamma(\frac{p+2}{2})}$ $\displaystyle=$
$\displaystyle-\frac{x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}}{\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p+2}{2})}+1,$

and thus,

$1-v_{\alpha,p}=\frac{x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}}{\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p+2}{2})}.$

Considering that $\chi_{p}^{2}$ is close to normal distribution $N(p,2p)$ when
$p$ is large, we have

$\xi_{\alpha,p}=p+\sqrt{2p}\cdot\Phi^{-1}(\alpha)+o(1),$

and thus,

$\displaystyle\frac{x^{\frac{p}{2}}\cdot
e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}}{p^{\frac{p}{2}}\cdot e^{-\frac{p}{2}}}$
$\displaystyle=$
$\displaystyle\left(1+\frac{\sqrt{2}\cdot\Phi^{-1}(\alpha)}{\sqrt{p}}+o\left(p^{-1}\right)\right)^{\frac{p}{2}}\cdot\exp\left\\{-\frac{\sqrt{p}\cdot\Phi^{-1}(\alpha)+o(1)}{\sqrt{2}}\right\\}$
$\displaystyle\rightarrow$
$\displaystyle\exp\left\\{\frac{\sqrt{p}\cdot\Phi^{-1}(\alpha)}{\sqrt{2}}\right\\}\cdot\exp\left\\{-\frac{\sqrt{p}\cdot\Phi^{-1}(\alpha)}{\sqrt{2}}\right\\}=1,$

i.e., $x^{\frac{p}{2}}\cdot
e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}=O\left(p^{\frac{p}{2}}\cdot
e^{-\frac{p}{2}}\right)$. On the other hand, because

$\Gamma\left(\frac{p}{2}+1\right)=\begin{cases}\frac{p}{2}!,&p\mbox{ is
even},\\\ \frac{p!\sqrt{\pi}}{2^{p}(\frac{p-1}{2})!},&p\mbox{ is
odd},\end{cases}$

and $p!=O(p^{p+\frac{1}{2}}\cdot e^{-p})$ based on the Stirling equation, we
have

$\Gamma\left(\frac{p}{2}+1\right)=O\left(\left(\frac{p}{2}\right)^{\frac{p+1}{2}}\cdot
e^{-\frac{p}{2}-1}\right).$

And, thus

$1-v_{\alpha,p}=\frac{x^{\frac{p}{2}}e^{-\frac{x}{2}}|_{x=\xi_{\alpha,p}}}{\alpha\cdot
2^{\frac{p}{2}}\cdot\Gamma(\frac{p+2}{2})}=\frac{O\left(p^{\frac{p}{2}}\cdot
e^{-\frac{p}{2}}\right)}{\alpha\cdot 2^{\frac{p}{2}}\cdot
O\left(\big{(}\frac{p}{2}\big{)}^{\frac{p+1}{2}}\cdot
e^{-\frac{p}{2}-1}\right)}=O\left(\frac{1}{\sqrt{p}}\right).$

Because $(1-v_{\alpha,1})$ is constant, we finally have

$\frac{\text{PRIASV}_{ReM}}{\text{PRIASV}_{ReO}}=\frac{1-v_{\alpha,p}}{1-v_{\alpha,1}}=O\left(\frac{1}{\sqrt{p}}\right).$

$\blacksquare$

###  A.6 Useful Lemmas for Theorem 3

To prove Theorem 3, we propose the following two lemmas.

######  Lemma 3.

Under Condition 1 and complete randomization, we have

$d_{\pi}\overset{\cdot}{\sim}\chi^{2}_{{\boldsymbol{\lambda}}}=\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2},$
(A.2)

where $Z_{1},...,Z_{p}$ are i.i.d. random variables from standard Gaussian
distribution, and vector
${\boldsymbol{\lambda}}=(\lambda_{1},\ldots,\lambda_{p})$ contains eigenvalues
of matrix
$\bm{P}=\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}\bm{\Lambda}_{\pi}\bm{V}_{xx}^{1/2}$
with $\bm{V}_{xx}^{1/2}$ being the Cholesky square root of $\bm{V}_{xx}$
satisfying $\bm{V}_{xx}^{1/2}\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}=\bm{V}_{xx}$.

######  Lemma 4.

Suppose random variables $X_{n}$ and $Y_{n}$ converge weakly to the same
distribution as $X$, i.e.,

$\begin{split}X_{n}\xrightarrow{d}X,\ Y_{n}\xrightarrow{d}X,\end{split}$

and the distribution function $F_{X}$ is continuous. Let $\xi_{Y_{n},\alpha}$
be the $\alpha$-quantile of the distribution of $Y_{n}$, then

$\begin{split}{\mathbb{P}}\left(X_{n}\leq\xi_{Y_{n},\alpha}\right)\rightarrow\alpha.\end{split}$

###  A.7 Proof of Lemma 3

Suppose the Cholesky decomposition of
${\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}$ is
${\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}={\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{1/2}\big{(}{\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{1/2}\big{)}^{T}$
and let
$\tilde{{\boldsymbol{Z}}}={\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{-1/2}{\boldsymbol{D}}$.
Then we have
${\boldsymbol{D}}={\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{1/2}\tilde{{\boldsymbol{Z}}}$
and

$d_{\pi}=N{\boldsymbol{D}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{D}}=N\tilde{{\boldsymbol{Z}}}^{T}\left(\big{(}{\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{1/2}\right)\tilde{{\boldsymbol{Z}}}=\tilde{{\boldsymbol{Z}}}^{T}\bm{P}\tilde{{\boldsymbol{Z}}},$

where
$\bm{P}=N\big{(}{\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{\Sigma}}_{\boldsymbol{D}}^{1/2}=\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{V}}_{xx}^{1/2}$
is a symmetric positive-definite real matrix with the following eigenvalue
decomposition

$\bm{P}={\boldsymbol{\Gamma}}diag\\{\lambda_{1},\cdots,\lambda_{p}\\}{\boldsymbol{\Gamma}}^{T}\quad\mbox{with}\quad{\boldsymbol{\Gamma}}^{T}{\boldsymbol{\Gamma}}={\boldsymbol{I}}_{p}.$

Central limit theorem of finite population implies that
$\sqrt{N}{\boldsymbol{D}}$ is asymptotically Gaussian under Condition 1:

$\begin{split}\sqrt{N}{\boldsymbol{D}}\overset{\cdot}{\sim}\bm{B},\
\bm{B}\sim\mathcal{N}(\bm{0},\bm{V}_{xx}).\end{split}$

Thus,

$\tilde{{\boldsymbol{Z}}}={\boldsymbol{\Sigma}}_{{\boldsymbol{D}}}^{-1/2}{\boldsymbol{D}}\overset{\cdot}{\sim}{\boldsymbol{V}}_{xx}^{-1/2}\bm{B}\sim\mathcal{N}(\bm{0},\bm{I}),\
\
{\boldsymbol{\Gamma}}^{T}\tilde{{\boldsymbol{Z}}}\sim\mathcal{N}(\bm{0},\bm{I}).$

Therefore,

$d_{\pi}=\tilde{{\boldsymbol{Z}}}^{T}\bm{P}\tilde{{\boldsymbol{Z}}}\overset{\cdot}{\sim}\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2},$

where $Z_{j}$’s are i.i.d. standard normal random variables and weights
$\lambda_{j}$’s are eigenvalues of matrix $\bm{P}$. $\blacksquare$

###  A.8 Proof of Lemma 4

Let $\xi_{X,a}$ be the $\alpha$-quantile of distribution of random variable
$X$, $F_{X_{n}}(\cdot),F_{Y_{n}}(\cdot),F_{X}(\cdot)$ be the distribution
function of $X_{n}$, $Y_{n}$ and $X$ respectively. Polya’s theorem implies

$\begin{split}&\lim_{n\rightarrow\infty}\sup_{x\in\mathcal{R}}\left|F_{X_{n}}(x)-F_{X}(x)\right|=0,\\\
&\lim_{n\rightarrow\infty}\sup_{x\in\mathcal{R}}\left|F_{Y_{n}}(x)-F_{X}(x)\right|=0.\end{split}$

Hence,

$\left|{\mathbb{P}}(X_{n}\leq\xi_{Y_{n},\alpha})-{\mathbb{P}}(X\leq\xi_{Y_{n},\alpha})\right|\rightarrow
0.$

Since

$\begin{split}\left|{\mathbb{P}}(X\leq\xi_{Y_{n},\alpha})-{\mathbb{P}}(X\leq\xi_{X,\alpha})\right|&=\left|F_{X}(\xi_{Y_{n},\alpha})-F_{X}(\xi_{X,\alpha})\right|\\\
&\leq\left|F_{X}(\xi_{Y_{n},\alpha})-F_{Y_{n}}(\xi_{Y_{n},\alpha})\right|+\left|F_{Y_{n}}(\xi_{Y_{n},\alpha})-F_{X}(\xi_{X,\alpha})\right|\\\
&=\left|F_{X}(\xi_{Y_{n},\alpha})-F_{Y_{n}}(\xi_{Y_{n},\alpha})\right|\rightarrow
0,\end{split}$

we have

$\begin{split}\left|{\mathbb{P}}(X_{n}\leq\xi_{Y_{n},\alpha})-{\mathbb{P}}(X\leq\xi_{X,\alpha})\right|&\leq\left|{\mathbb{P}}(X_{n}\leq\xi_{Y_{n},\alpha})-{\mathbb{P}}(X\leq\xi_{Y_{n},\alpha})\right|\\\
&+\left|{\mathbb{P}}(X\leq\xi_{Y_{n},\alpha})-{\mathbb{P}}(X\leq\xi_{X,\alpha})\right|\rightarrow
0.\end{split}$

That is,

$\begin{split}{\mathbb{P}}(X_{n}\leq\xi_{Y_{n},\alpha})\rightarrow\alpha.\end{split}$

$\blacksquare$

###  A.9 Proof of Theorem 3

First of all, Lemma 3 and Lemma 4 imply that the asymptotic acceptance rate of
$\phi_{\pi}$ is

$\begin{split}\gamma_{\phi_{\pi}}=\lim_{N\rightarrow\infty}{\mathbb{P}}\big{(}d_{\pi}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}=\alpha.\end{split}$

Hence, $\phi_{\pi}\in\Phi_{\alpha}$. Considering that
${\mathbb{E}}_{\pi}\left[(1-R^{2}_{\infty})V_{\tau\tau,\infty}\right]$ is
independent of $\phi$, minimizing
${\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}$ is
equivalent to minimizing

$\begin{split}{\mathbb{E}}_{\pi}\Big{[}{\mathbb{E}}\left(\bm{B}_{\infty}^{T}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}|\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right)\Big{]}&={\mathbb{E}}\left[\bm{B}_{\infty}^{T}{\mathbb{E}}_{\pi}\left({\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\right)\bm{B}_{\infty}\mid\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right]\\\
&={\mathbb{E}}\left[\bm{B}_{\infty}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{B}_{\infty}\mid\bm{B}_{\infty}\in\mathcal{B}_{\phi,\infty}\right],\end{split}$
(A.3)

where

$\bm{\Lambda}_{\pi}={\mathbb{E}}_{\pi}\left({\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\right)=\bm{\mu}_{\pi}\bm{\mu}_{\pi}^{T}+\bm{\Sigma}_{\pi}$

is the second moment of $\pi$ and is referred to as the characteristic matrix
of the prior distribution, with $\bm{\mu}_{\pi}$ and $\bm{\Sigma}_{\pi}$ being
the expectation vector and covariance matrix of $\pi$, respectively. Direct
application of Lemma 1 implies the asymptotic acceptance region of the optimal
criterion should take the form of

$\mathcal{B}_{\phi,\infty}=\big{\\{}\bm{\mu}:\bm{\mu}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{\mu}\leq
a\big{\\}},$ (A.4)

where the threshold $a$ is determined by solving

$\begin{split}{\mathbb{P}}\left(\bm{B}_{\infty}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{B}_{\infty}\leq
a\right)=\alpha.\end{split}$

Since
${\boldsymbol{B}}_{\infty}\sim\mathcal{N}({\boldsymbol{0}},{\boldsymbol{V}}_{xx,\infty})$,
we have $a=\xi_{\alpha,{\boldsymbol{\lambda}}_{\infty}}$, where
${\boldsymbol{\lambda}}_{\infty}$ is vector of eigenvalues of
$\big{(}\bm{V}_{xx,\infty}^{1/2}\big{)}^{T}\bm{\Lambda}_{\pi}\bm{V}_{xx,\infty}^{1/2}$,
and $Z_{1},...,Z_{p}$ being i.i.d. random variables from standard Gaussian
distribution. As the acceptance region of $\phi_{\pi}$ is

$\begin{split}\mathcal{B}_{\phi_{\pi}}=\left\\{\bm{\mu}:\bm{\mu}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{\mu}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right\\},\end{split}$

whose limit is consistent with the form of (A.4). Therefore,

$\phi_{\pi}={\operatorname{arg\,min}}_{\phi\in\Phi_{\alpha}}{\text{var}}_{\pi,a}\big{(}\sqrt{N}(\hat{\tau}-\tau)|\phi=1\big{)}.$

$\blacksquare$

###  A.10 Proof of Corollary 1

Based on Theorem 3, we have

$\displaystyle\phi_{\pi_{1}}=I\left(d_{\pi_{1}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}_{1}}\right),$
$\displaystyle\phi_{\pi_{2}}=I\left(d_{\pi_{2}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}_{2}}\right),$

where ${\boldsymbol{\lambda}}_{1}$ and ${\boldsymbol{\lambda}}_{1}$ are
vectors of eigenvalues of
$\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi_{1}}\bm{V}_{xx}^{1/2}$
and
$\big{(}\bm{V}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi_{2}}\bm{V}_{xx}^{1/2}$,
respectively. Since
${\boldsymbol{\Lambda}}_{\pi_{2}}=r{\boldsymbol{\Lambda}}_{\pi_{1}}$, we have

$d_{\pi_{2}}=N{\boldsymbol{D}}^{T}{\boldsymbol{\Lambda}}_{\pi_{2}}{\boldsymbol{D}}=rN{\boldsymbol{D}}^{T}{\boldsymbol{\Lambda}}_{\pi_{1}}{\boldsymbol{D}}=r\cdot
d_{\pi_{1}},\
\xi_{\alpha,{\boldsymbol{\lambda}}_{2}}=r\cdot\xi_{\alpha,{\boldsymbol{\lambda}}_{1}},$

and thus

$\phi_{\pi_{2}}=I\left(d_{\pi_{2}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}_{2}}\right)=I\left(r\cdot
d_{\pi_{1}}\leq
r\cdot\xi_{\alpha,{\boldsymbol{\lambda}}_{1}}\right)=\phi_{\pi_{1}}.$

The proof is complete. $\blacksquare$

###  A.11 Proof of Theorem 4

Under Condition 1,

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1&\overset{\cdot}{\sim}\epsilon+\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\\\
&\sim\epsilon+{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}\mid{\boldsymbol{B}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{B}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}},\end{split}$

where $\epsilon\sim\mathcal{N}(0,V_{\tau\tau}-\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau})$ is independent of
$\bm{B}\sim\mathcal{N}(0,{\boldsymbol{V}}_{xx})$. Since
${\boldsymbol{\Lambda}}_{\pi}={\boldsymbol{M}}{\boldsymbol{M}}^{T}$ is of full
rank, then ${\boldsymbol{M}}$ is reversible. Hence,

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1&\overset{\cdot}{\sim}\epsilon+{\boldsymbol{\beta}}^{T}({\boldsymbol{M}}^{-1})^{T}\bm{Q}\mid\bm{Q}^{T}\bm{Q}\leq\xi_{\alpha,{\boldsymbol{\lambda}}},\end{split}$

where
$\bm{Q}={\boldsymbol{M}}^{T}{\boldsymbol{B}}\sim\mathcal{N}(\bm{0},{\boldsymbol{M}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{M}})$,
and
${\boldsymbol{\beta}}^{T}({\boldsymbol{M}}^{-1})^{T}\bm{Q}\sim\mathcal{N}(0,\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau})$. Therefore,

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1&\overset{\cdot}{\sim}\sqrt{V_{\tau\tau}(1-R^{2})}\cdot\varepsilon_{0}+{\boldsymbol{\beta}}^{T}({\boldsymbol{M}}^{-1})^{T}\bm{Q}|\bm{Q}^{T}\bm{Q}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}.\end{split}$

By Eq. (15),

$r_{\phi_{\pi}}\triangleq\frac{{\mathbb{E}}\big{(}\bm{B}_{\infty}^{T}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\mid\bm{B}_{\infty}\in\mathcal{B}_{\phi_{\pi},\infty}\big{)}}{{\mathbb{E}}\big{(}\bm{B}^{T}_{\infty}{\boldsymbol{\beta}}_{\infty}{\boldsymbol{\beta}}_{\infty}^{T}\bm{B}_{\infty}\big{)}}=\lim_{N\rightarrow\infty}\frac{{\mathbb{E}}\big{(}\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\big{)}}{{\mathbb{E}}\big{(}\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}\big{)}},$

where $\bm{B}\sim N(\bm{0},\bm{V}_{xx})$ and

$\begin{split}\mathcal{B}_{\phi_{\pi}}=\left\\{{\boldsymbol{\mu}}:{\boldsymbol{\mu}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{\mu}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right\\},\end{split}$

where ${\boldsymbol{\lambda}}=(\lambda_{1},...,\lambda_{p})^{T}$ is the vector
of eigenvalues of
$\bm{P}=\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{V}}_{xx}^{1/2}$
and
$\bm{B}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{B}\sim\chi^{2}_{{\boldsymbol{\lambda}}}$.
For $\bm{Q}=\bm{M}^{T}\bm{B}$, we have

$\begin{split}{\mathbb{E}}\left(\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\right)&={\mathbb{E}}\left({\boldsymbol{\beta}}^{T}\bm{B}\bm{B}^{T}{\boldsymbol{\beta}}\mid\bm{B}^{T}{\boldsymbol{\Lambda}}_{\pi}\bm{B}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\\\
&={\boldsymbol{\beta}}^{T}\left(\bm{M}^{-1}\right)^{T}{\mathbb{E}}\left(\bm{Q}\bm{Q}^{T}\mid\bm{Q}^{T}\bm{Q}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\bm{M}^{-1}{\boldsymbol{\beta}},\end{split}$

where $\bm{Q}^{T}\bm{Q}\sim\chi^{2}_{{\boldsymbol{\lambda}}}$. Since
$\bm{Q}\sim N\left(\bm{0},\bm{M}^{T}\bm{V}_{xx}\bm{M}\right)$ and
$\bm{M}^{T}\bm{V}_{xx}\bm{M}$ has the same eigenvalues as $\bm{P}$, it can be
expressed as

$\begin{split}\bm{M}^{T}\bm{V}_{xx}\bm{M}=\bm{\Gamma}^{T}diag\\{\lambda_{1},...,\lambda_{p}\\}\bm{\Gamma},\end{split}$

where $\bm{\Gamma}$ is an orthogonal matrix. Let
$\bm{\Lambda}=diag\\{\lambda_{1},...,\lambda_{p}\\}$, then
$\tilde{\bm{Q}}=\bm{\Gamma}\bm{Q}\sim N\left(\bm{0},\bm{\Lambda}\right)$ and

$\begin{split}{\mathbb{E}}\big{(}\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\big{)}={\boldsymbol{\beta}}^{T}\left(\bm{M}^{-1}\right)^{T}\bm{\Gamma}^{T}{\mathbb{E}}\big{(}\tilde{\bm{Q}}\tilde{\bm{Q}}^{T}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}\bm{\Gamma}\bm{M}^{-1}{\boldsymbol{\beta}}.\end{split}$

Suppose $\bm{\tilde{}}{\bm{Q}}=\left(q_{1},...,q_{p}\right)^{T}$, for
$i\not=j$, symmetry implies

$\begin{split}{\mathbb{E}}\big{(}q_{i}q_{j}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}&={\mathbb{E}}\left({\mathbb{E}}\big{(}q_{i}q_{j}\mid
q_{j},\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\\\
&={\mathbb{E}}\left(q_{j}{\mathbb{E}}\big{(}q_{i}\mid
q_{j},\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\\\
&={\mathbb{E}}\big{(}q_{j}\times
0\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}=0.\end{split}$

Let
$c_{i}={\mathbb{E}}\left(\lambda_{i}Z_{i}^{2}\mid\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)$,
where $i=1,\ldots,p$ and $Z_{1},\ldots,Z_{p}$ are i.i.d. standard normal
distributed random variables, then

$\begin{split}{\mathbb{E}}\big{(}\tilde{\bm{Q}}\tilde{\bm{Q}}^{T}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}&=diag\left\\{{\mathbb{E}}\big{(}q_{i}^{2}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\big{\\}}_{i=1}^{p}\\\
&=diag\left\\{{\mathbb{E}}\left(\lambda_{i}Z_{i}^{2}\mid\sum_{j=1}^{p}\lambda_{j}Z_{j}^{2}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\right)\right\\}_{i=1}^{p}\\\
&=\bm{C}=diag\\{c_{1},...,c_{p}\\}.\end{split}$

Therefore,

$\begin{split}{\boldsymbol{\beta}}^{T}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}{\mathbb{E}}\big{(}\tilde{\bm{Q}}\tilde{\bm{Q}}^{T}\mid\tilde{\bm{Q}}^{T}\tilde{\bm{Q}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\big{)}\bm{\Gamma}\bm{M}^{-1}{\boldsymbol{\beta}}={\boldsymbol{\beta}}^{T}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}\bm{C}\bm{\Gamma}\bm{M}^{-1}{\boldsymbol{\beta}},\end{split}$

and thus,

$\begin{split}r_{\phi_{\pi}}&=\lim_{N\rightarrow\infty}\frac{{\mathbb{E}}\big{(}\bm{B}^{T}{\boldsymbol{\beta}}{\boldsymbol{\beta}}^{T}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\big{)}}{{\boldsymbol{\beta}}^{T}\bm{V}_{xx}{\boldsymbol{\beta}}}=\lim_{N\rightarrow\infty}\frac{{\boldsymbol{\beta}}^{T}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}\bm{C}\bm{\Gamma}\bm{M}^{-1}{\boldsymbol{\beta}}}{{\boldsymbol{\beta}}^{T}\bm{V}_{xx}{\boldsymbol{\beta}}}.\end{split}$

Since
${\boldsymbol{\beta}}^{T}\bm{V}_{xx}{\boldsymbol{\beta}}={\boldsymbol{\beta}}^{T}\big{(}\bm{M}^{-1}\big{)}^{T}\bm{\Gamma}^{T}\bm{\Lambda}\bm{\Gamma}\bm{M}^{-1}{\boldsymbol{\beta}}$,
we complete the proof. $\blacksquare$

###  A.12 Proof of Theorem 5

Under Condition 1,

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1&\overset{\cdot}{\sim}\epsilon+\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{B}\mid\bm{B}\in\mathcal{B}_{\phi_{\pi}}\\\
&\sim\epsilon+{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}\mid{\boldsymbol{B}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{B}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}},\end{split}$

where $\epsilon\sim\mathcal{N}(0,V_{\tau\tau}-\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau})$ is independent of
$\bm{B}\sim\mathcal{N}(0,{\boldsymbol{V}}_{xx})$. Since
${\boldsymbol{\Lambda}}_{\pi}={\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}$,
we have

$\begin{split}{\boldsymbol{\beta}}^{T}{\boldsymbol{B}}\mid{\boldsymbol{B}}^{T}{\boldsymbol{\Lambda}}_{\pi}{\boldsymbol{B}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}&\sim\epsilon_{\pi}+({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})^{-1}{\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{B}}\mid{\boldsymbol{B}}^{T}{\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{B}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}},\end{split}$

where $\epsilon_{\pi}\sim\mathcal{N}(0,\bm{V}_{\tau
x}\bm{V}_{xx}^{-1}\bm{V}_{x\tau}-({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})^{-1}({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}))$
is independent of ${\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{B}}$, and
$\xi_{\alpha,{\boldsymbol{\lambda}}}={\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi}\cdot\xi_{\alpha,1}$
is the $\alpha$-quantile of
${\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi}\cdot\chi^{2}_{1}$
distribution. Denote

$\begin{split}R^{2}_{\pi}=\frac{({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})^{-1}({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})}{V_{\tau\tau}},\end{split}$

we have

$\begin{split}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1&\overset{\cdot}{\sim}\epsilon+\epsilon_{\pi}+({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})^{-1}{\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{B}}\mid{\boldsymbol{B}}^{T}{\boldsymbol{\mu}}_{\pi}{\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{B}}\leq\xi_{\alpha,{\boldsymbol{\lambda}}}\\\
&\sim\sqrt{V_{\tau\tau}(1-R^{2}_{\pi})}\cdot\varepsilon_{0}+\sqrt{V_{\tau\tau}\cdot
R^{2}_{\pi}}\cdot L_{1,\xi_{\alpha,1}}.\end{split}$

Since $V_{\tau\tau}\cdot
R^{2}={\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}$,

$\begin{split}R^{2}_{\pi}=\frac{({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})^{-1}({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})}{{\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}}}R^{2}=\frac{({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})^{2}}{({\boldsymbol{\mu}}_{\pi}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\mu}}_{\pi})({\boldsymbol{\beta}}^{T}{\boldsymbol{V}}_{xx}{\boldsymbol{\beta}})}R^{2}.\end{split}$

Let $\theta$ be the angle between
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\mu}}_{\pi}$ and
$\big{(}{\boldsymbol{V}}_{xx}^{1/2}\big{)}^{T}{\boldsymbol{\beta}}$, then we
have

$\begin{split}R^{2}_{\pi}=(\cos\theta)^{2}\cdot R^{2}.\end{split}$

Moreover, denote the limit of $\theta$ and $R^{2}_{\pi}$ as $\theta_{\infty}$
and $R^{2}_{\pi,\infty}=(\cos\theta_{\infty})^{2}\cdot R^{2}_{\infty}$,
respectively, then

$\begin{split}\frac{{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{\pi}=1\big{)}}{{\text{var}}_{a}\big{(}\sqrt{N}(\hat{\tau}-\tau)\big{)}}&=\frac{V_{\tau\tau,\infty}(1-R^{2}_{\pi,\infty})+V_{\tau\tau,\infty}\cdot
R^{2}_{\pi,\infty}\cdot
v_{\alpha,1}}{V_{\tau\tau,\infty}}=1-(1-v_{\alpha,1})R^{2}_{\pi,\infty}.\end{split}$

By definition in Eq. (7),

$\begin{split}r_{\phi_{\pi}}&=1-(1-v_{\alpha,1})R^{2}_{\pi,\infty}/R^{2}_{\infty}=1-(1-v_{\alpha,1})(\cos\theta_{\infty})^{2}.\end{split}$

$\blacksquare$

###  A.13 Proof of Theorem 6

Take BCRD-ReB as an example and we suppose its acceptance probability in
second stage be $\alpha_{1}$. Since
$\hat{\tau}=\varrho\cdot\hat{\tau}_{1}+(1-\varrho)\cdot\hat{\tau}_{2}$, we
have

$\begin{split}{\text{var}}\left(\hat{\tau}\mid\phi_{BCRD\text{-}ReB}=1\right)&=\varrho^{2}\cdot{\text{var}}\left(\hat{\tau}_{1}\right)+(1-\varrho)^{2}\cdot{\text{var}}\left(\hat{\tau}_{2}\mid\phi_{\pi_{N_{1}}=1}\right).\end{split}$

Hence, for $N=N_{1}+N_{2}$, we have

$\begin{split}&\indent{\text{var}}_{a}\left(\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{BCRD\text{-}ReB}=1\right)\\\
&=\lim_{N_{1},N_{2}\rightarrow\infty}\left[\frac{\varrho^{2}}{N_{1}/N}\cdot{\text{var}}\left(\sqrt{N_{1}}(\hat{\tau}_{1}-\tau_{1})\right)+\frac{(1-\varrho)^{2}}{N_{2}/N}\cdot{\text{var}}\left(\sqrt{N_{2}}(\hat{\tau}_{2}-\tau_{2})\mid\phi_{\pi_{N_{1}}=1}\right)\right]\\\
&=\lim_{N_{1},N_{2}\rightarrow\infty}\left[\varrho\cdot{\text{var}}\left(\sqrt{N_{1}}(\hat{\tau}_{1}-\tau_{1})\right)+(1-\varrho)\cdot{\text{var}}\left(\sqrt{N_{2}}(\hat{\tau}_{2}-\tau_{2})\mid\phi_{\pi_{N_{1}}=1}\right)\right].\end{split}$
(A.5)

Denote $\hat{\tau}^{(2)}$ as the treatment effect estimator when we conduct
complete randomization in the second stage. For fixed $N_{1}$, as
$N_{2}\rightarrow\infty$, we have

$\begin{split}\lim_{N_{2}\rightarrow\infty}{\text{var}}\left(\sqrt{N_{2}}(\hat{\tau}_{2}-\tau_{2})\mid\phi_{\pi_{N_{1}}=1}\right)&=\lim_{N_{2}\rightarrow\infty}\left[1-(1-v_{\alpha_{1},\pi_{N_{1}}})R^{2}_{N_{2}}\right]\cdot{\text{var}}\left(\sqrt{N_{2}}(\hat{\tau}^{(2)}-\tau_{2})\right)\\\
&=\left[1-(1-v_{\alpha_{1},\pi_{N_{1}}})R^{2}_{\infty}\right]\cdot\lim_{N\rightarrow\infty}{\text{var}}\left(\sqrt{N}(\hat{\tau}-\tau)\right).\end{split}$

Moreover, under Condition 1, $\pi_{N_{1}}$ converges to the point mass at
${\boldsymbol{\beta}}_{\infty}$ as $N_{1}\rightarrow\infty$, and hence,
$v_{\alpha_{1},\pi_{N_{1}}}\rightarrow v_{\alpha_{1},1}$. Therefore,

$\begin{split}\lim_{N_{1},N_{2}\rightarrow\infty}{\text{var}}\left(\sqrt{N_{2}}(\hat{\tau}_{2}-\tau_{2})\mid\phi_{\pi_{N_{1}}=1}\right)&=\left[1-(1-v_{\alpha_{1},1})R^{2}_{\infty}\right]\cdot\lim_{N\rightarrow\infty}{\text{var}}\left(\sqrt{N}(\hat{\tau}-\tau)\right).\end{split}$

Since

$\begin{split}\lim_{N_{1}\rightarrow\infty}{\text{var}}\left(\sqrt{N_{1}}(\hat{\tau}_{1}-\tau_{1})\right)=\lim_{N\rightarrow\infty}{\text{var}}\left(\sqrt{N}(\hat{\tau}-\tau)\right),\end{split}$

for $\varrho\rightarrow\varrho^{*}$, we have

$\begin{split}\frac{{\text{var}}_{a}\left(\sqrt{N}(\hat{\tau}-\tau)\mid\phi_{BCRD\text{-}ReB}=1\right)}{{\text{var}}_{a}\left(\sqrt{N}(\hat{\tau}-\tau)\right)}&=\varrho^{*}+(1-\varrho^{*})\cdot\left[1-(1-v_{\alpha_{1},1})R^{2}_{\infty}\right].\end{split}$

Therefore, as $N_{1},N_{2}\rightarrow\infty$, we have

$\begin{split}\text{PRIASV}_{BCRD-
ReB}&=100\times\left\\{1-\left\\{\varrho^{*}+(1-\varrho^{*})\cdot\left[1-(1-v_{\alpha_{1},1})R^{2}_{\infty}\right]\right\\}\right\\}\\\
&=100\times(1-\varrho^{*})(1-v_{\alpha_{1},1})R^{2}_{\infty}.\end{split}$

PRIASV of ReM-ReB can be similarly derived. $\blacksquare$

##  Appendix B Additional simulation results

###  B.1 Comparison of ReO and ReM

To verify Theorem 2 and compare ReO and ReM, Figure B.1 shows average PRIV for
both methods across various covariate dimensions ($p$) with different
$\bar{R}^{2}$ values at $\rho=0$, $\sigma^{2}_{\boldsymbol{\beta}}=0.1$ and
$N=600$. As $\bar{R}^{2}$ increases, both methods show substantial PRIV
improvements. ReM’s average PRIV decreases with higher $p$, while ReO remains
stable at fixed $\bar{R}^{2}$.

Figure B.1: ReO versus ReM in terms of average PRIV for different $p$ and
$\bar{R}^{2}$ when $N=600$ and $\sigma^{2}_{\boldsymbol{\beta}}=0.1$.

###  B.2 Validation of acceptance probability

To verify the realized acceptance probability of the three re-randomization
schemes, running time of re-randomizing for each method in each setting is
recorded. Figure B.2 shows the reciprocal of the running time, which is
equivalent to acceptance rate for each method. We can see that as sample size
increases, the acceptance probability will get closer to 0.05, which is
exactly what we expect.

Figure B.2: The realized acceptance probability of ReO, ReB and ReM when (A)
$\sigma_{{\boldsymbol{\beta}}}^{2}=0.1$ and (B)
$\sigma_{{\boldsymbol{\beta}}}^{2}=1$. The black dotted line represents value
0.05.

###  B.3 Tuned acceptance rates for different re-randomization procedures and
their corresponding running time under linear settings

Table B.1 and Table B.2 show the numerically tuned acceptance rates of
different re-randomization procedures and their corresponding average running
time of sampling 1000 assignments under linear settings. We can see that the
computational cost of these methods are comparative under their acceptance
rates and thus, the comparisons of PRIV are fair.

Table B.1: Acceptance rates of different re-randomization procedures under linear settings. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.086 | 0.089 | 0.096 | 0.073 | 0.071 | 0.073 | 0.061 | 0.063 | 0.066
ReM-ReO | 0.082 | 0.086 | 0.087 | 0.069 | 0.071 | 0.071 | 0.063 | 0.062 | 0.062
BCRD-ReB | 0.052 | 0.047 | 0.045 | 0.042 | 0.040 | 0.045 | 0.046 | 0.038 | 0.034
BCRD-ReO | 0.053 | 0.048 | 0.049 | 0.046 | 0.040 | 0.040 | 0.043 | 0.040 | 0.035
| $p=5$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.093 | 0.087 | 0.092 | 0.069 | 0.069 | 0.071 | 0.063 | 0.068 | 0.058
ReM-ReO | 0.089 | 0.076 | 0.086 | 0.069 | 0.071 | 0.068 | 0.063 | 0.063 | 0.068
BCRD-ReB | 0.049 | 0.047 | 0.044 | 0.048 | 0.045 | 0.040 | 0.044 | 0.045 | 0.031
BCRD-ReO | 0.048 | 0.045 | 0.044 | 0.044 | 0.041 | 0.040 | 0.041 | 0.039 | 0.033
| $p=10$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.094 | 0.087 | 0.089 | 0.068 | 0.074 | 0.069 | 0.059 | 0.061 | 0.058
ReM-ReO | 0.092 | 0.085 | 0.082 | 0.067 | 0.071 | 0.067 | 0.059 | 0.060 | 0.059
BCRD-ReB | 0.052 | 0.047 | 0.045 | 0.048 | 0.079 | 0.040 | 0.041 | 0.039 | 0.039
BCRD-ReO | 0.047 | 0.044 | 0.039 | 0.043 | 0.041 | 0.038 | 0.041 | 0.038 | 0.036
| $p=20$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | $--$ | 0.090 | 0.087 | 0.063 | 0.067 | 0.066 | 0.058 | 0.056 | 0.059
ReM-ReO | $--$ | 0.088 | 0.081 | 0.065 | 0.063 | 0.068 | 0.058 | 0.056 | 0.058
BCRD-ReB | $--$ | 0.044 | 0.043 | 0.042 | 0.040 | 0.036 | 0.044 | 0.038 | 0.034
BCRD-ReO | $--$ | 0.039 | 0.037 | 0.042 | 0.038 | 0.035 | 0.041 | 0.040 | 0.033
Table B.2: Average running time (seconds) of sampling 1000 assignments using different re-randomization procedures under linear settings with tuned acceptance rates. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 0.852 | 0.867 | 0.848 | 1.581 | 1.571 | 1.579 | 2.653 | 2.642 | 2.656
ReM-ReB | 0.894 | 0.865 | 0.799 | 1.525 | 1.591 | 1.527 | 2.714 | 2.594 | 2.557
ReM-ReO | 0.919 | 0.892 | 0.873 | 1.586 | 1.568 | 1.560 | 2.554 | 2.644 | 2.662
BCRD-ReB | 0.850 | 0.919 | 0.907 | 1.750 | 1.689 | 1.402 | 2.510 | 2.783 | 2.796
BCRD-ReO | 0.826 | 0.887 | 0.820 | 1.581 | 1.652 | 1.537 | 2.664 | 2.591 | 2.713
| $p=5$
ReM | 0.923 | 0.950 | 0.917 | 1.744 | 1.721 | 1.721 | 2.897 | 2.912 | 2.908
ReM-ReB | 0.911 | 0.976 | 0.898 | 1.742 | 1.742 | 1.696 | 2.785 | 2.640 | 3.138
ReM-ReO | 0.918 | 1.090 | 0.923 | 1.724 | 1.695 | 1.746 | 2.762 | 2.833 | 2.626
BCRD-ReB | 0.957 | 1.018 | 0.980 | 1.653 | 1.629 | 1.670 | 2.844 | 2.729 | 3.392
BCRD-ReO | 0.957 | 1.004 | 0.934 | 1.783 | 1.760 | 1.654 | 2.958 | 2.900 | 3.091
| $p=10$
ReM | 1.057 | 1.037 | 1.035 | 1.951 | 1.956 | 1.971 | 3.360 | 3.345 | 3.360
ReM-ReB | 1.088 | 1.063 | 1.022 | 2.002 | 1.796 | 1.977 | 3.424 | 3.337 | 3.533
ReM-ReO | 1.076 | 1.064 | 1.069 | 1.997 | 1.882 | 2.006 | 3.420 | 3.370 | 3.461
BCRD-ReB | 1.007 | 1.065 | 1.040 | 1.846 | 1.024 | 1.853 | 3.495 | 3.346 | 3.022
BCRD-ReO | 1.072 | 1.033 | 1.098 | 1.999 | 1.913 | 1.924 | 3.451 | 3.410 | 3.250
| $p=20$
ReM | 1.339 | 1.320 | 1.311 | 2.473 | 2.459 | 2.459 | 4.272 | 4.259 | 4.284
ReM-ReB | $--$ | 1.350 | 1.324 | 2.652 | 2.469 | 2.445 | 4.433 | 4.445 | 4.275
ReM-ReO | $--$ | 1.328 | 1.343 | 2.558 | 2.543 | 2.333 | 4.326 | 4.363 | 4.329
BCRD-ReB | $--$ | 1.295 | 1.313 | 2.611 | 2.429 | 2.581 | 4.132 | 4.249 | 4.315
BCRD-ReO | $--$ | 1.333 | 1.326 | 2.492 | 2.510 | 2.468 | 4.227 | 3.908 | 4.319

###  B.4 Performance of two-stage ReB with various sample size proportion in
the first stage

To further study the effect of the proportion of sample in first stage on the
performance of two-stage ReB procedures, we conduct an additional experiment
on ReM, ReO, BCRD-ReO, BCRD-ReB, ReM-ReO and ReM-ReB with $r$ ranging from
0.05 to 0.8 when $N=1200$, $p=20$, $\bar{R}^{2}=0.5$ and $\rho=0$. We set
acceptance rate of ReM as $\alpha=0.05$ and tune the proper acceptance
probabilities of other re-randomization procedures numerically as before. The
running time ratio of different re-randomization procedures over ReM is
plotted in Figure B.3. It is observed that the running time of BCRD-ReB is
relatively unstable than the other two-stage ReB procedures. A conjecture is
that the instability arises from the use of Monte Carlo simulation during the
second stage for threshold estimation. The running time of ReM-ReB, however,
remains relatively stable, possibly because it constrains the allocation in
the first stage using ReM. Overall, the average running time of ReM and two-
stage procedures is comparable.

Average PRIVs of these methods under different specifications of $r$ are
plotted in Figure B.4, from which we can see the following facts immediately.
First, all two-stage methods achieve the best performance when $r$ falls into
the region of $[0.1,0.3]$. Second, ReM-ReB exhibits the best robustness to
different specifications of $r$ out of all the schemes considered. Thus, as a
rule of thumb, we recommend ReM-ReB with about 20% of samples in the first
stage as the primary method for practical usage.

Figure B.3: Running time ratio of re-randomization procedures over ReM under
different proportion of sample in first stage $r$ when total sample size
$N=1200$ and covariate dimensionality $p=20$ in linear setting. Figure B.4:
Average PRIV of ReO, ReM, ReM-ReB, ReM-ReO, BCRD-ReB and BCRD-ReO, with
different proportion of sample in first stage $r$ when total sample size
$N=1200$ and covariate dimensionality $p=20$ in linear setting.

###  B.5 Simulation results under non-linear settings

All above simulations are conducted under linear models. To evaluate the
performance of two-stage ReB under non-linear models, we repeated the above
numerical experiment with the linear regression model for potential outcomes
replaced by the non-linear regression model below:

$\begin{split}Y_{i}(1)&=a+{\boldsymbol{\beta}}_{1}^{T}\exp\\{{\boldsymbol{X}}_{i}\\}+\epsilon_{1,i},\
\ i=1,\ldots,N,\\\
Y_{i}(0)&={\boldsymbol{\beta}}_{0}^{T}\exp\\{{\boldsymbol{X}}_{i}\\}+\epsilon_{0,i},\
\ i=1,\ldots,N,\end{split}$

with all the other simulation settings unchanged. Under non-linear settings,
the acceptance rate of ReM is set to be $0.05$ and the acceptance
probabilities of other methods are tuned numerically as in experiments under
linear settings. The detailed acceptance rates as well as their corresponding
running time of sampling 1000 assignments are shown in Table B.3 and Table
B.4. The observations are similar to those in Section B.4.

Table B.3: Acceptance rates of different re-randomization procedures under non-linear settings. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.085 | 0.088 | 0.088 | 0.071 | 0.073 | 0.073 | 0.063 | 0.062 | 0.066
ReM-ReO | 0.096 | 0.086 | 0.083 | 0.068 | 0.071 | 0.072 | 0.063 | 0.061 | 0.061
BCRD-ReB | 0.055 | 0.053 | 0.049 | 0.042 | 0.040 | 0.042 | 0.044 | 0.040 | 0.033
BCRD-ReO | 0.062 | 0.047 | 0.046 | 0.044 | 0.045 | 0.043 | 0.043 | 0.039 | 0.036
| $p=5$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.096 | 0.085 | 0.096 | 0.068 | 0.069 | 0.071 | 0.061 | 0.063 | 0.063
ReM-ReO | 0.085 | 0.086 | 0.082 | 0.069 | 0.069 | 0.071 | 0.059 | 0.061 | 0.062
BCRD-ReB | 0.053 | 0.044 | 0.045 | 0.048 | 0.045 | 0.038 | 0.043 | 0.040 | 0.032
BCRD-ReO | 0.045 | 0.045 | 0.043 | 0.047 | 0.041 | 0.039 | 0.041 | 0.040 | 0.034
| $p=10$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | 0.094 | 0.087 | 0.092 | 0.065 | 0.067 | 0.069 | 0.059 | 0.061 | 0.059
ReM-ReO | 0.094 | 0.085 | 0.087 | 0.066 | 0.072 | 0.066 | 0.059 | 0.060 | 0.058
BCRD-ReB | 0.052 | 0.043 | 0.040 | 0.043 | 0.041 | 0.039 | 0.045 | 0.039 | 0.033
BCRD-ReO | 0.048 | 0.043 | 0.040 | 0.043 | 0.045 | 0.037 | 0.042 | 0.037 | 0.033
| $p=20$
ReM | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050 | 0.050
ReM-ReB | $--$ | 0.092 | 0.089 | 0.068 | 0.068 | 0.065 | 0.058 | 0.056 | 0.059
ReM-ReO | $--$ | 0.088 | 0.083 | 0.066 | 0.065 | 0.063 | 0.058 | 0.056 | 0.058
BCRD-ReB | $--$ | 0.044 | 0.045 | 0.046 | 0.039 | 0.036 | 0.079 | 0.038 | 0.032
BCRD-ReO | $--$ | 0.040 | 0.045 | 0.045 | 0.038 | 0.035 | 0.041 | 0.036 | 0.033
Table B.4: Average running time (seconds) of sampling 1000 assignments using different re-randomization procedures under non-linear settings with tuned acceptance rates. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 0.850 | 0.850 | 0.847 | 1.591 | 1.573 | 1.614 | 2.658 | 2.643 | 2.651
ReM-ReB | 0.910 | 0.866 | 0.867 | 1.568 | 1.524 | 1.591 | 2.565 | 2.652 | 2.553
ReM-ReO | 0.790 | 0.886 | 0.900 | 1.616 | 1.551 | 1.592 | 2.551 | 2.673 | 2.695
BCRD-ReB | 0.814 | 0.813 | 0.836 | 1.781 | 1.677 | 1.551 | 2.610 | 2.683 | 2.853
BCRD-ReO | 0.713 | 0.899 | 0.869 | 1.638 | 1.483 | 1.506 | 2.665 | 2.655 | 2.677
| $p=5$
ReM | 0.923 | 0.917 | 0.932 | 1.737 | 1.716 | 1.716 | 2.922 | 2.901 | 2.916
ReM-ReB | 0.873 | 0.957 | 0.852 | 1.742 | 1.728 | 1.718 | 2.923 | 2.844 | 2.825
ReM-ReO | 0.971 | 0.927 | 0.981 | 1.715 | 1.716 | 1.703 | 3.007 | 2.935 | 2.874
BCRD-ReB | 0.873 | 1.004 | 0.972 | 1.633 | 1.660 | 1.759 | 2.859 | 2.885 | 3.367
BCRD-ReO | 1.009 | 0.968 | 0.992 | 1.672 | 1.734 | 1.703 | 2.978 | 2.799 | 3.007
| $p=10$
ReM | 1.052 | 1.057 | 1.046 | 1.966 | 1.959 | 1.983 | 3.334 | 3.337 | 3.333
ReM-ReB | 1.066 | 1.093 | 1.002 | 2.098 | 2.037 | 1.979 | 3.410 | 3.294 | 3.446
ReM-ReO | 1.032 | 1.087 | 1.020 | 2.040 | 1.847 | 2.049 | 3.383 | 3.368 | 3.475
BCRD-ReB | 1.002 | 1.168 | 1.192 | 2.044 | 1.959 | 1.943 | 3.116 | 3.333 | 3.579
BCRD-ReO | 1.045 | 1.104 | 1.094 | 2.035 | 1.784 | 2.005 | 3.344 | 3.460 | 3.465
| $p=20$
ReM | 1.340 | 1.329 | 1.341 | 2.458 | 2.459 | 2.469 | 4.295 | 4.287 | 4.286
ReM-ReB | $--$ | 1.335 | 1.310 | 2.437 | 2.407 | 2.508 | 4.459 | 4.479 | 4.320
ReM-ReO | $--$ | 1.345 | 1.342 | 2.482 | 2.483 | 2.530 | 4.362 | 4.427 | 4.326
BCRD-ReB | $--$ | 1.299 | 1.275 | 2.395 | 2.550 | 2.605 | 2.313 | 4.355 | 4.543
BCRD-ReO | $--$ | 1.291 | 1.128 | 2.310 | 2.476 | 2.478 | 4.313 | 4.425 | 4.267

Table B.5 shows the average PRIV of different re-randomization methods with
different $p,r$ and $N$. It implies that two-stage ReB, especially ReM-ReB,
still outperforms ReM in most cases, as long as sample size is not very small,
even under nonlinear model. Figure B.5 shows the effect of $p$ on these
methods with fixed $r=0.2$ and $p=10$, and the effect of $N$ when $r=0.2$ and
$N=1200$. As sample size increases or dimension $p$ gets larger, the advantage
of two-stage ReB schemes gets more obvious. Also, ReM-ReB has the best
performance among the four schemes.

Table B.5: Average PRIV achieved by different re-randomization procedures under non-linear model. | $N=200$ | $N=600$ | $N=1200$
---|---|---|---
Scheme | $r=0.2$ | 0.3 | 0.4 | 0.2 | 0.3 | 0.4 | 0.2 | 0.3 | 0.4
| $p=2$
ReM | 48.28 | 48.28 | 48.28 | 47.53 | 47.53 | 47.53 | 47.28 | 47.28 | 47.28
ReM-ReB | 48.24 | 47.98 | 47.98 | 48.28 | 47.75 | 48.10 | 48.23 | 47.96 | 48.06
ReM-ReO | 46.86 | 46.92 | 48.04 | 47.44 | 47.73 | 48.35 | 48.19 | 47.91 | 47.85
BCRD-ReB | 38.74 | 34.97 | 30.40 | 39.34 | 33.54 | 29.90 | 38.65 | 33.99 | 28.57
BCRD-ReO | 37.17 | 33.06 | 29.23 | 38.42 | 33.14 | 29.48 | 38.72 | 33.65 | 28.88
| $p=5$
ReM | 40.55 | 40.55 | 40.55 | 41.24 | 41.24 | 41.24 | 40.86 | 40.86 | 40.86
ReM-ReB | 42.97 | 42.89 | 42.30 | 44.62 | 44.82 | 44.11 | 45.29 | 45.20 | 44.57
ReM-ReO | 38.66 | 40.92 | 41.63 | 43.79 | 44.21 | 43.74 | 45.03 | 44.95 | 44.23
BCRD-ReB | 36.53 | 32.38 | 27.04 | 37.71 | 32.40 | 28.51 | 37.97 | 33.66 | 29.13
BCRD-ReO | 31.18 | 29.06 | 25.75 | 36.61 | 32.05 | 28.08 | 37.59 | 32.55 | 28.87
| $p=10$
ReM | 32.88 | 32.88 | 32.88 | 32.77 | 32.77 | 32.77 | 33.36 | 33.36 | 33.36
ReM-ReB | 33.69 | 35.28 | 36.56 | 41.15 | 41.01 | 39.20 | 42.72 | 42.26 | 41.36
ReM-ReO | 25.57 | 30.78 | 33.74 | 38.58 | 39.11 | 38.86 | 42.44 | 41.87 | 41.09
BCRD-ReB | 30.53 | 28.91 | 25.86 | 35.59 | 31.04 | 26.90 | 36.85 | 32.70 | 28.57
BCRD-ReO | 20.64 | 23.56 | 21.80 | 33.11 | 29.80 | 26.71 | 35.33 | 31.95 | 27.99
| $p=20$
ReM | 23.93 | 23.93 | 23.93 | 24.08 | 24.08 | 24.08 | 25.66 | 25.66 | 25.66
ReM-ReB | $--$ | 24.50 | 27.46 | 34.30 | 34.40 | 33.88 | 39.75 | 38.68 | 37.01
ReM-ReO | $--$ | 17.66 | 24.07 | 31.05 | 32.52 | 32.63 | 38.16 | 37.60 | 36.39
BCRD-ReB | $--$ | 19.66 | 20.33 | 30.20 | 27.76 | 24.50 | 34.23 | 30.99 | 27.27
BCRD-ReO | $--$ | 11.64 | 14.60 | 25.60 | 25.36 | 23.16 | 33.21 | 30.21 | 26.90
Figure B.5: Performance of BCRD-ReO, BCRD-ReB, ReM-ReO and ReM-ReB compared
with ReM with $r=0.2$ under non-linear model. Sub-figure (A) investigates into
factor $N$ with covariate dimensionality $p$ fixed at 10, while sub-figure (B)
investigates into factor $p$ with $N$ fixed at 1200.

The effect of the proportion of sample in first stage on the performance of
two-stage ReB procedures can also be studied by conducting experiment on ReM,
BCRD-ReO, BCRD-ReB, ReM-ReO and ReM-ReB with $r$ ranging from 0.05 to 0.8 when
$N=1200$, $p=20$, $\bar{R}^{2}=0.5$ and $\rho=0$. Similar to the experiments
in linear setting, we set the acceptance rate of ReM be $\alpha=0.05$ and
calculate the acceptance probabilities of other re-randomization procedures
numerically. The running time ratio of re-randomization procedures over ReM is
shown in Figure B.6. Average PRIVs of these methods under different
specifications of $r$ are plotted in Figure B.7 and it is easy to see that
ReM-ReB still exhibits the best robustness to the proportion of sample in
first stage and similar conclusions as in linear settings can be drawn.

Figure B.6: Running time ratio of re-randomization procedures over ReM under
different proportion of sample in first stage $r$ when total sample size
$N=1200$ and covariate dimensionality $p=20$ under non-linear model. Figure
B.7: Average PRIV of ReM, ReM-ReB, ReM-ReO, BCRD-ReB and BCRD-ReO with
different proportion of sample in first stage $r$ when total sample size
$N=1200$ and covariate dimensionality $p=20$ under non-linear model.

###  B.6 Running time ratio in real-data analysis

Figure B.8 shows the running time ratio of re-randomization procedures against
ReM under different proportion of sample in first stage $r$ in real
application, i.e., Section 6.

Figure B.8: Running time ratio of re-randomization procedures over ReM under
different proportion of sample in first stage $r$ in real-data analysis.