I am conducting a meta-analysis using the brms in R with structure akin to the following (leaving out priors, etc. for simplicity):

brm(yi|se(sei) ~ 1 + (1 | eid), data = dat)

Here yi is the effect, sei is the standard error of that effect and eid is a unique identifier for each experiment. My understanding is that inclusion of a random effect by experiment (eid) addresses concerns related to dependencies in my data (e.g., if one experiment reports two effects from the same sample they would share the same eid). I have published several papers using this approach (sometimes with additional random effects, moderators, etc.) in the past without issue. I have just received a review stating that the use of a multilevel Bayesian model is insufficient to account for such dependencies, and I must instead apply some form of robust variance estimator (of the sort recommended by Hedges, Tipton, & Johnson, 2010). However, I am unable to find anything in the Bayesian meta-analytic literature discussing this issue. The reviewer has recommended abandoning my Bayesian approach so that I can use the robust function from the metafor package instead. I am not keen to do that (as much as I love the metafor package).


I was hoping someone could point me towards an article or two that either (a) discusses robust variance estimation in the context of Bayesian meta-analysis, or, (b) explains why robust variance estimation is unnecessary in the context of Bayesian meta-analysis.


The reviewer is correct in the sense that this model does not properly account for correlated sampling errors (when multiple effects are computed based on the same sample of subjects, then the sampling errors of the estimates are typically correlated). By adding random effects to your model, you are modeling potential correlation/dependencies in the underlying true effects, but not the sampling errors. Just like you are specifying the standard errors of the estimates (sei), one would have to specify the covariances of the estimates (i.e., an entire variance-covariance matrix of the sampling errors) to properly account for the latter. This is an example where this is done:


In practice, constructing this var-cov matrix ($V$) is difficult. One can construct an approximate var-cov matrix based on a 'guestimate' of the correlation among the sampling errors and then treat this as the 'working model'. Even if the $V$ matrix is somewhat misspecified, the estimates of the fixed effects should be approximately unbiased. But the standard errors might be off. Hence, for inferences, we can use cluster robust standard errors.

This aside, just adding a random effect at the study level implicitly assumes that the true effects within studies are homogeneous. That is not an assumption I would recommend to make a priori.

Here is an example illustrating the different approaches (note: esid is an 'effect size id', not the study id).


dat <- dat.assink2016
dat$sei <- sqrt(dat$vi)
head(dat, 10)

#    study esid id      yi     vi pubstatus year deltype        sei
# 1      1    1  1  0.9066 0.0740         1  4.5 general 0.27202941
# 2      1    2  2  0.4295 0.0398         1  4.5 general 0.19949937
# 3      1    3  3  0.2679 0.0481         1  4.5 general 0.21931712
# 4      1    4  4  0.2078 0.0239         1  4.5 general 0.15459625
# 5      1    5  5  0.0526 0.0331         1  4.5 general 0.18193405
# 6      1    6  6 -0.0507 0.0886         1  4.5 general 0.29765752
# 7      2    1  7  0.5117 0.0115         1  1.5 general 0.10723805
# 8      2    2  8  0.4738 0.0076         1  1.5 general 0.08717798
# 9      2    3  9  0.3544 0.0065         1  1.5 general 0.08062258
# 10     3    1 10  2.2844 0.3325         1 -8.5 general 0.57662813

# fit multilevel model (assuming no covariance among the sampling errors)
res1 <- rma.mv(yi, vi, random = ~ 1 | study/esid, data=dat)

# analogous brm() model
res2 <- brm(yi|se(sei) ~ 1 + (1 | study/esid), data = dat)
summary(res2, priors=TRUE)

# construct approx V matrix assuming r=0.6 for the sampling errors
V <- impute_covariance_matrix(dat$vi, cluster=dat$study, r=0.6)

# fit multilevel model with the approximate V matrix
res3 <- rma.mv(yi, V, random = ~ 1 | study/esid, data=dat)

# use cluster robust inference methods
robust(res3, cluster=dat$study)

# or use coef_test() from clubSandwich
coef_test(res3, vcov="CR2", cluster=dat$study)

You will find that the results in res1 and res2 are quite similar.$^1$ The results in res3 differ more by having specified the approximate $V$ matrix. Using cluster robust methods barely changes the results in this case, but could do so more in other cases.

Some relevant papers on the use of robust variance estimation (also known as the 'sandwich method') in the context of Bayesian analyses are https://doi.org/10.1214/10-AOAS362 and https://doi.org/10.3982/ECTA9097 but these are not focused on meta-analysis and I am not aware of any papers that do so.

If you want to stick with brms, then I would recommend to at least specify the random effects as above. There is some evidence that ignoring the sampling error covariances can also work (the study-level variance component then gets inflated to compensate for essentially assuming that the covariances are equal to zero, which we also see happening in this example). See https://doi.org/10.3758/s13428-012-0261-6 and https://doi.org/10.3758/s13428-014-0527-2 and https://doi.org/10.1080/13645579.2016.1252189 (none of these discuss Bayesian approaches though).

$^1$ The mean of the posterior for the study-level standard deviation is a tad higher in res2 than the REML estimate in res1, presumably because with 17 studies, the prior (a half student-t with df=3 and scale parameter 2.5) pulls the mean of the posterior up a bit.

  • 1
    $\begingroup$ Wolfgang to the rescue, as usual. Thanks! $\endgroup$ – jmfawcet Jun 22 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.