survrec: survival analysis for recurrent event data
In survrec: Survival Analysis for Recurrent Event Data

knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", message = FALSE,
  warning = FALSE, fig.width = 6.2, fig.height = 5,
  fig.align = "center"
)
set.seed(1)

Introduction

In many follow-up studies the event of interest is not terminal: a patient can be rehospitalized several times, a tumour can relapse, a machine fails and is repaired. The natural response is then the gap time, the time elapsed between consecutive occurrences, and each subject contributes as many gap times as observed events plus one final censored gap (from the last event to the end of follow-up).

Estimating the common survival function of the gap times is not a matter of pooling them into a Kaplan-Meier estimator: the last gap of every subject is always censored, longer gaps are more likely to be cut by the end of the study, and gap times within a subject may be correlated. survrec implements three estimators designed for this setting:

PSH, the generalized product-limit estimator of @pena2001, valid when the gap times of a subject are independent and identically distributed (a renewal process);
Wang-Chang, the estimator of @wang1999, which remains consistent when gap times are correlated within subjects;
MLE frailty, maximum likelihood under a gamma frailty model [@pena2001], where a subject-specific random effect \eqn{Z_i \sim \Gamma(\alpha, \alpha)} induces the correlation and the marginal survival is \eqn{S(t) = [\alpha/(\alpha + \Lambda_0(t))]^\alpha}.

The comparison of the three estimates on the same data is itself informative: agreement between PSH and the frailty estimate suggests independent gaps (large \eqn{\alpha}), while a Wang-Chang curve separated from the PSH one points to within-subject correlation.

library(survrec)

The data

We use the two datasets shipped with the package. MMC contains the times of the migratory motor complex, a cyclic intestinal motility pattern, for 19 healthy individuals [@husebye1990]:

data(MMC)
head(MMC)

colon records rehospitalizations after surgery in 403 colorectal cancer patients, with Dukes stage and chemotherapy as covariates

data(colon)
head(colon)

The response of every function in the package is a Survr object built from the subject identifier, the gap times and the event indicator (1 = event, 0 = the final censored gap of each subject):

x <- Survr(MMC$id, MMC$time, MMC$event)

Estimating the survival of the gap times

survfitr() is the formula interface; type selects the estimator:

fit <- survfitr(Survr(id, time, event) ~ 1, data = MMC, type = "wang-chang")
fit

summary() returns the estimated curve, and quantile() the survival times of any set of quantiles:

quantile(fit, probs = c(0.25, 0.5, 0.75))

The fitted curve can be drawn with autoplot() (a ggplot2 graphic; the classic base-graphics plot() method is also kept). Pointwise confidence bands use the log-minus-log transformation by default, so they always stay inside \eqn{[0, 1]}:

autoplot(fit)

Comparing the three estimators

plotEstimators() fits and overlays the three estimators:

plotEstimators(x)

The three curves are close, which suggests little within-subject correlation. The frailty fit makes this quantitative through \eqn{\alpha} (the larger \eqn{\alpha}, the weaker the association between gap times of a subject):

mle <- mlefrailty_fit(x, alpha.console = FALSE)
mle$alpha

The posterior frailty estimates of each subject are also returned; values spread away from 1 would indicate heterogeneity between subjects:

round(mle$frailties, 3)

Groups and covariates

A term on the right-hand side of the formula estimates one curve per group. For the colon data, Dukes stage:

fit.dukes <- survfitr(Survr(hc, time, event) ~ as.factor(dukes),
  data = colon, type = "pena"
)
autoplot(fit.dukes)

The estimated curves in tidy format, ready for any further processing:

head(as.data.frame(fit.dukes))

The mean cumulative function offers a complementary, calendar-time view of the process -- the expected number of events per subject up to each time point:

autoplot(mcf(Survr(hc, time, event) ~ as.factor(dukes), data = colon))

Comparing survival quantiles between groups

survdiffr() obtains bootstrap replicates of a survival quantile (the median by default) for each group. Three resampling schemes are available: nonparametric from the PSH or the Wang-Chang estimate of the gap-time distribution, and semiparametric under the fitted frailty model [@gonzalez2003]. Resampling uses R's random number generator, so a seed makes the analysis reproducible.

b <- survdiffr(Survr(hc, time, event) ~ as.factor(dukes),
  data = colon, q = 0.5, B = 199, boot.F = "WC", seed = 2026
)

summary() reports each group's observed median with a percentile bootstrap interval, and every pairwise difference with its interval and a two-sided bootstrap p-value:

summary(b)

autoplot(b)

Each group's element is a standard boot object, so the intervals of the boot package remain available:

boot::boot.ci(b$"1", type = c("norm", "basic", "perc"))

The re-estimation of the survival curve on each replicate runs in parallel when the package is compiled with OpenMP; survrecThreads() caps the number of threads (the resampling itself is serial on R's RNG, so results do not depend on the thread count).