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

[Zac Robinson]

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

	[[alternative HTML version deleted]]