I recall having heard that the hazard ratio, estimated in a Cox model, can be made robust against the parallel hazard functions assumption. The key to this is using a Huber-White, or Huber-Eicker-White, or Heteroscedasticity Consistent (HC) "sandwich" standard error. In this case the hazard ratio is interpreted as a failure-time-averaged hazard ratio.
Is this true? Is there a citation which I can refer to to claim this?
Apparently, this work came about through two different papers. First it was Struthers and Kalbfleisch who in 1986 published "Misspecified proportional hazard models" in Biometrika. They pointed out parameters from misspecified Cox models have a useful interpretation, but one that is defined implicitly. They use the example of accelerated failure time models as a probability model, an example where proportional hazard assumption would be violated.
The implicit expression is that the consistent estimator is the solution to $h(\beta) = 0$ where:
$$ h(\beta) = \int_{0}^{\infty} s^{(1)}(x) dx - \int_{0}^{\infty} \frac{s^{(1)}(\beta, x)}{s^{(0)}(\beta, x)} s^{(0)}(x) dx$$
and $S^{(j)}(x) = n^{-1} \sum_{i=1}^n Z_i^j Y_i(x) \lambda(x; Z_i)$ is the expected survival for a person having a covariate $Z$ differing by $j$ units. $s^{(j)}(x) = E \{ S^{(j)}(x)\}$ where the expectation is taken over the failure times.
A way they interpret that summary is as a time-averaged hazard ratio.
It is then Lin and Wei who in 1989 discuss the sandwich covariance estimator $\hat{A}(\beta)^{-1} \hat{B}(\beta) \hat{A}(\beta)^{-1}$ where $\hat{A}(\beta) = -n^{-1}\frac{\partial^2}{\partial \hat{\beta}^2} L(\hat{\beta})$ is the Hessian of the partial likelihood function evaluated at the MPLE and $\hat{B} = n^{-1} \sum U_i (\hat{\theta})U_i^{\prime}(\theta)$ is the (biased) MLE of the variance of the score function. They show that the misspecified $\beta$ estimated above is rigorous and it has an asymptotically normal distribution and thus tests and confidence intervals can be constructed about its value.
