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

-----Original Message-----
From: R-sig-meta-analysis <r-sig-meta-analysis-bounces at r-project.org> On Behalf
Of Yefeng Yang via R-sig-meta-analysis
Sent: Monday, March 4, 2024 00:41
To: r-sig-meta-analysis at r-project.org
Cc: Yefeng Yang <yefeng.yang1 at unsw.edu.au>
Subject: [R-meta] conversion between ranef() and blup() in rma() and rma.mv()

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?

blups$se.fitted.ranefs2 - blups$se.fitted.ranefs

Best,
Yefeng