Why Recurrent Events Need a Clear Risk-Set Definition
A person may experience multiple infections, hospital admissions, or disease flares during follow-up. Treating only the first event as the outcome discards subsequent events, but simply analyzing all events as independent observations ignores that events from one person share biology and history. With omics predictors, the modeling choice also determines whether an association describes an overall repeat-event rate, a rate conditional on reaching a particular event order, or a rate conditional on a person's latent susceptibility.
The three model families in this tutorial address these questions differently. Andersen-Gill (AG) uses a common event-time scale and a common baseline rate; Prentice-Williams-Peterson (PWP) models distinguish event order; frailty models represent unobserved differences in susceptibility between people. These are not interchangeable adjustments for “repeated data.” The risk set and interpretation must match the question.
For all three approaches, predictors used at time \(t\) must be available before the event at \(t\). A time-updated assay should therefore be aligned to a prespecified measurement window or lag, rather than using a measurement collected after an event to predict that event. The models also require an adequately specified relationship between predictors and event rate, and censoring must be non-informative conditional on modeled information. If follow-up ends because of a terminal event that is related to recurrent-event risk, treating that event as ordinary independent censoring may be inadequate; the recurrent-event models below do not by themselves resolve that problem.
Andersen-Gill: A Common Event-Time Scale
The AG model is an extension of the Cox proportional-hazards model to recurrent events. Each person can contribute several event times, with time measured from a common origin such as study entry. All events share a baseline hazard, and the same covariate effect applies at every event order:
Here \(\mathcal{F}_{i}(t-)\) denotes the person's observed history immediately before \(t\), \(\lambda_0(t)\) is the unspecified baseline rate, and \(\boldsymbol{\beta}\) contains log rate-ratio coefficients. A one-unit increase in predictor \(x_j\), holding other modeled covariates fixed, multiplies the conditional event rate by \(\exp(\beta_j)\), provided the proportional-rate model is appropriate. If \(x_j\) is standardized, the unit is one standard deviation.
The usual AG partial likelihood uses the risk set of people currently under observation, not a new risk set for each event order. Repeated event times within a person are not independent observations. A subject-clustered sandwich variance is commonly used to obtain standard errors robust to within-person dependence, while retaining independent sampling across people. This variance adjustment does not repair a misspecified mean intensity, an invalid risk set, or informative censoring.
Worked Example: An Omics Rate Ratio in an Andersen-Gill Model
Suppose a hypothetical study follows 80 patients for recurrent inflammatory flares. A prespecified, standardized omics score is measured before each relevant risk interval. The fitted AG model estimates a score coefficient of \(0.400\), with subject-clustered standard error \(0.180\). Assume the common baseline rate and proportional score effect are scientifically reasonable, predictor values precede events, and censoring is conditionally non-informative.
The estimated rate ratio for a one-standard-deviation increase in the score is
Thus, at the same follow-up time and covariate values, the fitted model estimates a 49.2% higher recurrent-flare rate per standard-deviation increase in score. This is an association in the modeled event rate, not evidence that the score causes flares.
Using a large-sample normal approximation, the test statistic for \(H_0:\beta=0\) is
The two-sided p-value is approximately \(0.0263\). A 95% confidence interval for the coefficient is
Exponentiating the endpoints gives a 95% confidence interval for the rate ratio of \((\exp(0.0472),\exp(0.7528))=(1.048,\ 2.123)\). Under the model and large-sample approximation, this procedure has 95% coverage in repeated sampling. The interval is compatible with rate increases from about 4.8% to 112.3% per standard deviation.
PWP: Conditioning on Event Order
The PWP approach recognizes that a person can be at risk for a second event only after experiencing a first event, and at risk for a third only after experiencing two. It fits a separate baseline hazard for each event order. In a total-time PWP model, time remains measured from study entry; in a gap-time PWP model, the clock resets after the preceding event. Either way, the risk set for event \(k\) contains people who have experienced event \(k-1\) and have not yet had event \(k\), while still under observation.
The event-order-specific baseline \(\lambda_{0k}(t)\) allows the baseline rate to differ for the first, second, and later events. If the same \(\boldsymbol{\beta}\) is used across orders, the model still assumes a common covariate rate ratio across those orders. The coefficient is conditional on being eligible for event \(k\), given the required preceding events. This conditional interpretation is often useful when recurrence itself changes who remains at risk.
Worked Example: Reading an Event-Order PWP Effect
In a hypothetical study of repeated disease flares, researchers fit a total-time PWP model with event-order-specific baseline hazards. For the same standardized omics score, suppose the common fitted coefficient is \(-0.300\), with standard error \(0.140\). Assume the risk sets correctly include only patients who have reached the relevant event order, the score is measured before each event interval, and the covariate effect is reasonably constant across event orders.
The fitted conditional rate ratio is
Among patients eligible for a given event order, a one-standard-deviation increase in the score is associated with an estimated 25.9% lower rate of that event, conditional on the model covariates. The estimate is not a comparison of all patients regardless of whether they experienced the preceding events.
The approximate test statistic is \(z=-0.300/0.140=-2.143\), giving a two-sided p-value of about \(0.0321\). The coefficient interval is \(-0.300\pm1.96(0.140)=(-0.5744,-0.0256)\). Exponentiating yields a rate-ratio interval of approximately \((0.563,\ 0.975)\). The interval excludes 1 under this approximation, but it does not remove the need to assess whether a shared score effect across event orders is credible.
Frailty Models: Representing Persistent Susceptibility
A frailty model adds a latent person-specific multiplicative effect. A common shared gamma frailty specification is
Here \(U_i\) is a positive random effect shared by all recurrent events for person \(i\); the gamma distribution is often chosen because it is positive and mathematically convenient. Conditional on \(U_i\) and the modeled covariates, the recurrent-event process is assumed to follow the specified intensity model. A positive \(\theta\) represents unobserved heterogeneity: some people have persistently higher event rates than others. The coefficient \(\beta_j\) gives a rate ratio conditional on the same frailty value, not automatically the population-average rate ratio after averaging over frailty.
A frailty term can be combined with an AG-style or event-order-specific risk structure. It answers a different modeling question from PWP stratification: PWP allows baseline rates to vary by event order, whereas frailty represents persistent person-level susceptibility. Both features may be needed, but extra complexity requires enough events and a scientifically justified model.
Worked Example: How Shared Frailty Connects Repeat Events
Consider a hypothetical gamma-frailty model with \(\theta=0.500\). For a patient with a fixed covariate profile, suppose the baseline cumulative rate over interval A is \(0.400\), and over a disjoint interval B it is \(0.600\). Conditional on the frailty, assume event counts in the two intervals are independent Poisson counts with means \(U_i(0.400)\) and \(U_i(0.600)\). These assumptions make the dependence calculation explicit.
The law of total covariance gives
The first term is zero by conditional independence. The second is
The positive covariance arises because the same high or low frailty affects both intervals. If the fitted score coefficient were \(0.250\), its conditional rate ratio would be \(\exp(0.250)=1.284\): among people with the same frailty and other covariates, a one-unit score increase corresponds to an estimated 28.4% higher rate. The value of \(\theta\) describes residual heterogeneity; it is not itself an omics effect.
Choosing and Checking a Model
Start by defining the scientific estimand and the risk set, rather than choosing a model from software defaults. AG is a natural option when the target is an overall event-rate association on a common time scale and event order is not central. PWP is preferable when eligibility for the next event depends essentially on prior events and the conditional rate at each order is the target. A frailty model is useful when persistent, unobserved patient-level susceptibility is substantively plausible and dependence should be represented parametrically. Robust AG variance and frailty are not equivalent: the former adjusts uncertainty, while the latter specifies a latent heterogeneity distribution and changes the likelihood and interpretation.
| Approach | Risk-set and baseline structure | Omics coefficient interpretation | Main caution |
|---|---|---|---|
| Andersen-Gill | All currently observed, at-risk people; common baseline on the chosen time scale | Rate ratio for recurrent events under a common proportional effect | Event-order differences are not modeled by separate baselines; use subject-clustered uncertainty when appropriate |
| PWP | Separate risk set and baseline for each event order | Rate ratio conditional on eligibility for that event order | Later-order estimates may be imprecise because fewer people reach those risk sets |
| Frailty | Risk structure depends on the chosen formulation; person-specific latent multiplier | Conditional rate ratio for people with the same frailty | Results depend on the frailty distribution and conditional model assumptions |
Define exactly what counts as a recurrence, when follow-up begins and ends, and whether elapsed time means time since entry or time since the preceding event.
Check which people can contribute to each event time or order. For PWP, verify that a person enters the kth-event risk set only after the first \(k-1\) events.
Use only predictor values available before the interval's event. Consider whether a common effect across time or event orders is plausible.
Inspect event counts by person and order, covariate overlap in risk sets, influential observations, and proportional-rate patterns. Compare model-based and subject-clustered standard errors where relevant, and examine whether frailty variance is supported rather than effectively zero.
Warning signs include a sharp decline in the number of people at risk at later event orders, unstable PWP estimates, or an AG coefficient that changes substantially when event order is represented. Check whether the omics effect appears time-varying; a single proportional rate ratio can conceal a changing association. For frailty models, a variance estimate near zero suggests little estimated residual heterogeneity under that model, while a large estimate may indicate strong heterogeneity or model misspecification. Neither result alone proves that a particular structure is correct.
Common alternatives answer different questions. A first-event Cox model is simpler but discards later recurrences. A robust AG model targets an event-rate association with a less restrictive variance calculation, but it does not estimate latent susceptibility. PWP explicitly conditions on event order but can leave few observations for late recurrences. Frailty accounts for unobserved heterogeneity but relies on a distributional assumption and yields a conditional effect. In every approach, a well-estimated association remains observational unless the study design and causal assumptions support a causal interpretation.
Common Mistake
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A fitted AG coefficient for a one-standard-deviation omics score is \(0.200\). What is the estimated rate ratio, and how should it be interpreted?
- In a PWP model for third events, who belongs in the risk set?
- For a gamma frailty model with \(\theta=0.400\), interval cumulative rates \(0.500\) and \(0.300\), and conditional independence given frailty, calculate the covariance of the two interval event counts.
- Which approach directly represents unobserved, persistent differences in patient susceptibility: robust AG variance, PWP event-order strata, or a frailty term?
Solutions. (1) \(\exp(0.200)=1.221\), so the estimated event rate is 22.1% higher per standard deviation, conditional on the AG model covariates and assumptions. (2) Patients who have experienced two events, remain under observation, and have not yet experienced their third event. (3) By the same total-covariance argument, \((0.500)(0.300)(0.400)=0.0600\). (4) A frailty term represents persistent latent susceptibility; robust AG variance changes the estimated uncertainty, and PWP stratifies by event order.