Violated non-proportional hazards - Cox regression model of time-dependent covariable

Question

I have been reading a lot on Cross Validated for a while now, but this is my first post here. Thanks for all the help I've received! I am currently working on a project where we are trying to estimate the effects of a drug on overall and disease-specific mortality. This is a prospective observational study, where patients are enrolled when they have a blood test, which happens at different time points. Patients are then followed until death, emigration, or the end of the study. To figure out if the drug impacts mortality, we have set up a Cox regression model where drug exposure is allowed to accumulate over time with the cumtdc function specified here.

The problem is that the proportional hazards assumption is severely violated - with p-values approaching 0 for most of the covariates. I am a beginner to this type of analysis and what I have read and tried so far is:

• Splitting continuous variables (such as age) into categories, which did not make a difference.
• Could interaction variables be helpful here - if so, how do I best implement them?
• Would this analysis be better suited to be performed with an accelerated failure time model?
• Is it possible that the analysis could be valid despite the non-proportional hazards? I am asking this since, intuitively, there are many individuals in the study (> 300 000) and many events (> 30 000), so even minimal violations might be registered as significant.

I am attaching the Schoenfeld plots with the respective p-values for the covariables.

I would be most appreciative of any suggestions to this problem. There seem to be many different approaches out there, making it a bit tricky to navigate for a beginner. Thank you!

ADDITON TO ORGINAL QUESTION

These are the log hazards for the variables age and X. I have now attempted both restricted cubic splines and categorization without meaningful changes to the Schoenfeld residuals.

ADDITION 2 TO ORIGINAL QUESTION

The output from cox.zph:

The requested plots from cox.zph for A and age:

Answer 1

With a data set this large, you are facing a problem similar to that with normality testing: any large real-world data set will tend to show a statistically "significant" deviation from an ideal. You thus have to engage closely with the data, using your knowledge of the subject matter to evaluate whether any deviations from proportional hazards (PH) are large enough to matter. This answer addresses the problem in the PH context, with suggestions for some graphical evaluations and for further reading.

With respect to your data set, it looks to me like most of your predictors don't deviate very far from constant hazards over time. The apparently categorical predictors (the last 6) seem quite flat. The coefficient for predictor A does, however, seem to show a tendency to increase over time. There also seem to be more high-positive residuals for X at early times than at late times. You might want to evaluate those more closely, based on your understanding of the subject matter. With continuous predictors (as A and X seem to be), an improper transformation might be the underlying problem, if their values aren't linearly associated with log hazards. Flexible modeling of those predictors with restricted cubic splines might fix the problem.

With the risk of survivorship bias, modeling with time-dependent covariates as you are doing is tricky enough. Be extra cautious if you think that you also need to move to time-dependent coefficients as in the vignette that you cite. Be warned that if you use the time-transform formalism to express coefficients as a function of time, the cox.zph() function can't be used to evaluate PH. This answer describes that and some other things that you might need to consider when using recent versions of the R survival package.

Answer 2

I have a very similar related question.

Not an expert, but based on what I've researched so far and what is highlighted on this page, here are some points that may be helpful in answer to your last question

Is it possible that the analysis could be valid despite the non-proportional hazards?

1. As others have mentioned, large sample size may be a factor leading to the violation of the statistical test of the PH assumption. As @EdM mentioned, it will depend on your context/knowledge of the subject as to whether the deviations matter.

From page 1461 of this paper

When the sample size is small, this method may lack power to detect deviations from PH; while for large sample sizes, hypothesis tests may be over sensitive to slight deviations from this assumption.

2. The answer to this question suggests that the effect of fitting a Cox model with non-proportional hazards is "slightly less power" which can be recovered with robust standard errors, leading to the hazard ratios being interpreted as the time-weighted average of the hazard ratio

3. This interesting paper (Why Test for Proportional Hazards?) was published recently highlighting that there are legitimate reasons to assume violation of PH, and that one of the effects is on the interpretation rather than invalidity of the results. They actually suggest that statistical testing of PH is unnecessary if it's expected that the hazard ratio varies over time. (I'd be interested to hear what others think of this paper)

In the section "How should hazard ratios be interpreted?" on page 1402

As a weighted average of the time-varying hazard ratios, the hazard ratio estimate from a Cox proportional hazards model is often used as a convenient summary of the treatment effect during the follow-up. However, a hazard ratio from a Cox model needs to be interpreted as a weighted average of the true hazard ratios over the entire follow-up period. The 95% confidence interval should be estimated using a valid method such as bootstrapping and also using inverse probability weighting to adjust for losses to follow-up.