Producing a regplot in metafor for multivariable model with a spline (marginal relation)

Question

I was trying to produce a regression plot with the regplot() in the metafor package, fitting a continuous predictor as a restricted cubic spline through the use of the rms package in R, in the context of a multivariable regression.

However, after carefully reviewing all the examples listed here: https://wviechtb.github.io/metafor/reference/regplot.html, I was unable to reproduce the plot desired (i.e., the regression plot of the continuous predictor modelled as a restricted cubic spline).

Here is a reproducible example:

library(metafor)
library(rms)
dat <- dat.bcg

#Fitting the meta-analysis
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat)

#Fit a multivariable meta-regression with ablat as a rcs with 3 knots
res <- rma(yi, vi, mods = ~ rcs(ablat,3)+year+alloc, data=dat)

The result of the model is reported below:

Model Results:

                     estimate       se     zval    pval      ci.lb    ci.ub   ​ 
intrcpt              -24.8731  38.8083  -0.6409  0.5216  -100.9361  51.1898    
rcs(ablat, 3)ablat     0.0043   0.0290   0.1497  0.8810    -0.0525   0.0611    
rcs(ablat, 3)ablat'   -0.0504   0.0466  -1.0799  0.2802    -0.1418   0.0410    
year                   0.0125   0.0196   0.6384  0.5232    -0.0259   0.0509    
allocrandom           -0.3525   0.4090  -0.8618  0.3888    -1.1541   0.4491    
allocsystematic       -0.2187   0.4832  -0.4526  0.6508    -1.1659   0.7284  

Now if I want to produce a regplot of this model (i.e., the marginal relationship between the spline predictor and the effect size), I must use:

regplot(res, mod="rcs(ablat, 3)ablat")

But the resulting plot (reported below) seems to me taking into account only the first term of the spline-modelled variable:

Consistently, I am able to produce another plot if I change my call of regplot to regplot(res, mod="rcs(ablat, 3)ablat'") (i.e., specifying as the moderator variable the non-linear term of ablat in the multivariable regression model).

I have tried to pre-specify prediction through the predict() function and then passing the argument to regplot, but I am still unable to produce the desired plot.

Answer 1

I am posting here the solution, that I have managed to arrange after the inestimable help of @Wolfgang in the comments, just in case someone else needs in the future:

res <- rma(yi, vi, mods = ~ rcs(ablat,3)+year+alloc, data=dat)

knots <- attr(rcs(model.matrix(res)[,2], 3), "parms")
knots
xs <- seq(10,60, length=500)
sav <- predict(res, newmods=cbind(rcspline.eval(xs, knots, inclx=TRUE),colMeans(model.matrix(res))[4],
                                  colMeans(model.matrix(res))[5],colMeans(model.matrix(res))[6]))
tmp <- regplot(res, mod=2, pred=sav, xvals=xs, las=1, digits=1, bty="l",
               psize=.20/sqrt(dat$vi), xlab="Predictor", main="Restricted Cubic Spline Model", transf=exp)
abline(v=knots, lty="dotted")
points(tmp)

And here's the resulting plot, with represent the marginal relation of rcs(ablat,3) with the effect size: