survrec: survival analysis for recurrent event data
Dolors Pelegrí-Sisó, Juan R. González, Edsel A. Peña and Robert L. Strawderman
2026-09-12
Abstract
survrec estimates the survival function of the time between occurrences of a recurrent event – repeated hospitalizations, tumour relapses, successive failures of a machine – from censored gap-time data. It implements the generalized product-limit estimator of Peña, Strawderman and Hollander (2001), the estimator of Wang and Chang (1999) for correlated inter-occurrence times, and maximum likelihood estimation under a gamma frailty model, together with bootstrap comparisons of survival quantiles between groups. This vignette works through a complete analysis of two real datasets.
Package
survrec 2.0.0
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 Peña et al. (2001), valid when the gap times of a subject are independent and identically distributed (a renewal process);
- Wang-Chang, the estimator of Wang and Chang (1999), which remains consistent when gap times are correlated within subjects;
- MLE frailty, maximum likelihood under a gamma frailty model (Peña et al. 2001), where a subject-specific random effect induces the correlation and the marginal survival is .
The comparison of the three estimates on the same data is itself informative: agreement between PSH and the frailty estimate suggests independent gaps (large ), while a Wang-Chang curve separated from the PSH one points to within-subject correlation.
library(survrec)2 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 (Husebye et al. 1990):
data(MMC)
head(MMC)
#> id time event group
#> 1 1 112 1 Males
#> 2 1 145 1 Males
#> 3 1 39 1 Males
#> 4 1 52 1 Males
#> 5 1 21 1 Males
#> 6 1 34 1 Malescolon records rehospitalizations after surgery in 403 colorectal
cancer patients, with Dukes stage and chemotherapy as covariates
(González et al. 2005):
data(colon)
head(colon)
#> hc time event chemoter dukes distance
#> 1 5634 24 1 2 3 1
#> 2 5634 433 1 2 3 1
#> 3 5634 580 0 2 3 1
#> 4 10767 489 1 1 2 1
#> 5 10767 693 0 1 2 1
#> 6 15843 15 1 1 2 1The 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)3 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
#> Survival for recurrent event data
#> n events mean se(mean) median recurrences: min max median
#> 19 80 106 12.7 95 1 9 4summary() returns the estimated curve, and quantile() the survival
times of any set of quantiles:
quantile(fit, probs = c(0.25, 0.5, 0.75))
#> 25% 50% 75%
#> 63 95 144The 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 :
autoplot(fit)
Figure 1: Wang-Chang estimate of the MMC gap-time survival
3.1 Comparing the three estimators
plotEstimators() fits and overlays the three estimators:
plotEstimators(x)
Figure 2: The three estimators on the MMC data
The three curves are close, which suggests little within-subject correlation. The frailty fit makes this quantitative through (the larger , the weaker the association between gap times of a subject):
mle <- mlefrailty_fit(x, alpha.console = FALSE)
mle$alpha
#> [1] 10.17569The posterior frailty estimates of each subject are also returned; values spread away from 1 would indicate heterogeneity between subjects:
round(mle$frailties, 3)
#> [1] 1.363 0.779 0.667 1.048 0.955 1.309 1.045 1.094 0.952 0.871 1.182 1.080
#> [13] 0.876 0.985 0.776 1.019 0.952 1.013 1.0354 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)
Figure 3: PSH estimates of the rehospitalization gap times by Dukes stage
The estimated curves in tidy format, ready for any further processing:
head(as.data.frame(fit.dukes))
#> group time n.event n.risk surv std.error
#> 1 1 2 1 299 0.9966555 0.003333296
#> 2 1 3 5 298 0.9799331 0.008050336
#> 3 1 4 3 293 0.9698997 0.009816343
#> 4 1 5 4 290 0.9565217 0.011714827
#> 5 1 6 2 285 0.9498093 0.012551139
#> 6 1 7 4 283 0.9363844 0.014032158The 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))
Figure 4: Mean cumulative number of rehospitalizations by Dukes stage
5 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
(González and Peña 2003). 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)
#> Bootstrap comparison of survival quantiles (95% percentile intervals)
#>
#> group estimate lower upper replicates
#> 1 1736 1188 Inf 199
#> 2 1028 487 1293 199
#> 3 199 161 335 199
#>
#> Pairwise differences:
#> contrast estimate lower upper p.value replicates
#> 1 - 2 708 159 Inf 0.0101 199
#> 1 - 3 1537 943 Inf 0.0000 199
#> 2 - 3 829 251 1123 0.0000 199autoplot(b)
Figure 5: Bootstrap distributions of the median rehospitalization-free time
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"))
#> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
#> Based on 199 bootstrap replicates
#>
#> CALL :
#> boot::boot.ci(boot.out = b$"1", type = c("norm", "basic", "perc"))
#>
#> Intervals :
#> Level Normal Basic Percentile
#> 95% ( 600, 3614 ) (1297, 3473 ) ( -1, 2175 )
#> Calculations and Intervals on Original Scale
#> Some basic intervals may be unstable
#> Some percentile intervals may be unstableThe 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).
6 Session information
sessionInfo()
#> R version 4.5.3 (2026-03-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.6.2
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: Europe/Madrid
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] survrec_2.0.0 BiocStyle_2.36.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_2.0.0 dplyr_1.2.1
#> [4] compiler_4.5.3 BiocManager_1.30.27 tinytex_0.57
#> [7] tidyselect_1.2.1 Rcpp_1.1.2 dichromat_2.0-0.1
#> [10] jquerylib_0.1.4 scales_1.4.0 boot_1.3-32
#> [13] yaml_2.3.12 fastmap_1.2.0 ggplot2_4.0.1
#> [16] R6_2.6.1 labeling_0.4.3 generics_0.1.4
#> [19] knitr_1.50 tibble_3.3.1 bookdown_0.43
#> [22] bslib_0.9.0 pillar_1.11.1 RColorBrewer_1.1-3
#> [25] rlang_1.3.0 cachem_1.1.0 xfun_0.52
#> [28] sass_0.4.10 S7_0.2.1 cli_3.6.6
#> [31] withr_3.0.3 magrittr_2.0.5 digest_0.6.39
#> [34] grid_4.5.3 rstudioapi_0.17.1 lifecycle_1.0.5
#> [37] vctrs_0.7.3 evaluate_1.0.4 glue_1.8.1
#> [40] farver_2.1.2 rmarkdown_2.29 tools_4.5.3
#> [43] pkgconfig_2.0.3 htmltools_0.5.8.1