Joint frailty models for recurrent and terminal events

In this post we’re going to take a look at joint frailty models, and how to fit them with our merlin command. Importantly, we’ll also discuss how to interpret the results.

Joint frailty models

An area of intense research in recent years is that of joint frailty models, which has become the commonly used name for a joint model for a recurrent event and a terminal event. We’re going to take a look at the most popular approach (Liu et al., 2004), and how to implement it in Stata.

In essence, we have a survival model for the recurrent event process, a survival model for the terminal event process, and we link them through a shared random effect. In other words, we have a random effect which accounts for the correlation between recurrent events, which is then included in the linear predictor for the terminal event model, with an accompanying coefficient to be estimated, which directly quantifies the strength of the association between the two processes. The nice property of this formulation is that if there is no association, then the models reduce to their separate versions. Let’s formalise it.

We have for the recurrent event model, the hazard function for the ith patient and the jth event,

h_{i j} (t) = h_{0} (t) \exp (X_{1 i j} β_{1} + b_{i})

where $h_{0} (t)$ is the baseline hazard rate, $X_{1 i j}$ is a vector of baseline covariates with associated log hazard ratios, $β_{1}$ , and finally a random intercept, $b_{i} ~ N (0, σ^{2})$ . So far, this is a standard frailty survival model (we’re going to use frailty and random effect interchangeably in this post), where each patient’s events share the same unobserved effect $b_{i}$ , which accounts for the correlation between events occurring in the same patient.

To bring in the terminal event process, we define the mortality rate for the ith patient,

λ_{i} (u) = λ_{0} (u) \exp (X_{2 i} β_{2} + α b_{i})

where $λ_{0} (u)$ is the baseline mortality rate, $X_{2 i}$ is a vector of baseline covariates with associated log hazard ratios, $β_{2}$ , and $α$ directly quantifies the association between the recurrent and terminal event processes. To be explicit, $\exp (α)$ represents the hazard ratio for a one-unit increase in $b_{i}$ … which is not that simple. The important way to look at it is if $α > 0$ , then those with a higher frailty (i.e. higher underlying recurrent event rate) have an increased mortality rate. The other way around, if $α < 0$ then those with a higher frailty have a reduced mortality rate.

Example

We illustrate these models with simulated data, loosely modelled on the readmission dataset that comes with the extensive frailtypack in R (Król et al., 2017), developed by Virginie Rondeau and her group. That dataset has information on re-hospitalisation times after surgery in patients diagnosed with colorectal cancer. Its covariates of interest include gender, Dukes’ tumour stage and comorbidity Charlson index, but to keep things as simple as possible, we’re just going to simulate gender, as a male binary indicator variable, and include it in our model below. We simulate 400 patients, about 60% of them male, and give each a frailty, $b_{i}$ , drawn from a normal distribution with $σ = 1$ .

Stata

. clear

. set seed 4219

. set obs 400                             // patients
Number of observations (_N) was 0, now 400.

. gen id = _n

. gen male = runiform() < 0.6

. gen b = rnormal(0, 1)                   // frailty, sd 1

Time is measured in years since surgery. First, we simulate each patient’s time to death with avalon standard, from a Weibull baseline hazard, $λ_{0} (u) = λ γ u^{γ - 1}$ , with $λ = 0.03$ and $γ = 1.1$ , a log hazard ratio of 0.2 for being male, and the frailty entering with $α = 1.5$ (covariates() takes each variable with its coefficient). Follow-up stops at five years.

Stata

. avalon standard stime death,            ///
>         distribution(weibull)           ///
>         lambda(0.03) gamma(1.1)         /// baseline
>         covariates(male 0.2 b 1.5)      /// alpha = 1.5
>         maxtime(5)                      //  five years' follow-up

Then the re-hospitalisations. On the clock-reset timescale, the gaps between them are, given the frailty, independent draws from a single distribution, here a Weibull with $λ = 0.25$ and $γ = 0.7$ , a log hazard ratio of 0.4 for being male, and the frailty entering with a coefficient of 1. avalon standard’s recurrence() option draws them one after another, restarting the clock at each re-hospitalisation, until the time we give it in maxtime(), here each patient’s own time of death or censoring, stime. It writes a pair of variables for each gap, stop1 and event1, stop2 and event2, and so on: stop# is the time since surgery at which the gap ended, and event# says whether it ended in a re-hospitalisation (1) or in death or censoring (0). maxevents(50) caps the number of gaps, far above the 17 that the busiest patient here needs.

Stata

. avalon standard stop event,             /// one stop# and event# per gap
>         distribution(weibull)           ///
>         lambda(0.25) gamma(0.7)         /// baseline, on the gap timescale
>         covariates(male 0.4 b 1)        /// frailty enters with coefficient 1
>         recurrence(maxevents(50))       /// the clock resets at each event
>         maxtime(stime)                  //  until death or censoring

reshape long turns the pairs into rows, one per gap, and we drop the empty ones after each patient’s last gap, checking that every patient’s last gap is censored, at death or at five years, whichever comes first. Each gap starts where the previous one stopped, so its length, the clock-reset time, is the difference between the two. The time to death and its indicator go on each patient’s last row only. Each patient therefore has as many rows as they had re-hospitalisations, plus one.

Stata

. reshape long stop event, i(id) j(j)
(j = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17)

Data                               Wide   ->   Long
-----------------------------------------------------------------------------
Number of observations              400   ->   6,800
Number of variables                  39   ->   8
j variable (17 values)                    ->   j
xij variables:
                 stop1 stop2 ... stop17   ->   stop
              event1 event2 ... event17   ->   event
-----------------------------------------------------------------------------

. drop if missing(event)                  // gaps after death or censoring
(5,867 observations deleted)

. bys id (j) : assert event[_N] == 0      // everyone's last gap is censored

. bys id (j) : gen double start = cond(_n == 1, 0, stop[_n-1])

. gen double time = stop - start          // clock reset: time since last event

. bys id (j) : replace stime = . if _n < _N
(533 real changes made, 533 to missing)

. bys id (j) : replace death = . if _n < _N
(533 real changes made, 533 to missing)

. drop b j

Let’s take a look at our dataset:

Stata

. list id time event stime death male if inlist(id,3,54)

     +---------------------------------------------------+
     | id        time   event       stime   death   male |
     |---------------------------------------------------|
  5. |  3   .29014997       1           .       .      0 |
  6. |  3   1.1419176       1           .       .      0 |
  7. |  3   3.5679325       0           5       0      0 |
134. | 54   1.1182214       1           .       .      1 |
135. | 54   .12372617       1           .       .      1 |
     |---------------------------------------------------|
136. | 54   1.9156641       0   3.1576117       1      1 |
     +---------------------------------------------------+

The gap times, since surgery or the previous re-hospitalisation, are stored in time, with corresponding event indicator event. This is in clock-reset formulation, i.e. each time a patient is re-hospitalised the clock is reset to zero, so we will be fitting a semi-Markov model in this example. Our overall survival time is stored in stime, with corresponding event indicator stored in death. Patient 3 was re-hospitalised twice and was still alive at five years; patient 54 was also re-hospitalised twice, the second time only 0.12 years after the first, and died 3.16 years after surgery.

We can fit such a model with merlin, adjusting for male,

Stata

. merlin (time                             /// rehosp. times
>              male                        /// male
>              M1[id]@1                    /// random intercept
>              , family(weibull,           /// distribution
>                         failure(event))) ///
>        (stime                            /// survival time
>              male                        /// male
>              M1[id]                      /// random effect & association
>              , family(weibull,           /// distribution
>                         failure(death))) //

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood = -1314.0523  (not concave)
Iteration 1:  Log likelihood =  -1291.855
Iteration 2:  Log likelihood = -1274.1574  (not concave)
Iteration 3:  Log likelihood = -1256.4583
Iteration 4:  Log likelihood = -1254.3132
Iteration 5:  Log likelihood = -1254.3009
Iteration 6:  Log likelihood = -1254.3007
Iteration 7:  Log likelihood = -1254.3007

Mixed effects regression model                             Number of obs = 933
Log likelihood = -1254.3007
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
time:        |
        male |   .2882978   .1411605     2.04   0.041     .0116283    .5649673
      M1[id] |          1          .        .       .            .           .
       _cons |  -1.329876   .1254899   -10.60   0.000    -1.575832   -1.083921
  log(gamma) |  -.4460228    .039422   -11.31   0.000    -.5232884   -.3687571
-------------+----------------------------------------------------------------
stime:       |
        male |   .4487432   .2461113     1.82   0.068    -.0336262    .9311125
      M1[id] |   1.317355   .2267665     5.81   0.000     .8729007    1.761809
       _cons |  -3.339901   .2978027   -11.22   0.000    -3.923584   -2.756218
  log(gamma) |  -.0430546   .0943663    -0.46   0.648    -.2280091    .1418998
-------------+----------------------------------------------------------------
id:          |
      sd(M1) |   .9331898   .0836885                      .7827701    1.112515
------------------------------------------------------------------------------

The first submodel is our model for the time to re-hospitalisation, where we’ve assumed a Weibull baseline, and a constant effect of male. The term M1[id] is the syntax required to specify a normally distributed random effect, with mean 0, in merlin. Its name must be M followed by a positive integer, which leaves room to add more (M2, M3, and so on). In square brackets we have to define the level at which the random effect applies – in our case we want to specify a random effect at the id level. Finally, the @1 notation constrains the coefficient on the random effect to be 1. If we didn’t specify this, it would by default have a coefficient which would be estimated. Keep reading.

The second submodel specification is our model for overall survival. We again adjust for male, and this time include the same random effect, M1, but without any constraint on its coefficient, so merlin estimates one. This corresponds to $α$ in the above formulation. Our syntax (based on Stata’s gsem command) provides a highly convenient way of linking random effects between outcome models. Let’s get to the results.

Our results show a highly positive estimate of association of 1.317 (95% CI: 0.873, 1.762) showing substantial association between the two event processes, i.e. a higher rate of re-hospitalisation also indicates a higher mortality rate. That is the direction we built into the data, and the one you’d expect with real re-hospitalisation and mortality. As a side note, we must remember to interpret our covariate effects conditional on the frailty.

Because we simulated the data, we can check the estimates against the truth. We set $α = 1.5$ , inside that interval, and $σ = 1$ , estimated at 0.933. The log hazard ratios for being male are estimated at 0.288 for re-hospitalisation and 0.449 for death (we simulated 0.4 and 0.2), each within about one standard error. The Weibull shape parameters are reported as log(gamma), and their estimates are −0.446 and −0.043, against the logs of the 0.7 and 1.1 we simulated, −0.357 and 0.095. The first is 2.3 standard errors away, which among eight estimates is not unusual.

Extensions

So what should we be thinking about next? Well, adjusting for more covariates, assessing proportional hazards, non-linear covariate effects, assessing the appropriateness of the Weibull baselines – all of these things can be incorporated very simply with merlin. Let’s fit Royston-Parmar models instead,

Stata

. merlin (time                             /// rehosp. times
>              male                        /// male
>              M1[id]@1                    /// random intercept
>              , family(rp, df(3)          /// distribution
>                         failure(event))) ///
>        (stime                            /// survival time
>              male                        /// male
>              M1[id]                      /// random effect & association
>              , family(rp, df(2)          /// distribution
>                         failure(death))) //
variables created: _rcs1_1 to _rcs1_3
variables created: _rcs2_1 to _rcs2_2

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood =  -1307.541  (not concave)
Iteration 1:  Log likelihood = -1278.7967
Iteration 2:  Log likelihood = -1257.1738
Iteration 3:  Log likelihood = -1253.9669
Iteration 4:  Log likelihood = -1253.4122
Iteration 5:  Log likelihood = -1253.3989
Iteration 6:  Log likelihood = -1253.3991

Mixed effects regression model                             Number of obs = 933
Log likelihood = -1253.3991
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
time:        |
        male |   .2848374   .1394276     2.04   0.041     .0115644    .5581104
      M1[id] |          1          .        .       .            .           .
       _cons |  -1.583851   .1290997   -12.27   0.000    -1.836882    -1.33082
-------------+----------------------------------------------------------------
stime:       |
        male |   .4572119   .2490982     1.84   0.066    -.0310116    .9454353
      M1[id] |   1.383109   .2463798     5.61   0.000     .9002135    1.866004
       _cons |  -2.167166   .2506775    -8.65   0.000    -2.658485   -1.675847
-------------+----------------------------------------------------------------
id:          |
      sd(M1) |   .9142385   .0866446                      .7592574    1.100855
------------------------------------------------------------------------------
    Warning: Baseline spline coefficients not shown - use ml display

The Royston-Parmar models tell the same story: $α$ is estimated at 1.383 and $σ$ at 0.914. The Weibull was the true model for both processes here, and the splines contain it as a special case, so they agree with it. Note that each submodel can be as different or similar as you like. Importantly, we now know how to formulate the crucial association structure in a joint frailty model, using a random intercept/frailty effect, so all of these extensions can be built on to our model.

In the above we are assuming a clock-reset timescale for the recurrent event process, which is how we simulated it. The clock-forward timescale, where time is measured from surgery and each interval between re-hospitalisations enters the model at the time the previous one ended (delayed entry), is supported from merlin 3.0.0, which makes the delayed-entry correction conditional on the frailty, as this model needs. Our data keep each interval’s start and end, on the time-since-surgery scale, in start and stop, so the clock-forward model would be merlin (stop male M1[id]@1, family(weibull, failure(event) ltruncated(start))) (stime male M1[id], family(weibull, failure(death))). For these data, simulated on the clock-reset timescale, the clock-reset model is the right one.

Joint frailty models can be delicate to fit. If the log-likelihood falls from one iteration to the next and merlin ends with “convergence not achieved”, with $α$ drifting upwards, try more quadrature points, for example intpoints(15). We have more to come on joint frailty models, so do subscribe or follow us on LinkedIn.

MethodsJoint longitudinal–survival models Recurrent events

Joint frailty models for recurrent and terminal events