{'fraction_non_alphanumeric': 0.02961135101789019, 'fraction_numerical': 0.009545764472872496, 'mean_word_length': 4.76455174995321, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New designs for Bayesian adaptive cluster-randomized trials

Adaptive approaches, allowing for more flexible trial design, have been proposed for individually randomized trials to save time or reduce sample size. However, adaptive designs for cluster-randomized trials in which groups of participants rather than individuals are randomized to treatment arms are less common. Motivated by a cluster-randomized trial designed to assess the effectiveness of a machine-learning based clinical decision support system for physicians treating patients with depression, two Bayesian adaptive designs for cluster-randomized trials are proposed to allow for early stopping for efficacy at pre-planned interim analyses. The difference between the two designs lies in the way that participants are sequentially recruited. Given a maximum number of clusters as well as maximum cluster size allowed in the trial, one design sequentially recruits clusters with the given maximum cluster size, while the other recruits all clusters at the beginning of the trial but sequentially enrolls individual participants until the trial is stopped early for efficacy or the final analysis has been reached. The design operating characteristics are explored via simulations for a variety of scenarios and two outcome types for the two designs. The simulation results show that for different outcomes the design choice may be different. We make recommendations for designs of Bayesian adaptive cluster-randomized trial based on the simulation results.


Introduction
Randomized controlled trials, which can ensure that subjects assigned to each treatment group are comparable with respect to all characteristics of interest to draw a causal conclusion, have played an essential role in evaluating the effectiveness of interventions [1]. In most trials designed to assess the effect of a drug or a treatment, individual participants are randomly allocated to each treatment arm. However, some interventions naturally operate at a group level, or target either a social network or physical environment. For example, in a trial assessing the clinical utility, safety, and potential effectiveness of a machine-learning based clinical decision support system (CDSS) developed by Aifred Health [2][3][4][5], the decision support tool is designed for physicians, naturally forming clusters of individual patients being treated by that physician [6,7]. For these types of interventions, a population-level effect is of more interests to researchers.
In addition, randomizing by individuals in this setting may lead to a contamination effect between trial arms, as physicians might struggle to treat patients from different treatment arms strictly differently. Contamination can cause dilution bias and affect the reliability and validity of the study. One way to reduce the possibility of contamination is to randomize by physicians which act as 'clusters' [8,9]. Thus, it may not be advisable to randomize individual participants to different treatment arms, and groups of subjects being treated by their clinicians can instead be randomly assigned to the treatment arms in what is known as a cluster-randomized trial [9][10][11][12].
Due to the difference in randomization unit, the design and analysis of cluster-randomized trials differ from individually randomized trials [13][14][15]. The independence assumption between participants' outcomes is violated in a cluster-randomized trial since subjects from the same cluster tend to have more similar responses than subjects from a different cluster. Correlation between subjects within the same cluster, as measured by intra-cluster correlation coefficient (ICC), should be considered carefully when analyzing data from cluster randomized trials. The correlation reduces the variability of responses in a clustered sample and thus reduces the statistical power to detect true differences between treatment arms relative to trials that randomize the same number of individuals [16,17].
A natural way to improve trial efficiency is to incorporate adaptive features into the design which allow for planned adjustments to the trial design after its initiation without undermining the validity and integrity of the trial [18]. In addition to flexibility and efficiency, adaptive designs are attractive to clinical scientists because they may reflect medical practice in the real world and they are ethical with respect to the need to determine and monitor efficacy as well as safety of the treatment [19]. Commonly used adaptive designs in clinical trials include, but are not limited to, adaptive randomization, group sequential design, and stopping rules [19,20].
Adaptive designs naturally fit into the Bayesian framework as that results or estimates are continually updated based on the accumulated information from interim data [21]. Adaptive designs based on Bayesian approaches have been extensively studied in recent years [22][23][24][25][26][27]. However, most adaptive designs, regardless of whether frequentist or Bayesian, focus on individually randomized trials. Adaptive designs for cluster-randomized trials are less common and the incorporation of adaptive features poses significant statistical challenges. Some specific adaptive features such as sample size re-estimation and group sequential design have been proposed in combination with clusterrandomized trials [28][29][30][31]. However, no formal statistical design of Bayesian adaptive cluster-randomized trials has been developed. We address this gap by proposing two Bayesian adaptive designs for cluster-randomized trials.
The organization this paper is as follows. The motivating trial is briefly described in section 2, followed by the development of two Bayesian adaptive designs and models for continuous and binary outcomes in section 3. Simulation studies are carried out in section 4, assessing the performance of the two proposed designs across a range of scenarios. The paper concludes in section 5.

Motivating setting
Aifred Health [2][3][4][5] has designed a machine-learning based clinical decision support system (CDSS) for physicians treating patients with depression. Patient characteristics of interests such as sociodemographic information, clinical information and medical history are input into the CDSS. The CDSS then, using a deep learning model, outputs the predicted efficacy for a number of possible treatments for that patient [6,7]. Treatments are ordered by efficacy and presented to the physician when they reach the treatment selection step of a clinical algorithm based on best practice guidelines [3,32]. Physicians with the CDSS can, on an individual patient basis, decide whether or not to use the information presented by CDSS as part of their medical decision-making.
An important step is to establish the clinical utility, safety, and potential effectiveness of a tool such as the CDSS. Practically, an intervention such as the CDSS must be delivered at the physician level, such that a cluster-randomized design would be appealing for the reasons described above. Participating physicians could be randomized, with patients recruited from physicians' usual practices in order to approximate the real-world clinical conditions and populations as closely as possible.
Typical outcomes in a study of depression in such a trial could be a continuous measure of symptoms (e.g., a visual analog scale, the Quick Inventory of Depressive Symptomatology [33], the 9-question depression scale from the Patient Health Questionnaire [34], the World Health Organization Disability Assessment Schedule 2.0 [35], etc.) or a binary measure of treatment response, minimum clinically significant change in symptoms, or remission, perhaps defined by a dichotomization of a standard depression score [36][37][38]. The importance of mental health treatment and the often relatively slow rate of patient accrual motivate the use of a Bayesian adaptive trial design to ensure adequate sample size and the possibility of early termination due to treatment effectiveness.

Two Bayesian adaptive cluster-randomized designs
Suppose that the maximum number of clusters and maximum cluster size are equal across the two treatment arms. Let K, n, m be the number of interim analyses (not including the final analysis), the maximum number of clusters for each treatment arm, and the maximum cluster size, respectively. For simplicity of exposition, we will assume that all clusters enroll the same number of participants. Two Bayesian designs, design 1 and design 2 in the remainder of this paper, are developed to sequentially enroll participants and analyze interim data in different ways.
In design 1, [ n K+1 ] clusters enter the trial at the start where [x] denotes the largest integer not exceeding x, and m individual participants will be enrolled for each cluster. At the subsequent analysis point, if the accumulated information up to the current analysis is sufficient to conclude the efficacy of the treatment or the final analysis has been reached, the trial will be terminated. Otherwise, another [ n K+1 ] new clusters will be enrolled, and m individual participants will be recruited for each new cluster. The trial then proceeds to the next analysis with new samples. This procedure is repeated until termination.
In design 2, at the beginning of the trial, all n clusters enter the trial, but only [ m K+1 ] individual participants are enrolled for each cluster. At the subsequent analysis point, if the trial is not to be terminated, another [ m K+1 ] individual participants are recruited for the same n clusters. The trial then proceeds to the next analysis. This procedure is repeated until termination.
The fundamental difference between the two designs lies in the way that participants are sequentially recruited. In design 1, the clusters are sequentially enrolled, and individual participants for each cluster are recruited all at once. In design 2, all clusters are enrolled at one time, but the individual participants for each cluster are sequentially enrolled. For illustration, consider an example for one treatment arm with K = 1, n = 4 and m = 6. That is, there is only one interim analysis planned (thus two analyses in total, including the final analysis). The maximum number of clusters is 4, and the maximum cluster size is 6. Figure 1 gives a graphical illustration. The black labelled circles represent different individual participants for the corresponding clusters. In design 1, clusters 1 and 2 will be first enrolled, and for each, six individual participants will be recruited. The interim analysis is based on data from clusters 1 and 2. If we decide to continue the trial, then we will further recruit clusters 3 and 4 and their respective six individual participants for the final analysis. The final analysis is based on data from all the four clusters. In design 2, all four clusters will be recruited at the start, but for each cluster, only three participants will be enrolled. If evidence based on the 12 individuals from 4 clusters is unable to conclude the efficacy of the intervention at the interim analysis, then an additional three individual participants for each cluster will be recruited, and the trial proceeds to the final analysis.

Early stopping at interim analysis
At each interim look, one should determine whether to stop the trial early or continue based on the interim result. The goal is to evaluate the efficacy of the treatment by testing the hypothesis where θ is the mean difference for continuous outcome and risk difference for binary outcome. The efficacy of the treatment can be concluded and the trial can be stopped early for efficacy at the k-th interim analysis if where δ is the minimal important difference, D k is all the available data up to the k-th analysis and U is the decision boundary.

Continuous outcomes
Based on the context of the motivating trial and without loss of generality, assume that a smaller value of the continuous outcome is preferred (as most depression and disability rating scales associate worse symptoms with larger total scores). Then, θ = µ c − µ t where µ c , µ t are the mean outcome for control and treatment groups, respectively. We further assume that at the kth analysis, k = 1, . . . , K + 1, there are n k clusters with m k observations in each cluster. Let Y ij be the continuous outcome of the i-th subject in the j-th cluster at current analysis point, j = 1, . . . , n k , i = 1, . . . , m k . The normality assumption in cluster-randomized trials states where µ j is the cluster-specific mean, µ is the population mean, σ 2 W and σ 2 B are within-and between-cluster variances respectively, and they can be related via the ICC, ρ = . It can be shown that the marginal distribution of Y ij is normal with mean µ and variance σ 2 . . , n k , and the covariance structure satisfies equal to µ and Σ is a block matrix of the form (2), and MVN indicates a multivariate normal distribution.
In the Bayesian framework, assume a normal prior for µ, where a is the prior mean and b 2 is the prior variance. Then, using Bayes' theorem, the posterior distribution for µ can be obtained: The above result applies to both µ t and µ c . Therefore the quantity of interest P (µ c − µ t > δ|Y ) can be estimated bŷ where µ post ci , µ post ti are sampled from the corresponding posterior distributions of µ c and µ t , and M is the number of Monte Carlo samples drawn from the posterior.
At any interim analysis, ifπ > U , the trial is stopped early for efficacy. Otherwise, the trial continues to enroll clusters/participants and proceeds to the next analysis where the prior mean and variance for the next analysis are updated with the posterior mean and variance for the current analysis. These steps are repeated until either the trial is stopped early for efficacy or reaches the final analysis.

Binary outcomes
For binary outcomes, assume that a larger proportion is preferred (e.g., a larger proportion of patients meeting the criteria for treatment response). Let θ = π t − π c where π c , π t are the population proportion for the control and treatment groups, respectively. For binary outcomes, the joint distribution of all observations cannot be obtained analytically and thus a tractable form of the posterior distribution of the parameters of interest is not available. Therefore, a hierarchical model is proposed as follows where r j is the number of events in the j-th cluster and π j is the clusterspecific proportion, j = 1, . . . , n k . Then, under the model, the cluster-specific proportion π j have the mean-variance relationship However, in cluster-randomized trials with binary outcomes, it is assumed that where π is the population proportion and ρ is the ICC. To make (3) and (4) consistent, define two transformed parameters where π is exactly the mean of the Beta distribution and v measures the information in the corresponding Beta distribution. Also, due to the consistency of (3) and (4), once ρ is fixed or can be estimated, v can also be determined through v = 1−ρ ρ . Thus, the only free parameter for the hierarchical model is π. In the Bayesian framework, assume a uniform prior for π, Then the posterior probability P (π t −π c > δ|R) can be approximated by drawing posterior samples using Markov Chain Monte Carlo (MCMC) implemented in, say RStan [39,40], where R = (r 1c , · · · , r n k ,c , r 1t , · · · , r n k ,t ) includes the number of events per cluster within each treatment group.

Simulation studies
The false positive rate and power cannot be obtained analytically in Bayesian adaptive trials since the sampling distributions of the test statistics (i.e., posterior probability statements (1) in section 3.2) are not known. Therefore simulation studies are required to specify the decision boundaries and other design characteristics [22]. For both outcomes, a single interim analysis was explored first; two interim analyses were then investigated. The minimal important difference was set as 0 in the simulation but it is straightforward to generalize to other values. For each scenario, 500 simulation replications were performed. The performance of designs was compared based on false positive rate and power. The false positive rate is estimated as False positive rate = Number of times the null hypothesis is falsely rejected Number of simulation runs , where falsely rejecting the null hypothesis means that for some k ≤ K + 1, forθ k the estimated mean difference for continuous outcome or estimated risk difference for binary outcome at the k-th analysis, θ is the true mean or risk difference, δ is minimal important difference, and δ = 0 in our simulation. D k is all the available data up to the k-th analysis, and U is the decision boundary as described in section 3.2.
Power is estimated as Power = Number of times correctly detecting the difference Number of simulation runs where correctly detecting the difference means that for some k ≤ K + 1, for θ 0 the value of θ under the alternative hypothesis. For continuous outcomes, the prior mean and variance for the population mean for both groups are fixed at 0 and 100, respectively. Various values of ICC were explored and the between-cluster variance σ 2 B was determined through To generate clustered continuous data, we first generate the n cluster-specific means from normal distributions with mean µ c (control group) and µ t (treatment group) and variance σ 2 B . Then, within each cluster, m samples are drawn from a normal distribution with mean equal to the clusterspecific mean and variance σ 2 W . The resulting m × n samples are expected to satisfy the preset correlation structure. For binary outcomes, clustered binary data are generated via Beta and binomial distributions; see Appendix C for details. All simulation parameters are summarized in Table 1.
Only results for continuous outcomes are presented here. The results for binary outcomes are summarized in Appendix C. Figures 2 and 3 show the false positive rate and power when a single interim analysis is planned. For power, only results for ρ = 0.2, 0.5, 0.8, representing low, moderate and high correlation, respectively, are displayed. In general, design 2 has higher power but also higher false positive rates compared with design 1.
However, in all scenarios, with U = 0.95, the false positive rates for both designs are well above 0.05 and thus are unsatisfactory. With U = 0.98, the false positive rates are reduced. Especially for design 1, the false positive rate can be controlled within 0.05 with U = 0.98. Therefore, U = 0.98 serves as a better decision boundary than U = 0.95. A larger decision boundary corresponds to a more conservative test since more evidence is required to conclude the treatment efficacy.
The effect of ICC on false positive rate is quite small but the effect on power is more direct. Under the same conditions, if the underlying ICC is higher, power is lower. For small ICC, if the underlying treatment effect is moderate to large, then only 20 clusters for each group will be sufficient to ensure the desired power, but if the underlying treatment effect is small, larger number of clusters, for example 40 or 60, should be considered. For moderate ICC and a small to moderate treatment effect, a larger number of clusters is recommended; for a large treatment effect only a small number of clusters, for example 20 clusters for each group may be enough. For large ICC, it is recommended to have more clusters than 60 as, even with a large treatment effect, power is still low.    Figures 4 and 5 show the false positive rate and power when multiple interim looks are built into the study design. Design 2 still has higher false positive rate and power compared with design 1. As was observed in Figures  2 and 3, a larger decision boundary can reduce false positive rate and the resulting reduction in power can be remedied by recruiting a larger number of clusters. Adding more interim analyses may increase the false positive rate and the power. Given a fixed decision boundary, with more interim analyses planned, there is higher chance of rejecting the null hypothesis.
As seen in the plot, for large ICC values, the addition of multiple interim analyses can bring obvious improvement in power under design 2. However,larger sample size is required for large ICC values to achieve sufficient power. For small to moderate ICC, regardless of design, a single interim analysis is preferred. The effect of cluster size on both false positive rate and power is very modest. The results for a larger cluster size are displayed in Figures B1  and B2 in Appendix B. To conclude, both designs perform better in terms of false positive rate with U = 0.98, and design 2 has a higher false positive rate and power. The choice between the two designs as well as the two decision boundaries depends on the research question, phase of the study, and the information we have about the treatment. For example, in some early phase clinical trials, if a relatively high false positive rate is acceptable, then design 2 with U = 0.95 may be recommended. However, for the situation we explored, a single interim analysis is preferred no matter which design or decision boundary is chosen.

Discussion
Motivated by a potential real-world cluster-randomized trial, we explore the statistical properties of Bayesian adaptive cluster-randomized trial designs. We explore stopping rules for efficacy in a cluster-randomized trial. Interim analyses may be planned over the course of the trial and at each interim analysis the trial may be stopped early if sufficient evidence is established to conclude efficacy otherwise the trial proceeds to the next analysis. We proposed two designs which sequentially enroll participants in different ways. The first design sequentially enrolls clusters, and individual participants for each cluster are recruited all together. The second recruits all clusters at the start of the trial, then sequentially enrolls batches of participants. Regardless of the design choice, the data analysis procedure is the same: the difference between the designs lies in the sequential enrollment of participants. Interim analyses of both continuous and binary outcomes are performed based on common models. For continuous outcomes, on the basis of the normality assumption, we obtained the analytical form of the posterior distribution of the population mean. Based on Monte Carlo simulation, the posterior probability of efficacy can be easily estimated by drawing random samples directly from the two posterior distributions, one for each treatment group. For binary outcomes, due to the complex correlation structure, the posterior distribution of population proportion cannot be obtained analytically. Instead, a hierarchical model was used and the posterior probability of efficacy is estimated via MCMC.
Through simulation, we found the design choice may be different for different outcomes. The preferred design may also depend on the research goal, phase of study, or feasibility considerations. For binary outcomes, within the parameter space we explored, design 2 is recommended based on design operating characteristics. Also, for both outcomes, one interim analysis may be recommended for the situations we investigated, as it makes the designs more feasible and the design performance with single interim analysis is satisfactory. However, a more general recommendation requires further exploration to account for practical issues or other design parameters we did not explore in our simulation. In addition, increasing cluster size may not necessarily bring much improvement in design performance. Keeping a small cluster size may not hurt design operating characteristics very much but it will save time and cost. In our motivating trial, a smaller cluster size is also more realistic, as it reduces the burden on the individual clinician to find suitable patients from their practice.
There are some limitations to our work. First, we only considered implementing a stopping rule for efficacy. There are many other adaptive features that we have not considered. One example is adaptive randomization. The response-adaptive randomization [41] has not been extensively investigated from Bayesian viewpoint. Second, we only focus on two-arm trials. Extension to multi-arm trials which may involve arm dropping can be considered in the future. Third, we only considered continuous and binary outcomes. However, survival endpoints are also prevalent in clinical trials. An extension to survival outcomes based on survival models can also improve the framework of Bayesian adaptive cluster-randomized trial. In addition, the effect of unequal cluster sizes has not been explored for our designs. However, the effect of varying cluster sizes on design operating characteristics in standard clusterrandomized trials has been extensively discussed [14,[42][43][44][45]. Typically, designs are most efficient for equal cluster sizes. Inflation of false positive rate will occur for imbalanced studies, and with the same sample size, imbalanced trials may be underpowered compared with their balanced counterparts. Similar impacts of unequal cluster sizes on power or false positive rate may also exist for our Bayesian adaptive cluster-randomized trials, and this would need to be determined in future work. Also, practical issues in planning a Bayesian adaptive cluster-randomized trial are not taken into account in our paper. For example, sometimes it may not be realistic to plan the preferred design recommended by design operating characteristics. In this case, feasibility may be the main reason for design choice. These context-specific issues require more exploration in the future.

Declaration of competing interest
D.B. is a director, shareholder, and employee of Aifred Health.

Appendix A Propositions
Proposition 1 For Y ij satisfying Y ij |µ j ∼ N (µ j , σ 2 W ) and µ j ∼ N (µ, σ 2 B ), the marginal distribution of Y ij is normal with mean µ and variance σ 2 Proof The marginal distribution of Y ij can be obtained from We can calculate the integral using the pdf of normal distribution, Therefore, the marginal distribution is

Proposition 2
The response vector Y = (Y T 1 , · · · , Y T n k ) T where Y j = (Y 1j , · · · , Y m k j ) T , j = 1, · · · , n k , and Y ij satisfies the conditions stated in proposition 1, follows a multivariate normal distribution with mean − → µ = (µ, . . . , µ) m k ×n k and variance-covariance matrix Σ as a block matrix of the form Proof We only need to prove that ∀ (a 11 , a 21 , · · · , amn), Z = a 11 Y 11 + · · · + amnYmn is univariate normal. Let M Z (t) be the moment generating function of Z. Therefore, This is the moment generating function of a normal random variable. Therefore, Z is univariate normal and Y is multivariate normal. Then the mean of Y is − → µ = (µ, . . . , µ) T and the variance-covariance structure satisfies N (a, b 2 ) for µ, under the conditions in proposition 1 and proposition 2, the posterior distribution of µ is

Proposition 3 Assume a prior distribution
Proof From proposition 2, the response vector Y ∼ MVN( − → µ , Σ). Therefore, we have That is,

Appendix C Simulation for binary outcomes
To generate clustered binary data, we first generate n cluster-specific proportions from Beta distributions with Beta parameters determined by the preset population mean and ICC so that all cluster-specific proportions are strictly between 0 and 1. Then, within each cluster, the number of events were generated from a binomial distribution with the cluster size at current stage and the cluster-specific proportions. Then the resulting m × n binary variables satisfy the predetermined correlation structure. The performance of the two designs in terms of false positive rate and power, when only one interim analysis is planned, is displayed in Figures C3  and C4. It can be observed that design 1 has higher false positive rate and both designs perform almost equally well in terms of power over the parameter space we explored. Also, a decision boundary of 0.98 helps control false positive rate compared with U = 0.95 while the decrease in power coming with a larger decision boundary is acceptable. However, for design 2, even with a smaller decision boundary, the false positive rates for most cases are acceptable considering that they are fluctuating around 0.05. In addition, as in the continuous case, when the underlying ICC is larger, power is relatively small. The design performance for multiple interim analyses is shown in Figures  C5 and C6. Design 1 still has higher false positive rate and the powers for the two designs are similar over the parameter space we explored. With multiple interim analyses, false positive rate will increase for most cases. The increase in power, when n = 20 and ∆ = 0.1, outweighs the increase in false positive rate. However, in this setting the powers for both designs are still unsatisfactory with each possible combination of ICC, decision boundary or baseline risk we explored. It suggests that when the effect size is 0.1, 20 clusters for each treatment arm may not be sufficient to ensure the desired power. When we have more clusters, such as 40 or 60, with multiple interim analyses, the increase in power is relatively slight compared with the increase in false positive rate for most cases. When the effect size is 0.2, a smaller sample size is sufficient to ensure the desired power. Therefore, the power is satisfactory for almost all scenarios we investigated. Similarly, with multiple interim analyses a larger boundary value can evidently help reduce false positive rate, especially for design 1. The resulting decrease in power is acceptable for a sufficient sample size. However, with a small sample size, a larger decision boundary may lead to insufficient power. The results for a larger cluster size are also displayed in Figures C7 and C8. Therefore, for binary outcomes, for the parameters explored in the simulation, design 2 is recommended since it has smaller false positive rate and performs comparably with design 1 in terms of power. Additionally, a single interim analysis is sufficient to control false positive rate while maintaining satisfactory power. Regarding the decision boundary, if design 2 with single interim analysis is planned, U = 0.95 may be conservative enough for obtaining a satisfactory false positive rate as well as power. However, if due to feasibility, design 1 is preferred or multiple interim analyses are of interest a larger decision boundary such as U = 0.98 may be recommended. Developing a more general set of recommendations requires further exploration.