I've been using Cox for survival analysis with the goal of identifying the effect of 5 different levels of a treatment on survival. My dataset is not small and has 130 covariates after dummy encoding.

About 15 of these - including one level of my treatment - violate the PH assumption.

I think in my situation the best extended Cox option would be adding time interactions to those that violate the PH assumption meaningfully. Another option, of course, could be an AFT, if the assumptions there are also satisfied.

My question is: If I were to try each of these approaches, how could I compare which these are preferred given my objective. That is, how would I know which is better? It's not like I'm making a prediction about survival in which case I could use CI or something more precise. Thus, comparing which is best and knowing how to know this is not obvious to me.


1 Answer 1


The most rational and linear approach would be to put Bayesian priors on the amount of non-PH for all variables, before making use of information you get from examining fit to the observed outcome variable. This will give you the correct uncertainties for "final" effects in the model.

Less satisfactory but simpler approaches are (1) use the unitness correlation $\rho$ that goes with the test of PH to gauge the extent of the violation, instead of using P-values, which can detect trivial violations, or (2) fit an overall PH model then fit an overall non-PH model (e.g., log-normal survival time model) and see which one "wins" by some global goodness of fit assessment. I'm less inclined to pick and choose covariate-by-covariate.

  • $\begingroup$ Thanks Professor Harrell. Can you explain what you mean by ρ here? When I print 'cox.zph' here (on R) I see chi-squared, df, and P-value. $\endgroup$
    – JED HK
    Commented Aug 10, 2022 at 13:15
  • 1
    $\begingroup$ Hmmm not sure why this was taken out of cox.zph. You can run Spearman's correlation of survival time vs r <- resid(fit, 'scaledsch') to get it. I think you can get the survival time from as.numeric(dimnames(r)[[1]]). $\rho$ gives you a unitness index of how correlated $\beta(t)$ is with $t$, which is a sample-size-less index of non-PH. $\endgroup$ Commented Aug 10, 2022 at 14:36
  • $\begingroup$ It was removed a while ago (cran.r-project.org/web/packages/survival/news.html). I followed your advice and now must interpret the output. Most ρ are around 0. The largest in magnitude away from 0 was 0.064 (though this variable was not significant according to cox.zph. This suggests correlation between β(t) and (t) is low, I presume. However according to Schoenfeld plots some hazards do seem to change quite a bit. How can I reconcile this? $\endgroup$
    – JED HK
    Commented Aug 10, 2022 at 15:49
  • $\begingroup$ I could not find mention of removal of $\rho$ in the link you provided. In terms of reconciling, check the scale on which the residual plot y-axis is drawn. Sometimes the scale exaggerates small effects. 0.064 is not very large. $\endgroup$ Commented Aug 10, 2022 at 16:00
  • $\begingroup$ Ah, sorry about that. I took it from EdM's answer here: stats.stackexchange.com/questions/144923/… Small enough to ignore, then? Feel free to link anything. You have one example in your book in Chapter 20 in which ρ is -0.23 and highly significant but you do not comment on what -0.23 means in terms of magnitude $\endgroup$
    – JED HK
    Commented Aug 11, 2022 at 8:01

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.