[R-meta] Restricted cubic spline behaviour in moderator variable

Dear Wolfgang,

Thank you for your quick reply. You are right this is to do with the NA values, thank you. Those subtle NA differences ever so slightly affected knot position. That wouldn?t be an issue but as you point out, we are then using a different knot when using  attr(rcs(model.matrix(res.rcs)[,2], 4), "parms")

That alone didn?t fully explain my very aberrant plot that I saw. I tried playing around with subtle changes to knot positions (e.g. moving + - 1) and that to the naked eye actually made almost no change. The differences in the plots were because:

predict(model1,  newmods=rcspline.eval(seq(1,100, length = 100), knots, inclx=TRUE))) was using knot positions that were slightly different position to those that were parameterised in the model (because of how I incorrectly attained them), but in the predict/plot calculation this then resulted in stranger predicted points. Hence strange tail.

I guess I didn?t think of this, as when manually providing the knots, and then using the same predict and plot before, there obviously wasn?t any difference as they were based on exactly the same knot positions.

Thanks for your help here, I was going to resort to manually providing knots as per Harrell, 2015 recommendations (which is what rcs does?), but this will save some hassle.

Thanks
Mick

Dear Mick,

This might have to do with missing values. In the example on the metafor website (https://www.metafor-project.org/doku.php/tips:non_linear_meta_regression<https://www.metafor-project.org/doku.php/tips:non_linear_meta_regression>), both of these lead to the exact same results:

res.rcs <- rma(yi, vi, mods = ~ rcs(xi, 4), data=dat)
res.rcs

knots <- attr(rcs(model.matrix(res.rcs)[,2], 4), "parms")
knots

res.rcs <- rma(yi, vi, mods = ~ rcs(xi, knots), data=dat)
res.rcs

But let's make some of the yi values missing:

dat$yi[1:10] <- NA

Then rerunning the code above will lead to somewhat different results.

This happens because:

rcs(dat$xi, 4)

is actually based on the entire dataset, while

attr(rcs(model.matrix(res.rcs)[,2], 4), "parms")

uses only the rows actually included in the analysis (i.e., removing the first 10 rows).

If you use:

knots <- attr(rcs(dat$xi, 4), "parms")

then

res.rcs <- rma(yi, vi, mods = ~ rcs(xi, knots), data=dat)
res.rcs

will give the same results again.

I have changed the code in the example to use 'knots <- attr(rcs(dat$xi, 4), "parms")' as this is more accurate as to what is happening (even though it makes no difference in the example on the website since there are no missing values). Also, I've expanded the first footnote which gives a hint in this direction.

It would be nice if you could confirm that this solves the issue. If not, we will have to go digging deeper.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org] On
Behalf Of Mick Girdwood via R-sig-meta-analysis
Sent: Thursday, 31 August, 2023 4:17
To: Gabriel Cotlier via R-sig-meta-analysis
Cc: Mick Girdwood
Subject: [R-meta] Restricted cubic spline behaviour in moderator variable

Hi everyone,

I am fitting an rma.mv model with a restricted cubic spline moderator. I have
followed the helpful examples from the metafor website, as well as looking
through this mailing list for previous examples. All of that went well. I have
however run into an interesting issue.

I guess this relates less to meta-analysis as such and more specifically to the
rcs fitting/package. Here is my model:

model1 <- rma.mv(yi, V,
mods = ~rcs(timepoint, 4),
data = data,
random = list(~ timepoint|cohort),
struct = c("CAR?))

I have been comparing the model fit for different moderator structures using
heuristics (AIC, BIC etc) as well as visually looking at predicted plots. The fit
of this example seemed quite strange, with the tail of the spline at the last
knot angling well away from any data points. I was curious and went and extracted
the knot positions with:

knots <- attr(rcs(model.matrix(model1)[,2], 4), "parms?)
(e.g. in my example: 3.00000, 6.00000, 12.00000, 51.64118)

I then refit another model, but specifying these exact same knot positions

model2 <- rma.mv(yi, V,
mods = ~rcs(timepoint, c(3.00000, 6.00000, 12.00000,
51.64118)),
data = data,
random = list(~ timepoint|cohort),
struct = c("CAR?))

This time the fit is completely different (not subtle, they are two completely
different angles on the 3rd spline/last knot) - i.e. in model 2 the tail of the
spline is perhaps what would be ?expected' of the data. What is going on here? To
my mind I am fitting the same model with the same knot positions (as I used the
knot positions from the first model to fit the second), but have just specified
them differently, why are they behaving differently? Is there some other
parameter of rcs that I am not using/missing/misunderstanding? I tried reading
into the help and source for rcs but didn?t notice anything. I?m guessing this
wouldn?t be metafor related... I tried testing what the attributes of the two
different rcs calls looked like but they were identical too, so I don?t
understand how it can fit them so differently? I also tried fitting the same
model with splines::ns and the fit was identical to model2. Any ideas what is
happening with model1?

Sorry if this question is perhaps slightly off topic.
Thank you for your help as always.

Mick Girdwood
La Trobe University | Australia

La Trobe University | TEQSA PRV12132 - Australian University | CRICOS Provider 00115M