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

[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

	[[alternative HTML version deleted]]