[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