Generalized Linear Mixed-Effects Meta-Analysis
František Bartoš
28th of April 2026
Source:vignettes/v13-metafor-parity-glmm.Rmd
v13-metafor-parity-glmm.RmdThis vignette illustrates how to use the RoBMA R package
(Bartoš & Maier,
2020) when the outcome is a count rather than an effect-size
estimate. We walk through the BCG vaccine dataset from the
metafor package (Viechtbauer, 2010), showing the
brma.glmm() call alongside its
metafor::rma.glmm() counterpart for log odds ratios. The
interface mirrors metafor::rma.glmm(): counts
(ai, bi, ci, di) are
passed in directly, and no escalc() step is needed.
All brma.glmm() calls keep the default prior
distributions; see the Prior
Distributions vignette for the prior definitions and
customization options. For the non-GLMM workflow showing more features
(also applicable to GLMM models) see the Bayesian
Meta-Analysis vignette.
brma.glmm() currently supports binomial outcomes (log
odds ratio, measure = "OR") and Poisson outcomes (log
incidence rate ratio, measure = "IRR").
Binomial Random-Effects Model
The BCG vaccine dataset from the metadat package
contains 13 randomized trials on tuberculosis prevention. Treated and
control event counts (tpos, tneg,
cpos, cneg) are passed in directly.
data("dat.bcg", package = "metadat")
head(dat.bcg)
#> trial author year tpos tneg cpos cneg ablat alloc
#> 1 1 Aronson 1948 4 119 11 128 44 random
#> 2 2 Ferguson & Simes 1949 6 300 29 274 55 random
#> 3 3 Rosenthal et al 1960 3 228 11 209 42 random
#> 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random
#> 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate
#> 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternatemetafor::rma.glmm() fits the random-effects logistic
model. We use model = "UM.FS" to obtain the unconditional
model with fixed study-specific intercepts, the formulation that
corresponds to the binomial GLMM in brma.glmm() (Model 4 in
Jackson et al. (2018)).
fit1_metafor <- metafor::rma.glmm(
ai = tpos, bi = tneg, ci = cpos, di = cneg, measure = "OR",
model = "UM.FS", data = dat.bcg
)
fit1_metafor
#>
#> Random-Effects Model (k = 13; tau^2 estimator: ML)
#> Model Type: Unconditional Model with Fixed Study Effects
#>
#> tau^2 (estimated amount of total heterogeneity): 0.2949
#> tau (square root of estimated tau^2 value): 0.5430
#> I^2 (total heterogeneity / total variability): 91.02%
#> H^2 (total variability / sampling variability): 11.14
#>
#> Tests for Heterogeneity:
#> Wld(df = 12) = 163.1649, p-val < .0001
#> LRT(df = 12) = 176.9544, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> -0.7450 0.1756 -4.2434 <.0001 -1.0891 -0.4009 ***
#>
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The matching Bayesian fit takes the same count columns and
measure = "OR".
fit1_brma <- brma.glmm(
ai = tpos, bi = tneg, ci = cpos, di = cneg, measure = "OR",
data = dat.bcg, seed = 1
)summary() reports posterior means, credible intervals,
and (suppressed here) MCMC convergence diagnostics for the pooled log
odds ratio mu and the heterogeneity tau.
summary(fit1_brma, include_mcmc_diagnostics = FALSE)
#>
#> Bayesian Random-Effects Model (k = 13)
#>
#> Estimates
#> Mean SD 0.025 0.5 0.975
#> mu -0.719 0.183 -1.089 -0.718 -0.353
#> tau 0.589 0.141 0.361 0.569 0.922The posterior mean for mu tracks the maximum-likelihood
estimate from the metafor package, with comparable
heterogeneity tau. pooled_effect() returns the
same summary backtransformed to the odds-ratio scale via
transform = "EXP":
pooled_effect(fit1_brma, transform = "EXP")
#>
#> Pooled Effect Size (odds ratio)
#> Mean Median 0.025 0.975
#> mu 0.495 0.488 0.337 0.703
#> Effect estimates are summarized on the odds ratio scale using EXP on the log odds ratio measure.Default Prior Distributions and the Auxiliary Parameter
The default prior distributions on mu and
tau follow the unit-information set-up as elsewhere in the
package, scaled by the known unit information standard deviation for the
chosen measure (sqrt(4) for OR
and IRR). Beyond mu and tau, GLMM
models include estimate-specific nuisance parameters for the per-study
baseline level. For binomial outcomes (measure = "OR"),
logit(pi[i]) is the midpoint of the two arm logits before
applying half of the study-specific log odds ratio to each arm. The
default prior_baserate = Beta(1, 1) is placed independently
on each pi[i]. For Poisson outcomes
(measure = "IRR"), phi[i] is the midpoint of
the two arm log incidence rates before applying half of the
study-specific log incidence rate ratio to each arm. Because this
midpoint log-rate is unbounded, there is no bounded flat default
analogous to Beta(1, 1), so the default
prior_lograte is a proper data-based unit-information
normal used independently for each phi[i]. Thus the default
nuisance-parameter prior distributions are exchangeable in the iid
sense, but they are not hierarchical prior distributions with a learned
common baseline-rate distribution. Both defaults can be overridden via
the corresponding prior_* argument; see the Prior Distributions
vignette. The Poisson model corresponds to Bagos
& Nikolopoulos (2009) and
uses event counts and exposure times (x1i,
t1i, x2i, t2i) in place of the
four-way count layout.
Other Inference Helpers
All inference helpers available for any other brma() fit
work the same way for a brma.glmm() fit. Meta-regression
via mods, multilevel structures via cluster,
location-scale models via scale, posterior summaries,
plots, predictions, residuals, influence, and LOO comparisons all carry
over. See the Bayesian
Meta-Analysis vignette for the full walkthrough and the Multilevel
Meta-Analysis vignette for the multilevel extension; the only
thing that changes here is the count-based input and the auxiliary
baseline-rate prior.