gcmrec: modelling recurrent events with effective age
Dolors Pelegrí-Sisó, Juan R. González, Elizabeth H. Slate and Edsel A. Peña
2026-08-31
Abstract
gcmrec fits the general class of semiparametric models for recurrent event data of Peña and Hollander (2004). Unlike standard survival models, it accounts for what an event does to the subject through an effective age function, for the effect of accumulating occurrences, and for unobserved heterogeneity through gamma frailties, so it applies to settings such as repeated hospital readmissions, cancer relapses or successive failures of a machine. Estimation is carried out in C++ by profile likelihood, with an EM algorithm for the frailty model, and the package provides descriptive, diagnostic and predictive tools for the fitted models. This vignette works through a complete analysis of a real cohort.
Package
gcmrec 2.0.0
1 Introduction
Survival analysis usually asks when an event happens. But many events happen more than once to the same subject: a patient is readmitted to hospital repeatedly, a machine breaks down again and again, a cancer relapses. Treating each occurrence as an independent observation ignores three things that make recurrent data different:
- The subject is not new after an event. A patient discharged after a readmission is not the same as a patient just out of surgery. How much the intervention “rejuvenates” the subject is itself part of the model.
- Occurrences accumulate. The tenth readmission may carry a different risk than the first, beyond anything the covariates explain.
- Subjects differ in ways we cannot measure. Some patients are frail and accumulate events; their inter-event times are not independent.
The model of Peña and Hollander (2004) puts all three in a single framework, and gcmrec fits
it. This vignette shows how, and — more importantly — how to choose its options
and read what it returns.
1.1 The model
You do not need the details to use the package, but the intuition helps interpret the output. The risk of the next event for subject \(i\) at time \(s\) is
\[\lambda_i(s) = \lambda_0\!\left[\mathcal{E}_i(s)\right]\; \rho\!\left(N_i(s^-); \alpha\right)\; \exp(X_i'\beta)\; Z_i ,\]
built from four pieces. The effective age \(\mathcal{E}_i(s)\) is how old the subject behaves, which resets (fully or partially) at each event — this is what encodes the effect of the intervention. The baseline \(\lambda_0\) is left unspecified (semiparametric, as in Cox regression) and is evaluated at the effective age rather than at calendar time. The function \(\rho(k;\alpha) = \alpha^k\) carries the effect of the \(k\) accumulated occurrences: \(\alpha > 1\) means risk grows with each event, \(\alpha < 1\) that it falls, \(\alpha = 1\) that occurrences leave the risk unchanged. Finally \(Z_i\) is an optional gamma frailty, a subject-specific multiplier that induces dependence between the inter-event times of the same subject.
2 Installation
Install the released version from CRAN:
install.packages("gcmrec")The package needs a C++ compiler at install time, which the standard R toolchain provides on all platforms, and depends on Rcpp/RcppArmadillo for the numerical core, survival for the model frame machinery and ggplot2 for the graphics. All of them are installed automatically.
Load it with:
library(gcmrec)3 Quick start
A complete analysis takes three lines: build the response with Survr(), fit
with gcmrec(), and read the result. Everything that follows in this vignette
expands on these steps.
data(readmission)
fit <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
data = readmission, s = 3000)
summary(fit) # hazard ratios with confidence intervals
plot(fit) # baseline survivor function4 Input data
gcmrec works on data in long format: one row per inter-event time, with a
subject identifier repeated across its recurrences. Three variables are
essential — id, time (the gap since the previous event, not calendar
time) and event (1 for an event, 0 for the censoring time) — plus any
covariates.
We use the readmission cohort (González, Fernández, et al. 2005):
861 records from 403 patients operated on for colorectal
cancer, followed for rehospitalisations.
data(readmission)
head(readmission, 4)
#> id enum t.start t.stop time event chemo sex dukes charlson
#> 1 1 1 0 24 24 1 2 2 3 3
#> 2 1 2 24 457 433 1 2 2 3 0
#> 3 1 3 457 1037 580 0 2 2 3 0
#> 4 2 1 0 489 489 1 1 1 2 0
length(unique(readmission$id)) # patients
#> [1] 403
table(table(readmission$id) - 1) # events per patient
#>
#> 0 1 2 3 4 5 6 8 9 10 11 16 22
#> 199 105 45 21 15 8 4 1 1 1 1 1 1Half of the patients are never readmitted, while a few accumulate many events — one reaches 22. That long tail is the kind of structure a recurrent model exploits and a single-event analysis would discard.
The response is built with Survr(), which plays the role that Surv() plays
in ordinary survival analysis:
head(Survr(readmission$id, readmission$time, readmission$event), 4)
#> id time event
#> [1,] 1 24 1
#> [2,] 1 433 1
#> [3,] 1 580 0
#> [4,] 2 489 14.1 Adding missing censoring times
Survr() requires each subject to end with a censored record (event = 0),
because the model needs to know how long the subject was observed after its
last event. Data sets where follow-up ends exactly at the last event fail with
Data doesn't match. addCenTime() repairs them by appending a row with time 0
and event = 0:
dat <- data.frame(id = c(1, 1, 2, 2),
time = c(5, 3, 7, 4),
event = c(1, 0, 1, 1)) # subject 2 ends on an event
addCenTime(dat)
#> id time event
#> 1 1 5 1
#> 2 1 3 0
#> 3 2 7 1
#> 4 2 4 1
#> 41 2 0 04.2 Other input formats
Historical gcmrec data sets are stored as a nested list (elements n and
subject) rather than a data frame. You can pass them to gcmrec() directly —
they are converted internally by as_gcmrec_data() — or convert them yourself
with List.to.Dataframe(). The hydraulic data set, times to failure of six
mining machines (Kumar and Klefsjö 1992), is in that format:
data(hydraulic)
head(List.to.Dataframe(hydraulic), 3)
#> id time event covar.1 covar.2
#> 1 1 327 1 0 0
#> 2 1 125 1 0 0
#> 3 1 7 1 0 0as_gcmrec_data() is a standard S3 generic, so support for further input
classes is a matter of adding a method. Tibbles and data.tables already work,
since they inherit from data.frame.
5 Exploring the data
Before fitting anything, look at the events. Two complementary views answer two different questions.
5.1 Event chart
graph.caltimes() draws one row per subject, a point at each recurrence and a
cross at the end of follow-up. By default subjects are sorted by length of
follow-up, which turns the chart into something readable:
graph.caltimes(readmission[readmission$id %in% 1:40, ])
Figure 1: Rehospitalisations of the first 40 patients
Each row is a patient, each dot a readmission, the cross the end of follow-up.
The picture already tells you what to expect: many patients with a single cross
and no event, a few with dense clusters of readmissions. Passing a variable to
var colours the subjects, and sortevents = "events" orders them by how many
recurrences they had — a quick way to see whether a covariate separates high-
from low-recurrence patients.
5.2 Mean cumulative function
The individual chart does not scale beyond a few dozen subjects. The mean cumulative function (MCF) summarises the same information for the whole cohort: it is the expected number of events accumulated by one subject up to each time, estimated without assuming any model.
m <- mcf(readmission, group = readmission$dukes)
m
#> Mean cumulative function
#> 1 136 event times, 2.06 events per subject by t=2175
#> 2 165 event times, 2.70 events per subject by t=2119
#> 3 122 event times, 5.90 events per subject by t=1572
plot(m)
Figure 2: Mean cumulative function by Dukes’ stage: the expected number of readmissions per patient
This single figure is the reason to fit a recurrent event model at all. A patient with a stage D tumour accumulates about 5.9 readmissions over the follow-up, against roughly 2.1 for stage A-B — and the separation appears early and grows steadily. The bands are pointwise confidence intervals; they widen at the right, where few patients remain under follow-up.
6 Fitting the model
gcmrec() is the main entry point. You give it a formula with a Survr()
response, the data, and a calendar time s up to which the analysis runs.
fit <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
data = readmission, s = 3000)
fit
#> Call:
#> gcmrec(formula = Survr(id, time, event) ~ as.factor(dukes) +
#> sex, data = readmission, s = 3000)
#>
#>
#> coef exp(coef) se(coef) z p
#> as.factor(dukes)2 0.398 1.488 0.112 3.54 4.0e-04
#> as.factor(dukes)3 0.991 2.695 0.133 7.44 1.0e-13
#> sex -0.361 0.697 0.102 -3.55 3.9e-04
#>
#> General class model parameter estimates
#> rho function: Alpha to k
#> alpha (s.e.): 1.12 (0.0142)
#>
#> log-likelihood=-2719.7
#> n= 403
#> n times= 861
#> number of iterations: 7 Newton-Raphson6.1 Arguments
A handful of arguments let you match the model to your problem:
| Argument | Default | What it controls |
|---|---|---|
s |
— | calendar time up to which subjects are followed. Events after s do not contribute. |
rhoFunc |
"alpha to k" |
the effect of accumulated occurrences: \(\rho(k;\alpha)=\alpha^k\), or "Identity" for \(\rho \equiv 1\) (occurrences carry no extra effect, so no \(\alpha\) is estimated). |
typeEffage |
"perfect" |
the effective age model: "perfect" repair resets the subject to as-good-as-new at every event; "minimal" leaves it as-bad-as-old. |
effageData |
NULL |
supply your own effective age per subject instead of generating it (see below). |
cancer |
NULL |
effective age driven by treatment response ("CR"/"PR"/"SD"), for the cancer model of González et al. (2005). |
Frailty |
FALSE |
add a gamma frailty per subject, fitted by EM. |
se |
"Information matrix" |
how standard errors are obtained; "Jacknife" for leave-one-out estimates. |
maxXi |
"Newton-Raphson" |
maximiser for the frailty parameter \(\xi\); "Brent" is a derivative-free alternative (Brent 1973). |
tol, maxit |
1e-6, 100 |
convergence tolerance and iteration cap. |
6.2 Interpreting the output
The print above has three parts. The coefficient table is read exactly as
in a Cox model: exp(coef) is a hazard ratio, so Dukes’ stage D patients have
2.69 times the readmission risk of stage A-B
patients, and women (sex = 2) have a lower risk than men. The alpha
line is what a Cox model cannot give you: here
\(\hat\alpha \approx 1.12\) with a small standard
error, meaning each readmission raises the risk of the next by about
12% — genuine event accumulation, on
top of the covariates. Last come the fit summaries: log-likelihood, number
of subjects, total records and iterations.
summary() returns the hazard ratios with confidence intervals as an object
with its own print method:
summary(fit)
#> hr 95% C.I.
#> as.factor(dukes)2 1.49 ( 1.19 - 1.85 )
#> as.factor(dukes)3 2.69 ( 2.08 - 3.50 )
#> sex 0.70 ( 0.57 - 0.85 )and plotForest() turns that table into the figure that usually goes into a
report, on a logarithmic axis with a reference line at 1:
plotForest(fit, labels = c("as.factor(dukes)2" = "Dukes C vs A-B",
"as.factor(dukes)3" = "Dukes D vs A-B",
"sex" = "Female vs male"))
Figure 3: Hazard ratios with 95% confidence intervals
Intervals crossing the dashed line are compatible with no effect.
The standard extractors work as with any other model in R, so the fit plugs into the usual tooling:
coef(fit)
#> alpha as.factor(dukes)2 as.factor(dukes)3 sex
#> 1.1204108 0.3975220 0.9912442 -0.3607120
sqrt(diag(vcov(fit))) # standard errors
#> alpha as.factor(dukes)2 as.factor(dukes)3 sex
#> 0.01424625 0.11223462 0.13322487 0.10160505
logLik(fit)
#> 'log Lik.' -2719.705 (df=4)
AIC(fit)
#> [1] 5447.4096.3 Testing model terms
The alpha line said that risk grows with each occurrence, but is that
difference real? anova() answers with a likelihood ratio test against the
same model with \(\rho \equiv 1\) (that is, \(\alpha = 1\): occurrences carry no
effect), which it refits internally:
anova(fit)
#> Likelihood ratio test for gcmrec models
#>
#> npar logLik Df Chisq Pr(>Chisq)
#> rho = Identity (alpha = 1) 3 -2748
#> rho = alpha^k 4 -2720 1 56.23 6.45e-14 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The test rejects \(\alpha = 1\) decisively, so the accumulation of events is not an artefact: a model that ignores it — an ordinary Cox model on the gap times, for instance — would be misspecified for these data.
Given two fits, the same function tests whether the extra covariates of the larger one are worth keeping:
fit.small <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
data = readmission, s = 3000)
anova(fit.small, fit)
#> Likelihood ratio test for gcmrec models
#>
#> npar logLik Df Chisq Pr(>Chisq)
#> fit.small 3 -2726
#> fit 4 -2720 1 13.05 0.000303 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A likelihood ratio test needs nested models fitted to the same data, so
anova() checks that and refuses the comparison otherwise:
| Comparison | Allowed? |
|---|---|
| One model against itself with \(\alpha = 1\) | yes, anova(fit) |
| Nested covariates, everything else equal | yes |
rhoFunc = "Identity" against "alpha to k" |
yes (this is the \(\alpha\) test) |
Different s, typeEffage, effageData or cancer |
no — not nested |
| Different data (different subjects or records) | no |
| Covariates that are not a subset of one another | no |
| One fit with frailties and one without | no — see section 6 |
| Two fits with frailties | approximate, with a warning |
6.4 Baseline functions
plot() draws the estimated baseline survivor function against effective
age (not calendar time), with a pointwise confidence band; type.plot = "hazard" draws the cumulative hazard instead.
plot(fit)
Figure 4: Baseline survivor function on the effective age scale, with its 95% confidence band
Every plotting function returns a ggplot object, so you can restyle it with
the usual syntax without the package getting in the way:
plot(fit, type.plot = "hazard", level = 0.99) +
ggplot2::labs(title = "Baseline cumulative hazard",
subtitle = "Colorectal cancer readmissions")
Figure 5: The same curve, restyled
6.5 Predictions for covariate profiles
The baseline describes a subject with all covariates at zero, which is rarely
an interesting patient. plotPredict() draws the curve implied by concrete
covariate profiles, and predict() returns the same numbers:
profiles <- data.frame(dukes = c(1, 2, 3), sex = 1)
predict(fit, profiles, type = "risk") # relative risk of each profile
#> 1 2 3
#> 0.6971798 1.0374959 1.8786102
plotPredict(fit, profiles,
labels = c("Dukes A-B", "Dukes C", "Dukes D"))
Figure 6: Predicted survivor function of the next readmission, by Dukes’ stage (men)
Read these as the probability of not yet having the next readmission as a function of effective age: at any point, a stage D patient is markedly more likely to have recurred than a stage A-B one.
7 Effective age models
The effective age is the heart of the model: it says what an event does to the subject. Two extremes come built in.
Perfect repair (the default) resets the effective age to zero at every event: the patient leaves the hospital as good as new. Minimal repair leaves the effective age untouched: the patient is exactly as old as before, and the event changes nothing.
mod.per <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
data = readmission, s = 3000, typeEffage = "perfect")
mod.min <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
data = readmission, s = 3000, typeEffage = "minimal")
rbind(perfect = coef(mod.per), minimal = coef(mod.min))
#> alpha as.factor(dukes)2 as.factor(dukes)3 sex
#> perfect 1.120411 0.3975220 0.9912442 -0.3607120
#> minimal 1.204088 0.4396272 1.2408367 -0.3771564The covariate effects are similar under both, but \(\alpha\) and the baseline are
not — because the two assumptions place the events on different time scales.
plotBaseline() takes a named list of fits and draws them together:
plotBaseline(list(perfect = mod.per, minimal = mod.min))
Figure 7: Baseline survivor function under perfect and minimal repair
7.1 Effective age from treatment response
Between those extremes, the intervention may repair the subject partially, and
by an amount the data tell you. The cancer argument implements the model of
González, Peña, et al. (2005) for cancer relapses, where the effective age after each
treatment is set by the response achieved: complete remission ("CR") resets it
to zero, partial remission ("PR") advances it by half the elapsed time, and
stable disease ("SD") — a null response — by the whole elapsed time.
The lymphoma data set carries exactly that variable:
data(lymphoma)
table(lymphoma$effage)
#>
#> CR PR SD
#> 101 9 2
mod.can <- gcmrec(Survr(id, time, event) ~ as.factor(distrib),
data = lymphoma, s = 1000, cancer = lymphoma$effage)
coef(mod.can)
#> alpha as.factor(distrib)1 as.factor(distrib)2 as.factor(distrib)3
#> 0.9684011 0.8395349 1.1647179 1.0022586If your effective age follows none of these schemes, compute it yourself and
pass it through effageData, a list with one entry per subject (intercepts,
slopes, lastperrep, perrepind, effagebegin, effage). GeneratedData
is an example in that format:
data(GeneratedData)
dat.gen <- List.to.Dataframe(GeneratedData)
mod.eff <- gcmrec(Survr(id, time, event) ~ covar.1 + covar.2,
data = dat.gen, effageData = GeneratedData, s = 100)
coef(mod.eff)
#> alpha covar.1 covar.2
#> 0.9946750 0.4980746 1.01806738 Frailty models
Two patients with identical covariates may still accumulate events at different
rates. A gamma frailty \(Z_i\) captures that unobserved heterogeneity, and
induces dependence among the inter-event times of the same subject. Setting
Frailty = TRUE fits the model by an EM algorithm:
mod.fra <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
data = readmission, s = 3000, Frailty = TRUE)
coef(mod.fra)
#> alpha as.factor(dukes)2 as.factor(dukes)3 sex
#> 1.0799902 0.4053905 1.1218867 -0.4430787
mod.fra$Xi # frailty parameter
#> [1] 2.511783The parameter \(\xi\) is read as an inverse measure of heterogeneity: the frailty has variance \(1/\xi\), so small \(\xi\) means very different subjects and \(\xi \to \infty\) recovers the model without frailties. Here \(\hat\xi \approx 2.51\) points to substantial heterogeneity between patients, which is why \(\hat\alpha\) drops compared with the model without frailties: part of what looked like event accumulation was really the frail patients contributing most of the events.
The estimated frailties themselves are returned, one per subject, and are informative in their own right:
summary(mod.fra$frailties)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.5660 0.7759 0.8842 1.0006 1.1155 4.7716
ggplot2::ggplot(data.frame(z = mod.fra$frailties),
ggplot2::aes(x = z)) +
ggplot2::geom_histogram(bins = 30, fill = "#0072B2", alpha = 0.85) +
ggplot2::labs(x = "Estimated frailty", y = "Patients") +
theme_gcmrec()
Figure 8: Estimated frailties
Patients to the right accumulate events faster than their covariates predict.
Note that anova() does not apply here: the log-likelihood of a frailty fit is
conditional on the estimated frailties, so it is not on the same scale as that
of a model without them. Judge the heterogeneity from \(\xi\) itself — with a
standard error if you fit with se = "Jacknife" — and from the spread of the
frailties, which here run from
0.57 to 4.77.
Frailty models are estimated by EM and are therefore slower than the plain fit, though the C++ core keeps a cohort of this size to a couple of seconds.
9 Standard errors
By default standard errors come from the inverse of the partial likelihood
information matrix. The alternative, se = "Jacknife", refits the model
leaving out each subject in turn: more robust, no distributional assumption,
but n refits instead of one. It is the option to use when the information
matrix is suspect — small samples, near-singular information — and the only way
to get standard errors for the frailty model.
sub <- readmission[readmission$id %in% unique(readmission$id)[1:60], ]
mod.info <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
data = sub, s = 3000)
mod.jack <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
data = sub, s = 3000, se = "Jacknife")
rbind(information = sqrt(diag(vcov(mod.info))),
jackknife = sqrt(diag(vcov(mod.jack))))
#> alpha as.factor(dukes)2 as.factor(dukes)3
#> information 0.1348344 0.2912920 0.3612783
#> jackknife 0.2208214 0.3885416 0.3331996The two agree closely here, which is reassuring; a large discrepancy would be a warning that the information matrix is not to be trusted. Since the leave-one-out fits are independent of each other, the jackknife is computed in parallel when the package is built with OpenMP support.
10 Migrating from version 1.x
If you used gcmrec 1.0-5, your analysis scripts still run: the function names,
arguments and the structure of the fitted object are unchanged. What changed is
underneath, plus a few additions:
| Task | Version 1.0-5 | Now |
|---|---|---|
| Numerical core | Fortran 77 | C++ (Rcpp/RcppArmadillo), 6-45x faster |
| Legacy list data | List.to.Dataframe() first |
accepted directly by gcmrec(); as_gcmrec_data() is the extension point |
summary() |
printed to the console | returns a summary.gcmrec object with a print method |
| Coefficients, covariance, log-likelihood | fit$coef, fit$var, fit$loglik |
also coef(), vcov(), logLik() (so AIC() works) |
| Plots | base graphics | ggplot2 objects you can restyle; shared look via theme_gcmrec() |
| Comparing two fits | plot(a); lines(b) |
plotBaseline(list(a = a, b = b)) (lines() is deprecated) |
| Maximum recurrences per subject | 200 (hard limit) | unlimited |
New in this version, with no equivalent before: mcf() (mean cumulative
function), plotForest(), predict() and plotPredict() for covariate
profiles, and anova() for likelihood ratio tests.
Two numerical results changed on purpose, both bug fixes. The baseline functions
under rhoFunc = "Identity" were computed from an incomplete coefficient vector
and are now correct; and the frailty jackknife used misaligned offsets in the
leave-one-out refits. If you need to reproduce the old numbers exactly, version
1.0-5 is available from the
CRAN archive and
tagged as v1.0-5 in the package repository.
Session info
#> 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] gcmrec_2.0.0 BiocStyle_2.36.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-4 gtable_0.3.6 jsonlite_2.0.0
#> [4] dplyr_1.2.1 compiler_4.5.3 BiocManager_1.30.27
#> [7] tinytex_0.57 tidyselect_1.2.1 Rcpp_1.1.2
#> [10] dichromat_2.0-0.1 jquerylib_0.1.4 splines_4.5.3
#> [13] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0
#> [16] lattice_0.22-9 ggplot2_4.0.1 R6_2.6.1
#> [19] labeling_0.4.3 generics_0.1.4 knitr_1.50
#> [22] tibble_3.3.1 bookdown_0.43 bslib_0.9.0
#> [25] pillar_1.11.1 RColorBrewer_1.1-3 rlang_1.3.0
#> [28] cachem_1.1.0 xfun_0.52 sass_0.4.10
#> [31] S7_0.2.1 cli_3.6.6 withr_3.0.3
#> [34] magrittr_2.0.5 digest_0.6.39 grid_4.5.3
#> [37] rstudioapi_0.17.1 lifecycle_1.0.5 vctrs_0.7.3
#> [40] evaluate_1.0.4 glue_1.8.1 farver_2.1.2
#> [43] survival_3.8-6 rmarkdown_2.29 tools_4.5.3
#> [46] pkgconfig_2.0.3 htmltools_0.5.8.1