# Evaluating Functional Form in Cox Regression using rcs

I was approaching Cox regression and need to evaluate whether is more suitable to model my independent variable as a non-linear prediction of the outcome.

Main questions is:

1. How do I check whether or not my independent variable should be modelled in a non-linear fashion?

I came into several methods, one of which was to plot the variable versus the martingale-based residual of the Cox regression (See for example www.sthda.com/english/wiki/cox-model-assumptions, last figures).

I therefore produce this type of plot with the following code, being var my independent variable:

time <- data$$time status <- data$$status
S <- Surv(time, status)
model1 <- coxph(S~1, data = data)
data\$residual <- residuals(model1, type = "martingale")
ggplot(data = data, mapping = aes(x = var, y = residual)) +
geom_point() +
geom_smooth()


Here is the result:

Here comes the other questions:

1. As you can see, in my plot martingale-based residuals seems to be a lot "clustered" at extreme values, different from other plots that I found online in which residuals were more or less "uniform". I don't know how to interpret this data - seems like the data are "splitted".
2. Based on the plot, should I assume a linear relationship? I was considering to use a restricted cubic spline with 4 knots, since my guess was that a non-linear relationship may be the case here.
3. What is the correct interpretation of this plot? Do anyone have a guidance paper to interpret this plots and/or more generally speaking a guide to Cox diagnostics?

[EDIT]

This is the result of the model and of the anova as requested by @EdM - actually to generate this result I have fitted model2 with cph.

Model:

      Coef    S.E.    Wald Z Pr(>|Z|)
var    4.5077  1.5089  2.99  0.0028
var'  -0.1702  4.8386 -0.04  0.9719
var'' -5.4636 12.5472 -0.44  0.6632


Anova:

                    Wald Statistics          Response: S

Factor     Chi-Square d.f. P
var        217.03     3    <.0001
Nonlinear  17.17     2    2e-04
TOTAL      217.03     3    <.0001


Instead of trying to guess from the shape of the plot of null-model martingale residuals against the value of your continuous predictor, why not let the regression itself tell you whether non-linear terms are necessary?

Harrell discusses restricted cubic splines in Section 2.4.5 of his class notes. If you fit with a restricted cubic spline, you get a set of associated coefficients, with the first being for the linear term and the rest being the coefficients of the non-linear terms. A joint test of the hypothesis that all of the non-linear term coefficients equal 0 is the test you need. See the end of that Section of the notes. If you use the rms package and its cph() command for Cox regressions, the anova() function applied to the cph object includes that test.

This regression approach has the advantage that you can apply it while taking other predictors into account. I can't say exactly why you have that particular distribution of residuals, but it might be due to other outcome-associated predictors that aren't accounted for in the null model used to generate the plot. See this answer and its link.

In response to comments and edited question:

The anova() results after modeling with the rcs() term suggest a weak overall non-linear association, dominated by the strong linear term (coefficient for var in the model results). That's consistent with the overall shape of the smoothed plot of martingale residuals.

Even though neither of the nonlinear coefficients (for var' or var") is individually distinguishable from 0 based on its own standard error, the multi-parameter Wald test used by anova() in this context also takes into account the covariance between the two coefficient estimates for non-linear terms. There's probably a very strong negative covariance between those estimates, leading to evidence that the two non-linear coefficients together have a "statistically significant" (p = 0.0002) association with outcome.

Whether that's sufficiently important in practice to devote the 2 extra degrees of freedom to non-linear modeling of var is up to you, based on your understanding of the subject matter and of your audience. With a large study like this I see no harm in keeping the non-linear terms, unless you have many other predictors to include in the model. The simple linear association might, nevertheless, be easier to explain. See Harrell's class notes, linked above, for further guidance.

Thoughts on the appearance of the martingale residual plot

Unlike some other residual types, martingale residuals are calculated for all cases, not just for cases with events. They are effectively the differences between observed and predicted event numbers for the cases. This means:

1. Martingale residuals for censored cases are necessarily negative (fewer "observed events" than predicted), so censoring patterns can influence what you see.

2. They add up overall to 0, setting constraints on the overall distribution of values (in part determined by the censored cases).

3. A null model like yours only "predicts" a combined survival curve, so you might find "extreme" residuals like this when there are many censored cases.

4. Some martingale residual plots in the literature and on line are based on multiple-regression models that include many predictors, with just one predictor of interest removed to try to gauge its functional form. As those models do a lot more "predicting" than a null model, you might find a more "uniform" distribution of residuals in those situations, when there is censoring.

• Thank you for your answer! I performed the anova and the output showed a p <0.0001 for "Non-Linear" term of "var", although the coefficients in the Cox Regression (var', var'') were highly non-significant. I don't know how this can happen but still, I deduced that a non-linear fit may be worthy. May 8, 2021 at 7:40
• @pankevedmo I'm a little concerned that there might be some confusion in how to interpret the anova() test on the Cox regression. It would help for you to edit your question to show those results. That will both give guidance to others who come across this question and help determine just why you are getting a significant "Non-Linear" anova() result with non-significant spline coefficients for var' and var''.
– EdM
May 8, 2021 at 16:19
• done! I've added the result of the anova. May 8, 2021 at 16:47
• @pankevedmo could you also include the result that suggests non-significant spline coefficients for var ' and var''? From the anova() output it looks like most of the response is from the linear term, with the "significance" of the non-linear term mostly due to the very large sample size.
– EdM
May 8, 2021 at 16:52
• there you go. Can you please help me how to interpret this output? May 8, 2021 at 19:44