An introduction to joint modelling of longitudinal and survival data

This post gives a gentle introduction to the joint longitudinal-survival model framework, and covers how to estimate them using our merlin command in Stata.

A joint model consists of a continuous, repeatedly measured (longitudinal) outcome, and a time-to-event, with the two models linked by random effects, or functions of them. Let’s formally define everything we need.

For the ith patient, we have the observed longitudinal outcome measured at time t, \(y_{i}(t)\),

$$y_{i}(t) = m_{i}(t) + \epsilon_{i}(t)$$

where

$$m_{i}(t) = X_{1i}(t)\beta_{1} + Z_{i}(t)b_{i}$$

and \(\epsilon_{i}(t)\) is our normally distributed residual variability (measurement error). We call \(m_{i}(t)\) our trajectory function, representing the true underlying value of the continuous outcome at time t, which is a function of fixed and random effects with associated design matrices, \(X_{1i}(t)\) and \(Z_{i}(t)\), and coefficients and random effects, \(\beta_{1}\) and \(b_{i}\), respectively. We assume normally distributed random effects,

$$b_{i} \sim N(0, \Sigma).$$

Our survival model can be defined in terms of the hazard function,

$$h_{i}(t) = h_{0}(t) \exp (X_{2i}\beta_{2} + \alpha m_{i}(t))$$

where \(h_{0}(t)\) is the baseline hazard function, \(X_{2i}\) is a vector of baseline covariates with associated log hazard ratios \(\beta_{2}\). We can then link the current value of the biomarker directly to survival, where \(\alpha\) is a log hazard ratio for a one-unit increase in the biomarker, at time t. This is just one way of linking the models…keep reading.

Let’s assume a random intercept and random linear trend for the biomarker trajectory and a further fixed effect of time squared to allow for more flexibility, i.e.

$$m_{i}(t) = (\beta_{10} + b_{0i}) + (\beta_{11}+b_{1i})t + \beta_{12}t^{2}$$

Our longitudinal model can be as simple or complex as we need, both in terms of fixed and random effects.

Joint model with the current value association structure

We’ll use the commonly utilised Primary Biliary Cirrhosis (PBC) example. In this post we will model the log of serum bilirubin (logb) over time, and investigate its association with survival, with survival times stored in stime. We load the dataset,

. use https://www.mjcrowther.co.uk/data/jm_example.dta, clear
(Example dataset for joint modelling, Michael Crowther 2011)

and let’s take a look:

. list id logb time trt stime died if inlist(id,1,2), sepby(id)
      +-------------------------------------------------------+
      | id        logb      time         trt     stime   died |
      |-------------------------------------------------------|
   1. |  1    2.674149         0   D-penicil   1.09517      1 |
   2. |  1    3.058707   .525682   D-penicil         .      . |
      |-------------------------------------------------------|
   3. |  2    .0953102         0   D-penicil   14.1523      0 |
   4. |  2   -.2231435   .498302   D-penicil         .      . |
   5. |  2           0   .999343   D-penicil         .      . |
   6. |  2    .6418539   2.10273   D-penicil         .      . |
   7. |  2    .9555114   4.90089   D-penicil         .      . |
   8. |  2    1.280934   5.88928   D-penicil         .      . |
   9. |  2    1.435084   6.88588   D-penicil         .      . |
  10. |  2    1.280934    7.8907   D-penicil         .      . |
  11. |  2    1.526056   8.83255   D-penicil         .      . |
      +-------------------------------------------------------+

Hopefully the dataset is fairly self-explanatory. We can fit the above joint model very simply using merlin, as follows

. merlin (logb                                      /// long. outcome
>                 fp(time, powers(1 2))             /// fractional polynomial
>                 time#M2[id]@1                     /// random slope on time
>                 M1[id]@1                          /// random intercept
>                 , family(gaussian)                /// distribution
>                 timevar(time))                    /// timevar
>        (stime                                     /// survival time
>                 trt                               /// baseline trt effect
>                 EV[logb]                          /// current value of 
>                 , family(weibull, failure(died))  /// distribution
>                 timevar(stime))                   /// timevar
>        , covariance(unstructured)                 //  VCV of random effects
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2
Fitting fixed effects model:
Fitting full model:
Iteration 0:   log likelihood = -3147.6761  (not concave)
Iteration 1:   log likelihood = -2295.5481  (not concave)
Iteration 2:   log likelihood = -2235.0456  (not concave)
Iteration 3:   log likelihood = -2145.8808  
Iteration 4:   log likelihood = -1946.6994  
Iteration 5:   log likelihood = -1914.2293  
Iteration 6:   log likelihood = -1911.9839  
Iteration 7:   log likelihood = -1911.9744  
Iteration 8:   log likelihood = -1911.9744  
Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -1911.9744
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |            
      fp():1 |   .1599175   .0152541    10.48   0.000     .1300201     .189815
      fp():2 |   .0040188    .001052     3.82   0.000     .0019569    .0060807
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5141206   .0580541     8.86   0.000     .4003366    .6279045
  sd(resid.) |   .3446222   .0066275                      .3318743    .3578598
-------------+----------------------------------------------------------------
stime:       |            
         trt |   .0428075   .1794715     0.24   0.811    -.3089501    .3945651
        EV[] |   1.246952   .0930087    13.41   0.000     1.064659    1.429246
       _cons |  -4.403106    .271847   -16.20   0.000    -4.935917   -3.870296
  log(gamma) |   .0054788   .0825863     0.07   0.947    -.1563873     .167345
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9940306   .0423725                      .9143567    1.080647
      sd(M2) |   .1891964   .0129363                      .1654672    .2163286
 corr(M2,M1) |   .4275544   .0714158                      .2780954    .5568003
------------------------------------------------------------------------------

Here, we use the EV[] element type to link the expected value of the response for the model for logb, our repeatedly measured biomarker, directly to survival. Because this is a time-dependent association structure, we must specify the timevar() options to make sure merlin knows which variables represent time in both submodels.

Looking at the results, we get an estimate of \(\alpha = 1.247\) (95% CI: 1.065, 1.429), showing a very strong association between the biomarker and survival, with a hazard ratio of 3.480, indicating a large increase in the mortality rate for a one-unit increase in the log of serum bilirubin, at time t.

Note, this is pretty simple trajectory model for our biomarker. We should investigate whether something more complex is required, by using more fractional polynomials or splines, but we leave that to you. Hint…it does matter (Crowther et al., 2016).

Random effects association structure

Alternatively, we can directly link individual random effects to survival, for example the random intercept,

$$h_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{2} + \alpha b_{0i})$$

This can be fitted with merlin as follows,

. merlin (logb                                       /// long. outcome
>                 fp(time, powers(1 2))              /// fractional polynomial
>                 time#M2[id]@1                      /// random sloope on time
>                 M1[id]@1                           /// random intercept
>                 , family(gaussian))                /// distribution
>        (stime                                      /// survival time
>                 trt                                /// baseline trt
>                 M1[id]                             /// random intercept
>                 , family(weibull, failure(died)))  /// distribution
>        , covariance(unstructured)                  //  VCV of random effects
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2
Fitting fixed effects model:
Fitting full model:
Iteration 0:   log likelihood = -3147.6721  (not concave)
Iteration 1:   log likelihood = -2277.6281  (not concave)
Iteration 2:   log likelihood = -2115.1192  
Iteration 3:   log likelihood = -2021.8568  
Iteration 4:   log likelihood = -1956.9546  
Iteration 5:   log likelihood = -1952.2007  
Iteration 6:   log likelihood = -1952.1805  
Iteration 7:   log likelihood = -1952.1805  
Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -1952.1805
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |            
      fp():1 |   .1600056   .0146768    10.90   0.000     .1312396    .1887716
      fp():2 |   .0028432   .0010681     2.66   0.008     .0007499    .0049366
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5160179   .0581287     8.88   0.000     .4020877     .629948
  sd(resid.) |   .3495502   .0068538                      .3363719    .3632448
-------------+----------------------------------------------------------------
stime:       |            
         trt |    .155547   .1762467     0.88   0.377    -.1898902    .5009842
      M1[id] |    1.32733    .114047    11.64   0.000     1.103802    1.550858
       _cons |  -3.857397   .2747575   -14.04   0.000    -4.395912   -3.318882
  log(gamma) |   .4230823   .0722191     5.86   0.000     .2815355    .5646292
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |    .994177   .0424516                        .91436    1.080961
      sd(M2) |   .1707932   .0118113                      .1491439    .1955851
 corr(M2,M1) |   .4819317   .0732525                        .32613    .6122441
------------------------------------------------------------------------------

We simply include the M1[id] element in both linear predictors, and estimate a coefficient for it in the survival model. Because it’s named the same, merlin knows there are still only two random effects in this model.

Our \(\alpha\) now represents the log hazard ratio for a one-unit increase in the patient-specific deviation from the mean intercept. Given we are now using a time-independent association structure in this example, we no longer require numerical integration to calculate the cumulative hazard function, and therefore don’t need to use the timevar() option. Such an association structure provides us with computational speed gains because of this.

Gradient or rate of change association structure

We can link the current gradient or rate of change of the biomarker directly to survival, for example

$$h_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{2} + \alpha m_{i}^{\prime}(t))$$

where

$$m_{i}^{\prime}(t) = \frac{\partial m_{i}(t)}{\partial t}$$

This can be fitted with merlin using the dEV[] element as follows

. merlin (logb                                      /// long. outcome
>                 fp(time, powers(1 2))             /// fractional polynomial
>                 time#M2[id]@1                     /// random slope on time
>                 M1[id]@1                          /// random intercept
>                 , family(gaussian)                /// distribution
>                 timevar(time))                    /// timevar
>        (stime                                     /// survival time
>                 trt                               /// baseline trt effect
>                 dEV[logb]                         /// current slope
>                 , family(weibull, failure(died))  /// distribution
>                 timevar(stime))                   /// timevar
>        , covariance(unstructured)                 //  VCV of random effects
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2
Fitting fixed effects model:
Fitting full model:
Iteration 0:   log likelihood = -3147.6761  (not concave)
Iteration 1:   log likelihood = -3122.1297  (not concave)
Iteration 2:   log likelihood = -2257.3203  
Iteration 3:   log likelihood = -2208.7113  (not concave)
Iteration 4:   log likelihood =  -2007.873  (not concave)
Iteration 5:   log likelihood = -2003.1171  (not concave)
Iteration 6:   log likelihood = -1990.9806  (not concave)
Iteration 7:   log likelihood = -1962.7313  (not concave)
Iteration 8:   log likelihood = -1955.4192  
Iteration 9:   log likelihood = -1934.4414  
Iteration 10:  log likelihood = -1928.9172  
Iteration 11:  log likelihood = -1928.8515  
Iteration 12:  log likelihood = -1928.8514  
Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -1928.8514
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |            
      fp():1 |   .1889818     .01636    11.55   0.000     .1569168    .2210468
      fp():2 |   .0049364    .001026     4.81   0.000     .0029255    .0069474
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5106389     .05561     9.18   0.000     .4016452    .6196325
  sd(resid.) |   .3519325   .0068616                      .3387378    .3656411
-------------+----------------------------------------------------------------
stime:       |            
         trt |   .0510815   .2218234     0.23   0.818    -.3836843    .4858473
       dEV[] |    11.1193   1.227478     9.06   0.000     8.713491    13.52512
       _cons |  -7.512335   .7239251   -10.38   0.000    -8.931202   -6.093468
  log(gamma) |   .6822592   .1013241     6.73   0.000     .4836676    .8808507
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9586501   .0405812                      .8823226     1.04158
      sd(M2) |   .2225962   .0152073                      .1946998    .2544896
 corr(M2,M1) |   .6767522   .0402024                      .5900355    .7480331
------------------------------------------------------------------------------

where \(\alpha\) now represents the log hazard ratio for a one-unit increase in the gradient of the biomarker, at time t. You can model the biomarker over time as simply or as flexibly as you need – dEV[] will take care of the differentiation for you.

Cumulative association structure

We can assess the cumulative effect of the biomarker directly to survival, for example

$$h_{i}(t) = h_{0}(t) \exp \left( X_{1i}\beta_{2} + \alpha \int_{0}^{t} m_{i}(u) \partial u \right)$$

where the integral is generally approximated using numerical integration. This can be fitted with merlin using the iEV[] element as follows

. merlin (logb                                      /// long. outcome
>                 fp(time, powers(1 2))             /// fractional polynomial
>                 time#M2[id]@1                     /// random slope on time
>                 M1[id]@1                          /// random intercept
>                 , family(gaussian)                /// distribution
>                 timevar(time))                    /// timevar
>        (stime                                     /// survival time
>                 trt                               /// baseline trt effect
>                 iEV[logb]                         /// current integral
>                 , family(weibull, failure(died))  /// distribution
>                 timevar(stime))                   /// timevar
>        , covariance(unstructured)                 //  VCV of random effects
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2
Fitting fixed effects model:
Fitting full model:
Iteration 0:   log likelihood = -3147.6761  (not concave)
Iteration 1:   log likelihood = -2298.6149  (not concave)
Iteration 2:   log likelihood = -2253.9988  
Iteration 3:   log likelihood = -2081.0664  
Iteration 4:   log likelihood = -1994.9057  
Iteration 5:   log likelihood = -1980.1636  
Iteration 6:   log likelihood = -1979.9388  
Iteration 7:   log likelihood = -1979.9382  
Iteration 8:   log likelihood = -1979.9382  
Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -1979.9382
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |            
      fp():1 |   .1549533   .0150315    10.31   0.000      .125492    .1844145
      fp():2 |   .0036447   .0010658     3.42   0.001     .0015558    .0057336
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5138816   .0579815     8.86   0.000       .40024    .6275231
  sd(resid.) |   .3463283   .0067158                      .3334126    .3597444
-------------+----------------------------------------------------------------
stime:       |            
         trt |  -.0997013   .1730687    -0.58   0.565    -.4389098    .2395072
       iEV[] |   .1585397   .0144384    10.98   0.000     .1302409    .1868385
       _cons |  -2.644635   .2010803   -13.15   0.000    -3.038745   -2.250525
  log(gamma) |  -.2274419   .1008161    -2.26   0.024    -.4250379    -.029846
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9919054   .0423277                      .9123192    1.078434
      sd(M2) |   .1797576   .0123527                      .1571064    .2056747
 corr(M2,M1) |   .4187535   .0750906                      .2614919    .5543573
------------------------------------------------------------------------------

where \(\alpha\) now represents the log hazard ratio for a one-unit increase in the cumulative value of the biomarker, at time t.

Combining association structures

In many situations, more than one aspect of the biomarker may be associated with survival, for example, we can link both the current value and slope,

$$h_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{2} + \alpha_{1} m_{i}(t) + \alpha_{2} m_{i}^{\prime}(t))$$

which can be estimated with,

. merlin (logb                                      /// long. outcome
>                 fp(time, powers(1 2))             /// fractional polynomial
>                 time#M2[id]@1                     /// random slope on time
>                 M1[id]@1                          /// random intercept
>                 , family(gaussian)                /// distribution
>                 timevar(time))                    /// timevar
>        (stime                                     /// survival time
>                 trt                               /// baseline trt effect
>                 EV[logb]                          /// current value
>                 dEV[logb]                         /// current slope
>                 , family(weibull, failure(died))  /// distribution
>                 timevar(stime))                   /// timevar
>        , covariance(unstructured)                 //  VCV of random effects
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2
Fitting fixed effects model:
Fitting full model:
Iteration 0:   log likelihood = -3147.6761  (not concave)
Iteration 1:   log likelihood = -3122.0457  (not concave)
Iteration 2:   log likelihood = -2250.5165  
Iteration 3:   log likelihood = -1948.5648  
Iteration 4:   log likelihood = -1919.9815  
Iteration 5:   log likelihood = -1906.7915  
Iteration 6:   log likelihood = -1905.8852  
Iteration 7:   log likelihood = -1905.8663  
Iteration 8:   log likelihood = -1905.8663  
Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -1905.8663
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |            
      fp():1 |   .1706864   .0158662    10.76   0.000     .1395893    .2017836
      fp():2 |   .0043629   .0010441     4.18   0.000     .0023166    .0064093
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5123052   .0577074     8.88   0.000     .3992008    .6254096
  sd(resid.) |   .3445854   .0066243                      .3318436    .3578165
-------------+----------------------------------------------------------------
stime:       |            
         trt |   .0351284    .189808     0.19   0.853    -.3368885    .4071453
        EV[] |   1.002625    .125386     8.00   0.000      .756873    1.248377
       dEV[] |   3.332876   1.010505     3.30   0.001     1.352324    5.313429
       _cons |  -5.189827   .4262606   -12.18   0.000    -6.025283   -4.354372
  log(gamma) |   .1376684   .0999516     1.38   0.168    -.0582331      .33357
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9886591   .0421755                      .9093581    1.074876
      sd(M2) |   .2002244   .0142931                      .1740819    .2302927
 corr(M2,M1) |   .4978334   .0658886                      .3581226    .6157387
------------------------------------------------------------------------------

We can clearly see that both the value of the biomarker, and where it’s going (increasing or decreasing, say) are important.

There are numerous extensions to consider, which we will come back to in later posts. If you want to find out more, take a look at our introductory video, and keep an eye out for updates on our forthcoming training course on joint models.