To fit a Cox model to recurrent event data (Andersen & Gill), using the R survival package, requires the user to cast the data into counting process format (see [1]).

For recurrent event analyses, a subject can occur multiple times in this dataset.

I understood that this is a programming trick that allows one to fit the recurrent model by means of a cox proportional hazard model for time to first event (potentially using a sandwich variance estimator, to compensate for intra-subject-correlations in event times).

The partial likelihood for the cox model, in the setting of time to first event, can be expressed in terms of risk sets. How will the risk sets be populated in a recurrent event scenario under this programming trick? Will the counting process format be used to shift times, viewing multiple at risk periods from a subject as actually multiple subjects with overlapping at-risk periods?

So for example, assume a subject occurs 3 times in the counting process format, (0,15], (15,40] and (40,100], and at each of end of an interval an event occurred.

Will that data be fitted as three different subjects with time to first events equal to 15, 25 and 60, and at time 15, the subject effectively occurs three times in the risk set? Again, correlations are compensated for through grouped jackknife or robust sandwich variance estimation?

I am trying to understand this programmatically and intuitively, what happens with the interval data frame above, which transformations will be used to view it programmatically as a time to first event model? Why does it work like that, intuitively?

Thank you.

[1] https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf


1 Answer 1


The counting process analysis might be better described as an analysis of "times to events" rather than of "time to first event." Therneau and Grambsch describe various approaches to multiple events of the same type per subject in Section 8.5.1; quotes below are from there.

In your scenario, the time scale is time since study entry:

a subject occurs 3 times in the counting process format, (0,15], (15,40] and (40,100], and at each of end of an interval an event occurred.

For this time scale, Therneau and Grambsch explain the difference between a typical Cox survival model and an Andersen-Gill model as follows (where $Y_i(t)$ is the at-risk indicator for individual $i$ at time $t$):

The difference lies in the definition of $Y_i(t)$; for survival data, the individual ceases to be at risk when an event occurs and $Y_i$ goes to zero, but for The Andersen-Gill model for recurrent events, $Y_i(t)$ remains one as events occur.

Your subject would thus appear in the risk set for all event times up to time = 100. "This model is ideally suited to the situation of mutual independence of the observations within a subject." Robust standard errors help account for things like unmeasured outcome-associated covariates that differ among individuals.

In this Andersen-Gill model the time for an individual is not reset to 0 at each event. The (unmodeled) baseline hazard can change over time, as in a Cox model, but it remains the same for all events and extends out to the very last event time.

There are alternatives that might be more appropriate for different situations.

You can, for example, model the inter-event intervals instead of the times to events. "The model assumption in this case is that the gap times form a renewal process." That no longer requires the counting process data format, but you have to re-work data of the type you show to so that the gap time (stopTime - startTime) is the time to event or censoring for each data row for an individual.

Therneau and Grambsch describe two other ways to handle recurrent events. One treats "the ordered outcome data set as though it were an unordered competing risks problem." That typically does not require the counting process form of data, and (unlike the above approaches) "allows a separate underlying hazard for each event."

A "conditional model," in contrast, assumes that the events have a strict order:

... in general, a subject is not at risk for the $k^{th}$ event until he/she has experienced event $k - 1$. To accomplish this, the counting process style of input is used, as in the AG model, but each event is assigned to a separate stratum.

The rest of Section 8.5 illustrates these approaches.

In response to comments

Quotes again are from Therneau and Grambsch, Section 8.5.1. If analysis proceeds along continuous time from "time since study entry," without resetting time = 0 at each event for an individual, then an Andersen-Gill model

can in fact be accurately approximated with Poisson regression software in the same manner in which Laird and Olivier [J. Am. Stat. Assoc., 76:231- 240, 1981] approximate an ordinary single-event Cox model.

That's effectively time transformation using the cumulative rate function, as indicated in a comment, a transformation that is the same for all individuals for each inter-event interval.

In that situation, the underlying

assumption is equivalent to each individual counting process possessing independent increments, i.e., the numbers of events in nonoverlapping time intervals are independent, given the covariates. Such processes are typically modeled as time-varying Poisson processes, hence the Poisson connection.

The alternative time scale that resets to time = 0 at each event, "time since entry or last event," will not in general provide the same results except in the special case of exponentially distributed gap times suitable for a simple Poisson model.

The model assumption in this case is that the gap times form a renewal process. The famed lack of memory of the exponential distribution implies that a renewal process with exponential gap times is also a counting process with independent increments; but in general counting processes cannot possess both independent increments and independent gap times.

One way to think about this: any time transformation of a gap time would only apply to the individual whose data had that particular gap. There would be no way in general to apply that time transformation to other individuals to identify the corresponding risk sets at (gap-end) event times.

Finally, I don't know whether it can be shown that the pseudo-likelihood is insensitive to correlations within subject event times. There might be some literature on that, but I haven't looked for it. My understanding is that the approach with recurrent events is similar to the working independence assumption used in other marginal modeling: you assume the independence to start, and then use robust covariance estimates for regression coefficients.

  • $\begingroup$ Thank you! Just to clarify: when I’m only interested in estimating the betas in a AG model: is it possible to estimate those betas for the AG using the inter-event intervals instead of the times to events. So avoid the counting process format, instead have multiple rows per subject for each inter-event interval, use a subject ID column, and add a term +cluster(ID), to have a variance estimate, that compensates for correlations between inter event times per subject? If so, isn’t it indeed possible that a subject appears multiple times in a risk set for the PL? $\endgroup$
    – ermeel
    Commented Jan 29 at 6:13
  • $\begingroup$ Further the AG is a non-homogenous Poisson process model, right? A non-homogenous Poisson process can be made homogenous by a time transformation using the cumulative rate function. This transformation preserves the ordering of times to event times. The pseudo likelihood only depends on the order of event times, not the actual event times. So the betas could be estimated by a homogenous Poisson process model with constant rate exp(beta X). For such Poisson processes the inter event times are identically distributed with exponential distr. having rate exp(beta X). $\endgroup$
    – ermeel
    Commented Jan 29 at 6:18
  • $\begingroup$ What is missing: can it be shown that the pseudo-likelihood is insensitive to correlations within subject event times, as long as one is interested only in MLE point estimates? For variance I would use grouped jackknife. $\endgroup$
    – ermeel
    Commented Jan 29 at 6:24
  • $\begingroup$ @ermeel with respect to the first comment, you can estimate betas from the counting-process approach or from the inter-event interval approach, but the betas you get might be different in general. That's due to they different assumptions about the baseline hazard. For the inter-event interval the baseline hazard resets at each event, for the counting-process model it doesn't and it only depends on total time since study start. Therneau and Grambsch discuss matters related to the second and third comments; I'll try to summarize later when I have a chance. $\endgroup$
    – EdM
    Commented Jan 29 at 8:28
  • $\begingroup$ All of my above remarks assume time independent covariates. $\endgroup$
    – ermeel
    Commented Jan 29 at 10:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.