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

[Viechtbauer, Wolfgang (NP)]

Dear Zac,

predict() currently does not provide predicition intervals for 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 <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Zac Robinson via R-sig-meta-analysis
> Sent: Sunday, March 10, 2024 16:26
> To: r-sig-meta-analysis using r-project.org
> Cc: Zac Robinson <zacrobinson2015 using gmail.com>
> Subject: [R-meta] Prediction Intervals for estimates of random slope 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=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=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=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