I've used pspline to add restricted cubic splines to certain continuous covariates that violated the PH assumption upon first analysis. I identified these covarites using cox.zph and also used Schoenfeld plots to plot the offending covariates (p<0.05).

After using pspline and re-running cox.zph I see that my offending variables are no longer significant. For completeness, however, I'd also like to show a new plot that shows the violation has been "fixed" because of the spline (i.e., that time-dependency has been accounted for).

However, when I make Schoenfeld plots for the new model (that also incorporates pspline) the plot looks hardly any different. I guess this is because the spline terms are all associated into one coefficient (fourth comment on the answer) but this then begs the question: how instead can I "show" graphically that the PH violation has been dealt with in a similar way that I "show" PH is not satisfied using Schoenfeld plots. I should be consistent with the reporting of my results.

Edit: Number of events in dataset = 1600, N=100k

  • $\begingroup$ The pspline() function doesn't provide restricted cubic splines. Instead it provides a penalized smoothing spline whose coefficients can't be readily interpreted on their own. If you want restricted cubic splines, try the rcs() function in Frank Harrell's rms package. Then specify terms=FALSE in the call to cox.zph() and you can get plots for each (unpenalized) spline coefficient separately. $\endgroup$
    – EdM
    Commented Aug 16, 2022 at 15:07

1 Answer 1


You seem to be mixing two effects: (1) overall effects of covariates modeled using spline functions, and (2) non-proportional hazards of covariates which can be thought of as interaction terms involving covariates and follow-up time (though not as simple as adding regular interaction terms to the model; the Cox model uses as special form of the likelihood function to incorporate these interactions). In your post you did not specify how (1) is tested. This needs to be done with multiple degree of freedom composite (aka "chunk") tests (Wald or likelihood ratio $\chi^2$ tests). For (2) you took the quick and appropriate approach of using a test of PH based on scaled Schoenfeld residuals (R survival package cox.zph function). cox.zph computes a multiple d.f. test of PH and also has an option to compute test statistics on overall spline "terms". Both of these are helpful when covariates are expanded into splines.

The general approach to the problem uses tensor splines: a spline in a covariate having all its terms interact with the spline components of follow-up time.

Please state the number of events in your dataset, as these analyses are often limited by that.

Don't use p<0.05 for making arbitrary decisions about modeling.

  • $\begingroup$ Thanks for the answer Professor Harrell. For (1), do you mean here comparing likelihood ratio of models with and without splines? If so, likelihood is lower with the splines. My "concern" I guess, Professor, is that before adding spline terms I use Schoenfeld to check also nonPH (i.e., I don't rely solely on significance of cox.zph. If, then, I modify terms in the model using splines, then I should also support that PH has been "fixed" with a similar visual method, otherwise I am being inconsistent $\endgroup$
    – JED HK
    Commented Aug 15, 2022 at 14:17
  • $\begingroup$ It is unclear why you would even look at a model that forced effects to be learn. I would have pre-specified a spline model because nature is seldom linear. $\endgroup$ Commented Aug 15, 2022 at 15:47
  • $\begingroup$ My understanding is we make the linear assumption for simplicity and accept that it's rarely true but close enough. In any case, I don't appear to have linearity issues with most of the predictors (based on Martingale plots). I just need to straighten out the PH violations and I think I should be basically done $\endgroup$
    – JED HK
    Commented Aug 15, 2022 at 15:57
  • $\begingroup$ By allowing for different possiblities (e.g., linear vs nonlinear) you are getting over-optimisitically small standard errors of the "final" parameter estimates, and possibly meaningless p-values. Best to allow for nonlinearity for effects not known to operate linearly. $\endgroup$ Commented Aug 15, 2022 at 19:56
  • $\begingroup$ Thanks Professor Harrell $\endgroup$
    – JED HK
    Commented Aug 15, 2022 at 20:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.