[R-meta] Metafor - regplot() for a categorical (and interaction) variable

[Viechtbauer, Wolfgang (NP)]

Interesting alternative. Speaking of alternatives, the orchaRd package also provides some nice functionality along those lines:

library(metafor)

remotes::install_github("daniel1noble/orchaRd")
library(orchaRd)

dat <- escalc(measure="RR", ai=tpos, bi=tneg,
              ci=cpos, di=cneg, data=dat.bcg)

res <- rma(yi, vi, mods = ~ alloc, data=dat)

orchard_plot(res, mod="alloc", group="trial", xlab="log(Risk Ratio)")

Best,
Wolfgang

> -----Original Message-----
> From: Reza Norouzian <rnorouzian using gmail.com>
> Sent: Saturday, September 7, 2024 14:40
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>; lelid26 using gmail.com
> Cc: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
> Subject: Re: [R-meta] Metafor - regplot() for a categorical (and interaction)
> variable
>
> Emanuel,
>
> Wolfgang's response is spot on. A shorter alternative might be to try the
> following.
>
> Reza
>
> dat <- escalc(measure="RR", ai=tpos, bi=tneg,
>               ci=cpos, di=cneg, data=dat.bcg)
>
> res <- rma(yi, vi, mods = ~ alloc, data=dat)
>
> grd <- emmprep(res) # check out metafor's documentation
>
> emmeans::emmip(grd, ~alloc, CIs=TRUE)  # check out emmeans' documentation
>
> On Sat, Sep 7, 2024 at 7:04 AM Viechtbauer, Wolfgang (NP) via R-sig-meta-
> analysis <mailto:r-sig-meta-analysis using r-project.org> wrote:
> This is one issue. However, there are more issues at hand here, since
> regplot(res, mod = ~ Continent) isn't gong to work either, as this will trigger
> the error "Can only specify a single variable via argument 'mod'." The reason is
> that the 'mod' argument is not meant to take a formula as input. However, when
> you do, then the length of 'mod' will be 2 (length(~ blah)), which triggers the
> error.
>
> @Emanuel: While I don't know what you really intend to do here, I suspect you
> would like to show all levels of 'Continent' on the x-axis. This is not so
> easily accomplished with the regplot() function, which is designed for
> visualizing quantitative/continuous moderators. For example:
>
> dat <- escalc(measure="RR", ai=tpos, bi=tneg,
>                             ci=cpos, di=cneg, data=dat.bcg)
>
> res <- rma(yi, vi, mods = ~ ablat, data=dat)
> res
>
> regplot(res, mod="ablat")
>
> A dichotomous moderator is also easily handled:
>
> dat$random <- ifelse(dat$alloc == "random", 1, 0)
>
> res <- rma(yi, vi, mods = ~ random, data=dat)
> res
>
> regplot(res, mod="random")
>
> We might want to make this look a bit nicer, maybe like this:
>
> regplot(res, mod="random", xlab="Method of Treatment Allocation", xaxt="n")
> axis(side=1, at=c(0,1), labels=c("Non-Random", "Random"))
>
> However, things become more tricky with a factor variable that has 3 or more
> levels:
>
> res <- rma(yi, vi, mods = ~ alloc, data=dat)
> res
>
> regplot(res, mod="alloc")
>
> This will not work, since there are two 'alloc' dummy variables, but regplot()
> is designed to place a single variable on the x-axis. One could do:
>
> par(mfrow=c(2,1))
> regplot(res, mod="allocrandom")
> regplot(res, mod="allocsystematic")
> par(mfrow=c(1,1))
>
> but this is showing the difference between random and not random (consisting of
> 'systematic' and 'alternate') in the first plot and the difference between
> systematic and not systematic (consisting of 'random' and 'alternate') in the
> second plot. Probably not how most people would want to visualize this. Instead,
> I suspect most would want to show three 'columns' of points, corresponding to
> the three levels, with lines connecting the fitted/predicted values for these
> levels.
>
> I had not quite considered that some may want to do something like this with
> this function. I am not sure how easy it would be to add this kind of
> functionality directly to regplot() given some of the internal intricacies.
> However, with a small trick, we can still accomplish this. A model with a
> categorical moderator with p levels can be represented as a polynomial
> regression model to the degree p-1 where the first term is the linear one (which
> we want to place on the x-axis). This, combined with the possibility to pass
> predicted values to regplot() via the 'pred' argument, we can place all levels
> on the x-axis as follows:
>
> dat$anum <- as.numeric(factor(dat$alloc))
> res <- rma(yi, vi, mods = ~ poly(anum, degree=2, raw=TRUE), data=dat)
> res
>
> pred <- predict(res, newmods=unname(poly(1:3, degree=2, raw=TRUE)))
> pred
>
> regplot(res, mod=2, pred=pred, xvals=c(1:3), xlim=c(1,3), xlab="Allocation
> Method", xaxt="n")
> axis(side=1, at=1:3, labels=levels(factor(dat$alloc)))
>
> Best,
> Wolfgang
>
> > -----Original Message-----
> > From: R-sig-meta-analysis <mailto:r-sig-meta-analysis-bounces using r-project.org>
> On Behalf
> > Of Michael Dewey via R-sig-meta-analysis
> > Sent: Friday, September 6, 2024 11:04
> > To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> > http://project.org>
> > Cc: Michael Dewey <mailto:lists using dewey.myzen.co.uk>
> > Subject: Re: [R-meta] Metafor - regplot() for a categorical (and interaction)
> > variable
> >
> > Dear Emanuel
> >
> > I think you will find that the parameter is named mods in rma.uni but
> > mod in regplot.
> >
> > Michael
> >
> > On 06/09/2024 00:56, Emanuel Schembri via R-sig-meta-analysis wrote:
> > > Hi
> > >
> > > I am trying to plot a bubble plot using the replot() function in
> > > Metafor. However, I cannot make it work for a categorical moderator
> > > (and an interaction between a categorical and numerical variable),
> > > while it does work when inputting a numerical or date.
> > >
> > > Any help would be greatly appreciated.
> > >
> > > Regards,
> > >
> > > Emanuel
> > >
> > > res <- rma(yi, vi, mods = ~ Continent, data= df)
> > > regplot(res, mods = ~ Continent)
> > >
> > > Mixed-Effects Model (k = 24; tau^2 estimator: REML)
> > >
> > > tau^2 (estimated amount of residual heterogeneity):     3.5128 (SE = 1.1529)
> > > tau (square root of estimated tau^2 value):             1.8743
> > > I^2 (residual heterogeneity / unaccounted variability): 99.79%
> > > H^2 (unaccounted variability / sampling variability):   477.76
> > > R^2 (amount of heterogeneity accounted for):            8.05%
> > >
> > > Test for Residual Heterogeneity:
> > > QE(df = 21) = 1150.4471, p-val < .0001
> > >
> > > Test of Moderators (coefficients 2:3):
> > > QM(df = 2) = 5.7609, p-val = 0.0561
> > >
> > > Model Results:
> > >
> > >                     estimate      se    zval    pval    http://ci.lb   ci.ub
> > > intrcpt              2.2738  1.3346  1.7037  0.0884  -0.3420  4.8897  .
> > > ContinentAsia      4.3152  1.8165  2.3756  0.0175   0.7549  7.8755  *
> > > ContinentEurope    2.6747  1.4053  1.9032  0.0570  -0.0797  5.4291  .
> > >
> > > ---
> > > Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> > >
> > > regplot(res, mods = ~ Continent)
> > >
> > > Error in regplot.rma(res, mods = ~Continent) :
> > >
> > >    Must specify 'mod' argument for models with multiple predictors.