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

[Viechtbauer, Wolfgang (NP)]

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

> -----Original Message-----
> From: Yefeng Yang <yefeng.yang1 using unsw.edu.au>
> Sent: Monday, March 4, 2024 13:23
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>; R
> Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-project.org>
> Subject: Re: conversion between ranef() and blup() in rma() and rma.mv()
>
> 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
> ________________________________________
> From: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
> Sent: 04 March 2024 20:07
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> Cc: Yefeng Yang <yefeng.yang1 using 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 using 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 using r-project.org
> > Cc: Yefeng Yang <yefeng.yang1 using 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