A Decision Framework for Recurrent-Event Studies
Recurrent-event analyses can produce different answers from the same data because they may target different quantities, measure time differently, or compare different risk sets. Choosing a model is therefore not a matter of selecting a familiar software option first. Begin with the scientific question: is the target the number of events accumulated by a horizon, an average event rate, or the instantaneous event hazard among people currently eligible to experience a particular recurrence? Then choose the clock and risk set that express that target.
The earlier material on mean cumulative functions and fixed-horizon burden contrasts distinguishes event burden from rate and event-hazard comparisons. This workshop brings those choices together with dependence-robust inference and reporting. It does not make the estimands interchangeable: a model’s coefficient answers only the question defined by its intensity, time scale, and eligible risk sets.
Choose the Estimand Before the Model
If the clinical question concerns total burden, a natural target is the mean cumulative function \(M(t)=E\{N_i(t)\}\): the expected number of events per subject accumulated through time \(t\) in a specified population. A contrast \(M_1(t)-M_0(t)\) is an absolute difference in expected events per subject by \(t\); a ratio \(M_1(t)/M_0(t)\), when its denominator is nonzero, is a relative burden contrast. These summaries do not by themselves describe an instantaneous hazard.
If the question concerns average event occurrence per unit of observed time, a marginal rate model may be suitable. If it concerns the event intensity among subjects who have reached a particular recurrence order and remain eligible, an event-specific hazard model may better match the question. A hazard ratio compares instantaneous rates among the risk sets defined by the model, conditional on the included covariates and model structure. It is not generally a ratio of cumulative event counts or of population-average rates.
Worked Example: Matching a Burden Question to an Estimand
Suppose a hypothetical randomized study follows 100 participants in each group for 12 months, with complete follow-up. At month 12, the estimated mean cumulative numbers of hospitalizations are 1.20 per participant in the intervention group and 1.60 in the control group. The scientific question is: “How many fewer hospitalizations per participant are expected by 12 months?”
The target is an absolute fixed-horizon burden contrast: \(M_1(12)-M_0(12)=1.20-1.60=-0.40\) hospitalizations per participant. Thus the intervention group has an estimated 0.40 fewer hospitalizations per participant by 12 months. The relative burden is \(1.20/1.60=0.750\), or 25.0% lower, but that is a different summary.
Neither result is automatically a hazard ratio. The burden contrast aggregates events over time, whereas a hazard ratio compares instantaneous event intensities in specified risk sets. If participants had unequal follow-up, the simple group means would also require appropriate handling of censoring and the target population; complete follow-up here makes that issue absent.
Choose the Clock and Risk Set Together
A total-time clock measures time from a common origin, such as randomization or study entry, regardless of earlier recurrences. It is appropriate when event intensity may depend on time since that origin, or when the aim is to compare recurrent-event experience over follow-up time. In an Andersen–Gill-style counting-process model, a subject may contribute repeatedly while under observation; the risk indicator \(Y_i(t)\) determines when that subject contributes. This formulation does not reset the clock after each event.
A gap-time clock resets time to zero after each event and measures time since the most recent recurrence. It is useful when the interval since the previous event is scientifically important. The risk set must still be defined: an event-specific conditional model, such as a Prentice–Williams–Peterson formulation, includes a subject in the risk set for event order \(j\) only after that subject has experienced event \(j-1\) and remains under observation. A total-time version instead indexes each event-order risk set by time since the common origin.
These choices change the comparison being made. A gap-time model for the third event concerns subjects who reached the third-event risk set, not all enrolled subjects. An event-order-specific model can also allow different baseline hazards, and possibly different covariate effects, by recurrence order. By contrast, a shared-effect model imposes a common covariate coefficient across the event contributions it includes. Such choices should follow the estimand, not be made solely because one model converges more easily.
| Analysis choice | Time and eligibility | Typical interpretation |
|---|---|---|
| Mean cumulative function | Follow-up from a common origin; population-level event counts | Expected cumulative events per subject by time \(t\) |
| Marginal rate model | Events per unit of at-risk time; specify how exposure time is counted | Rate ratio for the modeled population-average rate, subject to model assumptions |
| Andersen–Gill-style model | Usually total time; subjects contribute when under observation and at risk | Intensity ratio for the specified recurrent-event process and covariates |
| Event-specific conditional model | Risk set restricted to subjects who reached the required event order; total-time or gap-time clock | Hazard ratio within the event-order-specific risk set |
Worked Example: Selecting a Clock and Risk Set
In a hypothetical study of recurrent device failures, investigators record time from implantation to each failure. They want to know whether the interval since the last repair predicts the next failure, and they want the effect for the second failure among devices that experienced a first failure.
A gap-time clock matches the first question because it measures time from the most recent repair or failure, rather than time since implantation. For the second-failure analysis, the risk set includes only devices with a recorded first failure that remain under observation and eligible for a second failure. Devices with no first failure cannot contribute to that event-order risk set.
For example, suppose 30 of 80 devices experience a first failure; 4 of those are censored before a possible second failure, leaving 26 eligible devices at the start of second-event follow-up. The second-failure analysis begins with those 26, not all 80. This choice targets a conditional event-order comparison. It does not estimate the effect on overall failure burden in all implanted devices.
Dependence-Robust Inference and Its Scope
Repeated event contributions from one subject are generally dependent. A subject-level sandwich variance estimator accounts for this dependence in the estimated uncertainty when the model’s estimating equations have the appropriate mean structure and subjects are independent of one another. Let \(U_i(\hat\beta)\) be the contribution of subject \(i\) to the fitted estimating equation, and let \(A\) be the negative derivative of the total estimating equation with respect to the coefficient vector. A common sandwich estimate is:
The middle term aggregates subject-level score variation, so all event contributions from a subject stay together when estimating variability. The “bread” terms map that variation to coefficient scale. This is a variance estimator, not a change to the estimand or fitted coefficient. It does not repair an incorrect risk set, a misspecified mean or intensity model, informative censoring, confounding, or dependence between different subjects. With few independent subjects, the sandwich estimate can also be unreliable; a large number of events from a small number of people does not substitute for many independent clusters.
Worked Example: Reporting a Robust Rate-Ratio Estimate
In a hypothetical marginal rate analysis of 240 independent participants, the fitted intervention coefficient is \(\hat\beta=\log(1.50)=0.4055\). The subject-level robust standard error is \(SE(\hat\beta)=0.180\). Assume the model’s rate structure and censoring assumptions are judged reasonable, and the target is the intervention-to-control rate ratio.
The estimated rate ratio is \(\exp(0.4055)=1.50\). A large-sample 95% confidence interval on the coefficient scale is:
Exponentiating the endpoints gives the interval for the rate ratio: \((\exp(0.0527),\exp(0.7583))=(1.054,\ 2.135)\). The Wald statistic is \(z=0.4055/0.180=2.253\), giving a two-sided \(p\)-value of approximately 0.024. Under the model, the estimated intervention-group rate is 1.50 times the control-group rate; the 95% confidence interval ranges from 1.054 to 2.135. The result is inconsistent with a rate ratio of 1 at the 0.05 level. This is a rate comparison, not a claim that the intervention causes a 50% increase unless the design and causal assumptions support that interpretation.
Reporting the Analysis and Checking Its Fit
A useful report makes the analysis reproducible and the effect interpretable. State the population, event definition, follow-up origin, clock, risk-set rule, estimand, model structure, covariates, censoring approach, and variance estimator. Report the effect estimate and uncertainty in the units or scale of the estimand. If using a log-link rate model, exponentiate the coefficient and call it a rate ratio; if using a hazard model, identify the event order and risk set to which the hazard ratio applies.
Diagnostics should be tied to the model choices. Check event counts and follow-up by group and recurrence order; verify that event times and time resets were coded as intended; inspect whether event-order-specific effects or time-varying effects are needed; and assess whether influential subjects dominate the robust variance. Compare model-based and robust standard errors when available, but do not treat agreement as proof that assumptions hold. If the number of independent subjects is small, report this limitation and consider an inference method appropriate to the sampling structure.
The main alternatives answer different questions. Mean cumulative functions and fixed-horizon contrasts emphasize accumulated burden and are directly interpretable in events per subject. Marginal rate models summarize events per at-risk time. Recurrent-event hazard models describe intensity within their defined risk sets and may incorporate event order or time since the prior event. None is universally preferable. A sensitivity analysis using a scientifically defensible alternative clock or risk-set scheme can reveal whether conclusions depend on modeling choices, but it does not turn one estimand into another.
Common Mistakes and Limitations
- Calling every coefficient a hazard ratio. In a marginal rate model, exponentiated coefficients are rate ratios; the interpretation depends on the modeled rate.
- Ignoring event-order eligibility. A conditional second-event estimate applies to subjects who reached the second-event risk set, not automatically to the full enrolled population.
- Treating robust variance as a model cure. It addresses within-subject dependence in uncertainty under its assumptions; it does not fix a wrong clock, wrong risk set, biased coefficient, or informative censoring.
- Reporting only a p-value. Include the effect scale, estimate, confidence interval, units or context, and the risk set or population to which the result applies.
Recurrent-event models can also be difficult to interpret when follow-up differs substantially, covariates change over time, later event orders become sparse, or event occurrence affects subsequent eligibility. Robust inference does not remove those limitations. The report should state which individuals and event histories contribute to the estimate and avoid extending a conditional association to a broader population without justification.
Specify burden, rate, or event-specific hazard, including population and horizon where relevant.
Choose total time or gap time, then state precisely who is eligible to contribute and when.
Assess covariate structure, event-order effects, censoring, and the number of independent subjects.
Give the estimand, effect estimate, robust uncertainty, assumptions, and limitations.
Check Your Understanding
Answer each question, then compare with the solutions below.
- A study asks for the expected number of recurrent infections per participant by 18 months. Which estimand is the most direct match?
- For a conditional third-event model, which subjects belong in the risk set?
- A fitted log rate coefficient is 0.20. What is the corresponding rate ratio, approximately?
- Name one issue that a subject-level robust variance estimator does not fix.
Solutions. (1) The mean cumulative function \(M(18)\), or a group contrast in \(M(18)\) if comparing groups, directly targets cumulative event burden per participant by 18 months. (2) Subjects who have experienced the second event and remain under observation and eligible for a third event; the precise time scale must also be stated. (3) \(\exp(0.20)=1.221\), so the modeled rate is estimated to be about 1.22 times the reference rate, holding the model’s other terms fixed. (4) It does not fix, for example, an incorrect risk set, informative censoring, a misspecified mean structure, or confounding; any one of these answers is sufficient.