survrec: survival analysis for recurrent event data

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)

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 Males

colon 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        1

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
#> Survival for recurrent event data
#>    n events mean se(mean) median recurrences: min max median
#>   19     80  106     12.7     95                1   9      4

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))
#> 25% 50% 75% 
#>  63  95 144

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 :

autoplot(fit)
Wang-Chang estimate of the MMC gap-time survival.

Wang-Chang estimate of the MMC gap-time survival.

Comparing the three estimators

plotEstimators() fits and overlays the three estimators:

plotEstimators(x)
The three estimators on the MMC data.

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.17569

The 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.035

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)
PSH estimates of the rehospitalization gap times by Dukes stage.

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.014032158

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))
Mean cumulative number of rehospitalizations by Dukes stage.

Mean cumulative number of rehospitalizations by Dukes stage.

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        199
autoplot(b)
Bootstrap distributions of the median rehospitalization-free time.

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 unstable

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).

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] survrec_2.0.0    BiocStyle_2.41.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3         cli_3.6.6           knitr_1.52         
#>  [4] rlang_1.3.0         xfun_0.61           otel_0.2.0         
#>  [7] S7_0.2.2            jsonlite_2.0.0      labeling_0.4.3     
#> [10] glue_1.8.1          buildtools_1.0.0    htmltools_0.5.9    
#> [13] maketools_1.3.2     sys_3.4.3           sass_0.4.10        
#> [16] scales_1.4.0        rmarkdown_2.32      grid_4.6.1         
#> [19] evaluate_1.0.5      jquerylib_0.1.4     fastmap_1.2.0      
#> [22] yaml_2.3.12         lifecycle_1.0.5     BiocManager_1.30.27
#> [25] compiler_4.6.1      RColorBrewer_1.1-3  Rcpp_1.1.2         
#> [28] farver_2.1.2        digest_0.6.39       R6_2.6.1           
#> [31] bslib_0.12.0        withr_3.0.3         gtable_0.3.6       
#> [34] tools_4.6.1         boot_1.3-32         ggplot2_4.0.3      
#> [37] cachem_1.1.0

References

González, Juan R., Esteve Fernández, Víctor Moreno, et al. 2005. “Sex Differences in Hospital Readmission Among Colorectal Cancer Patients.” Journal of Epidemiology and Community Health 59 (6): 506–11. https://doi.org/10.1136/jech.2004.028902.
González, Juan R., and Edsel A. Peña. 2003. “Bootstrapping Median Survival with Recurrent Event Data.” IX Conferencia Española de Biometría (A Coruña, Spain).
Husebye, Einar, Viggo Skar, Odd O. Aalen, and Magne Osnes. 1990. “Digital Ambulatory Manometry of the Small Intestine in Healthy Adults.” Digestive Diseases and Sciences 35: 1057–65. https://doi.org/10.1007/BF01537575.
Peña, Edsel A., Robert L. Strawderman, and Myles Hollander. 2001. “Nonparametric Estimation with Recurrent Event Data.” Journal of the American Statistical Association 96 (456): 1299–315. https://doi.org/10.1198/016214501753381922.
Wang, Mei-Cheng, and Shu-Hui Chang. 1999. “Nonparametric Estimation of a Recurrent Survival Function.” Journal of the American Statistical Association 94 (445): 146–53. https://doi.org/10.1080/01621459.1999.10473831.