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

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