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

[Viechtbauer, Wolfgang (NP)]

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 <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-
> project.org>
> Cc: Michael Dewey <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    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.