Why Use an Andersen-Gill Model?
A person may experience the same type of event more than once: repeated hospital admissions, asthma exacerbations, or infections. A standard time-to-first-event analysis discards events after the first, even when those later events are part of the scientific question. The Andersen-Gill (AG) model uses the event history while the person remains at risk, representing each person’s follow-up with start-stop intervals.
The AG model is a counting-process form of the Cox model. Its regression coefficient describes a multiplicative association with the instantaneous recurrent-event rate among people currently at risk, conditional on the covariates in the model. Events from the same person may be dependent. A subject-level robust variance estimator adjusts uncertainty estimates for that within-person dependence; it does not change the fitted coefficient.
Here \(\mathcal{H}_i(t)\) denotes the subject’s observed history up to \(t\), \(\lambda_0(t)\) is an unspecified baseline intensity, and \(X_i(t)\) is a covariate vector known just before \(t\). The model uses the same baseline intensity and regression effects for each event order. A subject who has already had an event can remain at risk and contribute later events.
Representing Follow-Up in Start-Stop Form
Each row of the analysis data describes one interval \((\text{start},\text{stop}]\) during which a subject is at risk. The event indicator is 1 if the event occurs at the interval’s stop time and 0 otherwise. After a recurrent event, start a new interval at that event time if the subject remains under observation. The intervals for a subject must not overlap, and a subject’s covariates on a row must represent information available during that interval.
For example, suppose a hypothetical participant is followed from time 0 through time 8, has events at times 2 and 5, and is then censored at time 8. Their three rows are \((0,2]\) with event 1, \((2,5]\) with event 1, and \((5,8]\) with event 0. Under the usual start-stop convention, the risk set at event time \(t\) includes a row when \(\text{start}<t\leq\text{stop}\).
Worked Example: Turning Event Histories into Intervals
Consider three hypothetical participants. Participant A receives treatment and has events at times 2 and 5, then is censored at 8. Participant B is untreated and has events at 3 and 7, then is censored at 8. Participant C receives treatment and has an event at 4, then is censored at 8.
| id | treatment | start | stop | event |
|---|---|---|---|---|
| A | 1 | 0 | 2 | 1 |
| A | 1 | 2 | 5 | 1 |
| A | 1 | 5 | 8 | 0 |
| B | 0 | 0 | 3 | 1 |
| B | 0 | 3 | 7 | 1 |
| B | 0 | 7 | 8 | 0 |
| C | 1 | 0 | 4 | 1 |
| C | 1 | 4 | 8 | 0 |
There are \(3+3+2=8\) interval rows, of which \(2+2+1=5\) end in events. At time 5, A contributes an event from the row \((2,5]\); B contributes from \((3,7]\), and C from \((4,8]\). Thus all three participants are in the risk set just before time 5. A’s earlier event does not remove A from this risk set.
In real data, if treatment or another covariate changes during follow-up, split the record at the change time and assign each interval the covariate value then in effect. Do not assign a later value to an earlier interval: doing so uses future information.
Fitting the Model with coxph
In R, the start-stop response is specified as Surv(tstart, tstop, event). Include a subject identifier as a cluster term in the model formula to request a robust variance estimate clustered by subject. For example, for a binary treatment indicator coded 0 or 1, the call is:
R call: coxph(Surv(tstart, tstop, event) ~ treatment + cluster(id), data = agdat)
The interval rows identify who is at risk at each event time; cluster(id) identifies which rows belong to the same independent subject for variance estimation. The model’s partial likelihood uses the Cox risk sets at observed event times. It does not require specifying the shape of \(\lambda_0(t)\).
Worked Example: Interpreting an AG Treatment Estimate
Suppose a hypothetical analysis has 120 independent participants and 184 recurrent events. Treatment is coded 1 for treated and 0 for untreated, and a fitted model returns \(\hat{\beta}=0.405\) with subject-clustered robust standard error \(0.220\). These are stipulated illustrative software results, not observations from a real study.
For a one-unit increase in the treatment indicator, the estimated intensity ratio is
Under the model, the estimated event intensity while at risk is 1.50 times as high for the treated group as for the untreated group, at the same modeled covariate values. This is an intensity comparison, not a ratio of the probabilities of ever having an event.
A 95% Wald confidence interval on the log-intensity scale is \(\hat{\beta}\pm1.96SE\). Substituting the robust standard error gives \(0.405\pm1.96(0.220)=0.405\pm0.4312\), or \((-0.0262,0.8362)\). Exponentiating the endpoints gives
The interval includes 1, so these illustrative data do not provide statistically significant evidence of a treatment-associated intensity difference at the two-sided 5% level. The interval is compatible with an intensity slightly lower or substantially higher in the treated group; it is not a 95% probability statement about the fixed parameter.
The corresponding Wald statistic is \(z=0.405/0.220=1.84\), giving a two-sided \(p\)-value of approximately 0.066. The \(p\)-value is the probability, under the null hypothesis and the test’s assumptions, of a statistic at least this extreme; it is not the probability that the null hypothesis is true.
Why Clustered Robust Variance Is Used
The usual Cox variance calculation relies on a model-based information matrix. For recurrent events, the event contributions from one participant need not be independent. The robust, or sandwich, variance uses subject-level score contributions so that arbitrary within-subject dependence can be accommodated, while independent subjects remain the sampling units.
For clarity, let \(U_i(\beta)\) be subject \(i\)’s total partial-likelihood score contribution, and let \(A\) be the negative derivative of the total score with respect to \(\beta\), evaluated at the fitted coefficient. Ignoring finite-sample corrections, the robust variance estimator is
The middle term accumulates variability across subjects, not across interval rows. The reason is visible from a first-order expansion of the score equation. Since the fitted total score is zero, \(0\approx U(\beta_0)-A(\hat{\beta}-\beta_0)\), so \(\hat{\beta}-\beta_0\approx A^{-1}\sum_i U_i(\beta_0)\). Taking the variance across independent subjects gives the sandwich form: the variance of the summed subject scores is estimated by \(\sum_i U_iU_i^\top\), surrounded by the inverse information terms.
This adjustment protects standard errors and related inference against certain forms of within-person dependence. It does not repair a wrongly specified intensity model, make dependent subjects independent, or correct bias from informative censoring. If the scientific aim is to model event-order-specific risk sets or a subject-level random effect, those require different model structures.
Worked Example: Why the Event Rows Are Not Independent Clusters
Suppose a hypothetical data set contains 60 participants, each with one to four interval rows, for a total of 150 rows. The independent sampling units are the 60 participants, not 150 separate intervals. The clustered variance therefore groups rows by participant ID and sums each participant’s score contributions before forming the outer products.
If the treatment coefficient is \(0.30\) and the robust standard error is \(0.15\), the robust Wald statistic is \(0.30/0.15=2.00\). The corresponding two-sided normal-approximation \(p\)-value is about 0.0455. Treating all 150 rows as independent could produce a different standard error because it would fail to account for repeated contributions by the same 60 people. The coefficient and standard error values here are illustrative; the calculation shows why the cluster unit must match the independent subject.
Assumptions, Diagnostics, and Alternatives
The AG model and its robust inference require careful attention to the observation process and the data structure. In particular:
- Correct risk intervals: start and stop times must reflect actual time at risk. Intervals for the same participant must not overlap, and event indicators must mark the appropriate event-ending rows.
- Predictable covariates: each covariate value must be available just before the time it is used. Time-varying values need correctly split intervals.
- Independent clusters: the sandwich calculation allows dependence within a participant, but assumes independent participants (or independent clusters).
- Conditionally noninformative censoring: given the modeled covariates and history, censoring must not selectively remove participants in a way that invalidates estimation of the target intensity.
- Proportional intensity effects: unless modeled otherwise, a covariate’s coefficient is constant over the analysis time scale. The AG structure also assumes a common baseline intensity across event orders.
Check the interval construction before interpreting model output: look for negative or zero-length intervals, overlapping rows, events after censoring, and event counts that do not match the source histories. Confirm that a participant is still represented after a nonterminal event and that no intervals extend beyond their censoring time. Then assess proportional-hazards behavior, for example with scaled Schoenfeld residuals or covariate-by-time checks. A marked time pattern suggests the constant-coefficient form may be inadequate.
A robust standard error that differs from the model-based standard error is a useful reminder that within-subject dependence affects uncertainty. It is not, by itself, a diagnostic proving that either estimate is right. Very few independent participants can also make sandwich inference unstable, even when the number of interval rows or events is large. Tied event times should be handled with a documented tie method; the default Efron method in coxph is common, but the event-time resolution and tie frequency should be reviewed.
The AG model is convenient when the question concerns recurrent-event intensity over a shared time scale and event order does not define separate risk sets. A model stratified by event number, such as a Prentice-Williams-Peterson approach, instead makes event order central to the risk-set construction. A frailty model explicitly represents latent subject-level heterogeneity and requires additional distributional assumptions. These alternatives answer different questions; robust variance alone does not turn the AG model into a frailty model.
A final limitation is interpretive: an intensity ratio is not a cumulative probability, expected event count, or causal effect. A causal interpretation needs appropriate design and confounding assumptions in addition to a correctly specified model. When a terminal event ends the possibility of further recurrence, it also changes who can experience future events; treating that terminal event as ordinary independent censoring is not automatically justified.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A participant has events at times 4 and 9 and is censored at 12. Write the correct start-stop event rows.
- In an AG model, what does \(\exp(\beta)\) represent for a binary covariate?
- Why is the robust variance clustered by subject rather than by interval row?
- A fitted coefficient is 0.20 with robust standard error 0.10. Calculate the hazard-intensity ratio and its 95% Wald interval.
- Name one problem that a robust subject-level variance estimator does not fix.
Solutions. (1) The rows are \((0,4]\) with event 1, \((4,9]\) with event 1, and \((9,12]\) with event 0. (2) It is the ratio of conditional event intensities for covariate value 1 versus 0, with other modeled covariates held fixed. (3) Rows from the same participant may be dependent, while participants are treated as independent clusters; clustering sums score contributions within participant. (4) The ratio is \(\exp(0.20)=1.22\). The log-scale interval is \(0.20\pm1.96(0.10)=(0.004,0.396)\), so the intensity-ratio interval is approximately \((1.00,1.49)\). (5) Examples include informative censoring, a misspecified intensity or proportional-effects model, and dependence between participants.