# Recommendations for Non-Proportional Hazards

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.

• Have you tried adjusting for clinical/demographic variables with propensity scores instead of Cox regression? With this wealth of data propensity-score matching might be feasible.
– EdM
Mar 22, 2017 at 19:39
• @EdM Not for these data. My uncertainty regarding exactly how to propensity score match multiple category data (i.e., >2 categories) has prevented me from trying this method. However, in my experience, the results of the multivariable analysis will be very similar the results of the propensity score matched analysis (given that the purpose of both is to address selection bias). So, I would suspect that I'd end up with the same issue of non-proportionality. Mar 22, 2017 at 19:50
• That is for confounding not outcome heterogeneity. Mar 22, 2017 at 19:50

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

• Thanks for your insight; I hadn't considered parametric models. I do most of my modeling in SAS, so your course notes showing how to estimate the log-normal model using the rms package were incredibly useful. Although the fit of the log-normal model isn't ideal, I believe it's adequate and sufficient to report (Residuals). The PI will inevitably need to be reminded that a single study is by no means definitive. Thanks again. Mar 23, 2017 at 19:29

You certainly don't have marginal proportional hazards. That does not mean you don't have conditional proportional hazards!

To explain in more depth, consider the following situation: let's suppose we have group 1, which is very homogeneous and has constant hazard = 1. Now in group two, we have a heterogeneous population; 50% are at lower risk than group 1 (hazard = 0.5) and the rest are at higher risk than group 1 (hazard = 3). Clearly, if we knew whether everyone in group 2 was a higher or lower risk subject, then everyone would have proportional hazards. This is the conditional hazards.

But let's suppose we don't know (or ignore) whether someone in group 2 is at high or low risk. Then the marginal distribution for them is that of a mixture model: 50% chance they have hazard = 0.5, 50% they have hazard = 3. Below, I provide some R-code along with a plot of the two hazards.

# Function for computing the hazards from
# a 50/50 heterogenious population
mix_hazard <- function(x, hzd1 = 0.5, hzd2 = 3){
x_dens <- 0.5 * dexp(x, hzd1) + 0.5 * dexp(x, hzd2)
x_s    <- 1 - ( 0.5 * pexp(x, hzd1) + 0.5 * pexp(x, hzd2))
hzd    <- x_dens/x_s
return(hzd)
}

x <- 0:100/20
plot(x, mix_hazard(x),
type = 'l',
col = 'purple', ylim = c(0, 2),
xlab = 'Time',
ylab = 'Hazard',
lwd = 2)
lines(x, rep(1, length(x)), col = 'red', lwd = 2)

legend('topright',
legend = c('Homogeneous',
'Heterogeneous'),
lwd = 2,
col = c('red', 'purple'))


We see clearly non-proportional marginal hazards! But note that if we knew whether the subjects in group 2 were high risk or low risk subjects, we would have proportional hazards.

So how does this affect you? Well, you mentioned you have a lot of other covariates about these subjects. It is very possible that when we ignore these covariates, the hazards are non-proportional, but after adjusting for them, you may capture the causes of the heterogeneity in the different groups, and fix up your non-proportional hazards issue.

• Thank you for the reply! Your point is well-taken as the missing covariate problem can manifest as non-proportional hazards. I forgot to mention that, although the log-log survival curves were unadjusted, I tested non-proportionality using interactions with time as well as Schoenfield residuals after adjusting for all covariates of interest. I have edited my post to reflect this. Mar 22, 2017 at 18:58