Rerandomization and Regression Adjustment
Xinran Li Peng Ding






Rerandomization and Regression Adjustment

    Xinran Li and Peng Ding
[Xinran Li, Department of Statistics, University of Illinois at Urbana-Champaign, IL 61820 (E-mail: xinranli@illinois.edu). Peng Ding, Department of Statistics, University of California, Berkeley, CA 94720 (E-mail: pengdingpku@berkeley.edu).
]

    
================================================================================================================================================================================================================================================================================


Randomization is a basis for the statistical inference of treatment effects without strong assumptions on the outcome-generating process. Appropriately using covariates further yields more precise estimators in randomized experiments. R. A. Fisher suggested blocking on discrete covariates in the design stage or conducting analysis of covariance (ANCOVA) in the analysis stage. We can embed blocking into a wider class of experimental design called rerandomization, and extend the classical ANCOVA to more general regression adjustment. Rerandomization trumps complete randomization in the design stage, and regression adjustment trumps the simple difference-in-means estimator in the analysis stage. It is then intuitive to use both rerandomization and regression adjustment. Under the randomization-inference framework, we establish a unified theory allowing the designer and analyzer to have access to different sets of covariates. We find that asymptotically (a) for any given estimator with or without regression adjustment, rerandomization never hurts either the sampling precision or the estimated precision, and (b) for any given design with or without rerandomization, our regression-adjusted estimator never hurts the estimated precision. Therefore, combining rerandomization and regression adjustment yields better coverage properties and thus improves statistical inference. To theoretically quantify these statements, we discuss optimal regression-adjusted estimators in terms of the sampling precision and the estimated precision, and then measure the additional gains of the designer and the analyzer. We finally suggest using rerandomization in the design and regression adjustment in the analysis followed by the Huber–White robust standard error.

Keywords: covariate balance; experimental design; potential outcome; randomization








§ INTRODUCTION



In his seminal book Design of Experiments, <cit.> first formally discussed the value of randomization in experiments: randomization balances observed and unobserved covariates on average, and serves as a basis for statistical inference. Since then, randomized experiments have been widely used in agricultural sciences <cit.>, industry <cit.>, and clinical trials <cit.>. Recent years have witnessed the popularity of using randomized experiments in social sciences <cit.> and technology companies <cit.>. Those modern applications often have richer covariates. 



In completely randomized experiments, covariate imbalance often occurs by chance. <cit.> proposed to use analysis of covariance (ANCOVA) to adjust for covariate imbalance and thus improve estimation efficiency. <cit.>'s ANCOVA uses the coefficient of the treatment in the ordinary least squares (OLS) fit of the outcome on the treatment and covariates.  However, <cit.> criticized ANCOVA by showing that it can be even less efficient than the simple difference-in-means estimator under <cit.>'s potential outcomes framework.  <cit.>'s analyses allowed for treatment effect heterogeneity, in contrast to the existing literature on ANCOVA which often assumed additive treatment effects <cit.>. <cit.> proposed a solution to Freedman's critique by running the OLS of the outcome on the treatment, covariates, and their interactions. <cit.> showed the “optimality” of <cit.>'s estimator  within a class of regression-adjusted estimators. 



Aware of the covariate imbalance issue, <cit.> also proposed a strategy to actively avoid it in experiments. With a few discrete covariates, he proposed to use blocking, that is, to conduct completely randomized experiments within blocks of covariates. This remains a powerful tool in modern experiments <cit.>. Blocking is a special case of rerandomization <cit.>, which rejects “bad” random allocations that violate certain covariate balance criterion. Rerandomization can also deal with more general covariates. <cit.> demonstrated that rerandomization improves covariate balance. <cit.> further derived the asymptotic distribution of the difference-in-means estimator, and demonstrated that rerandomization improves its precision compared to complete randomization. 



Rerandomization and regression adjustment are two ways to use covariates to improve efficiency. The former uses covariates in the design stage, and the latter uses covariates in the analysis stage. It is then natural to combine them in practice, i.e., to conduct rerandomization in the design and use regression adjustment in the analysis. Several theoretical challenges remain. First, how do we conduct statistical inference? We will derive the asymptotic distribution of the regression-adjusted estimator under rerandomization without assuming any outcome-generating model. Our theory is purely randomization-based, in which potential outcomes are fixed numbers and the only randomness comes from the treatment allocation. 


Second, what is the optimal regression adjustment under rerandomization? The optimality depends on the criterion. We will introduce two notions of optimality, one based on the sampling precision and the other based on the estimated precision. Because our general theory allows the designer and analyzer to have different sets of covariates, it is possible that the estimated precision differs from the sampling precision asymptotically, even under the case with additive treatment effects. We will show that asymptotically, rerandomization never hurts either the sampling precision or the estimated precision, and <cit.>'s regression adjustment never hurts the estimated precision. Therefore, combining rerandomization and regression adjustment improves the coverage properties of the associated confidence intervals. Based on these findings, we suggest using <cit.>'s estimator in general settings, and show that the Huber–White variance estimator is a convenient approximation to its variance under rerandomization. Importantly, our theory does not rely on the linear model assumption. 


Third, how do we quantify the gains from the designer and analyzer? In particular, if the analyzer uses an optimal regression adjustment, what is the additional gain of rerandomization compared to complete randomization? If the designer uses rerandomization, what is the additional gain of using an optimal regression adjustment compared to the simple difference-in-means? Our theory can quantitatively answer these questions. 


This paper proceeds as follows. Section <ref> introduces the framework and notation. Section <ref> derives the sampling distribution of the regression-adjusted estimator under rerandomization. Section <ref> discusses optimal regression adjustment 
 in terms of
the sampling precision. Section <ref> addresses estimation and inference issues. Section <ref> discusses optimal regression adjustment 
in terms of
the estimated precision. 
Section <ref> quantifies the gains from the analyzer and the designer in both the sampling precision and the estimated precision. 
Section <ref> unifies the discussion and gives practical suggestions. 
Section <ref> uses  examples to illustrate the theory. Section <ref> concludes, and the Supplementary Material contains all the technical details. 






§ FRAMEWORK AND NOTATION


Consider an experiment on n units, with n_1 of them assigned to the treatment and n_0 of them assigned to the control. Let r_1 = n_1/n and r_0 = n_0/n be the proportions of units receiving the treatment and control. We use potential outcomes to define treatment effects <cit.>. For unit i, let Y_i(1) and Y_i(0) be the potential outcomes under the  treatment and control, and τ_i = Y_i(1) - Y_i(0) be the individual treatment effect. For this finite population of n units, the average potential outcome under treatment arm z (z=0,1) is Y̅(z) = n^-1∑_i=1^nY_i(z), and the average treatment effect is τ = n^-1∑_i=1^n τ_i = Y̅(1)-Y̅(0).
Let Z_i be the treatment assignment for unit i (Z_i=1 for the treatment; Z_i=0 for the control), and Z=(Z_1,Z_2, …,Z_n)' be the treatment assignment vector.  The observed outcome for unit i is Y_i = Z_iY_i(1)+(1-Z_i)Y_i(0).

 




 §.§ Regression adjustment in the analysis

In a completely randomized experiment (CRE), the probability that Z takes a particular value z=(z_1,…,z_n) is 
nn_1^-1, 
where ∑_i=1^n z_i = n_1 and ∑_i=1^n (1-z_i) = n_0 are fixed and do not depend on the values of covariates or potential outcomes. Equivalently, Z is a random permutation of a vector of n_1 1's and n_0 0's.  Let w_i=(w_i1, …, w_iJ)' be the J observed pretreatment covariates available to the analyzer. For descriptive convenience, we center these covariates at mean zero, i.e., n^-1∑_i=1^nw_i=0. Let 

    τ̂ = n_1^-1∑_i=1^n Z_i Y_i - n_0^-1∑_i=1^n (1-Z_i) Y_i,   τ̂_w = n_1^-1∑_i=1^n Z_i w_i- n_0^-1∑_i=1^n (1-Z_i) w_i

be the difference-in-means of the outcome Y and covariates w, respectively. Without covariate adjustment, τ̂ is unbiased for τ. After the experiment, the analyzer can improve the estimation precision for the average treatment effect by adjusting for the observed covariate imbalance τ̂_w. A general linear regression-adjusted estimator has the following equivalent forms: 

    τ̂(β_1, β_0)     =
    n_1^-1∑_i=1^n Z_i(
    Y_i - β_1' w_i) - 
    n_0^-1∑_i=1^n (1-Z_i)(
    Y_i - β_0' w_i)
    
        
    = τ̂ - ( r_0 β_1 + r_1 β_0)'τ̂_w 
    =  τ̂ - 'τ̂_w,

where β_1, β_0 and =  r_0 β_1 + r_1 β_0 are J dimensional coefficients. From (<ref>), τ̂(β_1, β_0) depends on (β_1, β_0) only through =  r_0 β_1 + r_1 β_0. Therefore, the choice of (β_1, β_0) is not unique to achieve the same efficiency gain. For simplicity, we will also call (<ref>) an adjusted estimator from now on. 


<cit.>'s ANCOVA chose β_1 =  β_0 to be the coefficient of w in the OLS fit of the observed outcome Y on the treatment Z and covariates w. <cit.> criticized ANCOVA because (a) the resulting estimator can be even less efficient than τ̂ and (b) the standard error based on the OLS can be inconsistent under the potential outcomes framework. <cit.> fixed (a) by choosing β_1 and β_0 to be the coefficients of w in the OLS fit of Y on w for treated and control units, respectively. The resulting adjusted estimator is numerically identical to the coefficient of Z in the OLS fit of Y on Z, w and Z×w. <cit.> fixed (b) by using the Huber–White robust standard error for linear models. Asymptotically, <cit.>'s estimator has smaller standard error and estimated standard error than τ̂. 



As a side note, <cit.>'s estimator also appeared in the semiparametric efficiency theory for the average treatment effect under independent sampling from a superpopulation <cit.>.  





 §.§ Rerandomization in the design


The above regression adjustment uses covariates in the analysis stage. We can also use covariates in the design stage to improve the quality of randomization and the efficiency of estimates. Before conducting the experiment, the designer collects K covariates x_i=(x_i1, …, x_iK)' for unit i. Similarly, we center the covariates at mean zero, i.e., n^-1∑_i=1^nx_i=0. Note that we allow x to be different from w. The CRE balances covariates on average, but an unlucky draw of the treatment vector can result in large covariate imbalance <cit.>. Therefore, it is sensible for the designer to check the covariate balance before conducting the experiment. 
Let 
τ̂_x  =n_1^-1∑_i=1^n Z_i x_i- n_0^-1∑_i=1^n (1-Z_i) x_i
be the difference-in-means of the covariates x between the treatment and control groups. It has mean zero under the CRE. However, imbalance in covariate distributions often occurs for a realized treatment allocation. We can discard those unlucky treatment allocations with large covariate imbalance, and rerandomize until the allocation satisfies a certain covariate balance criterion. This is rerandomization, which has the following steps: 



  (S1) collect covariate data and specify a covariate balance criterion;


  (S2) randomize the units into treatment and control groups;

  (S3)
if the allocation satisfies the balance criterion, proceed to (S4); otherwise, return to (S2);

  (S4) conduct the experiment using the accepted allocation from (S3). 
  

The balance criterion in (S1) can be a general function of the treatment assignment Z and the covariates (x_1, …, x_n).  

<cit.> focused on rerandomization using the Mahalanobis distance (ReM). ReM accepts a randomization if and only if M≡τ̂_x'{(τ̂_x)}^-1τ̂_x≤ a, where M is the Mahalanobis distance between the covariate means in two groups and a>0 is the predetermined threshold. <cit.> derived the asymptotic distribution of τ̂ under ReM, and showed that it is more precise than τ̂ under the CRE. They further showed that when a is small and x=w, the asymptotic variance of τ̂ under ReM is nearly identical to <cit.>'s adjusted estimator under the CRE. Therefore, rerandomization and regression adjustment both use covariates to improve efficiency of treatment effect estimates, but in the design and analysis stages, respectively. 
 







§ SAMPLING DISTRIBUTIONS OF REGRESSION ADJUSTMENT UNDER REM


Section <ref> shows that ReM trumps the CRE in the design stage and regression adjustment trumps the difference-in-means in the analysis stage. Therefore, it is natural to combine ReM and regression adjustment. Then a key question is how do we conduct statistical inference. This requires us to study the sampling distribution of the adjusted estimator τ̂(β_1, β_0) in (<ref>) under ReM. 
 




 §.§ Basics of randomization-based inference


To facilitate the discussion, we introduce some basic results from finite population causal inference. The first part describes fixed finite population quantities without randomness. The second part describes the repeated sampling properties and asymptotics under the CRE. 



  §.§.§ Finite population quantities, projections, and regularity conditions


For the treatment arm z (z=0,1), let S_Y(z)^2 = (n-1)^-1∑_i=1^n{Y_i(z)-Y̅(z)}^2 be the finite population variance of the potential outcomes, and S_Y(z),x=S_x,Y(z)' = (n-1)^-1∑_i=1^n{Y_i(z) - Y̅(z)}x_i' be the finite population covariance between the potential outcomes and covariates. Let S_x^2=(n-1)^-1∑_i=1^nx_i x_i' be the finite population covariance of the  covariates. We can similarly define S_Y(z), w, S^2_w, and other covariances. 

We introduce linear projections among these fixed quantities. For example, the linear projection of the potential outcome Y(z) on covariates w is Y̅(z) + β̃_z'w_i for unit i, with the coefficient

    β̃_z = _b∈ℝ^Jn^-1∑_i=1^n { Y_i(z) - Y̅(z) - b'w_i }^2 =  (S_w^2)^-1S_w, Y(z),    (z=0,1) .

The residual from this projection is Y_i(z) - Y̅(z) - β̃_z'w_i for unit i. Let S^2_Y(z)|w≡S_Y(z),w (S_w^2)^-1S_w, Y(z) denote the finite population variance of the linear projections, and S^2_Y(z)∖w≡ S^2_Y(z) - S^2_Y(z)|w the finite population variance of the residuals. We can similarly define S^2_τ|w, S^2_τ∖w,  S^2_Y(z)|x,  S^2_Y(z)∖x,S^2_τ|x, S^2_w|x and S^2_w∖x. 



 
The exact distributions of the estimators depend on unknown potential outcomes in general. We will use asymptotic approximations. Finite population asymptotics embeds the n units into a sequence of finite populations with increasing sizes. Technically, all quantities above depend on n, but we keep their dependence on n implicit for notational simplicity. Moreover, the sequence of finite populations must satisfy some regularity conditions to ensure the existence of the limiting distributions of the estimators. We use the regularity conditions motivated by the finite population central limit theorems <cit.>. 




 

As n→∞, 
the sequence of finite populations satisfies that, for z=0,1, 


  (i) r_z = n_z/n, the proportion of units receiving treatment z, has a positive limit; 

  (ii) 
the finite population variances and covariances, 
S_Y(z)^2, S^2_τ, S_x^2, S_w^2, S_Y(z),x, S_Y(z),w and S_x,w, have limiting values, 
and the limits of S_x^2 and S_w^2 are nonsingular;

  (iii) max_1≤ i≤ n|Y_i(z)-Y̅(z)|^2/n→ 0, max_1≤ i≤ nx_i_2^2/n → 0, and 
max_1≤ i≤ nw_i_2^2/n → 0.



In Condition <ref>, (i) and (ii) are natural, and (iii) holds almost surely if all the variables are independent and identically distributed (i.i.d.) draws from a superpopulation with more than two moments <cit.>. 
Throughout the paper, we assume the numbers of covariates K in the design and J in the analysis are both fixed and do not increase with the sample size n.

 



  §.§.§ Repeated sampling inference under the CRE





Under the CRE, over all nn_1 randomizations, n^1/2(τ̂-τ, τ̂_x',  τ̂_w') ' has mean 0 and covariance 

    V   ≡[ r_1^-1S_Y(1)^2 + r_0^-1S_Y(0)^2 - S_τ^2           r_1^-1S_Y(1),x+r_0^-1S_Y(0),x           r_1^-1S_Y(1),w+r_0^-1S_Y(0),w;          r_1^-1S_x, Y(1)+r_0^-1S_x,Y(0)                        (r_1r_0)^-1S_x^2                       (r_1r_0)^-1S_x, w;          r_1^-1S_w, Y(1)+r_0^-1S_w,Y(0)                       (r_1r_0)^-1S_w, x                        (r_1r_0)^-1S_w^2 ]
       ≡[ V_ττ V_τx V_τw; V_xτ V_xx V_xw; V_wτ V_wx V_ww ].


The finite population central limit theorem of <cit.> ensures that n^1/2(τ̂-τ, τ̂_x',  τ̂_w')' is asymptotically Gaussian with mean 0 and covariance matrix V under the CRE and Condition <ref>. We use  for two sequences of random vectors (or distributions) converging weakly to the same distribution. Therefore, 
n^1/2(τ̂ -τ , τ̂_x', τ̂_w')'
𝒩(0, V).



We define linear projections for random variables. We use 𝔼(·), (·) and (·) for mean, variance and covariance, and (·|·) and (·|·) for linear projections and corresponding residuals, exclusively under the CRE.  For example, the linear projection of τ̂ on τ̂_w is 
(τ̂|τ̂_w) = τ + 'τ̂_w,
with the coefficient 

    = _b∈ℝ^J𝔼(
    τ̂ - τ - b'τ̂_w)^2 
    = {(τ̂_w)}^-1(τ̂_w, τ̂)
     = V_ww^-1V_wτ .

The residual from this projection is 
(τ̂|τ̂_w) = τ̂ - (τ̂|τ̂_w) = τ̂ - τ - 'τ̂_w.
We can similarly define 
(τ̂_x|τ̂_w) and (τ̂_x|τ̂_w). 



Finally, the three linear projection coefficients β̃_1, β̃_0 and  defined in (<ref>) and (<ref>) have the following relationship.



r_0β̃_1+r_1 β̃_0 = .


Proposition <ref> is related to the non-uniqueness of the regression adjustment in (<ref>). It is important for the discussion below. 



 §.§ Asymptotic distribution of regression adjustment under ReM



Equipped with the tools in Section <ref>, we now can derive the asymptotic distribution of τ̂(β_1, β_0) under ReM. We first fix the coefficients β_1 and β_0, and will devote several sections to discuss the optimal choices of them. 


For unit i, let Y_i(z;β_z) ≡   Y_i(z) - β_z'w_i be the “adjusted" potential outcome under the treatment level z (z=0,1), 
and τ_i(β_1, β_0)≡τ_i- (β_1-β_0)'w_i be the “adjusted" individual treatment effect. The average “adjusted” treatment effect τ(β_1, β_0) ≡ n^-1∑_i=1^nτ_i(β_1, β_0) =  τ is identical to the average unadjusted  treatment effect because of the centering of w̅=0. The “adjusted" observed outcome is   Y_i(β_1,β_0)=Z_iY_i(1;β_1)+(1-Z_i)Y_i(0;β_0). The adjusted estimator (<ref>) is essentially the difference-in-means estimator with the “adjusted” potential outcomes. 
For z=0,1, let
S_Y(z;β_z)^2 and 
S_Y(z;β_z)|x^2 
be the finite population variances of 
Y_i(z;β_z) and its linear projection on x_i. Let
S_τ(β_1, β_0)^2 
and 
S_τ(β_1, β_0)|x^2
be the finite population variances of τ_i(β_1, β_0) and its linear projection on x_i. 
From Section <ref>, under the CRE, 
the variance of n^1/2{τ̂(β_1, β_0) - τ} is 

    V_ττ(β_1, β_0)
    = r_1^-1S_Y(1;β_1)^2 + r_0^-1S_Y(0;β_0)^2 - S_τ(β_1, β_0)^2,

and the squared multiple correlation between 
τ̂(β_1, β_0) and τ̂_x is  
<cit.>

    R^2_τ, x(β_1, β_0)
    = 
    {(τ̂(β_1, β_0) |τ̂_x)}/{τ̂(β_1, β_0) }
    =
    r_1^-1S_Y(1;β_1)|x^2+r_0^-1S_Y(0;β_0)|x^2- S_τ(β_1, β_0)|x^2/r_1^-1S_Y(1;β_1)^2 + r_0^-1S_Y(0;β_0)^2 - S_τ(β_1, β_0)^2
    .



The asymptotic distribution of τ̂(β_1, β_0) under ReM is a linear combination of two independent random variables ε and L_K,a, where ε∼𝒩(0,1) is  a standard Gaussian random variable and L_K,a∼ D_1|D'D≤ a is a truncated Gaussian random variable with D = (D_1,…, D_K)∼𝒩(0, I_K). Let ℳ denote the event M≤ a.




Under ReM and Condition <ref>, 

    n^1/2{τ̂(β_1, β_0)-τ}|ℳ  
     V_ττ^1/2(β_1, β_0)
    [
    {1-R^2_τ, x(β_1, β_0)}^1/2·ε + {R^2_τ, x(β_1, β_0)}^1/2·  L_K,a].



The ε component in (<ref>) represents the part of τ̂(β_1, β_0) that cannot be explained by τ̂_x and is thus unaffected by rerandomization. The L_K,a component in (<ref>) represents the part of τ̂(β_1, β_0) that can be explained by τ̂_x and is thus affected by rerandomization. Moreover, 
the asymptotic distribution (<ref>) is symmetric around zero, 
and the adjusted estimator is consistent for the average treatment effect, for any fixed values of the coefficients β_1 and β_0.  Theorem <ref> immediately implies the following two important special cases. 





  §.§.§ Special case: regression adjustment under the CRE



The CRE is a special case of ReM with a=∞. Therefore, Theorem <ref> implies that under the CRE,  n^1/2{τ̂(β_1, β_0)-τ} is asymptotically Gaussian with mean 0 and variance V_ττ(β_1, β_0). 


 
Under CRE and Condition <ref>, 
n^1/2{τ̂(β_1, β_0)-τ} 
V_ττ^1/2(β_1, β_0) ·ε .




Corollary <ref> is a known result from <cit.> and <cit.>. 



  §.§.§ Special case: no covariate adjustment under ReM

Using Theorem <ref> with β_1=β_0=0, we can immediately obtain the asymptotic distribution of τ̂≡τ̂(0, 0) under ReM. Let R^2_τ, x≡ R^2_τ, x(0, 0) be the squared multiple correlation between τ̂ and τ̂_x under the CRE: 

    R^2_τ, x    = 
    {(τ̂|τ̂_x)}/(τ̂)	
    
    = 
    V_τxV_xx^-1V_xτ/V_ττ
    = 
    r_1^-1S_Y(1)|x^2+r_0^-1S_Y(0)|x^2
    - S_τ|x^2
    /r_1^-1S_Y(1)^2 + r_0^-1S_Y(0)^2 - S_τ^2.

Then τ̂ has the following asymptotic distribution. 

Under ReM and Condition <ref>, 

    n^1/2( τ̂ - τ)|ℳ      
    V_ττ^1/2{( 1-R^2_τ, x)^1/2·ε + ( R^2_τ, x)^1/2· L_K,a}.



Corollary <ref> is a main result of <cit.>. 







§ 𝒮-OPTIMAL REGRESSION ADJUSTMENT


How do we choose the adjustment coefficients (β_1, β_0) or ? It is an important practical question. From Theorem <ref>, the adjusted estimator is consistent for any fixed coefficients β_1 and β_0. Therefore, it is intuitive to choose the coefficients to achieve better precision. A measure of precision is based on the quantile ranges of an estimator. 

We introduce the notion of the 𝒮-optimal adjusted estimator, using 𝒮 to emphasize its dependence on the sampling distribution. 



Given the design,
τ̂(β_1, β_0) is 𝒮-optimal if n^1/2{τ̂(β_1, β_0)-τ} has the shortest asymptotic 1-α quantile range among all adjusted estimators in (<ref>), for any α∈ (0,1). 



 
In general, quantile ranges are not unique. Importantly, the asymptotic distribution in (<ref>) is symmetric and unimodal around τ <cit.>. We consider only symmetric quantile ranges because they have the shortest lengths <cit.>.
The 𝒮-optimal adjusted estimator also has the smallest asymptotic variance among all estimators in (<ref>) <cit.>. Conversely, if all adjusted estimators in (<ref>) are asymptotically Gaussian, then the one with the smallest asymptotic variance is 𝒮-optimal. 
 


Theorem <ref> shows a complicated relationship between the coefficients (β_1, β_0) and the asymptotic distribution (<ref>). Below we simplify (<ref>). 
Let (τ̂_w|τ̂_x) ≡V_wxV_xx^-1τ̂_x be the linear projection of τ̂_w on τ̂_x, and 
(τ̂_w|τ̂_x)
≡τ̂_w - 

V_wxV_xx^-1τ̂_x be the residual from this linear projection. 
We further consider two projections. First, the linear projection of (τ̂|τ̂_x) on 
(τ̂_w|τ̂_x) has coefficient and squared multiple correlation

    _ 
    =  (
    V_wxV_xx^-1V_xw
    )^-1V_wxV_xx^-1V_xτ and  
    R^2_ .

Second, the linear projection of (τ̂|τ̂_x) on (τ̂_w|τ̂_x) has coefficient and squared multiple correlation

    _≡ (
    V_ww - V_wxV_xx^-1V_xw
    )^-1
    (V_wτ - V_wxV_xx^-1V_xτ )
     and  
    R^2_ .


Technically, the expressions for _ and _ above are well-defined only if the covariance matrices of  (τ̂_w|τ̂_x) and (τ̂_w|τ̂_x) are nonsingular. 
Otherwise, they are not unique. However, this will not cause any issues in the later discussions because the linear projections themselves



are always unique. 

Recall that we have defined V in (<ref>), and S^2_w|x and 
S^2_w∖x as the finite population covariances of the linear projections of w on x and the corresponding residuals. 
The following proposition shows the relationship among the three linear projection coefficients (, _, _). 


S^2_w∖x(-_) + 
S^2_w|x
(-_)  =  0.


The linear projection coefficients (,   _, _  ) are different in general.
However, if any two of them are equal, all of them must be equal with nonsingular S^2_w∖x and S^2_w|x. 
The following theorem decomposes the asymptotic distribution (<ref>) based on ( _, _). 



Under ReM and Condition <ref>, recalling that ≡ r_0 β_1 + r_1 β_0, we have

    n^1/2{τ̂(β_1, β_0)-τ}|ℳ
          {
    V_ττ(1-R^2_τ, x)
    (1-R^2_)
    + (r_1r_0)^-1(-_)'
    S^2_w∖x(-_)
    }^1/2·ε
           +  
    { V_ττR^2_τ, x(1 - R^2_) + (r_1r_0)^-1(-_)'
    S^2_w|x(-_)
    }^1/2·   L_K,a.






The asymptotic distribution (<ref>) has two independent components. The ε component in (<ref>) represents the part of τ̂(β_1, β_0) that is orthogonal to τ̂_x. The coefficient of ε attains its minimal value at  =_ with squared minimal value 
V_ττ   (1-R^2_τ, x)   (1-R^2_).
The first term V_ττ is the variance of n^1/2 (τ̂ - τ ). The second term 1-R^2_τ, x represents the proportion of the variance of τ̂ unexplained by τ̂_x. The third term 1-R^2_ represents the proportion of the variance of τ̂ unexplained by τ̂_w, after projecting onto the space orthogonal to τ̂_x. 
 

The L_K,a component in (<ref>) represents the linear projection of  τ̂(β_1, β_0) on τ̂_x with the ReM constraint. The coefficient of L_K,a attains its minimal value at  =_  
with squared minimal value 
V_ττ  R^2_τ, x  ( 1 - R^2_ ).
The first term V_ττ is again the variance of n^1/2(τ̂-τ). The second term R^2_τ, x represents the proportion of the variance of τ̂ explained by τ̂_x. The third term 1-R^2_ represents the proportion of the variance of τ̂ unexplained by τ̂_w, after projecting onto the space of τ̂_x.
 


Because _ and  _ are different in general, the coefficients of ε and L_K,a cannot attain their minimal values simultaneously. Consequently,
the adjusted estimator τ̂(β̃_1, β̃_0) may not be 𝒮-optimal under ReM, i.e., it may not have the shortest asymptotic 1-α  quantile range among (<ref>) for all α∈ (0,1). 
 
 

 
The 𝒮-optimal adjustment is complicated under ReM, especially when the designer and the analyzer have different covariate information. We will consider different scenarios based on the relative amount of covariate information used by the designer and the analyzer. 





 §.§ The analyzer has no less covariate information than the designer



We first consider the scenario under which the covariates w in the analysis can linearly represent the covariates x in the design. 


There exists a constant matrix B_1 ∈ℝ^K× J such that x_i=B_1 w_i for all unit i. 



Condition <ref> holds when the analyzer has access to all the covariates used in the design, and possibly collects more covariates after the experiment. For example, Condition <ref> holds if x is a subset of w. Under Condition <ref>, we can simplify the asymptotic distribution of the adjusted estimator under ReM. Analogous to (<ref>) and (<ref>), let R^2_τ, w be the squared multiple correlation between τ̂ and τ̂_w under the CRE <cit.>: 

    R^2_τ, w    = 
    {(τ̂|τ̂_w)}/(τ̂)	
    
    = 
    V_τwV_ww^-1V_wτ/V_ττ
    = 
    r_1^-1S_Y(1)|w^2+r_0^-1S_Y(0)|w^2
    - S_τ|w^2
    /r_1^-1S_Y(1)^2 + r_0^-1S_Y(0)^2 - S_τ^2.



Under Conditions <ref> and <ref>,

    _=_=,    
    R^2_ =   ( R^2_τ,w - R^2_τ,x ) / (1-R^2_τ,x) ,    
    R^2_=1,

and the asymptotic distribution of τ̂(β_1, β_0) under ReM is 

    n^1/2{τ̂(β_1, β_0)-τ}|ℳ
       {
    V_ττ(1-R^2_τ, w)
    + 
    (r_1r_0)^-1(-)'
    S^2_w∖x(-)
    }^1/2·ε
        +  
    {(r_1r_0)^-1(-)'
    S^2_w|x(-)
    }^1/2·   L_K,a.






From Corollary <ref>, the coefficients of ε and L_K,a attain minimum values at the same =  . We can then derive the 𝒮-optimal adjusted estimator and its asymptotic distribution. 


Under ReM and Conditions <ref> and <ref>, the 𝒮-optimal adjusted estimator is attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 , 
with the asymptotic distribution 

    n^1/2{τ̂(β̃_1, β̃_0)-τ}|ℳ      {V_ττ(1-R^2_τ,w)}^1/2·ε.




From Proposition <ref>, an optimal choice is (  β_1, β_0) = ( β̃_1, β̃_0 ) under Condition <ref>. An important feature of Theorem <ref> is that the limiting distribution in (<ref>) does not depend on covariates x and the threshold a of ReM.
Theorem <ref> has many implications, 
as discussed below. 



  §.§.§ Special case: 𝒮-optimal adjustment under the CRE



Theorem <ref> holds for the CRE (ReM with x = ∅ and a = ∞). 
It thus recovers the optimality property of <cit.>'s estimator under the CRE previously proved by <cit.>. 
Therefore, 
	 when the analyzer has no less covariate information than the designer, the
	  𝒮-optimal adjusted estimators under ReM and the CRE are the same and follow the same asymptotic distribution.
	 This implies that, with more covariates in the analysis, there is no additional gain from the designer through ReM as long as the analyzer performs the optimal  adjustment.
Section <ref> later contains related discussions.
	 


  §.§.§ Special case: the designer and analyzer have the same covariates


 Consider the case where the analyzer has the same covariates as the designer (x = w). Compare τ̂ under ReM to the 𝒮-optimal adjusted estimator τ̂(β̃_1, β̃_0) under the CRE. 
From Corollary <ref> and Section <ref>, the former has an additional independent component of (V_ττR^2_τ, x)^1/2 L_K,a in the asymptotic distribution. 
	When the threshold a is small, this additional component is approximately zero, and thus they have almost the same asymptotic distribution.
	Therefore, we can view rerandomization as covariate adjustment in the design stage <cit.>. 
	Moreover, the former has the following advantages. 
	First, rerandomization in the design stage does not use the outcome data. Second, τ̂ is simpler and thus provides a more transparent analysis <cit.>. Using τ̂ in rerandomization can thus avoid bias due to a specification search  of the outcome model (i.e., data snooping). 
	


Remark <ref> later contains related discussions.

 



 §.§ The analyzer has no more covariate information than the designer


We then consider the scenario under which the covariates x in the design can linearly represent
the covariates w in the analysis.


There exists a constant matrix B_2∈ℝ^J× K such that w_i=B_2x_i for all unit i. 



Condition <ref> is reasonable when the analyzer has access to only part of the covariates used in the design due to privacy or other reasons. 
For example, 
Condition <ref> holds if w is a subset of x. It also reflects the situation where the analyzer uses only the difference-in-means estimator with w  = ∅ 
even though 
the designer conducts ReM with x. Condition <ref> implies S^2_w∖x = 0, which further implies that the coefficient of ε in the asymptotic distribution (<ref>) does not depend on (β_1,β_0). We can then simplify the asymptotic distribution of the adjusted estimator. 
 


Under Conditions <ref> and <ref>, 

    _=,    R^2_ =  0,    
    	R^2_= R^2_τ, w/R^2_τ, x,
      S^2_w∖x=0,   S^2_w|x=S^2_w,

	and the asymptotic distribution of τ̂(β_1, β_0) under ReM is 

    n^1/2{τ̂(β_1, β_0)-τ}|ℳ
          {
    	V_ττ(1-R^2_τ, x)
    } ^1/2·ε   
    +  
    { V_ττ(R^2_τ, x - R^2_τ, w) + (r_1r_0 )^-1(-)'
    S^2_w(-)
    }^1/2·  L_K,a.




Under Condition <ref>, (τ̂_w|τ̂_x) = 0, and thus as discussed earlier, the projection coefficient _ is not unique. Nevertheless, Corollary <ref> does not depend on _.
Based on Corollary <ref>, we can obtain the 𝒮-optimal adjusted estimator and its asymptotic distribution under ReM. Let 

    ρ^2_τ, x∖w
    	= ( R^2_τ, x-R^2_τ, w)/(1 - R^2_τ, w)
    	∈ [0,1]

be the additional proportion of the variance of τ̂ explained by the covariates x in the design, after explained by the covariates w in the  analysis.


Under ReM and Conditions <ref> and <ref>, the 𝒮-optimal adjusted estimator is 
attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 ,  with the asymptotic distribution

    n^1/2{τ̂(β̃_1, β̃_0)-τ}|ℳ       
    V_ττ^1/2{( 1-R^2_τ, x)^1/2·ε + 
    ( R^2_τ, x-R^2_τ, w)^1/2·  L_K,a}
       ∼   { V_ττ(1-R^2_τ, w) }^1/2{( 1-ρ^2_τ, x∖w)^1/2·ε + 
    ( ρ^2_τ, x∖w)^1/2· L_K,a}.





From Theorem <ref>, although the analyzer has less covariate information than the designer of ReM, 
s/he can still obtain the 𝒮-optimal adjusted estimator using only the covariate information in the analysis. 


Theorems <ref> and <ref> give identical optimal coefficients, but different asymptotic distributions of the optimal estimators. 
When the designer and the analyzer have the same covariates (x = w), 
both Theorems <ref> and <ref> hold and give identical results. 
Specifically, ρ^2_τ, x∖w in (<ref>) reduces to zero, 
and the asymptotic distribution in (<ref>) simplifies to (<ref>), a Gaussian limiting distribution.


 

From Corollary <ref> and Theorem <ref>, under ReM, 
the asymptotic distribution of the 𝒮-optimal adjusted estimator  τ̂(β̃_1, β̃_0)
differs from that of τ̂ only in the coefficient of the truncated Gaussian random variable L_K,a. With a small threshold a, L_K,a is close to zero and thus the gain from adjustment is small. 
Similar to the discussion in Section <ref>, 
although τ̂ loses a little sampling precision compared to the optimal adjusted estimator, 
it does have the advantage of avoiding data snooping and improving transparency. 

 



 §.§ General scenarios


A practical complication is that the designer and analyzer may not communicate. Then it is possible that the designer and the analyzer do not use the same covariate information <cit.>. Consequently, the analyzer has part of the covariate information in the design and additional covariate information. Neither Condition <ref> or <ref> holds. Under general scenarios, unfortunately, the 𝒮-optimal adjusted estimator may not exist, in the sense that there does not exist an estimator among (<ref>) that has the shortest asymptotic 1-α quantile range for all α∈ (0,1). 



Some sub-optimal strategies exist. First, we can consider the adjusted estimator with the smallest asymptotic variance or the shortest asymptotic 1-α quantile range for a particular α∈ (0,1). The Supplementary Material gives the formulas for the former. However, explicit formulas  for the latter do not exist. 





Second, when a is small, L_K,a≈ 0, and the asymptotic distribution (<ref>) under ReM depends mainly on the ε component. The coefficient of ε attains its minimal value at _. 
Ignoring the L_K,a component, _ gives the 𝒮-optimal adjusted estimator. However, this result is not useful because it is infeasible for the analyzer to consistently estimate _ due to the incomplete information of the covariates  in the design. 
 


Third, we can still use τ̂(β̃_1, β̃_0) as a convenient adjusted estimator because we can easily obtain it via the OLS. Not surprisingly, this estimator is not 𝒮-optimal in general and it can be even worse than τ̂ under ReM. 
When a  is small, the ε components are the dominating terms in 
their asymptotic distributions under ReM. Therefore, we compare the coefficients of ε. 




Under ReM and Condition <ref>, 
the squared coefficient of ε is V_ττ(1-R^2_τ,w){1 - R^2_τ,x(β̃_1, β̃_0)} in the asymptotic distribution (<ref>) of τ̂(β̃_1, β̃_0), and is V_ττ(1 - R^2_τ,x) in the asymptotic distribution (<ref>) of τ̂. The former is smaller than or equal to the latter if and only if 

    R^2_τ,w+ (1-R^2_τ,w)· R^2_τ,x(β̃_1, β̃_0)≥  R^2_τ,x




A sufficient condition for (<ref>) is R^2_τ,w≥ R^2_τ,x, which holds under Condition <ref>. 
Recall that R^2_τ,x in (<ref>) measures the covariate information of the designer, and R^2_τ,w in (<ref>) measures the covariate information of the analyzer. From Theorem <ref>, when the analyzer has more covariate information, τ̂(β̃_1, β̃_0) is more precise than τ̂ if  the threshold a for ReM is small. 




A counterexample for (<ref>) is that the finite population partial covariance between Y(z) and w given x is zero for z=0,1. In this case,  the squared coefficient of ε 
for τ̂(β̃_1, β̃_0) is larger than or equal to that for τ̂. 
Intuitively, 

this is because 
the covariates in the analysis are unrelated to the potential outcomes after adjusting for the covariates in the design and using them only introduces additional variability. 
In the extreme case where x can linearly represent Y(1) and Y(0), 
the squared coefficient of ε for τ̂ is zero, 
while that for τ̂(β̃_1, β̃_0) is generally positive. See the Supplementary Material for more details.




 
Below we use a numerical example to illustrate the results above. It shows that   τ̂(β̃_1, β̃_0) can be superior or inferior to τ̂.


We choose n=1000, r_1=r_0=0.5, and generate the covariates and the potential outcomes using i.i.d. samples from the following model: 

    x,η, δi.i.d.∼𝒩(0,1),    
    w = x + η,    
    Y(0) = 2x + ρη + (1-ρ^2)^1/2δ,    
    Y(1) = Y(0) + 1.

Once generated, the covariates and potential outcomes are all fixed. 
We use ReM based on the covariate x, and choose the threshold a to be the 0.001th quantile of the χ^2_1 random variable. 
We then use regression adjustment based on the covariate w. 
Figure <ref>(a) shows the histograms of τ̂(β̃_1, β̃_0) and τ̂ under 
ReM when ρ=0.9.  In this case, regression adjustment increases the sampling precision. Figures <ref>(b) shows the histograms when ρ=0. In this case, regression adjustment decreases the sampling precision. The case with ρ = 0 reflects the scenario that the designer only gives a covariate with measurement error to the analyzer, possibly due to some privacy consideration. 

 











§ ESTIMATING THE SAMPLING DISTRIBUTIONS OF THE ESTIMATORS


Sections <ref> and <ref>  focused on the asymptotic distribution of the adjusted estimator and discussed optimal choices of the coefficients. In practice, we usually report the uncertainty of estimators in terms of confidence intervals when conducting frequentist inference. Confidence intervals are related to the quantile ranges of the estimated distributions of the corresponding estimators. Therefore, compared to 𝒮-optimality, a more practical definition of the optimal adjusted estimator should be based on the quantile ranges of the estimated distributions. This subtle issue does not exist in many other statistical inference problems, because usually consistent estimators exist for the true asymptotic distributions of the estimators. For example, in standard statistical problems, we can consistently estimate the variance of a Gaussian limiting distribution. Because of the possible miscommunication between the designer and analyzer, the analyzer may not be able to estimate all quantities based on the observed data in general. This is a feature of our framework. 



In this section, we discuss the estimation of the sampling distributions for fixed (β_1, β_0). 
In the next section, we will discuss the optimal choice of these coefficients. 
 
 


 §.§ The analyzer knows all the information in the design


We first consider the scenario under which the analyzer knows all the information of the designer.



The analyzer 
knows all the information in the design, 
including the covariates
x and the threshold a for ReM.

Condition <ref> implies Condition <ref>. 
However, Condition <ref> does not imply Condition <ref>, because the analyzer may not know which covariates are used in the design or which threshold a is chosen for ReM. 
From Theorem <ref>, 
the asymptotic distribution (<ref>) of τ̂(β_1, β_0) under ReM depends on 
V_ττ(β_1, β_0) and R^2_τ, x(β_1, β_0). Under treatment arm z (z=0,1), let
s^2_Y(z;β_z), s_Y(z;β_z),w
and s_Y(z;β_z),x
be the sample variance and covariances for the “adjusted" observed  outcome Y_i - β_z' w_i, covariates in the analysis w_i, and covariates in the design x_i; let s^2_Y(z;β_z)|x be the sample variance of the linear projection of Y_i - β_z' w_i on x_i. 
We estimate V_ττ(β_1, β_0) by 

    V̂_ττ(β_1, β_0)     =     
    r_1^-1s^2_Y(1;β_1) + r_0^-1s^2_Y(0;β_0)
          - 
    (s_Y(1;β_1),w-
    s_Y(0;β_0),w ) (S_w^2)^-1
    (s_w,Y(1;β_1)-s_w,Y(0;β_0) ),

 R^2_τ, x(β_1, β_0) by 

    R̂^2_τ, x(β_1, β_0)
        = 
    V̂_ττ^-1(β_1, β_0)
    {
    r_1^-1s_Y(1;β_1)|x^2+r_0^-1s_Y(0;β_0)|x^2
     . 
    
                       . 
     - (s_Y(1;β_1),x-
    s_Y(0;β_0),x ) (S_x^2 )^-1
    (s_x,Y(1;β_1)-s_x,Y(0;β_0) )
    },

and the asymptotic distribution of n^1/2{τ̂(β_1, β_0) - τ} by 

    V̂_ττ^1/2(β_1, β_0) 
    [
    {1-R̂^2_τ, x(β_1, β_0)}^1/2·ε + {R̂^2_τ, x(β_1, β_0)}^1/2·  L_K,a] .

The estimated distribution (<ref>) provides a basis for constructing confidence intervals for τ. However,  it is not convenient for theoretical analyses. Below we find the probability limit of (<ref>). 



Under ReM and Conditions <ref> and <ref>,  

the estimated distribution of τ̂(β_1, β_0) in (<ref>) has the same limit as 

    {
    V_ττ(1-R^2_τ, w)
    +
    S^2_τ∖w
    + 
    (r_1r_0)^-1(-)'
    S^2_w∖x
    	(-)
    }^1/2·ε
           + 
    { 
    (r_1r_0)^-1
    (-)'
    S^2_w|x
    (-)
    }^1/2· L_K,a.




From Corollary <ref> and Theorem <ref>, (<ref>) differs from the true asymptotic distribution (<ref>)  in S^2_τ∖w. We can not estimate S^2_τ∖w consistently using the observed data. Consequently, the probability limit has wider quantile ranges than the true asymptotic distribution, which results in conservative confidence intervals. This kind of conservativeness is a feature of finite population causal inference known ever since <cit.>'s seminal work. See the Supplementary Material for a rigorous proof of the conservativeness.
We  discuss two special cases of Theorem <ref> below. 



  §.§.§ Special case: regression adjustment under the CRE



Again, the CRE is ReM with a=∞ and x = ∅.  
Condition <ref> holds automatically under the CRE. 
Theorem <ref> immediately implies 
the following result. 



Under the CRE and Condition <ref>,  the estimated distribution of τ̂(β_1, β_0) in (<ref>) has the same limit as 

    { 
    V_ττ(1-R^2_τ, w) 
    +
    S^2_τ∖w
    + 
    (r_1r_0)^-1(γ-γ̃)'
    S^2_w
    	(γ-γ̃)
    }^1/2·ε.





 




  §.§.§ Special case: no covariate adjustment under ReM

Using Theorem <ref> with β_1 = β_0 = 0, 
we can immediately obtain the probability limit of the estimated distribution of τ̂ under ReM. 

Under ReM and Conditions <ref> and <ref>, the estimated distribution of τ̂ in (<ref>) has the same limit as 

    {V_ττ(1-R^2_τ,x) + S^2_τ∖w}^1/2·ε + 
    ( V_ττ R^2_τ,x)^1/2· L_K,a.






 §.§ General scenarios with partial knowledge of the design


We then consider scenarios without Condition <ref>. The analyzer either does not have all the covariate information used in the design or does not know the balance criterion for ReM. 
We can still estimate V_ττ(β_1, β_0) by V̂_ττ(β_1, β_0) in (<ref>). 
However, we cannot consistently estimate R^2_τ, x(β_1, β_0) due to incomplete information of the covariates used in the design stage. 

We can underestimate R^2_τ, x(β_1, β_0) by zero,
and then estimate the sampling distribution of τ̂(β_1, β_0) by 

    V̂^1/2_ττ(β_1, β_0) ·ε .


An important fact is that the lengths of quantile ranges of the asymptotic distribution (<ref>) are nonincreasing in R^2_τ, x(β_1, β_0). This fact guarantees that the estimated distribution
(<ref>) provides a conservative variance estimator of τ̂(β_1, β_0) and conservative confidence intervals for τ. See the Supplementary Material for a rigorous proof of the conservativeness.
Moreover, (<ref>) equals (<ref>) with a=∞, the estimated distribution pretending that ReM does not happen in the design stage. 
Consequently, the probability limit of (<ref>) equals (<ref>), as the following theorem states. 


 
 


Under ReM and Condition <ref>, 
the estimated distribution of τ̂(β_1, β_0) in (<ref>) has the same limit as (<ref>).















From Theorems <ref> and <ref>, under ReM, the limit of the estimated distribution of τ̂(β_1, β_0) in (<ref>) differs from the corresponding true asymptotic distribution in S^2_τ∖w and R^2_τ,x(β_1, β_0). Neither S^2_τ∖w nor R^2_τ,x(β_1, β_0) have consistent estimators based on the observed data. The difficulty in estimating  S^2_τ∖w comes from the fact that for each unit we can observe  at most one potential outcome, but the difficulty in estimating R^2_τ,x(β_1, β_0) comes from the incomplete information of the design.


Compared to (<ref>), 
the estimated distribution (<ref>) has unnecessarily wider quantile ranges due to the lack of information for ReM, including R^2_τ, x(β_1, β_0), K and a. 
In  (<ref>), we conduct conservative inference and consider the worse-case scenario, which is the CRE with x = ∅ and a = ∞. 
In practice, we can also conduct sensitivity analysis and check how the conclusions change as  R^2_τ, x(β_1, β_0), K and a vary. Without additional information, we still use (<ref>) to construct confidence intervals in the next section. 









§ 
OPTIMAL ADJUSTMENT BASED ON THE ESTIMATED DISTRIBUTION




𝒮-optimality is based on the uncertainty of the sampling distribution. As discussed in Section <ref>, we may not have consistent estimators for the sampling distributions, and often report conservative confidence intervals based on the estimated distributions. Because the lengths of confidence intervals provide important measures of uncertainty in frequentist inference, we focus on the optimal choice of (β_1, β_0) based on the estimated distributions proposed in Section <ref>. 

 
Under ReM, among the estimators in (<ref>), a choice of (β_1, β_0) is optimal  in terms of the estimated precision if the estimated distribution based on (<ref>) or (<ref>) has the
shortest asymptotic 1-α  quantile range for any α∈ (0,1). The corresponding estimated distribution then must have the smallest variance. Conversely, if the estimated distributions are Gaussian, then the one with the smallest estimated variance is optimal  in terms of the estimated precision. 

 
Apparently, the optimal adjusted estimator  in terms of the estimated precision
depends on the approach to constructing confidence intervals. 
In the ideal case, we want the probability limits of the estimated distributions to be identical to the true sampling distributions. Section <ref>, however, shows that this is generally impossible due to treatment effect heterogeneity or the information only known to the designer. 
In some cases, the confidence intervals based on (<ref>) and (<ref>) are asymptotically exact given the analyzer's observed information. For example, if 
the residual from the linear projection of individual treatment effect on covariates w is constant
across all units,  
or equivalently S^2_τ∖w=0,  
(<ref>) provides asymptotically exact confidence intervals; 
if further the designer conducts CRE (i.e., ReM with x = ∅ or a =∞), (<ref>) provides asymptotically exact confidence intervals. In general, these confidence intervals can be conservative, and can be improved <cit.>. 
We focus on (<ref>) and (<ref>) for their simplicity,

and divide this section to two subsections in parallel with Section <ref>. 
 
 



 §.§ When the analyzer knows all the information in the design stage


From Theorem <ref>, we can obtain the optimal adjusted estimator in terms of the estimated precision.  


Under ReM and Conditions <ref> and <ref>, based on Theorem <ref>, 
the optimal adjusted estimator in terms of the estimated precision
is attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 ,  with the estimated distribution having the same limit as 

    {
    V_ττ(1-R^2_τ, w)
    +
    S^2_τ∖w}^1/2·ε.
  




Again, from Proposition <ref>, an optimal choice is (  β_1, β_0) = ( β̃_1, β̃_0 ) under Condition <ref>. Importantly, Corollary <ref> does not depend on covariates x and the threshold a used in ReM, and thus also holds under the CRE.

 




 §.§ General scenarios without Condition <ref>

From Theorem <ref>, we can obtain the optimal adjusted estimator in terms of the estimated precision. 
 

Under ReM and Condition <ref>,  based on Theorem <ref>, the optimal adjusted estimator in terms of the estimated precision
is attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 ,  with  the estimated distribution having the same limit as (<ref>). 


The optimal estimators in terms of the estimated precision are identical in Corollaries <ref> and <ref>, no matter whether the analyzer knows all the information in the design or not. 
The optimal adjustment τ̂(β̃_1, β̃_0) can never hurt the estimated precision. In contrast, Section <ref> shows that  τ̂(β̃_1, β̃_0) may hurt the sampling precision in general. This is an important difference between the two notations of optimality in terms of the sampling precision and estimated precision. 
 

Below we give some intuition for Corollary <ref>. 
Under general scenarios without Condition <ref>, 
the analyzer does not know the information of the design. 
S/he pretends that the design was a CRE, and estimates the sampling distributions of the estimators under the CRE. Luckily, the resulting confidence intervals are still conservative. 
Dropping the term S^2_τ∖w, the estimated distribution converges to the sampling distribution under the CRE. Based on the discussion of 𝒮-optimality under the CRE in Section <ref>, 
the adjusted estimator 
τ̂(β̃_1, β̃_0)  is also optimal in terms of the estimated precision. 
 

 

example-1
[continued]
We revisit Example <ref>, and study the estimated distributions of τ̂(β̃_1, β̃_0) and τ̂ under ReM. 
Because the estimated distributions are Gaussian from Theorem <ref>, 
it suffices to report the estimated standard errors.  
Table <ref> shows the sampling standard errors and the average estimated standard errors.
On average, τ̂(β̃_1, β̃_0) results in shorter confidence intervals than τ̂ when ρ=0 or ρ=0.9. Interestingly, when ρ=0, τ̂(β̃_1, β̃_0) has less sampling precision as Figure <ref>(b) and Table <ref> show. 









§ GAINS FROM THE ANALYZER AND THE DESIGNER


In the design stage, we can use the CRE or ReM. In the analysis stage, we can use the unadjusted estimator τ̂ or the adjusted estimator τ̂(β̃_1, β̃_0). 
Based on the results in previous sections, we now study the additional gains of the designer and analyzer  in 
the sampling precision, the estimated precision, and the coverage probability. 




 §.§ Sampling precision


We first study the additional gains in the sampling precision from the analyzer and the designer, respectively. 
We measure the additional gain of the analyzer by comparing the asymptotic distributions of τ̂ and 
τ̂(β̃_1, β̃_0) under ReM. 
We measure the additional gain of the designer by comparing the asymptotic distributions of τ̂(β̃_1, β̃_0) under the CRE and ReM. 
Similar to Section <ref>, 
we consider different scenarios based
on the relative amount of covariate information for the analyzer and the designer.
Let v_K,a = P(χ^2_K+2≤ a)/P(χ^2_K≤ a)≤ 1 be the variance of L_k,a <cit.>, and q_1 - α/2(ρ^2) be the (1-α / 2)th quantile of (1-ρ^2)^1/2·ε+|ρ| · L_K,a. 





  §.§.§ The analyzer has no less covariate information than the designer


 

First, we measure the additional gain of the analyzer.



	Under ReM and Conditions <ref> and <ref>, compare τ̂(β̃_1, β̃_0) to τ̂.
	The percentage reduction in the asymptotic variance is 
	{  R^2_τ,w-(1-v_K,a)R^2_τ,x} / { 1-(1-v_K,a)R^2_τ,x} .
	For any α∈ (0,1), the percentage reduction in the length of the asymptotic 1-α  quantile range is 
	1 - (1-R^2_τ, w)^1/2·  q_1-α/2(0) / q_1-α/2(R^2_τ,x) .
	Both percentage reductions are nonnegative and nondecreasing in R^2_τ,w. 





From Corollary <ref>,  the gain from the analyzer through regression adjustment is nondecreasing in the analyzer's covariate information. Both percentage reductions in Corollary <ref> converge to 1 as R^2_τ,w converges to 1. 



Second, we measure the additional gain of the designer. Section <ref> demonstrates that τ̂(β̃_1, β̃_0) has the same asymptotic distribution under the CRE and ReM. Therefore, under Condition <ref>, the gain from the designer is zero. Nevertheless, this also implies that using ReM in the design will not hurt the sampling precision of τ̂(β̃_1, β̃_0). Moreover, under ReM, we can use the additional covariate information, in the same way as in the CRE, to improve the sampling precision. 



  §.§.§ The analyzer has no more covariate information than the designer


 

First, we measure the additional gain of the analyzer.
Both the asymptotic distributions of τ̂(β̃_1, β̃_0) and τ̂ are linear combinations of ε and L_K,a. The coefficients of ε are identical, but the coefficient of L_K,a for τ̂(β̃_1, β̃_0) is smaller than that for τ̂. 



	Under ReM and Conditions <ref> and <ref>, compare τ̂(β̃_1, β̃_0) to τ̂. 
	The percentage reduction in the asymptotic variance is
	v_K,aR^2_τ, w / { 1-(1-v_K,a)R^2_τ,x} .
	For any α∈ (0,1), the percentage reduction in the length of the asymptotic 1-α quantile range is
	1- (    1-R^2_τ, w )  ^1/2· 
	q_1-α/2(ρ^2_τ, x∖w) / q_1-α/2(R^2_τ, x) .
	Both percentage reductions are nonnegative and nondecreasing in R^2_τ,w. 



From Corollary <ref>, the improvement from regression adjustment is nondecreasing in the analyzer's covariate information. However, this improvment is small when the designer uses a small threshold a for ReM. Both percentage reductions in Corollary <ref> converge to 0 as a converges to 0. Intuitively, when the designer uses ReM with a small threshold, s/he has used more covariate information thoroughly in the design, and thus the analyzer has only a small additional gain through regression adjustment. 

Second, we measure the additional gain of the designer.



	Under Conditions <ref> and <ref>, compare τ̂(β̃_1, β̃_0) under ReM to
	that under the CRE. The percentage reduction in the asymptotic variance  is
	(1-v_K,a)ρ^2_τ, x∖w.
	For any α∈ (0,1), the percentage reduction in the length of the asymptotic 1-α quantile range is
	1-   q_1-α/2(ρ^2_τ, x∖w) / q_1-α/2(0) .
	Both percentage reductions are nonnegative and nondecreasing in R^2_τ, x.





From Corollary <ref>, the gain from the designer through ReM is nondecreasing in the designer's covariate information. The gain from the designer is substantial when R^2_τ, x is large and the threshold for ReM is small.
Both percentage reductions in Corollary <ref> converge to 1 as R^2_τ, x→ 1 and a→ 0.


Consider the special case where the designer and the analyzer have the same covariates (x = w). The additional gain from the analyzer is small given that the designer uses ReM with a small threshold a, and so is the additional gain from the designer given that the analyzer uses the 𝒮-optimal τ̂(β̃_1, β̃_0). 





  §.§.§ General scenarios



The complexity discussed in Section <ref> makes it difficult to evaluate the additional gains from the analyzer and the designer. Given any adjusted estimator, the designer can always use ReM to reduce the asymptotic variance and the lengths of asymptotic quantile ranges <cit.>. For the analyzer, in general, the performance of regression adjustment under ReM depends on the covariates used in the design, and thus the analyzer does not know the optimal adjusted estimator among (<ref>). For instance, without the covariate information used in the design, the analyzer is not sure whether τ̂(β̃_1, β̃_0) has smaller asymptotic variance than τ̂. Example <ref> shows two cases where τ̂(β̃_1, β̃_0) increases and decreases the sampling precision, respectively. 





 §.§ Estimated precision


We then study the additional gains in the asymptotic estimated precision from the analyzer and the designer, respectively. 
We measure the additional gain of the analyzer by comparing the estimated distributions of τ̂(β̃_1, β̃_0) and  τ̂ under ReM.
We measure the additional gain of the designer by comparing the estimated distributions of τ̂(β̃_1, β̃_0) under the CRE and ReM.
Similar to Section <ref>, we consider two scenarios based on whether the analyzer has full knowledge of the design or not.
Let 
κ=1+V^-1_ττS^2_τ∖w≥ 1,
which reduces to 1 when S^2_τ∖w = 0, i.e., the “adjusted” individual treatment effect τ_i(β̃_1, β̃_0) =  Y_i(1; β̃_1) - Y_i(0; β̃_0) is constant for all units.




  §.§.§ When the analyzer knows all the information in the design stage




First, we measure the additional gain of the analyzer.   


	Under ReM and Conditions <ref> and <ref>, compare the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) to that of τ̂ based on (<ref>). 
	The percentage reduction in the variance is
	{  R^2_τ,w - (1-v_K,a)R^2_τ, x} / {κ - (1-v_K,a)R^2_τ, x}.
	For α∈ (0,1), the percentage reduction in the length of the 1-α quantile range is 
	1 - ( 1-R^2_τ, w/κ)^1/2· 
	q_1-α/2(0) / q_1-α/2(R^2_τ,x/κ) .
	Both percentage reductions are nonnegative and nondecreasing in R^2_τ,w. 


From Corollary <ref>, 
the gain from regression adjustment is nondecreasing in the analyzer's covariate information. 
Both percentage reductions in Corollary <ref> converge to 1 as R^2_τ,w→ 1 and κ→ 1. 

Second, we measure the additional gain of the designer. From Corollary <ref> and the comment after it, the probability limits of
the estimated distributions of τ̂(β̃_1, β̃_0) are identical under both designs. 
Therefore, the gain from the designer is zero. 



	Consider the special case where the analyzer has the same covariate information as the designer and knows the balance criterion in the design.
As discussed above, the designer has no gain.	
	 Based on Corollaries <ref> and  <ref>, with a small threshold a, the estimated distributions of τ̂ and 
	 τ̂(β̃_1, β̃_0) 
	 have approximately the same probability limit. Therefore, the analyzer has small additional gain.





  §.§.§ General scenarios without Condition <ref>



First, 
we measure the additional gain of the analyzer. 


	Under ReM and Condition <ref>, compare the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) to that of τ̂ based on (<ref>). 
	The percentage reduction in the variance is R^2_τ, w/κ.  For any α∈(0,1), the percentage reduction in the length of the  1-α quantile range is 1-( 1-R^2_τ, w/κ)^1/2 . Both percentage reductions are nonnegative and nondecreasing in R^2_τ, w. 




Both percentage reductions in Corollary <ref> converge to 1 as R^2_τ, w→ 1 and κ→ 1, and they are larger than or equal to those in Corollary <ref>. 
 


Second, we measure the additional gain of the designer. From Corollary <ref> and Theorem <ref>, the estimated distributions of any adjusted estimator in (<ref>) have the same probability limit under both designs. Therefore, 
the gain from the designer  
is zero. Nevertheless, using ReM will not hurt the estimated precision of the treatment effect estimators, and we can use covariates in the analysis in the same way as in the CRE to improve the estimated precision. 





 §.§ Coverage probabilities


 

From Sections <ref> and <ref>,  (a) τ̂(β̃_1, β̃_0) is optimal in terms of the estimated precision no matter whether the analyzer knows all the information of the design or not, and (b) the designer provides no gain in the estimated precision of the τ̂(β̃_1, β̃_0). 
From (a), under ReM, the analyzer can never increase the asymptotic lengths of the confidence intervals by using τ̂(β̃_1, β̃_0) instead of τ̂. Therefore, we do not measure the additional gain of the analyzer in coverage probabilities. From (b), 
the asymptotic lengths of the confidence intervals based on τ̂(β̃_1, β̃_0) are the same under the CRE and ReM.   However, we will show shortly that  the designer can help to improve the coverage probabilities of the confidence intervals based on τ̂(β̃_1, β̃_0).  

 



  §.§.§ When the analyzer knows all the information in the design stage


From Sections 

<ref> 
and <ref>, 


τ̂(β̃_1, β̃_0) has the same sampling precision and estimated precision under ReM and the CRE.
Therefore,  the coverage probabilities of the associated confidence intervals  are asymptotically the same under the CRE and ReM.  
This implies that the designer provides no gain for the coverage probabilities of the confidence intervals based on τ̂(β̃_1, β̃_0).
We formally state the results as follows. 


	Under Conditions <ref> and <ref>, compare the confidence intervals based on τ̂(β̃_1, β̃_0) and (<ref>) under the CRE and ReM. 
	Their lengths are asymptotically the same after being scaled by n^1/2, and they have the same asymptotic coverage probability. 









  §.§.§ General scenarios without Condition <ref>



From Sections 

<ref>
and 
<ref>, 
τ̂(β_1, β_0) in (<ref>) with any β_1 and β_0
has better sampling precision under ReM than under the CRE, 
but 
it has the same estimated precision under ReM and the CRE. 
Therefore, 
the confidence intervals based on τ̂(β_1, β_0) under ReM have higher coverage probabilities than that under the CRE. 
We give a formal statement below. 





	Under Condition <ref>, compare the confidence intervals based on τ̂(β_1, β_0) and (<ref>) under the CRE and ReM. 
	Their lengths are asymptotically the same after being scaled by n^1/2. However, the asymptotic coverage probability under ReM is larger than or equal to that under the CRE. 


In Corollary <ref>, the confidence intervals under both ReM and the CRE are asymptotically valid and of the same length, 
but the one under ReM has higher coverage probabilities and is more conservative.
In particular,  for any α∈ (0,1) and any (β_1, β_0), 
as R^2_τ, x(β_1, β_0) → 1, S^2_τ∖w→ 0 and a → 0, 
the asymptotic coverage probabilities of the 1-α confidence intervals are 1 and 1-α under ReM and the CRE, respectively. 
Corollary <ref> holds for any adjusted estimator and thus holds for 
τ̂(β̃_1, β̃_0). Therefore, under general scenarios without Condition <ref>, the designer can provide a substantial gain in coverage probabilities of confidence intervals. This gives another justification for using ReM. 
 



 


§ UNIFICATION AND PRACTICAL SUGGESTIONS







 §.§ Unification


In total, there are four combinations in the design and analysis of experiments. Figure <ref> summarizes the sampling distributions and the probability limits of the estimated distributions for all combinations. 





<cit.> started the literature by discussing the property of τ̂ under the CRE. <cit.> showed that τ̂(β̃_1, β̃_0) improves τ̂ in terms of the sampling precision and estimated precision under the CRE. Arrow (i) in Figure <ref> illustrates this improvement. <cit.> showed that ReM improves the CRE in terms of the sampling precision and the estimated precision of τ̂. Arrow (ii) in Figure <ref> illustrates this improvement. Interestingly, τ̂(β̃_1, β̃_0) under the CRE and τ̂ under ReM have almost identical asymptotic sampling distributions and estimated distributions, if we use the same sets of covariates and a≈ 0 in ReM. 


However, both <cit.> and <cit.> compared sub-optimal strategies. We evaluated the additional gain from the analyzer given that the designer uses ReM.  Arrow (iii) in Figure <ref> illustrates this improvement. We also evaluated the additional gain from the designer given that the analyzer uses τ̂(β̃_1, β̃_0). 
Arrow (iv) in Figure <ref> illustrates this improvement. 

Table <ref> summarizes the results under all scenarios. 
We have the following conclusions. 

	
	
  * Compare the analyzer and the designer based on the sampling precision. 
	From the first two rows of Table <ref>, when one has more covariate information than the other, the one with more covariate information provides a substantial additional gain, while the other provides negligible additional gain.  
	
	
	
  * Compare the analyzer and the designer based on the estimated precision. 
	From the 6th and 8th
	columns of Table <ref>,
	the additional gain from the analyzer can be substantial, 
	while the additional gain from the designer is negligible in general. 
	
	
  * Consider the special case where the  analyzer has the same covariate information as the designer and knows the balance criterion in the design. 
	From the fourth row of Table <ref>,
	the additional gain from either the analyzer or the designer are negligible. 
	
	
  * From the last row of Table <ref>,
	the analyzer may hurt the sampling precision through regression adjustment, but can provide a substantial gain in the estimated precision. 
	The designer can improve sampling precision of any adjusted estimator, and does not hurt the estimated precision. 
	Therefore, although the designer cannot shorten the confidence intervals, s/he can increase the coverage probabilities. 
























 §.§ ReM, Lin's estimator, and the Huber–White variance estimator


Based on the summary in Section <ref>, we recommend using ReM in the design and  τ̂(β̃_1, β̃_0) in the analysis, which has better estimated precision and coverage property. 
However, some practical issues remain. 


First, we need to estimate the population OLS coefficients β̃_1 and β̃_0. Under both the CRE and ReM, we can use their sample analogues as consistent estimators, with β̂_1 and β̂_0 being the coefficients of w in the OLS fit of Y on w under the treatment and control, respectively. The corresponding adjusted estimator τ̂(β̂_1, β̂_0) numerically equals the coefficient of Z in the OLS fit of Y on Z, w and Z×w, i.e., it is <cit.>'s estimator. Replacing β̃_1 and β̃_0 by their sample analogues does not change the asymptotic distribution. Informally speaking, τ̂(β̃_1, β̃_0) and τ̂(β̂_1, β̂_0) have the same asymptotic behavior and optimality. The following corollary is a formal statement.



Under ReM and Condition <ref>, τ̂(β̃_1, β̃_0) and τ̂(β̂_1, β̂_0) have the same asymptotic distributions and the same probability limits of the estimated distributions. Thus, among (<ref>),
τ̂(β̂_1, β̂_0) is 𝒮-optimal under Condition <ref> or <ref>, 

and is
always optimal in terms of the estimated precision based on (<ref>) and (<ref>). 


 


Second, the estimated distributions (<ref>) and (<ref>) are identical with or without Condition <ref> because R̂^2_τ, x(β̂_1, β̂_0) based on (<ref>) equals zero under Condition <ref>. 
Moreover,
the variance estimator V̂_ττ(  β̂_1, β̂_0 ) based on (<ref>) is asymptotically equivalent to the Huber–White variance estimator V̂_ of the coefficient of Z from the OLS fit of Y on Z, w and Z×w. 





Under ReM and Condition <ref>, V̂_ττ(  β̂_1, β̂_0 ) - V̂_→ 0 in probability. 


Theorem <ref> extends <cit.>'s result for the CRE to ReM. It requires only Condition <ref>, but <cit.> requires higher order moment conditions. 
We finally construct the Wald-type confidence intervals based on a Gaussian approximation. 
The statistical inference based on τ̂(β̂_1, β̂_0), including variance estimation and confidence interval construction, is always the same no matter whether the design is a CRE or ReM and no matter whether the analyzer knows all the information of the design or not. 

From the above, using ReM and τ̂(β̂_1, β̂_0) enjoys the  
optimal estimated precision
and improves the coverage property, 
and the associated statistical inference can be conveniently implemented through the OLS fit and Huber–White variance estimate. 




§ ILLUSTRATION




 §.§ A simulation study

 
We conduct a simulation study to investigate the performance of the asymptotic approximation and the coverage properties of the confidence intervals in finite samples.  We generate the data in the same way as in Example <ref> with ρ=0 and vary the sample size n from 100 to 1000. 
For each simulated dataset, we generate ReM based on covariate x, where the the threshold a is the 0.001th quantile of the χ^2_1 random variable. 
Figure <ref>(a)–(b) show the histograms of τ̂ and τ̂(β̂_1, β̂_0) based on covariate w.
From Figure <ref>(a)–(b), the asymptotic approximation works fairly well. 
We then construct 95% confidence intervals for the average treatment effect, using the estimated distribution (<ref>) with 
either τ̂ or τ̂(β̂_1, β̂_0). 
From Figure <ref>(c)–(d), 
the confidence intervals based on estimator adjusted for w are shorter than that based on τ̂, 
and 
both confidence intervals are conservative with coverage probabilities larger than the nominal level, due to the analyzer's incomplete information of the design. 
We further consider τ̂(β̂_1, β̂_0) based on (x,w), assuming that the analyzer has access to the covariate x in the design. 
The corresponding confidence interval is even shorter and  becomes asymptotically exact, due to the additive treatment effects in the data generating process. 
From Figure <ref>(d),  the coverage probabilities are close to the nominal level as the sample size increases. 
Even when the sample size is small, the confidence interval works fairly well with coverage probability at least 94%. 
 










 §.§ The “Opportunity Knocks” experiment

The “Opportunity Knocks” experiment <cit.> aims at evaluating the impact of a financial incentive demonstration program on college students' academic performance. The experiment includes first- and second-year students who apply for the financial aid at a large Canadian commuter university. These students were randomly assigned to treatment and control groups. Students in the treated group have peer advisors and receive cash reward for attaining certain grades.


We use this dataset to illustrate rerandomization and regression adjustment. We consider the second-year students, 
and choose the outcome to be the average grade for the semester right after the experiment. 
We exclude students with missing outcomes or covariates, resulting in a treatment group of size 199 and a control group of size 369. 
We evaluate the repeated sampling properties of the adjusted estimators under rerandomization,  which depend on all the potential outcomes. 
However, half of the potential outcomes are missing from the observed data.  
To make the simulation more realistic, we impute all the missing potential outcomes based on a simple model fitting. 
We first fit a linear model of the observed outcome on the treatment indicator and covariates within each stratum classified by sex and high school GPA. We then impute the missing potential outcomes using the fitted linear model. 

We conduct ReM with two covariates  
and choose threshold a to be the 0.005th quantile of χ^2_2.
For the covariates in the design and the analysis, 
we consider the following two cases: 

	
  * the covariates in the design are sex and high school grade, and
	the  covariates in the analysis are 
	whether mother/father is a college graduate, whether correctly answer the first/second question in a survey, whether mother tongue is English, and  GPA in the previous year;
	
  * the covariates in the design and analysis are the same as in case (i), except that we switch the high school grade to the analysis stage and switch the GPA in the previous year to the design stage.

 
 


We first consider the sampling precision. 
Table <ref> shows the coefficients of ε in the asymptotic distributions. We omit the coefficients of L_K,a because the ε components are the dominating terms in the asymptotic distributions. 
Figure <ref> shows the histograms of τ̂(β̂_1, β̂_0) and τ̂ under ReM.
In Table <ref>,  compared to the second column, the reduction in coefficients in the first column shows the gain from the designer alone, and the reduction in the last column shows the gain from the analyzer alone. The magnitude of the reduction suggests the relative amount of the covariate information of the designer and analyzer.
In case (i), 
the first row of Table <ref> shows that τ̂(β̂_1, β̂_0) under the CRE is more precise than τ̂ under ReM. This holds because
R^2_τ,w = 0.65 ≥ R^2_τ,x = 0.31. 
Figure <ref>(a) shows that 
τ̂(β̂_1, β̂_0) outperforms τ̂, coherent with Theorem <ref>. 
In case (ii), 
the second row of Table <ref> shows that τ̂(β̂_1, β̂_0) under the CRE is less precise than τ̂ under ReM. This is because 
R^2_τ,w = 0.32 < R^2_τ,x = 0.65. 
Figure <ref>(b) shows that 
τ̂ outperforms τ̂(β̂_1, β̂_0) under ReM. 


We then consider the estimated precision. Because the estimated distributions in (<ref>) are Gaussian, it suffices to compare the estimated standard errors. Table <ref> shows the average estimated standard errors and the coverage probabilities of the 95% confidence intervals. 
In both cases, τ̂(β̂_1, β̂_0) has almost the same estimated precision under ReM and the CRE. So does τ̂. However, the 95% confidence intervals under ReM have higher coverage probability, which is coherent with Corollary <ref>.
Moreover, τ̂(β̂_1, β̂_0) always has higher estimated precision than τ̂, which is coherent with the 
optimality results in Corollary <ref>. 


































 




§ DISCUSSION



In sum, regression adjustment can improve the estimated precision but may hurt the sampling precision, and ReM can improve the sampling precision and  never hurts the estimated precision. 
The resulting adjusted estimator is optimal 
in terms of the estimated precision among all linearly adjusted estimators in (<ref>), 
has lower sampling variability under ReM than it would have had under the CRE, 
and the corresponding confidence intervals have higher coverage probabilities than that under the CRE.
Therefore, in practice, we recommend using ReM in the design and using <cit.>'s estimator in the analysis followed by the Huber–White robust standard error. Importantly, the analyzer should communicate with the designer, asking for detailed covariate information and assignment mechanism in the design stage. 



For the analysis, we focused on inferring the average treatment effect using regression adjustment. 
It is interesting to extend the discussion to covariate adjustment in more complicated settings, such as high dimensional covariates <cit.>, 
logistic regression for binary outcomes <cit.>, 
and adjustment using machine learning methods <cit.>. 
It is also important to consider covariate adjustment for general nonlinear estimands <cit.> and general designs <cit.>, 
such as blocking <cit.>, matched pairs <cit.>, and factorial designs <cit.>. 

 
For the design, we focused on rerandomization using the Mahalanobis distance. It is conceptually straightforward to extend the results to rerandomization with tiers of covariates <cit.>.  Recently, <cit.> discussed sequential rerandomization, and <cit.> discussed rerandomization in 2^K factorial experiments with tiers of both covariates and factorial effects.  It is important to discuss regression adjustment after these rerandomizations. 


 

The relationship between blocking and post-stratification for discrete covariates is analogous to the relationship between rerandomization and regression adjustment for general covariates. 
When the number of blocks is small compared to the sample size, our results suggest conducting post-stratification, or equivalently an OLS fit of the outcome on treatment, block indicator and their interactions, followed by the Huber–White robust standard error. 
	When the number of blocks is large, 
<cit.> showed that post-stratification can be worse than blocking, 
which sheds light on the possible advantage of rerandomization over regression adjustment with a large number of covariates. In this case, although deriving the asymptotic properties of rerandomization is challenging, it is still straightforward to conduct Fisher randomization tests. 






§ ACKNOWLEDGMENTS

We thank the Associate Editor and two reviewers for constructive comments. Peng Ding gratefully acknowledges financial support from the National Science Foundation (DMS grant # 1713152).

 
plainnat






































Supplementary Material












Appendix <ref> proves   the results related to the sampling distributions.
 
Appendix <ref> proves   the results related to 𝒮-optimality.

Appendix <ref> proves   the results related to the confidence intervals and the optimal adjusted estimators in terms of the estimated precision.

 

Appendix <ref> proves   the results related to the gains from the analyzer and the designer.
 

Appendix <ref> proves the asymptotic equivalence of τ̂(β̃_1, β̃_0) and τ̂(β̂_1, β̂_0) as well as the asymptotic equivalence of V̂_ττ(β̂_1, β̂_0) and the Huber–White variance estimator.









§ SAMPLING DISTRIBUTIONS OF ADJUSTED ESTIMATORS


 
 


By definition,  

    r_0β̃_1+r_1 β̃_0     = 
    r_0 ( S_w^2 )^-1S_w, Y(1) + r_1 (S_w^2 )^-1S_w, Y(0)
        = {(r_1r_0)^-1S^2_w}^-1{  
    r_1^-1S_w, Y(1) + r_0^-1S_w, Y(0)} 
     = V_ww^-1V_wτ = .





The regression adjustment coefficients β_1 and β_0 can depend on sample size n implicitly and have finite limits as n→∞. 
Recall that Y_i(z;β_z) = Y_i(z) - β_z' w_i is the “adjusted” potential outcome  under treatment z, and  τ_i(β_1, β_0) = τ_i - (β_1- β_0)'w_i is the “adjusted” individual treatment effect. 
Under Condition <ref>(ii),  
the finite population variances and covariances 
S^2_Y(z;β_z) = S^2_Y(z) + β_z' S^2_wβ_z - 2β_z'  S_w, Y(z),
S^2_τ(β_1, β_0) = S^2_τ + (β_1- β_0)' S^2_w (β_1- β_0) - 2 (β_1- β_0)' S_w,τ,
S_Y(z;β_z), x = S_Y(z),x - β_z'S_w, x,
and 
S^2_x
have finite limiting values. Under Condition  <ref>(iii), the maximum squared distances satisfy that as n→∞, 
max_1≤ i≤ nx_i_2^2/n → 0, and 

    n^-1max_1≤ i≤ n|Y_i(z;β_z) - Y̅(z;β_z)|^2 
        = n^-1max_1≤ i≤ n|Y_i(z) - Y̅(z) - β_z'w_i|^2 
       ≤ 
    n^-1   (  1+β_z'β_z ) max_1≤ i≤ n( |Y_i(z) - Y̅(z)|^2 + w_i_2^2 ) → 0,

where the inequality follows from the Cauchy–Schwarz inequality. 
Using <cit.>, we can show that, under ReM, n^1/2{τ̂(β_1, β_0)-τ} has the asymptotic distribution (<ref>). 




Corollary <ref> follows from Theorem <ref> with a=∞.



Corollary <ref> follows from Theorem <ref> with ( β_1, β_0 ) = (0,0).





By definition, Cov(τ̂_w, τ̂  - 'τ̂_w ) =0 under the CRE. We have 

    τ̂  - 'τ̂_w    = 
    (τ̂|τ̂_x) + 
    (τ̂|τ̂_x) - 
    '
    {( τ̂_w|τ̂_x) + 
    ( τ̂_w|τ̂_x)
    }
        =  
    {(τ̂|τ̂_x) - _'( τ̂_w|τ̂_x)} + 
    {(τ̂|τ̂_x)- _'( τ̂_w|τ̂_x)
    }
            
    - ( - _ )'( τ̂_w|τ̂_x) - 
    ( - _)'( τ̂_w|τ̂_x).

Thus, the covariances between τ̂_w and the four terms in (<ref>) sum to 0. Below we consider the four covariances separately. 

First, by definition, 
(τ̂|τ̂_x) - _'( τ̂_w|τ̂_x) is uncorrelated with ( τ̂_w|τ̂_x). 
Moreover, 
because (τ̂|τ̂_x) - _'( τ̂_w|τ̂_x) is a linear function of τ̂_x, it must also be uncorrelated with ( τ̂_w|τ̂_x). 
Thus, (τ̂|τ̂_x) - _'( τ̂_w|τ̂_x) is uncorrelated with 
( τ̂_w|τ̂_x) + ( τ̂_w|τ̂_x) = τ̂_w.

Second, 
by definition, 
(τ̂|τ̂_x)- _'( τ̂_w|τ̂_x) is uncorrelated with ( τ̂_w|τ̂_x). 
Moreover, because 
(τ̂|τ̂_x)- _'( τ̂_w|τ̂_x) is uncorrelated with τ̂_x, 
it must also be uncorrelated with (τ̂_w|τ̂_x). 
Thus, (τ̂|τ̂_x)- _'( τ̂_w|τ̂_x) is uncorrelated with 
( τ̂_w|τ̂_x) + 
( τ̂_w|τ̂_x) = τ̂_w.

Third, 
because ( τ̂_w|τ̂_x) is uncorrelated with ( τ̂_w|τ̂_x), 
we can simplify 
the covariance between  ( - _ )'( τ̂_w|τ̂_x) and τ̂_w = ( τ̂_w|τ̂_x) + ( τ̂_w|τ̂_x) as 
( - _ )' {( τ̂_w|τ̂_x) 
}. 
The covariance of ( τ̂_w|τ̂_x) has the following equivalent forms: 

    {
    (
    τ̂_w|τ̂_x
    )
    }    = n^-1V_wxV_xx^-1V_xw = (n r_1r_0)^-1S^2_w|x.














Fourth, because ( τ̂_w|τ̂_x) is uncorrelated with ( τ̂_w|τ̂_x), we can simplify the covariance between ( - _)'( τ̂_w|τ̂_x) and τ̂_w = ( τ̂_w|τ̂_x) + ( τ̂_w|τ̂_x) as 
( - _)' {( τ̂_w|τ̂_x)
}. 
The covariance of ( τ̂_w|τ̂_x) has the following equivalent forms: 

    {
    (
    τ̂_w|τ̂_x
    )
    }    = (τ̂_w) - {
    (
    τ̂_w|τ̂_x
    )
    }  = (n r_1r_0)^-1S^2_w∖x.













From the above, the zero covariance between τ̂_w and (<ref>) implies that 

    0    = -(nr_1r_0)^-1{( - _ )' S^2_w|x + 
    ( - _ )' S^2_w∖x}.

Therefore, Proposition <ref> holds. 





First, by the definitions of V_ττ(β_1, β_0)  in (<ref>) and R^2_τ,x(β_1, β_0) in (<ref>), the squared coefficients of L_K,a and ε in (<ref>) have the following equivalent forms:

    V_ττ(β_1, β_0) R^2_τ,x(β_1, β_0)     =
    n{(τ̂(β_1,β_0)|τ̂_x) }, 
    
    V_ττ(β_1, β_0){1- R^2_τ,x(β_1, β_0) }    = n{(τ̂(β_1,β_0)|τ̂_x) 
    }.


Second, 
because τ̂(β_1,β_0) = τ̂ - 'τ̂_w by (<ref>), the linear projection of τ̂(β_1,β_0) on τ̂_x under the CRE and the corresponding residual have the following equivalent forms: 

    (τ̂(β_1,β_0)|τ̂_x)
        = 
    (
    τ̂ - 'τ̂_w|τ̂_x) = 
    (
    τ̂|τ̂_x)
    - 
    '
    (
    τ̂_w|τ̂_x),
    
    (τ̂(β_1,β_0)|τ̂_x)
        = (
    τ̂ - 'τ̂_w|τ̂_x)
    = 
    (
    τ̂|τ̂_x)
    - 
    '
    (
    τ̂_w|τ̂_x).

Using the definitions of _ in (<ref>) and _ in (<ref>), we can express the above quantities as 

    (τ̂(β_1,β_0)|τ̂_x)
        = (
    τ̂|τ̂_x)
    - 
    _'
    (
    τ̂_w|τ̂_x) - 
    (
    
    -_)' 
    (
    τ̂_w|τ̂_x)
    
        =
    τ + 
    {(
    τ̂|τ̂_x
    ) |(
    τ̂_w|τ̂_x
    )
    } - 
    (
    
    -_
    )' 
    
    (
    τ̂_w|τ̂_x
    ), 
    
    (τ̂(β_1,β_0)|τ̂_x)
        = 
    (
    τ̂|τ̂_x)
    - 
    _'
    (
    τ̂_w|τ̂_x)
    - 
    ( - _)'
    (
    τ̂_w|τ̂_x)
    
        = 
    {(
    τ̂|τ̂_x) |(
    τ̂_w|τ̂_x)
    } - 
    ( - _)'
    (
    τ̂_w|τ̂_x).


Third, 
because the two terms in (<ref>) excluding the constant term τ are uncorrelated, the variance of 
the linear projection of τ̂(β_1,β_0) on τ̂_x is the summation of the variances of these two terms in (<ref>). 
Using (<ref>) and the definitions of R^2_τ,x and R^2_ in (<ref>) and (<ref>), we have 

    {(τ̂(β_1,β_0)|τ̂_x)
    }
        =    {(
    τ̂|τ̂_x)
    }·(
    1-R^2_) + 
    (
    
    -_)' 
    {(
    τ̂_w|τ̂_x)}(
    
    -_)
    
        =     
    n^-1V_ττ· R^2_τ,x(
    1-R^2_) + 
    (n r_1r_0)^-1(
    
    -_)' 
    S^2_w|x(
    
    -_).

Similarly, because the two terms in (<ref>) are uncorrelated, the variance of the residual of 
the linear projection of τ̂(β_1,β_0) on τ̂_x is the summation of the variances of these two terms in (<ref>). 
Using (<ref>) and the definitions of R^2_τ,w in (<ref>) and R^2_ in (<ref>), we have 

    {(τ̂(β_1,β_0)|τ̂_x)
    }
        =    {(
    τ̂|τ̂_x) 
    }·(
    1-R^2_) + 
    ( - _)'
    {(
    τ̂_w|τ̂_x)}( - _)
    
        =      
    n^-1V_ττ·(
    1 - R^2_τ,x)(
    1-R^2_) + 
    (n r_1r_0)^-1( - _)'
    S^2_w∖x( - _).


Fourth, 
using (<ref>), (<ref>), (<ref>) and (<ref>), we have   

    V_ττ(β_1, β_0) R^2_τ,x(β_1, β_0)     =
    V_ττ R^2_τ,x(
    1-R^2_) + 
    (r_1r_0)^-1(
    
    -_)' 
     S^2_w|x(
    
    -_), 
    
    V_ττ(β_1, β_0){1- R^2_τ,x(β_1, β_0) }    = 
    V_ττ(
    1 - R^2_τ,x)(
    1-R^2_) + 
    (r_1r_0)^-1( - _)'
    S^2_w∖x( - _).

These coupled with Theorem <ref> imply
Theorem <ref>. 





First, we prove that _=_= and R^2_=1.
Under Condition <ref>, τ̂_w can linearly represent τ̂_x as 
τ̂_x = B_1τ̂_w. Thus using the linearity of the projection operator, we have 

    (τ̂|τ̂_x) - τ -  '(τ̂_w|τ̂_x) = 
    (τ̂  - τ - 'τ̂_w|τ̂_x)
    = 
    {(τ̂|τ̂_w) |τ̂_x}.

By definition, (τ̂|τ̂_w) is uncorrelated with τ̂_w, and then it is also uncorrelated with τ̂_x = B_1τ̂_w. Therefore, 
(<ref>) equals zero, 
implying that 
(i) (τ̂|τ̂_x) = τ +  '(τ̂_w|τ̂_x) and  equals the linear projection coefficient of (τ̂|τ̂_x) on (τ̂_w|τ̂_x),  i.e., =_; 
(ii) the squared multiple correlation between (τ̂|τ̂_x) and (τ̂_w|τ̂_x) equals 1, i.e., R^2_=1.  
Moreover, (i) and Proposition <ref> imply that = _=_. 

Second, we prove that R^2_ =  (R^2_τ,w - R^2_τ,x)/(1-R^2_τ,x). 
Because _=, 
the residual from the linear projection of (τ̂|τ̂_x) on (τ̂_w|τ̂_x)
reduces to  

    (τ̂|τ̂_x) - _'(τ̂_w|τ̂_x)     = 
    (
    τ̂- 'τ̂_w|τ̂_x) = 
    {(τ̂|τ̂_w)
    |τ̂_x}.

Because, under Condition <ref>, τ̂_x = B_1 τ̂_w is uncorrelated with (τ̂|τ̂_w), 
(<ref>) reduces to (τ̂|τ̂_w ). 







Thus, the squared multiple correlation between 
(τ̂|τ̂_x )  and (τ̂_w|τ̂_x) reduces to 

    R^2_    =  
    1 - {(τ̂|τ̂_w )
    }}/{(τ̂|τ̂_x) }
     = 1 - 1-R^2_τ,w/1-R^2_τ,x = 
    R^2_τ,w-R^2_τ,x/1-R^2_τ,x.


Corollary <ref> then follows immediately from Theorem <ref>.




First, we prove that β̃_=β̃, R^2_= R^2_τ, w/R^2_τ, x, S^2_w|x=S^2_w, and S^2_w∖x=0. 
Under Condition <ref>, 
S^2_w|x=S^2_w,  S^2_w∖x=0, 
and 
τ̂_w = B_2 τ̂_x. 
Then using Proposition <ref>, we have _=. 
Moreover, 
(τ̂_w|τ̂_x) = τ̂_w, and thus the linear projection of 
(τ̂|τ̂_x) on (τ̂_w|τ̂_x) under the CRE reduces to  

    {(τ̂|τ̂_x) |(τ̂_w|τ̂_x)
    }
    = τ + _'
    (τ̂_w|τ̂_x) 
    = 
    τ + '
    τ̂_w = (τ̂|τ̂_w ).

Consequently, 
the squared multiple correlation between (τ̂|τ̂_x) and  (τ̂_w|τ̂_x) equals 

    R^2_    = 
    {(τ̂|τ̂_w )
    }/{(τ̂|τ̂_x)
    }
    = 
    (τ̂)R^2_τ, w/(τ̂)R^2_τ, x
    = 
    R^2_τ, w/R^2_τ, x.



Second, we prove that R^2_ = 0.
Under Condition <ref>, 
(τ̂_w|τ̂_x) = 0,
and thus the squared multiple correlation between (τ̂|τ̂_x) and 
(τ̂_w|τ̂_x) reduces to zero.

Corollary <ref> then follows immediately from Theorem <ref>.







§ 𝒮-OPTIMALITY





 §.§ Lemmas



Let ε∼𝒩(0,1), and 
L_K,a∼ D_1 |D'D≤ a, where D = (D_1,…, D_K)∼𝒩(0, I_K). Both ε and L_K,a are symmetric and unimodal around zero. 



It follows from <cit.>.



Let 
ζ_0,ζ_1 and ζ_2 be three mutually independent random variables. 
If 


  (1) ζ_0 is symmetric and unimodal around zero; 

  (2) ζ_1 and ζ_2 are symmetric around 0;

  (3) P(|ζ_1| ≤ c) ≥ P(|ζ_2 | ≤ c) for any c ≥ 0; 

then P( | ζ_0+ζ_1 | ≤ c) ≥ P( | ζ_0+ζ_2 | ≤ c) for any c ≥ 0.



It follows from <cit.>.



Let ε∼𝒩(0,1), 
L_K,a∼ D_1 |D'D≤ a, where D = (D_1,…, D_K)∼𝒩(0, I_K), and ε and L_K,a be mutually independent. 
For any nonnegative constants b_1≤ c_1, b_2≤ c_2, and any α∈ (0,1), 
the 1-α quantile range of b_1ε + b_2 L_K,a is narrower than or equal to that of c_1ε + c_2 L_K,a. 



From Lemma <ref>, b_1ε is symmetric and unimodal. Because b_2 ≤ c_2,
P(| b_2L_K,a | ≤ c) ≥ P( |c_2L_K,a| ≤ c) for any c≥ 0. Then from Lemma <ref>, 
P( | b_1ε + b_2 L_K,a | ≤ c) ≥ P( | b_1ε + c_2L_K,a | ≤ c) for any c≥ 0. 


From Lemma <ref>, c_2L_K,a is symmetric and unimodal. Because b_1 ≤ c_1, P( | b_1ε | ≤ c) ≥ P( | c_1ε | ≤ c) for any c≥ 0. Then from Lemma <ref>, 
P( | b_1ε + c_2L_K,a | ≤ c) ≥ P( | c_1 ε + c_2L_K,a | ≤ c) for any c≥ 0. 


From the above two results, for any c≥ 0,

    P( | b_1ε + b_2 L_K,a | ≤ c) ≥ 
    P( | b_1ε + c_2L_K,a | ≤ c) ≥ 
    P( | c_1 ε + c_2L_K,a | ≤ c),

which implies Lemma <ref>.



	For any α≥ 1/2, the αth quantile of (1-ρ^2)^1/2·ε+ | ρ | · L_K,a is nonincreasing in ρ^2.  



It follows from <cit.>. 







 §.§ Proofs



In the asymptotic distribution of τ̂(β_1, β_0) in Corollary <ref>, both coefficients of ε and L_K,a attain their minimum values at r_1β_1 + r_1β_0≡γ  = γ̃ . From Lemma <ref> and Corollary <ref>,
the 𝒮-optimal adjusted estimator is attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 , 
with the asymptotic distribution (<ref>). 
From Proposition <ref>, τ̂(β̃_1, β̃_0) is 𝒮-optimal.







From Corollary <ref>, 
in the asymptotic distribution of τ̂(β_1, β_0) under ReM, 
the coefficient of ε does not depend on , and the coefficient of L_K,a attains its minimum when =. 
From Lemma <ref> and Corollary <ref>, 
the 𝒮-optimal adjusted estimator is 
attainable when = or 
r_0β_1+r_1β_0 = r_0 β̃_1+r_1 β̃_0 ,  with asymptotic distribution (<ref>). 
Moreover, from Proposition <ref>, 
τ̂(β̃_1, β̃_0) is 𝒮-optimal. 










From Corollary <ref>, the squared coefficient of ε in the asymptotic distribution of τ̂ is 
V_ττ(1-R^2_τ,x).
From Theorem <ref>, the squared coefficient of ε in the asymptotic distribution of τ̂(β̃_1, β̃_0) is 
V_ττ(β̃_1, β̃_0) {
1-R^2_τ, x(β̃_1, β̃_0)
}.
Because 

    V_ττ(β̃_1, β̃_0)
    = n{τ̂(β̃_1, β̃_0) 
    } 
    = 
    n(
    τ̂ - τ̂_w
    )
    = 
    n{(τ̂|τ̂_w)}
    =
    V_ττ( 1-R^2_τ, w),

the squared coefficient of ε in the asymptotic distribution of τ̂(β̃_1, β̃_0) under ReM reduces to 
V_ττ(1-R^2_τ,w){1 - R^2_τ,x(β̃_1, β̃_0)}.
Therefore, under ReM, the squared coefficient of ε in the asymptotic distribution of τ̂(β̃_1, β̃_0)  is smaller than or equal to that of τ̂ if and only if 

    V_ττ(1-R^2_τ,w){1 - R^2_τ,x(β̃_1, β̃_0)}≤ 
    V_ττ(1 - R^2_τ,x)
       ⟺   (1-R^2_τ,w)
    - (1-R^2_τ,w)R^2_τ,x(β̃_1, β̃_0)
    ≤(1 - R^2_τ,x)
       ⟺    
    R^2_τ,w 
    + (1-R^2_τ,w)R^2_τ,x(β̃_1, β̃_0)
    ≥ R^2_τ,x.








 §.§ The adjusted estimator with the smallest asymptotic variance

From Theorem <ref>, the asymptotic variance of τ̂(β_1, β_0) under ReM is 

    {V_ττ( 1-R^2_τ, x)
    ( 1-R^2_)
    + (r_1r_0)^-1( -_)'
    S^2_w∖x( -_)}
          + 
    {
    V_ττR^2_τ, x( 1 - R^2_) + (r_1r_0)^-1( -_)'
    S^2_w|x( -_)
    } v_K,a
    
     =     
    V_ττ( 1-R^2_τ, x)
    ( 1-R^2_) + 
    V_ττR^2_τ, x( 1 - R^2_) v_K,a
           + 
    (r_1r_0)^-1( -_)'
    S^2_w∖x( -_) + 
    (r_1r_0)^-1( -_)'
    S^2_w|x( -_) v_K,a .

It is a quadratic form of . 
The derivative of (<ref>) with respect to  is   

    2 (r_1r_0)^-1S^2_w∖x( -_) + 
    2v_K,a (r_1r_0)^-1S^2_w|x( -_)
        =     2 (r_1r_0)^-1{( 
    S^2_w∖x + v_K,aS^2_w|x) 
     - 
    (
    S^2_w∖x_ +
    v_K,aS^2_w|x_)
    }.

Therefore, under ReM, τ̂(β_1, β_0) with the smallest asymptotic variance is attainable when 

    r_0 β_1 + r_1 β_0≡
    =
    (
    S^2_w∖x  + v_K,aS^2_w|x)^-1(
    S^2_w∖x_ + 
    v_K,aS^2_w|x_).

When a ≈ 0, the above coefficient is close to _.
 
 
 
 

 §.§ Technical details for Remark <ref>


We first give equivalent forms for the squared coefficients of ε in the asymptotic distributions of τ̂ and τ̂(β̃_1, β̃_0) under ReM. 
		From Corollary <ref> and the definition of R^2_τ, x in (<ref>), the squared coefficient of ε in the asymptotic distribution of τ̂ under ReM has the following equivalent forms:
		
    V_ττ (1-R^2_τ, x)     = n ·(τ̂) · (1-R^2_τ, x) = n ·{(τ̂|τ̂_x)}.

		From Theorem <ref> and the definition of R^2_τ, x(β_1, β_0) in (<ref>), the squared coefficient of ε in the asymptotic distribution of τ̂(β̃_1, β̃_0) under ReM has the following equivalent forms:
		
    V_ττ(β̃_1, β̃_0) {1-R^2_τ, x(β̃_1, β̃_0)}
        = 
    		n ·{τ̂(β̃_1, β̃_0)}·{1-R^2_τ, x(β̃_1, β̃_0)}
    		= n ·{(
    		τ̂(β̃_1, β̃_0) |τ̂_x)
    		}
        = 
    		n ·{(
    		τ̂  - ' τ̂_w|τ̂_x)
    		}
        = 
    		n ·{(
    		τ̂|τ̂_x)
    		-
    		' 
    		(
    		τ̂_w|τ̂_x)
    		}.

		
We then study the covariance between (
		τ̂|τ̂_x) and (
		τ̂_w|τ̂_x) under the CRE. 
		For z=0,1, let Y^⊥_i(z) be the residual from the linear projection of Y_i(z) on x_i, and w^⊥_i be the residual from the linear projection of w_i on x_i. 
		We have 
		
    (
    		τ̂|τ̂_x)     = n_1^-1∑_i=1^n Z_i Y^⊥_i(1) - n_0^-1∑_i=1^n (1-Z_i) Y^⊥_i(0), 
    (
    		τ̂_w|τ̂_x)     = n_1^-1∑_i=1^n Z_i w_i^⊥ - n_0^-1∑_i=1^n (1-Z_i) w_i^⊥,

		and thus 
		
    {(
    		τ̂|τ̂_x), 
    		(
    		τ̂_w|τ̂_x)
    		}    = 
    		n_1^-1S_Y(1), w|x  + n_0^-1S_Y(0), w|x ,

		where S_Y(z), w|x is the finite population covariance between Y^⊥(z) and w^⊥, or, equivalently, the finite population partial covariance between Y(z) and w given x, for z=0,1. 
		
We finally prove Remark <ref>. 
		When S_Y(z), w|x = 0 for z=0,1, from (<ref>), 
		(
		τ̂|τ̂_x)
		is uncorrelated with 
		(
		τ̂_w|τ̂_x). 
		This further implies that 
		
    {(
    		τ̂|τ̂_x)
    		-
    		' 
    		(
    		τ̂_w|τ̂_x)
    		}    = 
    		{(
    		τ̂|τ̂_x)
    		}
    		+ 
    		{' 
    		(
    		τ̂_w|τ̂_x)
    		}≥{(
    		τ̂|τ̂_x)
    		}.

		From (<ref>) and (<ref>), 
		the coefficient of ε 
		for τ̂(β̃_1, β̃_0) is larger than or equal to that for τ̂. 







§ ESTIMATED DISTRIBUTIONS AND THE ASSOCIATED OPTIMALITY






 §.§ Lemmas

For treatment arm z (z=0,1), 
let s^2_Y(z) be the sample variance of observed outcome, 
s^2_w,z be the sample covariance of covariates w, and s_Y(z),w be the sample covariance between observed outcome and covariates w. 

 

Under ReM and Condition <ref>, for z=0,1, 

    s^2_Y(z) - S^2_Y(z) = o_P(1),   s^2_w,z - S^2_w = o_P(1),   s_Y(z),w - S_Y(z),w = o_P(1), 
    
       s^2_x,z - S^2_x = o_P(1),   s_Y(z),x - S_Y(z),x = o_P(1),   s_w,x - S_w,x = o_P(1);

and for any β_1 and β_0 that can depend implicitly on sample size n but have finite limits, 

    s^2_Y(z;β_z) - S^2_Y(z;β_z) = o_P(1),       s_Y(z;β_z)|x^2  - S_Y(z;β_z)|x^2 = o_P(1), 
    
    s_Y(z;β_z),x - S_Y(z;β_z),x = o_P(1),      s_Y(z;β_z),w - S_Y(z;β_z),w = o_P(1).





First, 
we can view covariates w as “outcomes” unaffected by the treatment. 
Thus, 
(<ref>)
follows immediately from <cit.>. 


Second, 
the observed sample variances have the following equivalent forms:

    s^2_Y(z;β_z)  = s^2_Y(z)  + β_z's^2_w,zβ_z - 2β_z' s_w, Y(z),      s_Y(z;β_z),x  = s_Y(z),x - β_z's_w,x,
    
    s_Y(z;β_z),w  = s_Y(z),w - β_z's_w,z^2,       
    s_Y(z;β_z)|x^2 = 
    s_Y(z;β_z), x( s_x,z^2 )^-1
    s_x, Y(z;β_z).

From 
(<ref>),

    s^2_Y(z;β_z) - S^2_Y(z;β_z) = o_P(1),   s_Y(z;β_z),x - S_Y(z;β_z),x = o_P(1),   s_Y(z;β_z),w - S_Y(z;β_z),w = o_P(1).

These further imply that 
s_Y(z;β_z)|x^2 - S_Y(z;β_z)|x^2 = o_P(1). 
Thus,  
(<ref>) holds. 





Under ReM and Condition <ref>, 

    V̂_ττ(β_1, β_0) - Ṽ_ττ(β_1, β_0) = o_P(1),  R̂^2_τ, x(β_1, β_0) 
    - 
    R̃^2_τ, x(β_1, β_0) = o_P(1),

where 

    Ṽ_ττ(β_1, β_0)≡ V_ττ(β_1, β_0)+S^2_τ∖w,  R̃^2_τ, x(β_1, β_0)  ≡Ṽ_ττ^-1(β_1, β_0)
    	V_ττ(β_1, β_0)
    	R^2_τ, x(β_1, β_0).
 



From (<ref>), (<ref>) and Lemma <ref>, 

    V̂_ττ(β_1, β_0)      =  
    r_1^-1S^2_Y(1;β_1) + r_0^-1S^2_Y(0;β_0) - 
    S^2_τ(β_1, β_0) |w+ o_P(1) 
       = V_ττ(β_1, β_0) + S^2_τ(β_1, β_0) ∖w + o_P(1)
        = V_ττ(β_1, β_0) + S^2_τ∖w + o_P(1)
     = Ṽ_ττ(β_1, β_0) + o_P(1), 
    
    
    
    R̂^2_τ, x(β_1, β_0)
        = 
    Ṽ_ττ^-1(β_1, β_0)^-1(
    r_1^-1S_Y(1;β_1)|x^2+r_0^-1S_Y(0;β_0)|x^2
     - 
    S^2_τ(β_1, β_0) |x) + o_P(1)
        = 
    Ṽ_ττ^-1(β_1, β_0)
    V_ττ(β_1, β_0) ·
    V_ττ^-1(β_1, β_0)(
    r_1^-1S_Y(1;β_1)|x^2+r_0^-1S_Y(0;β_0)|x^2
     - 
    S^2_τ(β_1, β_0) |x)+ o_P(1)
        = Ṽ_ττ^-1(β_1, β_0)
    V_ττ(β_1, β_0)
    R^2_τ, x(β_1, β_0)+ o_P(1) 
    = 
    R̃^2_τ, x(β_1, β_0) + o_P(1).








V_ττ(β_1, β_0) = 
V_ττ(1-R^2_τ, w)
+ 
(r_1r_0)^-1(-)'
S^2_w
	(-).



It follows from Theorem <ref> and Corollary <ref> with a=∞ and x = ∅. 
We give a more direct proof below:

    V_ττ(β_1, β_0)     = 
    n(
    τ̂ - τ - 'τ̂_w)
    = 
    n{τ̂ - τ - 'τ̂_w - 
    (
     -  )'τ̂_w}
        = 
    n{(τ̂|τ̂_w) - 
    (
     -  )'τ̂_w}
    = n{(τ̂|τ̂_w) 
    }
    + 
    n(
     -  )'(
    τ̂_w)(
     -  )
        = V_ττ(1-R^2_τ, w)
    + 
    (r_1r_0)^-1(-)'
    S^2_w
    	(-).





 §.§ Proofs













From Lemma <ref>, under ReM, the probability limit of the estimated distribution of τ̂(β_1, β_0) is 

    Ṽ_ττ^1/2(β_1, β_0)
    [
    {1-R̃^2_τ, x(β_1, β_0) }^1/2·ε + {R̃^2_τ, x(β_1, β_0) }^1/2·  L_K,a].

By definition, it has the following equivalent forms:

    {Ṽ_ττ(β_1, β_0)-Ṽ_ττ(β_1, β_0)R̃^2_τ, x(β_1, β_0) }^1/2·ε + 
    {Ṽ_ττ(β_1, β_0)R̃^2_τ, x(β_1, β_0)}^1/2·  L_K,a
    ∼   {V_ττ(β_1, β_0) + S^2_τ∖w-V_ττ(β_1, β_0)R^2_τ, x(β_1, β_0) }^1/2·ε + {V_ττ(β_1, β_0)R^2_τ, x(β_1, β_0)}^1/2·  L_K,a
    ∼   {V_ττ(β_1, β_0){
    1 - R^2_τ, x(β_1, β_0)
    } + S^2_τ∖w}^1/2·ε + {V_ττ(β_1, β_0)R^2_τ, x(β_1, β_0) }^1/2·  L_K,a.

Because Condition <ref> implies Condition <ref>,
from Corollary <ref> and its proof, we can further write the probability limit of the estimated distribution of τ̂(β_1, β_0) as

    {
    V_ττ(1-R^2_τ, w)
    +
    S^2_τ∖w
    + 
    (r_1r_0)^-1(-)'
    S^2_w∖x
    	(-)
    }^1/2·ε   
     + 
    { 
    (r_1r_0)^-1
    (-)'
    S^2_w|x
    (-)
    }^1/2· L_K,a.





It follows  from Theorem <ref> with a=∞ and x = ∅. 







With β_1 = β_0 = 0,
Theorem <ref> implies that the the probability limit of the estimated distribution of τ̂ is 







    Ṽ_ττ^1/2(
    1-R̃^2_τ, x)^1/2·ε + 
    Ṽ_ττ^1/2(R̃^2_τ, x)^1/2·  L_K,a
    ∼    ( Ṽ_ττ
    -Ṽ_ττR̃^2_τ, x)^1/2·ε + 
    ( Ṽ_ττR̃^2_τ, x)^1/2·  L_K,a
    ∼    ( V_ττ + S^2_τ∖w
    -V_ττ R^2_τ, x)^1/2·ε + 
    ( V_ττ R^2_τ, x)^1/2·  L_K,a
    ∼    {V_ττ(1-R^2_τ,x) + S^2_τ∖w}^1/2·ε + 
    (V_ττ R^2_τ,x)^1/2· L_K,a.








From  Lemma <ref>, the probability limit of the estimated distribution of τ̂(β_1, β_0) in (<ref>) is Ṽ_ττ^1/2(β_1, β_0) ·ε. 
From Lemma <ref>, this probability limit has the following equivalent forms:

    Ṽ_ττ^1/2(β_1, β_0) ·ε   ∼{  V_ττ(β_1, β_0) +  S^2_τ∖w}^1/2·ε
       ∼{ V_ττ(1-R^2_τ, w) +  S^2_τ∖w
    + 
    (r_1r_0)^-1(-)'
    S^2_w
    	(-)  }^1/2·ε.





 


From Theorem <ref>, in the probability limit of the estimated distribution of τ̂(β_1, β_0) in (<ref>), both coefficients of ε and L_K,a attain their minimum values at r_0β_1+r_1β_0 ≡ =. 
Lemma <ref> then implies that the optimal
adjusted estimator among (<ref>) in terms of the estimated precision is attainable when r_0β_1+r_1β_0 ≡ =. The corresponding probability limit of the estimated distribution is 
{
V_ττ(1-R^2_τ, w)
+
S^2_τ∖w}^1/2·ε.





It follows immediately from Theorem <ref>. 





 §.§ Additional comments on the asymptotic conservativeness under ReM


First, we consider the scenario under Condition <ref>. From Lemma <ref>, Corollary <ref> and Theorem <ref>, the probability limit of the estimated distribution of τ̂(β_1, β_0) has larger variance and wider quantile ranges than the asymptotic distribution of τ̂(β_1, β_0). Therefore, both the variance estimator and confidence intervals are asymptotically conservative. 

Second, we consider general scenario without Condition <ref>.
Theorem <ref> implies that the probability limit of the estimated distribution of τ̂(β_1, β_0) has larger variance and wider quantile ranges than V^1/2_ττ(β_1, β_0) ·ε.
Lemma <ref> implies that V^1/2_ττ(β_1, β_0) ·ε has larger variance and wider quantile ranges than
the asymptotic distribution of τ̂(β_1, β_0) in (<ref>). 
Therefore, both the variance estimator and the Wald-type confidence intervals are asymptotically conservative. 
 












§ GAINS FROM THE ANALYZER AND THE DESIGNER



	First, we compare the asymptotic variances. 
	From Corollary <ref>, 
	the asymptotic variance of τ̂ is 
	V_ττ{
	1 - (1-v_K,a)R^2_τ, x}. From Theorem <ref>, the asymptotic variance of τ̂(β̃_1, β̃_0) is
	V_ττ(1-R^2_τ,w).
	Compared to τ̂, the percentage reduction in the asymptotic variance of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    	1-R^2_τ,w/
    		1 - (1-v_K,a)R^2_τ, x
    	= 
    	
    		R^2_τ,w - (1-v_K,a)R^2_τ, x/
    		1 - (1-v_K,a)R^2_τ, x.

	
	
	Second, we compare the asymptotic quantile ranges. From Corollary <ref>, the length of the asymptotic 1-α quantile range of 
	τ̂ is 2V_ττ^1/2· q_1-α/2(R^2_τ,x). From Theorem <ref>, the length of the asymptotic 1-α quantile range of τ̂(β̃_1, β̃_0) is 2V_ττ^1/2(1-R^2_τ,w)^1/2· q_1-α/2(0). Compared to τ̂, the percentage reduction in the length of the asymptotic 1-α quantile range of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    		2V_ττ^1/2(1-R^2_τ,w)^1/2· q_1-α/2(0)
    	/
    		2V_ττ^1/2· q_1-α/2(R^2_τ,x)
    	
    	= 1 - (1-R^2_τ,w)^1/2·q_1-α/2(0)/q_1-α/2(R^2_τ,x).

	
	Third, because τ̂(β̃_1, β̃_0) is 𝒮-optimal, both percentage reductions in the variance and the 1-α quantile range are nonnegative.  
	It is easy to verify that they are both nondecreasing in R^2_τ,w. 




	First, we compare the asymptotic variances. 
	From Corollary <ref>,  
	the asymptotic variance of τ̂ is V_ττ{
	1 - (1-v_K,a)R^2_τ, x}.
	From Theorem <ref>, the asymptotic variance of τ̂(β̃_1, β̃_0) is V_ττ{
	1-(1-v_K,a)R^2_τ, x - v_K,aR^2_τ, w}.
	Compared to τ̂, the percentage reduction in the asymptotic variance of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    		V_ττ{
    		1-(1-v_K,a)R^2_τ, x - v_K,aR^2_τ, w}/
    		V_ττ{
    		1 - (1-v_K,a)R^2_τ, x}
    	= 
    	v_K,aR^2_τ, w/
    		1 - (1-v_K,a)R^2_τ, x.

	
	
	Second, we compare the asymptotic quantile ranges. 
	From Corollary <ref>, 
	the length of the asymptotic 1-α  quantile range of τ̂ is 2 V_ττ^1/2· q_1-α/2(R^2_τ,x). 
	From Theorem <ref>, 
	the length of the asymptotic 1-α  quantile range of τ̂(β̃_1, β̃_0) is 2 V_ττ^1/2 (1-R^2_τ, w)^1/2· q_1-α/2(ρ^2_τ, x∖w). 
	Compared to τ̂, 
	the percentage reduction in the asymptotic 1-α quantile range of τ̂(β̃_1, β̃_0) is 
	
    1 - ( 1-R^2_τ, w)^1/2· 
    		q_1-α/2(ρ^2_τ, x∖w)
    /
    		q_1-α/2(R^2_τ,x).

	
	Third, because τ̂(β̃_1, β̃_0) is 𝒮-optimal, both percentage reductions in variance and 1-α quantile range are nonnegative.  
	It is easy to verify that the percentage reduction in the variance is nondecreasing in R^2_τ,w. 
	For the quantile range, from Lemma <ref>, 
	(1-R^2_τ, w)^1/2· q_1-α/2(ρ^2_τ, x∖w), the (1-α/2)th quantile of (1-R^2_τ, x)^1/2·ε + 
	(R^2_τ, x-R^2_τ, w)^1/2·  L_K,a, is nonincreasing  in R^2_τ,w. 
	Hence the percentage reduction in the 1-α quantile range is nondecreasing in R^2_τ,w.
	





	First, we compare the asymptotic variances. 
	From Section <ref>, 
	the asymptotic variance of τ̂(β̃_1, β̃_0) under the CRE is 
	V_ττ(1-R^2_τ, w).
	From Theorem <ref>, the asymptotic variance of τ̂(β̃_1, β̃_0) under ReM is 
	V_ττ(1-R^2_τ, w){
	1 - (1-v_K,a)ρ^2_τ, x∖w}.
	Compared to the CRE, the percentage reduction in the asymptotic variance under ReM is 
	
    1 - 
    		V_ττ(1-R^2_τ, w){
    		1 - (1-v_K,a)ρ^2_τ, x∖w}/
    		V_ττ(1-R^2_τ, w)
    	
    	= (1-v_K,a)ρ^2_τ, x∖w.

	
	
	Second, we compare the asymptotic quantile ranges. 
	From Section <ref>, 
	the length of the asymptotic 1-α quantile range of τ̂(β̃_1, β̃_0) under the CRE is 
	2 V_ττ^1/2(1 - R^2_τ,  w)^1/2· q_1-α/2(0).
	From Theorem <ref>, 
	the length of the asymptotic 1-α quantile range of τ̂(β̃_1, β̃_0) under ReM is 
	2 V_ττ^1/2(1-R^2_τ, w)^1/2· q_1-α/2(ρ^2_τ, x∖w). 
	Compared to the CRE,  
	the percentage reduction in the length of the asymptotic 1-α quantile range under ReM is 
	1-
	q_1-α/2(ρ^2_τ, x∖w)/
	q_1-α/2(0). 
	
	Third, 
	from Lemma <ref>, 
	both percentage reductions in variance and 1-α quantile range are nonnegative and nondecreasing in ρ^2_τ, x∖w = (R^2_τ, x-R^2_τ, w)/(1 - R^2_τ, w). 
	Consequently, both percentage reductions are nondecreasing in R^2_τ, x. 



	Recall that κ=1+V^-1_ττS^2_τ∖w≥ 1. 
	From Corollary <ref>, under ReM and Conditions <ref> and <ref>, 
	the probability limit of the estimated distribution of τ̂ is 
	
    {V_ττ(1-R^2_τ,x) + S^2_τ∖w}^1/2·ε + {V_ττR^2_τ,x}^1/2· L_K,a
       ∼   
    		V_ττ^1/2{(κ-R^2_τ,x)^1/2·ε + ( R^2_τ,x)^1/2· L_K,a}
       ∼   κ^1/2 V_ττ^1/2{( 1-R^2_τ,x/κ)^1/2·ε + ( R^2_τ,x/κ)^1/2· L_K,a}.

	From Corollary <ref>, under ReM and Conditions <ref> and <ref>, 
	the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	
    {
    	V_ττ(1-R^2_τ, w)
    	+
    	S^2_τ∖w}^1/2·ε∼ 
    	V_ττ^1/2(
    	κ - R^2_τ, w)^1/2·ε∼κ^1/2 V_ττ^1/2(
    	1 - R^2_τ, w/κ)^1/2·ε.

	
	First, we compare the variances. 
	The variance of the probability limit of the estimated distributions of τ̂ is
	κ V_ττ{
	1 - (1-v_K,a)R^2_τ,x/κ}.
	The variance of the probability limit of the estimated distributions of  τ̂(β̃_1, β̃_0) is 
	κ V_ττ(1 - R^2_τ, w/κ).
	Compared to τ̂, 
	the percentage reduction in variance of the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    		1 - R^2_τ, w/κ/
    		1  - (1-v_K,a)R^2_τ,x/κ
    	= 
    		R^2_τ, w - (1-v_K,a)R^2_τ,x/κ  - (1-v_K,a)R^2_τ,x.

	
	
	Second, we compare the quantile ranges. 
	The length of 1-α quantile range of the probability limit of the estimated distribution of  τ̂ is 
	2κ^1/2 V_ττ^1/2· q_1-α/2(R^2_τ,x/κ). 
	The length of 1-α quantile range of the probability limit of the estimated distribution of  τ̂(β̃_1, β̃_0) is 
	2κ^1/2 V_ττ^1/2
	(
	1 - R^2_τ, w/κ
	)^1/2· q_1-α/2(0). 
	Compared to τ̂, 
	the percentage reduction in 1-α quantile range of the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    		q_1-α/2(0) ( 1 - R^2_τ, w/κ)^1/2/
    		q_1-α/2(R^2_τ,x/κ)
    	
    	= 
    	1 - ( 1 - R^2_τ, w/κ)^1/2·q_1-α/2(0)/q_1-α/2(R^2_τ,x/κ).


	
	Third, the 	optimality of τ̂(β̃_1, β̃_0) 
in terms of the estimated precision 
	implies that both percentage reductions are nonnegative. It is easy to verify that they are both nondecreasing in R^2_τ, w.




	From
	Theorem <ref>, 
	the probability limit of the estimated distribution of τ̂ is
	Ṽ_ττ^1/2·ε∼ 
 (V_ττ + S^2_τ∖w )^1/2·ε.
	From Corollary <ref>, 
	the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	{V_ττ(1-R^2_τ, w)
	+
	S^2_τ∖w}^1/2·ε.
	Compared to τ̂, the percentage reduction in variance 
	of the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	
    1 - V_ττ(1-R^2_τ, w)
    		+
    		S^2_τ∖w/V_ττ + S^2_τ∖w
    	= 1 - κ - R^2_τ,w/κ =     R^2_τ,w/κ ,

	and the percentage reduction in length of the 1-α  quantile range of the probability limit of the estimated distribution of τ̂(β̃_1, β̃_0) is 
	
    1 - 
    		2q_1-α/2(0){V_ττ(1-R^2_τ, w)
    		+
    		S^2_τ∖w}^1/2/2q_1-α/2(0)
    		( V_ττ + S^2_τ∖w)^1/2
    	= 
    	1 - ( κ - R^2_τ,w)^1/2/κ^1/2
    	= 1 - (1 - R^2_τ,w/κ)^1/2.
 
	It is easy to show that both percentage reductions are nonnegative and nondecreasing in R^2_τ,w. 



In the following two proofs, we recall that q_1-α/2(0) is the (1-α/2)th quantile of a standard Gaussian distribution, and use the fact that under either the CRE or ReM,
		the 1-α confidence interval covers the average treatment effect if and only if
		
    q_1-α/2(0) ≤V̂^-1/2_ττ( β̃_1, β̃_0)  ×  n^1/2{τ̂(β̃_1, β̃_0) - τ}≤ q_1-α/2(0).

Therefore, the limit of the probability that (<ref>) holds is the asymptotic coverage probability of the confidence interval. 
		


		From 
		Corollary <ref> and the comment after it,
	 the probability limits of the estimated distributions of τ̂(β̃_1, β̃_0) are the same under  the CRE and ReM, and so are the lengths of confidence intervals after being scaled by n^1/2. 










		From Lemma <ref>, 
		V̂_ττ( β̃_1, β̃_0) in (<ref>) has the same probability under the CRE and ReM.  
		From Theorem <ref> and Section <ref>, 
		n^1/2{τ̂(β̃_1, β̃_0) - τ} in (<ref>)
		converges weakly to the same distribution under the CRE and ReM. 
		From Slutsky's theorem, 
		the quantity in the middle of (<ref>) converges weakly to the same distribution under the CRE and ReM.
		Therefore, for any α∈ (0,1), the limit of the probability that (<ref>) holds is the same under the CRE and ReM, and so is
the asymptotic coverage probability of the 1-α confidence interval.





	From Theorem <ref> and Corollary <ref>, the probability limits of the estimated distributions of τ̂(β_1, β_0) are the same under  the CRE and ReM, and so are the lengths of confidence intervals after being scaled by n^1/2. 










	Using Lemma <ref>, Theorem <ref> and Slutsky's theorem, we have that under ReM, the quantity in the middle of (<ref>) is asymptotically equal to
	
    Ṽ^-1/2_ττ(β_1, β_0)
    	V_ττ^1/2(β_1, β_0)
    	[
    	{1-R^2_τ, x(β_1, β_0)}^1/2·ε + {R^2_τ, x(β_1, β_0)}^1/2·  L_K,a].

	Using Lemma  <ref>, Corollary <ref> and Slutsky's theorem, we have that under the CRE, the quantity in the middle of (<ref>) is asymptotically equal to
	
    Ṽ^-1/2_ττ(β_1, β_0) V_ττ^1/2(β_1, β_0) ·ε.

	From Lemma <ref>, the distribution (<ref>) has shorter quantile ranges than (<ref>). 
	Therefore, for any α∈ (0,1), the limit of the probability that (<ref>) holds under ReM is larger than or equal to that under the CRE. 






§ Τ̂(Β̂_1, Β̂_0) AND VARIANCE ESTIMATORS UNDER REM




We need additional notation.
Let U_i = (1, Z_i, w_i', Z_i w_i')' ∈ℝ^2J +2. 
Let w̅_1 and w̅_0 be the averages of covariates, and Y̅_1 and Y̅_0 be the averages of observed outcomes in treatment and control groups. 
We can verify that in the OLS fit of Y on U, the coefficient of Z is τ̂(β̂_1, β̂_0), and the residual for unit i is ê_i =  Y_i - β̂_1' w_i - (Y̅_1 -  β̂_1' w̅_1 ) for treated units with Z_i = 1 and ê_i =Y_i - β̂_0' w_i - (Y̅_0 -  β̂_0' w̅_0 ) for control units with Z_i = 0.







For z=0,1, let 
σ̂^2_e,z =n_z^-1∑_i: Z_i=zê_i^2
be the average of squared residuals, and 
m^2_w,z = n_z^-1∑_i:Z_i=zw_iw_i' be the second sample moment of w. 
Define

    G     =
    n^-1∑_i=1^n U_i U_i'
     = 
    [ G_11 G_12; G_21 G_22 ]
        = 
    n^-1∑_i=1^n 
    [            1          Z_i         w_i'     Z_i w_i';          Z_i          Z_i     Z_i w_i'     Z_i w_i';          w_i      Z_i w_i     w_i w_i' Z_i w_i w_i';      Z_i w_i      Z_i w_i  Z_i w_iw_i'  Z_i w_iw_i' ]
    =
    (
    [            1          r_1          w̅'    r_1 w̅_1';          r_1          r_1    r_1 w̅_1'    r_1 w̅_1'; [2pt/2pt]
    w̅     r_1 w̅_1        S^2_w  r_1 m^2_w,1;     r_1 w̅_1     r_1 w̅_1  r_1 m^2_w,1  r_1 m^2_w,1 ])  
    
    
    
    H    = 
    n^-1∑_i=1^n ê_i^2 U_i U_i'
    ≡[ H_11 H_12; H_21 H_22 ]
       =n^-1∑_i=1^n[               ê_i^2          Z_i  ê_i^2          ê_i^2 w_i'     Z_i  ê_i^2 w_i';          Z_i  ê_i^2          Z_i  ê_i^2     Z_i  ê_i^2 w_i'     Z_i  ê_i^2 w_i';           ê_i^2 w_i      Z_i  ê_i^2 w_i      ê_i^2 w_i w_i' Z_i  ê_i^2 w_i w_i';      Z_i  ê_i^2 w_i       Z_i ê_i^2 w_i   Z_i ê_i^2 w_iw_i'  Z_i  ê_i^2 w_iw_i' ]
        = 
    (
    [                          σ̂^2_e                    r_1 σ̂^2_e,1           n^-1∑_i=1^nê_i^2 w_i'     n^-1∑_i=1^n Z_i  ê_i^2 w_i';                    r_1 σ̂^2_e,1                    r_1 σ̂^2_e,1     n^-1∑_i=1^n Z_i  ê_i^2 w_i'     n^-1∑_i=1^n Z_i  ê_i^2 w_i';  [2pt/2pt]
    n^-1∑_i=1^nê_i^2 w_i       n^-1∑_i=1^nZ_i  ê_i^2 w_i       n^-1∑_i=1^nê_i^2 w_i w_i' n^-1∑_i=1^n Z_i  ê_i^2 w_i w_i';      n^-1∑_i=1^n Z_i  ê_i^2 w_i       n^-1∑_i=1^n Z_i ê_i^2 w_i   n^-1∑_i=1^n Z_i ê_i^2 w_iw_i'  n^-1∑_i=1^n Z_i  ê_i^2 w_iw_i' ]) .


The Huber–White variance estimator for n^1/2{τ̂(β̂_1, β̂_0) - τ} is 
V̂_HW  = [G^-1HG^-1]_(2,2),
the (2,2)th element of G^-1HG^-1. 






 §.§ Lemmas





Under ReM and Condition <ref>, 
τ̂_w = r_0^-1w̅_1 = - r_1^-1w̅_0 = O_P(n^-1/2).


For 1≤ j≤ J, define pseudo potential outcomes (Ỹ(1), Ỹ(0)) = (W_j, W_j). We can verify that Condition <ref> also holds if we replace the original potential outcomes by the pseudo ones. Corollary <ref> implies that τ̂_W_j= O_P(n^-1/2) and thus τ̂_w = O_P(n^-1/2). 



	Under ReM and Condition <ref>, for z=0,1, we have β̂_z - β̃_z = o_P(1), and
	
    s^2_Y(z;β̂_z) - S^2_Y(z;β̃_z) = o_P(1),       s_Y(z;β̂_z)|x^2 - S_Y(z;β̃_z)|x^2 = o_P(1)
    s_Y(z; β̂_z),x - S_Y(z; β̃_z),x = o_P(1),      s_Y(z; β̂_z),w - S_Y(z;β̃_z),w = o_P(1).




The results follow directly from 
Lemma <ref>.





 


For any two matrices A and B, if  both A and A+B are nonsingular, then 

    (A + B)^-1 - A^-1 = 
     A^-1B(A + B)^-1BA^-1 - 
    A^-1BA^-1.




Lemma <ref> is known, but we give a direct proof for completeness. 
From 

    0    = A - A =  A(A+B)^-1(A+B) - A = A(A+B)^-1A +  A(A+B)^-1B - A,
    
    0    = B - B =   (A+B) (A+B)^-1B  - B = A(A+B)^-1B + B(A+B)^-1B  - B,

we have 
A(A+B)^-1A  - A =  B(A+B)^-1B  - B,
which further implies 

    (A+B)^-1 - A^-1    = 
    A^-1{A(A+B)^-1A  - A}A^-1 = 
    A^-1{B(A+B)^-1B  - B}A^-1
        = A^-1B(A+B)^-1BA^-1
    -A^-1BA^-1.




Under ReM and Condition <ref>, m_w,z^2 = S_w^2 + o_P(1) for z=0,1. 
Both G_11 and G_22 in (<ref>) converge in probability to nonsingular matrices. 



First, we consider m_w,z^2. 
By definition, 

    m_w,1^2  = 
    n_1^-1∑_i:Z_i=1w_i w_i' = 
    n_1^-1∑_i:Z_i=1( w_i - w̅_1 ) ( w_i - w̅_1 )' + 
    w̅_1 w̅_1' 
     = n_1^-1 (n_1-1) s^2_w,1 + w̅_1 w̅_1'.

From Lemmas <ref> and <ref>, 
s^2_w,1 = S^2_w + o_P(1) and w̅_1 w̅_1'= o_P(1). 
Thus, m_w,1^2 = S^2_w + o_P(1). 
Similarly, m_w,0^2 = S^2_w + o_P(1). 

Second, we consider G_11. 
By definition, G_11 has a limit as n→∞. 
Because

    G_11^-1 = 
    [   1 r_1; r_1 r_1 ]^-1
    = 1/r_1 r_0[  r_1 -r_1; -r_1    1 ],
 
the limit of G_11 is nonsingular.

Third, we consider G_22. From the above and by definition,  

    G_22 = 
    [       S^2_w r_1 m^2_w,1; r_1 m^2_w,1 r_1 m^2_w,1 ]
    = 
    [     S^2_w r_1 S^2_w; r_1 S^2_w r_1 S^2_w ] + o_P(1).

Thus, G_22 has a probability limit as n→∞. 
Because the limit of S^2_w is nonsingular, 
and 

    [     S^2_w r_1 S^2_w; r_1 S^2_w r_1 S^2_w ]^-1
    =
    ( G_11⊗S^2_w )^-1
    = G_11^-1⊗ ( S^2_w )^-1,

the probability limit of G_22 is nonsingular.






Under ReM and Condition <ref>,  

	
  (i) for z=0,1, 
	σ̂^2_e,z  = s^2_Y(z;β̂_z) + o_P(1) = S^2_Y(z;β̃_z) + o_P(1) = O_P(1), 
	and 
	H_11 = O_P(1); 
	
  (ii)
	n^-1∑_i=1^n Z_i  ê_i^2 w_i' = o_P(n^1/2),
	n^-1∑_i=1^n (1-Z_i)  ê_i^2 w_i'  = o_P(n^1/2),
	and 
	H_12 = H_21' = o_P(n^1/2); 
	
  (iii)
	n^-1∑_i=1^n Z_i  ê_i^2 w_iw_i' = o_P(n),
	n^-1∑_i=1^n (1 - Z_i ) ê_i^2 w_iw_i' = o_P(n),
	and 
	H_22 = o_P(n).




	
	
	First, we prove (i).  By definition and from Lemma <ref>,  
	
    σ̂^2_e,z
    	=n_z^-1  (  n_z-1 )  s^2_Y(z;β̂_z) 
    	= 
    	s^2_Y(z;β̂_z) + o_P(1) = S^2_Y(z;β̃_z) + o_P(1) = O_P(1),   (z=0,1).

	This further implies σ̂^2_e = r_1 σ̂^2_e,1 + r_0 σ̂^2_e,0 = O_P(1). 
	Thus, H_11 = O_P(1). 
	
	Second, we prove (ii). 
	For any 1≤ j≤ J,  
	
    | 
    	n^-1∑_i=1^n Z_i  ê_i^2 w_ij|
    	   ≤max_1≤ i≤ n |w_ij | ·
    	n_1^-1∑_i=1^n Z_i  ê_i^2 = 
    	max_1≤ i≤ n |w_ij | ·σ̂^2_e,1.

	Condition <ref> implies that max_1≤ i≤ n |w_ij |/n^1/2→ 0 and thus max_1≤ i≤ n |w_ij | = o(n^1/2). 
	Lemma <ref>(i) implies σ̂^2_e,1= O_P(1). 
	Thus, n^-1∑_i=1^n Z_i  ê_i^2 w_ij  = o(n^1/2) O_P(1) = o_P(n^1/2). 
	This further implies n^-1∑_i=1^n Z_i  ê_i^2 w_i'  = o_P(n^1/2). 
	Similarly, n^-1∑_i=1^n (1-Z_i)  ê_i^2 w_i'  = o_P(n^1/2). 
	By definition, 
	
    n^-1∑_i=1^nê_i^2 w_i'  = n^-1∑_i=1^n Z_i ê_i^2 w_i' + n^-1∑_i=1^n (1-Z_i) ê_i^2 w_i' = o_P(n^1/2).

	Thus, H_12 = o_P(n^1/2).  
	
	Third, we prove (iii). 
	For any 1≤ l,j≤ J,  
	
    | n^-1∑_i=1^n Z_i  ê_i^2 w_il  w_ij|    ≤ 
    	n^-1∑_i=1^n Z_i  ê_i^2 |w_il | |w_ij |
    	≤max_1≤ i≤ n |w_il | ·max_1≤ i≤ n |w_ij | · 
    	n_1^-1∑_i=1^n Z_i  ê_i^2
        = max_1≤ i≤ n |w_il | ·max_1≤ i≤ n |w_ij | ·σ̂^2_e,1.

	Condition <ref> implies max_1≤ i≤ n |w_il | = o(n^1/2) and max_1≤ i≤ n |w_ij |  = o(n^1/2). 
	Lemma <ref>(i) implies σ̂^2_e,1=O_P(1). 
	Thus, 
	n^-1∑_i=1^n Z_i  ê_i^2 w_il  w_ij = o(n^1/2) o(n^1/2) O_P(1) = o_P(n). 
	This then implies n^-1∑_i=1^n Z_i  ê_i^2 w_iw_i' = o_P(n). 
	Similarly, n^-1∑_i=1^n (1 - Z_i ) ê_i^2 w_iw_i' = o_P(n).
	Thus
	
    n^-1∑_i=1^nê_i^2 w_i w_i' = 
    	n^-1∑_i=1^n Z_i ê_i^2 w_i w_i' + n^-1∑_i=1^n (1-Z_i) ê_i^2 w_i w_i' = o_P(n).

	By definition, H_22 = o_P(n). 
	
































































































 §.§ Proofs



We first consider the asymptotic distributions. Let = r_0 β̂_1 + r_1 β̂_0, which satisfies  -  =o_P(1) by Lemma <ref>. From Lemma  <ref>, 
τ̂_w = O_P(n^-1/2).
Then


    τ̂(β̂_1, β̂_0) - τ̂(β̃_1, β̃_0)
        = 
    ( τ̂ - ' τ̂_w) - ( τ̂ - ' τ̂_w) = (
     - )' τ̂_w 
    =o_P(1) O_P(n^-1/2)
    = o_P(n^-1/2),










which implies that n^1/2{τ̂(β̂_1, β̂_0) - τ} has the same asymptotic distribution as n^1/2{τ̂(β̃_1, β̃_0) -τ}.


We then consider the probability limit of the estimated distribution.
From Lemma <ref>, 
we can show that 
V̂_ττ(β̂_1, β̂_0) - Ṽ_ττ(β̃_1, β̃_0) = o_P(1) and R̂^2_τ, x(β̂_1, β̂_0) 
- 
R̃^2_τ, x(β̃_1, β̃_0) = o_P(1)
under ReM. 
Therefore, under ReM and Condition <ref>
the estimated distributions of τ̂(β̂_1, β̂_0) and τ̂(β̃_1, β̃_0) have the same probability limit.  




 
Let  

    Λ = 
    [    G_11 0_2× 2J; 0_2J× 2    G_22 ],   Δ = 
    G - Λ = 
    [   0_2 × 2      G_12;      G_21 0_2J × 2J ],   Ψ = G^-1 - Λ^-1  .


First, we first find the stochastic orders of Δ and Ψ. 
From Lemma <ref>, G_12 = G_21' = O_P(n^-1/2), and thus Δ = O_P(n^-1/2). 
From Lemma <ref>,  

    Ψ    = (
    Λ + Δ)^-1 - Λ^-1
    = 
    Λ^-1Δ( 
    Λ + Δ)^-1ΔΛ^-1
    - 
    Λ^-1ΔΛ^-1.

From Lemma <ref>, the probability limit of Λ exists and is nonsingular. 
Thus, Ψ = O_P(n^-1/2). 

Second, we consider the difference between G^-1HG^-1 and  Λ^-1HΛ^-1: 

    G^-1HG^-1 - Λ^-1HΛ^-1    = 
    (
    Λ^-1 + Ψ) H(
    Λ^-1 + Ψ) - Λ^-1HΛ^-1
        = 
    ΨHΛ^-1 + Λ^-1HΨ + 
    ΨHΨ.

In (<ref>), we focus on the sub-matrix of the first two rows and the first two columns. We consider the corresponding submatrices of the
three terms in (<ref>). Below let [ ·]_(1-2, 1-2) denote the sub-matrix of the first two rows and the first two columns, and [ ·]_(2,2) denote the (2,2)th element of a matrix.
The first term in (<ref>) is 

    ΨHΛ^-1 
     = 
    [ Ψ_11 Ψ_12; Ψ_21 Ψ_22 ][ H_11 H_12; H_21 H_22 ][ G_11^-1       0;       0 G_22^-1 ]
    = 
    
    
    
    
    
    
    
    
    
    [ Ψ_11H_11G_11^-1 + Ψ_12H_21G_11^-1                                 *;                                                                   * ].

From Lemmas 
<ref> and <ref>,

    [ΨHΛ^-1]_(1-2, 1-2)  
    = 
    Ψ_11H_11G_11^-1 + Ψ_12H_21G_11^-1
    
    = 
    O_P(n^-1/2) O_P(1) O_P(1) + O_P(n^-1/2) o_P(n^1/2) O_P(1) .

For the second term in (<ref>), 
[Λ^-1HΨ]_(1-2, 1-2) = [ΨHΛ^-1]_(1-2, 1-2)' = o_P(1). 
For the third term in (<ref>), 

Lemma <ref> implies 
H = o_P(n) and thus ΨHΨ = O_P(n^-1/2) o_P(n) O_P(n^-1/2)  = o_P(1). 
Therefore, from (<ref>), 
[G^-1HG^-1]_(1-2, 1-2) - [Λ^-1HΛ^-1]_(1-2, 1-2) = o_P(1),
which further implies 

    V̂_HW - [Λ^-1HΛ^-1]_(2,2)
    = 
    [G^-1HG^-1]_(2,2) - [Λ^-1HΛ^-1]_(2,2) = o_P(1) .


Third, we consider the difference between [Λ^-1HΛ^-1]_(2,2) and V̂_ττ(β̂_1, β̂_0). 
Because

    [Λ^-1HΛ^-1]_(1-2, 1-2)    = 
    G_11^-1H_11G_11^-1 
    = (r_1r_0)^-2[  r_1 -r_1; -r_1    1 ][       σ̂^2_e r_1 σ̂^2_e,1; r_1 σ̂^2_e,1 r_1 σ̂^2_e,1 ][  r_1 -r_1; -r_1    1 ]
        = 
    [                  r_0^-1σ̂^2_e,0                - r_0^-1σ̂^2_e,0;                - r_0^-1σ̂^2_e,0 r_1^-1σ̂^2_e,1 + r_0^-1σ̂^2_e,0 ].

using Lemma 

<ref>, we have 

    [Λ^-1HΛ^-1]_(2,2)    = 
    r_1^-1σ̂^2_e,1 + r_0^-1σ̂^2_e,0
    = 
    
    
    r_1^-1s^2_Y(1;β̂_1) + r_0^-1s^2_Y(0;β̂_0) + o_P(1).

The property of OLS implies s_Y(z;β̂_z),w = 0 for z=0,1. 
Using the definition in (<ref>), we can then simplify V̂_ττ(β̂_1, β̂_0) as 
V̂_ττ(β̂_1, β̂_0)  =  
r_1^-1s^2_Y(1;β̂_1) + r_0^-1s^2_Y(0;β̂_0).
Therefore, 
[Λ^-1HΛ^-1]_(2,2) = V̂_ττ(β̂_1, β̂_0) + o_P(1).
From the above, we have 
V̂_HW = [Λ^-1HΛ^-1]_(2,2) + o_P(1) = V̂_ττ(β̂_1, β̂_0) + o_P(1). 


