Rejoinder on ‘Inference after covariate-adaptive randomization: aspects of methodology and theory’

I would like to thank all discussants for their insightful discussions on the topic of statistical inference after covariate-adaptive randomisation, especially for including reviews of some new results and references that are not in my review written more than a year ago. I hope these discussions together with my review will stimulate further studies in this important area having many applications particularly in clinical trials.

My rejoinder focuses on some main points from four separate groups of discussants.

1. The discussion by Drs. Ma, Zhang and Hu

Drs. Ma, Zhang, and Hu's discussion brings some new results or interesting directions for new studies. These include, but are not limited to, robust inference with the randomisation scheme considered in Hu and Hu (2012), inference after covariate-adaptive randomisation with covariate misclassification, unobserved covariates, missing data or non-compliance, and high dimensional covariates. Covariate-adaptive randomisation can also be combined with other adaptive designs such as sequential monitoring (Zhu & Hu, 2019), sample size re-estimation (Li et al., in press), and seamless phase II/III clinical trials (Ma et al., in press).

As I discussed in Section 7 of my review, the theoretical study of covariate-adaptive randomisation schemes such as Pocock and Simon's minimisation is not completed and still important, although some asymptotically valid inference procedures after these randomisation schemes have been derived. Let $\mathit{Z}$ be the discrete covariate used in covariate-adaptive randomisation, $z_1,\dots ,z_c$ be all possible categories of $\mathit{Z}$, ${\pi }_1,\dots ,{\pi }_k$ be the target assignment proportions in the trial, n be the total sample size of the trial, $n_a\left(l\right)$ be the number of patients in treatment a with $\mathit{Z}=z_l$, $a=1,\dots ,k$, $l=1,\dots ,c$, and $n\left(l\right)=n_1\left(l\right)+\cdots +n_k\left(l\right)$. A key result for studying the asymptotic validity of inference procedures after covariate-adaptive randomisation is

(R1)	$\sqrt{n}(\frac{n_a\left(l\right)}{n\left(l\right)}-{\pi }_a,a=1,\dots ,k,l=1,\dots ,c)|{\mathit{Z}}_1,\dots ,{\mathit{Z}}_n\to N(0,D)\text{in distribution},$

i.e. conditioned on n observed values of $\mathit{Z}$, ${\mathit{Z}}_1,\dots ,{\mathit{Z}}_n$, the $k\times c$ dimensional vector whose $(a,l)$th component is $\frac{n_a\left(l\right)}{n\left(l\right)}-{\pi }_a$ converges in distribution to a multivariate normal with mean 0 and some covariance matrix D. Result (R1) holds for type 1 or 2 covariate-adaptive randomisation schemes as described in Section 3 of my review. For Pocock and Simon's minimisation, however, it has not been rigorously shown in general whether or not (R1) holds, although (Ma et al., 2015) showed (R1) under restrictive conditions (a correct linear model between the response and $\mathit{Z}$ whose components are independent). Hu and Zhang (2020) derived the asymptotic normality for $\frac{n_a\left(l\right)}{n\left(l\right)}-{\pi }_a,a=1,\dots ,k$ with a single fixed l, but (R1) requires joint asymptotic normality of the entire vector over all $l=1,\dots ,c$.

Without (R1), asymptotic validity of inference procedures for most covariate-adaptive randomisation schemes including Pocock and Simon's minimisation can still be established under some conditions, which will be further explained in Section 4.

To construct asymptotic valid inference procedures, sometimes the explicit form of the covariance matrix D in (R1) is required.

I am less excited about balancing continuous covariates with covariate-adaptive randomisation. The reason is that covariate-adaptive randomisation is mainly use to enhance the credibility of the results of the trial (EMA, 2015) in terms of the balancedness of treatments across levels of some common discrete baseline prognostic factors such as institution, disease stage, prior treatment, gender, and age group. In fact, balancedness of marginal levels of these discrete prognostic factors is the main concern of agencies such as the EMA or FDA, and that is why Pocock and Simon's minimisation is popular. I do not see a clear motivation of balancing a continuous baseline covariate at the design stage. If efficiency is the concern, we may simply adjust for this continuous covariate in the inference procedure. It is much easier to construct a valid and efficient inference procedure, compared with to derive a valid inference procedure after balancing a continuous baseline covariate. So far, valid inference procedures after balancing a continuous baseline covariate are mostly model-based.

2. The discussion by Drs. Wang, Susukida, Mojtabai, Amin-Esmaeili and Rosenblum

My review mostly focuses on differences of sample means (or quantiles) and tests in survival analysis. The discussion by Drs. Wang, Susukida, Mojtabai, Amin-Esmaeili, and Rosenblum and their article (Wang et al., 2020) open up a wide range of robust inference methods to handle nonlinearity, various outcome types, repeated measures, missing outcomes, etc. Wang et al. (2020) also contain three examples of analyses of trial data, illustrating that the gain due to stratified permuted block randomisation and covariate adjustment could be as high as 36%.

In Section 5.3 of my review, empirical distribution estimators and related quantile estimators valid under covariate-adaptive randomisation are considered. Wang et al. (2020) established the asymptotic normality of Kaplan–Meier estimator under stratified permuted block or biased coin randomisation, which is important for survival analysis. Specifically, they showed that, in a trial with two arms (k = 2), $\left\{\sqrt{n}\right({\hat{S}}_n^{\left(a\right)}-S^{\left(a\right)}),t\in [0,\tau \left]\right\}$, where ${\hat{S}}_n^{\left(a\right)}$ is the Kaplan–Meier estimator of the survival function $S^{\left(a\right)}$, a is a fixed treatment, and $\tau >0$ is a fixed constant, converges weakly to a mean 0, tight Gaussian process with a covariance function $V^{\left(a\right)}(t,t^{\prime })$ explicitly given in the Supplementary Material of Wang et al. (2020). They also showed that $V^{\left(a\right)}(t,t)={\tilde{V}}^{\left(a\right)}(t,t)-\frac{U\left(t\right)}{{\pi }_1{\pi }_2},$ where ${\tilde{V}}^{\left(a\right)}(t,t^{\prime })$ is the covariance function under simple randomisation and $U\left(t\right)>0$ under stratified permuted block or biased coin randomisation. Again, we see the common phenomenon of reducing variance by applying covariate-adaptive randomisation compared with simple randomisation.

As commented by Wang et al. (2020), the result can be extended to the estimation of survival function with adjusted baseline covariates. Alternatively, one may consider a stratified version of Kaplan–Meier estimator along with the idea in formula (9) of my review.

3. The discussion by Dr. Liu

Dr. Liu's discussion provides useful details and references about asymptotic validity and efficiency of model-assisted inference procedures after covariate-adaptive randomisation and adjustment for covariates. The discussion about the efficiency gain in using ANOVA versus ANCOVA or ANCOVA with treatment-by-covariate interactions started as early as Yang and Tsiatis (2001), continued later by Freedman (2008), Lin (2013) and Wang et al. (2019), and studied under covariate-adaptive randomisation recently by Bugni et al. (2018), Bugni et al. (2019), Liu and Yang (2020), Ma et al. (2020b), Wang et al. (2020) and Ye et al. (in press).

I would like to emphasise two points here. The first one is, as pointed out by Dr. Liu, when there are only two treatment arms and equal allocation is used (k = 2 and ${\pi }_1={\pi }_2=1/2$), the use of ANCOVA with or without treatment-by-covariate interaction has the same asymptotic efficiency and is guaranteed to be more efficient than the use of ANOVA without adjusting for covariates. However, this phenomenon no longer exists once there are more than two treatment arms even if equal allocation is applied.

The second point is that, asymptotically, the most efficient estimator of the treatment difference $\theta =E(Y^{\left(a\right)}-Y^{\left(b\right)})$ defined in the beginning of Section 5.2 of my review is the ${\hat{\theta }}_A$ defined in Section 6.1 of my review, which adjusts for covariates by using a working linear model and the ordinary least squares estimator $\begin{array}{rl}{\widehat{\text{β}}}_a\left(z\right)& ={\left[\sum _{i\in L_a\left(z\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z\left)\right\}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z){\}}^{\mathrm{T}}\right]}^{-1}\\ & \times \sum _{i\in L_a\left(z\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z\left)\right\}{Y}_i\end{array}$ of covariate effect within each stratum $L_a\left(z\right)$ under treatment a and $\mathit{Z}=z$. As Dr. Liu pointed out, however, when there are small strata formed by levels of $\mathit{Z}$, the stratum-specific least squares estimator ${\widehat{\text{β}}}_a\left(z\right)$ might lead to inferior performance due to over-fitting (Liu & Yang, 2020). One modification is to combine ${\widehat{\text{β}}}_a\left(z\right)$ and ${\widehat{\text{β}}}_b\left(z\right)$ within each stratum level z, although they may estimate different quantities. Alternatively, utilising the fact that baseline covariates ${\mathit{U}}_i$'s have the same distribution over all treatment arms, Ye et al. (in press) recommended to replace the matrix inverse in ${\widehat{\text{β}}}_a\left(z\right)$ by an average over all treatment arms to remedy the issue of small strata, which leads to replace ${\widehat{\text{β}}}_a\left(z\right)$ by $\begin{array}{rl}& {\widetilde{\text{β}}}_a\left(z\right)\\ & ={\left[\frac{1}{n\left(z\right)}\sum _{a=1}^k\sum _{i\in L_a\left(z\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z\left)\right\}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z){\}}^{\mathrm{T}}\right]}^{-1}\\ & \times \frac{1}{n_a\left(z\right)}\sum _{i\in L_a\left(z\right)}\{{\mathit{U}}_i-{\overline{\mathit{U}}}_a(z\left)\right\}{Y}_i,\end{array}$ where $n_a\left(z\right)$ is the number of units in $L_a\left(z\right)$ and $n\left(z\right)=n_1\left(z\right)+\cdots +n_k\left(z\right)$. Note that stability issues related with dimensionality for not very large data sets are mainly in the inverses of estimated covariance matrices. Hence, using the inverse of an average may largely remedy the issue of small strata. Some simulation results in Ye et al. (in press) show that using ${\widetilde{\text{β}}}_a\left(z\right)$ in the estimation of θ leads to better finite-sample performance compared with using ${\widehat{\text{β}}}_a\left(z\right)$ or combining ${\widehat{\text{β}}}_a\left(z\right)$ and ${\widehat{\text{β}}}_b\left(z\right)$ when they actually estimate different quantities.

Finally, another way to handle many covariates is to apply high-dimensional technique as Dr. Liu commented (Ma et al., 2020a), or to use variable selection.

4. The discussion by Drs. Ye and Yi

In their discussion, Drs. Ye and Yi clearly described the working models behind estimators ${\hat{\theta }}_S$, ${\hat{\theta }}_A$ and ${\hat{\theta }}_B$ in Sections 5.2 and 6.1 of my review. This not only provides explanations about the asymptotic relative efficiencies among ${\hat{\theta }}_S$, ${\hat{\theta }}_A$ and ${\hat{\theta }}_B$, but also leads to a general working model (formula (1) in the discussion) that produces a class of model-assisted estimators of θ (formula (2) in the discussion) including ${\hat{\theta }}_S$, ${\hat{\theta }}_A$ and ${\hat{\theta }}_B$ as special cases.

In the beginning of Section 6.1 of my review, $\mathit{X}$ is considered to be the vector of all available baseline covariates, $\mathit{Z}$ is the discrete baseline covariate vector (part of $\mathit{X}$) used in covariate-adaptive randomisation, and $\mathit{U}$ is the vector of covariates not in $\mathit{Z}$ but in $\mathit{X}$ to be adjusted for efficiency in the analysis stage. I would like to point out that $\mathit{U}$ may contain some components which are interactions between $\mathit{Z}$ and covariates not in $\mathit{Z}$. Drs. Ye and Yi's discussion classifies $(\mathit{Z},\mathit{U})$ into two categories or vectors, $\mathit{W}$ and $\mathit{V}$, where $\mathit{W}$ contains covariates having treatment-by-covariate interaction in working model (1) in their discussion and $\mathit{V}$ has no treatment-by-covariate interaction. Note that either $\mathit{W}$ or $\mathit{V}$ could be empty. For example, for ANOVA without using any covariate, both $\mathit{W}$ and $\mathit{V}$ are empty; for classical ANCOVA without considering any treatment-by-covariate interaction, $\mathit{W}$ is empty but $\mathit{V}$ is not; as discussed by Drs. Ye and Yi, ${\hat{\theta }}_S$ in Section 5.2 of my review corresponds to $\mathit{W}=\mathit{Z}$ and empty $\mathit{V}$, ${\hat{\theta }}_A$ in Section 6.1 corresponds to $\mathit{W}=(\mathit{Z},\mathit{U})$ and empty $\mathit{V}$, and ${\hat{\theta }}_B$ in Section 6.1 corresponds to $\mathit{W}=\mathit{Z}$ and $\mathit{V}=\mathit{U}$.

In applications, a crucial question is, what is the minimum set of covariates to be included in $\mathit{W}$ or $\mathit{V}$ to ensure that the resulting model-assisted estimator of θ is asymptotically normal with mean θ and variance invariant to the covariate-adaptive randomisation schemes (including Pocock and Simon's minimisation)? As pointed out by Drs. Ye and Yi, a simple answer is that $\mathit{W}$ should contain the dummy variables for all joint levels of $\mathit{Z}$, and there is no requirement on $\mathit{V}$. In fact, $\mathit{V}$ is used to not let the dimension of $\mathit{W}$ become too high. Asymptotically, the estimator with $\mathit{V}$ being empty is most efficient, unless some components of $\mathit{W}$ are actually not related with the response. We must balance between adjusting for covariates and over-fitting, for which variable selection may be a useful solution.

Disclosure statement

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

Additional information

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.