# Cubic splines in Cox model

I have a question about the cubic splines used in the Cox model to test the linearity for the continuous variables. I read that usually the knots chosen are the quantiles. Can you find different results if you change the knots? I read that it is an important choice related to the number of knots. If you choose too many knots compared to the sample size you can have overfitting. I was wondering instead about the location of the knots. If I change the location of the knots, can this affect the results?

• Changing the positions is another axis of variation that can lead to overfitting. So the answer is yes. And aside from potential effects on over fitting, it will certainly have an effect on the results, that final question is certainly positive. Commented May 20 at 15:12

The question is about varying knot locations if the number of knots is fixed. I go into this in https://hbiostat.org/rmsc/genreg . With a restricted cubic spline (natural spline) the smoothness constraints are so strong that changing knot locations doesn’t have a very noticeable effect unless you put two knots very close together. In the link you will see links to interactive tools that will allow you to easily change knot locations and see how the fit changes. The online tool from Paul Lambert is one to look at first.

First, when you say something like

the cubic splines used in the Cox model to test the linearity for the continuous variables

you need to be careful about the type of spline being used.

For example, the pspline() function that's part of the R survival package, and seems to be preferred by Therneau and Grambsch, uses a penalized spline with knots spaced evenly along the range of predictor values. That starts with a larger number of knots and associated coefficients than restricted cubic splines for the same number of degrees of freedom, then penalizes the coefficient magnitudes to enforce the desired level of smoothness.

The restricted cubic splines implemented in Frank Harrell's rcs() function or the ns() function in the splines package don't involve penalization in that way. Smoothing is controlled by the number of knots and their placement. Default placement of inner knots is by quantiles of the predictor values. Unlike with pspline() you can adjust the placement of inner knots for those functions from the defaults.

Second, you have to watch out for the placement of the "boundary knots" near the limits of the range of predictor values. Linearity is enforced outside those boundary knots. A change in the locations of the boundary knots will affect the shape of the curve near the extremes.

The ns() function and the pspline() function by default place those boundary knots at the extremes of the predictor values, so there is no region of enforced linearity. In contrast, the default in the rcs() function places boundary knots within the range of predictor values. For example, with fewer than 100 non-missing values, the boundary knots are placed at the 5th largest and 5th smallest values. See the help page for rcspline.eval() in Frank Harrell's Hmisc package.

All those functions (including pspline()) allow you to specify different locations for those boundary knots and thus the shape of the curve near its extremes.

• Kudos for getting me to realize that placement of the outer knots in a restricted cubic spline is likely to be more important than placement of the interior knots. And I don’t like the default behavior of the other methods where the outer knots are placed at the extremes of the data. That’s not very robust to outliers. Commented May 21 at 12:22

So far, we have two answers, each from a very esteemed writer. Shawn says "of course it can" and demonstrates that. Frank says "it usually doesn't". Both are correct.

But I think the question is about a somewhat more specific aspect of using RCS: Testing linearity. AFAIK, this decision will always be somewhat subjective. Yes, you can compare a straight line to a RCS. Graphically, you could add a straight line fit the excellent plots in Shawn's answer. But then you have to decide: "Is the difference between the RCS and the straight line big enough to make it worthwhile?" because RCS is more complex to interpret.

The decision will also depend, in part, on things that can't be assessed statistically:

• What is the purpose of the model (e.g. prediction, interpretability, or both)?
• Who is the audience and how sophisticated are they, mathematically?
• Will the audience be interpreting the numerical output of the model, or will a graph be enough? Or perhaps they are only interested in the actual predictions?

Also, different though Shawn's three plots are, I think they might all result in the same answer to the question "is it worth it?" I'm just not sure what that answer will be!

• RCS = restricted cubic splines. Commented May 21 at 7:25

If I change the location of the knots, this can affect the results?

Certainly, this is true of any spline fit. You are changing where the basis functions are being constructed, so naturally this will shift where the line goes. Whether or not it is wise to do so depends on the data.

Can you find different results if you change the knots?

Yes. For example, if I fit 50 degrees of freedom using a standard cubic smoothing spline, I will overfit the data to extremity, but that won't necessarily mean it is right. The knot location and number of knots will again depend on the most parsimonious yet accurate fit to the data.

As an example, compare the three fits to the data below, one with $$\rm df = 4$$, another model with $$\rm df = 30$$, and one where the knots are changed with $$\rm df = 4$$ still intact but with different locations. You can see not only the results in the summaries change, but the fits change too. The first model is slightly curvilinear (and seems to model the data well enough), the second model is so overfit with knots that it simply chases the noise, and the third model with different locations interpolates poorly at the end of the scatterplot.

#### Fit Survival Model ####
library(survival)
fit1 <- coxph(Surv(futime, death) ~ bs(age, df=4), data=mgus)
fit1

#### Knot Number ####
fit2 <- coxph(Surv(futime, death) ~ bs(age, df=30), data=mgus)
fit2

#### Knot Location Change ####
fit3 <- coxph(Surv(futime, death) ~ bs(age, df=4,
knots = c(1,2,80,85)), data=mgus)
fit3

#### Plot ####
par(mfrow=c(1,3))
termplot(fit1, se = T, partial.resid = TRUE)
termplot(fit2, se = T, partial.resid = TRUE)
termplot(fit3, se = T, partial.resid = TRUE)