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
[R-meta] Partial dependence plots
4 messages · Cesar Terrer Moreno, Wolfgang Viechtbauer
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 at r-project.org] On Behalf Of Cesar Terrer Moreno Sent: Tuesday, 19 June, 2018 14:52 To: r-sig-meta-analysis at 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 _______________________________________________ R-sig-meta-analysis mailing list R-sig-meta-analysis at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
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 at 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 at r-project.org] On Behalf Of Cesar Terrer Moreno Sent: Tuesday, 19 June, 2018 14:52 To: r-sig-meta-analysis at 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
_______________________________________________ R-sig-meta-analysis mailing list R-sig-meta-analysis at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
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 at me.com] Sent: Tuesday, 19 June, 2018 19:38 To: Viechtbauer, Wolfgang (SP) Cc: r-sig-meta-analysis at 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 at 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 at r-project.org] On Behalf Of Cesar Terrer Moreno Sent: Tuesday, 19 June, 2018 14:52 To: r-sig-meta-analysis at 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