1
$\begingroup$

I am performing Cox regression modelling, searching for predictors for various clinical outcomes. Now, based on Schoenfeld residuals test, some of my binary predictors have p-value <0.05 with global p-value also <0.05. This means it violates the proportional hazard assumption / it is a time-dependent variable, am I right? I tried including tt(predictor) in my model but then cox.zph does not recognize it and I cannot test for Schoenfeld. I also tried including "predictor:time" in my model but then it spits 2 separate predictors "predictor0:time" and "predictor1:time" with no reference group. I don't want to do time splitting as it would not make sense clinically to interpret that variable in different time periods. In this circumstance, can I just use weighted Cox (cox.phw)? How can I test for goodness of fit with cox.phw? An example plot of weighted cox object:enter image description here Mind you; I am a clinician with no strong biostats background, just trying to learn and apply these methods, so please bear that in mind when reading and responding. Thank you!:)

$\endgroup$
2
  • $\begingroup$ Can you show the KM of survival for the separate covariates. It's easy for some to visualize "non-proportionality" on the survival scale. $\endgroup$
    – AdamO
    Dec 8, 2022 at 19:06
  • $\begingroup$ If the curves are only a "little" non-proportional, you can just use the robust=TRUE option to calculate sandwich based errors - i.e. asymptotically correct inference for a difference in survival. $\endgroup$
    – AdamO
    Dec 8, 2022 at 19:09

1 Answer 1

1
$\begingroup$

It's important in this case is to look at smoothed plots of scaled Schoenfeld residuals as a function of time for the predictors in question. That's true for several reasons.

First, with a large number of events you might find "statistically significant" violations of proportional hazards (PH) that are of small magnitude. You don't know how big a problem you have until you look.

Second, those plots suggest specific forms for the time dependence of the regression coefficients. The "time dependent" vignette for the R survival package illustrates that in Section 4.

Third, if you then choose to model a time dependence with the tt() functionality, as you note you can't then use cox.zph().* You can, however, plot the form of the resulting tt() function atop the original plot of scaled Schoenfeld residuals to see how well it works, as demonstrated in the last plot in Section 4 of the vignette. Just as you are unlikely in practice to get data that perfectly follow a normal distribution, you are unlikely to get perfect PH. With such a plot, you can apply your knowledge of the subject matter to asses whether the failure of PH is large enough to matter.

In terms of using the coxphw package, be careful that you aren't throwing away important information about the association of the covariate with outcome over time. In a standard Cox model that fails PH, you get event-weighted average hazard ratios that run a risk of bias. The coxphw() function seeks to remove that bias. If you aren't specifying an explicit function of time with that function, you are accepting the failure of PH (so there's no need to test for it) and getting a less-biased estimate of an "average" hazard ratio. But what if the change of hazard over time is clinically important, as with the Karnofsky score in Section 4 of the time-dependent vignette?

You can specify functions of time for covariates with coxphw(), and its vignette seems to provide code for a cox.zph() type of analysis for such models. It's based, however, on code adapted from Version 2.37-4 of the survival package and I don't know whether it fixes the error* with time-dependent coefficients.

Finally, don't feel limited to Cox PH models. Accelerated failure time (AFT) models (e.g., log-normal, fit by survreg() in the survival package) assume covariate-associated expansions/contractions of the time axis instead of proportional hazards. Such a model might fit well if PH doesn't hold in your Cox model. (Unless it's a Weibull model, a parametric AFT model necessarily can't have PH.) Aalen additive models (aareg() function) fit both the baseline hazard and coefficients as functions of time; whether that extreme flexibility is a blessing or a curse depends on your data and how carefully you interpret the model.


*According to the NEWS document of the survival package, you were able to use cox.zph() on models containing a tt() term in versions prior to 3.1-2, but the results were wrong. Keep that in mind if you find such plots in your searching for answers.

$\endgroup$

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.