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.

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     0

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       0

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

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, ])

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.

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)

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.

Arguments

A handful of arguments let you match the model to your problem:

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

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

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

The 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 ' ' 1

A likelihood ratio test needs nested models fitted to the same data, so anova() checks that and refuses the comparison otherwise:

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)

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

The same curve, restyled.

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

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.

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

If 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.0180673