Presenting results of a meta-analysis with multiple moderators?

Question

I wish to present the results of my meta-analysis using the best practices possible. I do not find, however, examples in articles similar to what my output is. Here's a simplification of my model and output (using function rma of the Metafor package, R).

rma (yi, vi, mods = ~ varA + varB)

For those unfamiliar with rma, this line code a meta-analysis mixed effect models, with 2 moderators, the variable A and variable B. Let's say the varA has 2 categorical levels and varB has 3 categorical levels.

The output would be similar to this, without any numbers:

intercept - estimate, statistic, p-value, 95% IC
varA.level2 - estimate, statistic, p-value, 95% IC
varB.level2 - estimate, statistic, p-value, 95% IC
varB.level3 - estimate, statistic, p-value, 95% IC

All those variable levels are significantly different from the intercept. I have no problem with varA. I can easily say that level 1 is different from level 2 and report the statistic and p-value for varA.level2.

But what about varB? level 1 differs from level 2, level 1 from level 3... but is level 2 different from 3? What would you report in a article in that case?

Thanks for insight or references to articles, books...

Answer 1

In essence, there are 6 different combinations of those two moderators. So, you could just compute and present the estimated/predicted effects for each of those combinations. Here is an example using a dataset from the metafor package that also happens to include two categorical moderators with 2 and 3 levels, respectively:

### load data
dat <- get(data(dat.mcdaniel1994))

### calculate r-to-z transformed correlations and corresponding sampling variances
dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat)

### fit mixed-effects meta-regression model (struct has 2 levels, type as 3 levels)
res <- rma(yi, vi, mods = ~ struct + type, data=dat)

### compute the estimated/predicted correlation for each combination
predict(res, newmods=rbind(c(0,0,0), c(0,1,0), c(0,0,1),
                           c(1,0,0), c(1,1,0), c(1,0,1)), transf=transf.ztor, addx=TRUE, digits=2)

The results are:

  pred ci.lb ci.ub cr.lb cr.ub X.intrcpt X.structu X.typep X.types
1 0.26  0.22  0.30 -0.08  0.54         1         0       0       0
2 0.19  0.01  0.36 -0.19  0.52         1         0       1       0
3 0.30  0.19  0.40 -0.05  0.58         1         0       0       1
4 0.20  0.13  0.26 -0.15  0.50         1         1       0       0
5 0.13 -0.04  0.29 -0.24  0.47         1         1       1       0
6 0.24  0.10  0.36 -0.13  0.54         1         1       0       1

One could also put these results into a forest plot (not showing the individual estimates, just these estimated outcomes):

slabs <- c("struct = s, type = j", "struct = s, type = p", "struct = s, type = s", "struct = u, type = j", "struct = u, type = p", "struct = u, type = s")
par(mar=c(4,4,1,2))
forest(sav$pred, sei=sav$se, slab=slabs, transf=transf.ztor, xlab="Correlation", xlim=c(-.4,.7))
text(-.4, 8, "Structure/Type",       pos=4, font=2)
text(.7,  8, "Correlation [95% CI]", pos=2, font=2)

That would look something like this:

There is a bit of redundancy in presenting all 6 combinations, since your model assumes that the influence of the two moderators is additive. But I think this makes the results quite clear.