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

Emanuel Schembri |e||d26 @end|ng |rom gm@||@com
Mon Sep 9 05:48:11 CEST 2024 

Dear Reza,

The orchaRd package has the orchard_plot() function for categorical
variables and the bubble_plot() function for numerical variables.

Regards,

Emanuel

On Sat, Sep 7, 2024 at 4:27 PM Reza Norouzian <rnorouzian using gmail.com> wrote:
>
> Thanks for sharing the github package, Wolfgang. It looks like orchard_plot--similar to emmprep()--can't handle rank deficiencies which commonly occur in meta-analytic datasets due to various missing variable combinations. It also seems that orchard_plot is only designed for categorical variables and can't be used with quantitative/continuous variables.
>
> I believe the issue that creates these limitations is that both emmprep() and orchard_plot() convert the rma objects to emmeans objects through emmeans::qdrg(). It might be possible to convert rma objects to one (I'm thinking lme::gls or something) that is readily acceptable in emmeans and likely both limitations may be lifted.
>
> Reza
>
>
> On Sat, Sep 7, 2024 at 8:15 AM Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>>
>> 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.

More information about the R-sig-meta-analysis mailing list