`vignettes/MedicineBMA.Rmd`

`MedicineBMA.Rmd`

Bayesian model-averaged meta-analysis allows researchers to seamlessly incorporate available prior information into the analysis (Bartoš, Gronau, et al., 2021; Gronau et al., 2017, 2021). This vignette illustrates how to do this with an example from Bartoš, Gronau, et al. (2021), who developed informed prior distributions for meta-analyses of continuous outcomes based on the Cochrane database of systematic reviews. Then, we extend the example by incorporating publication bias adjustment with robust Bayesian meta-analysis (Bartoš, Maier, et al., 2021; Maier et al., n.d.).

We illustrate how to fit the informed BMA (not adjusting for publication bias) using the `RoBMA`

R package. For this purpose, we reproduce the dentine hypersensitivity example from Bartoš, Gronau, et al. (2021), who reanalyzed five studies with a tactile outcome assessment that were subjected to a meta-analysis by Poulsen et al. (2006).

We load the dentine hypersensitivity data included in the package.

```
library(RoBMA)
data("Poulsen2006", package = "RoBMA")
Poulsen2006
#> d se study
#> 1 0.9073050 0.2720456 STD-Schiff-1994
#> 2 0.7207589 0.1977769 STD-Silverman-1996
#> 3 1.3305829 0.2721648 STD-Sowinski-2000
#> 4 1.7688872 0.2656483 STD-Schiff-2000-(2)
#> 5 1.3286828 0.3582617 STD-Schiff-1998
```

To reproduce the analysis from the example, we need to set informed empirical prior distributions for the effect sizes (\(\mu\)) and heterogeneity (\(\tau\)) parameters that Bartoš, Gronau, et al. (2021) obtained from the Cochrane database of systematic reviews. We can either set them manually,

```
fit_BMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
priors_effect = prior(distribution = "t", parameters = list(location = 0, scale = 0.51, df = 5)),
priors_heterogeneity = prior(distribution = "invgamma", parameters = list(shape = 1.79, scale = 0.28)),
priors_bias = NULL,
transformation = "cohens_d", seed = 1, parallel = TRUE)
```

with `priors_effect`

and `priors_heterogeneity`

corresponding to the \(\delta \sim T(0,0.51,5)\) and \(\tau \sim InvGamma(1.79,0.28)\) informed prior distributions for the “oral health” subfield and removing the publication bias adjustment models by setting `priors_bias = NULL`

\(^1\).

Alternatively, we can utilize the `informed_prior`

function that prepares informed prior distributions for the individual medical subfields automatically.

```
fit_BMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
priors_bias = NULL,
transformation = "cohens_d", seed = 1, parallel = TRUE)
```

The `name`

argument specifies the medical subfield name (use `print(BayesTools::prior_informed_medicine_names)`

to check names of all available subfields). The `parameter`

argument specifies whether we want prior distribution for the effect size or heterogeneity. Finally, the `type`

argument specifies what type of measure we use for the meta-analysis (only standardized mean differences `smd`

are available at the time of writing the vignette; see `?prior_informed`

for more details regarding the informed prior distributions).

We obtain the output with the `summary`

function. Adding the `conditional = TRUE`

argument allows us to inspect the conditional estimates, i.e., the effect size estimate assuming that the models specifying the presence of the effect are true and the heterogeneity estimates assuming that the models specifying the presence of heterogeneity are true\(^2\).

```
summary(fit_BMA, conditional = TRUE)
#> Call:
#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
#> transformation = "cohens_d", priors_effect = prior_informed(name = "oral health",
#> parameter = "effect", type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
#> parameter = "heterogeneity", type = "smd"), priors_bias = NULL,
#> parallel = TRUE, seed = 1)
#>
#> Robust Bayesian meta-analysis
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 2/4 0.500 0.995 217.517
#> Heterogeneity 2/4 0.500 0.778 3.511
#> Bias 0/4 0.000 0.000 0.000
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> mu 1.076 1.088 0.664 1.422
#> tau 0.231 0.208 0.000 0.726
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
#>
#> Conditional estimates:
#> Mean Median 0.025 0.975
#> mu 1.082 1.090 0.701 1.422
#> tau 0.297 0.255 0.075 0.779
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
```

The output from the `summary.RoBMA()`

function has 3 parts. The first one under the ‘Robust Bayesian Meta-Analysis’ heading provides a basic summary of the fitted models by component types (presence of the Effect/Heterogeneity/Publication bias). We can see that there are no models correcting for publication bias (we disabled them by setting `priors_bias = NULL`

). Furthermore, the table summarizes the prior and posterior probabilities and the inclusion Bayes factors of the individual components. The results show that the inclusion Bayes factor for the effect corresponds to the one reported in Bartoš, Gronau, et al. (2021) package, \(\text{BF}_{10} = 218.53\) and \(\text{BF}_{\text{rf}} = 3.52\) (up to an MCMC error).

The second part under the ‘Model-averaged estimates’ heading displays the parameter estimates model-averaged across all specified models (i.e., including models specifying the effect size to be zero). We ignore this section and move to the last part.

The third part under the ‘Conditional estimates’ heading displays the conditional effect size estimate corresponding to the one reported in Bartoš, Gronau, et al. (2021), \(\delta = 1.082\), 95% CI [0.686,1.412], and a heterogeneity estimate (not reported previously).

The `RoBMA`

package provides extensive options for visualizing the results. Here, we visualize the prior (grey) and posterior (black) distribution for the mean parameter.

`plot(fit_BMA, parameter = "mu", prior = TRUE)`

By default, the function plots the model-averaged estimates across all models; the arrows represent the probability of a spike, and the lines represent the posterior density under models assuming non-zero effect. The secondary y-axis (right) shows the probability of the spike (at value 0) decreasing from .50, to 0.005 (also obtainable from the ‘Robust Bayesian Meta-Analysis’ field in `summary.RoBMA()`

function).

To visualize the conditional effect size estimate, we can add the `conditional = TRUE`

argument,

`plot(fit_BMA, parameter = "mu", prior = TRUE, conditional = TRUE)`

which displays only the model-averaged posterior distribution of the effect size parameter for models assuming the presence of the effect.

We can also visualize the estimates from the individual models used in the ensemble. We do that with the `plot_models()`

function, which visualizes the effect size estimates and 95% CI of each specified model included in the ensemble. Model 1 corresponds to the fixed effect model assuming the absence of the effect, \(H_0^{\text{f}}\), Model 2 corresponds to the random effect model assuming the absence of the effect, \(H_0^{\text{r}}\), Model 3 corresponds to the fixed effect model assuming the presence of the effect, \(H_3^{\text{f}}\), and Model 4 corresponds to the random effect model assuming the presence of the effect, \(H_1^{\text{r}}\),). The size of the square representing the mean estimate reflects the posterior model probability of the model, which is also displayed in the right-hand side panel. The bottom part of the figure shows the model averaged-estimate that is a combination of the individual model posterior distributions weighted by the posterior model probabilities.

`plot_models(fit_BMA)`

We see that the posterior model probability of the first two models decreased to essentially zero (when rounding to two decimals), completely omitting their estimates from the figure. Furthermore, the much larger box of Model 4 (the random effect model assuming the presence of the effect) shows that Model 4 received the largest share of the posterior probability, \(P(H_1^{\text{r}}) = 0.77\))

The last type of visualization that we show here is the forest plot. It displays the original studies’ effects and the meta-analytic estimate within one figure. It can be requested by using the `forest()`

function. Here, we again set the `conditional = TRUE`

argument to display the conditional model-averaged effect size estimate at the bottom.

`forest(fit_BMA, conditional = TRUE)`

For more options provided by the plotting function, see its documentation using `?plot.RoBMA()`

, `?plot_models()`

, and `?forest()`

.

Finally, we illustrate how to adjust our informed BMA for publication bias with robust Bayesian meta-analysis (Bartoš, Maier, et al., 2021, p. @maier2020robust). In short, we specify additional models assuming the presence of the publication bias and correcting for it by either specifying a selection model operating on \(p\)-values (Vevea & Hedges, 1995) or by specifying a publication bias adjustment method correcting for the relationship between effect sizes and standard errors – PET-PEESE (Stanley, 2017; Stanley & Doucouliagos, 2014). See Bartoš et al. (n.d.) for a tutorial.

To obtain a proper before and after publication bias adjustment comparison, we fit the informed BMA model but using the default effect size transformation (Fisher’s \(z\)).

```
fit_BMAb <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
priors_bias = NULL,
seed = 1, parallel = TRUE)
```

```
summary(fit_BMAb, conditional = TRUE)
#> Call:
#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
#> priors_effect = prior_informed(name = "oral health", parameter = "effect",
#> type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
#> parameter = "heterogeneity", type = "smd"), priors_bias = NULL,
#> parallel = TRUE, seed = 1)
#>
#> Robust Bayesian meta-analysis
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 2/4 0.500 0.997 347.932
#> Heterogeneity 2/4 0.500 0.723 2.608
#> Bias 0/4 0.000 0.000 0.000
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> mu 1.045 1.052 0.705 1.344
#> tau 0.186 0.163 0.000 0.623
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
#>
#> Conditional estimates:
#> Mean Median 0.025 0.975
#> mu 1.048 1.053 0.720 1.344
#> tau 0.256 0.220 0.064 0.681
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
```

We obtain noticeably stronger evidence for the presence of the effect. This is a result of placing more weights on the fixed-effect models, especially the fixed-effect model assuming the presence of the effect \(H_1^f\). The increase in the posterior model probability of \(H_1^f\) is, in our case, the fact that it predicted the data slightly better after removing the correlation between effect sizes and standard errors (a result of using Cohen’s \(d\) for model fitting). Nevertheless, the conditional effect size estimate stayed almost the same.

Now, we fit the publication bias-adjusted model by simply removing the `priors_bias = NULL`

, which allows us to obtain the default 36 models ensemble called RoBMA-PSMA (Bartoš, Maier, et al., 2021).

```
fit_RoBMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
seed = 1, parallel = TRUE)
```

```
summary(fit_RoBMA, conditional = TRUE)
#> Call:
#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
#> priors_effect = prior_informed(name = "oral health", parameter = "effect",
#> type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
#> parameter = "heterogeneity", type = "smd"), parallel = TRUE,
#> seed = 1)
#>
#> Robust Bayesian meta-analysis
#> Components summary:
#> Models Prior prob. Post. prob. Inclusion BF
#> Effect 18/36 0.500 0.858 6.022
#> Heterogeneity 18/36 0.500 0.714 2.502
#> Bias 32/36 0.500 0.697 2.304
#>
#> Model-averaged estimates:
#> Mean Median 0.025 0.975
#> mu 0.722 0.880 0.000 1.283
#> tau 0.202 0.161 0.000 0.799
#> omega[0,0.025] 1.000 1.000 1.000 1.000
#> omega[0.025,0.05] 0.943 1.000 0.329 1.000
#> omega[0.05,0.5] 0.874 1.000 0.071 1.000
#> omega[0.5,0.95] 0.855 1.000 0.042 1.000
#> omega[0.95,0.975] 0.866 1.000 0.050 1.000
#> omega[0.975,1] 0.897 1.000 0.057 1.000
#> PET 0.931 0.000 0.000 4.927
#> PEESE 1.131 0.000 0.000 12.261
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
#> (Estimated publication weights omega correspond to one-sided p-values.)
#>
#> Conditional estimates:
#> Mean Median 0.025 0.975
#> mu 0.838 0.938 -0.035 1.297
#> tau 0.285 0.227 0.064 0.906
#> omega[0,0.025] 1.000 1.000 1.000 1.000
#> omega[0.025,0.05] 0.736 0.829 0.092 1.000
#> omega[0.05,0.5] 0.411 0.373 0.014 0.951
#> omega[0.5,0.95] 0.320 0.249 0.008 0.919
#> omega[0.95,0.975] 0.376 0.311 0.009 0.958
#> omega[0.975,1] 0.518 0.425 0.010 1.000
#> PET 2.909 3.136 0.171 5.614
#> PEESE 7.048 6.034 0.375 18.162
#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
#> (Estimated publication weights omega correspond to one-sided p-values.)
```

We notice that additional values in the ‘Components summary’ table in the ‘Bias’ row. The model is now extended with 32 publication bias adjustment models that account for 50% of the prior model probability. When comparing the RoBMA to the second BMA fit, we notice a large decrease in the inclusion Bayes factor for the presence of the effect \(\text{BF}_{10} = 6.02\) vs. \(\text{BF}_{10} = 347.93\), which still, however, presents moderate evidence for the presence of the effect. We can further quantify the evidence in favor of the publication bias with the inclusion Bayes factor for publication bias \(\text{BF}_{pb} = 2.30\), which can be interpreted as weak evidence in favor of the publication bias.

We can also compare the publication bias unadjusted and publication bias-adjusted conditional effect size estimates. Including models assuming publication bias into our model-averaged estimate (assuming the presence of the effect) slightly decreases the estimated effect to \(\delta = 0.838\), 95% CI [-0.035, 1.297] with a much wider confidence interval, as visualized in the prior and posterior conditional effect size estimate plot.

`plot(fit_RoBMA, parameter = "mu", prior = TRUE, conditional = TRUE)`

\(^1\) The additional setting `transformation = "cohens_d"`

allows us to get more comparable results with the `metaBMA`

R package since RoBMA otherwise internally transforms the effect sizes to Fisher’s \(z\) for the fitting purposes and the `seed = 1`

and `parallel = TRUE`

options grant us exact reproducibility of the results and parallelization of the fitting process.

\(^2\) The Model-averaged estimates that `RoBMA`

returns by default model-averaged across all specified models – a different behavior from the `metaBMA`

package that by default returns what we call “conditional” estimates in `RoBMA`

.

Bartoš, F., Gronau, Q. F., Timmers, B., Otte, W. M., Ly, A., & Wagenmakers, E.-J. (2021). Bayesian model-averaged meta-analysis in medicine. *Statistics in Medicine*. https://doi.org/10.1002/sim.9170

Bartoš, F., Maier, M., Quintana, D. S., & Wagenmakers, E.-J. (n.d.). Adjusting for publication bias in JASP & R – selection models, PET-PEESE, and robust Bayesian meta-analysis. *Advances in Methods and Practices in Psychological Science*. https://doi.org/10.31234/osf.io/75bqn

Bartoš, F., Maier, M., Wagenmakers, E.-J., Doucouliagos, H., & Stanley, T. D. (2021). *Robust Bayesian meta-analysis: Model-averaging across complementary publication bias adjustment methods*. PsyArXiv. https://doi.org/10.31234/osf.io/kvsp7

Gronau, Q. F., Heck, D. W., Berkhout, S. W., Haaf, J. M., & Wagenmakers, E.-J. (2021). A primer on Bayesian model-averaged meta-analysis. *Advances in Methods and Practices in Psychological Science*, *4*(3), 25152459211031256. https://doi.org/10.1177/25152459211031256

Gronau, Q. F., Van Erp, S., Heck, D. W., Cesario, J., Jonas, K. J., & Wagenmakers, E.-J. (2017). A Bayesian model-averaged meta-analysis of the power pose effect with informed and default priors: The case of felt power. *Comprehensive Results in Social Psychology*, *2*(1), 123–138. https://doi.org/10.1080/23743603.2017.1326760

Maier, M., Bartoš, F., & Wagenmakers, E.-J. (n.d.). Robust Bayesian meta-analysis: Addressing publication bias with model-averaging. *Psychological Methods*. https://doi.org/10.31234/osf.io/u4cns

Poulsen, S., Errboe, M., Mevil, Y. L., & Glenny, A.-M. (2006). Potassium containing toothpastes for dentine hypersensitivity. *Cochrane Database of Systematic Reviews*, *3*. https://doi.org/10.1002/14651858.cd001476.pub2

Stanley, T. D. (2017). Limitations of PET-PEESE and other meta-analysis methods. *Social Psychological and Personality Science*, *8*(5), 581–591. https://doi.org/10.1177/1948550617693062

Stanley, T. D., & Doucouliagos, H. (2014). Meta-regression approximations to reduce publication selection bias. *Research Synthesis Methods*, *5*(1), 60–78. https://doi.org/10.1002/jrsm.1095

Vevea, J. L., & Hedges, L. V. (1995). A general linear model for estimating effect size in the presence of publication bias. *Psychometrika*, *60*(3), 419–435. https://doi.org/10.1007/BF02294384