Bayesian Model-Averaged Meta-Analysis in Medicine

We outline a Bayesian model-averaged meta-analysis for standardized mean differences in order to quantify evidence for both treatment effectiveness $\delta$ and across-study heterogeneity $\tau$. We construct four competing models by orthogonally combining two present-absent assumptions, one for the treatment effect and one for across-study heterogeneity. To inform the choice of prior distributions for the model parameters, we used 50% of the Cochrane Database of Systematic Reviews to specify rival prior distributions for $\delta$ and $\tau$. The relative predictive performance of the competing models and rival prior distributions was assessed using the remaining 50\% of the Cochrane Database. On average, $\mathcal{H}_1^r$ -- the model that assumes the presence of a treatment effect as well as across-study heterogeneity -- outpredicted the other models, but not by a large margin. Within $\mathcal{H}_1^r$, predictive adequacy was relatively constant across the rival prior distributions. We propose specific empirical prior distributions, both for the field in general and for each of 46 specific medical subdisciplines. An example from oral health demonstrates how the proposed prior distributions can be used to conduct a Bayesian model-averaged meta-analysis in the open-source software R and JASP. The preregistered analysis plan is available at https://osf.io/zs3df/.


Introduction
Following Karl Pearson's first quantitative synthesis of clinical trials in 1904, meta-analysis gradually established itself as an irreplaceable method for statistics in medicine. 1 However, over a century later meta-analysis still presents formidable statistical challenges to medical practitioners, especially when the number of primary studies is low. In this case the estimation of across-study heterogeneity (i.e., across-study standard deviation) τ is problematic 2,3,4 ; moreover, these problematic τ estimates may subsequently distort the estimates of the overall treatment effect size δ. 2,5 The practical relevance of the small sample challenge is underscored by the fact that the median number of studies in a meta-analysis from the Cochrane Database of Systematic Reviews (CDSR) is only 3, with an interquartile range from 2 to 6. 6 One statistical method that has been proposed to address the small sample challenge is Bayesian estimation, either with weakly informative prior distributions, 7,8,9 predictive prior distributions based on pseudo-data, 10 or prior distributions informed by earlier studies. 11 These Bayesian techniques are well suited to estimate the model parameters when the data are scarce; however, by assigning continuous prior distributions to δ and τ , these estimation techniques implicitly assume that the treatment is effective and the studies are not homogeneous. 1 In order to validate these strong assumptions we may adopt the framework of Bayesian testing. Developed in the second half of the 1930s by Sir Harold Jeffreys,12,13 the Bayesian testing framework seeks to grade the evidence that the data provide for or against a specific value of interest such as δ = 0 and τ = 0 which corresponds to the null model of no effect and the fixed-effect model, respectively. Jeffreys argued that the testing question logically precedes the estimation question, and that more complex models (e.g., the models used for estimation, where δ and τ are free parameters) ought to be adopted only after the data provide positive evidence in their favor: "Until such evidence is actually produced the simpler hypothesis holds the field; the onus of proof is always on the advocate of the more complicated hypothesis." 14, p. 252 In the context of meta-analysis, Jeffreys's statistical philosophy demands that we acknowledge not only the uncertainty in the parameter values given a specific model, but also the uncertainty in the underlying models to which the parameters belong. Both types of uncertainty can be assessed and updated using a procedure known as Bayesian model-averaged (BMA) meta-analysis. 2 The BMA procedure applies different meta-analytic models to the data simultaneously, and draws inferences by taking into account all models, with their impact determined by their predictive performance for the observed data. 19,20,21 As in other applications of Bayesian statistics, BMA requires that all parameters are assigned prior distributions. However, in contrast to Bayesian estimation, Bayesian testing does not permit the specification of vague or "uninformative" prior distributions on the parameters of interest. Vague prior distributions assign most prior mass to implausibly large values, resulting in poor predictive performance. 13,21,22 In BMA, the relative impact of the models is determined by their predictive performance, and predictive performance in turn is determined partly by the prior distribution on the model parameters. In objective Bayesian statistics 23 so-called default prior distributions have been proposed; these default distributions meet a list of desiderata 24 and are intended for general use in testing. In contrast to this work, here we seek to construct and compare different prior distributions based on existing medical knowledge. 25,26,27 Specifically, we propose empirical prior distributions for δ and τ as applied to meta-analyses of continuous outcomes in medicine. To this aim we first used 50% of CDSR to develop candidate prior distributions and then used the remaining 50% of CDSR to evaluate their predictive accuracy and that of the associated models.
Below we first outline the BMA approach to meta-analyses and then present the results of a preregistered analysis procedure to obtain and assess empirical prior distributions for δ and τ for the medical field as a whole. Next we propose empirical prior distributions for the 46 specific medical subdisciplines defined by CDSR. Finally we demonstrate with a concrete example how our results can be applied in practice using the open-source statistical programs R 28 and JASP. 29

Bayesian Model-Averaged Meta-Analysis
The standard Bayesian random-effects meta-analysis assumes that a latent individual study effect θ i is drawn from a Gaussian group-level distribution with mean treatment effect δ and between-study heterogeneity τ . 30,31 Inference then concerns the posterior distributions for δ and τ . This estimation approach allows researchers to answer important questions such as "given that the treatment effect is nonzero, how large is it?" 3 and "given that there is between-study 1 Under a continuous prior distribution, the prior probability of any particular value (such as δ = 0, which represent the proposition that the treatment is ineffective) is zero. 2 A different approach is Bayesian model selection based on posterior (out-of-sample) predictive performance such as DIC/WAIC/LOO. 15,16 However, these approaches are unable to provide compelling support in favor of simple models. 17,18 3 More specific versions of this generic question are "given that the treatment effect is nonzero, is it positive or negative?" and "given that the treatment effect is nonzero, what is the posterior probability that it falls in the interval from a to b?" heterogeneity, how large is it?" Because the standard model assumes that the effect is nonzero, it cannot address the arguably more fundamental questions that involve a hypothesis test, 32,33 such as "how strong is the evidence in favor of the presence or absence of a treatment effect?" and "how strong is the evidence in favor of between-study heterogeneity (between-study standard deviation) versus homogeneity?" 34, p. 274 Here we outline a BMA approach that allows for both hypothesis testing and parameter estimation in a single statistical framework. 35 Our generic meta-analysis setup 19,20,36,37 (for the conceptual basis see Jeffreys) 13, p. 276-277 and p. 296 consists of the following four qualitatively different candidate hypotheses 4 : 1. the fixed-effect null hypothesis H f 0 : δ = 0 , τ = 0; 2. the fixed-effect alternative hypothesis H f 1 : δ ∼ g(·) , τ = 0; 3. the random-effects null hypothesis H r 0 : δ = 0, τ ∼ h(·); 4. the random-effects alternative hypothesis H r where δ represents the group-level mean treatment effect, τ represents the between-study standard deviation (i.e., the treatment heterogeneity), and g(·) and h(·) represent prior distributions that quantify the uncertainty about δ and τ , respectively. The four prior probabilities of the rival hypotheses are denoted by p(H f 0 ), p(H f 1 ), p(H r 0 ), and p(H r 1 ); these may or may not be set to 1 /4, reflecting a position of prior equipoise. The main advantage of this framework is that it does not fully commit to any single model on purely a priori grounds. Although in many situations the random-effects alternative hypothesis H r 1 is an attractive option, it may be less appropriate when the number of studies is small; in addition, as mentioned above, H r 1 assumes the effect to be present, whereas assessing the degree to which the data undercut or support this assumption may often be one of the primary inferential goals.
In our framework, after specifying the requisite prior distributions g(·) and h(·), the data drive an update from prior to posterior model probabilities, and pertinent conclusions are then drawn using BMA. 38,39 Specifically, the posterior odds of an effect being present, based on observed data y, is the ratio of the sum of posterior model probabilities for H f 1 and H r 1 over the sum of posterior model probabilities for H f 0 and H r 0 : Posterior odds for treatment effect = .
In model-averaging terms, this quantity is referred to as the posterior inclusion odds, as it refers to the post-data odds of 'including' the effect size parameter δ. As a measure of evidence, one may consider the change, brought about by the data, from prior inclusion odds to posterior inclusion odds. This change is known as the Bayes factor 22 , or one may quantify evidence by the change from prior to posterior inclusion odds:

BF rf
Inclusion Bayes factor for treatment heterogeneity Posterior inclusion odds for treatment heterogeneity Prior inclusion odds for treatment heterogeneity . ( An attractive feature of this framework is that it allows a graceful data-driven transition from an emphasis on fixedeffect models to random-effects models; with only few studies available, the fixed-effect models likely outpredict the random-effects models and therefore receive more weight. But as studies accumulate, and it becomes increasingly apparent that the treatment effect is indeed random, the influence of the random-effects models will wax and of the fixed-effect models will wane, until inference is dominated by the random-effects models. In addition, the Bayesian framework allows researchers to monitor the evidence as studies accumulate, without the need or want of corrections for optional stopping. 41 This is particularly relevant as the accumulation of studies is usually not under the control of a central agency, and the stopping rule is ill-defined. 42 Although theoretically promising, the practical challenge for our BMA meta-analysis is to determine appropriate prior distributions for δ and τ . Prior distributions that are too wide will waste prior mass on highly implausible parameter values, thus incurring a penalty for complexity that could have been circumvented by applying a more reasonably peaked prior distribution. On the other hand, prior distributions that are too narrow represent a highly risky bet; if the effect is not exactly where the peaked prior distribution guesses it to be, the model will incur a hefty penalty for predicting the data poorly, a penalty that could have been circumvented by reasonably widening the prior distribution.
There is no principled way around this dilemma: Bayes' rule dictates that evidence is quantified by predictive success, and predictions follow from the prior predictive distributions. 43 Thus, when the goal is to quantify evidence, the prior distributions warrant careful consideration. 44,45,46,47,48 Fortunately, the framework presented here contains only two key parameters, δ and τ ; moreover, a large clinical literature is available to help guide the specification of reasonable prior distributions. Our goal in this work is to use meta-analyses from the CDSR to create a series of informed prior distributions for both the effect size parameter δ and between-study variance parameter τ . 8,25,26,27 We will then assess the predictive adequacy of the various models in conjunction with the prior distributions on a hold-out validation set.

Candidate Prior Distributions
We developed and assessed prior distributions for the δ and τ parameters suitable for BMA of continuous outcomes using data from CDSR. 5 In the remainder of this work we adopt the terminology of Higgins et al. 49 : individual meta-analyses included in each Cochrane review are referred to as 'comparisons' and individual studies included in a comparison are referred to as 'studies'. All of the results were conducted using Cohen's d standardized mean differences (SMD). The analyses presented in this section were executed in accordance with a preregistration protocol (https://osf.io/zs3df/) unless explicitly mentioned otherwise.
In order to assess the predictive adequacy of the various prior distributions and models, we first randomly partitioned the data of the Cochrane reviews in a training and test set. The training set consisted of 3,092 comparisons with a total number of 23,333 individual studies, and the test set consisted of 3,091 comparisons with a total number of 22,117 individual studies. We used the training set to develop prior distributions for the δ and τ parameters and then assessed predictive accuracy using the test set.

Developing Prior Distributions Based on the Training Set
Left panel of Figure 1 outlines the data processing steps performed on the training set (further details are provided in the preregistration protocol, https://osf.io/zs3df/). First, in order to ensure that the training set yields estimates of τ that form a reliable basis for the construction of a prior distribution, we excluded comparisons with fewer than 10 studies. Second, we excluded comparisons for which at least one individual study was reported by the authors of the review to be "non-estimable" (i.e., the effect size of the original study could not be retrieved). Third, we transformed the reported raw mean differences to the SMD using the metafor R package. 50 Fourth, to ensure high consistency of the meta-analytic estimates, we re-estimated all comparisons using a frequentist random-effects meta-analytic model with restricted maximum likelihood estimator using the metafor R package. 50 These steps resulted in a final training set featuring 423 comparisons containing a total of 8,044 individual studies. The histograms and tick marks in Figure 2 display the δ and τ estimates from each comparison in the training set. 6 To develop candidate prior distributions for parameters δ and τ , we used the maximum likelihood estimator implemented in the fitdistrplus R package 51 to fit several distributions to the frequentist meta-analytic estimates from the training set. For the δ parameter, we considered normal and Student's t distributions fitted to the training set and compared them to an uninformed Cauchy distribution with scale 1 / √ 2 (a default choice in the field of psychology). 52 For the τ parameter, we considered half-normal, inverse-gamma, and gamma distributions fitted to the training set and compared 5 We identified systematic reviews in the CDSR through PubMed, limiting the period to Jan 2000 -May 2020. For that we used the NCBI's EUtils API with the following query: "Cochrane Database Syst Rev"[journal] AND ("2000/01/01"[PDAT]: "2020/05/31"[PDAT]). For each review, we downloaded the XML meta-analysis table file (rm5-format) associated with the review's latest version. We extracted the tables with continuous outcomes (i.e., MD and SMD) from these rm5-files with a custom PHP script. We labeled the tables with the Cochrane Review Group's taxonomy list for subfield analysis. 6 As specified in the preregistration protocol, we assumed that τ estimates lower than 0.01 are representative of H f · : τ = 0, and therefore these estimates were not used to determine candidate prior distributions for τ . Removing comparisons with at least one not estimable study (33 not estimable studies in 20 comparisons).
Transforming the reported raw mean differences to the standardized mean differences.
Re-estimating the comparisons using a random-effects meta-analytic model with restricted maximum likelihood.
Initial test data set: 3,091 comparisons (22,117  Removing comparisons with at least one not estimable study (156 not estimable studies in 85 comparisons).
Transforming the reported raw mean differences to the standardized mean differences.
Estimating the comparisons using Bayesian modelaveraged meta-analysis with different prior distribution configurations (10 comparisons were not estimable).  : Frequentist effect sizes estimates and candidate prior distributions from the training data set. Histogram and tick marks display the estimated effect size estimates (left) and between-study standard deviation estimates (right), whereas lines represent three associated candidate prior distributions for the population effect size parameter δ (left) and four candidate prior distributions for the population between-study standard deviation τ (right; see Table 1). Twelve effect sizes outside of the ±1.5 range are not shown and twenty-four τ estimates larger than 1 and sixty-eight τ estimates lower than 0.01 are not shown.
them to an uninformed uniform distribution on the range from 0 to 1. 45 The resulting distributions are summarized in Table 1 and their fit to the training set is visualized in Figure 2.

Assessing Prior Distributions Based on the Test Set
Right panel of Figure 1 outlines the data processing steps performed on the test set. Similarly to the training set, we removed non-estimable comparisons and transformed all effect sizes to SMD. However, in contrast to the training set, we retained all comparisons that feature at least 3 studies: there is no reason to limit the assessment of predictive performance to comparisons with at least 10 studies. These data processing steps resulted in a final test set consisting of 2,416 comparisons containing a total of 18,479 individual studies. The median number of studies in a comparison was 5 with an interquartile range from 3 to 9. Table 1: Candidate prior distributions for the δ and τ parameters as obtained from the training set. The inverse-gamma and gamma distributions follow the shape and scale parametrization and the Studen-t distributions follow the location, scale, and degrees of freedom parametrization. See Figure 2.
For each possible pair of candidate prior distributions depicted in Table 1 In the first analysis we evaluate the predictive performance associated with the different prior distribution configurations as implemented in H r 1 , that is the random-effects model that allows both δ and τ to be estimated from the data. Specifically, under H r 1 there are 3 × 4 = 12 prior configurations and each is viewed as a model yielding predictions. The prior probability of each prior configuration is 1 /12 ≈ 0.083 and the predictive accuracy of each prior configuration is assessed with the 2406 comparisons from the test set. Table 2 lists the 12 different prior configurations and summarizes the number of times their posterior probability was ranked 1, 2, . . . , 12. The results show that informed configurations generally outperformed the uninformed configurations (i.e., the Cauchy(0, 1 / √ 2) distribution on δ and the U(0, 1) distribution on τ ). The worst ranking performance was obtained with prior configuration 1 (i.e., uniformed distributions for both δ and τ ). Prior configurations 2, 3, 4, 5, and 9 feature an uninformed prior distribution on either δ or τ , and also did not perform well in terms of posterior rankings. The same holds for prior distribution configurations with the half-normal prior distribution for the τ parameter (i.e., prior configurations 2, 6, and 10). The best performing prior distribution configurations (i.e., numbers 7, 11, and 12) used more data-driven prior distributions for both δ (i.e., fitted normal or t-distributions) and τ (i.e., fitted inverse-gamma or gamma).  Figure 3 shows that, on average, the different prior distribution configurations perform similarly. As suggested by the posterior rankings from Table 2, configuration 1 predicted the data relatively poorly, resulting in an average posterior probability of 0.06; in contrast, configuration 11 predicted the data relatively well, resulting in an average posterior probability of 0.11. However, these posterior probabilities differ only modestly from the prior probability of 1 /12 ≈ 0.083, and the Bayes factor associated with the comparison of posterior probabilities of 0.11 and 0.06 is less than 2.
In sum, among the 12 prior configurations under H r 1 the best predictive performance was consistently obtained by data-driven priors, and the worst predictive performance was obtained by uninformed priors (cf. the posterior rankings from Table 2). However, the extent of this predictive advantage is relatively modest: starting from a prior probability of 1 /12 ≈ 0.083, the worst prior configuration has an average posterior probability of 0.06, and the best prior configuration has average posterior probability of 0.11 (cf. Figure 3). 8

Posterior Probability of the Four Model Types
In the second analysis we evaluate the predictive performance of the four meta-analytic model types (i.e., H f 0 , H f 1 , H r 0 , and H r 1 ) by model-averaging across all prior distribution configurations, separately for each of the 2406 comparisons. Table 3 shows the prior model probabilities obtained by first assigning probability 1 /4 to each of the four model types, and then spreading that probability out evenly across the constituent prior distribution configurations. 33, p. 47 For any particular comparison, a model type's model-averaged posterior probability is obtained by summing the posterior probabilities of the constituent prior distribution configurations. For example, the model-averaged posterior probability for H r 1 is obtained by summing the posterior probabilities for the 12 possible prior configurations, each of them associated with prior probability 1 /48 (cf. Table 3). Table 2: Ranking totals for each prior configuration based on the 2406 comparisons in the test set. The numbers indicate how many times a specific prior configuration attained a specific posterior probability rank amongst the 12 possible prior configurations. Rank '1' represents the best performance. Note that these rankings are conditional on assuming the meta-analytic model H r 1 (i.e., the posterior probabilities of the other meta-analytic models are not considered).

Rank
Prior δ  Table 2. Table 4 lists the four model types and summarizes the number of times their model-averaged posterior probability was ranked 1, 2, . . . , 4. The results show that, across all comparisons, complex models generally received more support than simple models. The model that predicted the data best was H r 1 , the random-effects model that assumes the presence of an effect; the model that predicted the data worst was H f 0 , the fixed-effect model that assumes the absence of an effect. However, even H f 0 outpredicted the other three model types in 662/2406 ≈ 28% of comparisons. Table 4 also shows the model-averaged posterior model probability across all comparisons. In line with the ranking results, the average probability for H f 0 decreased from 0.25 to 0.19, whereas that for H r 1 increased from 0.25 to 0.36. Nevertheless, the support for H r 1 across all comparisons is not overwhelming and does not appear to provide an empirical license to ignore H f 0 (or any of the other three model types) from the outset. Left panel of Figure 4 displays the model-averaged posterior probability for each model type across the 2406 comparisons. It is apparent that the posterior probability is highest for H r 1 . However, for a substantial number of comparisons (i.e., (662 + 573 + 406)/2406 ≈ 68%, cf. Table 4) a different model type performs better. For instance -and in contrast to popular belief-the fixed-effect models H f 1 and H f 0 together show the best predictive performance in a slight majority of comparisons (i.e., (662 + 573)/2406 ≈ 51%).

Inclusion Bayes Factors
In the third analysis we assess the inclusion Bayes factors for a treatment effect (cf. Eq. 1) and for heterogeneity (cf. Eq. 2); that is, we model-average across all prior distribution configurations and across two model types, separately for each of the 2406 comparisons. First, the inclusion BF 10 for a treatment effect quantifies the evidence that the data provide for the presence vs. the absence of a group-level effect, taking into account the model uncertainty associated with whether the effect is fixed or random. The left panel of Figure 5 displays a histogram of the log of the model-averaged BF 10 for the test set featuring 2406 comparisons. The histogram is noticeably right-skewed, which affirms the regularity that it is easier to obtain compelling evidence for the presence rather than the absence of an effect. 13, p. 196-197;9 Evidence for the presence of an effect was obtained in a small majority of the comparisons (i.e., 1336/2406 ≈ 55.5%).
Second, the inclusion BF rf for heterogeneity quantifies the evidence that the data provide for the presence vs. absence of between-study variability, taking into account the model uncertainty associated with whether the group-level effect is present or absent. The right panel of Figure 5 displays a histogram of the log of the model-averaged BF rf for the test set  .36 *Prior model probability **Average posterior model probability featuring 2406 comparisons. The right-skew again confirms the regularity: it is easier to find compelling evidence for heterogeneity than for homogeneity. Nevertheless, the data provide evidence for heterogeneity only in a slight majority of 1227/2406 ≈ 51.0% of the comparisons.
In sum, the inclusion Bayes factors revealed that in nearly half of the comparisons from the test set, the data provide evidence in favor of the absence of a treatment effect (i.e., 44.5%) and provide evidence in favor of the absence of heterogeneity (i.e., 49.0%). The distribution of the log Bayes factors is asymmetric, indicating that it is easier to obtain compelling evidence for the presence of a treatment effect (rather than for its absence) and for the presence of heterogeneity (rather than for homogeneity).

Exploratory Analysis: Model-Averaging Across Prior Distributions Under H r 1
To further investigate the predictive performance of the prior distributions, we performed one additional analysis that was not preregistered in the original analysis plan. We focused on H r 1 and evaluated the prior distributions for each parameter by model-averaging across the possible prior distributions for the other parameter. For instance, to obtain the model-averaged posterior probability for the Cauchy prior distribution on the δ parameter (i.e., δ ∼ Cauchy(0, 1 /  Table 4. In the right panel, the prior probability is 1 /3 for each prior distribution on δ, and 1 /4 for each prior distribution on τ , see also Table 5.  Table 3. This way we obtain an assessment of the relative predictive performance of a particular prior distribution, averaging over the uncertainty on the prior distribution for the other parameter. Table 5 lists the prior distributions and gives the number of times their model-averaged posterior probability attained a particular ranking. Consistent with the results reported earlier, the more data-driven prior distributions generally received more support than the prior distributions that are less informed. For the δ parameter, the best performing prior distribution was δ ∼ T (0, 0.33, 3); for the τ parameter, the best performing prior distribution was τ ∼ Inv-Gamma(1.26, 0.24). Although the preference for the data-driven prior distributions is relatively consistent, it is not particularly pronounced, echoing the earlier results. Specifically, Table 5 also shows the model-averaged posterior model probability across all comparisons. For the δ parameter, the t-prior has a model-averaged posterior probability of 0.39 (up from 1 /3), but the Cauchy prior retains a non-negligible probability of 0.25. For the τ parameter, the different prior distributions perform even more similarly; on average, the worst prior distribution is τ ∼ U (0, 1), and yet its model-averaged posterior model probability equals 0.23, down from 1 /4 but only a little. Likewise, the on-average best prior distribution is Table 5: Ranking totals for each prior distribution in H r 1 based on the 2406 comparisons in the test set. The numbers indicate how many times a specific prior distribution attained a specific posterior probability rank. Rank '1' represents the best performance. The rankings reflect predictive adequacy that is model-averaged across the possible prior distribution configurations of the other parameter.  Figure 4 displays the model-averaged posterior probability for each prior distribution across the 2406 comparisons. The figure confirms that the data-driven prior distributions perform somewhat better than the relatively uninformed prior distributions. The color band is darker red, on average, for the prior distributions with the highest posterior model probabilities, that is, δ ∼ T (0, 0.33, 3) and τ ∼ Inv-Gamma(1.26, 0.24).

Exploratory Analysis: Subfield-Specific Prior Distributions
Medical subfields may differ both in the typical size of the effects and in their degree of heterogeneity. In recognition of this fact we sought to develop empirical prior distributions for δ and τ that are subfield-specific. We differentiated between 47 medical subfields according to the taxonomy of the Cochrane Review Group. Based on their relatively good predictive performance detailed in the previous sections, we selected a t-distribution for the δ parameter (i.e., for subfield i, δ i ∼ T (µ = 0, σ i , ν i )) and an inverse-gamma distribution for the τ parameter (i.e., for subfield i, τ i ∼ Inverse-gamma(α i , β i )).
To estimate the parameters of these distributions separately for each subfield, we used the complete data set and proceeded analogously to the training set preparation: we removed comparisons with non-estimable studies, only used comparisons with at least ten studies, re-estimated the comparisons with a restricted maximum likelihood estimator, and removed comparisons with τ < 0.01 estimates. These frequentist estimates were used as input for constructing the data-driven subfield-specific prior distributions. However, since many subfields contain only a limited number of comparisons, we used Bayesian hierarchical estimation with weakly informative priors on the hyperparameters. The hierarchical aspect of the estimation procedure shrinks the estimated parameter values towards the grand mean, a tendency that is more pronounced if the estimated field-specific value is both extreme and based on relatively little information. 59,60,61 Specifically, we assumed that all field-specific parameters (i.e., σ i , ν i , α i , and β i ) are governed by positive-only normal distributions. For the t-distribution, we assigned positive-only Cauchy(0, k) prior distributions both to the across-field normal mean and to the across-field normal standard deviation, with k = 1 for parameter σ and k = 10 for parameter ν. For the inverse-gamma distribution, we assigned positive-only Cauchy(0, 1) prior distributions both to the across-field normal mean and to the across-field normal standard deviation for shape parameter α and scale parameter β. The hierarchical models were estimated using the rstan R package 62 that interfaces with the Stan probabilistic modeling language. 63 The Stan code is available alongside the supplementary materials at https://osf.io/zs3df/. Table 6 lists the 46 different subfields (the 47 th subfield "Multiple Sclerosis and Rare Diseases of the CNS" featured two comparisons, both of which were excluded based on the τ < 0.01 criterion), the associated number of comparisons and studies, and the estimated distributions for both δ and τ . The scale estimates for the δ parameter show considerable variation, ranging from 0.18 ("Developmental, Psychosocial and Learning Problems") to 0.60 ("Hepato-Biliary"). 10 A similar variation is present in the estimated distributions for the τ parameter. Figure 6 visualizes the prior distributions for each subfield.  Table 6.

Example: Dentine Hypersensitivity
We demonstrate BMA meta-analysis with an example from oral health. Poulsen et al. 64 considered the effect of potassium-containing toothpaste on dentine hypersensitivity. Five studies with a tactile outcome assessment were subjected to a meta-analysis. In their review, Poulsen et al. 64  We re-analyze the Poulsen et al. 64 comparison using the BMA meta-analysis implementation in the open-source statistical software package JASP (jasp-stats.org). 29,65,66,67,68 Appendix A provides the same analysis in R 28 using the metaBMA package. 53 Figure 7 shows the JASP graphical user interface with the left panel specifying the analysis setting and the right panel displaying the default output. After loading the data into JASP, the BMA meta-analysis can be performed by activating the "Meta-Analysis" module after clicking the blue "+" button in the top right corner, choosing "Meta-Analysis" from the ribbon at the top, and then selecting "Bayesian Meta-Analysis" from the drop-down menu. In the left input panel, we move the study effect sizes and standard errors into the appropriate boxes and adjust the prior distributions under the "Prior" tab to match the subfield-specific distributions given in Table 6. Specifically, for the "Oral Health" subfield the prior distributions are δ ∼ T (0, 0.51, 5) and τ ∼ Inv-Gamma(1.79, 0.28).
The JASP output panel displays the corresponding BMA meta-analysis results. The "Posterior Estimates per Model" table summarizes the estimates and evidence from the fixed-effect models, random-effects models, and finally the model-averaged results. The final row of the table shows an effect size estimate δ = 1.082, 95 % CI [0.686, 1.412] which is slightly lower and more conservative than the one provided by the frequentist random-effects meta-analysis, further quantified with extreme evidence for the presence of an effect, BF 10 = 218.53, and moderate evidence for the presence of heterogeneity, BF rf = 3.52. The JASP output panel also presents a forest plot that visualizes the observed effects size estimates from the individual studies, the overall fixed-effect and random-effects meta-analytic estimates and the corresponding model-averaged effect size estimate.
The JASP interface provides additional options not discussed here, such as (a) visualizing the prior and posterior distributions; (b) visualizing the estimated effect sizes from individual studies; (c) performing one-sided hypothesis tests; (d) updating evidence sequentially, study-by-study; (e) adding ordinal constraints 69 ; and (f) adjustments for publication bias. 70,71 Table 6: Subfield-specific prior distributions for 46 individual topics from the Cochrane database of systematic reviews estimated by hierarchical regression based on the complete data set. The t-distribution follows a location, scale, and degrees of freedom parametrization and the inverse-gamma distribution follows a shape and scale parametrization. See also Figure 6.

Concluding Comments
In this article, we introduced Bayesian model-averaged meta-analysis for continuous outcomes in medicine. The proposed methodology provides a principled way to integrate, quantify, and update uncertainty regarding both parameters and models. Specifically, the methodology allows researchers to simultaneously test for and estimate effect size and heterogeneity without committing to a particular model in an all-or-none fashion. In Bayesian model-averaged metaanalysis, multiple models are considered simultaneously, and inference is proportioned to the support that each model receives from the data. This eliminates the need for stage-wise, multi-step inference procedures that first identify a single preferred model (e.g., a fixed-effect model or a random-effects model) and then interpret the model parameters without acknowledging the uncertainty inherent in the model selection stage. The multi-model approach advocated here also decreases the potential impact of model misspecification.
Bayesian model-averaged meta-analysis comes with the usual advantages of Bayesian statistics -the ability to quantify evidence in favor or against any hypothesis (including the null hypothesis), the ability to discriminate between absence of evidence and evidence of absence, 72,73 the ability to monitor evidence as individual studies accumulate, 42 the straightforward interpretation of the results (i.e., probability statements that refer directly to parameters and hypotheses), 74 and the opportunity to incorporate historical information. 75,76 In this article, our goal was to take advantage of the existing medical knowledge base in order to propose and assess prior distributions that allow for more efficient inference.
Following a preregistered analysis plan, we fitted and assessed different prior distributions for both effect size δ and heterogeneity τ using comparisons of continuous outcomes from the Cochrane database of systematic reviews. We fitted prior distributions based on a training set of randomly selected comparisons, and then evaluated predictive performance based on a test set. The results showed that predictive performance on the test set was relatively similar for the different data-driven prior distributions. Moreover, and in contrast to popular belief and recommendations, 8,9,77,78 we did not find that the random-effects meta-analytic model provided a superior account of the data: the random-effects meta-analytic models outpredicted their fixed-effect counterparts in only 51.0% of comparisons. Although the random-effects alternative hypothesis H r 1 showed the best predictive performance on average, the data increased its model-averaged posterior probability from 0.25 to only 0.36, leaving 0.64 for the three competing model types (i.e., a model with no heterogeneity, a model without an effect, and a model without both).
Based on the outcome of our preregistered analysis, we used the data from Cochrane database of systematic reviews to develop empirical prior distributions for continuous outcomes in 46 different medical subfields. Finally, we applied Bayesian model-averaged meta-analysis with subfield-specific prior distributions to an example from oral health, using the free statistical software packages R and JASP. We believe that the proposed Bayesian methodology provides an alternative perspective on meta-analysis that is informed, efficient, and insightful.
fit_example <-meta_bma(y = study_effect_size, SE = study_se, data = data, d = prior("t", c(location = 0, scale = 0.51, nu = 5)), tau = prior("invgamma", c(shape = 1.79, scale = 0.28)), control = list(adapt_delta = .90)) To obtain the numerical summaries of the estimated model, we just execute the name of the object containing the fitted model: fit_example The resulting output corresponds to that given by JASP output (up to MCMC error): The inclusion Bayes factor quantifying the evidence in favor of heterogeneity can be obtained using Equation 2 and output from the "Model posterior probabilities" table.