[R-meta] Partial dependence plots

[Viechtbauer, Wolfgang (SP)]

Here is an example with a three-level factor. It's easier to just plug in the means for predictors we are averaging over.

library(metafor)

dat <- get(data(dat.bcg))
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
res <- rma(yi, vi, mods = ~ ablat + alloc, data=dat)
res

xs <- seq(min(dat$ablat), max(dat$ablat), by=1)
pred <- predict(res, newmods = cbind(xs, mean(res$X[,3]), mean(res$X[,4])))$pred
plot(xs, pred, type="o")

pred1 <- predict(res, newmods = cbind(mean(dat$ablat), 0, 0))$pred
pred2 <- predict(res, newmods = cbind(mean(dat$ablat), 1, 0))$pred
pred3 <- predict(res, newmods = cbind(mean(dat$ablat), 0, 1))$pred
plot(1:3, c(pred1, pred2, pred3), type="o", xaxt="n")
axis(side=1, at=1:3, labels=levels(factor(dat$alloc)))

There are more elegant ways of doing this, but I hope the idea is clear.

Best,
Wolfgang

-----Original Message-----
From: Cesar Terrer Moreno [mailto:cesar.terrer using me.com] 
Sent: Tuesday, 19 June, 2018 19:38
To: Viechtbauer, Wolfgang (SP)
Cc: r-sig-meta-analysis using r-project.org
Subject: Re: [R-meta] Partial dependence plots

Hi Wolfgang,

Many thanks for your quick response. 

How should I modify your code if e.g. C is a factor?

Best wishes,
César

> On 19 Jun 2018, at 17:15, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> 
> Hi César,
> 
> I do not know of an automated way of doing this. I looked at the 'pdp' package, but it might take some effort to make it work together with metafor. However, doing this manually shouldn't be too complicated. Here is an example using a relatively simple model with two predictors:
> 
> library(metafor)
> 
> dat <- get(data(dat.bcg))
> dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
> res <- rma(yi, vi, mods = ~ ablat + year, data=dat)
> res
> 
> xs <- seq(min(dat$ablat), max(dat$ablat), by=1)
> pred <- sapply(xs, function(x) mean(predict(res, newmods = cbind(x, dat$year))$pred))
> plot(xs, pred, type="o")
> 
> ### same as:
> pred <- predict(res, newmods = cbind(xs, mean(dat$year)))$pred
> lines(xs, pred, col="red")
> 
> xs <- seq(min(dat$year), max(dat$year), by=1)
> pred <- sapply(xs, function(x) mean(predict(res, newmods = cbind(dat$ablat, x))$pred))
> plot(xs, pred, type="o")
> 
> ### same as:
> pred <- predict(res, newmods = cbind(mean(dat$ablat), xs))$pred
> lines(xs, pred, col="red")
> 
> In your case:
> 
> M <- rma(yi, vi, dat, mods= ~ A*B + C + D)
> 
> xs <- seq(min(dat$C), max(dat$C), by=1) # or use an appropriate 'by' value
> pred <- sapply(xs, function(x) mean(predict(res, newmods = cbind(dat$A, dat$B, dat$A*dat$B, x, dat$D))$pred))
> plot(xs, pred, type="o")
> 
> xs <- seq(min(dat$D), max(dat$D), by=1) # or use an appropriate 'by' value
> pred <- sapply(xs, function(x) mean(predict(res, newmods = cbind(dat$A, dat$B, dat$A*dat$B, dat$C, x))$pred))
> plot(xs, pred, type="o")
> 
> Best,
> Wolfgang
> 
> -----Original Message-----
> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Cesar Terrer Moreno
> Sent: Tuesday, 19 June, 2018 14:52
> To: r-sig-meta-analysis using r-project.org
> Subject: [R-meta] Partial dependence plots
> 
> Dear all,
> 
> I have a mixed-effects model of the form:
> 
> M <- rma(yi, vi, dat, mods= ~ A*B + C + D)
> 
> How can I produce partial dependence plots of e.g. C and D (predicted effect sizes of C and D as a function of the value of each predictor variable). The idea is to show that once the interaction A*B is taken into account, C and D explain very little of the overall effect size.
> 
> Can this be coded into a function to make partial dependence plots for all variables?
> 
> Thanks
> César