[R-meta] Prediction Intervals for estimates of random slope rma.mv model with predict.rma()
Viechtbauer, Wolfgang (NP)
wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Thu Mar 14 14:22:39 CET 2024
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 using gmail.com>
> Sent: Monday, March 11, 2024 16:07
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using 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 using 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 using 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 using r-project.org
> > Cc: Zac Robinson <mailto:zacrobinson2015 using 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
More information about the R-sig-meta-analysis
mailing list