[R-meta] 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