I built a survival model with coxph. I have a large dataset (1 000 000 instances) and the global model is not proportional. It cannot be improved by introducing time-dependent coefficients or stratifying. My understanding is that estimating robust standard errors is in order in such a case. Coxph allows estimating robust standard error with robust = TRUE, but in the survival package documentation, it is always used together with cluster id option. I do not think that using clustering variables exists in my data. Does it make sense to estimate robust errors without cluster id in coxph? I would appreciate any help.

  • 1
    $\begingroup$ With such a large data set you are very likely to have proportional hazard (PH) tests fail, just as normality tests tend to fail with large data sets. Robust errors are fine, but please evaluate whether the deviation from PH is large enough to be important in practice. $\endgroup$
    – EdM
    Jul 29, 2021 at 15:36
  • $\begingroup$ Thank you. It is true, about large data, but I have never seen the theory about a threshold of ph deviation being too large or not, so I would not know how to assess it. $\endgroup$
    – Milo
    Jul 29, 2021 at 19:19
  • $\begingroup$ That’s not a statistical question about a “threshold.” It’s a practical consideration informed by your understanding of the subject matter. Examine, say, a loess-smoothed plot of scaled Schoenfeld residuals versus time, as you can get with the cox.zph function in R. Even if it’s not completely flat, the deviations might not be very large in comparison to the overall level of the hazard ratio/regression coefficient. $\endgroup$
    – EdM
    Jul 29, 2021 at 20:27
  • $\begingroup$ Thank you EdM@, it is very helpful. $\endgroup$
    – Milo
    Jul 30, 2021 at 8:56

2 Answers 2


The coxph model used to have the option robust=TRUE, it may still in fact, or it has been deprecated. This is because there's a connection between the robust error estimator and the frailty methods for repeated measures. The documentation says the following:

The cluster term is used to compute a robust variance for the model. The term + cluster(id) where each value of id is unique is equivalent to specifying the robust=TRUE argument.

So just create a unique row identifier and use it as a clustering variable to get robust variance. One way (in R) is to say mydat$id <- 1:nrow(mydat).

As the question: does it make sense? Sure! A parallel can be drawn to, say, comparing the generalized estimating equation with identity link. In the GEE one must specify the correlation structure as exchangeable or independent. The exchangeable structure has the residual-covariance for a particular cluster $\Sigma_{exch} = (1-\rho)I + \rho J$ where $I$ is the identity matrix, $J$ is the all-one matrix, and $\rho$ is the correlation. The identity structure is simply $\Sigma_{ident} = I$. The two are identical when the cluster size is 1.

The generalized estimating equation takes the form $0 = \sum_{i=1}^n D^T V^{-1} X^T\left( Y - g(X\beta)\right)$. The $A$ bread matrix is $\partial D(\mu) / \partial \beta$ and gives the usual model-based error and is related to the expected information. The $B$ meat matrix is $X^T\left( Y - g(X\beta)\right) \left( Y - g(X\beta)\right) ^T X / n$ and is related to the observed information. (The book by Stefanski and Boos will give a thorough treatment on the subject). But this is all a long way of saying that if you have an exchangeable correlation matrix with clusters of size 1, it will give the same result as with an independence correlation structure, that is the sandwich error estimate $\text{var}({\beta}) = A^{-1}B(A^{-1})^T$.

  • $\begingroup$ Thank you very much, Adam. This is very helpful. Coxph estimates the robust variance without including cluster-id at all. I assume it is equivalent to including cluster-id unique for every instance. Would you be kind enough to explain the last sentence in simpler terms? $\endgroup$
    – Milo
    Jul 29, 2021 at 15:39
  • $\begingroup$ @Milo hmm, sure, but what is the confusion? even I'm a little unclear how exactly the sandwich is calculated for a cox model specifically. Are you aware of the sandwich? Are you aware of the A, B matrices? I can try to link some other SE posts for background. $\endgroup$
    – AdamO
    Jul 29, 2021 at 15:51
  • $\begingroup$ I do not know A and B matrices and the term exchangeable correlation structure. I know that the sandwich estimator uses square residuals. This is probably not the subject that can be explained in one go. Thank you a lot. $\endgroup$
    – Milo
    Jul 29, 2021 at 20:35
  • $\begingroup$ @Milo see if the edit clarifies things at all. $\endgroup$
    – AdamO
    Jul 29, 2021 at 21:24
  • 2
    $\begingroup$ robust=TRUE is not deprecated. And it's a good thing that it isn't, because the robust variance estimator isn't the same thing as a frailty estimator -- not only are the assumptions and variances different, even the point estimates are different. $\endgroup$ Jul 29, 2021 at 23:01

Insofar as the problem is calculating a robust variance, the answer from AdamO (+1) solves the problem. You can calculate a robust variance estimate for a Cox model with only 1 data point per individual, if you specify robust = TRUE or include a cluster(id) term in the model with id values specified for the individuals.

For a Cox model, Therneau and Grambsch explain in Chapter 7 how the matrix product $D'D$ provides a robust "sandwich"-type estimator of the coefficient covariance matrix, where $D$ is the matrix of "dfbeta" residuals (score residuals scaled by the covariance of coefficient estimates). The "robust" variance is an adjustment like that with generalized estimating equations (GEE). (As Thomas Lumley notes in a comment, that's different from a frailty(id) term, which is like a mixed model.) You can do that directly with your data set.

The question is whether a robust variance estimate is what you should focus on here. Yes, as AdamO explains here, with a robust variance estimate you can do inference on the Cox model coefficients when the proportional hazards (PH) assumption is violated, and those coefficients (whose point estimates will be the same as in the usual Cox model) can be thought of as an average over all events--depending on the data, perhaps close to a time average.

There might, however, be even more to learn from investigating the reason for the apparent failure of PH.

First, it's possible that the massive scale of your data set gives a "statistically significant" violation of PH that doesn't really matter in practice. That's like the case with normality testing of large data sets--as hazards are never perfectly proportional and real data are never exactly distributed normally, large data sets can lead to "statistically significant violations" of model assumptions that have no practical significance. Look at the plots of scaled Schoenfeld residuals over time to see just how big you PH problem is. AdamO discusses the problem of PH in large data sets here.

Second, if the PH violation is large enough to be practically important, that might be a sign that something is wrong with your model. Therneau and Grambsch say on page 148:

There are several possible model failures that may appear as time-varying coefficients but would be dealt with more fruitfully by another approach. These include the omission of an important covariate, incorrect functional form for a covariate, and the use of proportional hazards when a different survival model is appropriate.

Omission of an important covariate will lead to PH violation even if the omitted covariate is balanced between groups. As the at-risk population changes over time with respect to the omitted covariate, "the operative hazard is the average hazard of those at risk at each time point, a mixture of hazards" (Therneau and Grambsch, page 150). There is also a bias toward 0 in the coefficient estimates of the included predictors.

On page 150, Therneau and Grambsch also say "a finding of lack of proportionality should lead one to check functional form diagnostics." Something as simple as a log transformation of a predictor might provide adequate PH.

It's also possible that the underlying survival phenomenon doesn't follow PH. As Therneau and Grambsch note, an accelerated failure time (AFT) or an Aalen additive model might be called for instead. (A Weibull model assumes both AFT and PH, so use a different AFT distribution if you want to deal with a PH violation.)

With this large a data set, you presumably have an opportunity to evaluate such possibilities in some detail. You might not have to limit yourself to accepting the violation of PH.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.