[R-meta] Plotting rma.mv as forest plots for individual moderator levels

[Viechtbauer, Wolfgang (NP)]

Dear Julian,

Please see below for my responses.

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On
>Behalf Of Julian Packheiser via R-sig-meta-analysis
>Sent: Saturday, 04 February, 2023 13:58
>To: r-sig-meta-analysis using r-project.org
>Cc: Julian Packheiser
>Subject: [R-meta] Plotting rma.mv as forest plots for individual moderator levels
>
>Dear members of the mailings list,
>
>Currently, we are conducting a multivariate multilevel meta-analysis with
>multiple outcomes coming from the same study. I use both study_id and
>effectsize_id for the calculation of V as well as as random effects in the model.
>We also added the laboratory in which the study was conducted as a random effect.
>The code thus looks like this:
>
>V <- vcalc(vi, cluster=study_id, obs=esid, data=mydata, rho=0.6)
>
># random-effects model
>res <- rma.mv(yi, V, mods = ~ outcome-1, random = ~ 1 | lab/study_id/esid,
>data=mydata)
>
>My first question is, should we use a sensitivity approach for rho values as they
>are unknown?

That's always a good idea when 'guestimating' a value. You could also consider using cluster-robust inference methods as a follow-up.

>My second question is about forest plots. Ideally, I would like to plot the
>results for each moderator separately. If I use forest.rma, I get all results in
>one plot. An overall result for this meta-analysis is however not meaningful
>since some outcomes are positively correlated and some are negatively correlated
>with the variable of interest.
>I assume that I cannot simply do a forest plot of a rma.mv result for each
>outcome individually as the results obviously differ from a multivariate
>approach. Could I simply re-arrange my outcomes and only plot them ordered by
>outcome? I assume the plotted values are yi.f and vi.f for the means and
>variances, but I could not find values for the fitted polygons.

With 'moderators', I assume you mean 'outcome' (as in the model above).

I would use forest() (forest.default to be precise) to draw a forest plot per outcome and then use addpoly() to add the estimated pooled effect from the model above to each plot. For example:

library(metafor)

dat <- dat.berkey1998
V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat)
res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial, struct="UN", data=dat)
res

with(dat, forest(yi, vi, subset=outcome=="AL", slab=paste0(author, year), header=TRUE, ylim=c(-1.5,8)))
addpoly(predict(res, newmods=c(1,0)), row=-1)

with(dat, forest(yi, vi, subset=outcome=="PD", slab=paste0(author, year), header=TRUE, ylim=c(-1.5,8)))
addpoly(predict(res, newmods=c(0,1)), row=-1)

I would recommend to add a note to these plots that the pooled estimate is based on the multivariate/multilevel model (as indeed pooling the estimates within each plot would lead to a slightly different result).

>Apologies if this question has already been asked, but I could not find it in the
>archives.
>
>Best
>Julian
>
>---
>Julian Packheiser
>Postdoc in the Social Brain Lab
>Netherlands Institute for Neuroscience