Why Clustered Variance Matters
A participant who experiences one recurrent event may be more likely than another participant to experience additional events. The event records from that participant are therefore not generally independent. A Cox partial likelihood can still provide an estimate of a covariate association under a specified recurrent-event intensity model, but the usual variance calculation treats the event contributions as if the working model fully described their dependence. When that assumption is doubtful, a subject-clustered sandwich variance allows arbitrary dependence among a subject’s event contributions while relying on independence across subjects.
The term marginal recurrent-event Cox model refers here to a model for the population-average event intensity, conditional on the recorded covariates and being at risk. Subject-clustered robust variance changes the estimated uncertainty, not the fitted partial-likelihood coefficient. It does not make a misspecified event-intensity model correct, nor does it remove dependence from the data.
Model and Data Structure
For subject \(i\), let \(N_i(t)\) count recurrent events through time \(t\), and let \(Y_i(t)\) indicate whether the subject is under observation and at risk just before \(t\). Let \(Z_i(t)\) be a vector of covariates, which may vary over time. A marginal proportional-intensity model is
Here \(\mathcal{F}_{t-}\) denotes the observed history just before time \(t\), \(\lambda_0(t)\) is an unspecified baseline intensity, and \(\beta\) is the vector of regression coefficients. Equivalently, the model specifies the event intensity while the subject is at risk. For a one-unit increase in a covariate, \(\exp(\beta_j)\) is the corresponding multiplicative intensity ratio, holding the other modeled covariates fixed. It is not a ratio of cumulative event probabilities.
In counting-process data, a subject’s follow-up is represented by intervals with start and stop times, an event indicator, covariates, and a subject identifier. Each interval contributes to a risk set while the subject is at risk; an event indicator of 1 marks an event at the interval endpoint, while 0 indicates no event there. Recurrent events require a new at-risk interval after an event when follow-up continues. The identifier must be the same for every record from the same person.
The model and variance calculation require independent subjects (or independent clusters), a correctly specified mean-intensity relationship for the coefficient interpretation sought, valid risk intervals and covariate histories, and censoring that is noninformative for the event process conditional on the modeled information. Robust inference also relies on having enough independent subjects for a large-sample approximation. Events within one subject may be dependent; that is the dependence the clustered variance is designed to accommodate.
Worked Example: Preparing Recurrent-Event Records
Suppose a hypothetical study follows four people for up to 12 months and records recurrent symptom episodes. Treatment is fixed at baseline. The records below show the event months; all subjects remain under observation through month 12.
| Subject | Treatment | Event months | Number of events |
|---|---|---|---|
| A | 0 | 3, 8 | 2 |
| B | 0 | 5 | 1 |
| C | 1 | 4, 7, 10 | 3 |
| D | 1 | 9 | 1 |
For an Andersen-Gill counting-process representation, subject A contributes an at-risk interval from month 0 to 3 with an event at 3, then an interval from 3 to 8 with an event at 8, and a final interval from 8 to 12 without an event. The same construction gives three event intervals for subject C and one final event-free interval after month 10. Subjects B and D each have one event interval followed by an event-free interval to month 12.
Treatment is repeated on every interval for its subject, and the subject identifier is repeated as well. The treatment coefficient compares event intensities under the model; it does not compare the proportion of subjects who ever had an event. With only four subjects, robust standard errors would be unreliable as large-sample inference, so these records illustrate data construction rather than support a dependable test.
Fitting and the Sandwich Variance
The Andersen-Gill (AG) approach fits a Cox model to these start-stop records using the risk sets defined by the chosen time scale. In software such as R’s survival package, a typical call is coxph(Surv(start, stop, event) ~ treatment + age + cluster(id), data = ...). The cluster(id) term requests a robust variance that aggregates contributions by subject. Without it, the same partial likelihood can be fitted with the usual model-based variance. Exact function options vary by software, so confirm that the reported variance is robust and that the cluster identifier is the subject identifier.
Let \(\hat{\beta}\) be the fitted coefficient vector, \(A\) the observed information matrix (the negative derivative of the partial-likelihood score), and \(U_i(\hat{\beta})\) the sum of the score contributions from all event records for subject \(i\). The uncorrected subject-clustered sandwich estimate is
where \(m\) is the number of independent subjects. The usual model-based variance is \(A^{-1}\). Some software applies a finite-sample correction; if so, the reported robust variance may differ slightly from this uncorrected expression.
The justification is a cluster-level version of the sandwich argument. At the fitted estimate, the total score is approximately zero, and a first-order expansion of the score around the target coefficient gives
If subjects are independent, the variance of the sum of their score contributions is the sum of their score variances. Estimating that sum by the outer products \(U_iU_i^\mathsf{T}\), then transforming by \(A^{-1}\) on both sides, yields the sandwich formula. Contributions from events within one subject are combined before taking the outer product, so their within-subject covariance is retained instead of being treated as independent event-level information.
Worked Example: Computing a Scalar Cluster-Robust Standard Error
For a hypothetical one-covariate fit, suppose the observed information at the fitted coefficient is \(A=20\), and the four subject-level score contributions are \(1.2\), \(-0.8\), \(0.5\), and \(-0.9\). Their sum is zero, as expected at the fitted value in this simplified illustration. Use the uncorrected sandwich formula.
The scalar “meat” is the sum of squared subject scores:
Since \(A^{-1}=1/20=0.05\), the robust variance and standard error are
The model-based variance would instead be \(1/20=0.05\), with standard error \(\sqrt{0.05}=0.224\). These illustrative values show that a sandwich estimate is not guaranteed to be larger: it uses the observed variability in cluster scores and can be smaller or larger. Four clusters are far too few to regard either large-sample standard error here as reliable.
Comparing a Robust Fit with an AG Fit
A common comparison is between (1) the AG partial-likelihood fit with its usual model-based variance and (2) that same AG fit with subject-clustered robust variance. If the records, risk sets, covariates, tie handling, and partial likelihood are identical, the coefficient estimates are identical: only the variance estimate changes. The AG model describes a common covariate effect on the event intensity across recurrences. Adding clustering does not create a different mean model.
This distinction matters for interpretation. The robust standard error supports inference that allows dependence among a subject’s events, under the independent-subject and model conditions above. It does not alter the fitted intensity ratio, change the risk sets, allow an event-order-specific coefficient, or model a subject-specific latent propensity. Other recurrent-event specifications can define different risk sets or effects, and their estimates need not agree with the AG estimate.
Worked Example: Comparing Coefficients and Inference
Suppose a hypothetical analysis uses the same AG records and partial likelihood for both fits. The estimated treatment coefficient is \(-0.405\) in each fit. The conventional model-based standard error is \(0.180\), while the subject-clustered robust standard error is \(0.260\). Compute the hazard-intensity ratio and compare the two Wald intervals.
The estimated intensity ratio is \(\exp(-0.405)=0.667\). Thus, under the fitted model, the treatment group’s recurrent-event intensity is estimated to be 0.667 times the comparison group’s intensity, conditional on the modeled covariates and being at risk. This is an estimated 33.3% lower intensity, not a claim of a 33.3% lower probability of any event.
Using the robust standard error, the Wald interval on the coefficient scale is
Exponentiating the endpoints gives an approximate 95% interval for the intensity ratio of \((\exp(-0.915),\exp(0.105))=(0.401,\ 1.111)\). Because this interval includes 1, the robust Wald test does not reject a null treatment coefficient of zero at the 0.05 level. Its test statistic is \(z=-0.405/0.260=-1.56\), with a two-sided normal-approximation p-value of about 0.119.
With the conventional standard error, the interval is \(-0.405\pm1.96(0.180)=(-0.758,-0.0522)\), or an intensity-ratio interval of approximately \((0.469,0.949)\). The corresponding test statistic is \(-2.25\), with a two-sided p-value of about 0.024. The difference in inference comes from the variance estimate, not from a change in the coefficient or its estimated intensity ratio. The robust interval is an approximate confidence interval; its repeated-sampling coverage is not a probability statement about this particular interval.
Diagnostics, Pitfalls, and Limitations
Before interpreting a fitted model, check that the data encode the event history correctly and that the variance targets the intended independent units. Useful checks include confirming that intervals for each subject are ordered and do not overlap improperly, event indicators mark the intended event times, and every record for a person carries the same cluster identifier. Check that subjects—not event records—are the independent units. Review the number of subjects, the spread of events across subjects, and whether a small number of subjects contribute a large share of the total score.
A sandwich variance is an asymptotic procedure. With few subjects, or highly unbalanced cluster sizes, its standard errors and Wald tests may be poorly calibrated even if the code runs without warnings. Robust variance also does not fix informative dropout, omitted covariates that distort the mean intensity, incorrect time-dependent covariates, or an inappropriate risk-set construction. These require substantive modeling decisions or additional methods, not merely a robust option. In particular, a robust variance should not be presented as evidence that the intensity model itself fits well.
A frequent mistake is to say, “The robust option accounts for dependence, so the fitted model no longer needs assumptions about the event process.” The correction is that robust variance permits within-subject dependence in the score contributions for inference, while the coefficient still comes from the specified partial likelihood and its intensity-model interpretation. If that mean model is wrong, the variance option does not repair the resulting estimand or bias.
The main alternative within this comparison is the conventional AG fit with model-based variance, which is more efficient when its working model fully captures the relevant dependence but can understate uncertainty otherwise. The subject-clustered sandwich trades some model-based efficiency for variance robustness across subjects, subject to adequate cluster count and a correctly specified mean model. A shared frailty model is a different approach that explicitly represents subject-level latent heterogeneity; it is not equivalent to clustering the variance and answers a different modeling question.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- In a recurrent-event data set with several rows per person, which identifier should define the robust cluster, and why?
- If the same AG partial likelihood is fitted with and without a cluster-robust variance, what should happen to the coefficient estimate?
- A treatment coefficient is \(0.20\). Interpret \(\exp(0.20)\) in the recurrent-event intensity model.
- Name two conditions that a subject-clustered sandwich variance does not remove.
Solutions. (1) Use the subject identifier, because all event records from one person may be dependent and must be combined into one cluster contribution. (2) The coefficient estimate should remain the same when risk sets, covariates, ties, and partial likelihood are identical; the estimated variance may change. (3) \(\exp(0.20)=1.221\), so the modeled event intensity is estimated to be about 1.22 times as high for a one-unit increase in that covariate, holding other modeled covariates fixed and conditional on being at risk. (4) Examples include a misspecified mean-intensity model, informative censoring or dropout, incorrect risk intervals, and having too few independent subjects for reliable large-sample inference.