1
$\begingroup$

I'm using cox regression to analyse my data. The explanatory variable is a congenital disease (X) and the outcome is an another disease (Y), which is a comorbidity of the congenital disease. I'm looking at the association between the two diseases. The individuals in my data that has a congenital disease, may be diagnosed with the other disease much later in life.

I have done a COX-PH, and I see a significant association between X and Y. But the proportional hazard assumption is violated (p is below 0.05). I have looked into this https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf which many of you have recommended almost every time someone asks about fixing the violation of proportional hazard assumption, and thank you for that!

I used the example in chapter 4 (page 16-19), and it fixed my problem. The proportional hazard assumption is no longer violated when dividing my data into time-intervals. In my case and in the example from chapter 4, the data is divided into 3 time-intervals. But this also means that I get 3 HR's (one for each time interval).

So what if I'm interested in only getting one HR? For my whole dataset? Is that really impossible?

and no, I can't use the tt() function, as my data is bigger than N=100.000. Every time I try to run tt(), then my R would run more than 20-30 minutes, and then I stop the process (I have never seen it stop by itself).

and as I have understood it, I can't just use the time interaction term either. according to the link I provided, the interaction term is most suitable for continuous explanatory variables? Maybe I'm wrong? Mine is categorical. But I have tried that too, and my proportional hazard assumption was still violated.

So if there is no option than having more than one HR, would I be able to make a metaanalysis, that would provide me one HR? Would that make sense?

Hope you can help me and clarify things!

$\endgroup$
4
  • 3
    $\begingroup$ With such a large data set, it's quite possible to have a "statistically significant" violation of PH that isn't of practical significance. Therneau and Grambsch say in Section 6.5: "A 'significant' nonproportionality may make no difference to the interpretation of a data set, particularly for large sampIe sizes." Could you please show a plot of smoothed scaled Schoenfeld residuals over time with associated confidence bands? (No need to show the tens of thousands of individual residuals.) $\endgroup$
    – EdM
    Nov 14, 2023 at 23:07
  • $\begingroup$ Thank you so much for your answer! I have been looking into that specific chapter you recommend. I have not fully understood it yet. But does this mean that I can ignore the p value for PH, since my data is so large? And I'm not allowed to share any data. Maybe you could tell me what I should look for on the Schoenfeld plot? $\endgroup$
    – Devi Sita
    Nov 15, 2023 at 10:39
  • 1
    $\begingroup$ This would be helped to see the kaplan-meier curves as well... although you're already showing the worst case scenario - whether treatment is better or worse depends on a few things - over what time frame are you considering? I show here that robust error estimation is needed for correct inference and the parameter is a time averaged hazard ratio stats.stackexchange.com/questions/341114/… but you can consider p-rho-gamma testing, or restricted mean survival testing. $\endgroup$
    – AdamO
    Nov 27, 2023 at 18:13
  • $\begingroup$ @AdamO what is exactly p-rho-gamma testing? Is there another name for it? and time-averaged hazard ratio, is that the same as the PH-ratio? Or are we taking about two different things? $\endgroup$
    – Devi Sita
    Nov 29, 2023 at 20:46

1 Answer 1

4
$\begingroup$

With such a large data set, it's possible that the "statistically significant" deviation from proportional hazards (PH) has no practical importance. This is similar to the problems with testing for normality in linear regression in a large data set. With enough data points, even a trivial deviation from the assumption might be "statistically significant."

The first thing to do is to examine the plot of smoothed, scaled Schoenfeld residuals over time. I'd recommend using the original plot functions from the standard R survival package, as at least prior versions of the survminer package had a coding error that made its plots pretty much useless. Therneau and Grambsch explain what to look for in Section 6.5, as explained on this page.

Fundamentally, you want to see how large the variation in the estimated Cox regression coefficient over time, $\hat\beta(t)$ from the smoothed plot, is with respect to the absolute value of the single value, $\hat\beta$, returned by the Cox model under the PH assumption. If the variation is relatively small (admittedly a judgment call), then you can assume that PH holds well enough, explain your reasoning in your report, and show the plot in supplemental information so that reviewers and readers can judge for themselves.

If the variation in $\hat\beta(t)$ is fairly large, particularly if the variation has a systematic trend with time, then either the model is mis-specified in terms of associations of continuous predictors with outcome or the PH assumption doesn't hold very well. If the PH assumption doesn't hold, then there is no single hazard ratio (HR) that can describe the predictor's association with outcome over all time periods. Having a change in HR over time might be clinically meaningful; apply your understanding of the subject matter in that case.

Even if PH doesn't hold, the single coefficient estimate $\hat\beta$ represents a sort of event-weighted average value that might be informative. AdamO shows the rationale on this page, with reference to the literature. Again, whether to accept that value is a judgment call that must be based on your understanding of the subject matter.

Alternatively, survival data that violate PH can sometimes be fit well with a parametric accelerated failure time (AFT) model. You should examine AFT models other than a Weibull model, as a Weibull model also effectively assumes PH. Log-normal or log-logistic models are examples worth examining.

AFT models can be a good choice, as some find AFT models to have a simpler intuitive interpretation than PH models. In AFT models, predictors effectively speed up or slow down the time scale. That's pretty easy to think about. The concept of hazard that underlies PH models can be trickier to grasp and often confuses those who are new to survival analysis. See this recent question, for example.

In response to new graph showing smoothed residuals

There is an extremely large drop-off in the time-varying coefficient up to about 10 years. It's far greater than the associated confidence intervals. The initial coefficient value of ~1.4 is equivalent to a hazard ratio of 4; the value of ~0.4 is equivalent to a hazard ratio of 1.5. It's reasonably constant thereafter.

The violation of PH for this variable is thus fairly substantial, not just representing your very large number of observations. There will be no single hazard ratio that applies over all times. The single coefficient reported by coxph() will provide an event-averaged estimate, if you wish some single value, but (in any event) follow the recommendation for robust error estimates made by @AdamO in a comment on your question.

I would recommend reporting your Cox model and then showing this graph to show how the estimated Cox regression coefficient (and thus the hazard ratio) changes over time. That's a direct way to show what you have found.

$\endgroup$
5
  • $\begingroup$ @SAphi11 look for a reasonably flat smoothed curve that doesn't deviate much from the coxph()coefficient value. Set resid=FALSE in the call to plot() on the cox.zph() output to remove individual residuals. You can show plots of Schoenfeld residuals without violating confidentiality: when you use plot() on the cox.zph() output, include a value for the ylab argument that doesn't include the actual name of the variable. There is then no way for anyone who sees the plot to figure out what your actual data and variable names are. $\endgroup$
    – EdM
    Nov 27, 2023 at 16:40
  • 1
    $\begingroup$ @SAphi11 I added a few more thoughts. You have a real, but potentially interesting and important, change in hazard ratio over time. Apply your understanding of the subject matter to try to come up with an explanation. $\endgroup$
    – EdM
    Nov 27, 2023 at 18:54
  • $\begingroup$ I really appreciate your time! $\endgroup$
    – Devi Sita
    Nov 27, 2023 at 19:03
  • $\begingroup$ @SAphi11 you don't ignore the violation of PH. You report it and illustrate it. The Cox model is still informative despite the violation of PH: this plot illustrates the model's estimate of how the Cox regression coefficient changes over time, as explained here. As AdamO notes in a comment on your question and in a link from this answer, the Cox coefficient and associated hazard ratio reported by the original model is a type of time-averaged value if you want a single value. $\endgroup$
    – EdM
    Nov 27, 2023 at 19:05
  • $\begingroup$ Let us continue this discussion in chat. $\endgroup$
    – Devi Sita
    Nov 27, 2023 at 19:35

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.