Multilevel Random-Effects Model
The school-calendar dataset from the metadat package
contains 56 standardized mean differences nested in 11 school districts.
Effect sizes yi and sampling variances vi are
already provided, so no escalc() step is needed.
data("dat.konstantopoulos2011", package = "metadat")
dat <- dat.konstantopoulos2011
head(dat)
#>
#> district school study year yi vi
#> 1 11 1 1 1976 -0.180 0.118
#> 2 11 2 2 1976 -0.220 0.118
#> 3 11 3 3 1976 0.230 0.144
#> 4 11 4 4 1976 -0.300 0.144
#> 5 12 1 5 1989 0.130 0.014
#> 6 12 2 6 1989 -0.260 0.014
metafor::rma.mv() fits a multilevel random-effects model
with the compound-symmetric structure
random = ~ school | district, parameterized by the total
heterogeneity \(\tau\) and the
within-district correlation \(\rho\).
fit1_metafor <- metafor::rma.mv(yi, vi, random = ~ school | district, data = dat)
fit1_metafor
#>
#> Multivariate Meta-Analysis Model (k = 56; method: REML)
#>
#> Variance Components:
#>
#> outer factor: district (nlvls = 11)
#> inner factor: school (nlvls = 11)
#>
#> estim sqrt fixed
#> tau^2 0.0978 0.3127 no
#> rho 0.6653 no
#>
#> Test for Heterogeneity:
#> Q(df = 55) = 578.8640, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> 0.1847 0.0846 2.1845 0.0289 0.0190 0.3504 *
#>
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The matching Bayesian fit replaces the random = ...
argument with cluster = district. brma()
parameterizes the multilevel model in the same compound-symmetric form:
a total heterogeneity \(\tau\) and a
cluster-level allocation \(\rho\) where
\(\tau_{between} = \tau\sqrt{\rho}\)
and \(\tau_{within} = \tau\sqrt{1 -
\rho}\). brma(..., cluster = ...) is the specialized
interface for this model. measure = "SMD" selects the
default prior distributions for standardized mean differences.
fit1_brma <- brma(
yi = yi, vi = vi, measure = "SMD",
cluster = district,
data = dat, seed = 1
)
summary() reports posterior means, credible intervals,
and (suppressed here) MCMC convergence diagnostics for mu,
tau, and rho.
summary(fit1_brma, include_mcmc_diagnostics = FALSE)
#>
#> Bayesian Multilevel Random-Effects Model (k = 56, clusters = 11)
#>
#> Estimates
#> Mean SD 0.025 0.5 0.975
#> mu 0.186 0.085 0.022 0.186 0.356
#> tau 0.318 0.056 0.232 0.310 0.450
#> rho 0.627 0.139 0.327 0.639 0.858
The posterior mean for the pooled effect tracks the REML estimate
from the metafor package, with comparable cluster-level
allocation rho and total heterogeneity
tau.
Heterogeneity: \(\tau\) and \(\rho\)
metafor::confint() reports profile-likelihood confidence
intervals for tau^2, tau, and
rho.
confint(fit1_metafor)
#>
#> estimate ci.lb ci.ub
#> tau^2 0.0978 0.0528 0.2398
#> tau 0.3127 0.2298 0.4897
#>
#> estimate ci.lb ci.ub
#> rho 0.6653 0.3282 0.8855
summary_heterogeneity() is the Bayesian counterpart,
reporting posterior summaries of the same quantities along with the
partitioned within- and between-cluster components.
summary_heterogeneity(fit1_brma)
#>
#> Heterogeneity Estimates:
#> Mean Median 0.025 0.975
#> tau 0.318 0.310 0.232 0.450
#> rho 0.627 0.639 0.327 0.858
#> tau [within] 0.187 0.185 0.131 0.253
#> tau [between] 0.253 0.245 0.147 0.405
#> tau2 0.104 0.096 0.054 0.203
#> tau2 [within] 0.036 0.034 0.017 0.064
#> tau2 [between] 0.068 0.060 0.022 0.164
#> I2 94.980 95.105 91.586 97.617
#> I2 [within] 35.298 34.226 13.779 63.148
#> I2 [between] 59.682 60.617 30.491 83.176
#> H2 22.097 20.430 11.886 41.970
plot() visualizes the posterior, and
prior = TRUE overlays the prior distribution so the shift
from prior to posterior is visible. Side-by-side panels for
tau and rho show the bulk of the posterior
moving away from the weakly informative defaults.
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.75)
plot(fit1_brma, parameter = "tau", prior = TRUE, main = "tau", xlim = c(0, 1))
plot(fit1_brma, parameter = "rho", prior = TRUE, main = "rho", xlim = c(0, 1))
The tau posterior concentrates around 0.3 with credible
interval comparable to the profile-likelihood bounds, while the
rho posterior favors a moderate-to-high cluster-level
allocation, broadly consistent with the REML point estimate of 0.67.
Other Inference Helpers
All inference helpers available for any other brma() fit
work the same way for the multilevel fit specified via
cluster. This includes meta-regression and location-scale
models, model comparison, and fit diagnostics. See the Bayesian
Meta-Analysis vignette for the full walkthrough; the only thing
that changes here is the addition of cluster = district.
For known sampling covariance matrices and richer random-effects
structures, see the Multivariate and
Multilevel Meta-Analysis vignette.