Why Account for Clusters?
A cause-specific Cox model can be fitted to subjects who belong to identifiable groups, such as patients treated at the same hospital, members of the same family, or repeated records from the same person. Outcomes within a group may be associated. If the analysis treats every subject as independent, its usual model-based standard error can be misleading, even when the fitted coefficient is useful. A cluster-robust variance estimate adjusts the estimated uncertainty by using the variation in score contributions across independent clusters.
The adjustment concerns variance estimation: it does not, by itself, change the fitted regression coefficient, the risk sets, or the cause-specific hazard being modeled. Its justification depends on having independent clusters and a suitably specified estimating model. It is not a general remedy for a misspecified hazard model, dependent clusters, or too few clusters.
The Sandwich Estimate
For cause \(k\), a subject remains in the cause-specific risk set until the first of the target event, a competing event, or loss to follow-up. A competing event therefore removes that subject from the risk set at its occurrence; it is not counted as a cause-\(k\) failure. This construction defines the cause-specific hazard and does not require treating competing event times as independent latent failure times. Valid estimation still requires appropriate handling of loss to follow-up, commonly conditional independence of censoring and the event process given the modeled covariates.
Suppose there are \(G\) independent clusters. Let \(U_g(\beta)\) be the contribution of cluster \(g\) to the Cox partial-likelihood score, and let \(\hat\beta\) solve the total score equation. With distinct event times, the score can be written as a sum over cause-\(k\) failures of the observed covariate vector minus its risk-set weighted mean. At a failure time \(t\), that weighted mean is
where \(R(t)\) contains those still event-free and under observation just before \(t\). The score contribution for a failure is \(X_i(t)-\bar X(\beta,t)\). Contributions are allocated to subjects and summed within their clusters. With tied failures, the score contributions must be calculated using the same tie-handling method used to fit the model.
Let \(A(\hat\beta)\) be the observed information, the negative derivative of the total score with respect to \(\beta^\top\), evaluated at \(\hat\beta\). The cluster-robust covariance estimate is
The middle term is often called the meat, and the two inverse-information factors are the bread. Since the total fitted score is zero, the cluster scores satisfy \(\sum_g U_g(\hat\beta)=0\), apart from numerical tolerance. For a single coefficient, the estimated robust standard error is the square root of the resulting scalar variance. For a coefficient \(\beta_j\), exponentiating gives a cause-specific hazard ratio; a Wald interval on the coefficient scale can be exponentiated to obtain an interval for that hazard ratio.
Why the Formula Works
The sandwich form follows from a first-order expansion of the score equation. Write \(\beta_0\) for the value targeted by the correctly specified model. Near \(\beta_0\), the fitted score is approximately the score at \(\beta_0\) minus the information times the estimation error:
If clusters are independent, the covariance of a sum of their scores is the sum of the cluster score covariances. The observed outer products \(U_g(\hat\beta)U_g(\hat\beta)^\top\) estimate those covariances without requiring subjects within a cluster to be independent. Substituting this estimate and the fitted information into the linearized expression gives the sandwich formula. The approximation is an asymptotic argument: it relies on having enough independent clusters and regularity conditions, not merely a large total number of subjects.
Worked Example: One-Coefficient Sandwich Calculation
Suppose a hypothetical cause-specific Cox fit estimates the coefficient for a one-unit increase in a treatment score as \(\hat\beta=0.40\). The observed information is \(A=4.00\). Four clusters have fitted score contributions \(1,-1,2,-2\). These contributions sum to zero, as expected at a fitted solution.
The meat is the sum of squared cluster scores:
Thus the robust variance and standard error are
The model-based standard error would be \(1/\sqrt{4}=0.5000\), under the additional working assumption that subjects contribute independent information. The robust 95% Wald interval for the coefficient is \(0.40\pm1.96(0.7906)=(-1.150,1.950)\). The estimated hazard ratio per one-unit increase is \(\exp(0.40)=1.492\), with robust interval \((\exp(-1.150),\exp(1.950))=(0.317,7.03)\). These data provide very imprecise evidence about the hazard ratio. With only four clusters, the usual large-cluster approximation is especially questionable; the calculation illustrates the formula, not a recommendation to trust this interval.
Worked Example: Two-Coefficient Sandwich Calculation
Consider a hypothetical fitted model with two covariates and observed information matrix \(A=\begin{pmatrix}4&0\\0&3\end{pmatrix}\). The four cluster score vectors, in coefficient order, are \((1,1)^\top\), \((-1,0)^\top\), \((0,-1)^\top\), and \((0,0)^\top\). The vectors sum to \((0,0)^\top\). Their outer products give
The off-diagonal element is \(1\), because the cross-products are \(1,0,0,0\); it is not zero. Since \(A^{-1}=\begin{pmatrix}1/4&0\\0&1/3\end{pmatrix}\), the robust covariance estimate is
The robust standard errors are \(\sqrt{0.1250}=0.3536\) and \(\sqrt{0.2222}=0.4714\). The estimated correlation between the coefficient estimates is \(0.08333/\sqrt{0.1250(0.2222)}=0.5000\). This nonzero covariance matters for joint inference, such as a Wald test of both coefficients; using only the two marginal standard errors would discard information about their estimated dependence.
Worked Example: Comparing Robust and Model-Based Inference
In another hypothetical one-covariate fit, let \(\hat\beta=0.50\), observed information \(A=9\), and six cluster scores be \(1.2,-0.8,0.6,-0.4,0.2,-0.8\). They sum to zero. The meat is
Therefore \(SE_{\mathrm{rob}}=\sqrt{3.28/9^2}=\sqrt{0.04049}=0.2012\). The model-based standard error is \(1/\sqrt{9}=0.3333\). Here the robust standard error is smaller, illustrating that the sandwich estimate need not inflate uncertainty in every sample. The robust coefficient interval is \(0.50\pm1.96(0.2012)=(0.1056,0.8944)\), or a cause-specific hazard-ratio interval of \((1.111,2.446)\), since the estimated hazard ratio is \(\exp(0.50)=1.649\). This difference from model-based inference does not establish that clustering is harmless: the cluster-score pattern and number of independent clusters determine the estimate, and six clusters remain few for reliable asymptotic inference.
Assumptions, Diagnostics, and Alternatives
The cluster sandwich is most defensible when clusters are independent of one another, the number of independent clusters is reasonably large, and the cause-specific Cox estimating model is appropriate for the scientific target. Within a cluster, the method permits dependence in outcomes and score contributions. It does not require specifying the exact correlation pattern. The Cox model must still have an appropriate covariate form and proportional cause-specific hazard structure for the coefficient to have its intended interpretation. Independent censoring assumptions for the cause-specific analysis also remain relevant.
Before reporting the estimate, check that cluster membership is unambiguous and that the clusters represent the independent sampling units. Inspect the number of clusters, their sizes, the number of target failures across clusters, and the sizes of the cluster score contributions. A few very large or highly influential clusters can dominate the meat. Compare robust and model-based standard errors as a diagnostic, not as a rule for selecting whichever is preferred. Large changes warrant checking the data structure, event coding, cluster labels, and influential clusters. A cluster-level deletion analysis can reveal whether one cluster drives the result.
With few clusters, the uncorrected sandwich estimate and normal-based Wald tests may be biased or poorly calibrated. Small-sample corrections and resampling methods designed for clustered data are possible, but their suitability depends on the design and cluster count. Reporting many subjects is not a substitute for having many independent clusters. Also ensure that the software computes cluster-level score contributions for the same cause-specific model and tie method that produced \(\hat\beta\).
The model-based Cox variance uses the information matrix and relies on its working independent-subject variance assumptions. The cluster-robust version relaxes the within-cluster independence assumption while retaining the fitted coefficient and model target. A frailty model instead specifies a probability model for shared cluster effects; it can answer a different, often conditional, question and depends on that distributional specification. Neither approach automatically fixes an incorrect covariate effect, nonproportional hazards, or informative censoring. The robust approach is attractive when the main concern is variance misspecification within known independent clusters, rather than a desire to model the dependence itself.
Check Your Understanding
Answer the questions, then compare your reasoning with the solutions below.
- In a cause-specific analysis for cause \(k\), what happens to a subject's risk-set membership at the time of a competing event?
- For a scalar coefficient, suppose \(A=5\) and the cluster scores are \(1,-1,1,-1\). Calculate the robust variance and standard error.
- Does a cluster-robust variance estimate change the fitted coefficient? What does it change?
- Name two conditions or diagnostics that become especially important when the number of clusters is small.
- Why is a frailty model not simply another way to calculate the same robust standard error?
Solutions. (1) The subject leaves the cause-\(k\) risk set at the competing event time and is not counted as a cause-\(k\) failure. (2) The meat is \(1+1+1+1=4\), so the variance is \(4/5^2=0.1600\) and the standard error is \(\sqrt{0.1600}=0.4000\). (3) No; it changes the estimated covariance, standard errors, and associated intervals or tests. (4) Check the number of independent clusters and whether a few clusters dominate the score contributions; few clusters also make usual large-sample Wald inference less reliable. (5) A frailty model specifies a distribution for shared cluster effects and may target a conditional coefficient, whereas the sandwich estimates variance from cluster score variation without specifying that dependence distribution.