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

[Viechtbauer, Wolfgang (NP)]

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), 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 using 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