Recurrent Events: Counting Every Admission Without Pretending They Are Independent

อ่านฉบับภาษาไทย (Thai version)

Abstract

Patients with heart failure are admitted again and again, yet many trials analyse only the time to the first admission. Recurrent-event models keep every admission, using a start-stop layout with one row per interval at risk. The Andersen-Gill model fits one Cox model to all admissions on a single clock from randomisation, with a standard error clustered on patient because one patient's admissions are correlated. The Prentice-Williams-Peterson model stratifies by admission number, with time measured from randomisation (total time) or from the previous admission (gap time). In a simulated trial of 1,500 patients with 1,313 admissions, the Andersen-Gill hazard ratio, comparing admission rates between arms, was 0.73 (95% CI 0.63 to 0.85). Clustering raised the standard error of its logarithm from 0.0555 to 0.0738 without moving the estimate. Death ends the count, so treating it as ordinary censoring (follow-up ending without the event) changes the question. This article shows how to build the rows, choose the clock and report deaths alongside admissions.

🔍 Click to enlarge Click anywhere to close

A trial that counts only the first admission

A steering committee is reviewing the analysis plan of a randomised heart-failure trial. The trial enrols 1,500 adults after a heart-failure admission and follows them for 24 to 36 months. The drug is fictional and the data are simulated.

The primary outcome is the time to the first heart-failure admission. Over follow-up, some patients are admitted again and again, a few of them four or more times, and most are never admitted. A cardiologist on the committee asks how to use every admission without treating one patient's repeated admissions as if they came from different people.

Admissions are recurrent events: events that can happen more than once in the same patient. The answer has three parts. Lay the data out with one row per interval at risk, and count every admission with a standard error that respects the patient. Then report deaths alongside, because a patient who has died can no longer be admitted.

What a first-event analysis leaves out

Take four patients followed for up to 24 months, a hand example rather than trial data. R1 was admitted at months 3 and 8 and was alive and in follow-up at month 24. R2 was admitted at month 5 and died at month 12. R3 was never admitted and was followed to month 24.

R4 was admitted at months 2, 4 and 10 and was lost to follow-up at month 18. The four patients had 6 admissions in all. An analysis of the time to the first event keeps one admission per patient: 3 of the 6.

Censoring means that follow-up ends without the event, so the patient is known to be event-free only up to that time. In a first-event analysis, R1, R2 and R4 leave at their first admission, and R3 is censored at month 24. R1's second admission and R4's second and third never enter the analysis. A patient admitted three times looks the same as one admitted once.

One row per interval at risk

Recurrent-event models take a counting process view: each patient carries a running count of admissions that steps up by one at every admission. The data go into a start-stop layout, with one row for each interval in which the patient is at risk of the next admission.

Each row has five entries. start and stop are the months since randomisation at which the interval opens and closes. event is 1 if the interval ends in an admission and 0 if it ends in censoring or death. enum is the number of the admission the patient is at risk for, and gap is stop minus start.

A patient with $k$ admissions therefore has $k + 1$ rows, the last ending at death or censoring. In the simulated trial analysed below, 1,500 patients with 1,313 admissions give 2,813 rows.

The ten start-stop rows of the hand example

Andersen-Gill: every admission on one clock

The Andersen-Gill model fits a Cox model to the start-stop rows, so every admission counts toward one hazard ratio (HR) [1].

Here $h_i(t)$ is the admission rate at month $t$ since randomisation for a patient $i$ who is still at risk, and $A_i$ is 1 for treatment and 0 for control. $h_0(t)$ is the control-arm rate, left unspecified as in any Cox model, $\beta$ is the log hazard ratio and $\exp(\beta)$ is the hazard ratio. One $h_0(t)$ serves the first, second and later admissions alike, so admissions are treated as interchangeable.

Time runs on one clock from randomisation, called total time. At each moment the risk set, the patients who could have the event then, holds every patient alive and in follow-up. When R4 is admitted for the second time at month 4, all four hand-example patients are in it.

Why the standard error is clustered on patient

Admissions from one patient are not independent. Some patients are more prone to admission than others and are admitted more often throughout, so one admission predicts more. A model-based, or naive, standard error treats the 1,313 admissions as if they came from 1,313 unrelated people. With positive dependence, as here, it comes out too small.

The robust (sandwich) standard error first adds up everything each patient contributes to the estimate. It then measures how much these patient totals vary, which keeps it valid when admissions cluster within patients. Clustering leaves the hazard ratio itself unchanged, and that estimate can be read as a ratio of average admission rates between the arms among patients still alive and in follow-up. The robust standard error keeps the interval valid under that reading, without assuming that one patient's admissions are independent [2].

A sandwich standard error repairs the standard error when the variance assumption is wrong. It does not repair a wrong mean model, and a biased coefficient keeps its bias.

In the simulated trial the Andersen-Gill hazard ratio is 0.73. The naive standard error of its logarithm, 0.0555, would give a 95% CI of 0.66 to 0.82. Clustering on patient raises the standard error to 0.0738 and widens the interval to 0.63 to 0.85. The estimate does not move; only its uncertainty does.

• Stata: stset stop, id(id) failure(event == 1) time0(start) exit(time .), then stcox arm, vce(cluster id)
• R: coxph(Surv(start, stop, event) ~ arm + cluster(id), data = rl, ties = "breslow"), where rl holds the start-stop rows

Prentice-Williams-Peterson: admissions in order

The Prentice-Williams-Peterson (PWP) model respects the order of admissions [3]. It stratifies the Cox model by admission number, so each admission number gets its own baseline rate. A patient enters the risk set for admission $k$ only after admission $k - 1$.

Here $k$ is the admission number, which defines the stratum, and $h_{0k}(t)$ is the control-arm rate for admission $k$ at time $t$ on the chosen clock. The hazard ratio $\exp(\beta)$ is shared across strata unless one is fitted per stratum. It compares the arms among patients with the same number of previous admissions.

PWP runs on either of two clocks. In total time (PWP-TT), $t$ counts months since randomisation and each row keeps its start and stop. In gap time (PWP-GT), $t$ counts months since the previous admission: the clock restarts at 0 after each admission, so each row runs from 0 to its gap.

In the hand example the four strata hold 4, 3, 2 and 1 rows. At month 4 in total time, only R1 and R4 are at risk of a second admission: R2 has not yet had its first, and R3 never has one. In gap time, R4's second admission comes 2 months after its first. The risk set is then every patient's interval at risk for a second admission that lasted at least 2 months: R1 (5 months), R2 (7) and R4 (2).

Total time suits a question about time since randomisation, for example when a treatment effect may build up or fade during the trial. Gap time suits a process that starts afresh after each admission, where the question is how long until the next one. The clock is best chosen in the analysis plan, before the data are seen.

Later admissions, selected patients

Later strata are small and selected. Few patients reach them, so the simulated analysis pools the fourth and later admissions into one stratum. Patients who reach a later admission also tend to be more prone to admission, and within that stratum the arms are no longer the randomised groups.

PWP-TT gives HR 0.78 (95% CI 0.69 to 0.87) and PWP-GT gives 0.77 (0.69 to 0.87), both with the clustered standard error. In a simulated trial of 100,000 patients both settle near 0.76, not 0.70. Treated patients who reach a second admission are, on average, more prone to admission than control patients who do, and that pulls the comparison toward 1. Clustering changes the PWP standard error far less than the Andersen-Gill one, probably because conditioning on the number of previous admissions already absorbs much of the dependence.

• Stata: stcox arm, strata(enum_c) vce(cluster id), after the start-stop stset for total time and after stset gap, failure(event == 1) for gap time
• R: coxph(Surv(start, stop, event) ~ arm + strata(enum_c) + cluster(id), data = rl, ties = "breslow"), with Surv(gap, event) for gap time; enum_c is the admission number with four and above pooled

Two other models: Wei-Lin-Weissfeld and shared frailty

The Wei-Lin-Weissfeld (WLW) model treats the first, second and third admissions as separate outcomes, each timed from randomisation [4]. Every patient counts as at risk for every admission number from the start, even for a third admission before a second. Because the risk set for a later admission does not wait for the earlier ones, its hazard ratio also reflects the effect on those earlier admissions, so WLW suits several distinct event types better than ordered admissions of one kind.

A shared frailty model gives each patient a random multiplier, the frailty, that scales all of that patient's admission rates. Its hazard ratio compares a treated and an untreated patient with the same frailty. The simulated trial was generated this way, with a frailty of mean one and variance 0.6 drawn from a gamma distribution (a distribution of positive values) and a rate ratio of 0.70.

Because the frailty is independent of arm, death and censoring, patients still in follow-up have the same average frailty in both arms at every month. The patient-level and population-level rate ratios are therefore both 0.70 here. They can differ when death depends on frailty, or when the risk set is limited to patients with a given history, as in a first-event or PWP analysis.

Death stops the count

Death is a terminal event: once a patient dies, no further admission can happen. Andersen-Gill and PWP code the last row of a patient who died as event = 0, like censoring, and the count models below simply end that patient's follow-up at death. Read as censoring, that assumes the patients who died would have been admitted at the same rate as similar patients still alive. In heart failure that is rarely credible, because the patients who die are often the sickest, and the sickest are admitted most.

Death therefore cannot be treated as independent censoring. Coding it as event = 0 is a convention of the layout that fixes the question: the hazard ratio describes the admission rate among patients still alive. A drug that left that rate unchanged but shortened life would still lower the admissions per randomised patient, because its patients would spend less time alive.

In the simulated trial the drug lowers the hazard of cardiovascular death (hazard ratio 0.70), so treated patients live longer and have more time to be admitted. With 100,000 simulated patients, the mean number of admissions per patient is 0.93 with control and 0.71 with treatment. That is a smaller relative difference than the rate ratio of 0.70. Death was generated independently of the frailty, which is why the Andersen-Gill value with 100,000 simulated patients is 0.70; real data give no such guarantee.

Three approaches keep death in view. A composite outcome counts cardiovascular deaths as events alongside admissions, so each such death counts against the drug, although one death can still remove several later admissions. The mean frequency function estimates the expected number of admissions per patient by each month, with patients who die adding no more admissions, so it is read beside the survival curves for death [5]. Joint models, which link admissions and death through a shared frailty, go further and are beyond this article.

For the first admission, death before admission is a competing event, one that stops the event of interest from ever happening. One minus the Kaplan-Meier curve then estimates the proportion who would be admitted if no one could die first, assuming that removing death would not change the admission rate: a hypothetical quantity. The cumulative incidence function gives the proportion actually admitted by a given month when death can come first; Part 4 of this series covers both.

Counting instead of timing: negative binomial regression

A simpler route ignores timing and counts each patient's admissions. Poisson regression models the count with follow-up time as the exposure, so it compares admissions per month of follow-up. The logarithm of follow-up enters as an offset, a term whose coefficient is fixed at one. Poisson regression assumes that the variance of a count equals its mean.

Frailty breaks that assumption, because a few patients have many admissions and most have none; this extra variation is called overdispersion. Negative binomial regression adds a patient-level multiplier with variance $\alpha$ that absorbs it. In the simulated trial $\hat\alpha$ is 0.67, where the hat marks an estimate, close to the generator's frailty variance of 0.6. The negative binomial rate ratio is 0.72 (95% CI 0.62 to 0.83).

The Poisson rate ratio is 0.73, and its naive and robust standard errors, 0.0555 and 0.0738, equal those of Andersen-Gill to four decimals. With one row per patient, the robust standard error of a Poisson model already respects the patient. Count models give a rate ratio and a mean count, not the timing of admissions, and they handle death no better than the Cox-type models. A separate article on Poisson regression covers rates in more depth.

• Stata: nbreg n_adm arm, exposure(fu_months) irr, where n_adm is the count and fu_months the follow-up in months
• R: glm.nb(n_adm ~ arm + offset(log(fu_months)), data = w3) from the MASS package

The simulated trial: every model side by side

The simulated trial has 1,500 patients, 750 per arm, and 1,313 admissions: 727 with control and 586 with treatment. At least one admission occurred in 364 control and 319 treated patients. The mean number of admissions per patient is 0.97 with control and 0.78 with treatment. The generator's admission rate ratio is 0.70.

A Cox model for the first admission only, with deaths before admission censored, uses 683 of the 1,313 admissions. It gives HR 0.77 (95% CI 0.66 to 0.89), and its value in a simulated trial of 100,000 patients is 0.75, not 0.70. This is the cause-specific Cox model of Part 5 of this series.

The population-level first-admission hazard ratio drifts from 0.70 at randomisation to 0.76 after one year and 0.80 after two years, while the patient-level ratio stays at 0.70. A single first-admission hazard ratio therefore averages a ratio that is not constant over follow-up. The frailest control patients are admitted early and leave the risk set faster than the frailest treated patients. So the control patients who remain are less prone to admission than the treated patients who remain.

Every model on the same simulated admissions

Stata: Andersen-Gill with a standard error clustered on patient

* Andersen-Gill: every admission counts, total time, robust (sandwich) SE clustered on patient
stset stop, id(id) failure(event == 1) time0(start) exit(time .)
stcox arm, vce(cluster id)

. * Andersen-Gill: every admission counts, total time, robust (sandwich) SE clustered on patient
. stset stop, id(id) failure(event == 1) time0(start) exit(time .)

Survival-time data settings

           ID variable: id
         Failure event: event==1
Observed time interval: (start, stop]
     Exit on or before: time .

--------------------------------------------------------------------------
      2,813  total observations
          0  exclusions
--------------------------------------------------------------------------
      2,813  observations remaining, representing
      1,500  subjects
      1,313  failures in multiple-failure-per-subject data
 29,562.808  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =  35.94684

. stcox arm, vce(cluster id)

         Failure _d: event==1
   Analysis time _t: stop
  Exit on or before: time .
        ID variable: id

Iteration 0:  Log pseudolikelihood = -9007.5576
Iteration 1:  Log pseudolikelihood = -8991.7416
Iteration 2:  Log pseudolikelihood = -8991.7416
Refining estimates:
Iteration 0:  Log pseudolikelihood = -8991.7416

Cox regression with no ties

No. of subjects =       1,500                           Number of obs =  2,813
No. of failures =       1,313
Time at risk    = 29,562.8078
                                                        Wald chi2(1)  =  17.80
Log pseudolikelihood = -8991.7416                       Prob > chi2   = 0.0000

                                 (Std. err. adjusted for 1,500 clusters in id)
------------------------------------------------------------------------------
             |               Robust
          _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         arm |      .7325    .054048    -4.22   0.000     .6338715    .8464749
------------------------------------------------------------------------------

R: the same model, naive and robust standard errors in one table

# Andersen-Gill: every admission counts, total time, robust (sandwich) SE clustered on patient
ag <- coxph(Surv(start, stop, event) ~ arm + cluster(id), data = rl, ties = "breslow")
print(summary(ag))

> ag <- coxph(Surv(start, stop, event) ~ arm + cluster(id),
+     data = rl, ties = "breslow")

> print(summary(ag))
Call:
coxph(formula = Surv(start, stop, event) ~ arm, data = rl, ties = "breslow",
    cluster = id)

  n= 2813, number of events= 1313

        coef exp(coef) se(coef) robust se     z Pr(>|z|)
arm -0.31129   0.73250  0.05555   0.07376 -4.22 2.44e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    exp(coef) exp(-coef) lower .95 upper .95
arm    0.7325      1.365    0.6339    0.8464

Concordance= 0.539  (se = 0.009 )
Likelihood ratio test= 31.63  on 1 df,   p=2e-08
Wald test            = 17.81  on 1 df,   p=2e-05
Score (logrank) test = 31.66  on 1 df,   p=2e-08,   Robust = 17.59  p=3e-05

  (Note: the likelihood ratio and score tests assume independence of
     observations within a cluster, the Wald and robust score tests do not).

Recurrent-event models compared

Choose by the question: Andersen-Gill for the admission rate among patients alive, PWP for the next admission given the number of earlier admissions, count models for admissions per month of follow-up. The table below lines up the models by clock, risk set and meaning. A published tutorial compares these models on epidemiological data [6], and a review applies them to recurrent heart-failure admissions in a large trial [7].

What to do in your own analysis

• Write the estimand, the precise quantity the analysis aims to estimate, before choosing a model: the admission rate among patients alive, the expected number of admissions per patient, or a composite with death.
• Build the start-stop rows from the admission dates, and check a few patients by hand against their records.
• When one patient contributes more than one event to Andersen-Gill, PWP or Wei-Lin-Weissfeld, cluster the standard error on patient; a frailty model instead handles the dependence through its random term.
• If the order of admissions matters, consider prespecifying PWP and its clock, and pooling sparse later strata.
• Consider deciding in advance whether days in hospital count as time at risk, since an inpatient cannot be admitted again.
• Report deaths by arm beside every admission result.

Related in the wiki: [[time-to-event-survival-analysis]]