Estimating robust variance without cluster id in survival model

Question

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.

Answer 1

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$.

Answer 2

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.