The following notations are consistently used throughout this article: 1.

y c k = ( y c k 1 , … , y c k n c k ) ′ , y t k = ( y t k 1 , … , y t k n t k ) ′ , X c k = ( x c k 1 , … , x c k n c k ) ′ , X t k = ( x t k 1 , … , x t k n t k ) ′ , y k = y c k y t k , X k = X c k X t k .

They denote data observations from the k-th basket.

We would like to compare the performance of the Bayesian methods with a separate analysis under the Bayesian framework. Under models adjusting for covariates, the linear model is (3.1) y c k = X c k β k + ϵ c k , k = 1 , … , K , y t k = X t k β k + τ k 1 n t k + ϵ t k , k = 1 , … , K , where the error terms ϵ c k and ϵ t k follow M V N n c k ( 0 , I n c k σ 2 ) and M V N n t k ( 0 , I n t k σ 2 ), respectively. In the k-th basket, the likelihood functions are given by (3.2) L ( β k , σ 2 | y c k , X c k ) = 1 ( 2 π σ 2 ) n c k exp ( − 1 2 σ 2 ‖ y c k − X c k β k ‖ 2 ) , L ( β k , τ k , σ 2 | y t k , X t k ) = 1 ( 2 π σ 2 ) n t k exp ( − 1 2 σ 2 | | y t k − τ k 1 n t k − X t k β k | | 2 ) . The complete likelihood function is (3.3) L ( β , τ , σ 2 | D ) = ∏ k [ L ( β k , σ 2 | y c k , X c k ) L ( β k , τ k , σ 2 | y t k , X t k ) ] . Assuming conjugate priors for parameters β, τ, σ 2 , (3.4) β k | σ 2 ∼ iid M V N r + 1 ( 0 , σ 2 Λ k − 1 ) , k = 1 , … , K , τ | σ 2 ∼ M V N K ( 0 , σ 2 Λ τ − 1 ) , σ 2 ∼ I G ( a 0 , b 0 ) . If no prior knowledge is available, non-informative priors are assumed, for example, by setting Λ k = diag r + 1 ( 1 10000 ), k = 1 , … , K, and Λ τ = diag K ( 1 10000 ).

The full posterior distribution is (3.5) p ( β , τ , σ 2 | D ) ∝ L ( β , τ , σ 2 | D ) ∏ k π ( β k | σ 2 ) π ( τ | σ 2 ) π ( σ 2 ) . The posterior for τ can be derived by integrating out the nuisance parameters: (3.6) p ( τ | D ) ∝ ∬ L ( β , τ , σ 2 | D ) ∏ k π ( β k | σ 2 ) π ( τ | σ 2 ) π ( σ 2 ) d β d σ 2 . Hence, the posterior distribution of τ is a K-dimensional multivariate t distribution with degrees of freedom ν = 2 a 0 + ∑ k n c k + ∑ k n t k . Several quantities are defined below in order to express the location parameters and the scale matrix of the posterior multivariate t distribution: (3.7) Λ D k = Λ k + X k ′ X k , μ k ∗ = X k ′ y k , c k = 1 n t k ′ 1 n t k − 1 n t k ′ X t k Λ D k − 1 X t k ′ 1 n t k , d k = y t k ′ 1 n t k − μ k ∗ ′ Λ D k − 1 X t k ′ 1 n t k ; and denote diagonal matrix C and vector d C = diag ( c 1 , … , c K ) , d = ( d 1 , … , d K ) ′ . The location parameter for the posterior of τ is (3.8) μ τ = ( C + Λ τ ) − 1 d . The scale matrix is (3.9) Σ = 2 B ν ( C + Λ τ ) , where B = b 0 + 1 2 ∑ k ( y k ′ y k − μ k ∗ ′ Λ D k − 1 μ k ∗ ) − 1 2 d ′ ( C + Λ τ ) − 1 d. Given the above forms, τ k ’s independently follow shifted and scaled t-distributions, and their posterior probabilities can be computed directly.

Berry et al. [4] and Thall et al. [36] introduced a Bayesian adaptive design with frequentist interim analyses and hierarchical modeling across the patient subgroups. In addition to model (3.1), a additional structure is built to model τ k ’s, i.e., τ k ∼ N ( μ τ , σ τ 2 ). Non-informative normal and inverse gamma priors are given to μ τ and σ τ 2 , respectively. Since all τ k ’s are sharing the same mean parameter μ τ in their priors, one can expect that each estimate of them is pulled towards μ τ . The shrinkage parameter σ τ 2 indicates the amount of information borrowing intensity between the baskets: larger σ τ 2 indicates less borrowing, and smaller σ τ 2 indicates stronger borrowing. To control the shrinkage parameter, we can define σ τ 2 = ϕ σ 2 and set a large value of ϕ to form a non-informative prior for τ k ’s.

The model is described in (3.1). The likelihood function is described in (3.2) and (3.3). The priors are assumed to be (3.10) β k ∼ ind M V N r + 1 ( 0 , σ 2 Λ k − 1 ) , k = 1 , … , K , τ k ∼ iid N ( μ τ , σ τ 2 ) , k = 1 , … , K , μ τ ∼ N ( 0 , ϕ σ τ 2 ) , σ τ 2 ∼ I G ( a 0 ′ , b 0 ′ ) , σ 2 ∼ I G ( a 0 , b 0 ) . To specify non-informative priors, set Λ k = diag r + 1 ( 1 10000 ) for k = 1 , … , K. The Gibbs collapsed sampling procedure is used to obtain the posterior samplers.

The challenge of achieving advantageous results for BHM to basket trials is that the number of baskets, K, is often small so the parameter σ τ 2 cannot be estimated robustly and the results are very sensitive for the prior setting. The detailed discussion can be found in [4]. To overcome this issue, Chu and Yuan [6] proposed a Calibrated Bayesian Hierarchical model (CBHM) for binary endpoints. Instead of giving a prior to σ τ 2 , they proposed to estimate σ τ 2 as a monotonic increasing function of a Chi-square test statistic T, where T is a quantity that measures the similarity of the treatment effects across baskets (i.e. Chi-square test statistics): (3.11) σ τ 2 = exp ( a + b log T ) , where a and b > 0 are pre-calculated constants obtained by simulation. The information is borrowed more if the treatment effects are similar across baskets, and less otherwise. To determine the values of a and b, one needs to pre-specify the two cases of σ τ 2 ’s representing the data being heterogeneous and homogeneous (e.g., σ τ 2 = 80 versus σ τ 2 = 1). In the simulation, one also needs to first decide on some scenarios where information should be borrowed (e.g., the treatment effects are the same for all baskets) across baskets and some other scenarios where information should not be borrowed (e.g., only one basket is truly efficacious), and then record the median of their Chi-square statistics. The constants a and b could be calculated by solving (3.11) under these two scenarios.

We extend this idea to the linear regression models, by replacing the Chi-square test statistic with the p-value of a F test statistic. The F test statistic is to test whether any treatment effect difference are heterogeneity from other baskets when the covariates are included in the model. Since the F test statistics is associated with different degrees of freedom, which depends on the sample size of the basket, number of covariates adjusted and the number of baskets being assessed, we use their p-values instead of the test statistics. See Section S.4 of the Supplementary Material for the formulae of the testing procedure. The main challenge may be to set a suitable function that relates the value of σ τ 2 to the p-value ( p v a l). Note that the p-value is always between 0 and 1, and σ τ 2 should monotonically decrease with the p-value. We simply apply the following relation: (3.12) σ τ 2 = a + b × ( 1 − p v a l ) , where a , b ≥ 0 are the tuning parameters. To find a robust function between σ τ 2 and p-values, some other functions are also tested, including the original exponential function linkage as well as the logit function of p v a l. The property of the exponential function may arise large variances in our setting and thus cause Type I error inflation. The logit function performs best on Type I error rate control, whereas the power improvement may not be as satisfactory. Compared with these functions, the above linear function performs quite robustly, and has good Type I error rate control and power improvement. A simplified procedure is used by adjusting the a and b directly instead of running simulations. Using these procedures, comparable results are produced with those by Chu and Yuan [6] in terms of power and Type I error rate comparison with other methods.

The model is described in (3.1). The likelihood function is described in (3.2) and (3.3). The priors for the parameters are given as (3.13) β k ∼ ind M V N r + 1 ( 0 , σ 2 Λ k − 1 ) , k = 1 , … , K , τ k ∼ N ( τ , σ ˆ τ 2 ) , k = 1 , … , K , τ ∼ N ( 0 , σ 0 2 ) , σ 2 ∼ I G ( a , b ) . Non-informative priors are specified by setting Λ k = diag r + 1 ( 1 10000 ) for k = 1 , … , K. The Gibbs collapsed sampling procedure to obtain the posterior samples is given in Algorithm 1.

Inferences from simple pooling or separate analyses are recognized to be inferior to hierarchical modeling methods, which allow adaptive borrowing of information across subgroups. However, the full Bayesian model bears the risk of too much shrinkage and excessive borrowing. In recent years, more and more attention is given to approaches based on local borrowing through mixture models. In such models, baskets can be grouped to clusters so that the information is borrowed within each cluster. Choosing or modeling the number of clusters M is critical when the mixture model is applied. It is even more important when applied to basket trials, as there are usually a limited number of baskets, and hence the overall model performance heavily relies on the model of M. We attempt to apply the Mixture of Finite Mixtures (MFM) approach [18], for its flexible and targeted control of modeling M.

The full MFM model can be specified as usual: (3.18) M ∼ p ( m ) : truncated Poisson p.m.f. on positive integer m with parameter λ , π 1 , … , π M ∼ D i r i c h l e t M ( γ , … , γ ) , P ( z k = m | M ) = π m , m = 1 , … , M , k = 1 , … , K , where z k , the latent variable, denotes cluster that the k-th basket belongs to, and π m is the corresponding probability P ( z k = m | M ). The truncated Poisson prior is chosen because of its convenience in sampling from the posterior distribution [18]. Several parameters remain unspecified in the above model. The value of λ indicates the prior choice of the number of clusters. Setting λ equal to K seems to help alleviate the excessive borrowing. γ affects the probability of assigning the basket to a distinct new cluster, hence influencing the borrowing strength across baskets.

The model is described in (3.1). The likelihood function is described in (3.2) and (3.3) by replacing τ k by τ z k . The priors for the parameters are given as follows: (3.19) β k ∼ iid M V N r + 1 ( 0 , σ 2 Λ k − 1 ) , k = 1 , … , K , τ m ∼ iid N ( 0 , ϕ σ 2 ) , m = 1 , … , M , σ 2 ∼ I G ( a 0 , b 0 ) , with Λ k = diag r + 1 ( 1 ξ , 1 10000 , … , 1 10000 ). For ϕ and ξ, they together control the amount of information to be borrowed across the baskets and therefore need to be carefully specified. A programmable algorithm to fit the MFM with covariates adjustment is given in Algorithm 2. In the Algorithm, as a priori the basket k is placed in (3.25) an existing cluster g ∈ G − k with probability ∝ | g | + γ , a new cluster with probability ∝ V K ( t + 1 ) V K ( t ) γ , where t = | G − k | is the number of clusters obtained by removing the k-th basket, and V K ( t ) needs to be precomputed as (3.26) V K ( t ) = ∑ m = 1 ∞ m ( t ) ( γ m ) ( K ) p M ( m ) . Here x ( t ) = x ( x + 1 ) … ( x + t − 1 ), and x ( t ) = x ( x − 1 ) … ( x − t + 1 ). By convention, x ( 0 ) = 1 and x ( 0 ) = 1. Meanwhile, p M ( m ) is the mass density function of the truncated Poisson distribution on {1, 2, …} with parameter λ.

The derivation of m ( D | β k , σ 2 ) in (3.24) is given in Section S.3 of the Supplementary Material.

Consider the case that the outcomes are influenced by two covariates, one continuous ( X 1 ) and one binary ( X 2 ). Following previous notation, the outcomes are modeled as (4.1) Y c k i = β k 0 + β k 1 X c k i 1 + β k 2 X c k i 2 + ϵ c k i , with i = 1 , … , n c k , Y t k i ′ = β k 0 + τ k + β k 1 X t k i ′ 1 + β k 2 X t k i ′ 2 + ϵ t k i ′ , with i ′ = 1 , … , n t k , where k = 1 , … , 4. The number of subjects in control and treatment arm are assumed to be the same ( n c k = n t k ≐ n k ). For baskets 1 to 4, they are assumed to be n 1 = 30, n 2 = 30, n 3 = 20, and n 4 = 20, and the same sample sizes are assumed for each stage (if the basket is available for analysis in that stage). Let β k 0 denotes the intercept in the control arm for the k-th basket. τ k denotes the difference in the intercept between treatment and control for the k-th basket. Each of the error terms ϵ c k i ( ϵ t k i ′ ) ∼ N ( 0 , 1 ). The coefficients are allowed to be different across baskets, and set to be ( β 11 , β 21 , β 31 , β 41 ) = ( 0.2 , 0.4 , 0.6 , 0.8 ). The continuous covariates X c k i 1 in the control arm are simulated from l o g N ( 0 , 0 . 9 2 ), and the covariates X t k i ′ 1 in the treatment arm are simulated from N ( 1.5 , 2.8 ). The binary covariates X c k i 2 in the control arm are simulated from B i n ( 0.4 ), and the covariates X t k i ′ 2 in the treatment arm are simulated from B i n ( 0.6 ). The values of the binary coefficients are assumed to differ by basket to be ( β 12 , β 22 , β 32 , β 42 ) = ( 0.4 , 0.3 , 0.2 , 0.1 ).

To evaluate the performance of these various methods, the data is generated assuming the control effects to be ( 0 , 0 , 0 , 0 ). Consider 6 scenarios of the true treatment effects, τ = ( τ 1 , … , τ 4 ), which are summarized in Table 1. Note that the first scenario corresponds to the global null (GN).

Under each scenario, 1000 trials are simulated. The following characteristics of all methods are evaluated under each scenario: 1.

FWER: The percentage of trials in which any ineffective basket is wrongly selected as efficacious.

For all methods, if not specified, the values of ξ and ϕ are 10,000, and the hyper-parameters a 0 = 0.5 and b 0 = 0.05, forming non-informative priors for the corresponding parameters. BHM, CBHM and MFM require sampling procedures. The MCMC sample size is 1300, with 300 burn-in. Trace plots are monitored to ensure convergence.

For CBHM, the tuning parameters in equation (3.12) are set to be a = 0 and b = 0.2. These values are picked based on the data observations and they yield the best overall performance.

For MFM, the truncated Poisson parameter is set to be λ = 4 and the Dirichlet parameter is set to be γ = 100. This reflects the original belief that the basket are separately analyzed, which essentially helps reduce FWER in the simulation. To determine the values of ϕ and ξ, we run the simulation for various combinations and select the following values that provide the best balance between FWER inflation and power: ϕ = 0.1 and ξ = 0.16. In particular, some settings actually provide even a larger power but those are not chosen because they also have a larger FWER. The results of the sensitivity analysis for other combinations of ϕ and ξ are reported in Tables S.4 to S.6 for models without covariates and Tables S.7 to S.9 for models adjusting for covariates in Section S.2 of the Supplementary Material.

All the methods are compared in One-Stage Design (i.e., no interim analysis) and two-stage designs (i.e., with one interim analysis). For two-stage design, we set q 2 = 1 − η ( 1 − q 1 ) as one of the conditions to determine q 1 and q 2 . When η = 1 2 , both early efficacy and futility rules are adopted; when η = 0, no early stopping for efficacy. For a fixed value of η, the value of q 1 (and hence q 2 ) is determined via simulation to control the FWER under the global null to be no more than 10%. The values of q 1 under different stage design are provided in Table 2 for Models adjusting for covariates; in Table S.2 for Models without covariates.

In this section, simulation results for the models adjusting for covariates when the covariates are correctly identified are presented. For simulation data without covariates, the simulation results are presented in Section S.1 of the Supplementary Material. The results for models with covariate adjustment and models without covariates are similar. This is because both models are “correctly specified”, i.e. the inclusion of covariates is consistent with how the data were generated. to assess the performance of these methods, the 6 scenarios summarized in Table 1 are applied. In Scenario 1, the treatment is ineffective in any of the baskets; In Scenario 2, the treatment is effective for basket 1 only; In Scenario 3, the treatment is effective for baskets 1 and 2; In Scenario 4, the treatment is effective for baskets 1, 2 and 3; In Scenario 5, the treatment is effective for all baskets. The effect size of 0.4 is equivalent for all baskets in Scenarios 2 to 5; In Scenario 6, the treatment is effective for baskets 2, 3 and 4 with effect sizes of 0.2, 0.4, and 0.6, respectively. Scenario 1 is used to determine the thresholds for all models to have the FWER of less or equal to 10% at the end of the final stage under the global null.

In the above simulations, we assume that the covariates are strong predictors of the outcomes and need to be adjusted in the model. The following tables reanalyze the simulated data in the presence of covariate effects, using a model with no adjustment for covariates. Since the effects of covariates are not included in the analytical model, the independent simulations used to determine the tuning parameters also assume no covariate effect, i.e. the tuning parameter values in these additional simulations are the same as those in Table 2.

This simulation is referred to “Mis-specified Model 1” and their results are summarized in Table 6. As expected, the performance of all methods is affected by the mis-specified analytical model. The proposed methods still provide significant power gains in most scenarios, but all methods have greatly increased FWER inflation. The FWER inflation for MFM is drastic. The effect on SEP seems to be minimal, however, it may also depend on the magnitude of the covariates.

Another type of Mis-specified Model is those adjusting covariates of which the outcome is independent. These are referred to as “Mis-specified Model 2” and their results are summarized in Table 7. Our simulation results show that including negligible covariates does not impact the performance of our methods. Therefore, it is advantageous to include at least some potential predictors as the covariates.

As an extension our current simulation results, additional simulation is done with the doubled sample size of the treatment arm. In one stage, sample sizes for treatment arm are assumed to be n t 1 = 60, n t 2 = 60, n t 3 = 40, and n t 4 = 40; sample sizes for control arm are assumed to be n c 1 = 30, n c 2 = 30, n c 3 = 20, and n c 4 = 20. Their results are summarized in Table 8, Table 9 and Table 10. A larger cross model performance difference is observed. As expected, there was a significant power increase for all methods. For CBH and BHM, the inflation of type I errors is moderate, if not severe, as a cost of increased power. Similar to previously reported simulation results, implementing the early efficacy rule helps Bayesian methods in terms of power increase and Type I error control.