Bayesian Meta-Analysis

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

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)

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")

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)

summary_heterogeneity(fit1_brma)

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

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

summary(fit2_brma, include_mcmc_diagnostics = FALSE)

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")

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)

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)

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)

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)
head(rstudent(fit3_brma)$z)

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)