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

Question

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.

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

3. 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.

Answer 1

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.

Answer 2

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.