Would you used prediction from Cox or stratified KM for drawing survival curves and testing them?

Question

Let's assume I have an experiment, where I track some event. I want to analyse the time to this event between two levels of a group G: g1 and g2

I want to do the comparative analysis in sub groups, for example:

• All
• age: <50yrs, >= 50yrs
• sex: male, female
• GroupA: a1, a2

crossed with GroupB: b1, b2, All

What I need:

1. KM curves. I want to plot them using ggplot2 as a grid of plots, like below: g1 = red, g2 = blue

2. measure of effect of difference between them. It will be HR and the restricted mean survival time (a standard where I work)

Of course for the non-proportional hazards the HR will be somehow "averaged" (no single HR over the entire trial, it varies) and probably meaningless in certain situations, but let's ignore it now.

3. test of significance of the difference. Let's focus now on the HR, assuming it makes sense. I want the log-rank or any other weighted test, depending on late, early, diminishing effects, but for now let's assume it's the classic log-rank.

I do not know nothing about the proportionality of hazards, cannot anticipate it.

Now, my questions:

1. To draw the curves, I can use either:

• the Kaplan-Meier stratified by the G={g1, g2}, which is a non-parametric estimator and does not assume proportional hazards equal baseline hazards

• prediction taken from Cox regression, where the G={g1, g2} is taken as a covariate. Cox will force the baseline hazard to be equal for both curves, so it makes stronger assumption than the K-M.

I already saw, that KM and Cox can produce different curves: Why do my survival curves generated by the Cox differ from Kaplan-Meier for the simplest model?

Which one would you suggest for drawing? Stratified KM vs. Cox? I think the curves should be maximally real, so no fake assumptions should be made, so the KM should be drawn. Would you agree?

2. To report the PH, I have to make the assumption on equal baseline hazards. Otherwise I cannot use Cox and calculate them.

If I use stratified KM + Cox-originated PH, I may get a discrepancy, but there's some cost to pay I'm afraid.

3. To compare the curves and assess significance of the PH, I need a test. I have two options:

• Mantel-Hanshel log-rank from the Kaplan-Meier, which does no assume equal baseline hazards
• log-rank score test from the Cox model.

I was told, that the log-rank is a special case of the Cox one, so they should provide exactly the same result. But I noticed, and not once, they differed. I'm wondering is this because the Cox assumes baseline hazards to be equal, and the Mantel-Hanshel does not? If so, why so many people say that the log-rank can be ignored and not taught any more (for example Professor Harrell; I saw it on this forum), because we have Cox, if Cox puts stronger assumptions and there are potentially differences?

TL;DR

To draw curves, report PH and test them, would you use the predictions from Cox, or stratified KM for curves, and the rest from Cox, even if there are potential discrepancies?

Answer 1

With Kaplan-Meier curves, each time you try to take another covariate into account by subsetting and stratifying, you lose events and power. That's one reason why so much emphasis is placed on survival models that include multiple predictors. That way you share information among all combinations of predictor values to get (1) higher-precision estimates of how the predictors together are associated with outcome and (2) reliable tests of significance--provided that the model assumptions hold.

I do not know nothing about the proportionality of hazards, cannot anticipate it.

Then let the data tell you. A useful approach is to start with a Cox proportional hazards (PH) model, and test whether PH holds adequately. That's typically done with smoothed plots of scaled Schoenfeld residuals for each predictor over time, for example with the cox.zph() function in R.

If PH holds reasonably well, and your model is well calibrated and not overfit, then all of the comparisons and tests that you want to perform are straightforward and much more powerful than the corresponding subset-based Kaplan-Meier log-rank comparisons would be. If PH don't seem to hold at first, there are several ways to adapt models to meet the assumption: e.g., transforming continuous predictors (to help meet the linearity assumption, too), stratifying by predictors that don't meet PH but aren't of primary interest, and modeling with time-dependent coefficients.

If you can't get PH to work, there are other types of models. Accelerated failure time (AFT) models sometimes work when PH doesn't hold, and if a fully parametric AFT model (e.g., from the survreg() function in R) fits well enough you can contemplate reporting mean survival rather than restricted mean survival. Or you can try an additive hazards model (aareg() in R) the doesn't depend on PH, allowing for time-varying coefficients of predictors as it estimates the baseline hazard.

The take-away here is that modeling done properly allows you to share information among cases in a way that provides superior power over the subsetting and stratifying needed to generate Kaplan-Meier curves. Cox PH models have advantages, but you can still get the power that comes from modeling with other approaches if you can't make PH work well enough.