Consider a clinical trial for assessing the causal effects from giving k treatments. Patients in a sample for the trial are randomized into the k treatment arms according to fixed proportions. Let Y (t) be a discrete or continuous potential response from a patient under treatment t (we assume each patient in the trial receives one and only one treatment). Let θ be the k-dimensional vector whose tth component is θ t = E(Y (t) ), the unknown mean potential response under treatment t, where E denotes population expectation. Based on data collected from the trial, we would like to make statistical inference on given functions of θ, such as a linear contrast θ t −θ s , a ratio θ t /θ s , or an odds ratio {θ t /(1−θ t )}/{θ s /(1−θ s )} between two treatment arms t and s.
In clinical trials, we typically observe some baseline covariates for each patient, which are measured prior to treatment assignments and, hence, are not affected by the treatment.
As emphasized in regulatory agency guidelines, the baseline covariates can be utilized in the following two ways. (i) In the design stage, covariate-adaptive randomization can be used to enforce the balancedness of treatments across levels of Z whose components are discrete baseline prognostic factors, such as institution, disease stage, prior treatment, gender, and age group. "Balance of treatment groups with respect to one or more specific prognostic covariates can enhance the credibility of the results of the trial" (EMA, 2015) .
(ii) In the analysis stage, a vector X of baseline covariates can be used to gain efficiency.
"Incorporating prognostic factors in the primary statistical analysis of clinical trial data can result in a more efficient use of data to demonstrate and quantify the effects of treatment" (FDA, 2019) . More specifically, the investigator is advised to identify "the covariates expected to have an important influence on the primary outcome" and to specify "how to account for them in the analysis in order to improve precision and to compensate for any lack of balance between groups" (ICH E9 (R1), 2019).
The well known analysis of variance (ANOVA) provides asymptotically valid inference on functions of θ for either discrete or continuous Y (t) when simple randomization is used, i.e., patient treatment assignments are randomized independently. However, ANOVA does not utilize the baseline covariates and can be inefficient.
An alternative is to model the dependence of Y on a baseline covariate vector X for efficiency gain. However, such a model-based approach is unreliable because the inference may be invalid if the model is misspecified (EMA, 2015) . In practice, it is often hard to specify the correct model or test any lack of fit. For this reason, a great deal of effort has been made to develop model-assisted approaches that gain efficiency by utilizing covariates in a working model and produce asymptotically valid inference even when the working model is misspecified. Development of covariate adjustment methods is encouraged by the most recent guidance for COVID-19: "To improve the precision of treatment effect estimation and inference, sponsors should ... propose methods of covariate adjustment" (FDA, 2020) .
One example of the model-assisted approach is the analysis of covariance (ANCOVA) which is widely recognized by researchers and regulatory agencies (EMA, 2015; FDA, 2019 ).
The working model used by the customary ANCOVA, originally proposed by Fisher (1935), does not include treatment-by-covariate interaction terms ( §2.3) and, thus, it is referred to as the homogeneous working model. In fact, treatment-by-covariate interaction terms are often ignored or even discouraged because of two perceptions: (i) even if the working model is misspecified, ANCOVA can still provide valid inference as it is model-assisted; (ii) a model without interaction terms have fewer coefficients to estimate and may have better finite sample properties. However, these perceptions only provide a partial picture. When the treatment effect is indeed heterogeneous, the ANCOVA estimator of θ without including treatment-by-covariate interaction in its working model may be even less efficient than the ANOVA estimator that uses no model assistance at all (Freedman, 2008; Lin, 2013) .
The confusion about how ANCOVA should be implemented can be seen from recent con-flicting recommendations by regulatory agencies: "The primary model should not include treatment by covariate interactions." (EMA, 2015) ; "The prespecified primary model can include interaction terms if appropriate" and "The closer the model approximates the true relationship between the outcome and the covariates, the greater the improvement in the power of significance tests and the precision of estimates..." (FDA, 2019) .
This leads to the first of our three principles presented in this article:
Principle 1 (Guaranteed efficiency gain). The working model should be chosen so that the model-assisted estimator often gains but never loses efficiency when compared to a benchmark estimator that does not adjust for any covariate.
Throughout this paper, the benchmark estimator is ANOVA. To fulfill Principle 1, we recommend a heterogeneous working model for ANCOVA that includes treatment-bycovariate interaction terms. We call this approach ANalysis of HEterogeneous COVAriance (ANHECOVA) to distinguish from ANCOVA that specifically refers to ANCOVA using the homogeneous working model. Our recommendation is motivated by the theoretical results in earlier work (Cassel et al., 1976; Lin, 2013; Liu and Yang, 2020; Li and Ding, 2020, among others), which we extend in §3.2- §3.3 to covariate-adaptive randomization and multiple treatment arms. We show that the ANHECOVA estimator of θ is consistent, asymptotically normal, and asymptotically more efficient than the ANCOVA estimator or the benchmark ANOVA estimator; in fact, it is asymptotically the most efficient estimator within a class of model-assisted estimators. In §3.1 we offer explanations of why the heterogeneous working model is generally preferable over the homogeneous working model.
A concern for adding treatment-by-covariate interaction terms in a working model is that the model may contain too many parameters and, hence, the finite sample performance of the resulting asymptotically efficient inference procedure may be worse than that from a model without interaction terms. For example, EMA (2015) recommends a compact working model with "no more than a few covariates". To address this issue, in §3.1 we develop a variant of ANHECOVA that improves its finite sample performance by using a pooled estimator of the covariance matrix of baseline covariate used in working model.
Covariate-adaptive randomization has been widely adopted in modern clinical trials.
According to a recent review of nearly 300 clinical trials published in 2009 and 2014, 237 of them used covariate-adaptive randomization (Ciolino et al., 2019) . Unlike simple randomization, covariate-adaptive randomization generates a dependent sequence of treatment assignments. Hence, conventional inference procedures constructed based on simple randomization are not necessarily valid under covariate-adaptive randomization (EMA, 2015) .
However, conventional methods are often applied in practice even when covariate-adaptive randomization is used. Starting from Shao et al. (2010) , various methods have been proposed for valid inference, which in general requires the randomization scheme to have a well-understood property (Shao and Yu, 2013; Ye, 2018; Bugni et al., 2018) . However, two limitations persist: (i) it is yet unclear how to obtain valid inference under complicated randomization schemes such as Pocock-Simon's minimization (Taves, 1974; Pocock and Simon, 1975; EMA, 2015) ; (ii) the standard error formula for inference varies with the randomization scheme and often has to be derived and proven under each particular scheme, which is a nuisance to the practice. These motivate our second principle.
Principle 2 (Validity and universality). The model-assisted inference should be asymptotically valid and can be universally applied to all commonly used randomization schemes.
In addition to validity, the universality means that a unified procedure can be used under all commonly used randomization schemes, including covariate-adaptive randomization and simple randomization. Our asymptotic theory in §3.3-3.4 shows that ANHECOVA achieves this principle as its asymptotic variance is invariant to the randomization scheme, as long as the randomization scheme satisfies a very mild condition (C2) stated in §2.2. However, this appealing property is not shared by ANCOVA. First, the asymptotic normality of the ANCOVA estimator of θ requires an additional condition (C3) on randomization, which is not satisfied by the popular Pocock-Simon's minimization method. Second, the asymptotic variance of the ANCOVA estimator varies with the randomization scheme.
Our third principle concerns robust variance estimation for asymptotically normally distributed estimators of functions of θ, another crucial step for inference.
Principle 3 (Robust standard error). The model-assisted inference should use robust standard errors that account for heteroscedasticity.
The asymptotic theory for heteroscedasticity-robust standard errors was developed decades ago (Huber, 1967; White, 1980) and has been widely used in econometrics. Alarmingly, its usage in clinical trials is scarce. This is particularly relevant for model-assisted inference because heteroscedasticity may arise from incorrect working models. Furthermore, the robust standard error also needs to take into account of the centering of baseline covariates before fitting working models, which is needed to identify the mean potential responses. Thus standard formulas for heteroscedasticity-robust standard errors do not directly apply to model-assisted inference for clinical trials.
The conclusions we draw in this paper are based on a new theoretical result on the joint asymptotic distribution of a class of covariate-adjusted estimators for θ, under simple or covariate-adaptive randomization. Our new result is quite general and subsumes previous results (mostly on the linear contrast θ t − θ s ) as special cases (Shao et al., 2010; Lin, 2013; Bugni et al., 2018 Bugni et al., , 2019 . Furthermore, our theory offers new insights on when ANCOVA may gain or lose efficiency over ANOVA. For example, under simple randomization with two treatment arms, Lin (2013) showed that ANCOVA with homogeneous working model satisfies Principle 1 (guaranteed efficiency gain over ANOVA) if the treatment allocation is balanced. However, our theory shows that this does not extend to trials with more than two arms, and is thus a peculiar property.
After introducing the notation, basic assumptions, and working models in §2, we present the methodology and theory in §3. Some numerical results are given in §4. The paper is concluded with recommendations and discussions for clinical trial practice in §5. Technical proofs can be found in the supplementary material.
2 Trial Design and Working Models
Recall that Y (t) represents the potential response under treatment t, t = 1, . . . , k. We use Z to denote the vector of discrete baseline covariates used in covariate-adaptive randomization and X to denote the vector of baseline covariates used in model-assisted inference. The vectors Z and X are allowed to share the same entries.
Suppose that a random sample of n patients is obtained from the population under inves-
, Z, and X, respectively. We impose the following mild condition.
. . , n, are independent and identically distributed with finite second order moments. The distribution of baseline covariates is not affected by treatment and the covariance matrix Σ X = var(X i ) is positive definite. Notice that no model between the potential responses and baseline covariates is assumed in (C1), and the potential responses can be either discrete or continuous.
Let π 1 , . . . , π k be the pre-specified treatment assignment proportions, 0 < π t < 1, and k t=1 π t = 1. Let A i be the k-dimensional treatment indicator vector that equals a t if patient i receives treatment t, where a t denotes the k-dimensional vector whose tth component is 1 and other components are 0. For patient i, the treatment A i is assigned after baseline covariates Z i and X i are observed; only one treatment is assigned according to
if and only if A i = a t . Once the treatments are assigned and the responses are recorded, the statistical inference is based on the observed
The simple randomization scheme assigns patients to treatments completely at random,
i , X i )'s and are independent and identically distributed with P (A i = a t ) = π t , t = 1, ..., k. It does not make use of covariates and, hence, may yield sample sizes that substantially deviate from the target assignment proportions across levels of the prognostic factors.
To improve the credibility of the trial, it is often desirable to enforce the targeted treatment assignment proportions across levels of Z by using covariate-adaptive randomization.
The three most popular covariate-adaptive randomization schemes are the stratified permuted block (Zelen, 1974) , the stratified biased coin (Shao et al., 2010; Kuznetsova and Johnson, 2017) , both of which use all joint levels of Z as strata, and Pocock-Simon's minimization (Taves, 1974; Pocock and Simon, 1975; Han et al., 2009) , which aims to enforce treatment assignment proportions across marginal levels of Z.
All these covariate-adaptive randomization schemes, as well as the simple randomization, satisfy the following mild condition (Baldi Antognini and Zagoraiou, 2015).
(C2) The discrete covariate Z used in randomization has finitely many joint levels in Z and satisfies (i) given {Z i , i = 1, . . . , n}, {A i , i = 1, . . . , n} is conditionally independent
i , X i ), i = 1, . . . , n}; (ii) as n → ∞, n t (z)/n(z) → π t almost surely, where n(z) is the number of patients with Z = z and n t (z) is the number of patients with Z = z and treatment t, z ∈ Z, t = 1, ..., k.
The benchmark ANOVA approach does not model how the potential responses Y 
where ϑ is a k-dimensional unknown vector and c T denotes the row vector that is the transpose of a column vector c. By Lemma 2 in the supplementary material, ϑ identifies
is the mean potential response under treatment t. In the classical exact ANOVA inference, the responses are assumed to have normal distributions. So a common perception is that ANOVA can only be used for continuous responses. As normality is not necessary in the asymptotic theory, the ANOVA and the other approaches introduced next can be used for non-normal or even discrete responses when n is large.
To utilize the baseline covariates, ANCOVA is based on the following homogeneous working model,
where ϑ and β / are unknown vectors having the same dimensions as A and X, respectively, and µ X = E(X i ). There is no treatment-by-covariate interaction terms in (2), which is incorrect if patients with different covariates benefit differently from receiving the same treatment, a scenario that often occurs in clinical trials. By Lemma 2 in the supplementary
i ). Thus, the ANCOVA estimator with working model (2) is modelassisted (Theorems 1 and 3 in §3). Then, what is the impact of ignoring the treatmentby-covariate interaction effect when it actually exists? The impact is that the ANCOVA estimator may be inefficient; in fact, Freedman (2008) gave examples in which the ANCOVA estimator is even less efficient than the benchmark ANOVA estimator.
To better adjust for X, we consider an alternative working model for ANCOVA that includes the treatment-by-covariate interactions:
where ϑ, β / 1 , . . . , β / k are unknown vectors and I(·) is the indicator function. We call model ( 
i ), i.e., ANCOVA with working model (3) is also model-assisted.
To differentiate the methods based on (2) and (3), we refer to the method based on (2) as ANCOVA and the one based on (3) as ANHECOVA.
As a final remark, both working models (2) and (3) use the centered covariate vector X − µ X . Otherwise, ANCOVA and ANHECOVA do not directly provide estimators of θ.
Centering is crucial, the only non-trivial exception is when homogeneous working model (2) is used and linear contrast θ t −θ s is estimated, as the covariate mean µ X cancels out. When fitting the working models (2) and (3) with real datasets, we can use the least squares with µ X replaced byX, the sample mean of all X i 's. In other words, we can center the baseline covariates before fitting the models. Since this step introduces non-negligible variation to the estimation, it affects the asymptotic variance of model-assisted estimator of θ and its estimation for inference. Thus, we cannot assume the data has been centered in advance and µ X = 0 without loss of generality (see §3.4).
We first describe the estimators of θ under (1)-(3). The traditional ANOVA estimator iŝ
whereȲ t is the sample mean of the responses Y i 's from patients under treatment t. As n → ∞,θ AN is consistent and asymptotically normal.
Using the homogeneous working model (2), the ANCOVA estimator of θ is the least squares estimator of the coefficient vector ϑ in the linear model (2) with (A i , X i ) as regressors. It has the following explicit formula,
whereX t is the sample mean of X i 's from patients under treatment t,X is the sample mean of all X i 's, and
is the least squares estimator of β / in (2). It is shown in Theorems 1 and 3 thatθ ANC is consistent and asymptotically normal as n → ∞ regardless of whether working model (2) is correct or not, i.e., ANCOVA is model-assisted.
The termβ T (X t −X) in (5) is an adjustment for covariate X applied to the ANOVA estimatorȲ t . However, it may not be the best adjustment in order to reduce the variance.
A better choice is to use heterogeneous working model (3). The ANHECOVA estimator of θ is the least squares estimator of ϑ under model (3),
whereβ
is the least squares estimator of β / t in (3) for each t. It is shown in Theorems 1-3 below that the ANHECOVA estimatorθ ANHC is not only model-assisted, but also asymptotically at least as efficient asθ AN andθ ANC , regardless of whether model (3) is correct or not.
The following heuristics reveal why the adjustmentβ T t (X t −X) in (7) is better than the adjustmentβ T (X t −X) in (5), and why ANHECOVA often gains but never hurts efficiency even if model (3) is wrong. As the treatment has no effect on X, bothX t andX estimate the same quantity and, hence,β T t (X t −X) is an "estimator" of zero. As n → ∞, β t converges to β t = Σ −1 X cov(X, Y (t) ) in probability, regardless of whether (3) is correct or not (Lemma 3 in the supplementary material). Hence, we can "replace"
Consequently,Ȳ t −β T t (X t −X) has a smaller asymptotic variance thanȲ t . Note that (9) does not hold with β t replaced by other quantities. This explains why the adjustment β T (X t −X) in ANCOVA may lose efficiency, asβ in (6) converges to π 1 β 1 + · · · + π k β k .
The variance reduction technique by (9) can be found in the generalized regression (GREG) approach in the survey sampling literature (Cassel et al., 1976; Cochran, 1977; Särndal et al., 2003; Fuller, 2009; Shao and Wang, 2014; Ta et al., 2020) . From the theory of GREG,β t in (7) may be replaced by any estimator that converges to β t in probability, without affecting the asymptotic distribution of the GREG estimator. This motivates the following potential improvement to (8), which utilizes the fact that X has the same covariance across treatments and estimates the covariance matrix of X using all patients,
where n t is the number of units under treatment t. This alternative estimator alleviates the concern of using an unstable inverse in (8) when the sample size is small.
We consider asymptotic theory under simple randomization for a general class of estimators of the formθ
whereb t 's have the same dimension as X and can either be fixed or depend on the trial data. Note that class (11) contains all estimators we have discussed so far:
Theorem 1. Assume (C1) and simple randomization for treatment assignment.
Then, as n → ∞,
.., β k ), the matrix with columns β 1 , ..., β k , and B = (b 1 , ..., b k ). In particular, (13) holds forθ AN ,θ ANC , andθ ANHC as described by (12).
We briefly describe the proof for part (ii) in Theorem 1 and defer other details to the supplementary material. Notice that
The positive semidefiniteness of this matrix follows from the following algebraic result with
Lemma 1. Let M be a matrix whose columns are m 1 , ..., m k , and π 1 , ..., π k be nonnegative constants with k t=1 π t = 1.
We would like to emphasize that Theorem 1(i) holds regardless of whether model (3) is correct or not. Theorem 1(ii) shows that ANHECOVA not only satisfies Principle 1 (guaranteed efficiency gain over ANOVA), but is also the most efficient estimator within the class of estimators in (11) as it attains the optimal V SR (B). Another consequence of Theorem 1(ii) is that adjusting for more covariates in ANHECOVA does not lose and often gains asymptotic efficiency, although adjusting for fewer covariates may have better performance when n is small.
For the important scenario of estimating the linear contrast θ t − θ s with fixed t and s, the corresponding model-assisted estimator is c T tsθ , whereθ is given by (11) and c ts is the k-dimensional vector whose tth component is 1, sth component is −1, and other components are 0. The following corollary provides a direct comparison of the asymptotic variances of ANOVA, ANCOVA, and ANHECOVA estimators of linear contracts, showing that the ANHECOVA estimator has strictly smallest asymptotic variance except for some very special cases.
Corollary 1. Assume (C1) and simple randomization.
(i) The difference between the asymptotic variances of c T tsθAN and c T tsθANHC is
which is always ≥ 0 with equality holds if and only if
(ii) The difference between the asymptotic variances of c T tsθANC and c T tsθANHC is
which is always ≥ 0 with equality holds if and only if
When k = 2, i.e., there are only two arms, (14) reduces to π 2 β 1 + π 1 β 2 = 0; (15) reduces to β 1 = β 2 or π 1 = π 2 = 1/2. The same conclusions were also obtained by Lin (2013) under a different framework that only considers the randomness in the treatment assignments.
Liu and Yang (2020) extended the results in Lin (2013) to stratified simple randomization.
We would like to emphasize that when there are more than two arms (k > 2), (14) or (15) 
Therefore, β t = β s is not enough for ANCOVA to be as efficient as ANHECOVA for estimating θ t − θ s . Moreover, even if treatment allocation is balanced, i.e., π 1 = · · · = π k , ANCOVA is generally less efficient than ANHECOVA when there are more than two arms; this is different from the conclusion in the case of two arms. Finally, in estimating θ t − θ s for all pairs of t and s, for ANOVA to have the same asymptotic efficiency as ANHECOVA, all β t 's need to be zero, i.e., X is uncorrelated with Y (t) for every t; for ANCOVA to have the same asymptotic efficiency as ANHECOVA, all β t 's must be the same, i.e., models (2) and (3) are the same.
It is worth to mention that when there are more than two treatment arms, the ANCOVA estimator can be either more efficient or less efficient than the ANOVA estimator even under balanced treatment allocation. This is also observed by Freedman (2008) 
We now consider the estimation of θ under covariate-adaptive randomization as described in §2.2. Specifically, we would like to provide answers to the following questions: Is there an estimator whose asymptotic distribution is invariant with respect to different randomization schemes, to attain the validity and universality principle? Is there an estimator that is asymptotically most efficient within the class of estimators of form (11) under any covariateadaptive randomization?
The answers to these two questions are affirmative, as established formally in Theorems 2 and 3, respectively. Importantly, the key to achieve the validity and universality principle as well as guaranteed efficiency gain is using the ANHECOVAθ ANHC in (7) with all the joint levels of Z included in the covariate X.
Theorem 2. (Validity and Universality of ANHECOVA). Assume (C1) and (C2). If heterogeneous model (3) is used with X containing the dummy variables for all the joint levels of Z as a sub-vector, then, regardless of whether working model (3) is correct or not and which randomization scheme is used, as n → ∞,
Comparing Theorem 1 with Theorem 2, we see that the ANHECOVA estimator including all dummy variables for Z has exactly the same asymptotic variance in simple randomization and any covariate-adaptive randomization satisfying (C2), which is reflected by the fact that V SR (B) is the same as V SR (B) in (13) with B = B. Therefore, this estimator fulfills Principle 2 (validity and universality) so that the same inference procedure can be used regardless of which randomization scheme is used. As we show next, however, this is not true for ANOVA or ANCOVA using model (2), or for ANHECOVA when Z is not fully included in the working model.
To answer the second question, we need a further condition on the randomization scheme, mainly for estimators not using model (3) or not including all levels of Z in X.
(C3) There exist k × k matrices Ω(z), z ∈ Z, such that, as n → ∞,
where D is a block diagonal matrix whose blocks are matrices Ω(z)/P (Z i = z), z ∈ Z.
Condition (C3) weakens Assumption 4.1(c) in Bugni et al. (2019) , where Ω(z) takes a more special form. For simple randomization, Ω(z) = diag(π t ) − ππ T for all z, where π = (π 1 , . . . , π k ) T . For stratified permuted block randomization and stratified biased coin randomization, Ω(z) = 0 for all z. Note that Pocock-Simon's minimization scheme does not satisfy (C3) because the treatment assignments are correlated across strata, despite some recent theoretical result has been obtained (Hu and Zhang, 2020). Thus, the following result does not apply to minimization. However, as our Theorem 2 does not require (C3), it applies to minimization.
The next theorem establishes the asymptotic distribution of the class (11) of estimators under covariate-adaptive randomization, based on which we show the optimality of ANHECOVA.
Theorem 3. Assume (C1), (C2), and (C3). Consider class (11) of estimators and, without loss of generality, we assume that all levels of Z are included in X, as the components of b t 's in (11) corresponding to levels of Z not in X may be 0.
where
The main technical challenge in the proof of Theorem 2 and Theorem 3 is that the treatment assignments A 1 , . . . , A n are not independent, so we cannot directly apply the classical Linderberg central limit theorem. Instead, we decomposeθ(b 1 , ...,b k ) − θ into four terms and then apply a conditional version of the Linderberg central limit theorem to handle the dependence. The details can be found in the supplementary material.
A number of conclusions can be made from Theorem 3.
1. While Theorem 2 answers the first question in the beginning of §3.3, i.e.,θ ANHC with all joint levels of Z included in model (3) achieves Principle 2 (validity and universality), the second question is answered by Theorem 3(ii) showing thatθ ANHC also attains Principle 1 (guaranteed efficiency gain); in fact, it is asymptotically the most efficient estimator compared with all estimators in class (11). (3) or not including all levels of Z in (3) is that the asymptotic validity of the resulting estimator requires condition (C3), which is not needed in Theorem 2. Furthermore, the resulting estimator not only is less efficient according to the previous conclusion, but also has a more complicated asymptotic covariance matrix depending on the randomization schemes (universality is not satisfied), which requires extra handling in variance estimation for inference; see, for example, Shao et al. (2010) and Bugni et al. (2018) .
3. Under covariate-adaptive randomization satisfying (C2)-(C3), it is still true that the ANCOVA estimator using model (2) may be asymptotically more efficient or less efficient than the benchmark ANOVA estimator. (18), the asymptotic covariance matrix V (B) is invariant with respect to ran- (18) is 0, which is the case when B = B, i.e.,θ ANHC is used with all levels of Z included in X. If R(B) is not 0, such as the case of ANOVA, ANCOVA, or ANHECOVA not adjusting for all joint levels of Z, then V (B) depends on randomization scheme and, the smaller the Ω(z), the more efficient the estimator is. Thus, the stratified permuted block or biased coin with Ω(z) = 0 for all z is preferred in this regard.
domization scheme if R(B) in
where V SR (0) is the asymptotic variance of ANOVA estimator under simple randomization. As we vary the randomization scheme and the working model, the change in the asymptotic variance is determined by two terms. The first term {V SR (B)−V SR (0)} arises from using a working model; the second term E [R(B){Ω SR − Ω(Z i )}R(B)] is the reduction due to using a covariate-adaptive randomization scheme, and the reduction also depends on the working model being used via R(B).
The last conclusion from Theorem 3 requires a further derivation, which leads to the following corollary.
Corollary 2. Assume (C1)-(C3) and that X only includes the dummy variables for all joint levels of Z. Then, for any B in (17)
A direct consequence from Corollary 2 is that, if Ω(z) = 0 for all z (e.g., stratified permuted block or biased coin randomization is used) and X only includes all joint levels of Z, then all estimators in class (11), including the benchmark ANOVA estimator, have the same asymptotic efficiency as the ANHECOVA estimator under any randomization. In other words, modeling with all joint levels of Z is equivalent to designing with Z.
For inference on a function of θ based on Theorems 1 to 3, a crucial step is to construct a consistent estimator of asymptotic variance. For ANOVA or ANCOVA, the customary linear model-based variance estimation assuming homoscedasticity can be inconsistent, as criticized by Freedman (2008). This motivates Principle 3-robust standard error. That is, we should use variance estimators that are consistent regardless of whether the working model is correct or not and whether heterpscedaticity is present or not.
Consider the ANHECOVA estimatorθ ANHC in (7) (using either (8) or (10)), where the adjustment covariates X include all dummy variables for Z that is used in the randomization. There exist formulas for heteroscedaticity-robust standard error (such as those provided in the sandwich package in R). However, those formulas cannot be directly applied here, because they do not account for the additional variation introduced by centering the covariate X as required by the identification of θ. In fact, the term B T Σ X B in the asymptotic variance V SR (B) in Theorem 2 arises from centering X.
Instead, we should use the robust variance estimator based on V SR (B), as described next. LetΣ X be the sample covariance matrix of X i based on the entire sample and S 2 t (β t ) be the sample variance of (Y i −β T t X i )'s based on the patients in treatment arm t. Then V SR (B) in (16) can be estimated bŷ
This variance estimator is consistent as n → ∞ regardless of whether the heterogeneous working model (3) or homoscedaticity (i.e., var(Y (t) i − β t X i | X i ) is a constant) holds or not, and regardless of which randomization scheme is used.
In most applications the primary analysis is about treatment effects in terms of the linear contrast θ t − θ s = c T ts θ for one or several pairs of (t, s). For large n, an asymptotic
and z α is the (1 − α) quantile of the standard normal distribution. The same form of confidence interval can be used for any linear contrast c T θ (the sum of components of c is 0) with c T tsθANHC and SE 2 ts replaced by c Tθ ANHC and SE 2 c = c TV c, respectively. Let C be the collection of all linear contrasts with dimension k. An asymptotic level (1 − α) simultaneous confidence band of c T θ, c ∈ C , can be obtained by Scheffé's method,
where χ 2 α,k−1 is the (1 − α) quantile of the chi-square distribution with (k − 1) degrees of freedom. Correspondingly, to test the hypothesis that all θ t 's are the same, i.e., H 0 : θ 1 = · · · = θ k , an asymptotic level α chi-square test rejects H 0 if and only if
where C is the (k − 1) × k matrix whose tth row is c T tk , t = 1, ..., k − 1. Inference procedures based on the ANOVA or ANCOVA estimator can be similarly obtained using Theorems 1 and 3. However, as they do not satisfy the validity and universality principle, special treatment is needed for each covariate-adaptive randomization scheme. For example, under the stratified permuted block or biased coin randomization, the ANOVA or ANCOVA estimator is asymptotically more efficient than the same estimator under simple randomization; thus, using variance estimators valid only under simple randomization may lead to conservative inference. To eliminate the conservativeness, modifications depending on covariate-adaptive randomization schemes have to be made (Shao et al., 2010; Bugni et al., 2018) . This is another reason why we recommend ANHECOVA over the other model-assisted estimators for the practice.
We perform a simulation study based on the placebo arm of 481 patients in a real clinical trial to demonstrate the finite-sample properties of the model-assisted procedures. We use the observed continuous response of these 481 patients as the potential response Y (1) under treatment arm 1, and a 2-dimensional continuous baseline covariate (U, W ). The empirical distribution of (Y (1) , U, W ) of these patients is the population distribution in simulations.
Notice that we do not know the true relationship between Y (1) and (U, W ) because they are from the real data. Thus, the working models (2) and (3) may be misspecified.
In all numerical results, we apply (10) for ANHECOVA.
We consider three simulation studies that differ in how the potential responses Y (2) and Y (3) of the other two treatment arms are generated, and how the treatment assignment is randomized. Our first simulation compares the standard deviations of the ANOVA, ANCOVA, and ANHECOVA estimators of θ 2 − θ 1 , with X = U for ANCOVA and AN-HECOVA. The two additional potential responses are generated according to
and with sample size 481 (all data points) and treatment being assigned by simple randomization with three different allocation proportions: 1:2:2, 1:1:1, and 2:1:1. Thus, the only randomness in the first simulation study is from treatment assignments. Here µ U is the mean of 481 U -values and ζ ranges from −1 to 1 with ζ = 0 corresponds to β 1 = β 2 = β 3 .
Although we only consider the estimation of θ 2 − θ 1 , data from the third arm is still used by ANCOVA and ANHECOVA.
Based on 10,000 simulations, all three estimators have negligible biases and their standard deviations are plotted for different values of ζ in Figure 1 . The simulation result shows that, as predicted by our theory, ANHECOVA is generally more efficient than the other two estimators, except when ζ is nearly 0 where ANCOVA is comparable to ANHECOVA.
Furthermore, the simulation with the 1:2:2 allocation (left panel in Figure 1 ) shows very clearly that there is no definite ordering of the variances of ANCOVA and ANCOVA. Our
Corollary 1 suggests that a balanced allocation does not guarantee the superiority of AN-COVA over ANOVA when there are multiple arms (in contrast with the case of two arms), which can be seen from the simulations with 1:1:1 allocation (middle panel in Figure 1 ).
The second and third simulation studies are intended to examine the performance of standard errors and the coverage of proposed 95% asymptotic confidence intervals for θ 2 −θ 1 and θ 3 − θ 1 . In the second simulation study, treatment assignments are still generated by simple randomization (with allocation 1:1:1 or 1:2:2), but a random sample of size n = 200 or 400 is drawn from the 481 subjects with replacement. The setup (20) is also changed to
For both ANCOVA and ANHECOVA, we include U but not W in their working models, so X = U . The simulation results based on 10,000 simulations are shown in Table 1 .
In the third simulation study, the treatment assignments are generated by stratified permuted block randomization with block size 6 and treatment allocation 1:1:1 or block size 10 and treatment allocation 1:2:2. The covariate Z used in randomization is a threelevel discretized W . Similar to the second simulation study, a random sample of size n = 200 or 400 is drawn, but the setup (21) is changed to
where µ W is the mean of all 481 W -values. For ANCOVA and ANHECOVA, we consider two working models with different choices of X. One model includes the dummy variables for Z and U but not W , motivated by the fact that Z is a discretization of W and carries some information from W . The other model includes the dummy variables for Z, U , and W (this is different from just including U and W because W is continuous and the working model is linear). The simulation results based on 10,000 simulations are shown in Table 2 .
Notice that in the third simulation that uses a covariate-adaptive randomization scheme, consistent standard errors for ANOVA or ANCOVA can be obtained using Theorem 3 but they are not readily available. For this reason, in the third simulation we use the standard error derived under simple randomization based on Theorem 1. In other words, standard errors for each estimator are computed using the same formula in the second and third simulation studies. Based on our theory, it is expected that confidence intervals based on ANOVA and ANCOVA are conservative.
The following is a summary of the simulation results in Tables 1-2. 1. All estimators have negligible bias compared to their standard deviation.
2. ANHECOVA has the smallest standard deviation in all our simulation settings. This is consistent with our asymptotic theory.
3. There is no unambiguous ordering of the standard deviations of ANCOVA and AN-COVA. In particular, in the third simulation study (Table 2) , ANCOVA adjusting for Z, U, W is better in estimating θ 2 − θ 1 but worse in estimating θ 3 − θ 1 . A similar observation can also be obtained from 
We further illustrate the different model-assisted procedures using a real data example. Chong et al. (2016) conducted a randomized experiment to evaluate the impact of low dietary iron intake on human capital attainment. They recruited students of age 11 to 19 in a rural area of Cajamarca, Peru, where many adolescents suffer from iron deficiency. The goal of this randomized trial is to quantify the causal effect of reduced adolescent anemia on school attainment. By using students' school grade as covariate Z with five levels, a stratified permuted block randomization with 1:1:1 allocation was applied to assign 219 students to one of the following three promotional videos:
Video 1: A popular soccer player is encouraging iron supplements to maximize energy;
Video 2: A physician is encouraging iron supplements for overall health;
Video 3: A dentist encouraging oral hygiene without mentioning iron at all. Chong et al. (2016) investigated whether showing different promotional videos to the students can improve their academic performance through increased iron intake. Video 3 is treated as a "placebo". After the treatment assignments, four students were excluded from the analysis for various reasons, which we also ignore in our analysis. The dataset is available at https://www.openicpsr.org/openicpsr/project/113624/version/V1/view. Chong et al. (2016) used various outcomes in their analysis; here we focus on one of their primary outcomes-the academic achievement-as an example. In this trial, the academic achievement is measured by a standardized average of the student's grades in math, foreign language, social sciences, science, and communications in a semester. For the model-assisted approaches, we use the baseline anemia status as the covariate in our working models (2) and (3), which is believed to moderate the treatment effect (Chong et al., 2016) . 
To improve its credibility and efficiency, we believe a clinical trial analysis can benefit from Principles 1 to 3 outlined in §1 and discussed throughout §2-3.
Our theoretical investigation shows that the only method that satisfies all the three principles is the ANHECOVA with all joint levels of Z included in heterogeneous working model (3), coupled with the robust variance estimator given by (19). Thus we believe it deserves wider usage in the clinical trial practice. Our theory shows that using any other estimator with the form (11), including ANOVA, ANCOVA using model (2), and ANHECOVA not adjusting for all joint levels of Z, suffers from invalidity, inefficiency, or non-universality in the sense that the asymptotic distribution of the estimator depends on a particular randomization scheme.
Our model-assisted asymptotic theory is given in terms of the joint asymptotic distribution in estimating θ, the vector of mean responses, with multiple treatment arms randomized with a wide range of covariate-adaptive randomization schemes. It can be readily used for inference about linear or nonlinear functions of θ, with either continuous or discrete responses. Although working models (2) and (3) In applications, one may be interested in comparing just two treatments, although the trial contains more than two treatment arms. A simple way of analysis is to ignore the data from other arms and apply inference procedures to the two arms of interest. For ANOVA, this is equivalent to using all the arms, since ANOVA does not borrow strength from other arms through using covariates. However, using data from two arms only is not recommended when ANHECOVA is applied, because it can utilize covariate data from other arms to gain efficiency. Regarding ANCOVA, there is no definite order of efficiency for using the whole dataset or data from two given arms, since using more covariate data in ANCOVA may increase or decrease efficiency.
Our Theorem 2 can also be applied to rerandomization schemes (Li et al., 2018; Li and Ding, 2020) using discrete covariates. Rerandomization also tries to balance the treatment assignments across levels of Z, but it randomizes the treatment assignments for all patients simultaneously. This is different from the covariate-adaptive randomization schemes that assign treatment sequentially. For two-armed trials, (Li et al., 2018 , Corollaries 1 and 2) have shown that rerandomization satisfies (C2), but it is yet unknown if it satisfies (C3). Similar results for model-assisted inference can also be found in (Li and Ding, 2020, Theorem 3).
As a final cautionary note, standard software does not produce asymptotically valid standard errors for model-assisted inference. We implement an R package called RobinCar to compute the model-assisted estimators and their robust standard errors, that is available at https://github.com/tye27/RobinCar.
The supplementary material contains all technical proofs. Table 2 : Bias, standard deviation (SD), average standard error (SE), and coverage probability (CP) of 95% asymptotic confidence interval under stratified permuted block randomization (block size = 6) and setup (22) based on 10,000 simulations
Allocation n Method (C1) and (C2) be given and suppose P (A i = a t | Z 1 , . . . , Z n ) = π t for all t = 1, . . . , k and i = 1, . . . , n.
(i) For any integrable function f , we have that
(ii) Let θ = (E(Y (1) , ..., E(Y (k) )) be the potential response mean vector, β = k t=1 π t β t ,
The condition P (A i = a t | Z 1 , . . . , Z n ) = π t for all t and i holds for most covariateadaptive randomization schemes. Note that it does not exclude the possibility that the set of random variables {A i , i = 1, . . . , n} is dependent on {Z i , i = 1, . . . , n}, which is indeed the case for covariate-adaptive randomization schemes. We impose this condition only in Lemma 2 to facilitate understanding the working models. This additional assumption is not needed for our asymptotic theory in Section 3, as condition (C2) (ii) is sufficient. 
i , X i , A i ) are independent and identically distributed. For covariate-adaptive randomization, we remark that the property of conditional independence (Dawid, 1979, Lemma 4 .3), (C2) (i) and
the third condition in Lemma 2 imply that
In the equalities above, we have used, respectively, the consistency of potential responses, the law of iterated expectation, (C2) (i), and the remark above.
(ii) We only prove the first result. The second result can be proved similarly. Notice that (θ, β) satisfies the following estimation equations:
for any t (S1)
From Lemma 2, (S1) implies that for any t, Proof. We prove the result forβ t . The proofs forβ are analogous and omitted. Notice that
becomes an average of independent random variables once conditional on {A, F}. From the existence of second moment of XY (t) , and the weak law of large numbers for independent random variables, we conclude that, for any > 0,
From the bounded convergence theorem, the above equation also holds unconditionally. In other words,
This together with the fact that n t /n = z n t (z)/{ z n(z)} = π t + o p (1), we have
Similarly, we can show the result with X i Y i replaced by X i or Y i also holds, i.e.,
The denominator ofβ t can be treated similarly, which leads to
The proof is completed by using the definition of β t .
2 Technical Proofs 2.1 Proof of (9)
Under simple randomization, A = (A 1 , . . . , A n ) are independent with other variables. We
where the last equality is from cov(X i , X j | A) = cov(X i , X j ) = 0 for i = j. Similarly, we can show thatȲ t andX −t are uncorrelated conditional on A.
Then,
where the second equality is from cov(X −t ,Ȳ t | A) = 0, E(Ȳ t | A) = E(Y (t) ) and the identity that cov(X,
where the second equality uses the identity that var(X) = E{var(X | Z)}+var{E(X | Z)},
For any fixed k-dimensional vector = ( 1 , . . . , k ) , we have
where tr denotes the trace of a matrix, Q denotes a p-dimensional random vector that takes value π −1 t t m t with probability π t , t = 1, . . . , k, and the last equality follows from the fact that the covariance matrix var(Q) = E(QQ ) − E(Q)E(Q ) is positive semidefinite.
and we analyze the two terms separately.
Consider M 1 , where the tth term equals 
which also holds unconditionally from the bounded convergence theorem, i.e.,
Then, consider M 2 , which can be reformulated as
From the central limit theorem, we can easily derive that
Finally, consider cov(M 1 , M 2 ), where the (t, s) element equals Combining the above results, we conclude that √ n{θ(b 1 , . . . , b k ) − θ} is asymptotically normal with mean 0 and variance V SR (B),
Thus, we have
The rest follows from applying Lemma 1 with M = Σ 1/2 X (B − B).
From Lemma 3, we know thatβ = β + o p (1) andβ t = β t + o p (1), t = 1, . . . , k. Let σ 2 A , σ 2 B , σ 2 U respectively be the asymptotic variance of c tsθANHC , c tsθANC and c tsθAN , where from Theorem 1,
The results in Corollary 1 (i) follows from
where the second equality follows from β t = Σ −1 X cov(X, Y (t) ). This also proves that σ 2 A ≤ σ 2 U , because Σ X is positive definite and π t + π s ≤ 1. If σ 2 A = σ 2 U , then we must have π s β t + π t β s = 0 and (1 − π t − π s ){β t − β s } = 0.
To show the results in Corollary 1 (ii), notice that
where the second equality is because
In order to show that σ 2 A − σ 2 B ≤ 0, we prove a stronger statement: it is true that for anyβ,
As a consequence, settingβ as β = k t=1 π t β t , the statement in (S3) also holds. This proves
In what follows, we prove the claim in (S3). Note that the gradient of the left hand side of (S3) is
which equals zero whenβ = {π s β t + π t β s }/(π t + π s ). This is also the unique solution from the positive definiteness of Σ X . It is also easy to see that the Hessian of the left hand side of (S3) is negative definite, which means thatβ = {π s β t + π t β s }/(π t + π s ) is the global and unique maximizer of the left hand side of (S3). The statement in (S3) is true because when evaluated atβ = {π s β t + π t β s }/(π t + π s ), the left hand side of (S3) equals
This completes the proof for σ 2 A ≤ σ 2 B , where the equality holds if and only if {β t − β s }(1 − π t − π s ) = 0 and k t=1 π t β t = {π s β t + π t β s }/(π t + π s ).
We analyze the class of estimators of form (11) under covariate-adaptive randomization.
First, fromX t −X = O p (n −1/2 ) andb t = b t + o p (1), we havê θ(b 1 , . . . ,b k ) =θ(b 1 , . . . , b k ) + {(X 1 −X)(b 1 −b 1 ), . . . , (X 1 −X)(b k −b k )} =θ(b 1 , . . . , b k ) + o p (n −1/2 ) Then, we still decomposeθ(b 1 , . . . , b k ) aŝ
For M 1 , its tth term equals
+o p (n −1/2 ) where the last equality is because n t /n = π t + o p (1) and the remaining terms are shown in the following to be O p (n −1/2 ).
Thus, the decomposition can be further refined aŝ
+o p (n −1/2 ) First, consider C 1 . Let A = {A 1 , . . . , A n } and F = {Z 1 , . . . , Z n }. Conditioned on A, F, the terms in C 1t , the tth element in C 1 , are independently distributed. From (C2) (i), we have E(C 1t | A, F) = 0. From Lindeberg's Central Limit Theorem, as n → ∞,
From the bounded convergence theorem, this result also holds unconditionally, i.e., C 1t var(C 1t | A, F) d − → N (0, 1).
Note that nvar(C 1t | A, F) = 1 nπ 2 and thus,
Therefore, when X contains all joint levels of Z and b t = β t , √ n(C 2 + C 3 ) = 0 a.s..
Consider M 2 , because M 2 does not involve treatment indicators, so M 2 follows the same asymptotic distribution as the M 2 in the proof of Theorem 1, with b t = β t , i.e.,
Finally, consider cov(C 1 , M 2 ), whose (t, s) element equals
where the first equality is from E(C 1t | A, F) = 0, the fourth equality is because cov(Y (t)
i − X i β t , X j β s | F) = 0 for i = j, the fifth equality is because E(Y (t) i −X i β t | Z i ) = θ t −µ X β t and, thus, cov{E(Y (t) i − X i β t | Z i ), E(X i β s | Z i )} = 0. The sixth equality is because
where the last equality is because n t (z)/n = P (Z i = z)π t + o p (1) and that n t (z)/n is bounded, hence, E{n t (z)}/n = P (Z i = z)π t + o(1).
In summary, we have shown that the asymptotic distribution of √ n{θ(b 1 , . . . , b k ) − θ} is invariant under randomization schemes satisfying (C2), and satisfies √ n{θ(b 1 , . . . , b k ) − θ} = √ n(C 1 + M 2 ) + o p (1) where π = (π 1 , . . . , π n ) . Therefore, we have that √ n(C 1 + C 2 + C 3 ) is asymptotically normal with mean zero and variance
As proved in Theorem 2, because M 2 does not involve treatment indicators, it is still
In the following, we deal with the covariance between C 1 , C 2 , C 3 and M 2 . Let M 2s be the sth element in M 2 , the covariance between C 1t and M 2s can be derived in a similar fashion as in (S6) and, thus,
Then, we show that (C 2 , M 2 ) converges in distribution to two independent random variables. Denote C 2 | F d − → ξ C , M 2 d − → ξ M , where ξ C is as defined in (S7). Note that P (C 2 ≤ t, M 2 ≤ s)
= E{P (C 2 ≤ t, M 2 ≤ s | F)} = E{P (C 2 ≤ t | F)P (M 2 ≤ s | F)} using (C2) (i)
where the last step follows from the bounded convergence theorem.
Moreover, cov{C 3t , M 2s } = cov{C 3t , (X − µ X ) b s } = 1 n 2