Generalized Linear Mixed-Effects Meta-Analysis

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 alternate

metafor::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 ' ' 1

The 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.922

The 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.