# Cox PH model: managing continuous variables and linearity assumption

In an epidemiological study, I'm using martingale plot to assess the linearity of continuous variables.

Here are the Martingale Residuals (from Null Model) using R's survminer::ggcoxfunctional() output for 2 variables, on which we see that the linearity assumption is violated. Please note that since I have left-truncated data, the timescale is age (start is age at inclusion and stop is age at event or censoring, see this or this).

ggcoxfunctional(Surv(age_start, age_stop, event) ~ x1 + x2, data=db)  These are zoomed-in images so user can see the scaled functionnal form of the curve, the lowess account for other points which you can see on unzoomed images here (x1) and here (x2). Indeed the curve form is quite different and non-linearity could be leaved unseen in the unzoomed picture.

Since these variables are for adjustment only, I don't need a rock-hard precision, so I would like to break it into a categorical factor. My problem is to decide which breaks to decide. Since there is no consensus or clinical evidence on where to break these variables, I see 2 main options:

• break around the median or some quantiles
• use the Martingale to decide where to break

Is deciding using the Martingale plot better ? If yes, how to decide ?

Side question: depending on the source, the reading of the martingale plot is quite different. I saw two interpretations: "the curve should be somehow linear", and "the curve should be somehow linear and parallel to the x axis". Which one is the right one ?

## EDIT :

Probably related (but unanswered) : What if linearity doesn't hold in a cox model?

I obviously should go for the spline method, but I couldn't find any ressource which explains it practically.

For instance, let's consider my case, with y my measured variable and x1, x2 and x3 my adjustment variables (confounding factors):

coxph(Surv(age_start, age_stop, event) ~ y + x1 + x2 + x3, data=db)


For what I witness with Martingale plots, x1 and x2 violate the linearity assumption, but y and x3 are fine (the proportional hasard assumption is OK for everyone). From what I understood, I should apply splines on my "guilty" variables, and so write something like this:

library(rms)
coxph(Surv(age_start, age_stop, event) ~ y + rcs(x1) + rcs(x2) + x3, data=db)


Unfortunately, the output is not easily understandable and the help (?rcs) is not providing much help. (Sorry Pr Harrell, this seems a great library but it is quite not beginner-friendly).

How can I correct my variables so the linearity assumption is not violated, without cheating by looking too much at my data ?

• I'm not familiar with survminer::ggcoxfunctional(), but if the red dots represent the individual residuals and the solid lines on your plot are supposed to be a lowess fit of them, then the lowess doesn't seem to fit the residuals very well. Also, the mean of martingale residuals over all cases should be 0, while none of the red dots seems to have a y value greater than 0. Might there be some extreme values beyond the plotted range that are distorting the apparent results? – EdM Aug 13 '18 at 14:30
• @EdM You are right, this is odd. I could be due to my large effective (n=30k) but I have no real clue. I use this function based on the explaination from this link. How could I track extreme Martingale residual values ? – Dan Chaltiel Aug 14 '18 at 12:16
• Think carefully about whether you really have "left-censored" data; it might be better to include age at study entry as a predictor rather than formulating the Surv() function in the way you show in your question. – EdM Aug 16 '18 at 14:25
• @EdM Maybe you misread my edit: my data are left-truncated, not left-censored. This point is actually a well thought decision with my supervisor, I think we are pretty sure about it. – Dan Chaltiel Aug 16 '18 at 15:03

Frank Harrell's recommendations get to the heart of the matter in general, but there may be issues with your dataset that are also giving you problems.

It seems that there must be outliers of x1 and x2 that are beyond the limits of the values displayed in your plots. I see no other way for almost all the martingale residuals to be at or below 0 and for the lowess fits to rise at both ends in the ways that they do. Your first task is to figure out what's going on with that issue and clean up the data if necessary. Examine the actual distributions of both those predictor variables and the martingale residuals (available with the residuals() function applied to a coxph object), making sure that you can see all the data points. That might indicate some problems with the data set that, once fixed, could remove the nonlinearity issue.

If you still want to allow for nonlinearity with respect to those predictors, restricted cubic splines are almost certainly the best choice. Although the math might be complicated, just think of a restricted cubic regression spline as the simplest smooth curve that best fits the data subject to some simple constraints: the numbers and locations of the knots on the x-axis, and linear rather than polynomial relations at the ends of the data range to prevent awkward behavior. Chapter 2 of Harrell's course notes discusses how to choose numbers and positions of knots. With 30,000 data points you can use many knots, but 5 are often quite adequate. Specifying the number of knots to rcs() and using the default positions typically works well.

If you are not directly interested in the relations of x1 or x2 to outcome and are simply controlling for them to analyze your predictor of main interest, then you don't have to worry much about the initially confusing display of p-values and hazard ratios for each of the terms associated with the splines. If you do care about their relations to outcome, take advantage of the facilities provided by anova() in the rms package to get an overall estimate of their significance including all those terms. Yes, there is a learning curve in starting to use that package. In particular, it's important to save the output from the datadist() function as run on your (cleaned) data frame (e.g., dd <- datadist(mydata), and then specify that as an option (e.g., options(datadist="dd") before running any of the linear modeling functions like cph(). I recall my initial hesitancy at starting to use that package and only wish that I had started even earlier than I did.

• Great answer, thank you very much ! Indeed, pictures were automatically zoomed-in. You can now see whole pictures links in my edit below them on which the lowess curve seems OK, but the fonctionnal form of the curve changes drastically. OK, I'll try to dive into this, this really seems worth it. One last question: shouldn't I test the linearity and th PH assumption on my rcs-transformed x1 and x2 variables ? How could I ? Also, if I'd care about their relation to outcome, how could I get a Hazard Ratio for the variable ? – Dan Chaltiel Aug 16 '18 at 13:34
• @DanChaltiel each of the extra spline terms is a predictor for which you can get HRs and check PH. The individual HR, however, are hard to interpret. You want to display the overall effects of x1 and x2; Harrell's class notes linked in my answer (e.g., chapter 18) show how to plot hazard versus value for spline fits and how to produce nomograms for that. For "linearity" the issue is whether the spline fits adequately adjusted for nonlinearity; you could check the significance of coefficients for extra knots. Validation and calibration of the overall model is probably much more important. – EdM Aug 16 '18 at 15:33
• OK, I get it, since the linearity is not valid, there is no individual HR for this variable. One last precision please: by adding rcs(x1, 5) in place of x1 in my formula, I get my outcome HR rightly adjusted, but the PH assumption is now violated (cox.zph() pvalue for all splines and for GLOBAL <<0.05). Is this a problem ? If yes, how to solve it ? – Dan Chaltiel Aug 20 '18 at 9:07
• @DanChaltiel the full plot of x1 suggests that the nonlinearity is dominated by the skew with a few cases having values of 4000 or 5000 while most are much lower.. See what happens when you try a log transformation of x1 instead of the spline fit. If PH is a problem for a continuous predictor, you can break up the cases into strata based on values of that predictor, allowing for different baseline hazards for each stratum, and treat it as a stratified predictor with strata() for coxph() fits or strat() for cph() fits. – EdM Aug 20 '18 at 9:37
• Sorry, I used x1 as an example, my question was broader, like "how to handle added PH problems when adjusting continuous variable with splines ?" Actually, in this dataset, the very high effective seems to remove the linearity problem (the full plot seems somehow linear and only the first spline member is significative for each variable), but I would like to understand what to do for other dataset. So, if the conclusion is to get back with strata, I have to get back to the original question: how to choose breaks when there is no clinical clue ? – Dan Chaltiel Aug 20 '18 at 9:44

Don't trust those plots were the data run out, and add confidence bands to the plots.

You are using a multi-step analysis for which model uncertainty will not be recognized in the final model fitting step. This will significantly distort confidence intervals and p-values.

A better approach is to just specify how many parameters you want to devote to that predictor, and fit that many using a regression spline. Stick with that result no matter how linear the fitted relationship looks. This will preserve all aspects of statistical inference. For details see my RMS course notes.

• Thanks for your answer. Unfortunately, my statistics level is not high enough to understand your course yet. I hope I will get there soon, but in the meantime, is there any easy way to do this in R ? I hope making it will help me to understand it. – Dan Chaltiel Aug 13 '18 at 13:35
• Also, are you saying that survminer::ggcoxfunctional() output is genuinely wrong or that Martingale plots should not be used this way ? If so, what is the right way ? – Dan Chaltiel Aug 13 '18 at 13:37
• @DanChaltiel I think the point is that if you use martingale plots this way then you will be selecting aspects of the model based on the data so that the usual assumptions for calculating p-values and confidence intervals won't hold. Splines are easy to incorporate into a model via the rcs() function provided in the rms package in R. For a continuous predictor they should provide much better control for nonlinearity than would choosing breakpoints. If you choose the number of spline knots without reference to the data, then p-values will be reliable. – EdM Aug 13 '18 at 15:19
• @EdM ok I understand the problem, this makes perfect sense. I've been reading a lot about splines since yesterday, they seems easy to use but not so much to understand. How should I choose the number of knots ? And more broadly, what algorithm should I use ? Let's say I compute my model with rcs(x1), what should I do next ? Thanks a lot for your help (BTW I think any answer to this comment could go as an answer and not a comment) – Dan Chaltiel Aug 14 '18 at 9:56
• I edited the title and the end of my question for clarity as I understand better what the problem really is. – Dan Chaltiel Aug 14 '18 at 10:29