Frailty modelling for multitype recurrent events in clinical trials

Abstract

Recurrent event outcomes are ubiquitous among clinical trial data which encourages a conventional approach to analysis. Yet a common feature of these data has received less attention, that is, survival times often comprise multiple types of events that may imply a disparity in cost and disease severity. Typically, we neglect this feature of the data by combining event-types or analyzing each type separately, thus ignoring any interdependence among them. This practice may reflect a dearth of readily available methods and software that more appropriately acknowledge the true data structure. We provide a review of the literature on multitype recurrent events and frailty modelling which reflects a renewed interest in the topic over the past decade and the emergence of software for estimation. Thus, a review of available methods seems timely, if not overdue.

Keywords

Frailty Mixed modelling multivariate random effects recurrent events

1 Introduction and motivating examples

In clinical trials of chronic diseases, we often have survival data including recurrent events that reflect disease progression and quality of life. The events may be upsetting and inhibiting for patients, such as migraines, hypoglycemic episodes in diabetes, hospital readmission in heart failure and so forth. These events may be classified by type, for example, according to severity for migraines, or whether hypoglycemia occurs during the daytime or nighttime, or whether the heart failure patient is admitted to the emergency department (ED) or hospital. Note that in all of these examples the occurrence of one type of event does not preclude the subsequent occurrence of another, that is, they are not competing risks, yet event-types may be interdependent. Recurrent events will often exhibit types (whether nominal or ordinal) in which case we can describe them as multitype recurrent events (MTREs). To further illustrate this feature of survival data, we describe several instances of MTREs we have encountered in our work.

In the mid-1990s, the Nambour Prevention Trial explored the effect of sunscreen use and beta-carotene on the incidence of nonmelanoma skin cancer in a high risk population in Northern Queensland, Australia (Green et al. 1999). About 1600 participants were randomized in a 2 × 2 factorial design, with new skin cancers detected at regular follow-up visits. Two types of skin cancers were noted: basal cell carcinoma (BCC) and squamous cell carcinoma (SCC). These are distinct histological types with different prognoses and implications (SCCs are a more serious and less common lesion) and, thus, we may not want to combine these event-types using a time-to-first BCC/SCC analysis. Thus, event-types are analyzed separately, using Cox regression, despite the very low failure rate (high censoring rate) especially for SCCs (>90%, possibly influenced by poor compliance), which renders statistical significance unattainable.

A meta-analysis of nine parallel-group and crossover studies (1674 patients) was carried out comparing rates of hypoglycemia between human insulin 30 and biphasic insulin aspart 30 (an insulin analog) in patients with Type 2 diabetes (Davidson et al. 2009). Using negative binomial regression, hypoglycemic episodes were analyzed overall and separately for major, minor and symptoms only, and also for nocturnal and daytime hypoglycemia. Patients may experience multiple episodes during the study period and the rate varies for the different event-types, for example, major episodes are less common than minor episodes. The inverse variance method was used to combine treatment effect estimates from the parallel-group and crossover studies. Although rates of overall episodes were not significantly different between the treatment groups, rates were considerably different for nocturnal and major episodes (these are particularly worrying for patients), although their incidence is low even in a meta-analysis of nine studies (14% and 2%, respectively). Overall, this approach seems ambivalent and simplistic.

Figure 1:

Example heart failure readmissions data

More recently, we examined the prognostic value of electrocardiogram (ECG) parameters with respect to heart failure related readmissions in a Canadian cohort of 900 patients with acute heart failure (Gouda et al. 2016). We have up to five years of readmissions data for patients who were recruited in the ED. Roughly 30% of patients had at least one ED visit and 20% had at least one hospital visit during the follow-up period (i.e., after discharge). The study design itself insists a more than nominal distinction between ED and hospital admissions, implying different costs, length of stay and severity. Using Cox regression to analyse time-to-events, patients with and without atrial fibrillation at the index visit showed no difference in ED visits ( $p = 0.40$ ), hospital readmissions ( $p = 0.55$ ) or when these event-types were combined ( $p = 0.23$ ). In Figure 1, we provide example readmissions data.

When confronted with data like those from these three studies, the statistician will typically either (a) consider event-types separately or (b) analyze them together. Regarding (a), if event-types are not independent, that is, if the occurrence of one type is suggestive of the risk of another type, then separate analyses are inadequate. Also, data may be sparse when separated out leading to equivocal results which may tempt us to consider analyzing types together. Regarding (b), although different event-types may all reflect disease progression, treating them as identical denies a potential disparity in cost, relationship to future outcomes and disease severity. Also, a strong treatment effect on one event-type will be diluted when combined with another event-type that shows a weaker effect. Thus, both approaches seem crude: we lose information if we analyze separately and we mask information if we combine them together. A more optimal modelling approach is needed.

Although a number of review articles of recurrent events analysis are available (Pandeya et al. 2005; Cook and Lawless 2007; Castaneda and Gerritse 2010; Rogers et al. 2014b) as far as we can tell no such overview exists for MTREs. A literature search reveals that models incorporating event-types are scarcely evident in medical journals reporting the results of clinical trials, with the topic confined to statistics journals. This implies limited uptake of the methods proposed; perhaps they are considered too esoteric or difficult to implement? Even among the statistics journals there appears to be a dearth of relevant research compared to the long, continuous publication history for univariate recurrent events. (Chen et al. 2012) stated that ‘statistical methods for handling multiple type recurrent events are relatively limited’ and (Zhu et al. 2010) noted that ‘there does not seem to exist an established method’. Although, with the publication of recent monographs (Cook and Lawless 2007) and programming code made available by authors (Rondeau et al. 2012), there now appears to be renewed interest in the topic (Li et al. 2015) and fewer obstacles to implementation. We focus our attention on proportional hazards frailty models because they are prevalent and a simple extension of more familiar models.

In the next section, we introduce the concept of frailty modelling and present an MTRE model by showing how it relates to common models for the analysis of recurrent events data. (The literature on recurrent events is extensive and we do not intend to provide an overview of that here.) In the third section, we describe multivariate frailty models for MTREs and then recent work regarding extensions of these models for particular scenarios. Finally, we cover software applications and use an example to illustrate how we would fit multivariate frailty model for MTREs in SAS.

2 From Cox regression to MTRE modelling

Survival data are ubiquitous in clinical trials and this encourages a conventional approach to analysis (Chun et al. 2012). For example, a Cox proportional hazards regression model (Cox 1972) is typically used to examine the influence of covariates such as treatment on survival outcomes.

λ (t; x) = λ_{0} (t) \exp (x β),

(2.1)

where

λ (t; x)

is the instantaneous probability of the occurrence of an event at time

t

, and we want to assess the effect of covariates (

x

) on this quantity.

λ_{0} (t)

is the baseline intensity corresponding to the case where covariates assume a value of 0, indicating the baseline or reference value. In the case of the heart failure study we may use a Cox model to analyse time-to-first readmission, or more commonly a composite of time-to-first readmission and death (Pocock et al. 2012). However, there is a growing demand to incorporate all hospital readmissions into the analysis because doing so entails greater statistical power and consequently a smaller required sample size or shorter follow-up (Anker and McMurray 2012; Rogers et al. 2014b). In this case, the Cox regression analysis of time-to-first becomes difficult to justify and readmissions may be analyzed using a popular extension of this model described by Andersen and Gill (1982).

Yet the Andersen–Gill model assumes events within a patient are independent, meaning that the likelihood of an event at time $t$ is unaffected by the event history and remains a function of the covariates and time $t$ only. Gamma random effects ( $θ_{i}$ for patient $i$ ) may be introduced to account for patient heterogeneity. These random effects have a multiplicative effect on the hazard leading to negative binomial regression (see Lawless 1987), analogous to the way Poisson regression is extended to negative binomial regression by including gamma random effects.

λ_{i} (t; x_{i}) = θ_{i} λ_{0} (t) \exp (x_{i} β) .

(2.2)

Despite these extensions to the original Cox model, an important feature of the event data is still ignored: readmissions are comprised of different types (e.g., ED and hospital visits), and the common approach of combining them should be questioned. Thus, we further extend our model by including random effects for event-types ( $ξ_{ij}$ for patient $i$ and event-type $j$ ), creating another mixed (i.e., fixed $(β)$ and random $(θ_{i}, ξ_{ij})$ ) effects model; see Abu-Libdeh et al. (1990).

λ_{ij} (t; x_{i}) = θ_{i} ξ_{ij} λ_{0} (t) \exp (x_{i} β) .

(2.3)

The relatedness of the recurrent events models given earlier is apparent from their form, that is, intensity-based, non-homogeneous Poisson process models (NHPP)—non-homogenous because the intensity is not constant in time. The $β$ are the regression coefficients for the covariates that correspond to, for example, the treatment effect, and we may specify a Weibull baseline intensity $λ_{0} (t) = ρ δ t^{δ - 1}$ because it is a popular choice and employed by both Lawless and Abu-Libdeh et al. (1990) noted earlier. Note if $δ = 1$ , then the intensity is constant and we have a time homogenous Poisson process. Note also that instead of the hazard we are speaking of the ‘intensity’ which is similar to the hazard, that is, it is the instantaneous probability of at least one event, although these terms may be used interchangeably. Other approaches could be noted (such as multi-state models, marginal models); however, the models given earlier are instructive because they lead us from the familiar Cox model to the model of interest with only small additions such as random effects. These random effects are referred to as ‘frailties’ and a number of monographs are available covering intensity-based frailty modelling, namely Wienke (2010), Aalen et al. (2008), Hougaard (2000), Andersen et al. (1993) and Duchateau and Janssen (2008). We will now introduce the concept of frailty modelling.

The assumption of independent events for the Andersen–Gill model implies that the baseline intensity is common for all events; thus a new event is unaffected by earlier events experienced by the patient (an assumption that is obviously violated in the context of heart failure hospitalizations). This differs from model (2.2) which includes random effects for patients ( $θ_{i}$ for patient i) to account for unexplained heterogeneity, that is, beyond that attributable to the covariates. (If we fail to account for heterogeneity then standard errors for $β$ will be underestimated and show bias that ‘grows with the correlation and the number of possible recurrences’; Villegas et al. (2013).) These random effects are termed ‘frailties’ because a higher value means a higher hazard, and events now occur according to each patient's own hazard, that is, events within a patient are correlated, with the independent frailties mimicking some probability distribution (typically gamma is assumed). The model is referred to as a ‘shared frailty’ model and is employed when we have clustered data such as animal experiments where a litter could be said to have a shared frailty, or in the heart failure study patients in the same site may assume a common frailty (as per O'Brien et al. 2014) with academic versus non-academic hospitals. Thus, we are concerned with describing associations between events. Since these scenarios are analogous to the situation of recurrent events clustered within a patient, we may adopt a shared frailty model for the analysis of recurrent events. An early and important paper making use of the shared frailty concept is Clayton (1978) who considered the risk of chronic disease in families (although the term ‘frailty’ was not used).

3 Multivariate frailty models for MTREs

The intensity-based frailty model (2.3) was described by Abu-Libdeh et al. (1990) (henceforth we will simply refer to as Abu-Libdeh) for the analysis of nonmelanoma skin cancer incidence in a randomized clinical trial of selenium supplements and includes random effects for patients $(θ_{i})$ and event-types ( $ξ_{ij}$ for patient i and event-type j). This is the earliest and most relevant paper concerning random effects modelling of multiple event-types; however, there is no apparent evidence that the model has since been employed in the recurrent clinical events setting. Incidentally, we may consider a more parsimonious model without the cluster effect for patients $(θ_{i})$ such that if event-types are equivalent, that is, $ξ_{ij} = ξ_{i}$ for $j = 1, \dots, J$ (where J is the number of types of events), then we have the simple shared frailty model. Although in this case the random effects for unexplained patient heterogeneity and dependence among events are confounded. Note that when patient-level effects enter the model it is natural to discuss random effects, rather than a fixed effect with very many levels, and we must choose a probability distribution that describes the random effects, and this distribution is defined by unknown parameters.

A number of probability distributions have been adopted in practice (see Lindeboom and Van Den Berg (1994) for some background). According to Cook and Lawless (2007), considerations for the choice of distribution include ‘tractability of the integral..., properties of the full intensity function, and the availability of software’ (see section ahead) and for event-types we desire a multivariate distribution. Abu-Libdeh assumed a gamma distribution for the patient effects (as per negative binomial regression) and a multivariate Dirichlet distribution for event-type effects (like the gamma distribution, the Dirichlet distribution is often used as a prior in Bayesian statistics). With two event-types (e.g., ED and hospitalization), the Dirichlet reduces to a beta distribution. However, it can be shown that, owing to the Dirichlet distribution (with $\sum_{j = 1}^{J} ξ_{ij} = 1$ for all i), identical treatment effect estimates would be obtained using the simpler shared frailty model that ignores event-type; Abu-Libdeh acknowledge as much when stating that ‘the parameters of the Dirichlet mixing distribution can be treated separately from the estimation of the other parameters’. The upside of this is that the Dirichlet enables the marginal likelihood to be derived analytically, which was likely a relevant concern when the paper was published in 1990. Nevertheless, an alternative distributional choice is desired.

A good option for our purpose is the multivariate log-normal distribution which is promoted and spelled out by Cook and Lawless (2007) for multiple event-types. It is more apt than the gamma when we have multivariate, correlated frailties. The random effects reflect correlation among events for a patient and the correlation between event-types is deduced from the multivariate distribution. Consider then, this alternative model:

λ_{ij} (t; x_{i}) = λ_{0 j} (t) \exp (x_{i} β_{j} + u_{ij}),

(3.1)

where

u_{ij} = log ξ_{ij}

and the vector

[u_{i 1}, . . ., u_{iJ}]

has a multivariate normal distribution, thus allowing positive and negative frailties (see, e.g., Chen and Cook 2009 and Chen et al. 2005). The approach is analogous to the so-called correlated lognormal frailty model which does not restrict units within a cluster to have a shared frailty (perfect positive correlation) but rather correlated frailties (Wienke 2010; Chapter 5). For example, in twin studies there will be separate random effects for each set of twins (clusters), yet we may desire that the individual frailties within a set of twins should be different but correlated, that is, via the multivariate lognormal distribution (Locatelli et al. 2004). Similarly, in the case of MTRE we do not want events within patients to have a shared frailty but wish to specify separate but dependent frailties for the event-type processes.

Regarding the formulation of model, as in model (2.3), we may have a common $β$ or separate regression coefficients for each event-type ( $β_{i}$ ) in the regression component depending on whether the assumption of a common effect across event-types is deemed reasonable. Likewise, we may have an event-type specific baseline hazard. Also note that compared to model (2.3), the patient-level random effects $θ_{i}$ are dropped in model (3.1) and thus patient and event-type variation is accounted for in the random effects $u_{ij}$ . Using a multivariate distribution for frailties $u_{ij}$ , we get a sense of the correlation among event-types which we would not otherwise obtain, thus providing new insights while accounting for the interdependence among event-types. More to the point, a model that ignored these associations may be inadequate or inefficient.

4 Other issues in MTREs: Terminal event and other developments

In this review we have so far neglected the likely possibility of a terminal event (competing risk) which adds further complexity to our model, that is, the follow-up time cannot be treated as independent of the recurrent event process. In the examples cited in the ‘Introduction’, it may be reasonable to treat deaths as uninformative censoring in one case, for example, nonmelanoma skin cancer incidence (thus, the competing risk of death is absent in Abu-Libdeh's analysis), but not in another, for example, repeat heart failure related hospitalizations (see Figure 1 where death follows readmissions). Such terminal events preclude further observations of the event of interest and in the latter example we should not treat these as censored survival times since the likelihood of death increases with the number of ED and hospital visits (Lee et al. 2009); thus these processes are not independent and we would like to account for the association between them (to avoid biased estimates). Other examples of terminal events may include consequences of ineffective treatment such as dropout due to adverse events or a switching of treatment. Liu and Huang (2008) state that ‘over the past decade, there has been a growing research interest in modelling correlated failure times in the presence of informative dropout or a dependent terminal event such as death’ which leads us again to the concept of frailty and random effects.

To introduce the ‘joint frailty’ model, let us first consider the simple case without multiple event-types. Joint frailty models are becoming increasingly popular as a means of accounting for terminal events in heart failure (e.g., Rogers et al. 2014b who promote its use, and Greenberg et al. 2014), and thus frailty as a concept and means of associating events is becoming familiar to clinicians. As Rogers et al. explain: ‘a common frailty term, which can be thought of as an unmeasured indication of the severity of illness that affects both hospitalisation rate and hazard for [cardiovascular] death, induces an association between the two processes’ (2014a). Effectively, this means the processes are jointly modelled and thus our clinical understanding of the recurrent events process and mortality is enhanced. Bear in mind we are now accounting for dependence in three ways: between recurrent events within patients, between different types of recurrent events and between recurrent events and the terminal event. Rogers et al. employed a model for which ‘the individual-specific frailties are assumed to affect the rate of heart failure hospitalisations and the hazard for cardiovascular death in the same way’ and this assumption could be challenged (Rogers et al. 2014b). However, a number of authors have described approaches that allow the frailty to differ between these processes (Mazroui et al. 2013). For the joint and shared frailty, both gamma (Mazroui et al. 2012; Rogers et al. 2014b) and log-normal (McGilchrist and Aisbett 1991; Belot et al. 2014) distributions are typically assumed; perhaps gamma is more common because it yields ‘relatively tractable likelihoods’ (Cook and Lawless 2007), although with the extra complexity of this model the marginal likelihood does not have a closed form and thus ‘using other distributions for the frailty, such as log normal... will not induce more difficulties’ (Rondeau et al. 2007). Suffice to say, the choice tends to be a mathematical one rather than a biological one (Wienke et al. 2003). In any case, authors have confirmed that results are robust to a misspecification of the frailty distribution (Pickles and Crouchley 1995; Mazroui et al. 2012; Belot et al. 2014). Note that if the processes are independent (or rather dependence is entirely captured by the covariates), then the distributional parameter will be close to zero and we may adopt the simpler model. Also, if the terminal event rate is low, we may opt for a more parsimonious model because, for example, the results can be more easily presented.

For the case of MTREs and a terminal event see Cook and Lawless (2007). An MTRE model like (3.1) which incorporates a terminal event as a joint frailty may be written as follows:

\begin{matrix} λ_{ij} (t; x_{i}) = λ_{0 j} (t) \exp (x_{i} β_{j} + u_{ij}) \\ d_{i} (t; x_{i}) = d_{0} (t) \exp (x_{i} α + γ_{1} u_{i 1} + γ_{2} u_{i 2}), \end{matrix}

(4.1)

where two event-types are assumed and

γ_{j}

determines the relationship between event-types and the terminal event. Several illustrations are noteworthy, for example, a multivariate lognormal random effects model incorporating a terminal event has been described by Mazroui et al. (2013) for breast cancer data. Zhu et al. (2010) considered the incidence of two infection types (bacterial, the most common, and fungal or viral) in acute myeloid leukaemia with relapse, transplant or death as the terminal event. Unlike Mazroui et al.’s model ‘the covariance effects on the terminal event were left arbitrary’ as was the dependence between recurrent events and death. Zeng et al. (2014) analyzed bleeding and transfusion events in patients with myelodysplastic syndrome; separate random effects for event-types by patient and patients were specified, the latter linking recurrent events and the informative censoring. Zhao et al. (2012) used a semi-parametric approach leaving the baseline intensity unspecified and treating the patient-level frailties as nuisance parameters and thus no distributional form need be assumed (covariates were event-type specific with no event-type frailties); they analyzed fever and other reaction rates (i.e., two event-types) after platelet transfusion among haematology/oncology patients. And Lin et al. (2017) analyzed cardiovascular events (coronary heart disease, stroke and heart failure) with all-cause mortality as the terminal event. To simplify matters, we could use a crude approach such as treating death as the last event (Rogers et al. 2014a; Rogers et al. 2014b), although when dealing with multiple types of events we must ask: which type of event? And equating death with an event will only make sense in certain contexts and if recurrent events are rare. Such analyses may be designated as supportive or sensitivity analyses.

Further, recent elaborations of the MTREs model have appeared in the literature such as the simultaneous modelling of longitudinal outcomes (Musoro et al. 2015), time-varying coefficients (Mazroui et al. 2016), ratios of intensity functions for types (Ning et al. 2017), interval censored data (Chen et al. 2005), gap time analysis (Lawless et al. 2001), a Bayesian approach (Lin et al. 2017) and handling of missing event-types (Chen and Cook 2009; Dewanji and Sengupta 2003; Lin et al. 2013). We may also require random effects for clustering levels such as patients within centres (Schaubel and Cai, 2005). Typically, such multi-level models lead to intractable likelihoods with high dimensional integration, that is, the frailty terms cannot be integrated out to obtain a closed form of the marginal likelihood (see (Duchateau and Janssen 2008) for multi-level frailty models). For these more complex frailty modelling approaches, we will want to know whether software is available or the statistician must code their analysis from scratch according to their specific requirements.

5 Software for MTREs

Developments described earlier, which have appeared since the publication of Abu-Libdeh's original random effects model, coincide with the increasing availability of more able statistical software required for the estimation of parameters, although ‘the current estimation methods for frailty proportional hazards models are not very satisfactory’ (Liu and Huang 2008) and ‘available software is limited’ (Rondeau and Gonzalez 2005). A majority of the papers cited earlier emerged in the recent literature, that is, within the last five years, although the topic of MTREs has seen renewed interest over the past decade: the statistical monograph of time-to-event data by Cook and Lawless (2007) (cited earlier) with a dedicated chapter on the topic is from 2007. Statistical software, naturally, lags behind. However, in a simple case such as the heart failure study where there are two event-types (ED and hospital), and with the baseline intensity specified (e.g., Weibull), we can use standard likelihood methods to obtain the parameter estimates and inference is straightforward, as follows.

The MTRE model (2.3) is a conditional model, events are independent given the random effects (since the random effects explain the dependence between events), and thus the conditional likelihood may be written out by hand as if the frailties were observed. The marginal likelihood is obtained by integrating, that is, averaging over the random effects (see Collett 2003, Section 11.3) to obtain a likelihood equation that contains the distributional parameters for the random effects and the parameters of interest (i.e., $β$ and the parameters of the Weibull). We then obtain maximum likelihood estimates of these parameters using optimization software (e.g., proc nlp in SAS), employing numerical methods such as Newton–Raphson (Lee and Wang 2013) and standard errors for the estimates are obtained from the inverse of the Hessian of the marginal likelihood. Abu-Libdeh use maximum likelihood estimation and provide the score vector and sample information matrix in their appendix, making their method easy to implement (Natarajan et al. 1994). We can obtain estimates of the random effects using Empirical Bayes which may be useful for calculating residuals to assess model fit and in the production of individual patient survival curves (with confidence intervals derived using the delta method; Collett 2003).

In cases that are not so straightforward, for example, when we have more event-types, or multivariate lognormal frailties, or incorporating a terminal event, the marginal likelihood contains high dimensional integrals, analytical integration is not possible and other more sophisticated methods (numerical integration or approximations) are required. The Bayesian approach seems a natural choice, a key feature of which is that parameters are unknown and described by a probability distribution (a prior distribution) like the random effects in the frailty model. Chen et al. (2005) postulate a model for interval censored data with multivariate lognormal frailties and implemented according to a Gibbs sampling algorithm. Likewise Lin et al. (also cited earlier) with three event-types and a terminal event used a Bayesian approach (Lin et al. 2017) (flat priors are assumed, e.g., $β$ parameters). In this case, WinBUGS software is favoured (Spiegelhalter et al. 2004) or proc mcmc in SAS. An alternative is to use quadrature (Naylor and Smith 1982) to approximate integrals and then the likelihood is maximized. (Liu and Huang 2008) applied Gaussian quadrature estimation using SAS proc nlmixed (Lu and Liu 2008) which we illustrate in the next section. A number of optimization algorithms are available in nlmixed for maximizing the likelihood, although the random statement is limited to normally distributed random effects. Mazroui et al. developed software that is capable of fitting the MTRE model, that is, an R package named frailtypack (Rondeau et al. 2012; Mazroui et al. 2013) with the penalized likelihood estimation for inference. Ideally, simulations should be exploited to compare the performance of the alternative estimation methods mentioned earlier. Liu et al. compared quadrature with Monte Carlo EM and penalized partial likelihood approaches and concluded that quadrature is preferable (Liu and Huang 2008). Further simulations are needed to clarify these results.

6 Illustration: Heart failure readmissions

In this section, we describe the thought process for planning and implementing an MTRE analysis. Consider the heart failure study described in the ‘Introduction’ with patients recruited when they presented at the ED. Heart failure related ED and hospital visits subsequent to discharge are the event-types of interest. ED visits are more subjective and we can expect they will show a higher incidence. We examined the predictive ability of atrial fibrillation at the index visit for subsequent ED and hospital visits. A simple analysis of such data using Cox regression of time-to-first readmission or death has been presented elsewhere (Gouda et al. 2016).

The model for repeat visits is as follows:

\begin{matrix} λ_{ij} (t; x_{i}) = λ_{0 j} (t) \exp (x_{i} β_{j} + u_{ij}) \end{matrix}

The random effects, or ‘frailties’, for event-types ( $u_{ij}$ : $i$ = patient, $j$ = hospital/ED visit) are bivariate normal, that is, N([0,0],[ $σ_{1}$ , ${cov}_{12}$ , $σ_{2}$ ]). (The multivariate normal distribution being a popular choice.) $β$ are the regression coefficients for the covariates (i.e., atrial fibrillation in this case) with hazard ratio $e^{β}$ . And we adopt a Weibull form for the baseline intensity as per Abu-Libdeh, that is, $λ_{0 j} (t) = ρ_{j} δ_{j} t^{δ_{j} - 1}$ . Estimates of the Weibull and $β$ parameters will be event-type specific. We could assume a common effect for event-types or test $H_{0} : β_{1} = β_{2}$ . We could also consider alternatives for the baseline intensity, for example, Liu and Huang (2008) prefer a piecewise hazard, although the coding is long-winded and thus more prone to errors.

For all-cause mortality, a joint frailty model is assumed:

\begin{matrix} d_{i} (t; x_{i}) = d_{0} (t) \exp (x_{i} α + γ_{1} u_{i 1} + γ_{2} u_{i 2}), \end{matrix}

where

γ_{j}

determines the relationship between Hospital/ED visit and death. The maximum likelihood estimates of the model parameters can be obtained using proc nlmixed in SAS. See Table 1 for a full summary of the results. Implementation is straight forward because we need only write out the integrand within nlmixed using if–then clauses for the different event-types with the distribution for the random effects specified in the random statement (SAS code for this example is available at the following link: https://papers-and-programs.blogspot.com). Initial estimates are important for convergence and could be obtained by fitting a simple shared frailty model, that is, without the event-type random effects. There is a lack of consensus regarding the use of adaptive quadrature and non-adaptive quadrature (this is the method of integral approximation and is specified in the proc statement in nlmixed). The statistician may opt to use adaptive quadrature by default (this is more computationally intensive) and resort to non-adaptive quadrature with limited q-points (e.g., 10), if difficulties with convergence arise (although it is important to understand the cause of non-convergence).

The correlation $(ρ)$ between event-types is calculable from the covariance matrix as ${cov}_{12} / \sqrt{(σ_{1} \times σ_{2})}$ . Using the estimate statement in nlmixed we obtain an estimate of 0.88 with 95% confidence interval 0.80–0.97. Thus, patients with a high risk of ED visits tend to also show a high risk of hospital readmissions. We can also examine the association between mortality and each event-type via $γ_{1}$ and $γ_{2}$ . The estimates are 1.14 (95% CI 0.26, 2.01) for hospital readmissions and $-$ 0.07 (95% CI $-$ 0.65, 0.51) for ED visits (Table 1) indicating that patients with a high risk of hospital readmission have more severe disease, leading to death. This justifies handling ED and hospital visits as different event-types and is an important finding consistent with previous research. For example, (Collins et al. 2015) state that ‘compared with outpatients with a similar degree of cardiac dysfunction, hospitalised patients have significantly worse outcomes’. Insight into the relatedness of event-types is a unique and important feature of the MTRE approach.

Using additional estimate statements we may test the null hypothesis of no association between mortality and the event types, for example, $γ_{1} = 0$ . Estimates of the random effects (or individual patient frailties) are produced using the Empirical Bayes method which is specified in the random statement. These estimates enable the statistician to calculate Martingale residuals and evaluate the model fit (we would hope to see the residuals for each event-type centred around zero). They may also allow us to produce prediction curves for patient survival; nlmixed has a predict statement which can be used for the purpose. These additional insights are an advantage of the MTRE analysis. (See SAS code at the weblink mentioned earlier for the full code including these estimate statements.)

Table 1:

Hazard ratio estimates and $p$ -values from the MTRE analysis

Analysis method	Hospital readmission	ED visit	Mortality
MTRE model	1.24 (0.78, 1.70) $p = 0.2590$	1.47 (1.00, 1.93) $p = 0.0183$
Mortality, $γ$	1.14 (0.26, 2.01) $p = 0.0112$	$-$ 0.07 ( $-$ 0.65, 0.51) $p = 0.8144$	1.23 (0.78, 1.68) $p = 0.2731$

Note: The correlation for bivariate frailties is 0.88 (0.80–0.97). The model includes A. Fib. covariate only. See model specification mentioned earlier.

The heart failure data are simple in that we have only two event-types. We may wish to extend our analysis of readmissions, for example, by including other primary diagnoses rather than restricting to heart failure (since, e.g., acute myocardial infarction increases the risk of heart failure (Velagaleti et al. 2008) and heart failure may comprise less than half of the total visits experienced (Babayan et al. 2003)). Alternatively, readmissions may be extended to three classifications: ED, ED to hospital (transferred to hospital on the same day or the following day) and hospital. In cardiovascular research, in general, we may encounter more than two event-types, for example, transient ischemic attacks or stroke may be classified by location. Likewise, Abu-Libdeh contemplated incorporating site of lesion into their analysis of BCCs and SCCs. We would, therefore, need to extend the model for more than two event-types, and perhaps considerably more. Data could become thin for these additional types and convergence may become difficult to attain.

7 Conclusion

The simple modelling approach, that is, a Cox regression of time-to-first event, remains a popular choice in clinical trials of chronic diseases. However, there is a recent push to include all events in the analysis and this has prompted a comparison of modelling strategies to identify which among them is optimal (Rogers et al. 2014b). Models for MTREs have been absent from this discussion even though recurrent events often imply some categorization of events into types. There may be extraneous or practical reasons that dissuade statisticians from allocating an MTRE model as the primary analysis. For example, despite the assistance offered by recent publications and statistical programmers, the statistician is required to invest more time in the early stages of the analysis. Also, a time-to-first analysis provides a simple basis for a power calculation; an MTRE model on the other hand requires data simulations. The time-to-first analysis, by amalgamating outcomes, obviates the issue of multiple testing and the adjustment of the nominal significance level (although an MTRE model could produce a weighted average across outcomes to estimate the net effect too). There may also be some resistance to unfamiliar, esoteric results when hoping to reach a broad audience.

Nevertheless, with the increasing relevance of recurrent events and MTREs in various clinical contexts, and with implementation using standard software now a straightforward matter, it will become increasingly difficult to neglect MTRE modelling. By providing a more complete characterization of the data and an increase in statistical efficiency, MTRE modelling could alter the conclusions drawn from clinical trial or registry data. If it is difficult to pre-specify a distribution for the random effects for event-types, or convergence is a concern, then the MTRE model could provide a supportive or sensitivity analysis, bearing in mind the trade-off between statistical rigour and cogency of results. Ultimately, of course, it will depend on the particular circumstances (e.g., the event rate, extent of follow-up, reliability of data collection, etc.) and the purpose of the analysis. A simulation study evaluating how MTRE models perform under differing circumstances would certainly be informative and aid this decision.

In the introduction, we quoted a Lancet paper where nonmelanoma skin cancers (SCC and BCC) were analyzed separately, yet we see Abu-Libdeh analyzing nonmelanoma skin cancer incidence using an MTRE model. This inconsistency implies more than a mere preference for one statistical approach over another. It implies a contradictory understanding of the clinical condition. Today, with increasing attention given to the associations among event-types, the justification for analyzing separately may no longer be readily accepted. A 2002 paper by the eminent cardiologist Dr Salam Yusuf was titled ‘Choice of Clinical Outcomes in Randomized Trials of Heart Failure Therapies: Disease-specific or Overall Outcomes?’ (Yusuf and Negassa 2002). One concern being that when all hospitalizations are combined it can mask a treatment effect that would be seen on disease-specific outcomes. Although it seems there is a third option we are neglecting, that is, the MTRE model.

References

Aalen

Borgan

Gjessing

(2008) Survival and Event History Analysis . New York: Springer.

Abu-Libdeh

Turnbull

Clark

(1990) Analysis of multi-type recurrent events in longitudinal studies; application to a skin cancer prevention trial. Biometrics , 46, 1017–34.

Andersen

Borgan

Gill

Keiding

(1993) Statistical Models Based on Counting Processes . New York: Springer.

Andersen

Gill

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

Anker

McMurray

(2012) Time to move on from ‘time-to-first’: Should all events be included in the analysis of clinical trials? European Heart Journal , 33, 2764–5.

Babayan

McNamara

Nagajothi

Kasper

Armenian

Powe

Baughman

Lima

(2003) Predi- ctors of cause-specific hospital readmission in patients with heart failure. Clinical Cardiology , 26, 411–8.

Belot

Rondeau

Remontet

Giorgi

(2014) A joint frailty model to estimate the recurrence process and the disease-specific mortality process without needing the cause of death. Statistics in Medicine , 33, 3147–66.

Castaneda

Gerritse

(2010) Appraisal of several methods to model time to multiple events per subject: Modelling time to hospitalizations and death. Revista Colombiana de Estadistica , 33, 43–61.

Chen

Cook

(2009) The analysis of multivariate recurrent events with partially missing event types. Lifetime Data Analysis , 15, 41–58.

10.

Chen

Cook

Lawless

Zhan

(2005) Statistical methods for multivariate interval-censored recurrent events. Statistics in Medicine , 24, 671–91.

11.

Chen

Wang

Cai

Shankar

(2012) Semiparametric additive marginal regre- ssion models for multiple type recurrent events. Lifetime Data Analysis , 18, 504–27.

12.

Chun

Wijeysundera

Austin

Wang

Levy

Lee

(2012) Lifetime analysis of hospitalizations and survival of patients newly admitted with heart failure. Circulation: Heart Failure , 5, 414–21.

13.

Clayton

(1978) A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika , 65, 141–51.

14.

Collett

(2003) Modelling Survival Data in Medical Research, second edition. Florida, US: Chapman and Hall/CRC Texts in Statistical Science. Chapman and Hall/CRC.

15.

Collins

Storrow

Albert

Butler

Ezekowitz

Felker

Fermann

Fonarow

Givertz

Hiestand

Hollander

Lanfear

Levy

Pang

Peacock

Sawyer

Teerlink

Lenihan

(2015) Early management of patients with acute heart failure: State of the art and future directions. A consensus document from the society for academic emergency medicine/heart failure society of america acute heart failure working group. Journal of Cardiac Failure , 21, 27–43.

16.

Cook

Lawless

(2007) The Statistical Analysis of Recurrent Events. (Statistics for Biology and Health), first edition. New York: Springer.

17.

Cox

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

18.

Davidson

Liebl

Christiansen

Fulcher

Ligthelm

Brown

Gylvin

Kawamori

(2009) Risk for nocturnal hypoglycemia with biphasic insulin aspart 30 compared with biphasic human insulin 30 in adults with type 2 diabetes mellitus: A meta-analysis. Clinical Therapeutics , 31, 1641–51.

19.

Dewanji

Sengupta

(2003) Estimation of competing risks with general missing pattern in failure types. Biometrics , 59, 1063–70.

20.

Duchateau

Janssen

(2008) The Frailty Model, first edition. (Statistics for Biology and Health). New York: Springer.

21.

Gouda

Brown

Rowe

McAlister

Ezekowitz

(2016) Insights into the importance of the electrocardiogram in patients with acute heart failure. European Journal of Heart Failure , 18, 1032–40.

22.

Green

Williams

Nale

Hart

Leslie

Parsons

Marks

Gaffney

Battistutta

Frost

Lang

Russell

(1999) Daily sunscreen application and betacarotene supplementation in prevention of basal-cell and squamous-cell carcinomas of the skin: A randomised controlled trial. The Lancet , 354, 723–29.

23.

Greenberg

Yaroshinsky

Zsebo

Butler

Felker

Voors

Rudy

Wagner

Hajjar

(2014) Design of a phase 2b trial of intracoronary administration of aav1/serca2a in patients with advanced heart failure: the cupid 2 trial (calcium up-regulation by percutaneous administra- tion of gene therapy in cardiac disease phase 2b). JACC: Heart Failure , 2, 84–92.

24.

Hougaard

(2000) Analysis of Multivariate Survival Data . New York: Springer.

25.

Lawless

(1987) Regression methods for poisson process data. Journal of the American Statistical Association , 82, 808–15.

26.

Lawless

Wigg

Tuli

Drake

Lamberti-Pasculli

(2001) Analysis of rep- eated failures or durations, with application to shunt failures for patients with paediatric hydrocephalus. Journal of the Royal Statistical Society, Series C (Applied Statistics) , 50, 449–65.

27.

Lee

Austin

Stukel

Alter

Chong

Parker

(2009) ‘Dose- dependent’ impact of recurrent cardiac events on mortality in patients with heart failure. The American Journal of Medicine , 122, 162–69 e1.

28.

Lee

Wang

(2013) Statistical Methods for Survival Data Analysis, fourth edition. New Jersey, US: Wiley.

29.

Liu

(2015) A generalized semiparametric mixed model for analysis of multivariate health care utilization data. Statistical Methods in Medical Research . doi 10.1177/0962280215615159.

30.

Lin

Cai

Fine

Lai

(2013) Nonparametric estimation of the mean function for recurrent event data with missing event category. Biometrika , 100, 727–40.

31.

Lin

Luo

Chen

Davis

(2017) Bayesian analysis of multi-type recurrent events and dependent termination with nonparametric covariate functions. Statistical Methods in Medical Research , 26, 2869–84.

32.

Lindeboom

Van Den

Berg G

(1994) Heterogeneity in models for bivariate survival: The importance of the mixing distribution. Journal of the Royal Statistics, Society B , 56, 4960.

33.

Liu

Huang

(2008) The use of gaussian quadrature for estimation in frailty proportional hazards models. Statistics in Medicine , 27, 2665–83.

34.

Locatelli

Lichtenstein

Yashin

(2004) The heritability of breast cancer: A bayesian correlated frailty model applied to swedish twins data. Twin Research , 7, 182–91.

35.

Liu

(2008) Analysis of correlated recurrent and terminal events data in sas. NESUG . URL http://www.lexjansen.com/nesug/nesug08/sa/sa16.pdf

36.

Mazroui

Mathoulin-Pelissier

Macgrogan

Brouste

Rondeau

(2013) Multivariate frailty models for two types of recurrent events with a dependent terminal event: application to breast cancer data. Biometrical Journal , 55, 866–84.

37.

Mazroui

Mathoulin-Pelissier

Soubeyran

Rondeau

(2012) General joint frailty model for recurrent event data with a dependent terminal event: Application to follicular lymphoma data. Statistics in Medicine , 31, 1162–76.

38.

Mazroui

Mauguen

Mathoulin-Pelissier

MacGrogan

Brouste

Rondeau

(2016) Time-varying coefficients in a multivariate frailty model: Application to breast cancer recurrences of several types and death. Lifetime Data Analysis , 22, 191–215.

39.

McGilchrist

Aisbett

(1991) Regression with frailty in survival analysis. Biometrics , 47, 461–6.

40.

Musoro

Geskus

Zwinderman

(2015) A joint model for repeated events of different types and multiple longitudinal outcomes with application to a follow-up study of patients after kidney transplant. Biometrical Journal , 57, 185–200.

41.

Natarajan

Turnbull

Slate

Wells

Clark

Abu-Libdeh

(1994) A computer program for the statistical analysis of repeated event data using a mixed effects regression model. Computer Methods and Programs in Biomedicine , 42, 283–94.

42.

Naylor

Smith

AFM

(1982) Applications of a method for the efficient computation of posterior distribution. Journal of the Royal Statistical Society. Series C (Applied Statistics) , 31, 214–25.

43.

Ning

Rahbar

Choi

Piao

Hong

Del Junco

Rahbar

Fox

Holcomb

Wang

(2017) Estimating the ratio of multivariate recurrent event rates with application to a blood transfusion study. Statistical Methods in Medical Research , 26, 1969–81.

44.

O'Brien

Subherwal

Roe

Holmes

Thomas

Alexander

Wang

Peterson

(2014) Do patients treated at academic hospitals have better longitudinal outcomes after admission for non-st-elevation myocardial infarction? American Heart Journal , 167, 762–9.

45.

Pandeya

Purdie

Green

Williams

(2005) Repeated occurrence of basal cell carcinoma of the skin and multifailure survival analysis: Follow-up data from the nambour skin cancer prevention trial. American Journal of Epidemiology , 161, 748–54.

46.

Pickles

Crouchley

(1995) A comparison of frailty models for multivariate survival data. Statistics in Medicine , 14, 1447–61.

47.

Pocock

Ariti

Collier

Wang

(2012) The win ratio: A new approach to the analysis of composite endpoints in clinical trials based on clinical priorities. European Heart Journal , 33, 176–82.

48.

Rogers

Jhund

Perez

Bohm

Cleland

Gullestad

Kjekshus

van Veldhuisen

Wikstrand

Wedel

McMurray

Pocock

(2014a) Effect of rosuvastatin on repeat heart failure hospitalizations: The corona trial (controlled rosuvastatin multinational trial in heart failure). JACC: Heart Failure , 2, 289–97.

49.

Rogers

Pocock

McMurray

Granger

Michelson

Ostergren

Pfeffer

Solomon

Swedberg

Yusuf

(2014b) Analysing recurrent hospitalizations in heart failure: A review of statistical methodology, with application to charm-preserved. European Journal of Heart Failure , 16, 33–40.

50.

Rondeau

Gonzalez

(2005) Frailtypack: A computer program for the analysis of correlated failure time data using penalized likelihood estimation. Computer Methods and Programs in Biomedicine , 80, 154–64.

51.

Rondeau

Mathoulin-Pelissier

Jacqmin- Gadda

Brouste

Soubeyran

(2007) Joint frailty models for recurring events and death using maximum penalized likelihood estimation: Application on cancer events. Biostatistics , 8, 708–21.

52.

Rondeau

Mazroui

Gonzalez

(2012) Frailtypack: An r package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametrical estimation. Journal of Statistical Software , 47, 1–28.

53.

Schaubel

Cai

(2005) Analysis of clustered recurrent event data with application to hospitalization rates among renal failure patients. Biostatistics , 6, 404–19.

54.

Spiegelhalter

Abrams

Myles

(2004) Bayesian Approaches to Clinical Trials and Health-care Evaluation, first edition. Chichester, England: John Wiley and Sons.

55.

Velagaleti

Pencina

Murabito

Wang

Parikh

D'Agostino

Levy

Kannel

Vasan

(2008) Long-term trends in the incidence of heart failure after myocardial infarction. Circulation , 118, 2057–62.

56.

Villegas

Julia

Ocana

(2013) Empirical study of correlated survival times for recurrent events with proportional hazards margins and the effect of correlation and censoring. BMC Medical Research Methodology , 13, 95.

57.

Wienke

(2010) Frailty models in survival analysis, first edition . Florida, US: Chapman and Hall/CRC.

58.

Wienke

Arbeev

Locatelli

Yashin

(2003) A simulation study of different correlated frailty models and estimation strategies. URL http://www.demogr.mpg.de/en/projects_publications/publications_1904/mpidr_working_papers/all.htm

59.

Yusuf

Negassa

(2002) Choice of clinical outcomes in randomized trials of heart failure therapies: Disease-specific or overall outcomes? American Heart Journal , 143, 22–8.

60.

Zeng

Ibrahim

Chen

Jia

(2014) Multivariate recurrent events in the presence of multivariate informative censoring with applications to bleeding and transfusion events in myelodysplastic syndrome. Journal of Biopharmaceutical Statistics , 24, 429–42.

61.

Zhao

Liu

(2012) Analysis of multivariate recurrent event data with time-dependent covariates and informative censoring. Biometrical Journal , 54, 585–99.

62.

Zhu

Sun

Tong

Srivastava

(2010) Regression analysis of multivariate recurrent event data with a dependent terminal event. Lifetime Data Analysis , 16, 478–90.

Frailty modelling for multitype recurrent events in clinical trials

Abstract

Keywords

1 Introduction and motivating examples

Figure 1:

Example heart failure readmissions data

6 Illustration: Heart failure readmissions

Table 1:

Hazard ratio estimates and p -values from the MTRE analysis

References

Hazard ratio estimates and $p$ -values from the MTRE analysis