Prediction intervals from predict_response(type="random") for meta-analysis using nlme

[zacrobinson5]

I am working on a dose-response multilevel meta-regression from which I'd like to extract the marginal effects with both confidence and prediction intervals. Traditionally, I've used the metafor and orchaRd packages for this, however, getting prediction intervals for random slope models has proven to be challenging.

I've since realized you can fit mixed-effect meta-analytic models using the nlme package as long as the weights and sigma arguments have been specified appropriately. Thus, I've been exploring options to move my post-fitting workflow to the numerous packages that support nlme models - which unfortunately metafor normally doesn't have access to.

With the preamble out of the way - I've now starting working with ggeffects to try to solve my issue. After reading the documentation, it seems the random effect variance for the prediction intervals is calculated using the work by Johnson et al., which is exactly what I was looking for; accounting for both random intercepts and random slopes.

That said, performing some preliminary tests with only random intercept multilevel meta-regressions, it seems like the prediction intervals are much larger than I'd expect compared to some of the other packages I'm used to. Please let me know if I am messing something up, or misunderstanding in any way. You can see my code below.

Thank you very much for the package and the thorough documentation - much appreciated!

#load packages
library(metafor)
library(orchaRd)
library(nlme)
library(tidyverse)
library(ggeffects)
library(emmeans)

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

#fit metafor
meta=rma.mv(
  data = df,
  yi ~ year + factor,
  random = list(~1|study/district/school),
  V = vi, 
  method = "REML"
)

#fit lme
mod=lme(
  data = df,
  yi ~ year + factor,
  random = ~1|study/district/school,
  weights = varFixed(~vi),
  control = lmeControl(sigma = 1),
  method = "REML"
)

#extract effects (metafor w/ orchaRd)
mod_results(meta,
            mod = "year",
            group = "study")%>%
  bubble_plot(group = "study")+
  ylim(-2.5,2.5)+
  theme_classic()+
  theme(legend.position = "none")+
  labs(title = "")

## prediction intervals look good

#extract effects (metafor w/ emmprep)
emmprep(meta,
        at=list(year=seq(1976,2000,length.out=100)))%>%
  emmeans(~year,
          weights = "prop",
          sigma=sqrt(sum(meta$sigma2)))%>%
  emmip(~year,CIs = TRUE,PIs = TRUE)+
  ylim(-2.5,2.5)+
  theme_classic()+
  theme(legend.position = "none")+
  labs(title = "")+
  geom_point(data = df,
             aes(x=year,
                 y=yi,
                 size=1/vi))

## prediction intervals same as orchaRd

#extract effects (lme)
eff.1=predict_response(mod,
                 terms = "year[1976:2000]",
                 margin = "marginalmeans",
                 weights = "prop",
                 type = "fixed")
eff.2=predict_response(mod,
                       terms = "year[1976:2000]",
                       margin = "marginalmeans",
                       weights = "prop",
                       type = "random")%>%
  rename(pred.low=conf.low,
         pred.high=conf.high)

eff=cbind(eff.1,eff.2[,c("pred.low","pred.high")])

ggplot(data = eff,
       aes(x=x,
           y=predicted))+
  geom_ribbon(aes(ymin=pred.low,
                  ymax=pred.high),
              linetype="dashed",
              color="black",
              alpha=0)+
  geom_ribbon(aes(ymin=conf.low,
                  ymax=conf.high),
              alpha=0.25)+
  geom_line(linewidth=1)+
  ylim(-2.5,2.5)+
  theme_classic()+
  theme(legend.position = "none")+
  labs(title = "")+
  geom_point(data = df,
             aes(x=year,
                 y=yi,
                 size=1/vi))

## prediction intervals much wider