This is an issue that has plagued me for a long time and I have found no good answers in textbooks, Google, or Stack Exchange.
I have data set of >100,000 patients for which four treatments are being compared. The research question is whether survival is different between these treatments after adjusting for a bunch of clinical/demographic variables. The unadjusted KM curve is below.
Non-proportional hazards were indicated by every method I used (e.g., unadjusted log-log survival curves as well as interactions with time and the correlation of Schoenfield residuals and ranked survival time, which were based on adjusted Cox PH models). The log-log survival curve is below. As you can see, the form of non-proportionality is a mess. Although none of the two-group comparisons would be too difficult to handle in isolation, the fact that I have six comparisons is really puzzling me. My guess is that I won't be able to handle everything in one model.
I'm looking for recommendations on what to do with these data. Modeling these effects using an extended Cox model is likely impossible given the number of comparisons and differing forms of non-proportionality. Given that they are interested in treatment differences, an overall stratified model isn't an option because it won't allow me to estimate these differences.
So, feel free to rip me apart, but I was thinking about initially estimating a stratified model to get the effects of the other covariates (testing the no-interaction assumption, of course), and then re-estimating separate multivariable Cox models for each two-group comparison (so, 6 total models). This way, I can address the form of non-proportionality for each two-group comparison and get a less wrong estimated HRs. I understand that the standard errors would be biased, but given the sample size, everything will likely be "statistically" significant.
Best Answer
Fantastic question fantastic answers. I'll add that you should consider a model making much different assumptions such as the lognormal survival model. Use the normal inverse function for the y_axis instead of log-log. Still need to covariate adjust. So also look at normality of residuals stratified by treatment. This is covered in a case study near the end of my course notes at http://biostat.mc.vanderbilt.edu/rms