Tutorials › Nonlinear Mixed-Effects Models (NLME) › Recurrent-Event Survival Models

Categorical, Count and Time-to-Event Responses · Tutorial 446 of 797

Recurrent-Event Survival Models

Specify recurrent-event survival models on appropriate time scales and interpret how event history and shared subject effects shape repeated-event dependence.

Advanced 12 min read

What You'll Learn

  • Define recurrent-event counting processes, risk indicators, and conditional intensities.
  • Specify Andersen–Gill models on the time-since-origin scale.
  • Formulate gap-time and event-order models and distinguish their clocks.
  • Interpret shared subject effects as a source of repeated-event dependence.
  • Construct recurrent-event likelihood contributions and check their assumptions.
  • Compare model choices and recognize common specification failures.

Why Recurrent Events Need Their Own Model

A standard time-to-event analysis often focuses on the first event. In many pharmacometric and clinical settings, however, the same subject can experience multiple events: repeated exacerbations, hospital visits, or episodes of rescue medication use. Discarding events after the first loses information; treating each event as if it came from a different person ignores the fact that events from one subject may be related.

A recurrent-event model represents the sequence of events over a subject’s observed follow-up. Its central modeling choices are the time scale used for risk, how the event number enters the intensity, and whether a subject-specific effect captures persistent differences between subjects. These choices determine the scientific meaning of the model, not just its computational form.

Definition: For subject \(i\), let \(N_i(t)\) be the number of recurrent events observed by time \(t\), and let \(Y_i(t)\) indicate whether the subject is under observation and at risk immediately before \(t\). The conditional intensity \(\lambda_i(t\mid\mathcal{F}_i(t-))\) is the instantaneous event rate at time \(t\), conditional on the observed history \(\mathcal{F}_i(t-)\) just before \(t\). It is defined so that, for a small interval \(dt\), the conditional probability of an event in \([t,t+dt)\) is approximately \(\lambda_i(t\mid\mathcal{F}_i(t-))dt\), with probability of multiple events of smaller order than \(dt\).

The event times are denoted \(T_{i1},T_{i2},\ldots\), in increasing order. Follow-up ends at a censoring time \(C_i\); events after that time are not observed. The intensity is zero when \(Y_i(t)=0\). A recurrent-event analysis therefore requires a clear definition of when a subject is observed and at risk, as well as a defensible account of censoring and any periods when events cannot be recorded.

Andersen–Gill: Time Since a Common Origin

The Andersen–Gill formulation uses a common time scale, usually time since study entry or another defined origin. At every event, the subject remains in the risk process if follow-up continues. A basic proportional-intensity specification is

$$ \lambda_i(t\mid\mathcal{F}_i(t-),b_i) = Y_i(t)\lambda_0(t) \exp\{\beta^\top X_i(t)+b_i\}. $$

Here \(\lambda_0(t)\) is the baseline intensity on the chosen time scale, \(X_i(t)\) is a covariate vector, \(\beta\) is its coefficient vector, and \(b_i\) is a subject-specific effect when included. In the basic Andersen–Gill structure, the same baseline intensity and covariate coefficients apply across event orders. Event history can still enter through \(X_i(t)\), if it is represented by appropriate history-dependent covariates, but event order is not automatically given a separate baseline.

A shared \(b_i\) represents persistent subject-level differences in event rate. Conditional on \(b_i\) and the modeled history, the intensity describes the event process for that subject. When \(b_i\) varies across subjects, integrating over that unobserved variation can induce dependence among their repeated events: subjects with higher underlying rates tend to contribute more events throughout follow-up. This is one way to represent repeated-event dependence, rather than assuming that each event is an independent observation.

Worked Example: Interpreting an Andersen–Gill Intensity

Suppose, hypothetically, the conditional baseline intensity is constant at \(0.0800\) events per month. A subject has a binary covariate \(X_i=1\), coefficient \(\beta=\ln(1.5)\), and subject effect \(b_i=\ln(1.2)\). Assume the subject remains under observation and at risk for 10 months, with no time-varying covariates.

The conditional intensity is

$$ \lambda_i = 0.0800\exp\{\ln(1.5)(1)+\ln(1.2)\} = 0.0800(1.5)(1.2) = 0.144 \text{ events per month}. $$

Thus, for this subject and these specified effects, the instantaneous event rate is 0.144 events per month. The integrated intensity over 10 months is \(0.144\times10=1.44\) events. Under this constant conditional intensity, the expected event count over the 10-month period, conditional on the covariate and subject effect, is 1.44. This is a model-based conditional mean, not a claim that the subject will experience exactly 1.44 events.

Gap-Time and Event-Order Formulations

A gap-time model resets the clock after each event. If \(U_{ij}=T_{ij}-T_{i,j-1}\) is the interval between event \(j-1\) and event \(j\), the hazard for the next event is modeled as a function of elapsed time \(u\) since the previous event. This scale is useful when the waiting time after an event has a direct scientific interpretation. It differs from Andersen–Gill time: two subjects at the same time since study entry can have different gap times because their most recent events occurred at different times.

A gap-time model must specify how event number is handled. One may use a common gap-time hazard for every event, or allow the hazard to differ by event order. Censoring before the next event contributes a censored final gap rather than an observed event time. The analysis must also define the origin of the first gap, commonly study entry or the start of the at-risk period.

Worked Example: A Gap-Time Hazard

Consider a hypothetical exponential gap-time model with conditional baseline hazard \(0.100\) per week and a binary covariate coefficient \(\beta=\ln(2)\). For a subject with \(X_i=1\), calculate the hazard and the probability that the next event-free interval exceeds three weeks, assuming the subject remains observed and at risk.

The hazard is constant over gap time and equals

$$ \lambda(u\mid X_i=1) = 0.100\exp\{\ln(2)\} = 0.200 \text{ events per week}. $$

For a constant hazard \(\lambda\), the probability of no event by gap time \(u\) is \(S(u)=\exp(-\lambda u)\). Therefore,

$$ P(U_{ij}>3\mid X_i=1) = \exp(-0.200\times3) = \exp(-0.600) \approx 0.5488. $$

The fitted model assigns approximately 54.88% probability to a gap longer than three weeks, conditional on the specified covariate and model. The hazard ratio for \(X_i=1\) versus \(X_i=0\) is \(\exp\{\ln(2)\}=2\), meaning twice the instantaneous event rate at a given gap time under this model—not half the expected gap for every possible model.

An event-order formulation allows risk to differ for the first, second, and later events. One general representation uses \(k=N_i(t-)+1\), the order of the next event, and an order-specific baseline:

$$ \lambda_i(t\mid\mathcal{F}_i(t-),b_i) = Y_i(t)\lambda_{0k}(t) \exp\{\beta_k^\top X_i(t)+b_i\}, \qquad k=N_i(t-)+1. $$

This formula uses time \(t\) since a common origin; an event-order model can instead be defined on a gap-time scale. The essential feature is that event order changes the modeled hazard, for example through \(\lambda_{0k}\), \(\beta_k\), or both. It is therefore important to state both the order structure and the clock. A subject who has not yet had event \(k\) is not at risk for the \(k\)th event-order stratum in formulations that model each order-specific risk period separately.

Worked Example: Event-Order-Specific Rates

Suppose, hypothetically, a subject has events at months 2 and 5 and is censored at month 7. For illustration, assume event-order-specific constant rates of \(0.100\), \(0.200\), and \(0.300\) per month for the first, second, and third event, respectively. The first event is observed after two months, the second after a further three months, and the subject remains event-free for two months in the third-event risk period.

For constant rates on these successive risk periods, the event-time likelihood contribution is the product of the rates at the two observed events multiplied by the probability of no third event during the final two months:

$$ L_i = (0.100)(0.200) \exp\{-(0.100\times2+0.200\times3+0.300\times2)\} = 0.0200\exp(-1.40) \approx 0.004932. $$

The likelihood has units of month\(^{-2}\), because two exact event times contribute two rate factors. The censoring at month 7 does not create an event-rate factor; it contributes the survival term for the observed event-free time in the third-event risk period. This calculation depends on the stated order-specific constant-rate model and complete observation of the event history up to censoring.

Likelihood and a Short Justification

For a simple recurrent-event process observed through \(C_i\), with no simultaneous events and a correctly specified conditional intensity, the likelihood conditional on the subject effect is proportional to the product of the intensities at observed event times multiplied by the exponential of minus the integrated intensity. This is the recurrent-event version of the event-time likelihood: each observed event contributes its instantaneous rate, while all observed event-free time contributes through the integrated rate.

$$ L_i(b_i) = \left\{\prod_{j:T_{ij}\le C_i} \lambda_i(T_{ij}\mid\mathcal{F}_i(T_{ij}-),b_i)\right\} \exp\left\{-\int_0^{C_i} \lambda_i(t\mid\mathcal{F}_i(t-),b_i)\,dt\right\}. $$

This expression presumes that the observation and censoring process is non-informative conditional on the modeled information, and that the intensity accounts for the relevant risk periods and covariate history. If \(b_i\) is random, a marginal likelihood is obtained by integrating this conditional likelihood over the distribution of \(b_i\). The fact that events share a subject does not justify multiplying separate event likelihoods while ignoring the shared effect or event history.

The basic mean relationship follows from the compensator of a counting process. If the intensity is integrable over follow-up, the difference between the event count and its integrated intensity has mean zero. Consequently,

$$ E\{N_i(\tau)\} = E\left\{\int_0^\tau \lambda_i(t\mid\mathcal{F}_i(t-))\,dt\right\}. $$

This identity is unconditional. It can also be stated conditional on suitable information available at the start of follow-up, when the corresponding conditional intensity and integrability conditions hold. It should not be conditioned on the event history through \(\tau\): that history includes the count being predicted and can change the conditional expectation. In the first worked example, the intensity was constant conditional on the specified covariate and subject effect, so the corresponding conditional mean count was the rate multiplied by the 10-month exposure.

Choosing a Formulation and Checking It

FormulationClock and event-order treatmentUseful whenKey limitation
Andersen–GillTime since common origin; typically a common baseline across event ordersRisk is naturally indexed by follow-up time and a common intensity structure is plausibleDoes not automatically allow first and later events to have different baseline hazards
Gap-timeTime since the previous event; event order may be common or modeled separatelyRecovery or waiting time after each event is scientifically centralResets the clock and may not represent risk that changes with time since study entry
Event-order-specificOrder-specific hazards; clock must be stated separatelyFirst and later events plausibly have different risk patternsLater-order strata can have sparse data and imprecise estimates

A practical specification should be checked against the event process and study design. Plot event counts and event times by follow-up time and by event order; inspect whether later events appear to follow a different pattern from first events; and verify that risk periods, event dates, and censoring dates were encoded correctly. Compare observed and model-implied event counts across follow-up and event order, and examine whether a small number of subjects account for a large share of events. These checks can reveal a poor baseline shape, an omitted event-order pattern, an incorrect risk interval, or unmodeled subject heterogeneity.

A common alternative to a subject-specific dependence model is to fit a marginal recurrent-event intensity and use a robust variance estimate to account for within-subject correlation in inference. This can protect standard errors against some forms of misspecified within-subject dependence when the marginal mean model is appropriate and there are enough independent subjects. It does not itself model subject-level heterogeneity, repair an incorrect mean intensity, or make event observations independent. Conversely, a shared subject effect models persistent heterogeneity but requires assumptions about its distribution and how it enters the intensity.

A frequent mistake is to treat every recurrent event as an independent first event, resetting the clock without stating that choice and ignoring subject membership. The correction is to define the risk process and time origin explicitly, choose a formulation that matches the scientific question, and represent within-subject dependence through the modeled history, a shared subject effect, or an appropriate marginal inference strategy. No formulation is universally best: a gap-time model can obscure follow-up-time effects, an Andersen–Gill model can conceal event-order changes, and event-order models can become unstable when few subjects reach later events.

Key takeaway: Recurrent-event analysis requires explicit choices about the intensity, the clock, event order, and dependence within subjects. Interpret every fitted rate conditional on those choices and on the stated risk and observation process.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions below.

  1. In an Andersen–Gill model, what does the baseline intensity use as its time scale, and does the basic model automatically give the second event a different baseline?
  2. A gap-time model estimates a constant hazard of \(0.0500\) per day. What is the model probability of no event during the next 4 days?
  3. Why can a shared subject effect induce dependence among repeated events, and why does a robust variance estimate not have the same role?
  4. What does censoring before the next recurrent event contribute to the event-time likelihood?

Solutions. (1) Andersen–Gill time is measured from a common origin such as study entry; its basic form uses a common baseline across event orders unless order is explicitly modeled. (2) \(S(4)=\exp(-0.0500\times4)=\exp(-0.200)\approx0.8187\), so the probability is approximately 81.87%. (3) Subjects with different shared effects have persistently different conditional intensities; averaging over that unobserved variation can make their repeated events statistically dependent. A robust variance estimate adjusts inference for within-subject correlation but does not itself specify a subject-level intensity effect. (4) The observed event-free interval contributes survival through the integrated intensity, not an event-rate factor at the censoring time.