# Bayesian Meta-Analysis

#### František Bartoš

#### 28th of April 2026

This vignette is a starting point for using the `RoBMA` R package (Bartoš & Maier, 2020). We walk through the BCG vaccine dataset from the `metafor` package (Viechtbauer, 2010), showing each `brma()` call alongside its `metafor::rma()` counterpart and building from a simple random-effects model up to meta-regression with a continuous and a categorical moderator. The interface is intentionally close to the `metafor` package: the same effect-size columns, moderator formulas, and standard plotting and diagnostic functions apply.

All `brma()` calls keep the default prior distributions; see the Prior Distributions vignette for the prior definitions and customization options.

## Random-Effects Model

The BCG vaccine dataset from the `metadat` package contains 13 randomized trials on tuberculosis prevention. `escalc()` computes log risk ratios `yi` and their sampling variances `vi`.

    data("dat.bcg", package = "metadat")

    dat <- metafor::escalc(
      ai = tpos, bi = tneg, ci = cpos, di = cneg,
      measure = "RR", data = dat.bcg
    )

`metafor::rma()` fits the random-effects model.

    fit1_metafor <- metafor::rma(yi, vi, data = dat, method = "REML")
    fit1_metafor
    #>
    #> Random-Effects Model (k = 13; tau^2 estimator: REML)
    #>
    #> tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
    #> tau (square root of estimated tau^2 value):      0.5597
    #> I^2 (total heterogeneity / total variability):   92.22%
    #> H^2 (total variability / sampling variability):  12.86
    #>
    #> Test for Heterogeneity:
    #> Q(df = 12) = 152.2330, p-val < .0001
    #>
    #> Model Results:
    #>
    #> estimate      se     zval    pval    ci.lb    ci.ub
    #>  -0.7145  0.1798  -3.9744  <.0001  -1.0669  -0.3622  ***
    #>
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The matching Bayesian fit reuses the same effect-size columns. `measure = "RR"` tells `brma()` to use the default prior distributions for log risk ratios. The default prior distributions are weakly informative and provide slight regularization toward zero.

    fit1_brma <- brma(
      yi = yi, vi = vi, measure = "RR",
      data = dat, seed = 1
    )

`summary()` reports posterior means, credible intervals, and (suppressed here) MCMC convergence diagnostics. Convergence diagnostics are included by default; chain-level visual diagnostics are available via `plot_diagnostic_trace()`, `plot_diagnostic_density()`, and `plot_diagnostic_autocorrelation()`.

    summary(fit1_brma, include_mcmc_diagnostics = FALSE)
    #>
    #> Bayesian Random-Effects Model (k = 13)
    #>
    #> Estimates
    #>       Mean    SD  0.025    0.5  0.975
    #> mu  -0.688 0.186 -1.056 -0.682 -0.292
    #> tau  0.563 0.136  0.334  0.552  0.858

The posterior mean for the log risk ratio is close to the REML estimate from the `metafor` package; the slight pull toward zero reflects the default `Normal(0, 1)` prior distribution. `pooled_effect()` returns the same summary, backtransformed to the risk-ratio scale via `transform = "EXP"`:

    pooled_effect(fit1_brma, transform = "EXP")
    #>
    #> Pooled Effect Size (risk ratio)
    #>     Mean Median 0.025 0.975
    #> mu 0.511  0.506 0.348 0.747
    #> Effect estimates are summarized on the risk ratio scale using EXP on the log risk ratio measure.

`plot()` visualizes the posterior distribution, and `prior = TRUE` overlays the prior distribution so the shift from prior to posterior is visible.

    par(mar = c(4, 4, 1, 1))
    plot(fit1_brma, parameter = "mu", prior = TRUE, xlim = c(-3, 3))

## Heterogeneity

`metafor::confint()` reports profile-likelihood confidence intervals for `tau^2`, `tau`, `I^2`, and `H^2`. `summary_heterogeneity()` is the Bayesian counterpart, returning posterior summaries of the same quantities.

    confint(fit1_metafor)
    #>
    #>        estimate   ci.lb   ci.ub
    #> tau^2    0.3132  0.1197  1.1115
    #> tau      0.5597  0.3460  1.0543
    #> I^2(%)  92.2214 81.9177 97.6781
    #> H^2     12.8558  5.5303 43.0680

    summary_heterogeneity(fit1_brma)
    #>
    #> Heterogeneity Estimates:
    #>        Mean Median  0.025  0.975
    #> tau   0.563  0.552  0.334  0.858
    #> tau2  0.336  0.305  0.112  0.735
    #> I2   91.153 92.030 80.883 96.532
    #> H2   13.712 12.547  5.231 28.832

## Funnel Plot

`funnel()` plots observed effect sizes against the fitted sampling distribution. We turn off the `sampling_heterogeneity` argument, since the `RoBMA` R package shows the full sampling distribution under the model by default.

    par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
    metafor::funnel(fit1_metafor, main = "metafor", xlim = c(-2, 1), ylim = c(0, 0.8))
    funnel(fit1_brma, main = "RoBMA", xlim = c(-2, 1), ylim = c(0, 0.8), sampling_heterogeneity = FALSE)

## Meta-Regression with a Continuous Moderator

A continuous moderator can be added with the standard formula syntax via the `mods` argument. `ablat` is the absolute latitude of the trial, frequently used to explain BCG effect-size variation.

    fit2_metafor <- metafor::rma(yi, vi, mods = ~ ablat, data = dat)
    fit2_metafor
    #>
    #> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
    #>
    #> tau^2 (estimated amount of residual heterogeneity):     0.0764 (SE = 0.0591)
    #> tau (square root of estimated tau^2 value):             0.2763
    #> I^2 (residual heterogeneity / unaccounted variability): 68.39%
    #> H^2 (unaccounted variability / sampling variability):   3.16
    #> R^2 (amount of heterogeneity accounted for):            75.62%
    #>
    #> Test for Residual Heterogeneity:
    #> QE(df = 11) = 30.7331, p-val = 0.0012
    #>
    #> Test of Moderators (coefficient 2):
    #> QM(df = 1) = 16.3571, p-val < .0001
    #>
    #> Model Results:
    #>
    #>          estimate      se     zval    pval    ci.lb    ci.ub
    #> intrcpt    0.2515  0.2491   1.0095  0.3127  -0.2368   0.7397
    #> ablat     -0.0291  0.0072  -4.0444  <.0001  -0.0432  -0.0150  ***
    #>
    #> ---
    #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    fit2_brma <- brma(
      yi = yi, vi = vi, measure = "RR",
      mods = ~ ablat,
      data = dat, seed = 1
    )

    summary(fit2_brma, include_mcmc_diagnostics = FALSE)
    #>
    #> Bayesian Mixed-Effect Model (k = 13)
    #>
    #> Common Estimates
    #>      Mean    SD 0.025   0.5 0.975
    #> tau 0.303 0.131 0.088 0.289 0.587
    #>
    #> Meta-Regression
    #>             Mean    SD  0.025    0.5  0.975
    #> intercept  0.217 0.275 -0.362  0.237  0.712
    #> ablat     -0.028 0.008 -0.043 -0.028 -0.011

`regplot()` produces the familiar bubble plot with the fitted regression line.

    par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
    metafor::regplot(fit2_metafor, main = "metafor", xlim = c(10, 60), ylim = c(-2, 0.5))
    regplot(fit2_brma, main = "RoBMA", xlim = c(10, 60), ylim = c(-2, 0.5))

`predict()` returns posterior summaries for user-specified moderator values.

    new_data <- data.frame(ablat = c(15, 35, 55))
    predict(fit2_brma, newdata = new_data, type = "terms")
    #>
    #> Location Term Posterior Prediction:
    #>         Mean Median  0.025  0.975
    #> mu[1] -0.196 -0.182 -0.587  0.123
    #> mu[2] -0.746 -0.743 -0.994 -0.506
    #> mu[3] -1.297 -1.307 -1.722 -0.835

## Comparing Models with Bayes Factors

Marginal likelihoods enable model comparison via Bayes factors (Kass & Raftery, 1995). After attaching marginal likelihoods to both fits with `add_marglik()`, `bf()` returns the Bayes factor for one model against another.

    fit1_brma <- add_marglik(fit1_brma)
    fit2_brma <- add_marglik(fit2_brma)

    bf(fit2_brma, fit1_brma)
    #> Bayes factor in favor of Model 1 over Model 2: 14.76119

The Bayes factor quantifies how much the data shift the relative predictive support between the two models. Here we obtain strong evidence in favor of the model including the `ablat` predictor.

## Adding a Categorical Moderator

`alloc` is a three-level factor describing the trial allocation method (random, alternate, systematic). The same formula syntax extends the meta-regression.

    fit3_metafor <- metafor::rma(yi, vi, mods = ~ ablat + alloc, data = dat)

    fit3_brma <- brma(
      yi = yi, vi = vi, measure = "RR",
      mods = ~ ablat + alloc,
      data = dat, seed = 1
    )

    summary(fit3_brma, include_mcmc_diagnostics = FALSE)
    #>
    #> Bayesian Mixed-Effect Model (k = 13)
    #>
    #> Common Estimates
    #>      Mean    SD 0.025   0.5 0.975
    #> tau 0.344 0.144 0.100 0.330 0.674
    #>
    #> Meta-Regression
    #>                     Mean    SD  0.025    0.5  0.975
    #> intercept          0.287 0.355 -0.468  0.312  0.929
    #> ablat             -0.026 0.009 -0.042 -0.027 -0.008
    #> alloc[random]     -0.249 0.251 -0.749 -0.248  0.251
    #> alloc[systematic]  0.003 0.276 -0.489 -0.017  0.580

The `RoBMA` R package also allows direct visualization of categorical moderators in the bubble plot (the `metafor` package only supports a single dummy coefficient at the moment).

    par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
    metafor::regplot(fit3_metafor, mod = "allocrandom", main = "metafor: random dummy", ylim = c(-2, 0.5))
    regplot(fit3_brma, mod = "alloc", main = "RoBMA: alloc factor", ylim = c(-2, 0.5))

`marginal_means()` computes estimated marginal means for each moderator term, averaging over the other moderators in the fit. Continuous moderators (here `ablat`) are evaluated at \(-1\), \(0\), and \(+1\) standard deviations of the predictor; categorical moderators (here `alloc`) at each level.

    emm <- marginal_means(fit3_brma)
    summary(emm)
    #>
    #> Marginal Means:
    #>                     Mean    SD  0.025    0.5  0.975
    #> intercept         -0.592 0.217 -1.061 -0.585 -0.180
    #> ablat[-1SD]       -0.212 0.247 -0.737 -0.197  0.234
    #> ablat[0SD]        -0.592 0.217 -1.061 -0.585 -0.180
    #> ablat[1SD]        -0.971 0.254 -1.494 -0.971 -0.429
    #> alloc[alternate]  -0.592 0.217 -1.061 -0.585 -0.180
    #> alloc[random]     -0.841 0.177 -1.192 -0.833 -0.516
    #> alloc[systematic] -0.588 0.214 -0.999 -0.605 -0.154

The `plot()` function visualizes the posterior distribution of the estimated marginal means.

    par(mar = c(4, 4, 1, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
    plot(emm, parameter = "alloc", xlim = c(-2, 0.5), lwd = 2)

## Model Diagnostics

The `RoBMA` R package implements a wide set of model fit diagnostics. Here, we highlight only a few.

`metafor::rstudent()` returns frequentist studentized residuals computed by deletion; `RoBMA::rstudent()` returns LOO-PIT residuals, which propagate posterior uncertainty and leverage through leave-one-out prediction. The two definitions are comparable but not identical; the `RoBMA` R package’s LOO-PIT residuals first require LOO computation via `add_loo()`.

    fit3_brma <- add_loo(fit3_brma)

    head(rstudent(fit3_metafor)$z)
    #> [1]  0.4369227 -0.2024807 -0.3084436 -0.1187827 -0.4726194  0.4726194
    head(rstudent(fit3_brma)$z)
    #> [1]  0.3515262 -0.3483678 -0.3940142 -0.4052585 -0.5379945  0.3420329

`qqnorm()` checks the residuals against a standard normal reference.

    par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
    qqnorm(fit3_metafor, main = "metafor: studentized", xlim = c(-1.5, 1.5), ylim = c(-3, 3))
    qqnorm(fit3_brma, main = "RoBMA: LOO-PIT", xlim = c(-1.5, 1.5), ylim = c(-3, 3))

`influence()` collects DFFITS, Cook’s distance, COVRATIO, leverage, and DFBETAS where the diagnostic is available.

    inf_brma <- influence(fit3_brma)
    head(inf_brma$inf)
    #>           rstudent      dffits      cook.d    cov.r   tau.del        hat
    #> Study 1  0.3515262  0.09742742 0.003051713 1.352210 0.3611943 0.08676324
    #> Study 2 -0.3483678 -0.16031362 0.009344464 1.672075 0.3683854 0.19988035
    #> Study 3 -0.3940142 -0.09201147 0.002596289 1.247964 0.3566583 0.06776189
    #> Study 4 -0.4052585 -0.25448939 0.028780571 2.209083 0.3775556 0.42974702
    #> Study 5 -0.5379945 -0.24916951 0.024724592 1.683551 0.3457254 0.53991804
    #> Study 6  0.3420329  0.21231431 0.022986373 2.516078 0.3785678 0.69164774

## References

Bartoš, F., & Maier, M. (2020). RoBMA: An R package for robust Bayesian meta-analyses. https://CRAN.R-project.org/package=RoBMA
Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773–795. https://doi.org/10.1080/01621459.1995.10476572
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03