Sage Journals: Discover world-class research

Abstract

Recurrent events analysis plays an important role in many applications, including the study of chronic diseases or recurrence of infections. Historically, many models for recurrent events have been variants of the Cox model. In this article we introduce and describe the application of the piece-wise exponential Additive Mixed Model (PAMM) for recurrent events analysis and illustrate how PAMMs can be used to flexibly model the dependencies in recurrent events data. Simulations confirm that PAMMs provide unbiased estimates as well as equivalence to the Cox model when proportional hazards are assumed. Applications to recurrence of staphylococcus aureus and malaria in children illustrate the estimation of seasonality, bivariate non-linear effects, multiple timescales and relaxation of the proportional hazards assumption via time-varying effects. The R package pammtools is extended to facilitate estimation and visualization of PAMMs for recurrent events data.

Keywords

Flexible survival models non-linear survival analysis piece-wise exponential models recurrent events time-varying effects

1 Introduction

A recurrent events setting in survival analysis is defined by repeated observations of an event of interest over the course of the observation period. In medical, clinical and biological research, such data play an important role, for example, in the context of chronic illnesses and infectious diseases. Some concrete examples include recurrence of hospitalization for cardiovascular events (Gibson et al., 2017; Varma et al., 2021), recurrent diseases and infections like COVID-19 (Dos Santos et al., 2021), pneumonia (Ramjith et al., 2021; Barakat et al., 2021), pneumococcus (Hernstadt et al., 2020), staphylococcus aureus (cf. Section 5.1) (Akinboyo et al., 2018; Abdulgader et al., 2019), asthma (Chan et al., 2021; Feng et al., 2021) and malaria (cf. Section 5.2) (Lawpoolsri et al., 2019; Ghosh et al., 2020; Kakuru et al., 2020).

Typical for this kind of data are the correlations, or more generally dependencies, between event recurrences. Such dependencies may be induced by an individual-specific unmeasured frailty, common across all events of the individual, or by dependence of the occurrence of an event on an individual’s history. In the latter case the hazard of an event could depend on the timing (e.g., time since last event) or the frequency (number of previous events). Ignoring either type of dependency, if present, leads to biased estimates or overoptimistic standard errors (Box-Steffensmeier and De Boef, 2006; Jahn-Eimermacher, 2008). To choose an appropriate model for the analysis of recurrent event data, assumptions on this dependency structure must be made as well as a choice on the timescale for the analysis of events: gap time (also referred to as waiting time, clock-reset or renewal approach) or calendar time (clock-forward). See Section 2.1 for more discussion about the timescale.

Several recurrent events models for calendar and gap times have been suggested in the literature. These models attempt to account for the effects of within-subject correlation either by adjusting the variances of the parameter estimators (the variance-adjustment models (Andersen and Gill, 1982; Prentice et al., 1981)) or by including a subject-specific random effect in the model—the shared frailty model (SFM). Variance-adjustment models assume independence among events when estimating the effect of the covariates, while adjusting the standard errors of the parameter estimates for the correlation between recurrent events by calculating robust variance through sandwich estimation (Huber, 1967; Freedman, 2006). SFMs (Hougaard, 2012) account for correlated events within an individual (and thus heterogeneity in the sample) by adding individual specific random effects, which act in a multiplicative way on the hazard function. In this work, we focus on the frailty (or random effects) model. In a standard SFM it is assumed that the hazards are proportional over time (conditional on the random effect) and the baseline hazard function is unspecified (in the Cox SFM), but a parametric form may be chosen as well. For interpretation of the estimated effects, an estimate of the baseline hazard function is important. Putter et al. (2005) mentioned that not knowing the baseline implies that the risk of observing an event at particular points in time remains unknown. The importance of knowing the absolute risk (or baseline hazard) when interpreting relative risks is further highlighted by Noordzij et al. (2017). They suggest that readers tend to overestimate or over-interpret relative risks when information about the absolute risks is missing. This is even strengthened if the hazard ratio is expected to be time-dependent. An understanding of the baseline hazard over time is important for clearer interpretation of time-varying effects (Ruhe, 2018).

Many publications on recurrent events analysis focus on variants of the Cox proportional hazards (CPH) model (Cox, 1972). This includes the SFM, mentioned earlier, and the variance-adjustment models, specifically the Andersen-Gill (AG) model, the Prentice-Williams-Peterson (PWP) models, and the Wei, Lin and Weissfeld (WLW) model (Andersen and Gill, 1982; Prentice et al., 1981; Wei et al., 1989). The AG model estimates a common baseline hazard for all events and covariate effects are usually estimated across events as well. In the PWP and WLW models the baseline hazards are stratified for the different events, and covariate estimates are either estimated across events or event-specific. The PWP models are constructed in either a calendar time approach or a gap time approach (see Section 2.1), while the WLW model is constructed using marginal times (where the time at risk for each event is counted from the study start). In the AG and PWP models, individuals only remain at risk for the event following their last observed event, while for the WLW model, individuals remain at risk for the maximum number of observed events. Studies (Kelly and Lim, 2000; Therneau and Hamilton, 1997) have shown a carry-over effect when using the WLW model and recommend that it can be used for modelling multiple event types rather than recurrent events. There are several tutorials available with applications and comparisons of these methods (Kelly and Lim, 2000; Therneau and Hamilton, 1997; Amorim and Cai, 2015; Ullah et al., 2014), and a recent tutorial on SFMs is provided by Balan and Putter (2020).

Parametric SFMs can be a useful alternative to the semi-parametric SFM. In the parametric SFM the baseline hazards can be modeled through a choice of many parametric distributions each with its own characteristics (Munda et al., 2012). Misspecification of the baseline hazards can lead to bias in estimated effects and possibly misleading results for the shape of the hazard function (Hutton and Monaghan, 2002; Kwong and Hutton, 2003; Li et al., 1996). In a recurrent events setting where the shape of the baseline hazard differs, for example between first and recurrent events, a mixture of parametric models must be used to appropriately capture the baseline hazards (Duchateau et al., 2003). This can be a complicated task.

More flexible models like the generalized survival models (GSMs) (Liu et al., 2017) are becoming increasingly more popular in modelling the baseline hazards. These models use a smoothing splines approach and allow more complexities such as time-dependent effects, nonlinear effects, and their interactions. The GSM has recently been extended to model correlated survival data: SFMs with time-constant covariates for gap time formulations or an AG variance-adjustment type model for calendar time formulations. However, GSMs still needs to be extended to SFMs in a calendar time framework.

Here we introduce the Piece-wise exponential Additive Mixed Model (PAMM) (Bender et al., 2018; Argyropoulos and Unruh, 2015) for modelling possibly time-varying effects, nonlinear effects and more complex interactions in a recurrent events setting where it is possible to consider multiple time scales, stratified baseline hazards and time-varying covariates. As the Cox model, PAMMs model the hazard rate without parametric assumptions about the distribution of event times and multiplicative covariate effects. Thus, if model specification is equivalent, both models lead to similar results and interpretation.

While conceptually similar to the Cox model, PAMMs offer an alternative method for survival analysis that has recently gained popularity. This is mainly due to the fact that PAMMs reformulate many different survival tasks to regression tasks via a data augmentation step and availability of software that facilitates practical application (Bender and Scheipl, 2018). While this allows the model to be estimated by any algorithm that can optimize the Poisson likelihood (Bender et al., 2021), the Generalized Additive Mixed Model (GAMM) framework has proven to be particularly useful, due to its flexibility, availability of robust and established estimation techniques (e.g., Hastie and Tibshirani (1986); Wood (2017)) and interpretability of covariate effects. This allows great versatility in modelling choices, as discussed in more detail in Section 6. In particular, PAMMs facilitate seamless estimation of multivariate, non-linear, time-varying effects via penalized splines with different basis functions, cyclic splines and splines with shape constraints (Pya and Wood, 2015). Recent examples include spline based stratified baseline hazard estimation in plant biology (Panel et al., 2020), cumulative effects in critical care (Bender et al., 2019), treatment effects with non-proportional hazards in patients with coronary stent restenosis (Giacoppo et al., 2020), estimation of bivariate effect surfaces in the context of randomized controlled trials (Argyropoulos and Unruh, 2015) and progression-free survival in patients with gastroenteropancreatic neuro-endocrine tumors (Carmona-Bayonas et al., 2019), and time-varying hazard ratios in post-transplant outcome assessment (Wey et al., 2020). An application to the recurrence of childhood pneumonia was recently presented in Ramjith et al. (2021). In general, however, detailed instructions on the application of PAMMs to survival tasks, especially tasks that go beyond the single-event setting with right-censoring as well as software tools that facilitate such analyses, particularly recurrent events analyses, are scarce.

In this article we give an introduction to PAMMs for recurrent events analysis. This includes data transformation specific to recurrent events and discussion on the choice of timescale. We show that the model performs well and generates unique insights by means of simulation studies and real data examples. The article is accompanied by the R package pammtools that was extended to facilitate application of PAMMs to recurrent events data. The pammtools package takes care of the required data transformation and greatly facilitates model estimation, visualization and interpretation of the estimated effects.

In the remainder of this work, we summarize hazard based model specification and estimation for recurrent events analysis in Section 2. In Section 3, we show how such models can be estimated using PAMMs, including a step by step guide on the necessary data augmentation. A comparison of the SFM and the PAMM when proportional hazards are assumed (conditional on the frailty) is made through a simulation study in Section 4. Section 5 concludes with examples of PAMMs applied to different recurrent events data sets and illustrates effect estimation of different complexities, that is, stratified baseline hazards, non-linear, cyclic effects of seasonality and multiple timescales. We conclude with a discussion and an outlook in Section 6.

2 Models for recurrent events

2.1 Timescale

The choice between the two timescales (gap time or calendar time) is usually driven by the research setting in which the data are collected, and the research question at hand, as the choice of the timescale affects the interpretation of the results (see Duchateau et al., 2003; for a thorough discussion). On both timescales, the study start time needs to be clearly defined and meaningful, for example, time from birth, time from randomization, time from diagnosis. Calendar time is appropriate when the interest is in the full course of the recurrent event process. On this timescale an individual’s time starts when entering the study and stops when leaving the study. For the first event, the individual’s contribution to time at risk is counted from entering the study. However, the individual’s contribution to time at risk for the second event is only counted from the first event, the start time of the risk interval for the second event. A calendar time data set can be viewed as a succession of (internally) left-truncated data sets for each of the events, as a subject’s time at risk is left-truncated at the event time of the previous event.

In the gap time approach one is interested in the time between events. Essentially, that means that after every event ‘the clock’ is reset; after each event the time since the previous event is used to define time intervals that an individual is at risk for an event. Thus, in contrast to the calendar time approach, a subject is at risk for all of his events at time t = 0, however, time to event is redefined as time from last event to next event. Historically, the gap time approach facilitates analysis, because the data could be evaluated with standard techniques for survival analysis without the need to take into account left-truncation. Assuming a full renewal process, all observations are considered independent. Dependence induced by number of previous events can be modeled by stratifying the hazards by event number. Dependence introduced by within-subject correlation can be modeled by introducing a frailty term.

In the PAMM framework, the choice of timescale is less important, as we can always add additional time-dependent covariates to describe the past. In calendar time, for example, we could include a variable for the gap time of the previous event or the number of previous events. In the gap time approach, we can include dependence on the total time since entering the study, and so on. More generally, PAMMs support dependence on multiple timescales (Iacobelli and Carstensen, 2013), as any dependence on time is modeled through covariate effects of different representations of time.

2.2 Model speciftcation

In the following we consider models for the (conditional) hazard function for the k − th event for subject i given by

λ_{i, k} (t ∣ z_{i}) : = λ (t ∣ x_{i} (t), z_{i}, k) = λ_{0, k} (t) \exp (g (x_{i} (t), t, k) + z_{i}), i = 1, \dots, n,

(2.1)

where, t is the time of interest, $x_{i} (t) = {(x_{i, 1} (t), x_{i, 2} (t), \dots x_{i, p} (t))}^{⊤}$ a vector of potentially time-varying covariates, z_i a subject specific unobserved random effect and g(x _i (t), t, k) is a general function of potentially non-linear, time-varying and event-specific covariate effects, not to be confused with a transformation function of the covariates. Further, $λ_{0, k} (t)$ are event-specific baseline hazards. In the CPH model and its extensions, baseline hazards are not specified and estimated non-parametrically. Standard models discussed in the literature follow from 2.1 by relaxing some of the dependencies. For example, with linear effects and under the proportional hazards assumption (conditional on frailty) we get the stratified SFM (Box-Steffensmeier and De Boef, 2006):

λ (t ∣ x_{i}, z_{i}, k) = λ_{0, k} (t) \exp (x_{i}^{⊤} β + z_{i}), i = 1, \dots, n .

Assuming a common baseline hazard for all events, the model could be further simplified by replacing λ_0,k with λ₀, resulting in the so-called unstratified (or unconditional on episode order) SFM (Hougaard, 2012).

The conditional likelihood for model 2.1 in calendar time (given the individual frailty terms) is proportional to

L (t ∣ z) = \prod_{i = 1}^{n} L_{i} (t_{i, 1}, \dots, t_{i, K_{i}} ∣ z_{i}) \propto \prod_{i = 1}^{n} \prod_{k = 1}^{K_{i}} λ_{i, k} {(t_{i, k} ∣ z_{i})}^{δ_{i, k}} \exp (- \int_{t_{i, k - 1}}^{t_{i, k}} λ_{i, k} (s ∣ z_{i}) d s),

(2.2)

where K_i is the number of events for which subject i was at risk, t_i,k the event time for event k for k < K_i, t_i,Ki is the time subject i was censored, and t_i,₀ the entering time for subject i which is usually equal to 0. Furthermore, δ_i,k is the event specific status indicator, so δ_i,k = 1 for k < K_i and $δ_{i, K_{i}} = 0$ (see also Duchateau et al., 2003). The likelihood in (2.2) is often given in the simplified form

L (t ∣ z) \propto \exp (- \sum_{i = 1}^{n} \int_{0}^{t_{i, K_{i}}} λ_{i, k} (s ∣ z_{i}) d s) \prod_{i = 1}^{n} \prod_{k = 1}^{K_{i} - 1} λ_{i, k} (t_{i, k} ∣ z_{i}),

(2.3)

The likelihood functions given so far are conditional the frailty variable. The marginal (unconditional) likelihood function is proportional to

L (t) \propto \prod_{i = 1}^{n} (\int_{z_{i}} \prod_{k = 1}^{K_{i}} λ_{i, k} {(t_{i, k} ∣ z_{i})}^{δ_{i, k}} \exp (- \int_{t_{i, k - 1}}^{t_{i, k}} λ_{i, k} (s ∣ z_{i}) d s) f_{Z} (z_{i}) d z_{i}),

(2.4)

where f_Z is the density function of the frailty variable.

The likelihood for the gap time model (conditional on frailty) is equivalent to 2.2 with t_i,k replaced by $d_{i, k} = t_{i, k} - t_{i, k - 1}$ and the lower integration limit set to 0:

L (t ∣ z) = \prod_{i = 1}^{n} L_{i} (t_{i, 1}, \dots, t_{i, K_{i}} ∣ z_{i}) \propto \prod_{i = 1}^{n} \prod_{k = 1}^{K_{i}} λ_{i, k} {(d_{i, k} ∣ z_{i})}^{δ_{i, k}} \exp (- \int_{0}^{d_{i, k}} λ_{i, k} (s ∣ z_{i}) d s) .

(2.5)

3 PAMMs for recurrent events analysis

3.1 Data augmentation and estimation

As described in detail elsewhere (Bender et al., 2018), application of PAMMs to time-to-event data requires a particular data augmentation step. We refer to the resulting data as piece-wise exponential data (PED). In principle, the data augmentation required for recurrent events analysis is equivalent, however, some specifics need to be taken into account, that also depend on the timescale of the analysis.

As for the single-event PED, the follow-up has first to be discretized into J intervals with cutpoints $0 = τ_{0} < τ_{1} < \dots < τ_{J}$ . The j^th interval is defined as $(τ_{j - 1}, τ_{j}]$ and τ_J as the maximum time, which can be set arbitrarily. In practice, the maximum observed event time is a typical choice, since no information about the event process is contained in the data beyond that point. The other cutpoints could be set to the unique observed event times. In the gap time analysis, cutpoints could also be set arbitrarily or at a subset of unique event times. In the calendar time approach, cut-points must be set at event times. Note that despite the discretization of the follow up, the exact time-to-event information is used for estimation, thus the PAMM is a method for continuous time-to-event data. Technical details about the data-transformation are provided below. A more intuitive illustration of the data transformation procedure is given in Section 3.2.

Let ${\{τ_{j}\}}_{j = 1}^{J}$ be the set of interval borders that partition the follow-up as before. We define the event-, interval- and subject-specific status and ‘time at risk’ variables as

δ_{i j k} = \{\begin{array}{l} 1 if t_{i, k} \in (τ_{j - 1}, τ_{j}] and δ_{i, k} = 1 \\ 0 else \end{array} and t_{i j k} = \{\begin{array}{l} t_{i, k} - τ_{j - 1} if t_{i, k} \in (τ_{j - 1}, τ_{j}) \\ τ_{j} - τ_{j - 1} else, \end{array}

where i = 1, ..., n the subject identifier as before, j = 1, ..., J denotes the interval and k = 1, ..., K the event number. Thus, δ_ijk indicates whether subject i, experienced an event of type k in interval j (1=yes, 0=no) and t_ijk is the time for which subject i was at risk for the k^th event in interval j. Note that for a particular subject i we only calculate these variables for time points (intervals) at which they were at risk for the particular event number. Doing so also naturally accounts for left-truncation at the beginning of the follow up as well as left-truncation for recurrent events on the calendar timescale (induced by the recurrent events process).

Consider a specific event k. Assuming piece-wise constant hazards

λ_{i, k} (t ∣ x_{i} (t), z_{i}, k) = λ_{j, k} \exp (g (x_{i} (t), t, k) + z_{i}) = λ_{i j k} \forall t \in (τ_{j - 1}, τ_{j}]

(3.1)

in each interval, the likelihood contribution of subject i (conditional on the frailty variable) in 2.2 can be rewritten as

λ_{i, k} {(t_{i, k} ∣ z_{i})}^{δ_{i, k}} \exp (- \int_{t_{i, k - 1}}^{t_{i, k}} λ_{i, k} (s ∣ z_{i}) d s) = (\prod_{j \in J_{i k}} λ_{i j k}^{δ_{i j k}}) \exp (- \sum_{j \in J_{i k}} λ_{i j k} t_{i j k}),

(3.2)

where $J_{i k}$ is the set of intervals for which i is at risk for event k. Here, left-truncation can be accounted for by defining $J_{i k}$ such that it excludes intervals prior to the left-truncation time. The right-hand side of 3.2 can be recognized to be proportional to likelihood contributions to a Poisson likelihood under the working assumption $δ_{i j k} \overset{i i d}{\sim} P o (μ_{i j k} = λ_{i j k} t_{i j k})$ and conditional on the frailties z_i. Making use of this, it can be seen that the conditional likelihood function $L (t ∣ z)$ is proportional to the conditional likelihood function that is found under the working assumption just mentioned:

L (t ∣ z) \propto L_{P o} (t ∣ z) = \prod_{k = 1}^{K} \prod_{i = 1}^{n} (\prod_{j \in J_{i k}} \frac{μ_{i j k}^{δ_{i j k}}}{δ_{i j k}!}) \exp (- \sum_{j \in J_{i k}} λ_{i j k} t_{i j k}) .

(3.3)

Similar to 2.4, the conditional likelihood in 3.3 can be expressed as a marginal likelihood by specifying a frailty distribution:

L_{P o} (t) = \prod_{i = 1}^{n} \int_{z_{i}} \prod_{k = 1}^{K} (\prod_{j \in J_{i k}} \frac{μ_{i j k}^{δ_{i j k}}}{δ_{i j k}!}) \exp (- \sum_{j \in J_{i k}} λ_{i j k} t_{i j k}) f_{Z} (z_{i}) d z_{i},

(3.4)

where $μ_{i j k} = λ_{i j k} t_{i j k}, λ_{i j k} = λ_{i, k} (t ∣ x_{i} (t), z_{i}, k) \forall t \in (τ_{j - 1}, τ_{j}], t_{i j k}$ is the time for which subject i was at risk for the k^th event in interval j and δ_ijk indicates whether subject i, experienced an event of type k in interval j. Thus, parameters of model 2.1 (i.e., the hazard) can be estimated by fitting a mixed-effects Poisson regression model to outcomes δ_ijk and offset o_ijk = log(t_ijk) with subject-specific random intercepts. In a generalized linear mixed-effects modelling framework, the integral can be approximated using numerical techniques such as adaptive Gauss-Hermite quadrature (Bates et al., 2015). In the GAMM framework, the integral is either numerically approximated or the link between penalized splines and random effects models (Lin and Zhang, 1999; Wand, 2003; Wood, 2004, 2017) is used to treat the frailties as penalized splines where the frailty variance is estimated as part of the penalty.

The reformulation of recurrent events models into Poisson regression models is interesting especially due to flexibility of GAMMs and readily available software. In PAMMs we parameterize the interval- and event-specific hazard rates 3.1 via (non-linear) covariate effects. An exemplary model specification is given below:

λ (t ∣ x_{i} (t), z_{i}, k) = \exp (β_{0, k} + f_{0, k} (t_{j}) + \sum_{p = 1}^{P} f_{p} (x_{i, p} (t_{j}), t_{j}) + z_{i}) \forall t \in (τ_{j - 1}, τ_{j}] .

(3.5)

In 3.5, $β_{0, k} + f_{0, k} (t_{j})$ denotes the event-specific baseline log-hazard rate and t_j is simply a covariate that represents time (e.g. t_j = τ_j) that is, used to make the hazard in interval j depend on time, and f_p(x_i,p(t_j), t_j) are covariate effects, that might be linear or non-linear and time-dependent. Non-linear functions are represented through basis functions and coefficients, for example, $f_{0, k} (s) = \sum_{l = 1}^{L} γ_{l, k} b_{l, k} (s)$ where γ_l,k are unknown regression parameters to be estimated and $b_{l, k} (s)$ are known basis functions, with several types of basis functions available (Wood, 2017). Popular choices for basis functions include B-Splines (Eilers and Marx, 1996), cubic regression splines and others. The concrete choice of basis functions, however, is usually not important. The number of basis functions L (or equivalently maximum degree of freedom) are also not important, as long as it is large enough to allow enough flexibility. Overfitting is avoided by enforcing a penalty on the coefficient estimates, that penalizes large deviations between (neighboring) coefficients, leading to a smooth function. Both of these aspects are illustrated in Supplement. As discussed elsewhere (Bender et al., 2018), penalized splines based estimation of the baseline hazard is particularly useful in this context, as it makes the model robust to the choice of interval cutpoints, enforces similarity of neighboring baseline hazards and reduces the number of parameters to estimate ( $L ≪ J$ ). Depending on the assumptions about the data generating process, 3.5 could be simplified, for example, by dropping the dependence on event number or by removing the frailty term.

From an estimation and inference perspective, the model is a Poisson regression task with additive and random effects. Therefore, the parameters of the model can be estimated by optimizing the penalized Poisson likelihood 3.6 using the GAMM framework (e.g., Wood, 2017).

L_{P o} (t ∣ z) + \sum_{m = 1}^{M} v_{m} γ_{m}^{⊤} S_{m} γ_{m} .

(3.6)

In 3.6, M is the total number of smooth and random effects in the model and ν_m the respective penalty parameter estimated from the data; γ _m is the coefficient vector associated with the m-th smooth or random effect and S _m the corresponding penalty matrix that will depend on the type of basis functions chosen for the m-th model term. Gaussian random effects can be subsumed in Equation (3.6) because of the equivalence of smooths and random effects (see for example, Wood, 2017, Sec. 5.8 and references therein). For example, if the m-th term is a random effect, $γ_{m} = {(γ_{1, m}, \dots, γ_{n, m})}^{⊤}$ represents the random intercepts, where n is the number of subjects in the data set. The corresponding basis functions in this case are 0/1 coded dummy variables for each subject (i.e. n dummy coded columns in the design matrix) and the penalty matrix S _m for these terms is an n × n identity matrix. The tuning parameter ν_m controls the shrinkage of the random effects term and the variance of the random effects is the inverse of ν_m.

Joint optimization of all parameters in 3.6 can be done in different ways. In the application section, we use the methods implemented in package mgcv, which provides different methods for optimization, including generalized cross validation, maximum likelihood and restricted maximum likelihood (see, Wood, 2011 for a comparison).

3.2 An example of recurrent events data in PED format

Exemplary recurrent time-to-event data for two individuals, $i \in {1, 2}$ both with a 3 months follow-up is given in Table 1 and will be used to illustrate the data augmentation step discussed in the previous section. Individual i = 1 has the first event at 0.5 months since the start of the study and does not experience a second event for the remainder of the follow-up. Individual i = 2 experiences the first event at 0.8 months since the start and the second event at 1.2 months since the start (or 0.4 months since the first event) and does not experience the third event in the remainder of the follow-up. As seen in Table 1, the events that each individual is at risk for are captured in different rows. Columns ‘entry’ and ‘exit’ define the time span (calendar time) in which the subjects are at risk for event k.

Table 1

Recurrent event data in ‘standard’ format, where each individual i has one row per event k he was at risk.

i	entry (t_i,k−1)	exit (t_i,k)	gap (d_i,k)	status (δ_i,k)	event (k)
1	0	0.5	0.5	1	1
1	0.5	3.0	2.5	0	2
2	0	0.8	0.8	1	1
2	0.8	1.2	0.4	1	2
2	1.2	3.0	1.8	0	3

3.2.1 Gap time

In Table 1 the unique event-specific gap times, ordered with respect to length, are given by 0.4, 0.5, 0.8 (values from column d_i,k in Table 1 with δ_ik = 1). The maximum gap time is 2.5, however, it usually makes sense to set interval border of the last interval J to the maximum event time, which is 0.8 in this example as there is no information about event times beyond that point. Therefore we will use cutpoints (0,0.4,0.5,0.8) and the intervals in this example will be defined as

{\{(τ_{j - 1}, τ_{j}]\}}_{j = 1}^{J} = {(0, 0.4], (0.4, 0.5], (0.5, 0.8]} .

The transformation of Table 1 to the PEM data format based on these intervals is shown in Table 2.

Table 2

Recurrent event data in PEM data format on the gap timescale.

i	gap (t_i,k)	k	j	(τ_j−1, τ_j]	δ_ijk	t_ijk	o_ijk = log (t_ijk)
1	0.5	1	1	(0,0.4]	0	0.4	log(0.4)
1	0.5	1	2	(0.4,0.5]	1	0.1	log(0.1)
1	2.5	2	1	(0,0.4]	0	0.4	log(0.4)
1	2.5	2	2	(0.4,0.5]	0	0.1	log(0.1)
1	2.5	2	3	(0.5,0.8]	0	0.3	log(0.3)
2	0.8	1	1	(0,0.4]	0	0.4	log(0.4)
2	0.8	1	2	(0.4,0.5]	0	0.1	log(0.1)
2	0.8	1	3	(0.5,0.8]	1	0.3	log(0.3)
2	0.4	2	1	(0,0.4]	1	0.4	log(0.4)
2	1.8	3	1	(0,0.4]	0	0.4	log(0.4)
2	1.8	3	2	(0.4,0.5]	0	0.1	log(0.1)
2	1.8	3	3	(0.5,0.8]	0	0.3	log(0.3)

3.2.2 Calendar time

On the calendar timescale we consider the ordered time points at which an event took place (or end of study) 0.5, 0.8, 1.2. Once again, we do not consider time points beyond the last observed event time. Next, we split the follow-up into intervals using these time points as cutpoints:

{\{(τ_{j - 1}, τ_{j}]\}}_{j = 1}^{J} = {(0, 0.5], (0.5, 0.8], (0.8, 1.2]} .

The resulting data in the piece-wise exponential format is given in Table 3.

Table 3

Recurrent event data in PEM data format on the calendar timescale.

i	t_{i k} (cal)	k	j	(τ_j−1, τ_j]	δ_ijk	t_ijk	o_ijk = log (t_ijk)
1	0.5	1	1	(0,0.5]	1	0.5	log(0.5)
1	3.0	2	2	(0.5,0.8]	0	0.3	log(0.3)
1	3.0	2	3	(0.8,1.2]	0	0.4	log(0.4)
2	0.8	1	1	(0,0.5]	0	0.5	log(0.5)
2	0.8	1	2	(0.5,0.8]	1	0.3	log(0.3)
2	1.2	2	3	(0.8,1.2]	1	0.4	log(0.4)

4 Simulation studies

4.1 Comparison of PAMM and SFM in a proportional hazards setting

Under the assumption of proportional hazards the proposed PAMM and the SFM should give similar estimates of the regression parameters and the frailty variance. In order to check whether this is true an extensive simulation study was performed. The simulation settings as well as the simulation results are given in Supplement 1. In summary, it is shown that under the proportional hazards assumption the PAMM and the SFM give (almost) identical results for the estimated regression parameters in both gap- and calendar timescales, using the stratified and/or unstratified models. Furthermore, both models tend to estimate the regression parameters and the frailty variance more accurately when there is a longer follow-up time or there are more recurrences. It must be noted that when the frailty variances are large both models underestimate the frailty variance, but this underestimation is lower for longer follow-up times. Specifically, in the gap time scenarios, we see that the underestimation is worse in the SFM, and in the calendar time scenarios the underestimation is worse in the PAMM. The consequences of these biases of the estimated frailty variance on the estimates of the regression coefficients are minimal. It should be noted that the varying results for the frailty variances could depend on model implementation, and not necessarily on the model formulation, where mgcv uses restricted maximum likelihood to maximize the penalized log likelihood and coxme uses a Laplace approximation to calculate an integrated penalized likelihood.

4.2 Comparison of PAMM and SFM in distinguishing frailty and time-varying effects

Balan and Putter (2019) showed in simulation studies that in certain scenarios of clustered survival data – including recurrent events, it is difficult to distinguish the presence of a frailty effect from the presence of a time-varying effect. We replicated their simulation study for recurrent events to compare how the PAMM and the SFM distinguishes between a frailty effect and a time-varying effect. The only difference is that we used log-normal frailties for the Cox-like SFM, since the implementation of PAMM using mgcv is limited to Gaussian distributed random effects only. The results are shown in Supplement 4. In summary, the results indicate that the PAMM distinguishes no worse (and in some cases better) between a frailty effect and a time-varying effect than the Cox-like SFM with Gaussian distributed random effects for at least two events per subject.

5 Motivating examples

In this section we apply the PAMM for recurrent events to two real-life data sets using R (R Core Team, 2020; RStudio Team, 2020). We use the R packages pammtools (Bender and Scheipl, 2018) for the data transformation and the visualization of the results and mgcv (Wood, 2017, 2003, 2004, 2011; Wood et al., 2016) for model estimation. The aim of this section is to show the application of several types of recurrent events models using the PAMM framework.

5.1 Example 1: Staphylococcus aureus infection data

The staphylococcus aureus (SA) infection data in Abdulgader et al. (2019) contains times at which 137 children were colonized with a new staphylococcus aureus infection in the nasopharynx in their first year of life in Paarl, South Africa. In Abdulgader et al. (2019), the gap time PWP model (Prentice et al., 1981) was used to model the survival curves for the different infections (i.e., first, second, third, etc.). A standard CPH model was also used to investigate the effects of risk factors on the time to the first infection. In this section, we aim to repeat their analysis but include recurrent infections using the stratified PAMM (see equation (3.5)) in gap time. We consider a different baseline hazard for the first event and for recurrent events (combined). An advantage of using the PAMM here is that we can easily model, estimate and visualize the baseline hazards stratified for the first and recurrent events. We use the HIV exposure variable, where some children are HIV exposed and uninfected (HEU) and the rest are HIV unexposed (HU). First we look at the baseline hazards and then we look at a proportional hazards model for HIV exposure where the effect of HIV exposure is different for first and recurrent events.

5.1.1 Baseline hazards

We first fit the model without the HIV covariate but including the frailty term (i.e., λ(t) = λ₀ _,k (t) exp (z_i)) to have an understanding of the functional form of the baseline hazard. The estimated frailty variance is 0.066 (p = 0.249). This is a small frailty variance and so we will exclude the frailty term and reduce our model to a simpler model.

The estimated degrees of freedom, explained in Supplement 2, for the evolution of the hazard rate over time is different for first and recurrent SA infection, indicating that the evolution of the baseline hazards for first and recurrent events may be different (note: we are not testing for a difference in the baseline hazards). We visualize the estimated baseline hazard in Figure 1, where the hazard rates are re-scaled to correspond with the number of events per child-year. For the first SA event, the hazard rate is highest after birth and decreases quickly until around 5 months, thereafter a very slow decline is seen. For recurrent SA infections, we see an increase in the hazard over time since the previous infection until about three months. Thereafter we see an almost constant hazard for a recurrent infection and a slow decline starting from around 9 months.

Figure 1

The hazard rates for (A) first and (B) recurrent SA infections respectively. 95% confidence intervals are shaded in gray.

5.1.2 PH model

The hazard function of the model we are interested in, is

λ (t) = λ_{0, k} (t) \exp (β_{1} \times HIV)

where $λ_{0, k} (t) = \exp (β_{0, k} + f_{0, k} (t))$ and t is the time in the gap timescale. The case k = 1 indicates the first SA infection, and k = 2 indicates recurrent SA infection. The results indicate that the HEU children had a higher hazard of SA infection over time than HU children (HR = exp(β₁) = 1.31, 95% CI : 0.98 − 1.74).

In the paper by Abdulgader et al. (2019), the univariable analysis of HIV exposure on the time to first SA infection showed a HR effect close to 1. Next, we model an interaction of HIV exposure and the recurrence indicator to allow a separate HIV exposure effect for first and recurrent events, that is,

λ (t) = λ_{0, k} (t) \times \exp (β_{1, k} \times HIV)

(5.1)

where $λ_{0, k} (t), k$ and t are as defined above. Here we find that the HR for HIV exposure for the first SA infection is 0.98 (95% CI : 0.64 − 1.50)) (similar to what was reported in their paper) and for recurrent SA infections is 1.72 (95% CI : 1.16 − 2.54)).

5.2 Example 2: Childhood malaria data

To illustrate application of PAMMs on calendar timescale, we reanalyze data from a study of childhood malaria undertaken by (Kakuru et al., 2020) in Uganda to study the effects of two intermittent preventive treatments, with Sulfadoxine-pyrimethamine (SP) and Dihydroartemisininpiperaquine (DP) on the incidence of childhood malaria from birth until children are one year old. Therefore, the timescale for events and recurrences coincides with the age of the children in this example. In the original analysis, treatment effect was investigated with respect to the time to first episode of malaria, and it was found that there may be larger treatment differences between boys but not between girls. Here we analyze time to first event and additionally recurrences and add preterm birth and gravidity (number of previous pregnancies) as additional risk factors in the model. We further illustrate incorporation of cyclic splines and bivariate smooth functions by incorporating a seasonal effect and later its interaction with child age. An advantage using PAMMs, similar to Example 5.1, is that the baseline hazards are flexibly estimated, while simultaneously estimating a non-linear effect of seasonality, all in a SFM setting. Furthermore, it can easily be investigated, through more complex interactions, if there is any indication of the non-linear seasonality effect changing over time). Additionally in Supplement 3 we show the effects of different choices of basis functions for the baseline hazards, as well as degrees of freedom.

5.2.1 Model specification

We fit the following model

\begin{matrix} \log (λ (t ∣ x_{i}, z_{i})) = β_{0} + f_{0} (t_{j}) + f ({doy}_{i}) \\ + β_{1} {DP}_{i} + β_{2} {female}_{i} + β_{3} {DP}_{i} \cdot {female}_{i} \\ + β_{4} {preterm}_{i} + β_{5} {gravidity}_{i} + z_{i} \end{matrix}

(5.2)

where $β_{0} + f_{0} (t_{j})$ is the log-baseline hazard as defined in Section 3 and Equation 3.5 (here equivalent to the smooth effect of child age). We model the effect of seasonality via variable ‘day of year’ (doy) and set $f (doy) = \sum_{m = 1}^{10} γ_{m} b_{m} (doy)$ , where b_m are cyclic cubic regression spline bases that are defined such that the function value is the same at the beginning and end of the year (see Wood, 2017, Section 5.3.2).

The results show a significant frailty variance ( $({\hat{σ}}^{2} = {0.69}^{2} = 0.48, p < 0.001)$ ), which indicates that the frailty term is needed to account for the correlation between events in a child’s longitudinal profile.

The results of the analysis show a significant treatment effect for boys, where DP is associated with a lower hazard of malaria over time compared with SP ( $H R = \exp ({\hat{β}}_{1})$ ) (95% CI): 0.74 (0.58, 0.95)), but not for girls ( $H R = \exp ({\hat{β}}_{1} + {\hat{β}}_{3})$ (95% CI): 0.98 (0.78, 1.24)). The model also shows a significant effect of gravidity (HR (95% CI): 1.05 (1.005, 1.09)) indicating that children whose mothers had one more prior pregnancy had an approximately 5% higher hazard over time.

The estimated non-linear effects of age and seasonality are depicted in Figure 2. The upper row shows the effect of the child’s age as log-hazard contribution of the non-constant baseline term f₀(t_j) (left panel) and as hazard/incidence rates (episodes per child year) (right panel). The hazard rates were calculated by using varying time values while holding all other covariate values constant.

Figure 2

Upper row: The non-constant part of the log-baseline (f₀(t)) including confidence intervals (left panel) and respective hazard/incidence rates (episodes per child year, other covariate held constant; right panel) across age. Bottom row: Effect of seasonality as contribution to the log-hazard (f (doy), left panel) and hazard/incidence rates (episodes per child year, given child of age 100 days, other covariates as before; right panel).

Specifically, doy=1, treatment=DP, sex=Female, preterm=no and gravidity=1. From both graphs, we can see that the hazard of contracting malaria rises quickly as children get older until approximately 3 months of age, where the effects starts to level off. The bottom row of Figure 2 depicts the seasonality of malaria infection, where once again, the left panel depicts the log-hazard contribution of the term f (doy) and the right panel the hazard/incidence rates (episodes per child year, for children of varying levels of seasonality, 100 days of age, other covariates fixed as before). They show that the hazard rates go up from the start of the rainy seasons (January) and start to go down from the start of the dry seasons, with the second rainy season of the year (around June) showing a larger effect on the hazard of malaria. Per construction, the effect is the same at the end of December as in the beginning of January.

5.2.2 Modelling age-dependent effects of seasonality:

As children get older, they might be left unattended for longer periods of time and potentially spent more time outside or otherwise unprotected (e.g., by nets). We therefore additionally investigate a potential interaction between child age and seasonality when it comes to the hazard of malaria infection, for example such that seasonality does not affect the hazard for very young children as much as for older children. Since both variables are continuous, we model the interaction via bivariate penalized splines represented by tensor products, such that Equation 5.2 becomes

\begin{matrix} \log (λ (t ∣ x_{i}, z_{i})) = β_{0} + f_{0} (t_{j}, d o y_{i}) + β_{1} {DP}_{i} + β_{2} s e x_{i} + β_{3} {DP}_{i} \cdot s e x_{i} \\ + β_{4} {preterm}_{i} + β_{5} {gravidity}_{i} + z_{i} \end{matrix}

(5.3)

Note that in Equation 5.3 the non-linear part of the baseline hazard is now a function of both, age (t_j) and day of year (doy). This term can thus be viewed as a time-varying effect of the day of year, as it depends on t_j. In this case, however, it can also be viewed as an example for the application of multiple timescales, as the hazard depends on age (time since origin) and the day of year (calendar time).

The estimated bivariate effect f₀(t_j, doy) is depicted in Figure 3 (as before, similar representations could be obtained for the hazard/incidence rate where the other covariates are fixed at selected values). The left panel shows the surface of the bivariate function, with brighter colours indicating larger hazards. The middle panel shows vertical slices through the surface for fixed values of age, while the right panel shows horizontal slices through the surface for fixed values of day of year (season). These results are in line with the results obtained from model 5.2 and in this case indicate, that there is no substantial interaction between child age and seasonality effect. This conclusion could also be obtained more formally, using a ANOVA type decomposition of the bivariate effect into main and interaction effects (Wood, 2006), as illustrated in Suppplement.

Figure 3

Note: Brighter colors indicate increased hazards. Horizontal and vertical dashed lines indicate slices through the surface shown in the mid and right panel respectively. Mid panel: Slices through the bivariate surface for child ages 1, 4 and 10 months. Right panel: Slices through the bivariate surface at 1 March, 1 June and 1 November for the day of the year.

6 Discussion

In this article we introduced PAMMs as an alternative modelling approach for recurrent events data. These models utilize the generalized additive mixed modelling framework for estimation and inference. This facilitates specification and estimation of models for recurrent event data with complex covariate effects, including (multivariate) non-linear effects, time-varying effects and covariates, multiple timescales as well as random effects. We illustrated the application of PAMMs and specification of different covariate effects on real data sets.

In simulation studies we have shown that in the proportional hazards setting, SFM PAMMs give similar results to the Cox SFM with respect to the estimated effects and frailty variance. In certain scenarios of clustered survival data (including recurrent events) it is difficult to distinguish between unobserved heterogeneity and time-varying effects (non-proportional hazards) in the SFM, especially for smaller cluster sizes (where a cluster size of one implies univariate survival data) (Balan and Putter, 2019). Such identifiability issues result from the data generating process and will persist independently of the model class chosen for the analysis of such data. However, from our simulation study in Supplement 4, we saw that the PAMM has a lower type 1 error rate for detecting a frailty effect in the presence of time-varying effects, even for clusters as small as size two, compared with the Cox SFM. This difference is much more pronounced for larger sample sizes. Also notable in the simulation results, the PAMM may be subject to inflated type 1 errors when used for testing time-varying effects, regardless of whether or not there is unobserved heterogeneity. This simulation study was however limited to a single binary covariate, monotonic shapes of the baseline hazards and without stratification. A considerable amount of work, possibly in terms of a much larger simulation study with a focus on recurrent events, must be carried out to assess these preliminary findings more carefully.

Compared to parametric SFMs, SFMs estimated by PAMMs do not make any assumption about the distribution of event times. Compared to Cox based SFMs, PAMMs directly model the (baseline) hazard via penalized splines. Thus, while for Cox models (and piece-wise exponential models for that matter) the hazard can have implausible jumps between two subsequent time points, PAMMs enforce a penalty on large differences between hazards from neighboring time points (intervals), leading to stable and more plausible hazard estimates.

While PAMMs offer a flexible framework for recurrent events analysis, some readers might be concerned with the expansion of the data set that results from the data augmentation step and, consequently, the computation time. In our experience this is rarely an issue. Especially for small and medium sized data sets, the computation time is barely noticeable (even when using all unique event times as cut-points). For example, on a standard laptop with Intel core-i7 CPU with 1.90GHz, 16GB RAM, fitting the interaction model in Section 5.1 (Equation (5.1)) takes about 1 second; the model estimating incidences of Malaria infection (with bivariate effects, Equation (5.3)), in Section 5.2, takes about 10 seconds. For high-dimensional data efficient estimation methods exist (Wood et al., 2017; Greven and Scheipl, 2020). For models on the gap timescale, computation time could be further improved without sacrificing predictive performance by reducing the number of interval cut points (cf. complexity analysis in Bender et al., 2021).

The R package pammtools has been extended to support recurrent events analysis, including the required data augmentation on the gap as well as calendar timescale. The convenience functions in pammtools allow the calculation of covariate effects, hazards rates, hazard ratios, cumulative hazards and survival probabilities according to the user’s specifications. This means that these quantities, including confidence intervals, can be estimated and visualized in a few simple steps, even in case of complex models (cf. Section 5, Supplement 3).

Another advantage of the PAMM framework is that it does not require any particular implementation in order to estimate the model. While pammtools provides many convenience functions, once the data is transformed, any software that can maximize the Poisson likelihood can be used for model estimation. However, some methods and implementations will be more suitable than others. In this article, we estimated PAMMs using the R package mgcv, which offers great support for complex modelling of covariate effects. Nevertheless, in some settings, it might have limitations, that can be compensated by using a different software for estimation. For example, while mgcv offers support for random effects estimation, it is inefficient when random effects have high cardinality. In such cases, specialized mixed modelling software might be preferred, for example, lme4 (Bates et al., 2015), nlme (Pinheiro et al., 2020) or gamm4 (Wood and Scheipl, 2020). The mixed modelling framework implemented in these packages can further be exploited to estimate hierarchical models with nested or crossed random effects or more general hierarchical generalized additive models, that can be estimated via empirical Bayes (Pedersen et al., 2019) or fully Bayesian techniques (Bu¨rkner, 2017).

A current limitation of these techniques and respective implementations is that they are limited to modelling Gaussian distributed random effects, while Gamma distributed random effects are popular in the context of survival analysis (Balan and Putter, 2020). However, through simulation studies by Gasparini et al. (2019); Liu et al. (2017), it has been shown that the estimates of the regression coefficients are quite robust to misspecification in the choice of frailty distributions regardless of sample size or the amount of heterogeneity as determined by the frailty variance. Gasparini et al. (2019) showed that it is more important to prioritize modelling the baseline hazard correctly. Alternatively, package hglm could be used to estimate Gamma distributed random effects, however, it has no native support for penalized spline estimation. Furthermore, the gamlss package (Rigby and Stasinopoulos, 2005), offers estimation of non-parametric random effects (Coupé, 2018). Both options, however, need to be explored further in the context of survival and recurrent events analysis.

This work highlights the usefulness and possibilities of extending PAMMs to recurrent events data.In future research, this model class could be further extended to settings with multivariate recurrent events (i.e., when multiple time-to-event outcomes can (re)occur on the same observation unit), recurrent events data with a terminal or competing event or for recurrent events data with non-zero time to recovery.

7 Availability of data and materials

Data for the Staphylococcus aureus example are contained within the pammtools package and are also available from the publisher’s webpage as Supplement 5. The malaria data is not available for public sharing but can be accessed as indicated in the data availability section in the manuscript Kakuru et al. (2020). We did however simulate a malaria dataset with 200 subjects from the original data where we changed the peaks of the seasonality effect. This data is available in Supplement 6. Additionally, the analysis presented here is accompanied by an online supplement that illustrates the data structure and code for data transformation, model fitting and visualization (cf. Supplements 2 and 3) in detail. Further, a vignette for recurrent events PAMMs can be found at https://adibender.github.io/pammtools/articles/recurrent-events.html. The R code and supplementary materials for the simulation study in Supplement 1 can be found in https://github.com/jordache-ramjith/PAMM and doi.org/10.6084/m9.figshare.14638353.v1 respectively. The R code for the simulation data in Supplement 4, can be requested from the authors Balan and Putter (2019) directly.

Supplementary material

Supplementary material is available online.

Supplemental Material for Recurrent events analysis with piece-wise exponential additive mixed models by Jordache Ramjith, Andreas Bender, Kit C. B. Roes, Marianne A. Jonker, in Statistical Modelling

Footnotes

Acknowledgements

We thank Shima Mohammed, Mark Nicol and Heather Zar for access and support to use the staphylococcus aureus infection data from the Drakenstein Child Health Study. We also thank Abel Kakuru, Grant Dorsey and the rest of the authors in the Uganda malaria birth cohort study for access and support to use the data for secondary analysis. We would also like to thank Theodore Balan for sharing his simulation code, so that we could replicate the simulations for our own comparisons. Lastly we would like to thank the reviewers for their thorough feedback and suggestions.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

AB was funded by the German Federal Ministry of Education and Research (BMBF) under Grant No. 01IS18036A.

References

Abdulgader

, Robberts

, Ramjith

, Nduru

, Dube

, Gardner-Lubbe

, Zar

, and Nicol

(2019) Longitudinal population dynamics of staphylococcus aureus in the nasopharynx during the first year of life. Frontiers in Genetics , 10, 198.

Akinboyo

, Voskertchian

, Gorfu

, Betz

, Ross

, Carroll

, and Milstone

(2018) Epidemiology and risk factors for recurrent staphylococcus aureus colonization following active surveillance and decolonization in the nicu. Infection Control and Hospital Epidemiology , 39, 13349.

Amorim

, and Cai

(2015) Modelling recurrent events: a tutorial for analysis in epidemiology. International Journal of Epidemiology , 44, 324–33.

Andersen

, and Gill

(1982) Cox’s regression model for counting processes: a large sample study. The Annals of Statistics , 10, 1100–20.

Argyropoulos

, and Unruh

(2015) Analysis of time to event outcomes in randomized controlled trials by generalized additive models. PLoS ONE , 10, e0123784. https://doi.org/10.1371/journal.pone.0123784.

Balan

, and Putter

(2020) A tutorial on frailty models. Statistical Methods in Medical Research , 29, 3424–54.

Balan

, and Putter

(2019) Nonproportional hazards and unobserved heterogeneity in clustered survival data: When can we tell the difference? Statistics in Medicine , 38, 3405–20.

Barakat

, Hussein

, Fouda

, Zoair

, and Abd El-Razek

(2021) The underlying causes of recurrent pneumonia in children: A two-center study. Journal of Advances in Medicine and Medical Research , 33, 62–9.

Bates

, Ma¨chler

, Bolker

, and Walker

(2015) Fitting linear mixed-effects models using lme4. Journal of Statistical Software , 67, 1–48. https://doi.org/10.18637/jss.v067.i01.

10.

Bender

, and Scheipl

(2018) pammtools: Piece-wise exponential Additive Mixed Modeling tools . https://doi.org/10.48550/arXiv.1806.01042.

11.

Bender

, Groll

, and Scheipl

(2018) A generalized additive model approach to time-to-event analysis. Statistical Modelling , 18, 299–321.

12.

Bender

, Scheipl

, Hartl

, Day

, and Ku¨chenhoff

(2019) Penalized estimation of complex, non-linear exposurelag-response associations. Biostatistics , 20, 315–31. ISSN 1465-4644. https://doi.org/10.1093/biostatistics/kxy003.

13.

Bender

, Ru¨gamer

, Scheipl

, and Bischl

(2021) A general machine learning framework for survival analysis. In Machine Learning and Knowledge Discovery in Databases , pages 158–173. Springer International Publishing. https://arxiv.org/abs/2006.15442.

14.

Box-Steffensmeier

, and De Boef

(2006) Repeated events survival models: the conditional frailty model. Statistics in Medicine , 25, 3518–33.

15.

Bu¨rkner

P-C

(2017) brms: An r package for bayesian multilevel models using stan. Journal of Statistical Software , 80, 1–28. ISSN 1548-7660. https://doi.org/10.18637/jss.v080.i01.

16.

Carmona-Bayonas

, Jiménez-Fonseca

, Lamarca

, Barriuso

, Castan˜o

, Benavent

, Alonso

, Riesco-Martínez

MdC

, Alonso-Gordoa

, Custodio

, SÁnchez CÁnovas

, Hernando Cubero

, López

, Lacasta

, FernÁndez Montes

, Marazuela

, Crespo

, Escudero

, Diaz

JÁ

, Feliciangeli

, Gallego

, Llanos

, Segura

, Vilardell

, Percovich

, Grande

, Capdevila

, Valle

, and García-Carbonero

(2019) Prediction of progression-free survival in patients with advanced, well-differentiated, neuroendocrine tumors being treated with a somatostatin analog: The GETNE-TRASGU study. Journal of Clinical Oncology , 37, 2571–80.

17.

Chan

, Ko

FWS

, Ling

, Cheung

, Chan

, Lo

, Ng

, MM-s

Lui

, Yee

KSW

, . (2021) A territory-wide study on the factors associated with recurrent asthma exacerbations requiring hospitalization in hong kong. Immunity, Inflammation and Disease , 9, pages 1–13.

18.

Coupé

(2018) Modeling linguistic variables with regression models: Addressing non-gaussian distributions, non-independent observations, and non-linear predictors with random effects and generalized additive models for location, scale, and shape. Frontiers in psychology , 9, 513.

19.

Cox

(1972) Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological) , 34, 187–202.

20.

Dos Santos

, de Góis Filho

, Silva

AMF

, Santos

JVG

, Santos

, Aquino

, de Jesus

, Almeida

MLD

, da Silva

, Altmann

, . (2021) Recurrent covid-19 including evidence of reinfection and enhanced severity in thirty brazilian healthcare workers. Journal of Infection , 82, 399–406.

21.

Duchateau

, Janssen

, Kezic

, and Fortpied

(2003) Evolution of recurrent asthma event rate over time in frailty models. Journal of the Royal Statistical Society: Series C (Applied Statistics) , 52, 355–63.

22.

Eilers

PHC

, and Marx

(1996) Flexible smoothing with B-splines and penalties. Statistical Science , 11, 89–121. https://doi.org/10.1214/ss/1038425655.

23.

Feng

, He

, Liu

, Chen

, and He

(2021) Effect of congenital heart disease on the recurrence of cough variant asthma in children. BMC Cardiovascular Disorders , 21, 1–9.

24.

Freedman

(2006) On the so-called “huber sandwich estimator” and “robust standard errors”. The American Statistician , 60, 299–302.

25.

Gasparini

, Clements

, Abrams

, and Crowther

(2019) Impact of model misspecification in shared frailty survival models. Statistics in medicine , 38, 4477–502.

26.

Ghosh

, Olaniyi

, and Obabiyi

(2020) Mathematical analysis of reinfection and relapse in malaria dynamics. Applied Mathematics and Computation , 373, 125044.

27.

Giacoppo

, Alfonso

, Xu

, Claessen

BEPM

, Adriaenssens

, Jensen

, Pérez-Vizcayno

, D-Y

Kang

, Degenhardt

, Pleva

, Baan

, Cuesta

, D-W

Park

, Kukla

, Jiménez-Quevedo

, Unverdorben

, Gao

, Naber

, S-J

Park

, Henriques

JPS

, Kastrati

, and Byrne

(2020) Drug-coated balloon angioplasty versus drug-eluting stent implantation in patients with coronary stent restenosis. Journal of the American College of Cardiology , 75, 2664–78.

28.

Gibson

, Pinto

, Chi

, Arbetter

, Yee

, Mehran

, Bode

, Halperin

, Verheugt

, Wildgoose

, . (2017) Recurrent hospitalization among patients with atrial fibrillation undergoing intracoronary stenting treated with 2 treatment strategies of rivaroxaban or a dose-adjusted oral vitamin k antagonist treatment strategy. Circulation , 135, 323–33.

29.

Greven

, and Scheipl

(2020) Comments on: Inference and computation with generalized additive models and their extensions. TEST , 29, 343–50.

30.

Hastie

, and Tibshirani

(1986) Generalized Additive Models. Statistical Science , 1, 297–310. ISSN 0883-4237. https://doi.org/10.1214/ss/1177013604.

31.

Hernstadt

, Cheung

, Hurem

, Vasilunas

, Phuong

, Quinn

, Agrawal

, Daley

, Cole

, and Gwee

(2020) Changing epidemiology and predisposing factors for invasive pneumococcal disease at two australian tertiary hospitals. The Pediatric infectious disease journal , 39, 1–6.

32.

Hougaard

(2012) Analysis of multivariate survival data . Springer Science & Business Media.

33.

Huber

(1967) The behavior of maximum likelihood estimates under nonstandard conditions. In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability , volume 1, pages 221–233. University of California Press.

34.

Hutton

, and Monaghan

(2002) Choice of parametric accelerated life and proportional hazards models for survival data: asymptotic results. Lifetime Data Analysis , 8, 375–93.

35.

Iacobelli

, and Carstensen

(2013) Multiple time scales in multi-state models. Statistics in Medicine , 32, 5315–27. https://doi.org/10.1002/sim.5976.

36.

Jahn-Eimermacher

(2008) Comparison of the andersen–gill model with poisson and negative binomial regression on recurrent event data. Computational Statistics & Data Analysis , 52, 4989–97.

37.

Kakuru

, Jagannathan

, Kajubi

, Ochieng

, Ochokoru

, Nakalembe

, Clark

, Ruel

, Staedke

, Chandramohan

, . (2020) Impact of intermittent preventive treatment of malaria in pregnancy with dihydroartemisinin-piperaquine versus sulfadoxine-pyrimethamine on the incidence of malaria in infancy: a randomized controlled trial. BMC Medicine , 18, 1–11.

38.

Kelly

, and LL-Y

Lim

(2000) Survival analysis for recurrent event data: an application to childhood infectious diseases. Statistics in Medicine , 19, 13–33.

39.

Kwong

, and Hutton

(2003) Choice of parametric models in survival analysis: applications to monotherapy for epilepsy and cerebral palsy. Journal of the Royal Statistical Society: Series C (Applied Statistics) , 52, 153–68.

40.

Lawpoolsri

, Sattabongkot

, Sirichaisinthop

, Cui

, Kiattibutr

, Rachaphaew

, Suk-Uam

, Khamsiriwatchara

, and Kaewkungwal

(2019) Epidemiological profiles of recurrent malaria episodes in an endemic area along the thailand-myanmar border: a prospective cohort study. Malaria Journal , 18, 1–11.

41.

Y-H

, Klein

, and Moeschberger

(1996) Effects of model misspecification in estimating covariate effects in survival analysis for small sample sizes. Computational Statistics & Data Analysis , 22, 177–92.

42.

Lin

, and Zhang

(1999) Inference in generalized additive mixed modelsby using smoothing splines. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 61, 381–400.

43.

Liu

, Pawitan

, and Clements

(2017) Generalized survival models for correlated time-to-event data. Statistics in Medicine , 36, 4743–62.

44.

Munda

, Rotolo

, and Legrand

(2012) parfm: Parametric frailty models in r. Journal of Statistical Software , 51, 1–20.

45.

Noordzij

, van Diepen

, Caskey

, and Jager

(2017) Relative risk versus absolute risk: one cannot be interpreted without the other. Nephrology Dialysis Transplantation , 32, ii13–18.

46.

Panel

ADC

, Pen

, Pannebakker

, Helsen

HHM

, and Wertheim

(2020) Seasonal morphotypes of drosophila suzukii differ in key life-history traits during and after a prolonged period of cold exposure. Ecology and Evolution , 10, 9085–99.

47.

Pedersen

, Miller

, Simpson

, and Ross

(2019) Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ , 7, https://doi.org/10.7717/peerj.6876.

48.

Pinheiro

, Bates

, DebRoy

, Sarkar

, and Core Team

(2020) nlme: Linear and Nonlinear Mixed Effects Models . URL https://CRAN.R-project.org/package=nlme. Prentice

, Williams

and Peterson

(1981) On the regression analysis of multivariate failure time data. Biometrika , 68, 373–9.

49.

Putter

, Sasako

, Hartgrink

, Van de Velde

, and Van Houwelingen

(2005) Long-term survival with non-proportional hazards: Results from the dutch gastric cancer trial. Statistics in Medicine , 24, 2807–21.

50.

Pya

, and Wood

(2015) Shape constrained additive models. Statistics and Computing , 25, 543–59. https://doi.org/10.1007/s11222-013-9448-7.

51.

Core Team

(2020) R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

52.

Ramjith

, Roes

, Zar

, and Jonker

(2021) Flexible modelling of risk factors on the incidence of pneumonia in young children in south africa using piece-wise exponential additive mixed modelling. BMC Medical Research Methodology , 21. https://doi.org/10.1186/s12874-020-01194-6.

53.

Rigby

, and Stasinopoulos

(2005) Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics) , 54, 507–54.

54.

Team

RStudio

(2020) RStudio: Integrated Development Environment for R . RStudio, PBC., Boston, MA. URL http://www.rstudio.com/.

55.

Ruhe

(2018) Quantifying change over time: Interpreting time-varying effects in duration analyses. Political Analysis , 26.

56.

Therneau

, and Hamilton

(1997) rhdnase as an example of recurrent event analysis. Statistics in medicine , 16, 2029–47.

57.

Ullah

, Gabbett

, and Finch

(2014) Statistical modelling for recurrent events: an application to sports injuries. British Journal of Sports Medicine , 48, 1287–93.

58.

Varma

, Bourge

, Stevenson

, Costanzo

, Shavelle

, Adamson

, Ginn

, Henderson

, Abraham

and CHAMPION Investigator Group (2021) Remote hemodynamic-guided therapy of patients with recurrent heart failure following cardiac resynchronization therapy. Journal of the American Heart Association , 10, e017619.

59.

Wand

(2003) Smoothing and mixed models. Computational statistics , 18, 223–49.

60.

Wei

L-J

, Lin

, and Weissfeld

(1989) Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of the American Statistical Association , 84, 1065–73.

61.

Wey

, Hart

, Salkowski

, Skeans

, Kasiske

, Israni

, and Snyder

(2020) Posttransplant outcome assessments at listing: Long-term outcomes are more important than short-term outcomes. American Journal of Transplantation , 20, 2813–21.

62.

Wood

(2003) Thin-plate regression splines. Journal of the Royal Statistical Society (B) , 65, 95–114.

63.

Wood

(2004) Stable and efficient multiple smoothing parameter estimation for generalized additive models. Journal of the American Statistical Association , 99, 673–86.

64.

Wood

(2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) , 73, 3–36.

65.

Wood

, and Scheipl

(2020) gamm4: Generalized Additive Mixed Models using ‘mgcv’ and ‘lme4’ . URL https://CRAN.R-project.org/package=gamm4.

66.

Wood

(2006) Low-Rank Scale-Invariant Tensor Product Smooths for Generalized Additive Mixed Models. Biometrics , 62, 1025–36.

67.

Wood

(2017) Generalized additive models: an introduction with R . CRC press.

68.

Wood

, Li

, Shaddick

, and Augustin

(2017) Generalized additive models for gigadata: modeling the uk black smoke network daily data. Journal of the American Statistical Association , 112, 1199–1210.

69.

Wood

, Pya

, and Sa¨fken

(2016) Smoothing parameter and model selection for general smooth models (with discussion). Journal of the American Statistical Association , 111, 1548–75.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.48 MB

0.27 MB

0.75 MB

0.10 MB

0.01 MB

0.02 MB