# Which method of handling non-proportional hazards (here: crossed KM curves in 2+ points) would you prefer?

I have an experiment, where I compare two groups over time in frames of the time to event analysis. The Kaplan-Meier curves cross 3 times.

It means, that the hazards ratios switch to the opposite side multiple times. So, the problem is not even in non-constant HR, but even worse - the interpretation varies from period to period. Like on the picture below.

I am interested in 2 analyses: the Cox regression (with just single 2-level categorical group) to get the HR and the CI for it and the curve comparison test.

Normally, I would run the Cox model and get the HR and the p-value (equivalent to long-rank in this case).

But here I have multiple HRs. What should I do?

1. Run separate Cox / KM analysis in each period determined by the moment of crossing? But within each period the curves do not diverge from each other, as the PH suggests, but rather first diverge and then converge (to cross each other). So the PH assumption is still violated!

I could "manually" find the moments of crossing, R has such function and split the analysis by those periods.

1. employ time dependent covariate? I mean group x time interaction? But how to interpret such model?

2. regarding the statistical test, I can imagine a test, that is not fooled by the turning hazards ratios and reporting that "in overall the two curves differ, regardless of the nature of the difference". But what kind of test can do that?

Is there any other, better method to handle such analysis? I know the AFT model doesn't need the PH assumption, but still - how will it handle the "race" between the two curves?

Please note, I am aware, that missing covariate could cause that, but I have no additional information in my data set. Only the basic survival data I need to summarize somehow.

Let's assume the difference reported naively with the log rank is stat. significant. So saying "if there is no significant difference, your problem diminishes". The curves diverge far from each other to then meet each other and cross.

• If these are survival curves, then your annotations with respect to HR aren't correct. The hazard is the ratio of the slope of a survival curve to the current survival value. One survival curve higher than another over a period of time doesn't mean that its HR is lower over that period of tine. In fact, at the end of period 1 where the overall survival is the same for both groups, the steeper slope for the blue group means that its hazard it greater than that for the red group at that point in time.
– EdM
Mar 26, 2021 at 17:54
• For each of the first 3 periods of time, overall survival is the same for both groups at each endpoint. Thus for each of those 3 periods one definition of an average hazard--the integral of the hazard over time (cumulative hazard within the period), divided by the elapsed time--is necessarily the same for the 2 groups. See the Wikipedia entry for example. From your curves, a difference in cumulative or average hazard (as defined above) is only seen in period 4.
– EdM
Mar 26, 2021 at 18:07

In small samples under the null hypothesis, we get curves that cross often; this is just an artifact of randomness, and a null result is expected. In large samples, however, this can be challenging because a nuance is added to the initial hypothesis, and it is dangerous to change your hypothesis based on observed data.

In textbooks, they will tell you non-proportional hazards reduce your power. That may be true if the hypothesis is whether the survival curves differ (at all). But that's a misinterpretation of the log rank test: the test is actually whether the hazard ratio is different from 1, just like in a Cox model. In fact, I think it's rather unimportant to know whether the survival curves differ on the whole, they almost always do.

A time-by-covariate interaction in a Cox model is a special case of an accelerated failure time model. If the crossing curves are because the hazard ratio varies by time and by covariate, this can solve the muddle, but you have to be either very lucky, or actually understand the science of what's happening. The competing risks package has some implementation for these models. And of course, interpreting an interaction is very difficult in any linear model, let along a Cox model.

There are 3 solutions I can think of:

1. Use a robust error estimate. A surprising finding is that when the hazard ratio is non-proportional, the Cox model estimates a failure-time weighted average hazard ratio over the duration of follow-up.

2. Use a $$g\rho \gamma$$ family of estimators: these can weight earlier or later survival if those are of greater interest and it's simply the tails of the analysis where high variability is causing curves to cross. For instance, if a cancer kills all patients within 5 years, but the median survival is 2 on drug A and 1 on drug B, we might value those first two years of survival greatly

3. Use a restricted mean estimate of survival. This is the only estimator that estimates average survival time from censored data. It is just the area under the Kaplan-Meier curve. The "restriction" is just the duration of follow-up that should be curtailed to an appropriate time frame as with approach 2, so the area under the first two years, or four years, will give you an appropriate estimand and test that does not rely on a proportional hazard assumption.

• Thank you very much for your answer. But shouldn't I focus on the separate periods, since the HR get opposite? We do not have just "one kind of differences between the two phenomena", because one time the first curve "is above the other" (lower risk), and then the situation turns opposite. It tells me,that, over time, the nature of the phenomenon varies. Maybe I shouldn't ignore this and pretend these are just "fluctuations" over time (that's what we do by using a single method, ignoring this switching) but rather split it and analyse separately? And, within period,I could use a robust method? Mar 24, 2021 at 16:21
• No you shouldn't focus on the periods defined by where the KM curve crossed because those are incidental findings, even if the true hazards cross at least once, you would need to define some estimator for the crossing point, or for the total number of crossings. The required sample size to estimate these things reliably are way out of scope for any standard analysis. Any effort to combine results over the different periods would give you an estimate similar to the overall estimate. Mar 24, 2021 at 16:58
• I run the analysis and got very different HRs per each time point, so I think I will have to end up with the time varying covariates. My client won't accept that so highly crossing curves are unimportant (highly I mean the difference in HR is really big, like 0.3 to 2.5). In the medical area, where I work, it's common that over time the nature of the phenomenon changes and one method may get much worse than the other. I have no other information about the process that could justify it, but I have to handle it and report separate HRs anyway. It comes just from the domain knowledge. Thank you. Mar 24, 2021 at 17:16
• Yes, the estimator would tell me if this is statistically detectable. This is problematic, as I this is not a designed confirmatory experiment, where the statistical inference makes sense. It's rather ad hoc, undesigned observatory experiment, where statistical inference has very poor justification. That's why there only the descriptive analysis makes sense, and any kind of inference is dangerous. That's why I thought I could use the observed crossing points. Something like "turning/changing point analysis". I'm now searching for R/SAS packages/macros. Maybe this will help there. Mar 24, 2021 at 17:19

Consider whether a proportional hazards model is appropriate for your data in the first place.

To my eye, the red curve is very similar to the blue curve except that it is stretched out in time. Each of the inflection points on the red curve happens later than on the blue curve. If those are survival curves, then that's just what you would expect with an accelerated failure time (AFT) situation: time with respect to experiencing this event goes "slower" for those in the red group than for those in the blue group. The "race" that you see with this plot shown as a function of absolute time might disappear when evaluated with time axes differentially stretched for the 2 groups, as in AFT. The difficulty is that the survival functions don't immediately suggest any single parametric form for an AFT model that can capture what seem to be four biologically distinct stages with respect to survival.

If these curves represent your data and you have an adequately sized study, apply your understanding of the subject matter. What underlying processes might distinguish those 4 separate time periods: first a very slow event rate, then highly accelerated, then slow again, finally accelerated again. That functional interpretation of your data might be much more important than trying to force your data into a proportional hazard framework.