Balancing Treatment Allocation over Continuous Covariates: A New Imbalance Measure for Minimization

Abstract

In many clinical trials, it is important to balance treatment allocation over covariates. Although a great many papers have been published on balancing over discrete covariates, the procedures for continuous covariates have been less well studied. Traditionally, a continuous covariate usually needs to be transformed to a discrete one by splitting its range into several categories. Such practice may lead to loss of information and is susceptible to misspecification of covariate distribution. The more recent papers seek to define an imbalance measure that preserves the nature of continuous covariates and set the allocation rule in order to minimize that measure. We propose a new design, which defines the imbalance measure by the maximum assignment difference when all possible divisions of the covariate range are considered. This measure depends only on ranks of the covariate values and is therefore free of covariate distribution. In addition, we developed an efficient algorithm to implement the new procedure. By simulation studies we show that the new procedure is able to keep good balance properties in comparison with other popular designs.

1. Introduction

Balanced allocation among treatment groups is often desirable in many clinical trials. A well-balanced design enhances the trials by increasing the credibility of the trial, precision of subgroup or interim analysis, and robustness to model misspecification [1]. Moreover, ignorance of balance at the design stage may lead to loss of statistical efficiency, especially in small trials [2]. In the existence of covariates (or prognostic factors), well-balanced allocation does not only mean similar group sizes but also similar distributions of covariate values across treatment groups [3–6]. With discrete covariates, a great many papers have been published, which include stratified permuted block designs, marginal procedures (or minimization) [7, 8], and hierarchical models [9–11]. Y. Hu and F. Hu [12] are among the few who explored the theoretical properties of such procedures.

With respect to balancing over continuous covariates, a traditional way is to discretize a continuous covariate by splitting its range into several small intervals, each of which defines a distinct category [4, 13]. However, such methods have several potential problems. First, there is no standard way of defining such intervals: usually the boundaries are set by experts’ opinion from a nonstatistical perspective, or they are heavily dependent on the model assumption of the covariate distribution. Weir and Lees [13] studied the randomization procedure in an acute stroke clinical trial in which 12 covariates were involved. The five continuous covariates were discretized before the traditional minimization was applied. Based on the historical information of the covariate distribution, the covariates age and mean arterial pressure were discretized by their quintiles; plasma glucose level was first discretized into two groups by the critical value 8 mmol/L and the one with larger values was further split by quartiles; the remaining two covariates, time of delay from stroke onset and Glasgow coma scale, were categorized whose critical values were set by clinical experience. Therefore, we see that discretization methods are rather covariate-specific and may depend on the distribution of the covariates. Second, there is a dilemma as to the selection of the number of categories (denoted by 𝑚). While a small 𝑚 may lead to loss of information, a large 𝑚 can result in a severe overall imbalance, as is in the case of stratified designs when the number of strata is large [7].

An alternative that saves the step of discretization is to define a certain “imbalance measure” for continuous covariates among different treatment groups and then allocate the treatments sequentially to minimize that specific measure [5, 6, 14–16]. With two treatment groups 𝐴 and 𝐵, many authors defined imbalance measures as the mean difference of covariate values between the two groups: some of the measures may further be adjusted for standard deviations within the two groups [5, 14], others normalized the mean difference by calculating the two-sample 𝑡-statistic and the corresponding 𝑃 value [15], that is, a larger 𝑃 value indicates more balanced allocation. An alternative to the above mean comparison methods is model-based approach [17–20], which relates different allocation to different design matrices in a linear model, which in turn influence the maximum likelihood estimator of the treatment effect 𝛼_𝑛; it then sequentially favors an assignment that would lead to a smaller var(𝛼_𝑛). This approach only implies balance under a homogeneous linear model.

The methods discussed above, however, usually depend on the assumption of normal distributions of the covariates. Among the methods that apply to more general distributions, Stigsby and Taves [16] considered the difference of rank sums as their imbalance measure; Su [6] used a weighted average of the difference in overall patient numbers and the difference in quartiles of the continuous covariates. Rosenberger and Sverdlov [21] found the empirical distributions of the covariates in the two groups and calculated their Kolmogorov-Smirnov distance to measure the imbalance for a series of randomization procedures.

In this paper, we propose a new minimization procedure which carefully examines the assignment differences over all intervals and set their maximum value as the actual imbalance measure. The new procedure utilizes only the rank information of the covariates and is therefore robust against various distributions. Moreover, we developed an efficient algorithm whose computation is even less intensive than discretization. By simulation studies, we show that our new procedure leads to well-balanced allocation in terms of overall patient numbers and covariate distributions. For simplicity, we first consider only one continuous covariate and then extend the idea to more general cases.

In Section 2, by giving a motivating example, we analyze the downside of simple discretization and then highlight the idea underlying the new procedure. In Section 3, details of the procedure are described, and the key part of the algorithm is explained. Simulation studies are carried out in Section 4, comparing our procedure with discretization methods as well as some other designs. In Section 5, we conclude the paper by discussing possible extension of the current procedure.

2. Motivation

We will first focus on two treatments, 𝐴 and 𝐵 and only one continuous covariate 𝑍. For simplicity, assume 𝑍 takes value over interval [0,1]; any other type of 𝑍 can be transformed to [0,1] by a linear function if 𝑍 is bounded, or a nonlinear function such as 𝑒^𝑥/(𝑒^𝑥+1) otherwise.

To see the difficulty in defining “balance’’ for a continuous covariate, we use an example as shown in Figure 1. The first 8 patients have been randomized. The 9th has just arrived, who is tentatively assigned to treatment 𝐴 (a) and 𝐵 (b). The figure shows the assignment as well as the covariate values of the patients. First, we point out that discretization could cause complications since different ways of splitting may lead to different assignments. For example, if the covariate range [0,1] in Figure 1 is split into 𝑚 intervals of equal lengths, then the choices of 𝑚=1, 𝑚=2, and 𝑚=4 result in different preferences of the assignment: if 𝑚=1, that is, only overall patient numbers are considered, then the 9th patient could be assigned to treatment 𝐴 or treatment 𝐵, since in either way the absolute difference is 1; if 𝑚=2, then treatment 𝐴 would be favored, producing balanced patient numbers in the second category, that is, over the interval [0.5,1.0] where the 9th patient belongs; if 𝑚=4, then treatment 𝐵 would be preferred instead, since it produces allocation of 2 : 1 in the third category, that is, over the interval [0.5,0.75), rather than 3 : 0 when treatment 𝐴 is assigned. Second, we emphasize that balance in the mean values does not necessarily lead to balance in distributions. Taking the upper panel in Figure 1, for instance, the assignment of treatment 𝐴 would not cause a significant mean difference in the two treatment groups. However, in terms of covariate distribution, a severe imbalance would exist, since the 5 patients (7th, 6th, 9th, 4th, and 8th) in group 𝐴 have covariate values at the center of the range and the 4 patients (2nd, 3rd, 5th, and 1st) in group 𝐵 at the two ends.

(a)

(b)

Figure 1 

Two allocation sequences for the same covariate values.

In light of the above discussion, we propose a new imbalance measure, which is defined as the maximum (absolute) assignment difference when all possible ways of discretization are considered. By doing so, we want to make sure that the assignment difference over any interval does not become too extreme. Thus, if the 9th patient is tentatively assigned to 𝐴, we need to examine the differences of patient numbers in the two groups in all neighborhoods of 𝑍₉ and use the maximum of them as the resulting degree of imbalance. The corresponding imbalance caused by the assignment of 𝐵 can be calculated in the same manner. Then with a probability greater than 1/2, the 9th patient should be assigned to a treatment that would cause a smaller maximum difference.

3. Procedure

Consider a sequential trial with two treatments 𝐴 and 𝐵 (control and test). For the first 𝑛 patients that have arrived, let 𝐓_𝑛=(𝑇₁,…,𝑇_𝑛) be the allocation sequence with 𝑇_𝑖=𝐴 if the 𝑖th patient is assigned to treatment 𝐴 and 𝑇_𝑖=𝐵 otherwise. Suppose we need to balance treatment allocation over one single continuous covariate 𝑍 with cumulative distribution function 𝐹, whose range 𝑈 is an interval on the real line. The two endpoints of 𝑈 can be finite or infinite. Let 𝐙_𝑛=(𝑍₁,…,𝑍_𝑛), where 𝑍_𝑖 is the covariate value of the 𝑖th patient. Suppose 𝑍₁, …, 𝑍_𝑛 are independent. Define Δ𝐼,𝐙_𝑛,𝐓_𝑛=𝑁_𝐴𝐼,𝐙_𝑛,𝐓_𝑛−𝑁_𝐵𝐼,𝐙_𝑛,𝐓_𝑛(3.1) as the difference of patient numbers in treatment groups 𝐴 and 𝐵 over the interval 𝐼⊂𝑈, given the covariate information 𝐙_𝑛 and allocation sequence 𝐓_𝑛. Suppose the first 𝑛 patients have been randomized and the (𝑛+1)th patient has just arrived, that is, 𝐙_𝑛+1=(𝑍₁,…,𝑍_𝑛+1) and 𝐓_𝑛=(𝑇₁,…,𝑇_𝑛) are known and 𝑇_𝑛+1 is to be determined. Define 𝐷^(𝑘)𝑛+1=sup_𝐼||Δ𝐼,𝐙_𝑛+1,𝐓_𝑛+1||∶𝑍_𝑛+1∈𝐼⊂𝑈,𝑇_𝑛+1=𝑘,(3.2) where 𝑘=𝐴 or 𝐵. Thus, 𝐷^(𝑘)𝑛+1 is the potential maximum absolute difference it would cause if the new patient was assigned to treatment 𝑘. Note that the interval 𝐼 in (3.2) needs to contain the new covariate value 𝑍_𝑛+1, because the difference over any other interval will not be affected by the the arrival of 𝑍_𝑛+1 and is therefore not of our interest. Note also that 𝐷^(𝑘)𝑛+1 is a function of (𝐙_𝑛+1,𝐓_𝑛,𝑘).

Then, assign the (𝑛+1)th patient to treatment 𝐴 with the following probability: 𝑃𝑇_𝑛+1=𝐴∣𝐙_𝑛+1,𝐓_𝑛=⎧⎪⎨⎪⎩𝑝,if𝐷^(𝐴)𝑛+1𝑞,if𝐷_(𝐴)𝑛+1<𝐷^(𝐵)𝑛+1,0.5,otherwis>e,𝐷^(𝐵)𝑛+1,(3.3) where 0<𝑞<𝑝<1 and 𝑝+𝑞=1.

We will show how the 9th patient is randomized according to our new procedure. The critical part lies in the calculation of 𝐷^9(𝐴) and 𝐷^9(𝐵). For the former, that is, the 9th patient is temporarily assigned to treatment 𝐴 (see the upper panel in Figure 1), then the maximum absolute difference is 5, attained over intervals which exclusively contain 𝑍₇, 𝑍₆, 𝑍₉, 𝑍₄, and 𝑍₈, that is, 5 patients over these intervals are assigned to 𝐴 and 0 is to 𝐵, so 𝐷^9(𝐴)=5. Similarly, 𝐷^9(𝐵)=3. Therefore, since 𝐷^9(𝐴)>𝐷^9(𝐵), with probability 𝑝>1/2 the 9th patient will be assigned to treatment 𝐵.

We would like to point out that in order to calculate the imbalance 𝐷^(𝑘)𝑛+1, which is defined as a supremum, it is sufficient to examine intervals whose endpoints belong to the set {𝑍₁,…,𝑍_𝑛+1}. Hence the total number of such intervals has the order 𝑂(𝑛²). Moreover, the difference of patient numbers over any of these intervals is only related to the ranks of {𝑍₁,…,𝑍_𝑛+1}, whose joint distribution places an equal probability 1/[(𝑛+1)!] on any permutation of [1,2,…,(𝑛+1)]. To support the above argument, we can reexamine the upper panel in Figure 1: for any interval 𝐼=[𝑎,𝑏], where 𝑍₂<𝑎≤𝑍₇ and 𝑍₅≤𝑏≤𝑍₁, the difference of patient numbers over 𝐼 is exactly the same as that over interval [𝑍₇,𝑍₅], which is +3=5−2; furthermore, so long as the relative positions of {𝑍₁,…,𝑍_𝑛+1} remain the same, this difference of +3 does not change. Nor does 𝐷^(𝑘)𝑛+1. Therefore, we come to the conclusion that the new procedure is free of the underlying distribution 𝐹.

In fact, the computation time of 𝐷^(𝑘)𝑛+1 can be reduced by examining an even smaller number of intervals, that is, (𝑛+2) instead of 𝑂(𝑛²). Before demonstrating this, we need a few more notations and definitions. Since our new procedure is distribution-free, simply assume that the covariate 𝑍 is from uniform [0,1]. Suppose 𝐙_𝑛+1 and 𝐓_𝑛 have been observed. Δ(𝐼,𝐙_𝑛,𝐓_𝑛) in (3.1), defined as difference of patient numbers in groups 𝐴 and 𝐵 over interval 𝐼, will simply be written as Δ_𝑛𝐼. For the ease of representation, let 𝑍_𝐿=0 and 𝑍_𝑅=1. Define two sets 𝒮_𝐿 and 𝒮_𝑅 as 𝒮𝒮_𝐿=Δ_𝑛𝑍_𝑖,𝑍_𝑛+1∶𝑍_𝑖≤𝑍_𝑛+1,𝑖=𝐿,1,…,𝑛_𝑅,=Δ_𝑛𝑍_𝑛+1,𝑍_𝑗∶𝑍_𝑛+1≤𝑍_𝑗,𝑗=𝑅,1,…,𝑛.(3.4) That is, 𝑍_𝑖 is any point from the set {0,𝑍₁,…,𝑍_𝑛}, that is, to the left of or equal to 𝑍_𝑛+1, and [𝑍_𝑖,𝑍_𝑛+1) is a left-closed and right-open interval. For instance, in Figure 1,  𝑛=8 and 𝒮_𝐿={0,𝑍₂,𝑍₇,𝑍₆} from left to right. The interpretation of 𝑍_𝑗 and (𝑍_𝑛+1,𝑍_𝑗] is similar.

Proposition 3.1. Let 𝐶_𝐿1=min𝒮_𝐿, 𝐶_𝐿2=max𝒮_𝐿, 𝐶_𝑅1=min𝒮_𝑅, and 𝐶_𝑅2=max𝒮_𝑅. Then, 𝐷𝐷^(𝐴)𝑛+1=max||𝐶_𝐿1+𝐶_𝑅1+1||,||𝐶_𝐿1+𝐶_𝑅2+1||,||𝐶_𝐿2+𝐶_𝑅1+1||,||𝐶_𝐿2+𝐶_𝑅2+1_(𝐵)𝑛+1||,=max||𝐶_𝐿1+𝐶_𝑅1−1||,||𝐶_𝐿1+𝐶_𝑅2−1||,||𝐶_𝐿2+𝐶_𝑅1−1||,||𝐶_𝐿2+𝐶_𝑅2−1||.(3.5)

The proof of Proposition 3.1 is given in the Appendix. Proposition 3.1, together with the definitions of 𝒮_𝐿 and 𝒮_𝑅, suggests the following.(1)The calculation of 𝐷^(𝐴)𝑛+1 and 𝐷^(𝐵)𝑛+1 only requires the examination of (𝑛+2) intervals, instead of 𝑂(𝑛²).(2)For two consecutive intervals [𝑍_𝑖,𝑍_𝑛+1) and [𝑍_𝑖,𝑍_𝑛+1) in 𝒮_𝐿, where 𝑍_𝑖<𝑍_𝑖, 𝑖≠𝐿 and 𝑖^≠𝐿, we have Δ_𝑛[𝑍_𝑖,𝑍_𝑛+1)=Δ_𝑛[𝑍_𝑖,𝑍_𝑛+1)±1, depending on the assignment 𝐴 or 𝐵 for the patient at 𝑍_𝑖 (the same argument applies to intervals in 𝒮_𝑅).

The above two observations form the basis of the algorithm, which was developed for the new procedure in the simulation studies (Section 4). We found that the computation time of the new procedure was even less than discretization methods.

4. Simulation Studies

Suppose 𝑁=60 patients enrolled. As mentioned in Section 1, for continuous covariates, it is desirable to keep similarity between treatment groups in two aspects: the group sizes and the distributions of the covariates. Therefore, the new procedure was first compared with several other procedures in terms of the following two criteria: (1) the mean absolute difference 𝐸|𝑁_𝐴(𝑁)−𝑁_𝐵(𝑁)| of all patient numbers in the two groups, shown as 𝐸𝐷_all; (2) the mean Kolmogorov-Smirnov distance (K-S) between the empirical distributions of covariate 𝑍 in groups 𝐴 and 𝐵, shown as 𝐸𝐷_ks, which basically measures the similarity between two distributions. In addition, we used a new criterion: (3) the “maximum imbalance” defined by us as: sup_𝐼||Δ𝐼,𝐙_𝑁,𝐓_𝑁||∶𝐼⊂[0,1],(4.1) shown as 𝐸D_max, which is the maximum absolute difference over all possible intervals after all patients have been assigned to a treatment. We will show that criterion (3) acts as a compromise between criterion (1) and criterion (2).

Since the procedures we compared are all distribution-free, the independent covariate values 𝑍₁,…,𝑍_𝑁 were simply generated from Unif(0,1). All procedures use the strategy of minimization, but each has a different imbalance measure. More specifically, under a certain imbalance measure 𝐷, we calculate 𝐷^(𝐴) or 𝐷^(𝐵), defined as the imbalance that would occur if the new patient was assigned to treatment 𝐴 or 𝐵. Depending on whether 𝐷^(𝐴)−𝐷^(𝐵) is positive, negative, or zero, the allocation probability toward the treatment 𝐴 is 𝑞, 𝑝, or 1/2, where 0<𝑞<𝑝<1 and 𝑝+𝑞=1. In the simulation, we used 𝑝=2/3 and 𝑝=1, with the latter corresponding to deterministic allocation unless there is a tie.

The following procedures were studied.(1) Efron’s design (EFRON) [22]. Let 𝑁_𝐴 and 𝑁_𝐵 be the patient numbers in the two groups. Define the imbalance measure 𝐷 as |𝑁_𝐴−𝑁_𝐵|. This method solely focuses on the balance of patient numbers.(2) Kolmogorov-Smirnov measure (K-S). Let 𝐹_𝑘 be the empirical distribution function of the covariate values in group 𝑘, 𝑘=𝐴, 𝐵. Define imbalance measure 𝐷 as the Kolmogorov-Smirnov distance between 𝐹_𝐴 and 𝐹_𝐵. This method solely focuses on the balance of distributions. The above two methods are rarely used as a way of balancing over a continuous covariate, since each of them is designed to meet only one criterion. In our simulations, they served as two controls to evaluate other procedures.(3) Discretization (DSCRT). In practice, in order to discretize a continuous covariate 𝑍 with cumulative distribution function 𝐹, the range is often split by the quantiles of 𝐹 at probabilities 1/𝑚,2/𝑚,…,(𝑚−1)/𝑚. This is equivalent to splitting [0,1] into 𝑚 intervals of equal length 1/𝑚 for 𝑍∼ unif[0,1]. In our simulations, we tried 𝑚=2,4,8. Within each category Efron’s design was applied.(4)The new procedure (MAX-IMB).(5)Stigsby and Taves’ rank sum (RANK-SUM) [16]. Let (𝑅₁,…,𝑅_𝑁) be the ranks of (𝑍₁,…,𝑍_𝑁). Suppose patients 𝑖₁,…,𝑖_𝑁1 are in group 𝐴 and patients 𝑗₁,…,𝑗_𝑁2 are in group 𝐵. The imbalance measure 𝐷 is defined by |∑^𝑁1_𝑘=1𝑅_𝑖𝑘−∑^𝑁2_𝑘=1𝑅_𝑗𝑘|.(6) Su’s weighted average (WGT-AVE) [6]. Let 𝑄^𝑘1,𝑄^𝑘2,𝑄^𝑘3 be the quartiles of the covariate values in group 𝑘, 𝑘=𝐴, 𝐵. The “qualitative’’ imbalance measure 𝐷_𝐿 is defined by 𝑤₀𝐼||𝑁_𝐴−𝑁_𝐵||>𝑐₀+𝑤₁𝐼max_𝑖=1,2,3||𝑄^𝐴𝑖−𝑄^𝐵𝑖||max𝑄^𝐴𝑖,𝑄^𝐵𝑖>𝑐¹,(4.2)where 𝑤₀ and 𝑤₁ are two weights placed on the two items, 𝑐₀ and 𝑐₁ are two upper limits, and 𝐼 is the indicator function; the “quantitative’’ imbalance measure 𝐷_𝑇 is defined by 𝑤₀||𝑁_𝐴−𝑁_𝐵||𝑐₀+𝑤¹max_𝑖=1,2,3||𝑄^𝐴𝑖−𝑄^𝐵𝑖||/max𝑄^𝐴𝑖,𝑄^𝐵𝑖𝑐₁.(4.3) Therefore, let 𝐷^𝐿(𝑘) be the the qualitative imbalance resulted from the tentative assignment of treatment 𝑘 to the new patient, 𝑘=𝐴,𝐵. If 𝐷^𝐿(𝐴)−𝐷^𝐿(𝐵) is positive (negative), the allocation probability toward the treatment 𝐴 will be 𝑞(𝑝); if 𝐷^𝐿(𝐴)−𝐷^𝐿(𝐵)=0, use measure 𝐷_𝑇 to determine the probability 𝑝, 𝑞, or 1/2.

𝑤₀, 𝑤₁, 𝑐₀, and 𝑐₁ can be changed freely from a subjective point of view. In the simulations, we fixed 𝑤₀=𝑤₁=1 and 𝑐₁=10%, but tried 𝑐₀=2,4,6.

The results for 𝑝=2/3 and 𝑝=1 under 5000 repetitions are shown in Tables 1 and 2.

Table 1 

Comparison of 𝐸𝐷_all, 𝐸𝐷_ks, and 𝐸𝐷_max for 𝑝=2/3.

Table 2 

Comparison of 𝐸𝐷_all, 𝐸𝐷_ks, and 𝐸𝐷_max for 𝑝=1.

We first focus on Table 1 (𝑝=2/3). The 1st and 2nd columns suggest that the “best’’ 𝐸𝐷_all and the “best” 𝐸𝐷_ks that can be achieved are 1.28 by EFRON and 0.137 by K-S, respectively, at the expense of large imbalance under the other criterion. For DSCRT when 𝑚 increases from 2 to 4 and 8, 𝐸𝐷_all increases from 2.17 to 2.94 and 3.76, whereas 𝐸𝐷_ks decreases from 0.178 to 0.161 and 0.159. Therefore, we see that there is a trade-off between the balance of group sizes and the balance of covariate distributions. Similar trend can be observed for WGT-AVE when 𝑐₀ increases from 2 to 4 and 6. In terms of these two criteria, the new procedure (MAX-IMB), with 𝐸𝐷_all 2.36 and 𝐸𝐷_ks 0.159, has better performance than DSCRT with 𝑚=4 and 𝑚=8 and WGT-AVE with 𝑐₀=4 and 𝑐₀=6; it also has lower 𝐸𝐷_all and 𝐸𝐷_ks than RANK-SUM.

In terms of 𝐸𝐷_max (the 3rd column), MAX-IMB has the minimum value 7.38 since it sequentially minimizes this criterion. On the contrary, DSCRT minimizes the imbalance of patient numbers over the selected intervals, but ignores the imbalance over others. As a result, on average, the maximum imbalance under DSCRT is higher than that under MAX-IMB. In a sense, 𝐸𝐷_max serves as a tool which detects any allocation imbalance that is ignored by DSCRT. Since the new procedure examines both “global’’ imbalance, that is, over the whole range, and “local’’ imbalance, that is, over any small interval, it can be regarded as a compromise between achieving balance in overall group sizes and achieving balance in covariate distributions.

Similar conclusion can be drawn for 𝑝=1 (see Table 2), that is, the allocation is deterministic except the case of a tie. From 𝑝=2/3 to 𝑝=1, the decrease in 𝐸𝐷_all is most significant under DSCRT, from (2.17,2.94,3.76) to (0.49,0.93,1.45). This is because when 𝑝=1 only covariate values are random, and so long as the numbers of patients over the selected intervals are even (e.g., 34 over [0,0.5] and 26 over (0.5,1] with 𝑚=2), DSCRT can always achieve perfect overall balance. Moreover, even if the patient numbers are odd (e.g., 35 over [0,0.5] and 25 over (0.5,1]), there are still chances that the allocation differences are +1 and −1 or in the reversed way, again resulting perfect overall balance. Other procedures are more complex and the decrease in 𝐸𝐷_all is less significant. As a result, when 𝑝=1, MAX-IMB is only uniformly better than DSCRT with 𝑚=8, not 𝑚=4.

We also compared the above procedures by other commonly used measures including the mean absolute difference of sample means (𝐸𝐷_mean) and the mean absolute difference of sample standard deviations (𝐸𝐷_std) of the covariate values in the two treatment groups. Furthermore, Lin and Su [23] introduced another criterion, the area between the empirical cumulative distribution functions of the covariate values in the two treatment groups (normalized by the difference of the maximum and the minimum values), denoted as 𝐸𝐷_area, and pointed out that this criterion has better performance than Kolmogorov-Smirnov distance in capturing the difference in two distributions. We thus included this criterion in the simulation. Since the measurements of mean, standard deviation, and area under a distribution function depend on the underlying distribution of the covariate, we did simulation studies under a uniform distribution 𝑍∼Unif[0,1] and under a normal distribution 𝑍∼𝑁(0,1) and show the results in Tables 3 and 4, respectively.

Table 3 

𝑍∼Unif[0,1]. Comparison of 𝐸𝐷_area, 𝐸𝐷_mean, and 𝐸𝐷_std for 𝑝=2/3.

Table 4 

𝑍∼Norm[0,1]. Comparison of 𝐸𝐷_area, 𝐸𝐷_mean, and 𝐸𝐷_std for 𝑝=2/3.

From Table 3 under the uniform distribution, it is seen that K-S has the best performance, since its 𝐸𝐷_mean, 𝐸𝐷_std, and 𝐸𝐷_area (2.29, 1.77, and 5.02) are the lowest among all procedures. This is expected, since K-S solely minimizes the distance between the two distributions. Once the distributions are closest, so are the summary statistics of means and standard deviations as well as the area between the distributions. However, K-S is likely to produce severe imbalance of group sizes, as shown in Tables 1 and 2. EFRON has the worst performance under the three criteria since it completely ignores the covariate distributions.

Among the three choices of 𝑚 under DSCRT, roughly speaking 𝑚=4 performs best: its 𝐸𝐷_mean and 𝐸𝐷_area (3.22 and 5.56) are the lowest and its 𝐸𝐷_std (1.92) slightly higher than that under 𝑚=8. Moreover, under these three criteria, DSCRT with 𝑚=4 is uniformly better than MAX-IMB, RANK-SUM, and WGT-AVE. However, the good performance of DSCRT with 𝑚=4 is based on the correct identification of quartiles of the true covariate distribution, which may not be feasible before the collection of data. In contrast, other methods do not require such information.

Comparing MAX-IMB, RANK-SUM, and WGT-AVE, RANK-SUM has the highest values under the three criteria; MAX-IMB has comparable performance to WGT-AGE with 𝑐₀=4, with the former having slightly higher 𝐸𝐷_mean and the latter slightly higher 𝐸𝐷_std and 𝐸𝐷_area. Similar conclusion can be reached for the normal distribution (see Table 4). The result for 𝑝=1 under the three criteria 𝐸𝐷_mean, 𝐸𝐷_std, and 𝐸𝐷_area resembles that for 𝑝=2/3, the only difference being that the best choice of 𝑚 under DSCRT is 𝑚=8 instead of 𝑚=4. In fact, we also did simulations under different sample sizes (𝑁=30 and 150) and the results are quite consistent.

5. Discussion and Conclusions

In this paper, we propose a new minimization procedure that balances treatment allocation over continuous covariates. For any new patient, it examines the imbalances in the neighborhoods of his or her covariate value and bias the allocation probability towards the treatment that would result in a smaller value of the maximum imbalance. The new method only depends on the ranks of the covariates and is therefore distribution-free. Our simulation studies have shown that it is able to maintain relatively good balance in terms of group sizes and covariate distributions across treatment groups.

In addition, the new procedure does not require the specification of any critical values, which is usually needed for discretization methods in order to define categories. For the latter methods, if quantiles of the covariate distribution 𝐹 are used for the critical values, then lack of knowledge about 𝐹 may lead to wrong guesses of the quantiles. The new procedure saves this step by considering all possible divisions of the range. Nevertheless, only the assignment differences over (𝑛+2) intervals have to be examined to calculate the new imbalance measure, and the corresponding algorithm is computationally efficient.

Borrowing the idea of Pocock and Simon’s design [7], our method can easily be generalized to two or more continuous covariates or a mix of discrete and continuous covariates. Suppose that for a total of 𝐿 covariates 𝑍¹,…,𝑍^𝐿, the first 𝐿₁ is continuous and the rest are discrete. When the (𝑛+1)th patient is enrolled, for any continuous covariate 𝑍^𝑖, 𝑖=1,…,𝐿₁, we define 𝐷^{(𝑘)𝑛+1,𝑖}, 𝑘=𝐴,𝐵 by (3.2), which is the the maximum imbalance measure with respect to the 𝑖th covariate; for any discrete covariate 𝑍^𝑗, 𝑗=𝐿₁+1,…,𝐿, observe the category the new patient belongs to, tentatively assign him to treatment 𝑘, and define 𝐷^{(𝑘)𝑛+1,𝑗} as the absolute difference of patient numbers in the two treatment groups with respect to that specific category. For example, if the 𝑗th covariate is gender and the new patient is a male, then 𝐷^{(𝑘)𝑛+1,𝑗} is calculated among all males. Define 𝐷^{(𝑘)𝑛+1,}_mix=^𝐿1_𝑖=1𝑤_𝑖𝐷^{(𝑘)𝑛+1,𝑖}+^𝐿_𝑗=_𝐿1+1𝑤_𝑗𝐷^{(𝑘)𝑛+1,𝑗},(5.1) where 𝑤_𝑖’s and 𝑤_𝑗’s are the weights placed on the covariates and can be assigned by the importance of the different covariates. Depending on whether 𝐷^{(𝐴)𝑛+1,}_mix is greater than, less than, or equal to 𝐷^{(𝐵)𝑛+1,}_mix, assign the (𝑛+1)th patient to treatment 𝐴 with probability 𝑞, 𝑝 or 1/2. From (5.1), it is seen that 𝐷^{(𝑘)𝑛+1,}_mix is similar to Pocock and Simon’s weighted average of marginal imbalances. The only difference is that for those continuous covariates we redefine the marginal imbalances by the new measure proposed in the current paper, so that the negative effect caused by discretization can be mitigated. Since the marginal imbalances for discrete covariates in 𝐷^{(𝑘)𝑛+1,}_mix remain the same as in Pocock and Simon’s, it is expected that the good balance properties for discrete covariates in their design can be preserved when 𝐷^{(𝑘)𝑛+1,}_mix is used in the minimization. We did simulations for the case of two continuous covariates and the new procedure again showed improvement over other procedures.

In the case that an unequal allocation ratio such as 𝑟_𝐴∶𝑟_𝐵 is desired, one can easily generalize the proposed metric by redefining Δ(𝐼,𝐙_𝑛,𝐓_𝑛) in (3.1) as 𝑁_𝐴𝐼,𝐙_𝑛,𝐓_𝑛𝑟_𝐴−𝑁^𝐵𝐼,𝐙_𝑛,𝐓_𝑛𝑟_𝐵.(5.2) By doing so, it can be ensured that the allocation ratio over any interval is close to 𝑟_𝐴∶𝑟_𝐵. In practice, one can also modify the maximum imbalance measure by adding weights to different intervals. The weight of each interval can be assigned as a function of the number of patients within the interval, so that the procedure remains distribution-free. But the algorithm to implement such a procedure will be more complicated. We will leave these as future research topics.

Appendix

Proof of Proposition 3.1

We will use the basic fact that for a sequence {𝑎₁,…,𝑎_𝑛} and a constant 𝑏, max_{𝑖=1,…,𝑛}||𝑎_𝑖+𝑏||=max|||max_𝑖𝑎_𝑖+𝑏|||,|||min_𝑖𝑎_𝑖+𝑏|||.(A.1) Following the notations and argument in Section 3, 𝐷^(𝐴)𝑛+1=max_{𝑖=𝐿,1,…,𝑛,𝑗=𝑅,1,…,𝑛}||Δ𝑍_𝑖,𝑍_𝑗,𝐙_𝑛+1,𝐓_𝑛+1||∶𝑍_𝑖≤𝑍_𝑛+1≤𝑍_𝑗,𝑇_𝑛+1=𝐴.(A.2) But [𝑍_𝑖,𝑍_𝑗]=[𝑍_𝑖,𝑍_𝑛+1)∪{𝑍_𝑛+1}∪(𝑍_𝑛+1,𝑍_𝑗] and Δ({𝑍_𝑛+1},𝐙_𝑛+1,𝐓_𝑛+1)=1 if 𝑇_𝑛+1=𝐴. Thus, 𝐷^(𝐴)𝑛+1=max_{𝑖=𝐿,1,…,𝑛}max_{𝑗=𝑅,1,…,𝑛}||Δ_𝑛𝑍_𝑖,𝑍_𝑛+1+1+Δ_𝑛𝑍_𝑛+1,𝑍_𝑗||∶𝑍_𝑖≤𝑍_𝑛+1≤𝑍_𝑗.(A.3) Note that by definition of the notation Δ_𝑛𝐼, Δ_𝑛[𝑍_𝑖,𝑍_𝑛+1)=∶Δ([𝑍_𝑖,𝑍_𝑛+1),𝐙_𝑛,𝐓_𝑛), which further equals Δ([𝑍_𝑖,𝑍_𝑛+1),𝐙_𝑛+1,𝐓_𝑛+1) since the interval [𝑍_𝑖,𝑍_𝑛+1) does not contain 𝑍_𝑛+1. Now for any fixed 𝑖, apply (A.1) to the constant Δ_𝑛[𝑍_𝑖,𝑍_𝑛+1)+1 and the sequence Δ_𝑛(𝑍_𝑛+1,𝑍_𝑗], 𝑗=𝑅,1,…,𝑛, and we have 𝐷^(𝐴)𝑛+1==mmaxax_{𝑖=𝐿,1,…,𝑛}mmaaxx_{𝑖=𝐿,1,…,𝑛}||𝐶_𝑅2+Δ_𝑛𝑍_𝑖,𝑍_𝑛+1+1||,||𝐶_𝑅1+Δ_𝑛𝑍_𝑖,𝑍_𝑛+1+1||∶𝑍_𝑖≤𝑍_𝑛+1max_{𝑖=𝐿,1,…,𝑛}||𝐶_𝑅1+1+Δ_𝑛𝑍_𝑖,𝑍_𝑛+1||∶𝑍_𝑖≤𝑍_𝑛+1,=∶max{𝐼,𝐼𝐼}||.𝐶_𝑅2+1+Δ_𝑛𝑍_𝑖,𝑍_𝑛+1||∶𝑍_𝑖≤𝑍_𝑛+1(A.4) Applying (A.1) again to 𝐼 with constant 𝐶_𝑅1+1 and sequence Δ_𝑛[𝑍_𝑖,𝑍_𝑛+1), we have 𝐼=max||𝐶_𝐿1+𝐶_𝑅1+1||,||𝐶_𝐿2+𝐶_𝑅1+1||.(A.5) Similarly, 𝐼𝐼=max||𝐶_𝐿1+𝐶_𝑅2+1||,||𝐶_𝐿2+𝐶_𝑅2+1||.(A.6) Therefore,𝐷^(𝐴)𝑛+1=max||𝐶_𝐿1+𝐶_𝑅1+1||,||𝐶_𝐿2+𝐶_𝑅1+1||,||𝐶_𝐿1+𝐶_𝑅2+1||,||𝐶_𝐿2+𝐶_𝑅2+1||.(A.7) The derivation of 𝐷^(𝐵)𝑛+1 is similar.

Acknowledgments

This work was supported by Grants DMS-0907297 and DMS-0906661 from the National Science Foundation (USA).