Meta-Analysis with Few Studies and Binary Data: A Bayesian Model Averaging Approach

: In meta-analysis, the existence of between-sample heterogeneity introduces model uncertainty, which must be incorporated into the inference. We argue that an alternative way to measure this heterogeneity is by clustering the samples and then determining the posterior probability of the cluster models. The meta-inference is obtained as a mixture of all the meta-inferences for the cluster models, where the mixing distribution is the posterior model probabilities. When there are few studies, the number of cluster conﬁgurations is manageable, and the meta-inferences can be drawn with BMA techniques. Although this topic has been relatively neglected in the meta-analysis literature, the inference thus obtained accurately reﬂects the cluster structure of the samples used. In this paper, illustrative examples are given and analysed, using real binary data.

where for each study i ∈ {1, . . . , k}, the observed effect x i is assumed to be normally distributed with mean parameter θ i (the treatment effectiveness conditional on study i) and variance τ i . Similarly, θ represents the unconditional treatment effectiveness, i.e., the pooled effect size; τ is its variance; and [−, −] indicates a prior density to be assigned. When the effectiveness is measured by a discrete 0-1 random variable, the above hierarchical normal model is typically applied to the logit transformation of the data log (x i /(n i − x i )) with the reparametrisation log (θ i /(1 − θ i )) [3][4][5][6][7][8]. If the samples contain zeros, a fixed data continuity correction is applied, and a logit approximation is then used. However, this practice has been criticised by Sweeting et al. [9], who proposed an alternative empirical data correction. Furthermore, small studies present more heterogeneity than large ones [10]. Recently, Weber et al. (2020) [11] recommended that zero-cell corrections should be avoided due to the poor performance of their statistical properties. Obviously, continuity corrections are not the only way to proceed. A number of meta-analysis methods based, e.g., on binomial models were proposed in [12] and the references therein.
In the meta-analysis literature, determining statistical heterogeneity mainly consists of either estimating the between-study variation, characterised as the variance, or heterogeneity parameter, or otherwise testing the null hypothesis that the true treatment effects across studies are identical H 0 : θ 1 = θ 2 = . . . = θ k = θ [8,13]. A measure of evidence for heterogeneity is obtained by means of the Q test, a chi-squared test with k − 1 degrees of freedom, first defined by Cochran [14]. Alternatively, and in order to overcome the dependence on the number of studies considered for the meta-analysis, which is the case with the Q test, Higgings and Thompson [15] proposed the dimensionless I 2 index, which describes the percentage of variation across studies that is due to heterogeneity rather than chance. Mittlböck and Heinzl [16] compared these two relative measures of heterogeneity jointly with the index H 2 M = Q k − 1 − 1, through an intensive simulation study. Basically, I 2 and H 2 M have similar properties, although H 2 M is preferred when few studies are considered. As shown in [16], the number of studies affects the power of the heterogeneity test and the heterogeneity measures I 2 , but not H 2 M . It is well known that all tests of heterogeneity have relatively little power to detect heterogeneity when data are sparse and/or when meta-analyses are based on a small number of studies [1,9,15,17]. These situations are very common in daily practice, for instance with respect to rare or orphan diseases, where zero values are often observed in both the trial and the treatment arms considered.
Statistical heterogeneity is a characteristic of the studies that compose a meta-analysis and can be present at different levels of intensity, ranging from homogeneity to heterogeneity. There are also many intermediate situations that should be considered in order to correctly estimate the quantities of interest. Higgins et al. [18] proposed a tentative naive categorisation of values for I 2 , relating values of 25%, 50% and 75% with low, moderate and high levels of heterogeneity, respectively, and considering a significant degree of heterogeneity to be present when I 2 > 50%.
In the present analysis, we address the question from a different standpoint, arguing that between-sample heterogeneity is a clustering problem and that model uncertainty can be incorporated into the inference using a Bayesian procedure. Under this procedure, the posterior probabilities of the cluster models are computed, and the meta-inference is then obtained as a mixture of the meta-inferences obtained for the cluster models, using the posterior model probabilities as a mixing distribution. Another point of difference is that most previous studies analysed between-study heterogeneity for a quantity of interest in the comparison of two treatments, such as the mean difference, relative risk or odds ratio. In our proposal, the configuration of between-study heterogeneity is analysed for each treatment separately, which allows the heterogeneity structure to differ between the experimental treatment and the control.
The proposed methodology is known as Bayesian model averaging and has recently begun to be used in meta-analyses. For instance, Bayesian meta-analyses modelling averages over fixed and random effect models were considered in Gronau et al. [19] and Scheibehenne et al. [20].

A Motivating Example
Günhan et al. [12] suggested the use of weakly informative priors for the treatment effect parameter of a Bayesian meta-analysis model, to be applied in a paediatric transplant dataset. In the present study, and for illustrative purposes, we focus on a dataset corresponding to renal post-transplant lymphoproliferative diseases (PTLD). The results obtained are shown in Table 1. In this respect, Crins et al. [21] conducted Cochrane's Q test and found no evidence of heterogeneity between the trials. Moreover, Günhan et al. obtained low values for the estimated heterogeneity parameter under weakly informative or vague priors and maximum likelihood estimation. However, the structure of the partition clustering induced by the PTLD dataset in Table 1 shown in Figure 1 indicates that different forms of inter-trial variability could be present.
Homogeneity Type 2-Heterogeneity Heterogeneity There are |M| = 5 possible partitions (models or heterogeneity configurations), and for each of the cluster classes, there are 1, 3, 1 partitions, respectively. Thus, the selection of any one model (homogeneity, for example) among the five that are possible may ignore the uncertainty underlying model selection and misrepresent the uncertainty concerning the quantities of interest. In view of these considerations, we believe that Bayesian model averaging (BMA) is an appropriate means of accounting for model uncertainty and, thus, heterogeneity. Expectations and quantities of interest are obtained by weighted averages over the set of models, rather than by selecting a single model and drawing inferences as if it were the true model.
For each of the treatments considered (control and experimental), Table 2 shows the (posterior) probabilities of each of the five models induced by the clustering applied to the data in this example. These posterior probabilities constitute the weights to be used for averaging the set of models. Their expressions are given in Section 3. As an example, let us observe the control treatment. The case of homogeneity has the highest posterior probability of being true, 0.39. However, this option does not consider the remaining heterogeneity structures and therefore would result in misleading inferences being drawn, by rejecting all other possible models, which constitute a combined probability of 0.61. The same reasoning is applicable to the experimental treatment.
The BMA approach to meta-analysis involves averaging over all the possible models when making inferences about the quantities of interest. For example, for the dataset considered in this case study, if θ is the treatment effect under the control treatment, its posterior distribution, for the data given, is: i.e., the average of the posterior distributions under each of the models ("type of heterogeneity") considered, weighted by its posterior model probability. For instance, the posterior mean for θ is a weighted average of the posterior means in each model: where the posterior probability for model M m can be obtained by: where:

Summary
The rest of this paper is organised as follows. The binomial Bayesian models and the conditional distributions linking the experimental parameters θ i and the meta-parameter θ are presented in Section 2, where the Bayesian procedure for clustering the samples and the likelihood of the meta-parameter are also given. The Bayesian model averaging procedure for estimating the meta-parameters is described in Section 3. Three illustrative examples are provided in the following section, with real datasets. Finally, Section 5 summarises the main conclusions drawn and presents some concluding remarks.

The Bayesian Binomial Model
Assume a meta-analysis involves k studies that provide k independent discrete samples with the binomial distribution {Bin(x i |n i , θ i ), i = 1, . . . , k}, where θ i represents the treatment effectiveness, n i the number of patients and x i the number of successful treatments, conditional on the study i. Assume, moreover, that the prior information on the conditional treatment effectiveness θ i is weak. Accordingly, the uniform prior Unif(θ i |0, 1) is used in a context where the data contain zeros [22,23]. Therefore, for i = 1, . . . , k, we consider the Bayesian sampling models: where: and 1 A is the indicator function. The value of one is assigned to all elements of A, and zero elsewhere. In this analysis, it is convenient to consider a 0 − 1 latent variable x, the meta-variable, which is defined as the result obtained when a treatment with meta-effectiveness θ is applied to a patient in a virtual study, which is not affected by between-study variability. The distribution of this meta-variable x is of the same type as that of x i , and hence, we have the Bernoulli meta-model Ber(x|θ), where the meta-parameter θ represents the true (unconditional) treatment effect. This meta-model gives a precise meaning to the meta-parameter θ. The objective Bayesian meta-model M is then given by: The Bayesian meta-analysis is then based on the posterior distribution of the parameter θ, which is given by: The likelihood function in (9) cannot be obtained with the information on sample x i on study i, which is related to θ i , but not to θ. Therefore, further steps are required to derive an appropriate likelihood function for θ: 1. A distribution π(θ i |θ) is needed to link the experimental parameters θ i and the meta-parameter θ. This linking distribution must ensure there is coherence between the conditional and marginal distributions of the experimental parameter and the meta-parameter. Mathematically, this requires that the corresponding bivariate distribution belong to the class of bivariate distributions with given marginals. 2. It is clear that the likelihood of θ strongly depends on the cluster structure of the samples.
Therefore, before formulating the likelihood of θ for the samples x i , i = 1, . . . , k, we must address the important question of whether k, the dimension of the model, can be reduced by clustering homogeneous samples. 3. In other words, each cluster model indicates a different heterogeneity structure of the sampling model for x, and the posterior probability informs us about the uncertainty for this structure.
For the data x = (x 1 , . . . , x k ), we need to obtain the likelihood of θ conditional on a given cluster model. Finally, the likelihood of θ for the data x is obtained as a mixture of the above conditional likelihood functions.
We now proceed step-by-step.
Following Moreno et al. [24], we consider the conditional intrinsic linking distributions {π I (θ i |θ, t), t = 1, 2, . . .} obtained by comparing model M with model M i , which presents interesting properties as a linking distribution. For any positive integer t, the intrinsic method gives the conditional intrinsic prior as the following mixture of Beta distributions, From (11), we obtain: and hence, the hyperparameter t indicates how strongly the conditional distribution π I (θ i |θ) concentrates mass around θ. In practice, hyperparameter t is fixed. Hence, for the sake of simplicity in notation, we refer to the linking distribution π(θ i |θ) rather than π(θ i , |θ, t).
The bidimensional prior π I (θ i , θ) = π I (θ i |θ)1 (0,1) (θ) satisfies Equation (10) for any t, and therefore, the linking class of intrinsic distributions and the Bayesian models (6) and (8) are coherent. In the following, we assume that θ i , i = 1, . . . , k, are conditional independent given θ. Therefore, the linking distribution of θ 1 , . . . , θ k conditional on θ is given by: Finally, observe that the proposed linking intrinsic distribution of θ i , conditional on θ, does not require the use of the logit transformation for the sparse data as is the case with standard random-effect models.

Clustering the Experimental Samples
Following Moreno et al. [25], we first define what is meant by cluster. The samples x i and x j , i = j, from f (x|θ i ) and f (x|θ j ), respectively, are said to be in the same cluster if θ i = θ j . The between-study heterogeneity is then determined by the number of clusters and by the location of the samples x 1 , . . . , x k within these clusters. The posterior distribution of the clusters is needed in order to draw inferences from these quantities.
To cluster the samples, we adopt the product partition model approach proposed by Barry and Hartigan [26], together with a Bayesian model selection procedure based on Bayes factors for the intrinsic priors for the model parameters and on hierarchical uniform priors for the cluster models.

The Likelihood of θ
From (7), given a partition r p = (r 1 , . . . , r k ), the sampling distribution of x = (x 1 , . . . , x k ) is: where θ p = (θ r 1 , . . . , θ r k ) is an unknown parameter of dimension p and the component θ j in (13) corresponds to r i = j, m j = ∑ i:r i =j n i and s j = ∑ i:r i =j x i . For example, the cluster model given by the heterogeneity structure (1, 2, 2) in Table 3 has the sampling distribution: The heterogeneity partition r k = (1, 2, . . . , k) has the corresponding likelihood function given by: Now, following Moreno et al. [25], integrating out θ p with the intrinsic prior π(θ p |p, r p ) = π I (θ 1 , . . . , θ p |θ)1 (0,1) (θ) dθ, we obtain the likelihood of θ, conditional on the cluster model (p, r p ). After some algebra and using the assistance Wolfram Mathematica, the likelihood of θ, conditional on the cluster model (p, r p ), is given by: where 3 F 2 (v, w, z) denotes the generalised hypergeometric function with argument z and vector parameters v and w of dimensions two and three, respectively, a j = (−t, −t, s j + 1) and b j = (1, −m j − t + s j ).
To derive the necessary likelihood function of θ, we then integrate out (14) with respect to a discrete prior on (p, r p ). As recommended by Moreno et al., an appropriate uniform hierarchical prior is used, in which each factorised prior is uniform.
For example, for the cluster configuration (1, 2, 2) in Table 3, the uniform hierarchical prior is: The complexity in this task lies in knowing how many elements are in each heterogeneity subclass. According to Casella et al. [27], who presented a comprehensive description of these calculations, the first product factor in (15) corresponds to the inverse of the multinomial factor, where k 1 ≤ k 2 ≤ . . . ≤ k p are integers assigned to cluster 1, . . . , p such that k 1 + . . . + k p = k. The second product factor takes into account the redundancies in the number of samples k 1 , . . . , k p , and b(k, p) is the number of these possible sizes of k j that satisfy the recursive equation b(k, p) = b(k − 1, p − 1) + b(k − p, p), 1 ≤ p ≤ k with b(k, 1) = b(k, k) = 1. For each value given for k and p, this number is readily obtained with the Mathematica package PartitionFunctionP. The third factor corresponds to the discrete uniform distribution of the number of studies in the meta-analysis.
Finally, from (14) and (15), the (unconditional) likelihood of θ for the data x is given by:

Bayesian Model Averaging in the Meta-Analysis
As mentioned in Section 1, the BMA approach to meta-analysis requires us to average all possible models (heterogeneity structures) when making inferences on θ. In this case, the posterior probabilities in (4) correspond to those of any heterogeneity structure given by a pair (p, r p ), which is represented by: where m r p (x|p, r p ) = f (x|p, r p , θ p )π(θ p |p, r p ) dθ p is the marginal of the data x conditional on model (p, r p ), with f (x|p, r p , θ p ) and π(θ p |p, r p ) as in (13) and (15), respectively. These posterior model probabilities Pr(p, r p |x) are the weights for the meta inference. The posterior distribution in (2) becomes: where, since as in (8) The posterior distribution in (18) is computed numerically over the parametric space Θ = (0, 1) using Wolfram Mathematica, which has a huge library of ready-to-use functions and, moreover, provides a simple code. Furthermore, once the posterior distribution has been obtained, the command ProbabilityDistribution can be used to generate any sample of the posterior distribution of the parameter of interest and its transforms, as shown in the examples given in Section 4.

Illustrations
To illustrate the arguments developed in the preceding section, we now analyse two real datasets. For these case studies, we assume that t = 48, indicating that the intrinsic link distribution concentrates a considerable mass of probability around θ. Other values of t have also been used, and the results obtained are very robust. The Mathematica code for the case study datasets is available from the authors upon request.

Motivating Example Revisited
Returning to the dataset considered in Section 1, from Table 2, observe that we have different posterior probabilities for the cluster structures and thus for the Bayesian model average of the true treatment effect under treatment (T) or control (C), θ T and θ C . All the measures of interest can be obtained from the BMA posterior distribution in (18). Table 4 shows the posterior summaries of the PTLD rate for both cases (mean and highest density interval (HDI)).  Figure 2 (left panel) shows the BMA posterior density of θ C and θ T given by (18). It is apparent that the posterior density of θ under the experimental treatment is slightly shifted towards higher PTLD rates with respect to the control treatment. Observe, moreover, that the sparsity of the data yields reverse J-shaped posterior distributions with a mode at θ = 0.
On the other hand, since we have an explicit expression for the posterior distribution of θ C and θ T from (18), it is straightforward to simulate these independent distributions. The usual parameters of interest, such as the risks ratio or the odds ratio and their log, are obtained immediately using the corresponding transform of the above-simulated values. For comparison, Table 5 shows the estimations of the risk and odds ratio from Crins et al. [21] and from Günhan et al. [12], together with those from our BMA approach. These estimated values show some differences from those obtained by Crins et al. and Günhan et al. Although Crins et al. asserted that there is no evidence for between-trial heterogeneity for PTLD outcomes, this conclusion does not follow from Table 2. Certainly, the homogeneity model has the largest posterior probability, for both treatment options, but a large proportion of the uncertainty accumulated in the remaining heterogeneity structures must be considered, in line with the BMA approach. This, together with the very small number of informative studies conducted in this respect, no more than three, accounts for the wider intervals obtained by the BMA procedure.

Another Case with Few and Double-Zero Studies
The dataset in Table 6 is extracted from Cosmi et al. [28] and corresponds to the comparison of total mortality in four studies (k = 4) of patients treated with ticlopidine (T) versus oral anticoagulation for coronary stenting (C). Cosmi et al. found no evidence of heterogeneity, and I 2 was estimated at 0.0%. This fact occurs frequently for meta-analyses with a small number of studies [29]. Furthermore, the Q test for this situation possesses a very low power [30]. As a consequence, we have little information on homogeneity (heterogeneity) either way. In this case, there are 15 models. These are all averaged in the BMA posterior distributions, but for simplicity, the top cluster models (posterior probability greater than 0.05) for both treatments are shown in Table 7. The homogeneity cluster model is the best model for the control treatment, but only has a 0.221 mass of probability, while the remaining probability is distributed among the other 14 models. The outcome is similar for the experimental treatment, for which the best option is now the heterogeneity cluster model. Obviously, these between-study variabilities should be considered in any meta-analysis.
Cosmi et al. used the Mantel-Haenszel method under random-effects to estimate the risks ratio, for which a value of 0.73 (95% CI 0.25-2.18) was obtained. The main conclusion derived from these data was that ticlopidine plus aspirin versus oral anticoagulation did not affect all-cause mortality. Using the proposed BMA approach, a risks ratio of 0.959 is obtained, which is in line with the value reported by Cosmi et al. However, once again, a realistic and larger interval becomes apparent when all of the uncertainty is managed via the BMA perspective (95% HDI 0.02-41.86).

Conclusions
The question of statistical heterogeneity in a meta-analysis based on a relatively small number of studies is important and has been addressed by Friede et al. [17], IntHout et al. [10] and Pateras et al. [31], among many others. Moreover, if zero events occur in these small trials, the challenge is multiplied. There has been much debate about the presence of significant between-study heterogeneity in these cases, and it has been reported that the misuse of continuity correction procedures can lead to incorrect inferences being drawn [9]. In this respect, Veroniki et al. [32] catalogued various methods for estimating between-study heterogeneity. The examples given in the present paper also show that the usual measures employed to determine the presence or otherwise of heterogeneity, I 2 and τ, may not consider all the scenarios contained in the meta-analysis, and therefore, it is currently recommended to accompany the I 2 statistic and the estimate of τ with a confidence interval. Such a confidence interval will be very wide, especially in the case of a small number of effect sizes in a meta-analysis.
The Bayesian model averaging that we propose is an intuitively attractive approach to the problem of accounting for heterogeneity uncertainty, even in meta-analyses with zero events. This method involves averaging all the heterogeneity structures (models) obtained by clustering the observed sample when making inferences about the true treatment effect. The number of models (types of heterogeneity) considered in BMA rises in line with the number k of studies, according to Bell's number [27]. However, this number is easily manageable for the case of few studies (2 for k = 2, 5 for k = 3, 15 for k = 4, and 52 for the limit case k = 5), and therefore, the computation of the posterior probabilities of each model and the derivation of the posterior BMA distribution of the parameter of interest are immediate and computationally feasible. In this respect, Rhodes et al. [33] conducted an extensive review of the Cochrane Database of Systematic Reviews (CDSR), developed in Davey et al. [34], and pointed out that: "Of 22453 meta-analyses from CDSR, containing at least two studies, just under 75% contained five or fewer studies". We recommend our method when the meta-analysis has 10 or fewer studies.
The clustering procedure considered takes into account all uncertainty in the heterogeneity structures. This fact has important consequences for treatment comparison, due to the likelihood of θ being different for each heterogeneity configuration.
On the other hand, observed heterogeneity could be a starting point to investigate sources of heterogeneity. In this paper, every potential cluster model is estimated. However, this might also be a starting point to investigate sources of heterogeneity before the simple model averaging takes place.
From these considerations, we conclude that the BMA procedure proposed is not only more compelling in a conceptual sense, but also provides novel probabilistic clustering methods that are capable of resolving challenging meta-analysis scenarios, in which there are few studies and sparse data. Funding: Financial support for this study was provided in part by Grant ECO2017-85577-P (Ministerio de Ciencia, Innovación y Universidades, Agencia Estatal de Investigación, Spain).