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

________________________________________
From: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer at maastrichtuniversity.nl>
Sent: 04 March 2024 20:07
To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r-
project.org>
Cc: Yefeng Yang <yefeng.yang1 at unsw.edu.au>
Subject: RE: conversion between ranef() and blup() in rma() and rma.mv()

Dear Yefeng,

Please see my responses below.

Best,
Wolfgang

-----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.

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.

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?

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:

blups$se.fitted.ranefs2 - blups$se.fitted.ranefs
 [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

Best,
Yefeng