# 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)\mathbf{\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}\mathbf{\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(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\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(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\beta_{2}} + \alpha m^{\prime}_{i}(t))$

where $m^{\prime}_{i}(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(t) = h_{0}(t) \exp \left( X_{1ij}\mathbf{\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(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\beta_{2}} + \alpha_{1} m_{i}(t)+ \alpha_{2} m^{\prime}_{i}(t))$

which can be estimed 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.