Recurrent Events Analysis

Recurrent-event outcomes

are events that a single patient can experience more than once over follow-up: COPD/asthma exacerbations, heart-failure hospitalizations, sickle-cell pain crises, infections under immunosuppression, severe hypoglycemia, falls, seizures, and most healthcare-resource-utilization (HCRU) endpoints. The naive habit of analyzing only time to first event throws away the majority of the information and, more importantly, answers a different question than payers and clinicians are asking. A drug that does not delay the first exacerbation but halves the long-run exacerbation rate will look null in a first-event Cox model and clearly beneficial in a recurrent-event model. The choice of method is therefore an estimand decision, not a software preference.

Core estimand distinction

Recurrent-event analyses split into three estimand families that are NOT interchangeable.

1. Rate / intensity: how frequently events occur per unit person-time. The marginal rate (Lin-Wei-Yang-Ying, LWYY) targets the population-averaged event rate and is robust to unspecified within-person dependence; the intensity (Andersen-Gill, AG) conditions on the prior event history through the risk set. Both yield a rate/hazard ratio.
2. Gap-time / conditional: time between successive events, modeled with Prentice-Williams-Peterson (PWP) stratified by event number - this answers "given you have had k events, what is the effect on time to the (k+1)th?" and is appropriate when biological risk genuinely changes after each event (e.g., post-MI).
3. Absolute burden: the mean cumulative function (MCF) / expected cumulative number of events by time t (Nelson-Aalen-type estimator, or the Ghosh-Lin/Cook-Lawless MCF in the presence of death), which is the most communicable quantity for clinical and HTA audiences because it is on the natural scale of "events per patient." Crucially, AG and LWYY share the same point estimate of the rate ratio but differ in variance: LWYY uses a robust sandwich variance that does not require the AG independent-increments assumption, which is almost always violated in chronic disease (patients with one exacerbation are prone to more). For most RWE rate questions, LWYY (or a negative-binomial rate model) is the defensible default; reserve AG for when the event-history dependence is itself of interest, and PWP for ordered gap-time questions.

Pros, cons, and trade-offs

• vs time-to-first-event Cox (cox-ph-regression): Recurrent-event models use the entire event process and match burden/HCRU estimands; first-event Cox discards all subsequent events and can be badly underpowered or even sign-wrong when the treatment acts on rate rather than on time-to-first. Cost: more data engineering (counting-process intervals), a harder-to-explain effect measure, and explicit assumptions about within-person dependence and terminal events. Prefer recurrent-event methods whenever the second-and-later events carry clinical or economic weight.
• vs simple Poisson/negative-binomial count models (poisson-negative-binomial-count-models): A negative-binomial rate model with a log person-time offset IS a recurrent-event method and is often the right first choice for total burden - it absorbs overdispersion (the empirical signature of within-person clustering) and gives a clean rate ratio. Its limitation is that it collapses the process to a single count and cannot represent event timing, time-varying exposure, or the shape of risk over follow-up. Prefer AG/LWYY/PWP when timing, time-updated covariates, or the risk trajectory matter; prefer negative binomial for a transparent, payer-friendly total-burden summary.
• vs composite "first hospitalization or death" endpoints: The composite forces a single first event and treats a hospitalization as equivalent to death; recurrent-event-with-terminal-event methods (joint frailty, Ghosh-Lin, while-alive estimands) keep them distinct and let death act as the informative truncation it is. Cost: more complex modeling and stronger reliance on a complete death source.
• vs negative-binomial without competing-risk handling: Ignoring death is the central trap. If the more effective arm keeps frailer patients alive longer, those survivors accrue more events, and a naive rate model can make the better drug look worse. Prefer joint frailty or a while-alive / MCF-with-death estimand whenever mortality is non-trivial and plausibly differential by arm.

When NOT to use — and when it is actively misleading or dangerous

• The outcome is genuinely a first/terminal event (death, first stroke as a one-time terminal endpoint, first MI in a primary-prevention question where you only care about onset). Forcing a recurrent-event frame here invents a process that does not exist.
• Death is common and differential by arm but you use a naive rate/AG model. This is the dangerous case: the informative terminal event biases the rate ratio, often toward harm for the more effective drug via the survivor-accrual mechanism above. If you cannot model death jointly, you must at minimum present a while-alive estimand or restrict to a window where mortality is negligible - and say so.
• Within-person dependence is ignored. Using ordinary (model-based) Cox or Poisson standard errors on stacked intervals understates variance because events within a person are correlated; robust/sandwich variance (cluster on person_id) or a frailty term is mandatory, not optional.
• Events cannot be cleanly delimited. If your data source cannot separate one episode from its own follow-up claims (e.g., a hospitalization plus its readmission transfer, or a 30-day steroid taper recorded as daily fills), the "event count" is an artifact of coding, and any rate ratio inherits that artifact.

Data-source operational depth

• Claims (FFS): The workhorse substrate, but every event must be built from raw claims into clinical episodes. A single hospitalization generates a facility claim plus multiple professional and DME claims with different service dates; collapse them, and stitch inter-facility transfers (discharge-to-admit gap of 0-1 days) into one event or you will count a transfer as a "recurrent" admission. Apply a setting-specific clean window (e.g., a moderate COPD exacerbation requires a steroid/antibiotic burst with no qualifying event in the prior 14 days) so that the medication refills sustaining one episode are not read as new events. Person-time (the offset) must end exactly at disenrollment, death, or study end - extending it past disenrollment fabricates exposure with no observable events and dilutes the rate.
• Claims (Medicare Advantage): MA encounter data are notoriously incomplete and inconsistently submitted across plans, so events are differentially undercounted; restrict rate analyses to fee-for-service (Parts A/B with Part D for drug-defined events) and exclude MA-only person-time, or recurrence rates will be biased downward by missingness rather than by true clinical benefit.
• Competing risks in elderly claims: In older or sicker cohorts, death rates differ by exposure, so the censoring of the recurrent process is informative AND differential. A rate model that treats death as ordinary administrative censoring will mis-rank the arms. Carry a reliable mortality source (Medicare vital status, the limited Part D death flag, or linked NDI) and use a death-aware estimand.
• Immortal time in procedure/initiation studies: If follow-up (and thus the at-risk person-time for recurrence) starts at a landmark that the patient had to survive event-free to reach (e.g., counting readmissions only among those who survived the index surgery, with time zero set at discharge but exposure defined post-discharge), the interval before exposure is immortal and inflates the comparator's apparent event-free time. Align time zero to the exposure decision and start counting events from there.
• EHR: Events are encounters, labs, vitals, notes, or rescue orders. Visit frequency is itself outcome-correlated (sicker patients visit more), so raw event counts confound disease severity with capture intensity - this is informative observation/observation bias. Restrict to unambiguous severe events, model the visit process, or use inverse-intensity-of-observation weighting; never treat "more recorded events" as "more true events."
• Registry / linked: Registries give adjudicated, clean events (the numerator) but usually miss out-of-registry utilization and may have irregular assessment schedules that create panel-count rather than exact-time data. Link to claims for complete hospitalization burden and to a death index for the terminal event; reconcile registry visit dates against claim service dates before building intervals.

Worked claims example

Question: does drug A vs active comparator B reduce the rate of moderate-or-severe COPD exacerbations among new initiators in a 100% Medicare FFS sample (Parts A/B/D)?

1. Eligibility: age >=40, >=2 COPD diagnoses (J44.x), and 365 days of continuous A/B/D enrollment before index_date (first qualifying fill of A or B; new users of both).
2. Event definition: a severe exacerbation = an inpatient or ED claim with a COPD principal/first diagnosis; a moderate exacerbation = an outpatient/ED claim for COPD accompanied by a Part D fill of a systemic corticosteroid and/or a COPD-relevant antibiotic within +-5 days.
3. Episode cleaning: collapse facility + professional claims with overlapping or adjacent service dates into one event; stitch transfers (admit within 1 day of a prior discharge) into the same episode; impose a 14-day clean window so the steroid taper sustaining one episode and any immediate follow-up visit are not counted as new events.
4. Counting-process layout: for each person_id build start-stop rows from index_date with one row ending at each cleaned event date (event=1) and a final administrative row (event=0) ending at the earliest of disenrollment, death, or study end; carry days_supply-derived on-treatment status as a time-varying covariate if an as-treated estimand is wanted.
5. Person-time = sum of (tstop-tstart); it must terminate at the FFS-observable end, never at end-of-data for a patient who left FFS earlier.
6. Estimation: headline result = MCF of exacerbations by month, by arm, accounting for death (Ghosh-Lin); adjusted effect = LWYY marginal rate model (or a negative-binomial rate model with log person-time offset) with robust variance clustered on person_id and a high-dimensional propensity score or PS weights.
7. Because COPD patients with severe disease both exacerbate and die more, fit a joint frailty model or present a while-alive exacerbation rate as the primary death-aware sensitivity analysis, and report attrition at every step.

Interpreting the output

An LWYY marginal rate model of COPD exacerbations returns: rate ratio = 0.73 (95% CI 0.58–0.92) for treated vs untreated, with robust sandwich variance clustered on patient.

Formal interpretation. The rate ratio of 0.73 estimates that the marginal mean exacerbation rate in the treated group is approximately 73% of the rate in the untreated group, averaged over all patients and follow-up time. The robust variance accounts for within-person correlation — each patient's events are not independent, so standard Poisson or Cox standard errors would be anti-conservative. Terminal events (death) create informative censoring: patients who die can no longer exacerbate, compressing observed rates in higher-mortality groups. This is why a while-alive rate or joint frailty model is reported alongside as a sensitivity analysis. The rate ratio summarizes the full event burden across follow-up, not merely the time to first event.

Practical interpretation. A rate ratio of 0.73 means the treated group experiences roughly 27% fewer exacerbations per unit of follow-up time — a burden-of-disease summary more policy-relevant than a time-to-first-event HR when the disease is relapsing-remitting. Pair this with the mean cumulative function (MCF) plot by arm to visualize diverging event burden over time and confirm the rate ratio is not driven by a single early episode or differential dropout rather than a sustained reduction in recurrence.

Decision diagram

flowchart TD
  Q{What is the recurrent-event estimand?} -->|Total burden / rate| RATE[Marginal rate: LWYY or negative-binomial<br/>log person-time offset, robust variance]
  Q -->|Effect conditional on event history| AG[Andersen-Gill intensity<br/>counting-process, time-varying exposure]
  Q -->|Effect on time between ordered events| PWP[Prentice-Williams-Peterson<br/>gap-time, strata by event number]
  Q -->|Communicable absolute burden| MCF[Mean cumulative function<br/>expected events by time]
  RATE --> DEATH{Is death common and<br/>plausibly differential by arm?}
  AG --> DEATH
  PWP --> DEATH
  MCF --> DEATH
  DEATH -->|Yes| JOINT[Joint frailty / Ghosh-Lin /<br/>while-alive estimand]
  DEATH -->|No| REPORT[Report rate ratio + MCF<br/>with attrition and sensitivity]
  JOINT --> REPORT

gantt
  title Counting-process layout for one patient (COPD exacerbations, claims)
  dateFormat YYYY-MM-DD
  axisFormat %b %Y
  section Baseline
  365-day continuous FFS enrollment + washout :done, base, 2023-01-01, 2023-12-31
  section At-risk follow-up (start-stop rows)
  Interval 1 -> event at exacerbation 1 :active, i1, 2024-01-01, 70d
  Interval 2 -> event at exacerbation 2 :active, i2, 2024-03-11, 120d
  Interval 3 -> admin censor (disenroll / death / study end) :crit, i3, 2024-07-09, 90d

Worked example

Scenario

Patient 3041 has COPD and is followed for one full calendar year (January 1 through December 31, 2024) after starting a new inhaler. During that year the patient has three moderate-to-severe COPD exacerbations. We want to know the patient's annual exacerbation rate. A first-event-only analysis would record only the March flare-up and then stop watching — yielding a single event. The recurrent-events approach keeps watching and captures all three, giving an event rate of 3.0 exacerbations per person-year.

Dataset

Counting-process layout for patient 3041: one row per at-risk interval, each ending at the next exacerbation or the administrative end of follow-up. This is the format the Andersen-Gill model actually reads.

FIG. 1 — DESIGN TIMELINE

Steps

Result

3 exacerbations in 365 days of follow-up = 3.0 events per person-year. A first-event analysis using only the day-70 event would see 1 event in 0.19 person-years and miss two-thirds of this patient's disease burden.

Trade-offs

Runnable example

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Offset = log of observable person-time (in years here so the rate is events/person-year).
df["log_pt_years"] = np.log(df["pt_days"] / 365.25)

nb = smf.glm(
    formula="n_events ~ arm + age + cci",
    data=df,
    family=sm.families.NegativeBinomial(),      # absorbs within-person overdispersion
    offset=df["log_pt_years"],
    freq_weights=df["ps_weight"] if "ps_weight" in df else None,
).fit(cov_type="cluster", cov_kwds={"groups": df["person_id"]})

rr = np.exp(nb.params["arm"])
ci = np.exp(nb.conf_int().loc["arm"])
print(f"Rate ratio (arm) = {rr:.3f}  95% CI [{ci[0]:.3f}, {ci[1]:.3f}]")

from lifelines import CoxTimeVaryingFitter

ctv = CoxTimeVaryingFitter()
ctv.fit(
    long,
    id_col="person_id",
    start_col="tstart",
    stop_col="tstop",
    event_col="event",
    formula="arm + on_treatment + age + cci",
    robust=True,            # sandwich variance for within-person dependence (AG -> LWYY-style SE)
)
ctv.print_summary()         # exp(coef) on 'arm' is the recurrent-event rate/intensity ratio

library(survival)

# Andersen-Gill intensity model; cluster(id) -> robust (LWYY) variance, the RWE default for the rate ratio.
ag <- coxph(Surv(tstart, tstop, event) ~ arm + age + cci + cluster(id),
            data = long_df)

# Prentice-Williams-Peterson GAP-TIME: clock resets each event (use gap = tstop - tstart), stratified by event order.
long_df$gap <- long_df$tstop - long_df$tstart
pwp <- coxph(Surv(gap, event) ~ arm + age + cci + strata(event_number) + cluster(id),
             data = long_df)

# Mean cumulative function (Nelson-Aalen-type) of expected events by time, by arm, for the communicable headline.
mcf <- survfit(Surv(tstart, tstop, event) ~ arm, data = long_df, id = id)
summary(ag); summary(pwp); plot(mcf, cumhaz = TRUE, xlab = "Days", ylab = "Mean cumulative events")

/* Negative-binomial RATE model with a log person-time offset -> rate ratio = exp(arm estimate). */
proc genmod data=work.summary;
  class arm (ref='0') / param=ref;
  model n_events = arm age cci / dist=negbin link=log offset=log_pt;
  repeated subject=person_id / type=ind;   /* robust SE for residual within-person clustering */
  estimate 'Rate ratio (arm)' arm 1 -1 / exp;
run;

/* Andersen-Gill counting-process Cox; covs(aggregate)+id = robust (LWYY) variance for the rate/intensity ratio. */
proc phreg data=work.long covs(aggregate);
  class arm (ref='0') / param=ref;
  id person_id;
  model (tstart, tstop)*event(0) = arm age cci / rl;
run;

/* Prentice-Williams-Peterson: stratify the risk set by event number; gap-time uses tstop-tstart as the response. */
data work.long; set work.long; gap = tstop - tstart; run;
proc phreg data=work.long covs(aggregate);
  class arm (ref='0') / param=ref;
  id person_id;
  strata event_number;
  model gap*event(0) = arm age cci / rl;
run;

Citations

1. [1]Andersen PK, Gill RD. Cox's regression model for counting processes: a large sample study. Annals of Statistics. 1982;10(4):1100-1120.
2. [2]Amorim LDAF, Cai J. Modelling recurrent events: a tutorial for analysis in epidemiology. International Journal of Epidemiology. 2015;44(1):324-333.

1. [3]Prentice RL, Williams BJ, Peterson AV. On the regression analysis of multivariate failure time data. Biometrika. 1981;68(2):373-379.

1. [4]Lin DY, Wei LJ, Yang I, Ying Z. Semiparametric regression for the mean and rate functions of recurrent events. Journal of the Royal Statistical Society Series B. 2000;62(4):711-730.