Inference after covariate-adaptive randomisation: aspects of methodology and theory

Abstract

Covariate-adaptive randomisation has a more than 45 years of history of applications in clinical trials, in order to balance treatment assignments across prognostic factors that may have influence on the outcomes of interest. However, almost no theory had been developed for covariate-adaptive randomisation until a paper on the theory of testing hypotheses published in 2010. In this article, we review aspects of methodology and theory developed in the last decade for statistical inference under covariate-adaptive randomisation. We focus on issues such as whether a conventional procedure valid under the assumption that treatments are assigned completely at random is still valid or conservative when the actual randomisation is covariate-adaptive, how a valid inference procedure can be obtained by modifying a conventional method or directly constructed by stratifying the covariates used in randomisation, whether inference procedures have different properties when covariate-adaptive randomisation schemes have different degrees of balancing assignments, and how to further adjust covariates in the inference procedures to gain more efficiency. Recommendations are made during the review and further research problems are discussed.

Keywords:

• balancedness of assignments
• efficiency
• model-assisted approach
• model free inference
• stratification
• survival analysis

1. Introduction

In a clinical trial to compare $k\ge 2$ treatments, patients are typically randomised into treatment arms according to fixed treatment assignment proportions ${\pi }_1,\dots ,{\pi }_k$, where each ${\pi }_t$ is a known number strictly between 0 and 1 and $\sum _{t=1}^k{\pi }_t=1$. The simplest randomisation scheme assigns patients to treatments completely at random and, thus, is called complete randomisation or simple randomisation. However, simple randomisation may yield imbalance assignments, i.e. sample sizes not following the assignment proportions across some prognostic factors or covariates, e.g., institution, disease stage, prior treatment, gender and age, which are thought to have significant influence on the outcomes or responses of interest. For instance, a trial exhibiting a substantial imbalance in patient age or disease stage between two treatment arms may not pass a regulatory review even though a statistically significant treatment effect has been shown. The issue is more serious when patients are not all available for simultaneous assignment of treatments but rather arrive sequentially and must be treated immediately.

This leads to the development of covariate-adaptive randomisation (which is also referred to as dynamic allocation), i.e., treatment assignment of the ith patient is made dependent on the observed covariate value of this patient and the assignments and covariate values of all i−1 previously assigned patients. It should be emphasised that covariate-adaptive randomisation does not use any outcomes or responses from the i−1 previous patients when the ith patient is randomised to a treatment arm. Adaptive randomisation methods using outcomes or responses are not our focus and can be found, for example, in Hu and Rosenberger (2006), Zhang et al. (2007), Hu et al. (2009), Rosenberger and Lachin (2015), and the references therein. The oldest covariate-adaptive randomisation scheme is the minimisation proposed by Taves (1974) and its extensions in Pocock and Simon (1975). Other popular covariate-adaptive randomisation methods include the stratified permuted block randomisation method (Zelen, 1974), the stratified urn design (Wei, 1977; Zhao & Ramakrishnan, 2016) and the stratified biased coin method (Kuznetsova & Johnson, 2017; Shao et al., 2010). Summaries of different allocation schemes are given by Kalish and Begg (1985), Schulz and Grimes (2002), and Rosenberger and Sverdlov (2008).

How often are covariate-adaptive schemes applied in clinical trials? According to Taves (2010), from 1989 to 2008, over 500 clinical trials implemented the minimisation method to balance important covariates. In a recent review of nearly 300 clinical trials published in year 2009 and year 2014 (Ciolino et al., 2019), 237 of them used covariate-adaptive randomisation. Other examples can be found in van der Ploeg et al. (2010), Fakhry et al. (2015), Breugom et al. (2015), Stott et al. (2017), and Sun et al. (2018). In the 2018 New England Journal of Medicine, there are seven articles about covariate-adaptive randomisation (Horn et al., 2018; Jourdain et al., 2018; McKeever et al., 2018; Mehra et al., 2018; Myles et al., 2018; Ramirez et al., 2018; Zannad et al., 2018). Applications of covariate-adaptive randomisation are not limited to clinical trials, as they are relevant for randomised experiments with many interventions.

How should inference be carried out with data collected under covariate-adaptive randomisation? Unfortunately, tests and other inference procedures constructed based on simple randomisation, which will be called conventional tests and inference procedures, are often applied in practice after data are collected under covariate-adaptive randomisation. For example, in the seven articles cited previously in the 2018 New England Journal of Medicine, they all used conventional tests for treatment effect. On one hand, over 35 years between 1974 and 2009, there were many empirical results showing that some conventional tests could still control Type I errors in spite of using covariate-adaptive randomisation; see, for example Birkett (1985), Forsythe (1987), Aickin (2002), Weir and Lees (2003), Hagino et al. (2004), and Zhong and Kim (2008). On the other hand, the Committee for Proprietary Medicinal Product commended that ‘it remains controversial whether the analysis adequately reflects the randomisation scheme’ (Committee for proprietary medicinal products, 2004) and the European Medicines Agency 2015 guidelines stated that ‘possible implications of dynamic allocation methods [ minimisation] on the analysis, e.g., with regard to bias and Type I error control should be carefully considered,… conventional statistical methods do not always control the Type I error’ (EMA, 2015). Because a statistical inference procedure on treatment effects should be valid under the particular randomisation scheme used in data collection, the application of conventional inference procedures after covariate-adaptive randomisation has definitely raised concerns and controversies.

Why don't we always apply an inference procedure valid under a given covariate-adaptive randomisation scheme? In their review of covariate-adaptive randomisation, Rosenberger and Sverdlov (2008) stated:

Very little theoretical work has been done in this area, despite the proliferation of papers. The original source papers are fairly uninformative about theoretical properties of the procedures.

That is, the lack of theoretical work in developing valid inference procedures associated with covariate-adaptive randomisation schemes is probably the main reason why conventional procedures are applied in applications. Why is there so little theoretical work in this problem prior to 2008? Unlike simple randomisation, covariate-adaptive randomisation generates some dependence among treatment assignments, covariates and outcomes under which asymptotic distributions of treatment effect estimators (such as the difference of sample averages) could not be easily derived. Shao et al. (2010) initiated theoretical studies on the validity of statistical tests under covariate-adaptive randomisation. The following three issues are addressed in their paper:

1. Can we develop a test procedure valid under covariate-adaptive randomisation?

2. If we use covariate-adaptive randomisation and a conventional test procedure valid under simple randomisation, will the Type I error of the test be inflated?

3. Is a test under covariate-adaptive randomisation more powerful than it is under simple randomisation?

If we have affirmative answers to Questions (A)–(C), or at least Questions (A)–(B), then the concerns and controversies about using covariate-adaptive randomisation will be largely eliminated. As the first piece of theoretical work, the results in Shao et al. (2010) are limited to certain types of tests, randomisation schemes, and models between covariates and responses. Fortunately, significant progresses in the theory of this area have been made in the last decade, e.g., Hu and Hu (2012), Shao and Yu (2013), Ma et al. (2015), Bugni et al. (2018, 2019), Ye (2018), Ma et al. (2020), Ye and Shao (2020), and Ye et al. (2020). Another stream of results is based on permutation or re-randomisation methods, e.g., Simon and Simon (2011), Kaiser (2012), and Bugni et al. (2018).

The purpose of this article is to review aspects of methodology and theory for statistical inference after covariate-adaptive randomisation. We concentrate on Questions (A)–(C) previously stated and the main results in the last decade, some of which are very recent. It is our hope that this review will provide some guidance for clinical trialists about which valid inference procedures to use for various situations, and will shed light on further research and development in this important area.

2. Covariates, outcomes and treatment effects

First, let's describe covariates and outcomes or responses under a clinical trial. Consider a clinical trial with a total of n patients that are assigned to $k\ge 2$ treatment arms denoted by $a=1,\dots ,k$. From patient $i\in \{1,\dots ,n\}$, let ${\mathit{X}}_i$ be the vector of all observed covariates and let $Y_i^{\left(a\right)}$ be the potential outcome or response of interest under treatment assignment $a=1,\dots ,k$. $Y_i^{\left(a\right)}$ is called potential outcome because only one of $Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)}$ will be observed from patient i, as each patient receives only one treatment. Thus what we observe from patient i is $Y_i=Y_i^{\left(a\right)}$ if treatment a is assigned to patient i.

The outcome or response $Y_i^{\left(a\right)}$ could be continuous or discrete, or a survival time. In a survival trial, censoring is typically involved so that, for patient i, $Y_i^{\left(a\right)}=min(T_i^{\left(a\right)},C_i^{\left(a\right)})$ together with an indicator of the event $T_i^{\left(a\right)}\le C_i^{\left(a\right)}$ are observed, where $T_i^{\left(a\right)}$ is the potential survival or failure time and $C_i^{\left(a\right)}$ is the potential censoring time, under treatment $a=1,\dots ,k$.

Throughout, we assume the following minimal condition on covariates and outcomes.

(C1)	$(Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{X}}_i)$, $i=1,\dots ,n$, are independent and identically distributed.

Note that there is no assumption on the relationship between the covariates and potential outcomes. We allow arbitrary treatment effect heterogeneity, i.e., the effect of treatment and covariate interaction on potential outcomes.

For convenience, we use $Y^{\left(1\right)},\dots ,Y^{\left(k\right)},\mathit{X}$ to denote the variables from a generic patient. Under (C1), $(Y^{\left(1\right)},\dots ,Y^{\left(k\right)},\mathit{X})\sim (Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{X}}_i)$ for every i, where $\mathcal{X}\sim \mathcal{Y}$ means that $\mathcal{X}$ has the same distribution as $\mathcal{Y}$.

To assess treatment effect, we may be interested in the average treatment effect between any fixed pair of treatment arms, a and b, defined as $E(Y^{\left(a\right)}-Y^{\left(b\right)})$, where E is the population expectation and $E\left(Y^{\left(a\right)}\right)$ is assumed to be well defined. Another important measure in comparing treatments a and b is the quantile treatment effect defined as $q_{\tau }^{\left(a\right)}-q_{\tau }^{\left(b\right)}$ (Firpo, 2007; Zhang et al., 2020), where $q_{\tau }^{\left(a\right)}$ is the τth quantile of the distribution of $Y^{\left(a\right)}$ under treatment a and τ is a fixed fraction. Quantile treatment effect is more appropriate when potential outcomes are highly skewed and is more relevant and informative than the average treatment effect when some distributional impacts have to be assessed.

Both average treatment effect and quantile treatment effect are some characteristics of the distributions of potential outcomes. In some applications, we would like to assess the treatment effect on the entire distribution of a potential outcome or the entire conditional distribution of a potential outcome given covariates. For example, in a survival analysis we may be interested in testing whether the conditional distributions of $T^{\left(a\right)}$ given $\mathit{X}$ are the same for different a's.

Here we would like to make it clear that treatments may have effect not only on the marginal distributions of potential outcomes, but also the conditional distributions of $Y^{\left(a\right)}$ given $\mathit{X}$, although typically treatments may not have any effect on the marginal distribution of the covariate $\mathit{X}$.

So far we have not yet discussed the treatment assignment of patients. Suppose that treatment assignments are made according to some probabilistic mechanism. For patient i, let $A_i$ be the treatment assignment indicator vector, i.e., $A_i=e_a$ if patient i is assigned to treatment a, where $e_a$ is a vector whose ath component is 1 and rest components are 0's, $a=1,\dots ,k$. The observed outcome from patient i is $Y_i=Y_i^{\left(a\right)}$ if and only if $A_i=e_a$, $a=1,\dots ,k$, $i=1,\dots ,n$.

3. Covariate-adaptive randomisation schemes

We now introduce details about how $A_i$'s are generated according to a randomisation scheme, using or without using covariates ${\mathit{X}}_i$'s.

Under simple randomisation, $A_i$'s are independent of $(Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{X}}_i)$'s and, further, $A_i$'s are independent and identically distributed with $P(A_i=e_a)={\pi }_a$, where P denotes the probability under a given randomness mechanism. It should be emphasised that the independence between $A_i$'s and $(Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{X}}_i)$'s means that the treatment assignments are independent of potential outcomes and covariates, not that treatments have no effect on potential outcomes or conditional distribution of $Y^{\left(a\right)}$ given $\mathit{X}$ as discussed in Section 2.

Let $\mathit{Z}$ be a vector of discrete covariates with finitely many levels to be utilised in covariate-adaptive randomisation. Typically, components of $\mathit{Z}$ are some discrete components of $\mathit{X}$ and/or some discretised continuous components of $\mathit{X}$ that are thought to have significant influence on the potential outcomes. In the following, we describe some popular covariate-adaptive randomisation schemes for enforcing assignment allocation across at levels of $\mathit{Z}$. In a typical covariate-adaptive randomisation scheme, for the ith patient arrived with observed ${\mathit{Z}}_i$, the treatment assignment indicator $A_i$ is generated depending on not only the value of ${\mathit{Z}}_i$ but also the $\mathit{Z}$-values and assignments of the previous i−1 patients, $i=1,\dots ,n$.

The stratified permuted block randomisation method (Zelen, 1974) randomly assigns a block of size B patients into k arms each having $B{\pi }_a$ patients for every B sequentially arrived patients with $\mathit{Z}=\mathit{z}$, a particular level of $\mathit{Z}$, where B is appropriately chosen so that $B{\pi }_a$'s are integers and the last block is allowed to be incomplete. This method is called stratified permuted block randomisation since randomisation is carried out within each stratum (joint level of $\mathit{Z}$) to achieve balancedness of assignments across strata.

When k = 2 and ${\pi }_1={\pi }_2=1/2$, the stratified biased coin method (Shao et al., 2010) assigns patient i with ${\mathit{Z}}_i=\mathit{z}$ according to the biased coin randomisation in Efron (1971), $P(A_i=e_1)=\left\{\begin{array}{ll}p,& D_{i-1}\left(\mathit{z}\right)<0,\\ 1/2,& D_{i-1}\left(\mathit{z}\right)=0,\\ 1-p,& D_{i-1}\left(\mathit{z}\right)>0,\end{array}\right.$where p is a fixed constant satisfying $1/2<p<1$ and $D_{i-1}\left(\mathit{z}\right)$ is one half of the within $\mathit{z}$ stratum difference between the numbers of patients in treatment 1 and treatment 2 after i−1 assignments have been made. An extension of the stratified biased coin for general case of $k\ge 3$ can be found in Kuznetsova and Johnson (2017).

The stratified urn design (Wei, 1977, 1978a, 1978b) is the stratified biased coin randomisation with p depending on i. When k = 2 and ${\pi }_1={\pi }_2=1/2$, the fixed p in biased coin is replaced by a $p_i$ depending on $D_{i-1}\left(\mathit{z}\right)$. According to Wei (1977), the urn design would force balance at the beginning of treatment allocation, and approach simple randomisation as the size of trial increases. A stratified urn design for general situation of $k\ge 3$ can be constructed using the method described in Zhao and Ramakrishnan (2016).

The previous three stratified covariate-adaptive randomisation schemes enforce balancedness of treatment assignment allocation across all strata, i.e., joint levels of $\mathit{Z}$. However, the oldest covariate-adaptive randomisation scheme, the minimisation, is very different from these three methods.

First, consider k = 2 and ${\pi }_1={\pi }_2=1/2$. For each $i\in \{1,\dots ,n\}$, let $G_i^{\left(1\right)}$ be a weighted sum of squared or absolute differences between numbers of patients in two treatment arms over marginal levels of $\mathit{Z}$, where the calculation is based on i−1 previously assigned patients and the assumption that the ith patient i is assigned to treatment 1, and let $G_i^{\left(2\right)}$ be the same sum except that the ith patient is assumed to be in treatment 2. For a = 1 or 2, $G_i^{\left(a\right)}$ represents the ‘total amount of imbalance’ in treatment numbers across the marginal levels of $\mathit{Z}$ which exists if treatment a is assigned to the ith patient. Therefore, we would like to assign the ith patient by minimising $G_i^{\left(a\right)}$ over a = 1, 2, i.e., we assign the ith patient to treatment 1 if $G_i^{\left(1\right)}<G_i^{\left(2\right)}$, to 2 if $G_i^{\left(1\right)}>G_i^{\left(2\right)}$, and to 1 or 2 randomly if $G_i^{\left(1\right)}=G_i^{\left(2\right)}$. This is why the method is called the minimisation by Taves (1974). Pocock and Simon (1975) extended the minimisation by allowing minimisation with a given probability, i.e., $P(A_i=e_1)=\left\{\begin{array}{ll}p,& G_i^{\left(1\right)}<G_i^{\left(2\right)},\\ 1/2,& G_i^{\left(1\right)}=G_i^{\left(2\right)},\\ 1-p,& G_i^{\left(1\right)}>G_i^{\left(2\right)},\end{array}\right.$where $p>1/2$ is a fixed constant. Pocock and Simon's method is still referred to as the minimisation and Taves' minimisation is the special case with p = 1. For a general k and/or allocation, the minimisation can be similarly constructed (Han et al., 2009; Pocock & Simon, 1975).

If $\mathit{Z}$ is one dimensional, then the minimisation is the same as the stratified biased coin method. For a multivariate $\mathit{Z}$, the key distinction between the minimisation and the three previously described stratified randomisation methods is that enforcing treatment balancedness is at all joint levels of $\mathit{Z}$ for the latter but only at marginal levels of $\mathit{Z}$ for the former. For this reason, the minimisation is also called the marginal method in Ma et al. (2015) and Ye and Shao (2020). Enforcing treatment balance in marginal levels of $\mathit{Z}$ is sufficient in most applications.

Any of the previously introduced covariate-adaptive randomisation schemes satisfy

(D1)	$\{A_i,i=1,\dots ,n\}$ and $\{Y_i^{\left(1\right)},\dots ,Y_i^{\left(k\right)},{\mathit{X}}_i,i=1,\dots ,n\}$ are conditionally independent given $\{{\mathit{Z}}_i,i=1,\dots ,n\}$.

Actually, (D1) almost always holds for covariate-adaptive randomisation, because treatments, not their assignments, may affect the potential responses as we discussed earlier for simple randomisation, and given ${\mathit{Z}}_i$'s, the rest of ${\mathit{X}}_i$'s contain covariate information not used in randomisation.

Furthermore, all covariate-adaptive randomisation schemes considered so far satisfy the following condition (D2) (Baldi Antognini & Zagoraiou, 2015). In what follows, ⇒ denotes convergence in distribution as the sample size $n\to \mathrm{\infty }$, and $\Rightarrow 0$ is in fact convergence to 0 in probability.

(D2)

For every $i=1,\dots ,n$, $P(A_i=e_a|{\mathit{Z}}_1,\dots ,{\mathit{Z}}_n)={\pi }_a$ and, for every a and every level $\mathit{z}$ of $\mathit{Z}$, $\left\{n\right(\mathit{z}\left){\}}^{-1}{D}^{\left(a\right)}\right(\mathit{z})$ $\Rightarrow 0$, where $D^{\left(a\right)}\left(\mathit{z}\right)=n_a\left(\mathit{z}\right)-{\pi }_{a}n\left(\mathit{z}\right)$, $n\left(\mathit{z}\right)$ is the number of patients with ${\mathit{Z}}_i=\mathit{z}$, and $n_a\left(\mathit{z}\right)$ is the number of patients with ${\mathit{Z}}_i=\mathit{z}$ under treatment a.

Note that $D^{\left(a\right)}\left(\mathit{z}\right)$ is a measure of the assignment imbalance in stratum $\mathit{z}$. According to the asymptotic property of $D^{\left(a\right)}\left(\mathit{z}\right)$ in (D2), covariate-adaptive randomisation schemes can be classified into one of the following three types.

Type 1.	For every a and every $\mathit{z}$ of $\mathit{Z}$, $\left\{n\right(\mathit{z}\left){\}}^{-1/2}{D}^{\left(a\right)}\right(\mathit{z})\Rightarrow 0$.
Type 2.	For every a, $D^{\left(a\right)}\left(\mathit{z}\right)$'s with all different strata $\mathit{z}$'s are mutually independent and, for every $\mathit{z}$, $\left\{n\right(\mathit{z}\left){\}}^{-1/2}{D}^{\left(a\right)}\right(\mathit{z})\Rightarrow N(0,v_a)$, the normal distribution with mean 0 and a known variance $v_a>0$.
Type 3.	Methods not in Type 1 or 2.

The three types are in the order of the degree in enforcing the balancedness within every $\mathit{z}$ using the assignment imbalance measure $D^{\left(a\right)}\left(\mathit{z}\right)$. Type 1 is the strongest, requiring $D^{\left(a\right)}\left(\mathit{z}\right)$ diverging slower than the square root of within stratum $\mathit{z}$ sample size. Representatives of Type 1 covariate-adaptive randomisation methods are stratified permuted block and biased coin schemes. In fact, under stratified permuted block randomisation, $D^{\left(a\right)}\left(\mathit{z}\right)$ is bounded; for the stratified biased coin method, it follows from a result in Efron (1971) that $D^{\left(a\right)}\left(\mathit{z}\right)$ is bounded in probability for every $\mathit{z}$.

Type 2 is weaker than Type 1, as $\left\{n\right(\mathit{z}\left){\}}^{-1/2}{D}^{\left(a\right)}\right(\mathit{z})$ converges in distribution to $N(0,v_a)$, not 0. The stratified urn design is Type 2 with $v_a=1/12$ when k = 2, ${\pi }_1={\pi }_2=1/2$ (Wei, 1978a, 1978b). Simple randomisation treated as a special case of covariate-adaptive randomisation is also Type 2. Finally, the minimisation is Type 3, since it is neither Type 1 nor Type 2 (Ye & Shao, 2020). Specifically, under minimisation, $D^{\left(a\right)}\left(\mathit{z}\right)$ and $D^{\left(a\right)}\left({\mathit{z}}^{\prime }\right)$ with $\mathit{z}\ne {\mathit{z}}^{\prime }$ are not independent, and their relationship is complicated, because assignments are made according to marginal levels of $\mathit{Z}$.

4. Validity and conservativeness of tests

Testing a null hypothesis of no treatment effect on potential outcomes is the most utilised statistical inference procedure in clinical trials. For a given null hypothesis $H_0$ and a significance level $\alpha >0$, a test statistic $T$ is a function of observed $\{Y_i,{\mathit{X}}_i,i=1,\dots ,n\}$, which is constructed such that $H_0$ is rejected if and only if T is outside of the interval $[z_{\alpha /2},z_{1-\alpha /2}]$, where $z_r$ is the rth quantile of a known distribution, usually the standard normal distribution, in which case $H_0$ is rejected if and only if $|T|>z_{1-\alpha /2}$ as $z_{1-\alpha /2}=-z_{\alpha /2}$. Here, we consider two sided tests; the discussion for a one sided test is similar and omitted. T is said to be asymptotically valid (or valid for simplicity) if (1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) T is said to be asymptotically conservative (or conservative for simplicity) if (2) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)<\alpha .$(2)

4.1. Validity of conventional tests

As we discussed in Section 1, prior to 2010, there was almost no theoretical work and practitioners applied conventional tests developed under simple randomisation, which caused concerns about whether Type I error could be inflated. That is, if a conventional test T is applied after covariate-adaptive randomisation, does (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) still hold?

Forsythe (1987) concluded that a conventional test T still controls Type I error when $\mathit{Z}$ used in minimisation is also included in the construction of T. However, this conclusion was based on simulation results under certain models.

The first piece of theoretical work in this area obtained by Shao et al. (2010) is that, under covariate-adaptive randomisation, a conventional T is valid in the sense of (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) if both of the following hold:

1. The covariate $\mathit{Z}$ used in covariate-adaptive randomisation is a function of all covariates used to construct the test T.

2. T is valid in the sense of (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) under any fixed set of treatment allocation $A_1,\dots ,A_n$.

Note that (i) coincides with Forsythe's simulation discovery. But (ii) requires the validity of T under any deterministic allocation $A_1,\dots ,A_n$, which can be realistically achieved only when a correct statistical model is used in constructing T. However, correctly impose a model is difficult. Although mathematically, (i)–(ii) is only sufficient not necessary for the validity of a conventional test T under covariate-adaptive randomisation, we can easily find an example in which T is not valid under covariate-adaptive randomisation when either (i) or (ii) fails; e.g., Shao et al. (2010) and Shao and Yu (2013).

4.2. Conservativeness of conventional tests

Before we answer Question (A) in Section 1 regarding the development of a valid test according to (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) under covariate-adaptive randomisation, we would like to address Question (B) in Section 1, i.e., whether or not a conventional test T is conservative in the sense of (2(2) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)<\alpha .$(2) ). If the answer is yes, then at least the Type I error is not inflated by using conventional tests.

The first result of this kind was obtained by Shao et al. (2010) regarding the two sample t-test under a homogeneous one-way analysis of covariance model. The result is, the conventional two sample t-test for comparing two treatments (k = 2) is conservative according to (2(2) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)<\alpha .$(2) ) under the stratified biased coin randomisation. Following this work, results about the conservativeness of different conventional tests under different models and covariate-adaptive randomisation methods have been obtained by Hu and Hu (2012), Shao and Yu (2013), Ma et al. (2015), Bugni et al. (2018), Ye (2018), and Ye and Shao (2020). In particular, under Type 1 or 2 randomisation schemes described in Section 3, Ye and Shao (2020) proved the conservativeness of the conventional log-rank and score tests for survival analysis, which is a substantial advance in the theory of this area. Unfortunately, no result is available for the conservativeness of conventional tests under minimisation, except for some unrealistic cases. Furthermore, the available results are for particular conventional tests, i.e., no general result is available.

The reason why conventional tests become conservative under some covariate-adaptive randomisation schemes as well as why the result is not available for minimisation can be explained as follows. Many (if not most) conventional tests are ratios with numerators being statistics accessing the plausibility of the null hypothesis $H_0$ and denominators being standard errors estimating the standard deviations of the corresponding numerators. For example, the two sample t-test for testing effect between two treatments (k = 2) is (3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) where $n_a$ is the number of patients assigned to treatment a, ${\overline{Y}}_a$ and $S_a^2$ are the sample mean and sample variance, respectively, based on $Y_i$'s under treatment a; the numerator ${\overline{Y}}_1-{\overline{Y}}_2$ of T in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ) accesses the plausibility of the null hypothesis $H_0:E(Y^{\left(1\right)}-Y^{\left(2\right)})=0$, and the denominator $\sqrt{S_1^2/n_1+S_2^2/n_2}$ estimates the asymptotic standard deviation of ${\overline{Y}}_1-{\overline{Y}}_2$. Under a Type 1 or 2 covariate-adaptive randomisation scheme, it is usually true that the numerator of T in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ) still measures the plausibility of $H_0$, and the denominator of T in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ) is too large because the Type 1 or 2 covariate-adaptive randomisation scheme typically reduces the variation of numerator after enforcing the balancedness of treatment assignments. Specifically, ${\overline{Y}}_1$ and ${\overline{Y}}_2$ are independent under simple randomisation but are negatively correlated under Type 1 or Type 2 covariate-adaptive randomisation and, consequently, the variance of ${\overline{Y}}_1-{\overline{Y}}_2$ is smaller under covariate-adaptive randomisation and $S_1^2/n_1+S_2^2/n_2$ still estimates the variance of ${\overline{Y}}_1-{\overline{Y}}_2$ under simple randomisation. The reduction in variation together with the fact that the denominator of conventional test does not account for this reduction lead to the conservativeness of conventional test.

As we discussed in Section 3, the stratified permuted block and biased coin randomisation schemes are Type 1 and the stratified urn designs are Type 2. Hence conventional tests are conservative under these randomisation schemes.

The minimisation, however, is neither Type 1 nor Type 2 (Ye & Shao, 2020). The only available result on the asymptotic distribution of ${\overline{Y}}_1-{\overline{Y}}_2$ under minimisation is obtained (Ma et al., 2015) under a very restrictive and nearly unrealistic condition that not only the relationship between the observed response $Y_i$ and ${\mathit{Z}}_i$ is linear, but also all components of ${\mathit{Z}}_i$ are independent and there is no other covariate in the linear model. Because the minimisation only enforces the marginal balancedness of treatment assignments, its asymptotic properties are very complicated and a general result about the asymptotic distribution of a simple statistic like ${\overline{Y}}_1-{\overline{Y}}_2$ is not available. Some progress has been made in some recent work (Hu & Zhang, 2020), but the problem is not completely solved.

4.3. Development of valid tests

We now return to address Question (A) in Section 1. Although a conservative test controls the Type I error rate, it may lose power of the test and, thus, may not be appreciated by clinical trialists.

From the discussion in Section 4.1, a conventional test is valid according to Equation (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) if (i)–(ii) hold, but (ii) requires prefect modelling that may be unrealistic, since model misspecification often occurs especially when there are many covariates. The discussion in Section 4.2 actually suggests that we modify the denominator of a conventional test to develop a valid test under covariate-adaptive randomisation. The first result was also obtained by Shao et al. (2010) who proposed a bootstrap variance estimator for the two sample t-test with a component of re-generating treatment assignments in every bootstrap sample to account for the correct variation under the stratified biased coin randomisation. The resulting bootstrap test replaces the denominator of two sample t-test in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ) by the squared root of the bootstrap variance estimator and is valid according to (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ). This bootstrap method can be extended to modifying many other conventional tests, for Type 1 or 2 covariate-adaptive randomisation scheme (Shao & Yu, 2013; Ye & Shao, 2020).

With some effort on deriving the asymptotic distribution of the numerator of a conventional test under Type 1 or 2 covariate-adaptive randomisation, a valid test can also be constructed by correctly estimating the asymptotic variance of the numerator (Ye, 2018; Ye & Shao, 2020). For the conventional two sample t-test in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ), for example, Ye (2018) showed that a valid test under stratified biased coin randomisation can be obtained by replacing the denominator of the t-test by $2\sqrt{\sum _{\mathit{z}}n\left(\mathit{z}\right)S^2\left(\mathit{z}\right)}/n$, where $S^2\left(\mathit{z}\right)$ is the sample variance based on $Y_i$'s in stratum $\mathit{Z}=\mathit{z}$. Compared with the bootstrap, this approach does not require a large amount of computation and has another advantage to be discussed later.

Perhaps a better approach is to directly derive a valid test based on a given covariate-adaptive randomisation scheme or a general group of randomisation schemes. This will be discussed in Section 5 when we consider general inference procedures.

Another stream of methods is based on re-randomisation or permutation, e.g., Simon and Simon (2011), Kaiser (2012), and Bugni et al. (2018). In the rest of this section, we discuss in details about the re-randomisation approach in Simon and Simon (2011), which is somewhat similar to the bootstrap method. Consider k = 2 and $H_0:Y^{\left(1\right)}\sim Y^{\left(2\right)}$. Under $H_0$, $Y^{\left(1\right)}$ and $Y^{\left(2\right)}$ are exchangeable so that we create potential outcome ${\tilde{Y}}_i^{\left(1\right)}={\tilde{Y}}_i^{\left(2\right)}=Y_i$ for patient i. Let $\mathcal{A}=(A_1,\dots ,A_n)$ be the observed treatment assignments under the given covariate-adaptive randomisation scheme. Any test T can be written as $T(\mathcal{A},\mathcal{O})$, where $\mathcal{O}=\{{\tilde{Y}}_i^{\left(1\right)},{\tilde{Y}}_i^{\left(2\right)},{\mathit{X}}_i,i=1,\dots ,n\}$. Let $\mathcal{C}=(C_1,\dots ,C_n)$ be randomly generated treatment assignments under the same randomisation scheme, i.e., $\mathcal{C}\sim \mathcal{A}$ conditioned on $\mathit{Z}$, $T(\mathcal{C},\mathcal{O})$ be $T(\mathcal{A},\mathcal{O})$ with $\mathcal{A}$ replaced by $\mathcal{C}$, and let $F_{\mathcal{O}}$ be the cumulative conditional distribution function of $T(\mathcal{C},\mathcal{O})$ given $\mathcal{O}$. From the probability theory, $\begin{array}{rl}& P\left\{\phantom{|}T(\mathcal{C},\mathcal{O})<F_{\mathcal{O}}^{-1}\left(\alpha /2\right)\text{or}T(\mathcal{C},\mathcal{O})\right.\\ & \left.>F_{\mathcal{O}}^{-1}(1-\alpha /2)|\mathcal{O}\right\}\le \alpha .\end{array}$Hence, unconditionally, under $H_0$, $\begin{array}{rl}& P\{T(\mathcal{C},\mathcal{O})<F_{\mathcal{O}}^{-1}\left(\alpha /2\right)\text{or}T(\mathcal{C},\mathcal{O})>F_{\mathcal{O}}^{-1}(1-\alpha /2)\}\\ & \le \alpha \end{array}$and if we reject $H_0$ if and only if T is outside of the interval $\left[F_{\mathcal{O}}^{-1}\right(\alpha /2),F_{\mathcal{O}}^{-1}(1-\alpha /2\left)\right]$, then this T has Type I error rate $\le \alpha$ for every n.

Two issues remain to be considered. The first one is that the quantile $F_{\mathcal{O}}^{-1}\left(r\right)$ usually has no explicit form and approximation such as Monte Carlo is needed. The second issue is that this method may be conservative for every n, because $T(\mathcal{C},\mathcal{O})$ with random $\mathcal{C}$ is discrete. At this stage, it is still unknown whether result (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) holds for this method, since the previous argument shows that the left-hand side of (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) $\le \alpha$, but we cannot prove the equality in (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) holds, i.e., we cannot rule out the possibility that (2(2) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)<\alpha .$(2) ) actually holds so that the re-randomisation method is conservative.

4.4. Tests in survival analysis

We review some available theory for survival analysis, since covariate-adaptive randomisation has a long history of application in survival trials. In fact, all 7 articles in the 2018 New England Journal of Medicine cited in Section 1 are about survival trials.

For simplicity, we focus on the case of k = 2.

The data structure for survival analysis is described in the beginning of Section 2, where the potential outcome $Y^{\left(a\right)}=min(T^{\left(a\right)},C^{\left(a\right)})$, $T^{\left(a\right)}$ is the potential survival, and $C^{\left(a\right)}$ is the potential censoring, under treatment a. It is typically assumed that conditional on covariate $\mathit{X}$, $T^{\left(a\right)}$ and $C^{\left(a\right)}$ are independent and the ratio $P(C^{\left(1\right)}\ge t|\mathit{X})/P(C^{\left(2\right)}\ge t|\mathit{X})$ is a function of t only.

The most common analysis in survival trials is testing whether two treatments have different effect on the conditional distributions of $T^{\left(a\right)}$ given $\mathit{X}$. Let $\lambda (t,\mathit{x},a)$ be the underlying hazard function of $T^{\left(a\right)}$ given $\mathit{X}=\mathit{x}$, a = 1, 2. The null hypothesis of interest is $H_0:\lambda (t,\mathit{x},1)=\lambda (t,\mathit{x},2)$ for all possible t and $\mathit{x}$.

Without imposing any model, a conventional nonparametric test for $H_0$ is the log-rank test (4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) where $S_a\left(t\right)=\sum _{i=1}^{n}I(A_i=e_a)I(Y_i^{\left(a\right)}\ge t)$, $I\left(C\right)$ is the indicator of event C, $S\left(t\right)=S_1\left(t\right)+S_2\left(t\right)$, $N_i\left(t\right)=I(A_i=e_1)N_i^{\left(1\right)}\left(t\right)+I(A_i=e_2)N_i^{\left(2\right)}\left(t\right)$, and $N_i^{\left(a\right)}\left(t\right)=I(T_i^{\left(a\right)}\le C_i^{\left(a\right)})I(Y_i^{\left(a\right)}\le t)$, a = 1, 2. Similar to the two sample t-test in (3(3) $T=\frac{{\overline{Y}}_1-{\overline{Y}}_2}{\sqrt{S_1^2/n_1+S_2^2/n_2}},$(3) ), the log-rank test that is valid according to (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) under simple randomisation is conservative in the sense of (2(2) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)<\alpha .$(2) ) under Type 1 or 2 covariate-adaptive randomisation, because the denominator of T in (4(4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) ) is too large as a standard error for the numerator of T. A valid modified log-rank test is derived by replacing the denominator of T with the squared root of a stratified variance estimator given in Formula (20) of Ye and Shao (2020).

In survival analysis, the following Cox proportional hazard model is very popular: (5) $\lambda (t,\mathit{x},a)={\lambda }_0\left(t\right)\exp (\theta a+{\text{β}}^{\mathrm{T}}\mathit{x}),$(5) where θ is an unknown parameter, ${\text{β}}^{\mathrm{T}}$ is the transpose of a vector β of unknown parameters, and ${\lambda }_0\left(t\right)$ is an unspecified baseline hazard function. If the Cox model is correct, then the null hypothesis is the same as $H_0:\theta =0$, and a score test of $H_0$ can be derived using the partial likelihood under the Cox model. The idea is that the score test is more powerful than the log-rank test if the Cox model is correct. Even if the Cox model could be misspecified, it can be used as a working model under the model-assisted approach, i.e., a model is used to assist the derivation of an inference procedure that is efficient when the model is correct and is still asymptotically valid when the model is incorrect.

Under simple randomisation, a valid model-assisted score test was derived (DiRienzo & Lagakos, 2002; Kong & Slud, 1997; Lin & Wei, 1989), which is often more powerful than the log-rank test in (4(4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) ) without using any covariates. This conventional score test, however, is shown in Ye and Shao (2020) to be conservative under Type 1 or 2 covariate-adaptive randomisation, because of the same reason that the denominator of the score test is too large as a standard error. Again, we can obtain a valid score test by replacing the denominator with the squared root of a stratified variance estimator (Ye & Shao, 2020).

We can also apply the bootstrap or re-randomisation discussed in Section 4.3 to construct valid tests. However, the bootstrap or re-randomisation discussed in Section 4.3 is not correct in survival analysis, unless we assume $P(C^{\left(1\right)}\ge t|\mathit{X})=P(C^{\left(0\right)}\ge t|\mathit{X})$ for all t. The reason is that, to apply the bootstrap or re-randomisation, the observed $(Y_i,{\mathit{X}}_i)$'s have to be exchangeable across i under $H_0$. Under $H_0$, although $T_i^{\left(a\right)}$'s are exchangeable, $C_i^{\left(a\right)}$'s are not unless $P(C^{\left(1\right)}\ge t|\mathit{X})=P(C^{\left(2\right)}\ge t|\mathit{X})$ for all t. Even if the treatment has no effect on the potential survival time, it may have some effect on the potential censoring due to some practical reasons.

5. Valid inference

We have already discussed to some extent how to construct valid tests under covariate-adaptive randomisation. There are a few shortcomings in those available results reviewed in Section 4. First, an obvious one is that some results/methods rely on correct specification of a model. Second, all results/methods in Section 4 depend on covariate-adaptive randomisation schemes; in particular, Type 1 or 2 randomisation method is required, which excludes the minimisation. Third, only testing hypotheses is considered, not other inference such as confidence sets. Finally, all methods in Section 4 are modifications of conventional procedures.

In this section, we would like to address the following re-phrased Question (A) in Section 1:

1. Can we develop an inference procedure valid under covariate-adaptive randomisation with very little model assumption?

5.1. Testing in survival analysis

We begin with the log-rank test for survival data in the case of k = 2. The stratified log-rank test (Peto et al., 1976) is simply the log-rank test in (4(4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) ) stratified with all levels of the discrete covariate $\mathit{Z}$ utilised in covariate-adaptive randomisation: (6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) where $L\left(\mathit{z}\right)$ is the stratum of patients with ${\mathit{Z}}_i=\mathit{z}$, $S_a(t,\mathit{z})=\sum _{i\in L\left(\mathit{z}\right)}I(A_i=e_a)I(Y_i^{\left(a\right)}\ge t)$, and $S(t,\mathit{z})=S_1(t,\mathit{z})+S_2(t,\mathit{z})$. Although the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) exhibits nice empirical properties under covariate-adaptive randomisation (Lachin et al., 1988; Xu et al., 2016) and has been used for a long time, the first proof of its validity according to (1(1) $\underset{P\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{r}{H}_0}{sup}\underset{n\to \mathrm{\infty }}{lim}P\left(T\notin [z_{\alpha /2},z_{1-\alpha /2}]\right)=\alpha$(1) ) comes from Ye and Shao (2020) with some efforts. The proof actually shows that the stratified log-rank test is valid for any covariate-adaptive randomisation method, including the minimisation, as long as the minimal conditions (D1) –(D2) are satisfied.

Why does stratification make so much difference? Recall that in Section 4.1 we comment that a test will be valid if two conditions are satisfied: (i) $\mathit{Z}$ used in randomisation is also used in constructing the test and (ii) a correct model is used to derive the test. Note that the stratification with strata being levels of $\mathit{Z}$ can be viewed as a kind of modelling based on the discrete covariate $\mathit{Z}$, and such modelling is always correct. Thus (ii) has been met if we stratify using $\mathit{Z}$. To meet (i), we must fully stratify, i.e., use all strata defined by joint levels of $\mathit{Z}$, not partially stratify. It can be shown that if we combine some strata in the construction of the stratified log-rank test, then the resulting test is not valid.

The only issue with the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) is that it is not efficient if $\mathit{Z}\ne \mathit{X}$, i.e., $\mathit{X}$ contains more information than $\mathit{Z}$. In fact, we cannot definitely tell whether the stratified log-rank test is more powerful than the unstratified log-rank test in (4(4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) ) under simple randomisation, which is similar to the issue of a stratified sample mean may not be always more efficient than the unstratified sample mean in survey sampling. Ye and Shao (2020) showed by simulation that a modified log-rank test that replaces the denominator of T in (4(4) $\begin{array}{rl}T& =\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1\left(t\right)}{S\left(t\right)}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{i=1}^n{\int }_0^{\mathrm{\infty }}\frac{S_1\left(t\right)S_2\left(t\right)}{\left\{S\right(t){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(4) ) with a stratified standard error may be more powerful than the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ). The efficiency issue will be further considered in Section 6.

5.2. Inference on average or quantile treatment effect

Next, we consider inference on the population mean difference $\theta =E(Y^{\left(a\right)}-Y^{\left(b\right)})$ with any two fixed treatments a and b in a trial with $k\ge 2$ treatment arms. As the development of inference procedures often starts with finding estimators of the parameter of interest, we first review some available estimators of θ.

The simplest estimator of θ is the sample mean difference ${\overline{Y}}_a-{\overline{Y}}_b$, where ${\overline{Y}}_a$ is the sample mean of $Y_i$'s under treatment $a=1,\dots ,k$. Bugni et al. (2018) proposed another estimator called the strata fixed effect estimator in their Section 4.2. The asymptotic distributions of ${\overline{Y}}_a-{\overline{Y}}_b$ and the strata fixed effect estimator have been derived under Type 1 or 2 covariate-adaptive randomisation, but they are not available for Type 3 covariate-adaptive randomisation such as minimisation due to the lack of theory on Type 3 methods.

The following post-stratified estimator of θ, similar to the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ), is proposed by Bugni et al. (2019) and Ye et al. (2020), (7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) where ${\overline{Y}}_a\left(\mathit{z}\right)$ is the sample mean of $Y_i$'s from patients in post-stratum $L\left(\mathit{z}\right)$ under treatment $a=1,\dots ,k$. If the weight $n\left(z\right)/n$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) is replaced by the population weight $P(\mathit{Z}=\mathit{z})$, then ${\hat{\theta }}_S$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) is the stratified estimator in survey sampling. Since $P(\mathit{Z}=\mathit{z})$ is substituted by $n\left(z\right)/n$ and $L\left(\mathit{z}\right)$ is formed after $\mathit{Z}$ is observed, the estimator ${\hat{\theta }}_S$ is referred to as post-stratified estimator in survey sampling.

Applying different techniques, Bugni et al. (2019) and Ye et al. (2020) independently established that, if (C1) and (D1)–(D2) hold and the second order moments of $Y^{\left(a\right)}$ and $Y^{\left(b\right)}$ are finite, then (8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) where $\begin{array}{rl}{\sigma }_S^2& =E\left\{\mathrm{v}\mathrm{a}\mathrm{r}\left(Y^{\left(a\right)}|\mathit{Z}\right)/{\pi }_a+\mathrm{v}\mathrm{a}\mathrm{r}\left(Y^{\left(b\right)}|\mathit{Z}\right)/{\pi }_b\right\}\\ & +\mathrm{v}\mathrm{a}\mathrm{r}\left\{E\right(Y^{\left(a\right)}-Y^{\left(b\right)}|\mathit{Z}\left)\right\}.\end{array}$Result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) is model free, i.e., only (C1) and the second-order moments of the potential outcomes are required. It is applicable to any covariate-adaptive randomisation method satisfying (D1)–(D2), most noticeably the minimisation for which very little is known about its theoretical property, as the minimisation is neither Type 1 nor Type 2. Another interesting fact is that the limiting variance ${\sigma }_S^2$ is invariant with respect to randomisation methods. Hence, not only result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) holds for any covariate-adaptive randomisation method as long as the minimal (D1) –(D2) are satisfied, but also ${\hat{\theta }}_S$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) has the same asymptotic distribution and efficiency regardless of which randomisation scheme is used for treatment assignments. Such kind of result has not be seen in the literature except that Ye and Shao (2020) showed that the asymptotic distribution of the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) is invariant to the randomisation schemes. Existing results in the literature (Bugni et al., 2018; Ma et al., 2015; Shao & Yu, 2013; Shao et al., 2010) are typically dependent with randomisation methods and many of them are not applicable to Type 3 methods such as the minimisation.

When the covariate-adaptive randomisation scheme is Type 1, result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) also holds with ${\hat{\theta }}_S$ replaced by the strata fixed effect estimator in Bugni et al. (2018). In general, however, ${\hat{\theta }}_S$ is asymptotically more efficient than the strata fixed effect estimator or the simple estimator ${\overline{Y}}_a-{\overline{Y}}_b$.

For inference on θ under any type covariate-adaptive randomisation, if ${\hat{\theta }}_S$ is adopted to estimate θ, then all we need to do is to derive an estimator ${\hat{\sigma }}_S^2$ of ${\sigma }_S^2$ that is consistent, i.e., ${\hat{\sigma }}_S^2-{\sigma }_S^2\Rightarrow 0$ under any type covariate-adaptive randomisation. This is actually not difficult once we establish a result like (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ). It is shown in Ye et al. (2020) that a consistent estimator of ${\sigma }_S^2$ under any type covariate-adaptive randomisation is $\begin{array}{rl}{\hat{\sigma }}_S^2& =\frac{1}{n}\sum _{\mathit{z}}{n}^2\left(\mathit{z}\right)\left\{\frac{S_a^2\left(\mathit{z}\right)}{n_a\left(\mathit{z}\right)}+\frac{S_b^2\left(\mathit{z}\right)}{n_b\left(\mathit{z}\right)}\right\}\\ & +\frac{1}{n}\sum _{\mathit{z}}n\left(\mathit{z}\right){\left\{{\overline{Y}}_a\left(\mathit{z}\right)-{\overline{Y}}_b\left(\mathit{z}\right)\right\}}^2-{\hat{\theta }}_S^2,\end{array}$where $n_a\left(\mathit{z}\right)$ and $S_a^2\left(\mathit{z}\right)$ are the sample size and sample variance of $Y_i$'s, respectively, of the patients in stratum $\mathit{Z}=\mathit{z}$ and under treatment a.

Under any randomisation scheme satisfying (D1)–(D2), an asymptotically valid $(1-\alpha )\mathrm{\%}$ confidence interval for θ has limits ${\hat{\theta }}_S\pm z_{1-\alpha /2}{\hat{\sigma }}_S$, where $z_{1-\alpha }$ is the quantile of the standard normal distribution.

More estimators of the average treatment effect θ are considered in Section 6.

We now consider inference on another important parameter, the quantile treatment effect defined as $q_{\tau }^{\left(a\right)}-q_{\tau }^{\left(b\right)}$ in Section 2, where $q_{\tau }^{\left(a\right)}$ is the τth quantile of the distribution of $Y^{\left(a\right)}$ under treatment a and τ is a fixed fraction.

Unlike the means, for quantiles we cannot use differences as in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ). Instead, we estimate $q_{\tau }^{\left(a\right)}$ and $q_{\tau }^{\left(b\right)}$ separately, and then take a difference of estimates. Under treatment a, we estimate the marginal distribution of $Y^{\left(a\right)}$ at a fixed point y as (9) ${\hat{F}}^{\left(a\right)}\left(y\right)=\frac{1}{n}\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n_a\left(\mathit{z}\right)}\sum _{i\in L\left(\mathit{z}\right)}I(A_i=e_a)I(Y_i^{\left(a\right)}\le y),$(9) $a=1,\dots ,k$. Then, $q_{\tau }^{\left(a\right)}$ is estimated by ${\hat{q}}_{\tau }^{\left(a\right)}=$ the τth quantile of ${\hat{F}}^{\left(a\right)}$, and $q_{\tau }^{\left(a\right)}-q_{\tau }^{\left(b\right)}$ is estimated as ${\hat{q}}_{\tau }^{\left(a\right)}-{\hat{q}}_{\tau }^{\left(b\right)}$. For inference on quantiles, however, a simple estimator of the asymptotic variance of ${\hat{q}}_{\tau }^{\left(a\right)}$ may not be easily obtained. Methods such as the bootstrap or Woodruff's interval may be applied (Shao, 2003).

The stratification in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ), (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) or (9(9) ${\hat{F}}^{\left(a\right)}\left(y\right)=\frac{1}{n}\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n_a\left(\mathit{z}\right)}\sum _{i\in L\left(\mathit{z}\right)}I(A_i=e_a)I(Y_i^{\left(a\right)}\le y),$(9) ), together with the asymptotic theory, provides a solid foundation for valid and model free inference after covariate-adaptive randomisation and, thus, it largely eliminates the concern and controversy as discussed by regulatory agencies about the use of covariate-adaptive randomisation such as minimisation.

Combining the results and discussions in this section and Section 5.1, we reach a general conclusion that a valid inference procedure can be obtained as long as the covariate $\mathit{Z}$ utilised in covariate-adaptive randomisation is fully used in the construction of inference procedure. A simple way to do this is to use all joint levels of $\mathit{Z}$ as strata.

It can be seen that the conditions needed for this conclusion is much weaker than (i) and (ii) stated in Section 4.1, but (i)–(ii) in Section 4.1 are considered for the validity of a conventional test under covariate-adaptive randomisation.

5.3. Effect of types of randomisation schemes

Result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) about the asymptotic distribution of ${\hat{\theta }}_S$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) is invariant to any types of randomisation schemes described in Section 3. But this does not imply that the stratification in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) or in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) is the best way for inference, especially when problems other than the inference on average treatment effect are considered. An example is that the modified log-rank test in Ye and Shao (2020) may be more powerful than the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ), as discussed in the end of Section 5.1.

If an inference procedure is not invariant to different randomisation schemes, then it is interesting to find out which randomisation scheme, or which type, provides better inference procedures. For the modified log-rank test in Ye and Shao (2020), it is more powerful when a Type 1 randomisation scheme is used, rather than the Type 2 or 3. The same may be true for any inference procedure not invariant to different randomisation schemes. For different Type 1 methods, such as the stratified permuted block and the stratified biased coin methods, so far there is no result indicating that the inference procedures based on these two randomisation schemes have different performances.

6. Efficiency considerations

Question (C) in Section 1 is about whether a test under covariate-adaptive randomisation can be more powerful than it is under simple randomisation. Another question is, if $\mathit{Z}$ is used in randomisation and stratification as in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) or (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) and if $\mathit{X}$ contains more information than $\mathit{Z}$, can we obtain more powerful tests or more efficient estimators by utilising covariate information in $\mathit{X}$ that is not in $\mathit{Z}$? Note that $\mathit{X}$ may contain a component that is not in $\mathit{Z}$ but is related with the potential responses $Y^{\left(1\right)},\dots ,Y^{\left(k\right)}$, or some components of $\mathit{Z}$ are discretised components of $\mathit{X}$ and the remaining information after discretisation is still useful in predicting the potential responses.

6.1. Adjusting for covariates

We first consider the second question in the estimation of $\theta =E(Y^{\left(a\right)}-Y^{\left(b\right)})$ for two fixed treatments a and b. Let $\mathit{U}$ be a function of $\mathit{X}$ that we want to further utilise in improving the efficiency of ${\hat{\theta }}_S$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ). Since the information generated by $\mathit{Z}$ is not in that of $\mathit{U}$, we assume that $\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{U}|\mathit{Z}=\mathit{z})$ is positive definite for every $\mathit{z}$.

For model free estimation and inference, we do not want to impose any model between the potential responses and $\mathit{U}$. In fact, it is hard to find a correct model within each stratum $\mathit{Z}=\mathit{z}$, if we still apply stratification in estimating θ. How do we adjust for covariates without using a model? Ye et al. (2020) adopted the model-assisted generalised regression approach in survey sampling, first discussed in Cassel et al. (1976) and studied extensively in the literature, for example, Särndal et al. (2003), Shao and Wang (2014), and Ta et al. (2020).

In this section, we review some results from Ye et al. (2020). Let ${\mathit{U}}_i$ be the covariate $\mathit{U}$-value of patient i, and for each $\mathit{z}$, let ${\overline{\mathit{U}}}_a\left(\mathit{z}\right)$ be the sample mean of ${\mathit{U}}_i$'s of patients in stratum $L_a\left(\mathit{z}\right)=\{i:{\mathit{Z}}_i=\mathit{z}\text{ under treatment}a\}$, and $\begin{array}{rl}{\widehat{\text{β}}}_a\left(\mathit{z}\right)& ={\left[\sum _{i\in L_a\left(\mathit{z}\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}\left)\right\}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}){\}}^{\mathrm{T}}\right]}^{-1}\\ & \times \sum _{i\in L_a\left(\mathit{z}\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}\left)\right\}{Y}_i.\end{array}$Within treatment a and stratum $L\left(\mathit{z}\right)=\{\mathit{Z}=\mathit{z}\}$, ${\widehat{\text{β}}}_a\left(\mathit{z}\right)$ is the least squares estimator of the coefficient vector in front of $\mathit{U}$ under a linear model between $Y^{\left(a\right)}$ and $\mathit{U}$, but the model is not required to be correct. An estimator of θ following ${\hat{\theta }}_S$ but further adjusting for covariate $\mathit{U}$ is (Ye et al., 2020) $\begin{array}{rl}{\hat{\theta }}_A& =\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left[{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z})-\{{\overline{\mathit{U}}}_a\left(\mathit{z}\right)\\ & -\overline{\mathit{U}}\left(\mathit{z}\right){\}}^{\mathrm{T}}{\widehat{\text{β}}}_a\left(\mathit{z}\right)+\left\{{\overline{\mathit{U}}}_b\right(\mathit{z})-\overline{\mathit{U}}(\mathit{z}\left){\}}^{\mathrm{T}}{\widehat{\text{β}}}_b\right(\mathit{z}\left)\right],\end{array}$where $\overline{\mathit{U}}\left(\mathit{z}\right)$ is the sample mean of ${\mathit{U}}_i$'s of all patients in stratum $L\left(\mathit{z}\right)$.

An alternative estimator ${\hat{\theta }}_B$ of θ in Ye et al. (2020) is obtained by replacing both ${\widehat{\text{β}}}_a\left(\mathit{z}\right)$ and ${\widehat{\text{β}}}_b\left(\mathit{z}\right)$ in the definition of ${\hat{\theta }}_A$ with a combined estimator $\begin{array}{rl}& \widehat{\text{β}}\left(\mathit{z}\right)\\ & ={\left[\sum _{a=1}^k\sum _{i\in L\left(\mathit{z}\right),A_i=a}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}\left)\right\}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}){\}}^{\mathrm{T}}\right]}^{-1}\\ & \times \sum _{a=1}^k\sum _{i\in L\left(\mathit{z}\right),A_i=a}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(\mathit{z}\left)\right\}{Y}_i.\end{array}$When k>2, both $\overline{\mathit{U}}\left(z\right)$ and $\widehat{\text{β}}\left(\mathit{z}\right)$ involve data from patients in treatment arms other than a and b.

The following result parallel to result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) is established in Ye et al. (2020). If (C1) and (D1)–(D2) hold and the second order moments of $Y^{\left(a\right)}$ and $\mathit{U}$ are finite, then (10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) where $\begin{array}{rl}{\sigma }_A^2& =E\left[\mathrm{v}\mathrm{a}\mathrm{r}\right\{Y^{\left(a\right)}-{\mathit{U}}^{\mathrm{T}}{\text{β}}_a\left(\mathit{Z}\right)|\mathit{Z}\}/{\pi }_a\\ & +\mathrm{v}\mathrm{a}\mathrm{r}\{Y^{\left(b\right)}-{\mathit{U}}^{\mathrm{T}}{\text{β}}_b(\mathit{Z}\left)|\mathit{Z}\right\}/{\pi }_b]\\ & +E\left[\right\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right){\}}^{\mathrm{T}}\\ & \times \mathrm{v}\mathrm{a}\mathrm{r}\left(\mathit{U}|\mathit{Z}\right)\left\{{\text{β}}_a\right(\mathit{Z})-{\text{β}}_b(\mathit{Z}\left)\right\}]\\ & +\mathrm{v}\mathrm{a}\mathrm{r}\left\{E\right(Y^{\left(a\right)}-Y^{\left(b\right)}|\mathit{Z}\left)\right\},\\ {\sigma }_B^2& =E\left[\mathrm{v}\mathrm{a}\mathrm{r}\right\{Y^{\left(a\right)}-{\mathit{U}}^{\mathrm{T}}\text{β}\left(\mathit{Z}\right)|\mathit{Z}\}/{\pi }_a\\ & +\mathrm{v}\mathrm{a}\mathrm{r}\{Y^{\left(b\right)}-{\mathit{U}}^{\mathrm{T}}\text{β}(\mathit{Z}\left)|\mathit{Z}\right\}/{\pi }_b]\\ & +\mathrm{v}\mathrm{a}\mathrm{r}\left\{E\right(Y^{\left(a\right)}-Y^{\left(b\right)}|\mathit{Z}\left)\right\},\end{array}$${\text{β}}_a\left(\mathit{z}\right)=\left\{\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathit{U}|\mathit{Z}=\mathit{z}\left){\}}^{-1}\mathrm{c}\mathrm{o}\mathrm{v}\right(\mathit{U},Y^{\left(a\right)}|\mathit{Z}=\mathit{z})$, $a=1,\dots ,k$, and $\text{β}\left(\mathit{z}\right)=\sum _{a=1}^k{\pi }_a{\text{β}}_a\left(\mathit{z}\right)$.

Several conclusions can be made from result (10(10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) ). First, result (10(10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) ) is model free and invariant with respect to covariate-adaptive randomisation schemes, as long as the minimal (D1)–(D2) hold.

Second, from the definitions of ${\sigma }_S^2$ and ${\sigma }_A^2$, it is shown in Ye et al. (2020) that $\begin{array}{rl}& {\sigma }_S^2-{\sigma }_A^2\\ & =E\left[\left\{{\pi }_b{\text{β}}_a\right(\mathit{Z})+{\pi }_a{\text{β}}_b(\mathit{Z}){\}}^{\mathrm{T}}\right.\\ & \times \left.\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{U}\mid \mathit{Z})\left\{{\pi }_b{\text{β}}_a\right(\mathit{Z})+{\pi }_a{\text{β}}_b(\mathit{Z}\left)\right\}\phantom{{}^{\mathrm{T}}}\right]\\ & \times \left\{{\pi }_a{\pi }_b\right({\pi }_a+{\pi }_b){\}}^{-1}\\ & +E\left[{\left\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right)\right\}}^{\mathrm{T}}\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{U}\mid \mathit{Z})\left\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right)\right\}\phantom{{\left\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right)\right\}}^{\mathrm{T}}}\right]\\ & \times \left\{\right({\pi }_a+{\pi }_b{)}^{-1}-1\}\end{array}$and, hence, adjusting covariate $\mathit{U}$ always gains efficiency, i.e., ${\hat{\theta }}_A$ is asymptotically more efficient than ${\hat{\theta }}_S$, unless (11) $\begin{array}{rl}& {\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)=0\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\text{for every}\mathit{z}\text{,}\end{array}$(11) in which case ${\hat{\theta }}_S$ and ${\hat{\theta }}_A$ have the same asymptotic efficiency. When there are more than two treatments, $1-{\pi }_a-{\pi }_b>0$ and, consequently, (11(11) $\begin{array}{rl}& {\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)=0\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\text{for every}\mathit{z}\text{,}\end{array}$(11) ) holds only when ${\text{β}}_a\left(\mathit{z}\right)={\text{β}}_b\left(\mathit{z}\right)=0$ for every $\mathit{z}$, i.e., $\mathit{U}$ is uncorrelated with the potential responses $Y^{\left(a\right)}$ and $Y^{\left(b\right)}$ after conditioning on $\mathit{Z}$ so that adjusting for $\mathit{U}$ is unnecessary. When there are only two treatments, (11(11) $\begin{array}{rl}& {\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)=0\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\text{for every}\mathit{z}\text{,}\end{array}$(11) ) also holds if ${\pi }_a={\pi }_b=1/2$ and ${\text{β}}_a\left(\mathit{z}\right)=-{\text{β}}_b\left(\mathit{z}\right)$ for every $\mathit{z}$.

Third, from the definitions of ${\sigma }_A^2$ and ${\sigma }_B^2$, it can be shown (Ye et al., 2020) that $\begin{array}{rl}& {\sigma }_B^2-{\sigma }_A^2\\ & =E\left[\left\{{\text{β}}_a\right(\mathit{Z})-\text{β}(\mathit{Z}\left){\}}^{\mathrm{T}}\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathit{U}\mid \mathit{Z}\left)\right\{{\text{β}}_a\left(\mathit{Z}\right)-\text{β}\left(\mathit{Z}\right)\}\right]{\pi }_a^{-1}\\ & +E\left[\left\{{\text{β}}_b\right(\mathit{Z})-\text{β}(\mathit{Z}\left){\}}^{\mathrm{T}}\mathrm{v}\mathrm{a}\mathrm{r}\right(\mathit{U}\mid \mathit{Z}\left)\right\{{\text{β}}_b\left(\mathit{Z}\right)-\text{β}\left(\mathit{Z}\right)\}\right]{\pi }_b^{-1}\\ & -E\left[{\left\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right)\right\}}^{\mathrm{T}}\mathrm{v}\mathrm{a}\mathrm{r}(\mathit{U}\mid \mathit{Z})\left\{{\text{β}}_a\left(\mathit{Z}\right)-{\text{β}}_b\left(\mathit{Z}\right)\right\}\right].\end{array}$and, hence, ${\hat{\theta }}_A$ is asymptotically more efficient than ${\hat{\theta }}_B$ unless (12) $\begin{array}{rl}\text{β}\left(\mathit{z}\right)& =\frac{{\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)}{{\pi }_a+{\pi }_b}\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\\ & \text{for every}\mathit{z}\text{,}\end{array}$(12) in which case ${\hat{\theta }}_B$ and ${\hat{\theta }}_A$ have the same asymptotic efficiency.

Note that $\widehat{\text{β}}\left(\mathit{z}\right)$ used in ${\hat{\theta }}_B$ ignores the fact that $\mathrm{c}\mathrm{o}\mathrm{v}(\mathit{U},Y^{\left(a\right)}\mid \mathit{Z}=\mathit{z})$ may depend on treatment a. That is why ${\hat{\theta }}_B$ is asymptotically not as efficient as ${\hat{\theta }}_A$ in general, and ${\sigma }_B^2={\sigma }_A^2$ when these covariances are the same for every a and every $\mathit{z}$, i.e., ${\text{β}}_1\left(\mathit{z}\right)=\cdots ={\text{β}}_k\left(\mathit{z}\right)$ so that (12(12) $\begin{array}{rl}\text{β}\left(\mathit{z}\right)& =\frac{{\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)}{{\pi }_a+{\pi }_b}\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\\ & \text{for every}\mathit{z}\text{,}\end{array}$(12) ) holds. If (12(12) $\begin{array}{rl}\text{β}\left(\mathit{z}\right)& =\frac{{\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)}{{\pi }_a+{\pi }_b}\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\\ & \text{for every}\mathit{z}\text{,}\end{array}$(12) ) holds, ${\hat{\theta }}_B$ may have better finite sample performance than ${\hat{\theta }}_A$, although two estimators are asymptotically equivalent. An exceptional case for ${\sigma }_A^2={\sigma }_B^2$ is when k = 2 and ${\pi }_1={\pi }_2=1/2$, in which we even do not need ${\text{β}}_a\left(\mathit{z}\right)={\text{β}}_b\left(\mathit{z}\right)$.

In general, ${\hat{\theta }}_B$ may be asymptotically less efficient than ${\hat{\theta }}_S$, i.e., covariate adjustment with only the main effects may hurt efficiency, a perspective in Freedman (2008) and Lin (2013). For example, there are scenarios in which (11(11) $\begin{array}{rl}& {\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)=0\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\text{for every}\mathit{z}\text{,}\end{array}$(11) ) holds but (12(12) $\begin{array}{rl}\text{β}\left(\mathit{z}\right)& =\frac{{\pi }_b{\text{β}}_a\left(\mathit{z}\right)+{\pi }_a{\text{β}}_b\left(\mathit{z}\right)}{{\pi }_a+{\pi }_b}\text{and}\\ & \left\{{\text{β}}_a\right(\mathit{z})-{\text{β}}_b(\mathit{z}\left)\right\}(1-{\pi }_a-{\pi }_b)=0\\ & \text{for every}\mathit{z}\text{,}\end{array}$(12) ) does not.

Finally, inference about θ can be carried out based on (10(10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) ) and the availability of consistent estimators of ${\sigma }_A^2$ and ${\sigma }_B^2$. Some model free consistent variance estimators under any covariate-adaptive randomisation schemes are derived in Ye et al. (2020), which are similar to ${\hat{\sigma }}_S^2$ in Section 5.2.

6.2. Can covariate-adaptive randomisation boost efficiency?

We now address Question (C) raised in Section 1 and the beginning of this section: Can a test (or an inference procedure) under covariate-adaptive randomisation be more efficient than it is under simple randomisation?

For the types of covariate-adaptive randomisation schemes described in Section 3, the answer is no, assuming that exactly the same test is used under simple randomisation or under covariate-adaptive randomisation without adjusting for conservativeness. This answer is based on the first-order asymptotic property. With a fixed n, the test or inference procedure under covariate-adaptive randomisation may perform slightly better due to the balancedness of treatment assignments.

In our previous discussions, a conventional procedure may be conservative under covariate-adaptive randomisation, and a valid procedure can often be constructed by modifying the conventional procedure. This modified procedure can be more efficient than the conventional procedure, but the comparison is not fair because the modified procedure makes some adjustment typically depending on the covariate $\mathit{Z}$.

Then, what is the advantage of applying covariate-adaptive randomisation? It is applied mainly for balancing treatment assignments across prognostic factors, which may be important for reviewing clinical results and other practical considerations.

There is a stream of developments and results in balancing discrete or continuous covariates and increasing estimation efficiency at the same time (Atkinson, 1982, 1999, 2002; Baldi Antognini & Zagoraiou, 2011; Rosenberger & Sverdlov, 2008; Senn et al., 2010). The approaches are typically model-based and the gains in efficiency may be from the second-order asymptotics.

Boosting efficiency can also be achieved by adjusting covariates under simple randomisation with less effort compared with applying covariate-adaptive randomisation, which is discussed next.

6.3. Designing versus modelling

Utilising covariate $\mathit{Z}$ in randomisation can be viewed as a kind of designing for better quality of data, although this is not the same as what in the traditional experiment design, because in clinical trials we typically cannot control covariate values of patients. Adjusting for covariates, either model-based or model-assisted, fits into the general framework of modelling. In this section, we address the issue of designing versus modelling.

First, consider inference on the average treatment effect $\theta =E(Y^{\left(a\right)}-Y^{\left(b\right)})$. If $\mathit{Z}$ is the only covariate, i.e., $\mathit{X}=\mathit{Z}$, then the conclusion is that designing and modelling (adjusting for covariate) can achieve the same efficiency asymptotically. In this case, ${\hat{\theta }}_S={\hat{\theta }}_A={\hat{\theta }}_B$ and it has the same asymptotic normal distribution under simple randomisation and under any other covariate-adaptive randomisation satisfying (D1)–(D2). The stratification in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) serves the purpose of modelling under simple randomisation, but it is essential for obtaining easy inference under covariate-adaptive randomisation including minimisation.

Consider next the situation where $\mathit{X}\ne \mathit{Z}$ and the covariate $\mathit{U}$ as discussed in Section 6.1 together with $\mathit{Z}$ are available for modelling (the entire covariate $\mathit{X}$ may still contain more information than that from $\mathit{U}$ and $\mathit{Z}$). The conclusion is, modelling with $\mathit{Z}$ and $\mathit{U}$ achieves more efficiency than designing with $\mathit{Z}$ only, and is the same as designing with $\mathit{Z}$ plus an additional modelling with $\mathit{U}$ (adjusting for $\mathit{U}$). This directly comes from result (10(10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) ). Under simple randomisation, ${\hat{\theta }}_A$ in Section 6.1 is the estimator after modelling with $\mathit{Z}$ and $\mathit{U}$, in view of the fact that $\mathit{Z}$ is discrete so that stratification is the same as modelling with $\mathit{Z}$, and its limiting variance is ${\sigma }_A^2$ in (10(10) $\begin{array}{rl}& \sqrt{n}({\hat{\theta }}_A-\theta )\Rightarrow N(0,{\sigma }_A^2)\text{and}\\ & \sqrt{n}({\hat{\theta }}_B-\theta )\Rightarrow N(0,{\sigma }_B^2),\end{array}$(10) ). On the other hand, designing with $\mathit{Z}$ only leads to the estimator ${\hat{\theta }}_S$ in (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ), which has limiting variance ${\sigma }_S^2$ in (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ) regardless of which covariate-adaptive randomisation is applied, and ${\sigma }_A^2\le {\sigma }_S^2$. Finally, designing with $\mathit{Z}$ plus an additional modelling with $\mathit{U}$ leads to estimator ${\hat{\theta }}_A$.

Similar conclusions can be obtained for testing in survival analysis as discussed in Section 4.4. Consider the situation of $\mathit{X}=\mathit{Z}$. Since $\mathit{Z}$ is discrete, the Cox model given by (5(5) $\lambda (t,\mathit{x},a)={\lambda }_0\left(t\right)\exp (\theta a+{\text{β}}^{\mathrm{T}}\mathit{x}),$(5) ) is always correct. Modelling with $\mathit{Z}$ produces the score test under simple randomisation, whereas designing with $\mathit{Z}$ leads to the stratified log-rank test in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ). It is shown in Ye and Shao (2020) that the two tests have the same Pitman's asymptotic efficiency. If $\mathit{X}\ne \mathit{Z}$ and the Cox model (5(5) $\lambda (t,\mathit{x},a)={\lambda }_0\left(t\right)\exp (\theta a+{\text{β}}^{\mathrm{T}}\mathit{x}),$(5) ) with $\mathit{X}$ is correct, then it is shown in Ye and Shao (2020) that the score test under simple randomisation is more efficient than the stratified log-rank test based on designing and stratification with $\mathit{Z}$, in terms of Pitman's asymptotic efficiency. In this case, designing with $\mathit{Z}$ plus an additional modelling leads to the score test. Unlike the case of inference on average treatment effect, however, in survival testing all results for the situation of $\mathit{X}\ne \mathit{Z}$ relies on the correctness of Cox model (5(5) $\lambda (t,\mathit{x},a)={\lambda }_0\left(t\right)\exp (\theta a+{\text{β}}^{\mathrm{T}}\mathit{x}),$(5) ). If model (5(5) $\lambda (t,\mathit{x},a)={\lambda }_0\left(t\right)\exp (\theta a+{\text{β}}^{\mathrm{T}}\mathit{x}),$(5) ) is wrong, then the score test can be less powerful than the unstratified log-rank test.

7. Further research work

We end this review with the following discussion of further research topics in this area.

1. Although some estimation and inference procedures previously discussed have asymptotic distributions invariant to covariate-adaptive randomisation schemes, it may be still important to study and understand the Type 3 randomisation methods such as the minimisation whose properties are unclear at this stage. In particular, the asymptotic property of $D^{\left(a\right)}\left(\mathit{z}\right)$ defined in (D2). Efforts should be made to establish the joint asymptotic distribution of $D^{\left(a\right)}\left(\mathit{z}\right)$ with $\mathit{z}$ being all levels of $\mathit{Z}$. A different direction is to develop more and better covariate-adaptive randomisation schemes. For example, in Section 5.3 we point out that a Type 1 randomisation scheme may produce more efficient inference procedures than a Type 2 or 3 randomisation scheme. Hu and Hu (2012) modified Pocock and Simon's approach and proposed to use an imbalance measure that is a weighted sum of the overall imbalance, marginal imbalance, and strata imbalance. Some effort should be made to study the implementation of this scheme for practical uses.

2. To utilise covariates, we considered the model-assisted generalised regression approach for the estimation of average treatment effect and score test under a working Cox model for testing hypotheses in survival analysis. It is interesting to develop other model-assisted approaches to gain efficiency without relying on models.

3. From result (8(8) $\sqrt{n}({\hat{\theta }}_S-\theta )\Rightarrow N(0,{\sigma }_S^2),$(8) ), if ${\mathit{Z}}^{\prime }$ is another covariate such that the σ-field of ${\mathit{Z}}^{\prime }$ contains the σ-field of $\mathit{Z}$, then the ${\hat{\theta }}_A$ using ${\mathit{Z}}^{\prime }$ in randomisation is asymptotically more efficient than the ${\hat{\theta }}_A$ using $\mathit{Z}$ in randomisation. That is, utilising more covariate information in randomisation can increase asymptotic efficiency. On the other hand, using a $\mathit{Z}$ with too many levels may cause sparsity of data. Some guidance on this may be useful for practical users.

4. The stratification in (6(6) $\begin{array}{rl}T& =\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\left\{A_i-\frac{S_1(t,\mathit{z})}{S(t,\mathit{z})}\right\}\mathrm{d}{N}_i\left(t\right)\\ & \times {\left[\sum _{\mathit{z}}\sum _{i\in L\left(\mathit{z}\right)}{\int }_0^{\mathrm{\infty }}\frac{S_1(t,\mathit{z})S_2(t,\mathit{z})}{\left\{S\right(t,\mathit{z}){\}}^2}\mathrm{d}{N}_i\left(t\right)\right]}^{-1/2},\end{array}$(6) ) or (7(7) ${\hat{\theta }}_S=\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n}\left\{{\overline{Y}}_a\right(\mathit{z})-{\overline{Y}}_b(\mathit{z}\left)\right\},$(7) ) uses all levels of $\mathit{Z}$ as strata. In applications, it is possible that some strata contain very few number of patients or even no patient. Some methods of handling this scenario should be developed to produce asymptotically valid or at least conservative inference procedures, such as combining some strata with small sizes.

5. The result and discussion on inference about quantile treatment effects are very limited. In survival analysis, due to the presence of censoring, the distribution function estimator in (9(9) ${\hat{F}}^{\left(a\right)}\left(y\right)=\frac{1}{n}\sum _{\mathit{z}}\frac{n\left(\mathit{z}\right)}{n_a\left(\mathit{z}\right)}\sum _{i\in L\left(\mathit{z}\right)}I(A_i=e_a)I(Y_i^{\left(a\right)}\le y),$(9) ) has to be replaced by the Kaplan–Meier product-limit type estimator. Furthermore, how to adjust for covariates has not been considered.

6. The bootstrap, re-randomisation and permutation methods described in Section 4.3 are promising alternative tools to the approach of asymptotic distribution plus variance estimation for statistical inference. Two issues have to be addressed. The first one is that the re-randomisation and permutation methods are naturally developed for testing. Applying these tools for inference on parameters other than the average treatment effect requires further development. The other issue is what we discussed in the end of Section 4.4, i.e., the development of bootstrap or re-randomisation methods when censoring distributions conditioned on $\mathit{X}$ can be different under the null hypothesis that the survival distributions conditioned on $\mathit{X}$ are identical.

Acknowledgements

Our research was supported by the National Natural Science Foundation of China (11831008) and the U.S. National Science Foundation (DMS-1914411).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Our research was supported by the National Natural Science Foundation of China (11831008) and the U.S. National Science Foundation (DMS-1914411).

Notes on contributors

Jun Shao

Dr Jun Shao holds a PhD in statistics from the University of Wisconsin-Madison. He is a Professor of Statistics at the University of Wisconsin-Madison. His research interests include variable selection and inference with high dimensional data, sample surveys, and missing data problems.