[R-meta] Prediction Intervals for estimates of random slope rma.mv model with predict.rma()

Dear Zac,

Please respond to the mailing list, not just to me personally.

Just summing up all the variances isn't correct - the correct computation
of the prediction interval for such a random-effects structure is
unfortunately more complex.

Best,
Wolfgang

-----Original Message-----
From: Zac Robinson <zacrobinson2015 at gmail.com>
Sent: Monday, March 11, 2024 16:07
To: Viechtbauer, Wolfgang (NP) <

wolfgang.viechtbauer at maastrichtuniversity.nl>

Subject: Re: [R-meta] Prediction Intervals for estimates of random slope

rma.mv

model with predict.rma()

Hey Wolfgang,

Thank you very much for the reply. I actually have played around with

things a

bit more and think I've found a solution. I may not be specifying the

sigma

argument perfectly for a random slope model - but I think in concept

this works!

Looks like emmeans / predict() have cross compatibility that allows for
prediction intervals when sigma is specified

Let me know what you think,
Zac

# SETUP

-------------------------------------------------------------------

#load packages and data
library(metafor)
library(tidyverse)
library(emmeans)

data("dat.konstantopoulos2011")
df=dat.konstantopoulos2011%>%
  as.data.frame()%>%
  mutate(factor=sample(c("A","B"),56,replace=TRUE))

# RANDOM SLOPE (x1)

-------------------------------------------------------

#fit single random slope model
m2= rma.mv(yi ~ year+factor, vi,
          random = list(~year|study, ~1|district, ~1|school),
          data = df, method = "REML",struct = "GEN",
          control = list(optimizer="bobyqa"))

#extract estimates with emmmprep + emmeans + predict
est.ci=emmprep(m2,
        at=list(year=1976:2000))%>%
  emmeans(~year,
          weights = "prop",
          sigma=sqrt(sum(
            m2$sigma2,
            m2$tau2,
            m2$gamma2
          )))

pi= est.ci%>%
  predict(interval="prediction")%>%
  select(lower.PL,upper.PL)

preds=cbind(est.ci,pi)
preds # prediction intervals are now added

ggplot(data = preds,
       aes(x=year,
           y=emmean))+theme_classic()+
  geom_point(data=df,
             aes(y=yi,
                 size=1/vi))+
  geom_ribbon(aes(ymin=lower.PL,
                  ymax=upper.PL),
              linetype="dashed",
              alpha=0,
              color="black")+
  geom_ribbon(aes(ymin=asymp.LCL,
                  ymax=asymp.UCL),
              alpha=0.25)+
  geom_line()

On Mon, Mar 11, 2024 at 2:06?AM Viechtbauer, Wolfgang (NP)
<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
Dear Zac,

predict() currently does not provide predicition intervals for

http://rma.mv

models with a GEN random effects structure. Not sure when I will get

around to

implementing that.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis <mailto:

r-sig-meta-analysis-bounces at r-project.org>

On Behalf

Of Zac Robinson via R-sig-meta-analysis
Sent: Sunday, March 10, 2024 16:26
To: mailto:r-sig-meta-analysis at r-project.org
Cc: Zac Robinson <mailto:zacrobinson2015 at gmail.com>
Subject: [R-meta] Prediction Intervals for estimates of random slope

http://rma.mv

model with predict.rma()

Dear All,

I am working on extracting the adjusted effects from a dose-response
multilevel meta-regression. One issue I've run into is getting the
appropriate prediction intervals on my estimates when the model
includes random slopes (specifically when struct = "GEN" is included).

From

what I gather from the help files, this comes down to appropriately
specifying the `tau2.levels` and `gamma2.levels` arguments within
predict.rma() - however, I am unsure how to get this to work.

In the code below I've tried to demonstrate my issue and explore the

other

solutions I'm aware of through the orchaRd package). My end goal

requires

adjusting for both continuous and categorical moderators so ideally the
solution stays within the emmprep() / emmeans() framework (which makes
proportional weights and setting non-focal continuous variables at the

mean

easy), but I realize this may not be possible.

Thank you very much in advance!

Zac

# SETUP

-------------------------------------------------------------------

#load packages and data
library(metafor)
library(tidyverse)
library(orchaRd)
library(emmeans)

data("dat.konstantopoulos2011")
df=dat.konstantopoulos2011%>%
  as.data.frame()%>%
  mutate(factor=sample(c("A","B"),nrow(df),replace=TRUE))

# RANDOM INTERCEPT

--------------------------------------------------------

#fit random intercept model
m1=http://rma.mv(yi ~ year+factor, vi,
          random = list(~1|study, ~1|district, ~1|school),
          data = df, method = "REML")

#extract estimates with predict.rma
predict(m1,
        newmods =cbind(rep(c(1976:2000),2),
                       rep(c(0,1),each=25)),
        addx = TRUE) # works fine

#extract estimates with emmprep
emmprep(m1,
        at=list(year=1976:2000))%>%
  emmeans(~year,
          weights = "prop")%>%
  pred_interval_esmeans(model=m1,
                        mm=.,
                        mod="year") # works fine

#extract estimates with orchaRd
mod_results(m1,
            mod = "year",
            group = "study") # works fine

# RANDOM SLOPE (x1)

-------------------------------------------------------

#fit single random slope model
m2=http://rma.mv(yi ~ year+factor, vi,
          random = list(~year|study, ~1|district, ~1|school),
          data = df, method = "REML",struct = "GEN",
          control = list(optimizer="bobyqa"))

#extract estimates with predict.rma
predict(m2, newmods =cbind(rep(c(1976:2000),2),
                           rep(c(0,1),each=25)),
        addx = TRUE) # no prediction intervals

predict(m2, newmods =cbind(rep(c(1976:2000),2),
                           rep(c(0,1),each=25)),
        addx = TRUE,
        tau2.levels = levels(year), # try adding tau2.levels
        ) # still no prediction intervals

#extract estimates with emmprep
emmprep(m2,
        at=list(year=1976:2000))%>%
  emmeans(~year,
          weights = "prop")%>%
  pred_interval_esmeans(model=m2,
                        mm=.,
                        mod="year") # works but gives an error

#extract estimates with orchaRd
mod_results(m2,
            mod = "year",
            group = "study") # works but gives an error


# RANDOM SLOPE (x2)

-------------------------------------------------------

#fit double random slope model
m3=http://rma.mv(yi ~ year+factor, vi,
          random = list(~year|study, ~year|district, ~1|school),
          data = df, method = "REML",struct = c("GEN","GEN"),
          control = list(optimizer="bobyqa"))

#extract estimates with predict.rma
predict(m3, newmods =cbind(rep(c(1976:2000),2),
                           rep(c(0,1),each=25)),
        addx = TRUE) # no prediction intervals

predict(m3, newmods =cbind(rep(c(1976:2000),2),
                            rep(c(0,1),each=25)),
        addx = TRUE,
        tau2.levels = levels(year),# try adding tau2.levels
        gamma2.levels = levels(year)# try adding gamma2.levels
        ) # still no prediction intervals

#extract estimates with emmprep
emmprep(m3,
             at=list(year=1976:2000))%>%
  emmeans(~year,
          weights = "prop")%>%
  pred_interval_esmeans(model=m3,
                        mm=.,
                        mod="year") # works but gives an error

#extract estimates with orchaRd
mod_results(m3,
            mod = "year",
            group = "study") # works but gives an error