Can it be said, that in case of non-proportional hazards, Cox returns average hazard ratio over the entire period?

Question

I have read two opposite statements about the hazard ratios obtained from the Cox procedure in presence of non-proportional hazards. One is that the HR for a certain covariate, that fails this assumption, is averaged over the entire period.

So, if initially it's 0.8 and then it grows to 2, we can say it's 1.4 in average. Which hides the fact, that the hazards were opposite and swapped.

The others say, that it cannot be simply averaged and instead other measures should be used, like the restricted mean survival time. Can you name any source resolving or giving more insights on that? Or simply explain, how to interpret Cox in case of crossing curves (resulting in swapped hazards) or varying hazards in general?

Answer 1

There is a very large statistics literature on this topic. In the face of non-proportional hazards (non-PH), the hazard ratio is a kind of average, and we know the weights that are used as shown in some of the references and R software appearing here. The weighting is not as simple as what you wrote, but you have the right idea.

When there is strong non-PH such as when curves cross, it is not very fruitful to interpret Cox regression coefficients. And switching to mean survival time brings in other problems and leaves us with a hard-to-interpret result (e.g., the mean life length for subjects who died within 3 years was 2.4) and is too dependent on the choice of time horizon. Instead I recommend flexibly modeling what you don't know. See Section 20.7 of this especially the Breslow et al reference. Fit non-PH directly and get time-specific hazard ratios. You can also get cumulative incidence at a specific time estimated from this extended Cox model.

If the hazard ratios of all the predictors converge to 1.0 as $t \rightarrow \infty$ you should switch to an accelerated failure time model instead. The above link to the RMS course notes goes into detail about that.

Answer 2

According to

Schemper, Michael, Samo Wakounig, and Georg Heinze. "The estimation of average hazard ratios by weighted Cox regression." Statistics in medicine 28.19 (2009): 2473-2489

the parameter estimated by Cox regression model is
$$\frac{\int \frac{h_1(t)}{h(t)}f(t)dt}{\int \frac{h_0(t)}{h(t)}f(t)dt}$$ where $h_1(t)$ and $h_0(t)$ are the hazard functions in the two groups, $h(t)=h_0(t)+h_1(t)$ and $f(t)$ is the density of events. This is their formula (3) with the weight function $w(t)=1$.

The paper refers to Kalbfleisch, John D., and Ross L. Prentice. "Estimation of the average hazard ratio." Biometrika 68.1 (1981): 105-112.

If that is true, then it is not a weighted average of $h_1(t)/h_0(t)$.

Consider the specific case where the hazards are piecewise constant with $h_1(t)$ equal 0.1 for $t \in [0,3]$ and 0.3 for $t \in (3,6]$ and with $h_0(t)$ equal 0.2 for $t \in [0,6]$ and where everyone without an event by time 6 is censored. Also, assume 1:1 allocation. The probability of an event before time $t$ for $t<3$ is $F(t)=\frac{1-e^{-0.1*t}+1-e^{-0.2*t}}{2}$. The probability of an event before time $t$ for $t \in [3,6]$ is $F(t)=\frac{(1-e^{-0.1*3}e^{-0.3*(t-3)})+1-e^{-0.2*t}}{2}$. The density, $f(t)$, is the derivative of $F(x)$.
The formula above then becomes
$$\frac{ \int_0^{3} \frac{0.1}{0.3}f(t)dt+ \int_3^{6} \frac{0.3}{0.5}f(t)dt} {\int_0^{3} \frac{0.2}{0.3}f(t)dt+ \int_3^{6} \frac{0.2}{0.5}f(t)dt}$$ $$\approx \frac{ \frac{0.1}{0.3} 0.355+\frac{0.3}{0.5}0.344} {\frac{0.2}{0.3} 0.355+\frac{0.2}{0.5}0.344}\approx 0.868$$

I also confirmed by simulation that the average estimated HR is about 0.87 and the exp of the average estimated logHR is about 0.86.

However, the formula in the linked slides is $HR1^{p1}HR2^{p2}$ where p1 and p2 are the proportion of events in the two periods and HR1 and HR2 are the hazard ratios in the two periods and in this case, they calculated it to be 0.86. Both formulas come up with similar answers. So, the question is, after rearranging the formula above with some algebra, does $$HR1^{p1}HR2^{p2}=\frac{ \frac{HR1}{HR1+1} p1+\frac{HR2}{HR2+1} p2} {\frac{1}{HR1+1} p1+\frac{1}{HR2+1}p2}$$ in general? It can't be. So, one or both of them is wrong. There is an outline of the argument on slide 8 and some references. Some approximations are used there.

One thing you could do is set up an example where the two are very different and then use simulation to check which one is closer. I tried the following example: Arm 1: constant hazard rate of 0.05 in first 3 units of time, constant hazard rate of 20 afterwards; Arm 0: constant hazard rate of 1. Everyone without an event censored at time 6. p1=0.545 and p2=0.455. The formulas gave close results (0.76 and 0.85). But, when I simulated data (1000 per arm, 10,000 replications), the mean hazard ratio was 0.277. The formulas are both wrong! In this example, most of the people in Arm 0 have events before time 3, so I am thinking Cox regression puts more weight in the first time period. I don't know how to fix the formulas. So, the paper

Uno, Hajime, et al. "Moving beyond the hazard ratio in quantifying the between-group difference in survival analysis." Journal of clinical Oncology 32.22 (2014): 2380.

is right when it says

"When the PH assumption is violated (ie, the true hazard ratio is changing over time), the parameter actually being estimated by the Cox procedure may not be a meaningful measure of the between-group difference; it is not, for example, simply an average of the true hazard ratio over time.[6]"

Reference [6] is the above referenced paper by Kalbfleisch and Prentice. The notation is too complicated there for me to easily read, so I don't know whether they have a formula there that would be correct for a scenario like this.