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)
# 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)
# 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
# or use coef_test() from clubSandwich
coef_test(res3, vcov="CR2", cluster=dat$study)
You will find that the results in
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.