I have a question on the cox proportional hazard model, in particular the proportionality assumption. I use cox.zph() function in R to check whether the proportionality assumption is satisfied. The global test of cox.zph() has p-value much greater than 0.05, but one of the covariates has p-value close to 0.01. In this situation, should I do a remodelling (can be stratification or time interaction), or should i say that the assumption is still not violated since globally its not significant? Thank you
2 Answers
The global test of proportional hazards is not well-calibrated. You haven't controlled for multiple comparisons. It's difficult to gauge power of the test. $\alpha=0.05$ is probably too lax in most sample sizes. The test is arbitrarily powerful in large sample sizes. It's possible that the covariate you identify is a spurious finding, and that it arises from natural variability in the observation of time-to-event data.
Even if the hazards were not proportional, altering the model to fit a set of assumptions fundamentally changes the scientific question. As Tukey said, "Better an approximate answer to the exact question, rather than an exact answer to the approximate question." If you were to fit the Cox model in the presence of non-proportional hazards, what is the net effect? Slightly less power. In fact, you can recover most of that power with robust standard errors (specify robust=TRUE
or cluster = ~id
). In this case the interpretation of the (exponentiated) model coefficient is a time-weighted average of the hazard ratio--I do this every single time.
When the actual hazard ratio over-time is of interest, there are flexible methods of estimating its value. You may create a flexible, polynomial representation of time using basis splines and fit their interaction with the covariate(s) to estimate a hazard ratio time function. The power of the Cox model may be compromised by this. Using a parametric exponential survival model with spline adjustment for time can approximate the semi-parametric inference of the Cox model very well, and is better powered to detect interactions of time with one or more covariates.
If hypothesis testing is your main goal, you should not do anything at all and stick with the model you had anticipated using before seeing the data.
Changing the pre-defined model based on the results of the cox.zph() test can lead to biased estimates and invalid p-values.
See SiM paper: https://harlanhappydog.github.io/files/SiM.pdf
p < 0.05
in the PH test violates the proportionality of hazards assumption. Hence, as you said, you need to do either stratified analysis (if the variable that violates the assumtion is categorical -- i.e., nominal or ordinal) or run a Cox model with time-dependent covariate (i.e., time interaction). $\endgroup$