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

[Viechtbauer, Wolfgang (NP)]

My preference is to set the knot positions manually anyway. That way, one is always going to use the same positions in all calls of rcs(), rcspline.eval(), and so on. And one can then avoid that knots are set at values like 51.64118 (I would just round this to 52).

Best,
Wolfgang

>-----Original Message-----
>From: Mick Girdwood [mailto:M.Girdwood using latrobe.edu.au]
>Sent: Thursday, 31 August, 2023 15:12
>To: Viechtbauer, Wolfgang (NP)
>Cc: R Special Interest Group for Meta-Analysis
>Subject: Re: 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
>
>On 31 Aug 2023, at 18:21, Viechtbauer, Wolfgang (NP)
><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>
>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
>
>La Trobe University | TEQSA PRV12132 - Australian University | CRICOS Provider
>00115M