---
title: "Posterior Densities and Parameter Tests"
author: "František Bartoš"
date: "25th of October 2026"
output:
  rmarkdown::html_vignette:
    self_contained: yes
bibliography: "`r normalizePath('../inst/REFERENCES.bib', winslash = '/', mustWork = TRUE)`"
csl: "`r normalizePath('../inst/apa.csl', winslash = '/', mustWork = TRUE)`"
link-citations: true
vignette: >
  %\VignetteIndexEntry{Posterior Densities and Parameter Tests}
  %\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     = "v22-posterior-densities-parameter-tests",
  objects  = c(
    "fit_simple",
    "fit_simple_mu_null",
    "fit_simple_mu_negative",
    "fit_simple_mu_positive",
    "fit_meta_reg"
  ),
  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  = vignette_fig_retina()
)
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
)

prior_effect_normal   <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1))
prior_effect_null     <- prior(distribution = "spike",  parameters = list(location = 0))
prior_effect_negative <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = -Inf, upper = 0))
prior_effect_positive <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = 0, upper = Inf))

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
)

prior_effect_normal   <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1))
prior_effect_null     <- prior(distribution = "spike", parameters = list(location = 0))
prior_effect_negative <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = -Inf, upper = 0))
prior_effect_positive <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = 0, upper = Inf))

fit_simple <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_normal,
  data = dat, seed = 1
)
fit_simple <- add_marglik(fit_simple)

fit_simple_mu_null <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_null,
  data = dat, seed = 2
)
fit_simple_mu_null <- add_marglik(fit_simple_mu_null)

fit_simple_mu_negative <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_negative,
  data = dat, seed = 3
)
fit_simple_mu_negative <- add_marglik(fit_simple_mu_negative)

fit_simple_mu_positive <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_positive,
  data = dat, seed = 4
)
fit_simple_mu_positive <- add_marglik(fit_simple_mu_positive)

fit_meta_reg <- brma(
  yi = yi, vi = vi, measure = "RR",
  mods = ~ ablat + alloc,
  set_contrast_factor_predictors = "meandif",
  prior_effect = prior_effect_normal,
  data = dat, seed = 5
)

vignette_cache_save(cached_fits)
```

Bayesian meta-analysis produces posterior distributions for the average effect, the between-study heterogeneity, and any meta-regression coefficients.
Evidence for or against specific values of these parameters can be quantified via Bayes factors.

The posterior distributions are usually summarized via posterior means and credible intervals, and visualized via (prior and) posterior density plots, usually obtained via kernel density estimation (KDE).
The KDE visualization can be, however, quite "bumpy" even with a sizable number of posterior samples.
In the first part of this vignette we illustrate how the posterior density visualization can be improved via q-grid conditional marginal density estimator (qCMDE) implemented in the `RoBMA` R package [@RoBMA].

In the second part of the vignette, we illustrate how to compute Bayes factors based on posterior distribution in the `RoBMA` R package (which takes advantage of the qCMDE estimator when testing a point null hypothesis).
This posterior samples route allows us to re-use a single model fit with a Savage--Dickey density ratio [@dickey1971weighted; @wagenmakers2010bayesian] to compute Bayes factor for point-null hypotheses and encompassing-prior odds [@klugkist2005bayesian] for inequality and region hypotheses.
It adds to the existing bridge sampling route [@gronau2017tutorial] which compares marginal likelihoods of separately fitted constrained and unconstrained models (see "Comparing Models with Bayes Factors" in [*Bayesian Meta-Analysis*](v02-bayesian-meta-analysis.html)),
and the product-space route [@lodewyckx2011tutorial] which fits multiple models simultaneously in a single MCMC run and estimates the Bayes factor from the posterior inclusion probabilities (see [*Bayesian Model Averaging*](v20-bayesian-model-averaging.html)).

We use the BCG vaccine dataset from the `metafor` package [@metafor] to walk through the density estimation and Bayes factor computation.
We start with a simple random-effects meta-analysis, compare KDE and qCMDE on the average effect and heterogeneity, and use `hypothesis()` to test three hypotheses on the average effect with results matched against bridge sampling.
We then add a continuous and a categorical moderator and combine `marginal_means()` with `hypothesis()` for tests on estimated marginal means.
All effect sizes are on the log risk-ratio scale, so the no-effect value is 0.

## Simple Meta-Analysis

The BCG vaccine dataset 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
)
```

We first fit a random-effects meta-analysis with a standard normal prior distribution on the mean effect.

```{r fit-simple, eval = FALSE}
prior_effect_normal <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1))

fit_simple <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_normal,
  data = dat, seed = 1
)
```

```{r summary-simple}
summary(fit_simple)
```

## qCMDE Posterior Density Estimation

By default, `plot.brma()` visualizes marginal posterior distributions with kernel density estimation (KDE) on the sampled MCMC draws.
KDE summaries are convenient and generic, but they might not look smooth even with a sizable number of posterior samples, and the density value at a specific point, for example `mu = 0`, inherits the same Monte Carlo noise.

qCMDE is a model-aware alternative that leverages information in the joint posterior distribution.
It provides a smoother density curve and a more stable estimate of the marginal posterior density at any user-specified value.
The construction extends the Rao--Blackwellized marginal density estimator of @gelfand1992bayesian and the importance-weighted marginal density estimator of @chen1994importance.

Passing `density_method = "qCMDE"` to `plot.brma()` switches the density route from the default KDE to qCMDE.
We compare the two estimators (KDE in black, qCMDE in blue) for the average effect (`mu`) and the between-study heterogeneity (`tau`); the grey curves overlay the corresponding prior distribution.

```{r simple-mu-tau-densities, fig.height = 6}
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)

plot(fit_simple, parameter = "mu",  prior = TRUE, density_method = "KDE",   main = "mu: KDE",   xlim = c(-1.8, 0.6))
plot(fit_simple, parameter = "mu",  prior = TRUE, density_method = "qCMDE", density_control = list(samples = 1000L), main = "mu: qCMDE", xlim = c(-1.8, 0.6))
plot(fit_simple, parameter = "tau", prior = TRUE, density_method = "KDE",   main = "tau: KDE",   xlim = c(0, 1.2))
plot(fit_simple, parameter = "tau", prior = TRUE, density_method = "qCMDE", density_control = list(samples = 1000L), main = "tau: qCMDE", xlim = c(0, 1.2))
```

The qCMDE curves are smoother on both panels, but otherwise track the KDE summary in shape and location.

## Testing the Mean Effect

Bayes factors compare two models by their prior predictive performance, usually through the ratio of marginal likelihoods [@kass1995bayes].
For point-null hypotheses on a single parameter, the same Bayes factor can be obtained from the Savage--Dickey density ratio [@dickey1971weighted; @wagenmakers2010bayesian], which needs only the prior and posterior density ordinate at the tested value.
For inequality and region hypotheses, the encompassing-prior approach computes the Bayes factor as the change from prior odds to posterior odds for the restricted region [@klugkist2005bayesian].
The Savage--Dickey and encompassing-prior routes are both available from a single fit through `hypothesis()`; the rest of this section walks through three example hypotheses on the average effect and contrasts both routes:

- $\text{BF}_{10}$: an unrestricted average effect against a point null at 0;
- $\text{BF}_{-0}$: a negative average effect against a point null at 0;
- $\text{BF}_{-+}$: a negative average effect against a positive average effect.

```{r bf-helper, include = FALSE}
compare_bf <- function(posterior_bf, fit, fit_null, fit_negative,
                       fit_positive) {

  test_labels <- c("$\\text{BF}_{10}$", "$\\text{BF}_{-0}$", "$\\text{BF}_{-+}$")
  bridge_bf <- c(
    bf(fit, fit_null)[["bf"]],
    bf(fit_negative, fit_null)[["bf"]],
    bf(fit_negative, fit_positive)[["bf"]]
  )
  posterior_raw   <- posterior_bf[["BF"]]

  out <- data.frame(
    Test                     = test_labels,
    `Posterior BF`           = format_bf_table(posterior_raw),
    `Bridge BF`              = format_bf_table(bridge_bf),
    check.names              = FALSE
  )
  row.names(out) <- NULL

  return(out)
}
format_bf_table <- function(x) {

  out <- formatC(
    round(as.numeric(x), 1),
    format = "f",
    digits = 1
  )
  out[!is.finite(as.numeric(x))] <- as.character(as.numeric(x)[!is.finite(as.numeric(x))])

  return(out)
}
```

The posterior-based route uses `hypothesis()`.
With `density_method = "qCMDE"`, the point-null hypotheses use the qCMDE ordinate at 0 which improves stability of the estimate.
The region test $\text{BF}_{-+}$ uses prior-posterior odds and does not need a density ordinate.

```{r hypothesis-mu-posterior}
mu_posterior_bf <- hypothesis(
  fit_simple, density_method = "qCMDE",
  hypothesis     = c(
    "mu != 0",
    "mu < 0 vs mu = 0",
    "mu < 0 vs mu > 0"
  ),
)
mu_posterior_bf
```

The bridge sampling route fits one constrained model per hypothesis side, attaches marginal likelihoods with `add_marglik()`, and compares them with `bf()`.

```{r bridge-model-syntax, eval = FALSE}
prior_effect_null     <- prior(distribution = "spike",  parameters = list(location = 0))
prior_effect_negative <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = -Inf, upper = 0))
prior_effect_positive <- prior(distribution = "normal", parameters = list(mean = 0, sd = 1), truncation = list(lower = 0, upper = Inf))

fit_simple <- add_marglik(fit_simple)

fit_simple_mu_null <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_null,
  data = dat, seed = 2
)
fit_simple_mu_null <- add_marglik(fit_simple_mu_null)

fit_simple_mu_negative <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_negative,
  data = dat, seed = 3
)
fit_simple_mu_negative <- add_marglik(fit_simple_mu_negative)

fit_simple_mu_positive <- brma(
  yi = yi, vi = vi, measure = "RR",
  prior_effect = prior_effect_positive,
  data = dat, seed = 4
)
fit_simple_mu_positive <- add_marglik(fit_simple_mu_positive)
```

```{r bridge-bf-syntax, eval = FALSE}
bf_10         <- bf(fit_simple, fit_simple_mu_null)
bf_minus_0    <- bf(fit_simple_mu_negative, fit_simple_mu_null)
bf_minus_plus <- bf(fit_simple_mu_negative, fit_simple_mu_positive)
```

The two routes answer the same questions through different computations.
Bridge sampling compares separately fitted constrained models, while the posterior-based route reuses one fitted model and only needs the qCMDE ordinate at the tested value.
The table reports the resulting Bayes factors and the corresponding error estimates in percentages.

```{r hypothesis-mu-comparison, echo = FALSE}
knitr::kable(
  compare_bf(
    posterior_bf = mu_posterior_bf,
    fit          = fit_simple,
    fit_null     = fit_simple_mu_null,
    fit_negative = fit_simple_mu_negative,
    fit_positive = fit_simple_mu_positive
  ),
  caption   = "Bayes factor estimates for the mean-effect hypotheses.",
  align     = c("l", "r", "r"),
  escape    = FALSE,
  row.names = FALSE
)
```

The point-null Bayes factors agree closely across the two routes.
The posterior-based value for $\text{BF}_{-+}$ is limited by the number of posterior draws as only a very small proportion of samples is at the positive side.
Consequently, the encompassing prior estimate is much more variable than the corresponding bridge sampling approach (see the the large `error%(BF)` reported earlier).
In general, Bayes factors larger than 100 (lower than 1/100) are likely to have limited precision when using the posterior based or product-space approaches.
(This usually does not matter for the interpretation, i.e., the data provide overwhelming evidence for one of the hypotheses.)

## Meta-Regression

We next add absolute latitude (`ablat`) and allocation method (`alloc`) as moderators.
For the categorical moderator we use mean-difference contrasts via `set_contrast_factor_predictors = "meandif"`.
It specifies the prior distribution such that the allocation levels are treated as exchangeable;
the prior distribution does not depend on an arbitrary reference level and is better suited to hypothesis tests comparing allocation methods among each other.

```{r fit-meta-reg, eval = FALSE}
fit_meta_reg <- brma(
  yi = yi, vi = vi, measure = "RR",
  mods = ~ ablat + alloc,
  set_contrast_factor_predictors = "meandif",
  prior_effect = prior_effect_normal,
  data = dat, seed = 5
)
```

```{r summary-meta-reg}
summary(fit_meta_reg)
```

qCMDE is again available directly through `plot.brma()`.
When interpreting categorical moderators, estimated marginal means (EMMs) are a clearer way to report the results.
The `marginal_means()` allows us to compute Bayes factors using `bf = TRUE` and use the qCMDE via `density_method = "qCMDE"` arguments.
The density estimate is directly saved in the EMMs object, so the corresponding hypotheses can be tested later with `hypothesis()`.

```{r meta-reg-marginal-means}
emm_meta_reg_default <- marginal_means(fit_meta_reg)
emm_meta_reg_qcmde   <- marginal_means(fit_meta_reg, bf = TRUE, density_method = "qCMDE")
```

The four panels below compare KDE and qCMDE for the `ablat` coefficient and for the per-level posterior of the three `alloc` marginal means.

```{r meta-reg-densities, fig.height = 5}
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)

plot(fit_meta_reg, parameter = "ablat", prior = TRUE, density_method = "KDE",   main = "KDE",   xlim = c(-0.08, 0.03))
plot(fit_meta_reg, parameter = "ablat", prior = TRUE, density_method = "qCMDE", main = "qCMDE", xlim = c(-0.08, 0.03))
plot(emm_meta_reg_default, parameter = "alloc", prior = TRUE, main = "alloc means: KDE",   xlim = c(-1.6, 0.3))
plot(emm_meta_reg_qcmde,   parameter = "alloc", prior = TRUE, main = "alloc means: qCMDE", xlim = c(-1.6, 0.3))
```

Hypotheses on the categorical moderator can then be stated directly in terms of named allocation-method marginal means.
The bracket syntax (`alloc[alternate]`) selects a named level from the allocation marginal posterior.

```{r hypothesis-alloc}
hypothesis(
  emm_meta_reg_qcmde,
  hypothesis = c(
    "alloc[alternate] > alloc[random]",
    "alloc[alternate] > alloc[systematic]",
    "alloc[random] > alloc[systematic]"
  )
)
```


## References