Skip to content

[R-meta] conversion between ranef() and blup() in rma() and rma.mv()

5 messages · Wolfgang Viechtbauer, Yefeng Yang

Yefeng Yang
# 
Dear MA community,

I am testing how to calculate the so-called empirical Bayes using metafor package. I found there are three ways of doing it. Theoretically, those ways should return the same values, but I found numerical differences. Please take a look at my illustration below.

I use dat.bcg embedded in metafor package as an example.

# load metafor
library(metafor)
# calculate effect size and sampling variance
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

The first way is to apply blup() to a rma object (fitted via rma()):
# use rma() to fit a random-effects model:
res <- rma(yi, vi, data=dat)

# use blup() to calculate empirical Bayes, which is the sum of the fitted/predicted values based on the fixed effects and the estimated random effects
blups <- as.data.frame(blup(res))

The second way is to sum up ranefs() and fitted():

# firstly, use ranef() to get the blups of the random effects
ranefs <- as.data.frame(ranef(res))

# secondly, sum up the the fitted fixed effect and the the blups of the random effects
blups$pred.fitted.ranefs <-  as.numeric(fitted(res)) + ranefs$pred # get the point estimate
blups$se.fitted.ranefs <- sqrt(res$se^2 + ranefs$se^2) # get the error

We see the point estimates (blups$pred vs.  blups$pred.fitted.ranefs) match well, while standard error (blups$se vs. blups$se.fitted.ranefs) does not.

The third way is to use rma.mv() to fit a RE model, and then sum up ranefs() and fitted():
# use rma.mv() to fit a random-effects model:
obs <- 1:nrow(dat)
res2 <- rma.mv(yi, vi, random = ~ 1 | obs, data = dat)

# firstly, use ranef() to get the blups of the random effects
ranefs2 <- as.data.frame(ranef(res2))

# secondly, sum up the the fitted fixed effect and the the blups of the random effects
# to compare with those derived from rma(), we add the values to the data frame blups
blups$pred.fitted.ranefs2 <- as.numeric(fitted(res2)) + ranefs2$obs.intrcpt # get the point estimate
blups$se.fitted.ranefs2 <- sqrt(res2$se^2 + ranefs2$obs.se^2) # get the error

While the values (both point estimate and error) from the three ways are very close, there are numerical differences. Anyone can comment on whether I am doing wrong?

Best,
Yefeng
Wolfgang Viechtbauer
# 
Dear Yefeng,

Please see my responses below.

Best,
Wolfgang
That's because sqrt(res$se^2 + ranefs$se^2) assumes independence between the fitted values and the BLUPs of the randome effects, which is not the case.
Nothing. The internal implementations ranef() is slightly different for rma.uni and rma.mv objects, which can lead to minor numerical differences in the SEs. But they are practically identical:
[1] 8.999273e-08 7.447359e-08 9.906757e-08 6.278849e-08 6.076396e-08 6.553534e-08 7.801192e-08
 [8] 6.639436e-08 6.081443e-08 6.143914e-08 6.418636e-08 1.089931e-07 6.135178e-08
Yefeng Yang
# 
Hi Wolfgang,

Brilliant. I appreciate your clarification. Would you like to point me the literature/formula underpinning the blup() (in rma()) and ranef() (in both rma() and rma.mv())?

I only know the following ref providing the formula for calculating study-specific effect. It is a kind of weighted mean of the overall effect and observed effect size.
Higgins J P T, Thompson S G, Spiegelhalter D J. A re-evaluation of random-effects meta-analysis[J]. Journal of the Royal Statistical Society Series A: Statistics in Society, 2009, 172(1): 137-159.

Best,
Yefeng
Wolfgang Viechtbauer
# 
I would have to start digging again, but these are useful references:

Raudenbush, S. W., & Bryk, A. S. (1985). Empirical Bayes meta-analysis. Journal of Educational Statistics, 10(2), 75-98. https://doi.org/10.3102/10769986010002075

van Aert, R. C. M., Schmid, C. H., Svensson, D., & Jackson, D. (2021). Study specific prediction intervals for random-effects meta-analysis: A tutorial. Research Synthesis Methods, 12(4), 429-447. https://doi.org/10.1002/jrsm.1490

One of the best references in my opinion is:

Searle, S. R., Casella, G., & McCullough, C. E. (1992). Variance components. Hoboken, NJ: Wiley.

It is not focused on meta-analysis, but the standard MA models are just mixed-effects models and this book provides one of the most thorough coverages of the underlying theory.

Best,
Wolfgang
Yefeng Yang
# 
Hi Wolfgang,

Thanks for spending time on my question and sending refs.

Best,
Yefeng