Skip to content

[R-meta] Partial residual plots and meta-analysis

4 messages · Viechtbauer Wolfgang (STAT), Cesar Terrer Moreno, Michael Dewey

Original
Cesar Terrer Moreno
# 
Dear All,

Using our beloved metafor and model selection, I got to the conclusion that my best model was of the form:

mymodel <- rma(es, var, data=dat , mods= ~ 1 + A + B*C)

I want to show in a figure (to be included in a paper) the relationship between es and the individual moderators, including the interaction, while taking into account the effect of the rest of the moderators. Although I don?t have experience in this matter, it seems that partial residual plots is what I am looking for. I have tried to use the package visreg for this purpose, but apparently it cannot handle a rma output:

visreg(mymodel)
Error in formula.default(fit) : invalid formula

Do you know how to visualise the individual effects of the moderators in a meta-regression? 
Thanks
C?sar
Viechtbauer Wolfgang (STAT)
# 
Hi C?sar,

To show how the expected value of 'es' changes as a function of a moderator, just use the predict() function, computing predicted values for varying values of the moderator. An example for the simplest case (with a single moderator) is shown here:

http://www.metafor-project.org/doku.php/plots:meta_analytic_scatterplot

In models with multiple moderators, you can do the same thing. "Taking into account the effect of the rest of the moderators" simply means that you hold the other moderators constant. A common practice is to hold continuous moderators constant at their mean (although one can choose any other sensible value). Using the same example but with two moderators:

library(metafor)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
res <- rma(yi, vi, mods = ~ ablat + year, data=dat)
predict(res, newmods=cbind(13:55, mean(dat$year)))      # holding year constant at the mean
predict(res, newmods=cbind(mean(dat$ablat), 1948:1980)) # holding ablat constant at the mean

For models with interactions, you do the same thing. For example:

res <- rma(yi, vi, mods = ~ ablat + year + ablat:year, data=dat)
predict(res, newmods=cbind(13:55, mean(dat$year), 13:55 * mean(dat$year)))
predict(res, newmods=cbind(mean(dat$ablat), 1948:1980, mean(dat$ablat) * 1948:1980))

Since the slope of one moderator changes as a function of the value of the other moderator when the two moderators interact, one may also want to obtain predicted values when holding the other moderator constant at a couple different values. For example:

predict(res, newmods=cbind(13:55, 1948, 13:55 * 1948)) # holding year constant at lower bound
predict(res, newmods=cbind(13:55, 1966, 13:55 * 1966)) # holding year constant at the mean (rounded)
predict(res, newmods=cbind(13:55, 1980, 13:55 * 1980)) # holding year constant at upper bound

One can then plot the three lines to show how the expected value of 'es' changes as a function of ablat when year is held constant at its lower bound, mean, and upper bound.

Best,
Wolfgang
Cesar Terrer Moreno
# 
Hi Wolfgang,

Many thanks, I like this approach too.

One minor question: when I try to plot the result in ggplot, it complains about the class of the output from pred:

MAPpred <- predict(ECM, newmods=cbind(645:1750, mean(ecmdat$MAT), 300,  mean(ecmdat$MAT) * 300),  addx = T)

ggplot(MAPpred, aes(X[,"MAP"], make_pct(pred))) + ?

Error: ggplot2 doesn't know how to deal with data of class list.rma

I tried to do ?as.data.frame?, but it doesn?t work. Any suggestion?
Thanks
C?sar
Michael Dewey
# 
Dear Cesar

You could try data.frame(unlist(MAPpred)) although I am not absolutely 
sure it gives you what ggplot2 expects (I am not a ggplot2 user)

Michael
On 15/07/2017 08:46, Cesar Terrer Moreno wrote: