When should the robust error function in R be used in Cox proportional hazard regressions?

Let's consider a scenario where I compute a HR that appears to violate the proportional hazard assumption according to the cox.zph test, but when plotting the smoothed scaled Schoenfeld residual plot, the violation is found to be not "real", as determined by examining the residual plot to assess the deviation of the coefficient through time, and in this case, the deviation is minor.

I’ve been studying this issue based on a previous Cross Validated discussion in this link. Despite the assumption being violated, the hazard ratio (HR) can still be used for interpretation as long as the function robust=T is applied in the Cox model. Then, the HR will be referred to as the failure-time-averaged hazard ratio.

However, while I'm not a statistician, I've attempted to understand the concept of when to use robust=T.

My questions are as follows:

1. From my understanding, robust errors are employed to calculate standard deviations in heteroscedastic data. Therefore, if the proportional hazard assumption is violated, is it then assumed that the data used for the Cox analysis is heteroscedastic? Is that why it's recommended to calculate the robust error in this case? In other words, does the violation of the proportional hazard assumption indicate heteroscedasticity, or does it imply it for certain?

2. Should the robust=T function always be included in a cox model, even when the proportional hazard assumption is not violated?

3. Should the time-averaged HR always be referred to as the “failure-time-averaged hazard ratio”, or only when the robust=T function is implied in the Cox model?

Question 1. There are several types of robust standard errors that can provide useful estimates of the (co)variances of regression coefficient estimates when model assumptions aren't met. Besides the methods for dealing with heteroscedasticity in linear regression, there are forms of robust standard errors that deal with autocorrelation and clustering; see the Wikipedia page on clustered standard errors for some links. That said, I'm not sure how to define "heteroscedasticity" in a survival model.

Question 2. The sandwich estimator for Cox models is "robust to several possible misspecifications in the Cox model including the lack of proportional hazards, incorrect functional form for the covariates, and omitted covariates." (Therneau and Grambsch, p. 160). A potential problem with that robust sandwich estimator is its own variability in small data sets. Nevertheless, Therneau and Grambsch conclude (p. 161):

we still recommend examination of the sandwich estimate in cases where the Cox model assumptions are suspect.

Question 3. If the proportional hazards assumption isn't met, then any coefficient estimate will necessarily be an event-averaged estimate. The robust estimator does not change the coefficient estimates themselves. It provides an estimate of the variance-covariance matrix of the coefficient estimates that can be used more reliably for inference.

• thanks for your answer! Q2: I have a large dataset, more than 100.000 observations, so I applied the robust=T to all my cox models to be on the "safe side", but is it incorrect to think like that? Because I do know that the ph assumption can be violated by the smallest deviations in large sample sizes. Some of my cox models met the ph assumption, and some did not. For those models that didn't violate the ph assumption, should they be computed without the robust=T then? Commented Apr 29 at 18:03
• and for Q3: I thought that the HR should be referred to as the event-averaged HR when the robust=T is included in the Cox model. However, I now understand that when the proportional hazards (ph) assumption isn't met, it is referred to as an event-averaged HR, regardless of whether robust=T is included or not. This may sound trivial, but what is a HR that does not violate the ph assumption referred to as, if not event-averaged HR? @EdM Commented Apr 29 at 18:08
• @Aphi11 robust=T also helps deal with omitted covariates and improper functional forms for associations of continuous predictors with outcome. It thus can be a good, cautious choice, particularly with a very large data set in which the variability of the robust estimator (seen in small data sets) shouldn't be a problem. If PH isn't violated, then a hazard ratio is conventionally just called a hazard ratio. That's not really any different from a linear regression slope, which is effectively a type of average over all the observations, when linearity holds.
– EdM
Commented Apr 29 at 18:17
• Thank you for clarifying things! @EdM Commented Apr 29 at 18:19
• @Aphi11 one caution: a "large" data set for survival analysis mainly is a function of the number of events, not the total number of observations. I'm assuming that a reasonable fraction of your 100,000 observations are of event times.
– EdM
Commented Apr 29 at 18:20