On estimation of overall treatment effects in multiregional clinical trials under a discrete random effects model

2.1. DREM with binary endpoints

To address the heterogeneous treatment effect across regions in an MRCT, the DREM was proposed by Lan and Pinheiro. ³ In this section, we present the application of the DREM to binary endpoints ⁵ by taking the OR as an example.

We consider the patient population to be partitioned into K disjointed clinical regions. The regional weight of region k is represented by w k , k = 1 , … , K , with ∑ k = 1 K w k = 1 , where w k is a parameter rather than a random variable. ³ In a randomized control trial, patients are randomly assigned to either the treatment group or the control group. Suppose that the patient-allocation ratio of the treatment group and the control group is T : C and let π = T / ( T + C ) , which represents the overall treatment allocation ratio. The treatment allocation ratio for region k is denoted as π k , k = 1 , … , K , where π k is a parameter rather than a random variable. ³ Note that the overall treatment allocation ratio π can be expressed as the weighted sum of π k , i.e., π = ∑ k = 1 K w k π k .

When the primary endpoint is binary, to compare the treatment effect in the treatment group and the control group, the overall treatment effect is estimated by the event rates of the two groups. The Bernoulli random variable X i j k is used to indicate whether the event occurs or not for the jth patient in the kth region for group i = T , C . We assume that

Under the DREM, the treatment effect for a randomly selected patient in the population is a random variable with the distribution defined by the possible values { θ 1 , θ 2 , … , θ K } and their respective probabilities { w 1 , w 2 , … , w K } . The mean is θ and the between-region variance ( τ 2 ) is expressed as

OPEN IN VIEWER

Table 1.

Between-region variance τ 2 for event rates p T k and p C k , given ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) and ( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 ) .

p T 1	p T 2	p T 3	p C 1	p C 2	p C 3	τ 2
0.50	0.60	0.60	0.65	0.70	0.70	0.0071
0.55	0.60	0.65	0.70	0.75	0.80	0.0024
0.40	0.40	0.45	0.60	0.60	0.60	0.0088

When a patient from the population is randomly selected, the variance of his/her treatment effect should not only consider the between-region variance τ 2 , but also take into account the within-region variance. For a patient from region k in group i, the within-region variance is the variance of Bernoulli ( p i k ) , which is denoted as σ i k 2 = p i k ( 1 − p i k ) . As a result, we define the variance of the random variable X i j k as

So far, we have defined the overall treatment effect θ (Eq. (3)) and have obtained the naïve estimator θ ^ (Eq. (5)). We will now discuss the potential bias of θ ^ and how to quantify it. As a common protocol is used in an MRCT, we generally assume that the overall treatment allocation ratio π is applied across all K regions, meaning that π = π 1 = … = π K . However, in practice, at the final analysis stage, the treatment allocation ratios π ^ k ′ s and the regional weights w ^ k ′ s may differ from those predetermined at the design stage. This reduces the accuracy of estimating the overall treatment effect. In this article, we take into account the gap between the design stage and the actual implementation results with the following assumptions:

(1)

the observed regional weights w ^ k ′ s are possibly different from the predetermined regional weights w k ′ s in the protocol, and

(2)

the observed treatment allocation ratios π ^ k ′ s are possibly different across regions.

Although efforts are made to enroll patients from each region to achieve predetermined regional weights w k ′ s and a predetermined treatment allocation ratio π , they are often different in practice. We recognize this “real-world” challenge and attempt to fully understand its impact on the overall treatment effect. We anticipate that the observed naïve estimator θ ^ would not be equal to θ ^ based on the predetermined patient assignments.

One approach to empirically quantify the bias of θ ^ is to use conditional inference with ancillary statistics for θ .^6,9 There is a vast amount of literature on how to adjust the estimator based on baseline covariate information under unconditional inference.^10–13 In this study, we focus on the conditional distribution of θ ^ − θ given appropriate ancillary statistics, as it is more informative than unconditional inference to study the stochastic behavior of θ ^ .^7,8 Here, both the observed regional weights { w ^ k , k = 1 , … , K } and the observed treatment allocation ratios { π ^ k , k = 1 , … , K } are ancillary statistics.

We consider the naïve estimator θ ^ to be generated through random sampling, but { ( w ^ k , π ^ k ) , k = 1 , … , K } of each sample would be the same as the observed patient assignment. We also assume random sampling of patients from each region and then their random assignment to one of two groups. Within the framework of the DREM, we extend the method proposed by Tian et al. ⁶ By the Central Limit Theorem, it can be readily shown that the conditional distribution of ( θ ^ − θ ) | { ( w ^ k , w k , π ^ k , π ¯ ) , k = 1 , … , K } is approximately normal, with the mean of

b θ = 1 ( 1 − p T ) p T ∑ k = 1 K p T k ( w ^ k π ^ k π ¯ − w k ) − 1 ( 1 − p C ) p C ∑ k = 1 K p C k ( w ^ k ( 1 − π ^ k ) 1 − π ¯ − w k ) ,

and the variance of

σ π w 2 = ( 1 ( 1 − p T ) p T ) 2 ∑ k = 1 K w k 2 V a r ( p ^ T k ) + ( 1 ( 1 − p C ) p C ) 2 ∑ k = 1 K w k 2 V a r ( p ^ C k ) .

(6)

It is worth noting that the variance of p ^ i k accounts for both between-region and within-region variations. According to Eq. (4), the variance of p ^ i k can be derived as

Var ( p ^ i k ) = σ i k 2 + τ 2 N i k , i = T , C .

(7)

Next, let p i * = ∑ k = 1 K w k p ^ i k be another estimator of p i (Eqs. (1) and (2)). Note that if w ^ k = w k and π ^ k = π ¯ , k = 1 , … , K , then the estimator p ^ i is equal to the estimator p i * . Following Tian et al., ⁶ define p T * = ∑ k = 1 K w k p ^ T k and p C * = ∑ k = 1 K w k p ^ C k . Using the aforementioned notation, we construct an estimator of b θ as

b ^ θ = 1 ( 1 − p T * ) p T * ∑ k = 1 K p ^ T k ( w ^ k π ^ k π ¯ − w k ) − 1 ( 1 − p C * ) p C * ∑ k = 1 K p ^ C k ( w ^ k ( 1 − π ^ k ) 1 − π ¯ − w k ) ,

and an estimator of σ π w 2 as

σ ^ π w 2 = ( 1 ( 1 − p T * ) p T * ) 2 ∑ k = 1 K w k 2 σ ^ T k 2 + τ ^ 2 N T k + ( 1 ( 1 − p C * ) p C * ) 2 ∑ k = 1 K w k 2 σ ^ C k 2 + τ ^ 2 N C k ,

where σ ^ i k 2 = p ^ i k ( 1 − p ^ i k ) , i = T , C . By adjusting θ ^ directly, we can obtain

θ ^ π w = θ ^ − b ^ θ .

(8)

The bias-adjusted estimator θ ^ π w is the sum of θ ^ and a function of ( w ^ k , w k , π ^ k , π ¯ ) , k = 1 , … , K . Note that θ ^ π w is a consistent estimator of θ . Moreover, the conditional distribution ( θ ^ π w − θ ) | { ( w ^ k , w k , π ^ k , π ¯ ) , k = 1 , … , K } under the null hypothesis is approximately normal with mean 0 and variance σ π w 2 .

3. Power analysis and sample size determination under the DREM

3.2. Sample size calculation for the naïve estimator

In this section, we introduce the distribution of the naïve estimator θ ^ and derive the required sample size. The estimator θ ^ follows an approximately normal distribution with the mean of θ and the variance of Var ( θ ^ ) . To calculate the variance Var ( θ ^ ) , let g ( p ) = log ( p / ( 1 − p ) ) , which is the log odds function of the event rate p. Then, through Taylor's expansion, g ( p ^ ) ≈ g ( p ) + ( p ^ − p ) ( 1 / ( p ( 1 − p ) ) ) . Taking the variance on both sides, we obtain

Var ( g ( p ^ ) ) ≈ ( 1 p ( 1 − p ) ) 2 Var ( p ^ ) .

(11)

Under the DREM, substituting Eq. (7) into Eq. (11) and considering the log(OR) of the treatment and control groups, we obtain

V a r ( θ ^ ) = 1 N [ 1 π ( σ T 2 + τ 2 p T 2 ( 1 − p T ) 2 ) + 1 ( 1 − π ) ( σ C 2 + τ 2 p C 2 ( 1 − p C ) 2 ) ] ,

where σ i 2 = p i ( 1 − p i ) , i = T , C . Therefore, the test statistic for the naïve estimator θ ^ is given by Z n e = ( θ ^ − θ ) / var ( θ ^ ) , and the required sample size for θ ^ is

N n e = ( z 1 − α + z 1 − β δ ) 2 [ 1 π ( σ T 2 + τ 2 p T 2 ( 1 − p T ) 2 ) + 1 ( 1 − π ) ( σ C 2 + τ 2 p C 2 ( 1 − p C ) 2 ) ] .

(12)

We have thus derived the sample sizes N π w (Eq. (10)) and N n e (Eq. (12)) for θ ^ π w and θ ^ , respectively.

With α = 2.5 % , β = 20 % , ( p T 1 , p T 2 , p T 3 ) = ( 0.50 , 0.60 , 0.60 ) , ( p C 1 , p C 2 , p C 3 ) = ( 0.65 , 0.70 , 0.70 ) , and ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) , along with τ 2 = 0.0071 , Table 2 shows the sample sizes N n e and N π w required for given observed patient assignments ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) . Although a treatment allocation of T : C = 1 : 1 is common in an MRCT, there may be a higher proportion in the treatment arm in some clinical trials.^14,15 Therefore, we considered two cases: π k ′ s = 0.5 and π k ′ s = 0.7 , as shown in the top and bottom halves of Table 2, respectively. In both cases, the first row shows that the observed patient assignments were equal to the designed ones, i.e., w ^ k = w k and π ^ k = π k for all regions k, k = 1 , 2 , 3 . Since there are infinitely many combinations of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } , the subsequent rows show some randomly generated { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } as described below. The observed regional weights w ^ k ′ s were generated from a uniform distribution in the range [ 0.1 , 0.7 ] with ∑ k = 1 3 w ^ k = 1 . For π k ′ s = 0.5 and π k ′ s = 0.7 , the observed treatment allocation ratios π ^ k ′ s were randomly generated from the ranges [ 0.3 , 0.7 ] and [ 0.55 , 0.85 ] , respectively.

OPEN IN VIEWER

Table 2.

The required sample size for θ ^ and θ ^ π w with ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) using α = 2.5 % and 1 − β = 80 % .

w ^ 1	w ^ 2	w ^ 3	π ^ 1	π ^ 2	π ^ 3	N n e	N π w
( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 )
0.350	0.350	0.300	0.500	0.500	0.500	554	551
0.462	0.318	0.220	0.692	0.328	0.354	554	651
0.235	0.216	0.549	0.674	0.308	0.430	554	789
0.352	0.387	0.261	0.376	0.454	0.514	554	566
( π 1 , π 2 , π 3 ) = ( 0.7 , 0.7 , 0.7 )
0.350	0.350	0.300	0.700	0.700	0.700	677	674
0.526	0.260	0.214	0.606	0.555	0.599	677	640
0.216	0.591	0.193	0.672	0.694	0.672	677	805
0.339	0.310	0.351	0.661	0.683	0.778	677	696

In the case of π k ′ s = 0.5 , while w ^ k = w k and π ^ k = π k for all k, the required sample sizes are N n e = 554 and N π w = 551 . Moreover, regardless of the combination of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } , the sample size N n e remains the same. On the other hand, depending on the given { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } , the sample size N π w will be different accordingly. For instance, given ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) = ( 0.462 , 0.318 , 0.220 , 0.692 , 0.328 , 0.354 ) , we obtained N π w = 651 for the desired power level of 80 % . In contrast to the proposed estimator θ ^ π w , θ ^ does not take into account the difference in patient assignments between the design stage and the final analysis stage. Therefore, when calculating the sample size N n e , it will not change with { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } . In the case of π k ′ s = 0.7 , a similar result was observed. As shown in Table 2, N π w may be greater or less than N n e , depending on the combination of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } . When considering the predetermined patient assignments at the design stage and the possible observed patient assignments, it becomes evident that conducting a trial with the sample size N n e may not guarantee adequate statistical power. This deficiency in power can be attributed to the bias of the naïve estimator θ ^ . This study addresses the issue of loss of statistical power by hypothetically adjusting the required sample size at the design stage.

To evaluate the performance of the proposed approach under the DREM, we compared the bias-adjusted estimator θ ^ π w with the naïve estimator θ ^ . For simplicity, we considered three regions. Let the event rates of the treatment group and the control group be { p T k } = ( 0.50 , 0.60 , 0.60 ) and { p C k } = ( 0.65 , 0.70 , 0.70 ) , respectively. The regional treatment effects are ( θ 1 , θ 2 , θ 3 ) = ( − 0.619 , − 0.442 , − 0.442 ) . At the design stage, we set the patient assignments as ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) and ( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 ) . The between-region variance τ 2 was 0.0071. To detect the expected treatment effect θ < 0 in favor of the treatment group at a significance level of α = 2.5 % and power of 1 − β = 80 % , the required sample size was N n e = 554 . It is important to note that both θ ^ π w and θ ^ were evaluated using a sample size of 554 for comparison. The true parameter θ is − 0.504 .

In order to investigate various scenarios that may occur in practice, we generated different combinations of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } . The observed regional weights w ^ k ′ s were randomly generated from a uniform distribution within the range of [ 0.1 , 0.7 ] with a constraint on ∑ w ^ k = 1 . This wide range allows us to cover a diverse set of possible values for the observed regional weights. The observed treatment allocation ratios π ^ k ′ s were randomly generated from a uniform distribution with the mean of 0.5 and the range of [ 0.3 , 0.7 ] . By bounding the observed treatment allocation ratio between this range, we avoid values that are too large or too small, as this indicates a very poor randomization scheme, which is rarely the case in practice. For each combination of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } , we performed 10,000 simulations to assess the performance of the two estimators. The simulation results shown in Table 3 include empirical bias, standard deviation, root mean squared error (RMSE), empirical test size, and empirical power. To interpret the size of the biases, we considered the standardization of b ^ θ and defined it as

The value of b ^ θ * is also provided in Table 3. Note that the results presented in Table 3 are sorted by the absolute values of b ^ θ * .

We started with scenarios where the value of b ^ θ * is negative, taking the first row of Table 3 as an example. It can be observed that θ ^ has a bias of − 0.134 , while θ ^ π w is almost unbiased. This indicates that our approach successfully corrected the bias. Considering a significance level of α = 2.5 % , the empirical test size of θ ^ is 6.4 % , and the empirical test size of θ ^ π w is 2.5 % . The inflated test size of θ ^ is primarily due to its large negative bias. In terms of standard deviation, θ ^ π w has a slightly higher value of 0.193 compared to 0.177 for θ ^ . The RMSE of θ ^ π w is 0.193, which is less than the value of 0.222 for θ ^ . The estimator θ ^ has the empirical power of 95.2 % , while θ ^ π w has the empirical power of 74.2 % . To explain the unusually high empirical power of θ ^ , we need to look at its corresponding test size and standard deviation. When the test size is as high as 6.4 % , such high power is not meaningful.

Next, let's turn our attention to the scenarios where the value of b ^ θ * is positive, taking the fifth row of Table 3 as an example. The naïve estimator θ ^ is overestimated with a bias of 0.095, resulting in the empirical test size of 0.9 % , which is well below 2.5 % . The excessively small test size of θ ^ leads to the empirical power of 62.6 % , which is well below the desired power of 80 % . In contrast, θ ^ π w is nearly unbiased with the empirical test size of 2.4 % and the empirical power of 76.1 % .

Across all the combinations of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } presented in Table 3, regardless of the value of b ^ θ * , it is evident that our proposed estimator θ ^ π w is almost unbiased. When the observed patient assignments { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } remarkably deviate from the designed assignments { ( w k , π k ) , k = 1 , 2 , 3 } , indicated by b ^ θ * being far from 0, the proposed estimator clearly outperforms the naïve estimator. The empirical test size of θ ^ π w is always close to 2.5 % compared to that of θ ^ . The empirical power of θ ^ π w is more meaningful when the test size is stable and close to 2.5 % compared to that of θ ^ . However, the standard deviation of θ ^ π w is slightly greater than that of θ ^ . This is because our approach sacrifices some precision in the process of bias adjustment. Therefore, when there is almost no bias ( b ^ θ * is close to 0), the RMSE of θ ^ π w will be slightly larger than that of θ ^ .

These results can also be clearly seen in Figure 1, where the x-axis represents b ^ θ * and each dot represents a combination of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } . A total of 500 combinations were compared between the two estimators, and each combination represents the result of 10,000 simulations. Figure 1 depicts the relationships between b ^ θ * and the empirical bias, standard deviation, test size, and power. When the size of bias b ^ θ * is far from 0, the empirical test size of θ ^ is far from 2.5 % and the empirical power of θ ^ is far from 80 % . On the other hand, since the empirical test size of the biased-adjusted estimator θ ^ π w is stable around 2.5 % , the empirical power of θ ^ π w is statistically informative. Figure 2 displays a histogram comparing the RMSEs of the two estimators for the 500 combinations. Since the RMSE is affected by the bias and standard deviation and our approach sacrifices some precision in the process of adjusting the bias, our approach will show a greater advantage when | b ^ θ * | is large.

OPEN IN VIEWER

Figure 1.

Comparisons of the empirical bias, standard deviation, test size, and power of θ ^ and θ ^ π w , given ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) , ( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 ) , and 500 combinations of ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) .

OPEN IN VIEWER

Figure 2.

Comparisons of the RMSEs of θ ^ and θ ^ π w when ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) and ( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 ) . The absolute values of b ^ θ * were obtained from 500 combinations of ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) . The darker area means that the RMSE value of θ ^ is greater than or equal to θ ^ π w ; otherwise, it is represented by the lighter area.

In the second set of predetermined patient allocation settings, we have ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) and ( π 1 , π 2 , π 3 ) = ( 0.7 , 0.7 , 0.7 ) . We conducted 10,000 simulations for each combination { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } , considering a total of 500 combinations for analysis. For this set, the observed regional weights w ^ k ′ s were randomly generated from a uniform distribution within the range of [ 0.1 , 0.7 ] with a constraint on ∑ w ^ k = 1 . The observed treatment allocation ratios π ^ k ′ s were randomly generated from a uniform distribution with the mean of 0.7, and the range of [ 0.55 , 0.85 ] . The required sample size N n e was determined to be 667. The results of the simulation study are presented in Figures 3 and 4. Figure 3 illustrates the relationships between b ^ θ * and the bias, standard deviation, test size, and power for the two estimators θ ^ π w and θ ^ . Figure 4 presents a histogram comparing the RMSEs of the two estimators for the 500 combinations of { ( w ^ k , π ^ k ) , k = 1 , 2 , 3 } . Similar to the previous analysis, the findings shown in Figures 3 and 4 align with the conclusions drawn from Figures 1 and 2.

OPEN IN VIEWER

Figure 3.

Comparisons of the empirical bias, standard deviation, test size, and power of θ ^ and θ ^ π w , given ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) , ( π 1 , π 2 , π 3 ) = ( 0.7 , 0.7 , 0.7 ) , and 500 combinations of ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) .

OPEN IN VIEWER

Figure 4.

Comparisons of the RMSEs of θ ^ and θ ^ π w when ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) and ( π 1 , π 2 , π 3 ) = ( 0.7 , 0.7 , 0.7 ) . The absolute values of b ^ θ * were obtained from 500 combinations of ( w ^ 1 , w ^ 2 , w ^ 3 , π ^ 1 , π ^ 2 , π ^ 3 ) . The darker area means that the RMSE value of θ ^ is greater than or equal to θ ^ π w ; otherwise, it is represented by the lighter area.

When the observed regional weights w ^ k ′ s and treatment allocation ratios π ^ k ′ s closely align with the design, the bias of b ^ θ * is close to 0, and our proposed estimator θ ^ π w demonstrates its advantages by reducing bias, maintaining a stable test size close to 2.5 % , and providing informative empirical power. To demonstrate this, we adopted the same designed set as the first set in this section, and assumed the observed regional weights w ^ k ′ s were randomly generated from a uniform distribution within the range of [ 0.25 , 0.40 ] , the observed treatment allocation ratios π ^ k ′ s were randomly generated from a uniform distribution within the range of [ 0.45 , 0.55 ] . Such assumptions about the range of observed values w ^ k ′ s and π ^ k ′ s are not exceptional but rather represent a common scenario. The results are summarized in Figure 5, which illustrates the relationships between b ^ θ * and the bias, standard deviation, test size, and power for the two estimators θ ^ π w and θ ^ . The stability is demonstrated by the consistent performance of our method despite the variability of b ^ θ * . Therefore, our method demonstrates better performance.

OPEN IN VIEWER

Figure 5.

Comparisons of the empirical bias, standard deviation, test size, and power of θ ^ and θ ^ π w , given ( w 1 , w 2 , w 3 ) = ( 0.35 , 0.35 , 0.30 ) , ( π 1 , π 2 , π 3 ) = ( 0.5 , 0.5 , 0.5 ) , and 500 combinations of ( w ^ k , π ^ k ) , k = 1 , 2 , 3 , where 0.25 ≤ w ^ k ≤ 0.40 and 0.45 ≤ π ^ k ≤ 0.55 .