# Can one compare Hazard Ratio and Average Hazard Ratio?

I would like to compare three predictors of survival by comparing three cox regression models. One model violates the PH assumption which brought me to weighted cox regression as implemented by coxphw package by Dunkler, Ploner, Schemper and Heinze (2018).

However, I was wondering whether one can compare Hazard Ratios from cox proportional hazard models and average Hazard Ratios from cox weighted hazard models (i.e. "predictor a was associated with a greater risk of ... (HR = ...) than predictor b (AHR = ...").

# Edit:

Following the sugesstion of @EdM, I elaborate my plans.

• The goal is to evaluate the predictive power of three predictors for occurence of a specific disease (i.e. baseline test A, a novel baseline test B and a combination of these C).
• Three different models are to be set up in order to adjust for covariates. I would like to establish a statement like "Predictor C was associated with a greater risk of occurence of the disease than Predictor A and B (HR A = ; HR B = ...)".
• However, one model does not meet the PH assumption (mainly due to predictor B, which is categorical), this led me to weighted cox regression and to the question whether one can just compare averaged HR with HR obtained by proportional hazards cox regression. Or is this like comparing apples and pears?
• Could you please say more about what you mean by "comparing three Cox regression models"? Are you evaluating 3 predictors separately, or are you using the same predictors but including interactions, nonlinear transformations, etc differently among the models? Please provide that information by editing the question, as comments are easy to overlook and can get deleted.
– EdM
Jul 6, 2021 at 16:46
• @EdM, I edited my post, hope this can shed light on my question Jul 7, 2021 at 8:17

This type of model comparison might actually be leading to your apparent loss of proportional hazards (PH).

There are problems when you omit any outcome-associated predictor from a Cox model, as there are in logistic regression. At the least, that can lead to bias in the coefficient estimates of the retained predictors. The situation can get exacerbated in Cox models evaluating events over time, particularly if predictors are correlated.

Say, for example, that your predictor A is associated with a faster development of disease than is predictor B while A and B are associated with each other. Then if you omit A from a model that still contains B, those with the A marker will preferentially develop disease at earlier times and be lost from the risk set at later times. So in a model that omits A it might seem that the association of B with outcome decreases over time as the higher-risk cases (marked by A) leave the risk set: an apparent loss of PH. If you keep A in your model, the model can take the A status into account and you might find that the combination of A and B maintains PH adequately.

In general, the more outcome-associated predictors you include in your model without overfitting, the better predictive performance will be. As your predictor C is the combination of A and B, one would expect the C model to be superior.

Comparing hazard ratios (HRs), however, is not the best way to compare models. HRs tell you how quickly the disease develops on average, but they don't necessarily demonstrate the model's overall ability to discriminate event times among cases or to reliably gauge the probability of disease development by any given time. You might, for example, have a predictor that works very well with a large HR on a small subset of cases but doesn't work well over the entire population of interest.

Provided that the models are fit to the same data set, overall model comparison can be done by anova() for nested models, or by the Akaike Information Criterion (AIC). There are tools for evaluating Cox model discrimination and calibration in, for example, the rms package in R.

If PH doesn't hold, then there aren't well-defined HRs in any event--they necessarily change over time. Without PH, comparisons among HRs thus make even less sense. The best way to proceed to to grapple directly with PH, as the lack of proportionality might be of substantive interest. For example, use time-varying coefficients, as explained in the R survival package time-dependence vignette. Something as simple as a step function could both solve the PH problem and indicate different associations of a marker with earlier versus later events. Or try an accelerated failure time model that doesn't follow PH, like a log-normal or log-logistic model.

If you nevertheless still want to compare HRs among Cox models that don't meet PH, then you have to choose what type of "average" HR you want to compare. Absent PH, the standard Cox model provides a type of event-averaged HR. That average thus depends on the distributions of event and censoring times. The weighting procedures in coxphw try to to provide more of a time-averaged HR less dependent on the specific event and censoring times. According to the article describing the package, the procedures

weight the contributions of each event time to the partial likelihood according to the product of the survivor function and the inverse cumulative probability of follow-up at that time.

Absent PH, however, those HR values from coxphw aren't time-constant HRs in the sense envisioned in a Cox model. Those HR values can be interpreted in terms of odds of concordance, a useful measure of a single model's performance, but concordance is not sensitive for distinguishing among models. That continues to argue against using HR values to compare among models.

Absent PH, comparing HRs between standard Cox models and those from coxphw is "comparing apples to pears," as you fear. So choose one or the other type of model for such comparisons. If all models are fit on the same data, my sense is that the coxphw correction will be similar for all models, so that inter-model comparisons should be similar with either standard or coxphw models even if the absolute HR values differ between those model types. One might argue that the coxphw approach for all models will be the most generally reliable, but I have no practical experience with it.

• Thanks for your elaborate answer, @EdM. Since the violating variable does not change over time bc it refers to a baseline category (to be more exact, the df consists of one row per subject, i.e. one observation per subject) I applied a step function as described in the vignette by Therneau et al. This yielded convergence warnings, huge coefficients and NA values. Applying accelerated failure time models (weibull dist) worked out but gave out higher AIC values in all 3 models than in Cox models (dif AIC ca. 80). Would you suggest weighted cox regression for all is therefore the best choice? Jul 8, 2021 at 9:35
• @a.henrietty Look at the plot of smooth scaled Schoenfeld residuals over time from cox.zph(). Are the deviations from PH large enough to matter? If they are, try to match the shape over time with a step or a tt() function. Convergence warnings for a step model suggest that data weren't reformatted properly with survSplit() to match the step time. Check that. I'm not sure that AIC values are comparable between Cox and AFT models. Examine model validation and calibration. I see weighted Cox as a last resort.
– EdM
Jul 8, 2021 at 12:24
• @a.henrietty Just noticed that you tried a Weibull for your AFT model. Weibull is the only distribution family that is both AFT and PH. So if PH isn't working for your Cox model it won't work for a Weibull parametric model either. Try non-PH AFT models like log-normal or log-logistic. For AIC comparison, see my answer to your related question. I suspect that the AIC values from partial-likelihood Cox models can't be compared directly against those from full-likelihood parametric models, but I'm not completely sure.
– EdM
Jul 8, 2021 at 14:19