---
title: "Bayesian Meta-Analysis"
author: "František Bartoš"
date: "28th of April 2026"
output:
  rmarkdown::html_vignette:
    self_contained: yes
bibliography: ../inst/REFERENCES.bib
csl: ../inst/apa.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{Bayesian Meta-Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown_notangle}
---

```{r child = "_vignette-nowrap.md", echo = FALSE, eval = TRUE}
```

```{r setup, include = FALSE}
source("vignette-cache.R", local = knitr::knit_global())

cached_fits <- vignette_cache(
  name     = "v02-bayesian-meta-analysis",
  objects  = c("fit1_brma", "fit2_brma", "fit3_brma"),
  packages = c("metafor", "metadat")
)

knitr::opts_chunk$set(
  collapse    = TRUE,
  comment     = "#>",
  eval        = vignette_cache_eval(cached_fits),
  message     = FALSE,
  warning     = FALSE,
  fig.width   = 7,
  fig.height  = 3.5,
  dev         = "png",
  fig.retina  = 3
)
if (.Platform$OS.type == "windows") {
  knitr::opts_chunk$set(dev.args = list(type = "cairo"))
}
```

```{r load-models, include = FALSE}
library("metafor")
library("RoBMA")

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

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

vignette_cache_load(cached_fits)
```

```{r precompute-models, include = FALSE, eval = FALSE}
library("metafor")
library("RoBMA")

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

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

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

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

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

vignette_cache_save(cached_fits)
```

This vignette is a starting point for using the `RoBMA` R package [@RoBMA].
We walk through the BCG vaccine dataset from the `metafor` package [@metafor], 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*](v01-prior-distributions.html) 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`.

```{r escalc}
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.

```{r fit1-metafor}
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.

```{r fit1-brma, eval = FALSE}
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()`.

```{r fit1-brma-summary}
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"`:

```{r fit1-brma-pooled}
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.

```{r fit1-brma-plot, fig.width = 5, fig.height = 3.5}
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.

```{r heterogeneity-confint}
confint(fit1_metafor)
```

```{r heterogeneity-bayes}
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.

```{r funnel-plot, fig.width = 7, fig.height = 3.2}
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.

```{r fit2-metafor}
fit2_metafor <- metafor::rma(yi, vi, mods = ~ ablat, data = dat)
fit2_metafor
```

```{r fit2-brma, eval = FALSE}
fit2_brma <- brma(
  yi = yi, vi = vi, measure = "RR",
  mods = ~ ablat,
  data = dat, seed = 1
)
```
```{r fit2-brma-summary}
summary(fit2_brma, include_mcmc_diagnostics = FALSE)
```

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

```{r regplot-continuous, fig.width = 7, fig.height = 3.2}
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.

```{r fit2-predict}
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 [@kass1995bayes].
After attaching marginal likelihoods to both fits with `add_marglik()`, `bf()` returns the Bayes factor for one model against another.

```{r add-marglik, eval = FALSE}
fit1_brma <- add_marglik(fit1_brma)
fit2_brma <- add_marglik(fit2_brma)
```
```{r bayes-factor}
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.

```{r fit3-metafor}
fit3_metafor <- metafor::rma(yi, vi, mods = ~ ablat + alloc, data = dat)
```

```{r fit3-brma, eval = FALSE}
fit3_brma <- brma(
  yi = yi, vi = vi, measure = "RR",
  mods = ~ ablat + alloc,
  data = dat, seed = 1
)
```
```{r fit3-brma-summary}
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).

```{r regplot-categorical, fig.width = 7, fig.height = 3.2}
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.

```{r marginal-means}
emm <- marginal_means(fit3_brma)
summary(emm)
```

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

```{r marginal-means-plot, fig.width = 4.5, fig.height = 3.25}
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()`.

```{r add-loo, eval = FALSE}
fit3_brma <- add_loo(fit3_brma)
```
```{r residuals}
head(rstudent(fit3_metafor)$z)
head(rstudent(fit3_brma)$z)
```

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

```{r qqplot, fig.width = 7, fig.height = 3.2}
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.

```{r influence}
inf_brma <- influence(fit3_brma)
head(inf_brma$inf)
```

## References