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

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 at gmail.com>
Sent: Saturday, September 7, 2024 14:40
To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at r-
project.org>; lelid26 at gmail.com
Cc: Viechtbauer, Wolfgang (NP) <

wolfgang.viechtbauer at 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 at 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 at 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 at r-
http://project.org>
Cc: Michael Dewey <mailto:lists at 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.