Regression-Based Multiple Treatment Effect Estimation under Covariate-Adaptive Randomization

Abstract

Covariate-adaptive randomization methods are widely used in clinical trials to balance baseline covariates. Recent studies have shown the validity of using regression-based estimators for treatment effects without imposing functional form requirements on the true data generation model. These studies have had limitations in certain scenarios; for example, in the case of multiple treatment groups, these studies did not consider additional covariates or assumed that the allocation ratios were the same across strata. To address these limitations, we develop a stratum-common estimator and a stratum-specific estimator under multiple treatments. We derive the asymptotic behaviors of these estimators and propose consistent nonparametric estimators for asymptotic variances. To determine their efficiency, we compare the estimators with the stratified difference-in-means estimator as the benchmark. We find that the stratum-specific estimator guarantees efficiency gains, regardless of whether the allocation ratios across strata are the same or different. Our conclusions were also validated by simulation studies and a real clinical trial example.

1 Introduction

Randomization is considered to be a gold standard in clinical trials and other intervention studies, such as online A/B tests and experiments in economics. As one of the basic randomization methods, simple randomization allocates experimental units into different treatment groups with a fixed probability. However, simple randomization may cause an imbalance in baseline covariates. In a regulatory guidance, European Medicines Agency (2015) recommended an adjustment for baseline covariates in randomized clinical trials. Covariate-adaptive randomization is arguably the most commonly used method to adjust covariate imbalance in the design stage of randomized clinical trials.

Covariate-adaptive randomization aims to balance treatment allocations among baseline covariates. For example, stratified block randomization (Zelen, 1974) defines sets of strata based on covariates and allocates units in each stratum by using block randomization. It is one of the most commonly used randomization methods in clinical trials, as according to recent surveys, it is implemented in nearly 70% of trials (Ciolino et al., 2019; Lin et al., 2015). Recently, stratified block randomization has also been used in COVID-19 vaccine or treatment trials (e.g., Baden et al., 2021; Wang et al., 2020). Minimization (Pocock & Simon, 1975; Taves, 1974) has been used for balancing covariates over their margins. This scheme has been generalized to control other types of imbalance measures, which may include overall and within-stratum imbalance measures (Hu & Hu, 2012; Hu et al., 2023). Details of above randomization methods and other methods such as stratified biased coin design (Efron, 1971; Shao et al., 2010) and model-based approaches (Atkinson, 1982; Begg & Iglewicz, 1980) are described in Rosenberger and Lachin (2015).

A critical challenge after randomization is to obtain a valid statistical inference for treatment effects. One of the practical strategies applied to address this challenge is the use of linear regression. The robustness of regression estimators under simple randomization has been discussed in several significant papers (Freedman, 2008; Lin, 2013; Yang & Tsiatis, 2001). However, under covariate-adaptive randomization, the challenge becomes complex. Various valid tests have been proposed for different allocation methods under model assumptions of the true data generation process (e.g., Ma et al., 2015, 2020; Wang & Ma, 2021; Shao et al., 2010). A model-assisted approach with regression adjustment for stratification covariates has recently been proposed in Bugni et al. (2018, 2019). This approach produces a valid inference by using a working model between responses and covariates, regardless of whether the working model is correct or not. To consider the regression adjustment for baseline covariates in addition to stratification covariates, stratum-common estimators and stratum-specific estimators have been developed, mainly for the case in which the allocation ratios are the same across strata (Liu et al., 2023; Ma et al., 2022; Ye et al., 2022a, 2022b). However, little attention has been paid to the case in which the allocation ratios are different across strata, especially when additional baseline covariates are included, although different allocation ratios are commonly used in practice and are more flexible (Angrist et al., 2014; Chong et al., 2016).

To fill this gap, we develop stratum-common and stratum-specific estimators under multiple treatments when allocation ratios vary across strata. We justify the asymptotic behaviors of these two estimators and propose nonparametric consistent estimators for the asymptotic variances. Moreover, we demonstrate that for any randomization scheme fulfilling our assumptions, the stratum-specific estimator can guarantee efficiency gains over the other estimators that are considered in this study. The results of our study are the most general obtained so far, with respect to regression adjustment for covariate-adaptive randomization. That is, our methods can handle multiple treatments, adjust for additional baseline covariates to improve efficiency, and allow different allocation ratios across strata.

The remainder of the paper is organized as follows. In Section 2, we introduce the framework and assumptions for covariate-adaptive randomization. In Section 3, after a brief review of the benchmark estimator, we describe the stratum-common estimator, the stratum-specific estimator, and their asymptotic properties. We present the simulation studies in Section 4. Section 5 provides a real example of a clinical trial and Section 6 consists of our recommendations and a discussion. All the estimators considered in this paper are implemented in the R package caratMULT, available at GitHub.

2 Framework and Assumptions

Consider a covariate-adaptive randomization procedure with n units. Let us suppose that treatments ⁠. Let be the control group and ⁠. Let be the indicators, such that indicates that the ith unit is assigned to the ath treatment. Let be the number of units in each treatment group. Let denote the stratum label, which takes value in ⁠. Here, K represents the total number of strata, which is fixed and finite. Let be the p-dimensional additional baseline covariates that are not used in the randomization procedure. Let be the target proportion of stratum k and be the target proportion of treatment a in stratum k, so that is the target proportion of treatment a. Let be the number of units in stratum k and be the number of units in stratum k and treatment group a. Let be the estimated proportion of stratum k.

We use the Neyman–Rubin model to define potential outcomes and treatment effects (Neyman, 1923; Rubin, 1974). Let {} be the potential outcomes and be the observed outcomes. Our goal is to estimate the treatment effect for all ⁠. That is, we will estimate based on the observed data ⁠. The sample means of the potential outcomes are defined as and ⁠. The population variance of a transformed outcome ⁠, such as ⁠, is denoted as ⁠. We define as the population variance of a transformed outcome in stratum k.

The following assumptions are made for the data generation and covariate-adaptive randomization procedures. We denote the sets of random variables with at least one positive stratum-specific variance as ⁠. Let be independent and identically distributed samples from the population distribution ⁠.

Assumption 1. and ⁠, for all ⁠.

Assumption 2. Conditional on ⁠, is independent with ⁠.

Assumption 3. as ⁠, for all and ⁠.

According to Assumption 2, given the strata, the potential outcomes and the additional covariates are conditionally independent of the assignment procedure. Assumption 3 is the same as Assumption 2.2(b) in Bugni et al. (2019). Several well-known randomization schemes satisfy this assumption, such as simple randomization, stratified permuted block randomization (Zelen, 1974) and stratified biased coin randomization (Kuznetsova & Johnson, 2017; Shao et al., 2010). In the special case in which are identical across all strata for each treatment group, Pocock and Simon's minimization (Pocock & Simon, 1975) also satisfies this assumption (Hu et al., 2023).

3 Regression-based Multiple Treatment Effect Estimation

3.1 Benchmark Estimator

First, we introduce the benchmark estimator—an estimator for without using additional covariates. This estimator was proposed in Bugni et al. (2019) and was obtained by considering the estimation of the linear regression

by the ordinary least squares (OLS) method. We denote as the OLS estimator of ⁠. From calculations,

Note that we use the estimated proportion rather than the target proportion in our estimators. This estimator can also be considered as a plug-in estimator (Liu & Yang, 2020). We define

as the treatment effect in stratum k. Therefore, the treatment effect can be obtained from

It was shown that is an unbiased and consistent estimator of under Assumptions 1–3, and the plug-in estimator is a consistent estimator of (Ma et al., 2022). Let be the estimator of τ. The asymptotic behavior of was also derived in Bugni et al. (2019), implying that

We provide the expanded form of later in the proof in the Web Appendix A. We seek to determine whether adjusting the imbalance of the additional covariates will improve the efficiency compared with ⁠.

Remark 1. Regression (1) is not the original regression used in Bugni et al. (2019); however, both regressions yield the same results when estimating ⁠. In our formulation, the estimator for can be directly obtained. Regression (1) is an extension of the third regression in Ma et al. (2022) to multiple treatments.

In addition to stratification indicators, the covariates may contain additional information that can improve estimation of the treatment effect. Regression is a common strategy for adjusting additional covariate imbalance. Below, we discuss two regression models with the inclusion of ⁠. Any covariates that can be linearly represented by stratification indicators should be removed from the additional covariates. For simplicity, we continue to use as the additional covariates without loss of generality.

Before formally specifying the models and the asymptotic results, we define several population-level regression coefficients. Let be the covariance matrix for any two random vectors R and Q, be the covariance matrix for R and Q given stratum k, and ⁠. We assume that X has a finite second moment and ⁠, ⁠, and are positive definite. Two types of coefficients are used for additional covariates: stratum-common coefficients and stratum-specific coefficients ⁠, which are defined as

Here, are the population-level weighted regression coefficients for regressing on with weights ⁠, and are the population-level regression coefficients for regressing on in stratum k. We also define some overall and stratum-specific sample-level means: ⁠, ⁠, ⁠; and some treatment-specific and stratum-specific sample-analog covariance matrices: ⁠, ⁠, ⁠, ⁠.

3.2 Stratum-common Estimator

We first consider the stratum-common estimator. The regression is

Let be the OLS estimator for τ and be the ath element of ⁠. As shown in our proof, regression (2) is equivalent to running regressions separately in the treatment and the control groups. Then, we obtain

where and are the estimators for the regression coefficients.

Theorem 1. Suppose that Assumptions 1–3 hold. Let ⁠. Then,

where⁠,

Here, 1_A is an A-dimensional vector of ones and .

Remark 2. The asymptotic result is not restricted to the regression coefficients ⁠. For any population coefficients ⁠, if is a consistent estimator for ⁠, then retains this asymptotic behavior. In the special case in which for all ⁠, the benchmark estimator and its asymptotic behavior can be obtained.

An area of concern is whether the asymptotic variance will decrease when the additional covariates are added. Results from Ma et al. (2022) show that not all regression models can achieve efficiency gain due to the additional covariates. In general, the stratum-common estimator may have a more adverse impact on the efficiency than ⁠. In the special case in which for all ⁠, can improve the efficiency. Let denote that matrix V is negative semi-definite.

Corollary 1. When for all and ⁠, we obtain

Remark 3. When there is only one treatment group and one control group, our model reduces to the case in Ma et al. (2022), the regression model with interaction and additional covariates. Furthermore, our estimator and the asymptotic variance are identical to those derived in Ma et al. (2022) for this particular case.

Remark 4. The use of stratum-common coefficients was recommended by Ye et al. (2022a). In their work, the coefficients were achieved by the analysis of covariance using a heterogeneous working model (ANHECOVA) and the estimator derived by these coefficients for potential outcomes was referred to ANHECOVA estimator. With respect to the considerations in their study, the ANHECOVA estimator achieved a guaranteed efficiency gain over other estimators that they considered, not only under simple randomization but also under covariate-adaptive randomization. In general, when different allocation ratios are used across strata, the stratum-common estimator may adversely affect the efficiency. As shown in Corollary 1, when for all strata, the stratum-common estimator can improve the efficiency, as the ANHECOVA estimator does. To guarantee efficiency gains when using the stratum-common estimator, different weights can be assigned for each stratum; however, a more efficient estimator can be constructed, as described in the next section.

3.3 Stratum-specific Estimator

To develop a regression model for the stratum-specific estimator, we introduce the following linear model.

for all ⁠. Here, we use to indicate the units in stratum k. Let denote the OLS estimator for ⁠, denote the estimator for ⁠, and denote the estimator for τ. The regression can also be split into regression models, giving

where and are the estimators for the regression coefficients.

Remark 5. We do not provide an overall regression model because of the complicated expression arising from the three-way interaction between additional covariates ⁠. Such expression may not be welcomed for application reasons.

Theorem 2. Suppose that Assumptions 1–3 hold. Let ⁠; then

where with

and has the same definition as in Theorem 1.

Furthermore, is more efficient than and due to the following corollary.

Corollary 2. Suppose that Assumptions 1–3 hold.

In addition,

Remark 6. The use of a stratum-specific estimator under covariate-adaptive randomization has been discussed by Liu and Yang (2020) and Ye et al. (2022b). Liu and Yang (2020) showed that when there are two treatments, a stratum-specific estimator can achieve guaranteed efficiency gains over a stratum-common estimator. We generalize their results to multiple treatments. Ye et al. (2022b) also claimed efficiency gains for a stratum-specific estimator over other estimators when the same allocation proportion was applied across strata. Our results are more general, justifying that our stratum-specific estimator is more efficient than the other estimators considered in our study, regardless of whether the target allocation ratios across strata are the same or different.

3.4 Consistent Variance Estimators

For all of the estimators mentioned above, a crucial step when drawing a valid inference is to construct a consistent asymptotic variance estimator. Specifically, we consider the difference in the potential outcomes for any two treatments b and c. Let be an A-dimensional vector, such that all elements are zero, except for the bth element, which is 1, and the cth element, which is −1. Let ⁠, be the real difference in potential outcomes, and be the estimator for generated from the three estimators above, for ⁠. Based on simple calculation,

where

The matrix is the asymptotic variance for ⁠, and denotes the element in the ith row and jth column of ⁠, for ⁠. By using as the transformed outcomes again and as an estimator of ⁠, ⁠. Let ⁠. For the three estimators above, ⁠, where ⁠, respectively, for ⁠. We define

and

Let and denote the overall sample covariance matrix and stratum-specific sample covariance matrix, respectively. Let be the estimator of for ⁠.

Theorem 3. Suppose that Assumptions 1–3 hold. The variance estimators

are consistent estimators of for ⁠, respectively.

4 Simulation Studies

We ran simulations for several models to compare the finite-sample behavior of the proposed estimators. The potential outcomes are generated by the following model:

where is a discrete variable, takes a value from with probabilities (0.2,0.2,0.3,0.3), and is used for stratification such that ⁠. is a continuous variable serving as additional covariates. is a random error. We also ran additional simulations when was not normally distributed. The results are presented in Section B.1 in the Web Appendix. We considered ⁠, that is, two treatment groups and one control group. We estimated two treatment effects, τ₁ and τ₂. The number of units n was 600. Details of the outcome models are as follows:

For coefficients, we considered ⁠, and as an arithmetic sequence with 12 numbers, with 1 being the first number and 6 being the last number, first ordered by strata and then by treatments. Three types of randomization methods were considered: simple randomization, stratified block randomization, and Pocock and Simon's minimization. For Tables 1 and 2, the treatment allocations were the same across strata, with an equal ratio (1:1:1) and an unequal ratio (1:2:2), respectively. For Table 3, we assigned different treatment allocation ratios to different strata. When and ⁠, the allocation ratio was 1:2:2, and when and ⁠, the allocation ratio was 2:3:5. Note that we only considered simple randomization and stratified block randomization for Table 3, because minimization requires the same allocation ratios across strata. For stratified block randomization, the block size was 6 in Table 1 and 10 in Tables 2 and 3. For Pocock and Simon's minimization, range was used as the marginal imbalance metric, and no overall balance was considered. The allocation probabilities were determined by the procedure given in Han et al. (2009); in particular, when the treatment with the smallest allocation ratio was preferred, the probability for this treatment was 0.9. As the allocation ratios may vary step by step under unequal allocation ratio for the procedure in Han et al. (2009) in Table 2, we also considered the minimization procedure proposed by Kuznetsova and Tymofyeyev (2012), which preserves the allocation ratio at every step. The additional simulation results can be found in Section B.2 in the Web Appendix. The bias, standard deviation (SD) of the estimators, standard error (SE) estimators, and empirical coverage probabilities (CP) of 95% confidence intervals were evaluated based on 1,000 simulation runs. The results are summarized below.

• (1)

  Almost all estimators had negligible bias (⁠⁠). The results of Model 3 are an exception, with a maximum bias of approximately 3%.

• (2)

  Most of the proposed variance estimators matched the corresponding simulated variance, and the expected CP of 95% was reached.

• (3)

  In Tables 1 and 2, for all models, was more efficient than ⁠, which supports Corollary 1. Specifically, in Models 1 and 2, reduced the standard errors by approximately 40%–50%.

• (4)

  For all models we considered, was more efficient than ⁠, as we showed in Corollary 2. The efficiency gain was greater when was used than when was used. As expected, in Model 1, the standard errors were similar for and ⁠. In Model 2, led to an approximately 10% reduction in standard errors compared with ⁠. These results also confirm that the efficiency gain is guaranteed, regardless of whether the target proportions across strata are the same or different.

• (5)

  For two linear models (Models 1 and 2), estimators with additional covariates considerably improved the efficiency compared with ⁠. However, the improvement is not apparent in Model 3, which is nonlinear and may not benefit most from the regression-based adjustment.

5 Clinical Trial Example

In this section, we first analyze real data from a clinical trial that compared nefazodone, the cognitive behavioral analysis system of psychotherapy (CBASP), and their combination for the treatment of chronic depression (Keller et al., 2000). We consider the combination as the control group indexed as 0, nefazodone as the treatment indexed as 1, and the CBASP as the treatment indexed as 2. The trial randomly assigned 681 patients to three treatments with an allocation ratio of 1:1:1, and the outcome of interest was the final score of the 24-item Hamilton rating scale for depression. To evaluate the performance of our proposed estimators, we consider five covariates in the estimation: AGE, HAMD17, HAMD24, HAMD_COGNID, and GENDER, which are detailed in Table 4. The data are stratified using GENDER and HAMD17, where HAMD17 is discretized into two levels: less than 18 and greater than or equal to 18. AGE, ⁠, HAMD17, HAMD24, and HAMD_COGNID are treated as additional covariates. Table 5 presents the estimates along with the 95% confidence intervals and the variance reductions. The results show that our proposed estimators reduce the variances by 5.83% to 10.16%. Note that the proposed methods lead to slightly increased estimates. The increase is mainly due to the adjustment of the imbalance in HAMD17, which is highly predictive of the outcomes.

Furthermore, to illustrate the effectiveness of proposed estimators under more general settings of allocation ratio, we also analyze the synthetic data generated based on the trial. To be specific, we use simple randomization or stratified block randomization to allocate the patients with different allocation ratios. The unobserved potential outcomes are imputed by fitting the real trial data with a nonparametric spline using the function bigssa in the R package bigspline. The covariates we choose to fit the data are the same with the covariates in Table 4. The treatment effects and variances are estimated based on 500 replicates. The results are presented in Table 6 and show that our proposed estimators have a variance reduction compared to the benchmark estimator. The stratum-specific estimator reduces variance by 9.07% to 13.26%, and the stratum-common estimator reduces variance by 4.96% to 7.97%. The efficiency gain of the proposed estimators is observed under both simple randomization and stratified block randomization.

6 Discussion

In covariate-adaptive randomization, efficiency could be improved by adjusting baseline covariates. In this study, we consider stratum-common and stratum-specific treatment effect estimators for general scenarios, including multiple treatments and different allocation ratios across strata. We also develop consistent estimators for asymptotic variance. We found that when compared with the benchmark estimator and the stratum-common estimator, the stratum-specific estimator can ensure efficiency gains for all cases under our settings. Therefore, we recommend the use of the stratum-specific estimator when certain conditions hold.

A vital condition for the validity of the stratum-specific estimator is that the coefficient estimators obtained from regression must be consistent, which requires a considerable number of units in each stratum. This condition may not be satisfied in clinical trials that have many small strata. Future work could focus on addressing such issues. In addition, future studies could explore the extension of our results to generalized linear models, such as logistic regression for binary outcomes or Poisson regression for count data (Guo & Basse, 2023). Finally, another potential area of interest is a generalization to the high-dimensional settings under multiple treatments.

Data Availability Statement

The data that support the findings in this paper are available in the Supporting Information of this paper.

Supporting Information

Web Appendices A and B, referenced separately in Sections 3 and 4, and R codes to reproduce the simulation results are available with this paper at the Biometrics website on Wiley Online Library. An R package caratMULT, which implements all the estimators in this paper, is available at GitHub (https://github.com/guyujia98/caratMULT).

Acknowledgments

This work was supported by National Natural Science Foundation of China (Grant Nos. 12171476, 12071242).

References