# Multivariate joint longitudinal-survival models

The do file with all the code from this tutorial is available to download here

Joint longitudinal-survival models have been widely developed, but there are many avenues of research where they are lacking in terms of methodological development, and importantly, accessible implementations. Hopefully, merlin fills a few gaps.

A current topic is the ability to model multiple longitudinal outcomes, and link each to survival. Check out Graeme Hickey’s recent review (Hickey et al. (2017)) (Graeme developed joineRML in R which models multiple continuous biomarkers and the Cox survival model, and maintains joineR which is the original implementation of the Henderson et al. (2000) paper). For simplicity, I’ll start by concentrating on an example with two continuous biomarkers, where I want to look at the association between aspects of the biomarker and the risk of event, using either the random effects (time-independent association) or the current value association structures, to show some flexibility in model choice. I’ll use the standard Primary Biliary Cirrhosis (PBC) example.

I model the survival outcome using a Weibull model, so we have:

$$T \sim Weib(\lambda,\gamma)$$ $$Y_{1} \sim N(\mu_{1},\sigma^{2}_{1})$$ $$Y_{2} \sim N(\mu_{2},\sigma^{2}_{2})$$

Now each model for the biomarkers can be as flexible as we like. We can use fractional polynomials with the fp() elements, or restricted cubic splines with the rcs() elements. For simplicity let’s just assume a random intercept and fixed linear slope model for each of them, and work from there. I’m going to link each random intercept to the survival linear predictor.

$$\log (\lambda) = X \beta_{0} + b_{01i} \alpha_{1} + b_{02i} \alpha_{2}$$

We’ll start by loading the data and having a look at a particular patient’s data:

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


. list id logb logp time trt stime died if id==3, noobs

+-----------------------------------------------------------------+
| id       logb       logp      time         trt     stime   died |
|-----------------------------------------------------------------|
|  3   .3364722   2.484907         0   D-penicil   2.77078      1 |
|  3   .0953102   2.484907   .481875   D-penicil         .      . |
|  3   .4054651   2.484907   .996605   D-penicil         .      . |
|  3   .5877866   2.587764   2.03428   D-penicil         .      . |
+-----------------------------------------------------------------+



where logb and logp are the log of serum bilirubin and log of prothrombin index, respectively, time is the time at which the biomarkers were measured, trt is treatment group (either placebo or D-penicillamine), stime is the observed event time, and died is the event indicator. Patient 3 had four repeated measures of both serum bilirubin and prothrombin index, measured at the same time points, recorded in the variable time. They were on treatment and died after 2.77 years.

We fit our model with merlin as follows,

. merlin  (logb  time M1[id]@1              , family(gaussian))                 ///
>         (logp  time M2[id]@1              , family(gaussian))                 ///
>         (stime trt M1[id] M2[id]  , family(weibull, failure(died)))   ///
>           , covariance(unstructured)  restartvalues(M2 0.1)

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -1653.9599  (not concave)
Iteration 1:   log likelihood = -1636.8538  (not concave)
Iteration 2:   log likelihood =  -1620.449  (not concave)
Iteration 3:   log likelihood = -705.13898  (not concave)
Iteration 4:   log likelihood = -596.86757
Iteration 5:   log likelihood = -454.38452
Iteration 6:   log likelihood = -428.44064
Iteration 7:   log likelihood = -426.98648
Iteration 8:   log likelihood = -426.98094
Iteration 9:   log likelihood = -426.98094

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -426.98094
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |
time |   .0977578   .0043158    22.65   0.000      .089299    .1062165
M1[id] |          1          .        .       .            .           .
_cons |   .5801571   .0651011     8.91   0.000     .4525613    .7077529
sd(resid.) |    .491377   .0085925                      .4748214    .5085098
-------------+----------------------------------------------------------------
logp:        |
time |    .012967   .0007165    18.10   0.000     .0115626    .0143713
M2[id] |          1          .        .       .            .           .
_cons |   2.368182   .0050595   468.06   0.000     2.358266    2.378099
sd(resid.) |   .0843888   .0014806                      .0815362    .0873413
-------------+----------------------------------------------------------------
stime:       |
trt |   .0215182   .1801246     0.12   0.905    -.3315196    .3745559
M1[id] |   .8489245   .1324585     6.41   0.000     .5893105    1.108538
M2[id] |   9.413099   2.061308     4.57   0.000      5.37301    13.45319
_cons |  -4.096368   .2912854   -14.06   0.000    -4.667277   -3.525459
log(gamma) |   .5006391   .0723945     6.92   0.000     .3587485    .6425297
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   1.109998   .0471555                      1.021318    1.206378
sd(M2) |   .0740779   .0040597                      .0665335    .0824777
corr(M2,M1) |   .7083502   .0376808                      .6265019    .7747462
------------------------------------------------------------------------------



where we assume a random intercept and fixed linear trend for each of the biomarkers. We specify an unstructured variance-covariance structure through the covariance() option, which lets us estimate the correlation between the two random intercepts. Our parameters _cmp_3_2_1 and _cmp_3_3_1 estimate the association between a one-unit increase in the subject specific deviations from the mean intercepts for log serum bilirubin and log prothrombin, respectively. Both show higher values increase the mortality rate.

We can conduct a sensitivity analysis for the above model by using multivariate $t$-distributed random effects, but first I’ll store the estimates from the previous model, and use those as starting values in this model (I’m estimating the same number of parameters) simply using the from() option, as follows

. mat startvals = e(b)

. set seed 9741308

. merlin  (logb  time M1[id]@1              , family(gaussian))                 ///
>         (logp  time M2[id]@1              , family(gaussian))                 ///
>         (stime trt M1[id] M2[id]          , family(weibull, failure(died)))   ///
>         , covariance(unstructured) intmethod(mcarlo) intpoints(300)           ///
>           redist(t) df(3) from(startvals)

Fitting full model:

Iteration 0:   log likelihood = -472.27797
Iteration 1:   log likelihood = -455.59535
Iteration 2:   log likelihood = -450.51196
Iteration 3:   log likelihood = -443.82075
Iteration 4:   log likelihood = -443.11608
Iteration 5:   log likelihood = -443.08982
Iteration 6:   log likelihood =  -443.0898

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood =  -443.0898
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |
time |    .096978   .0043078    22.51   0.000     .0885349    .1054211
M1[id] |          1          .        .       .            .           .
_cons |   .3301169   .0513636     6.43   0.000     .2294462    .4307877
sd(resid.) |   .4925842   .0086487                      .4759214    .5098303
-------------+----------------------------------------------------------------
logp:        |
time |   .0126674   .0007143    17.73   0.000     .0112675    .0140674
M2[id] |          1          .        .       .            .           .
_cons |   2.356792    .004528   520.49   0.000     2.347918    2.365667
sd(resid.) |   .0847224    .001493                      .0818461    .0876997
-------------+----------------------------------------------------------------
stime:       |
trt |   .0240181   .1763539     0.14   0.892    -.3216292    .3696653
M1[id] |   .8744751   .1286834     6.80   0.000     .6222603     1.12669
M2[id] |   8.316883   1.890573     4.40   0.000     4.611428    12.02234
_cons |  -4.332625    .296188   -14.63   0.000    -4.913143   -3.752107
log(gamma) |   .4792559   .0713973     6.71   0.000     .3393197    .6191921
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   .8244423   .0252386                      .7764304    .8754231
sd(M2) |   .0540816   .0032575                      .0480596    .0608583
corr(M2,M1) |    .717381   .0464357                      .6136249    .7967836
------------------------------------------------------------------------------



where we specify the redistribution(t) and choose a degrees of freedom with df(). If using $t$-distributed random effects, then we have to use Monte Carlo integration. Given the log-likelihood, there doesn’t seem to be evidence that the more robust model is needed, but we should investigate an increasing number of intpoints().

We can now explore alternative association structures, such as the commonly used current value form. In merlin we use the EV[] syntax to link the expected value of an outcome within the linear predictor of another. They are referred to by using the appropriate response variable name within the square brackets, or by the model number (indexed by the order they are specified). Since we are linking a time-dependent expected value to survival, we must explicitly use the fp() or rcs() (for restricted cubic splines) functions when using time, as it must be numerically integrated out in the survival outcome model. For each Gaussian response, we must specify the variable which holds the time of measurements, which does not have to be the same between outcomes, as in this case. Further inbuilt functions include the first and second derivatives of the expected value, with respect to time, using dEV[] and d2EV[], respectively, and the integral of the expected value, using iEV[].

My survival linear predictor is now,

$$\log (\lambda) = X \beta_{0} + E[y_{1}(t)|\mu_{1}(t)] \alpha_{1} + E[y_{2}(t)|\mu_{2}(t)] \alpha_{2}$$

With continuous outcomes, and identity link functions, we have,

$$\log (\lambda) = X \beta_{0} + \mu_{1}(t) \alpha_{1} + \mu_{2}(t) \alpha_{2}$$

We fit our model with,

. merlin  (logb  fp(time,pow(1)) M1[id]@1,                       ///
>                family(gaussian) timevar(time))                 ///
>         (logp  fp(time,pow(1)) M2[id]@1,                       ///
>                family(gaussian) timevar(time))                 ///
>         (stime trt EV[logb] EV[logp],                          ///
>                family(weibull, failure(died)) timevar(stime))  ///
>         , covariance(unstructured) restartvalues(M2 0.1)
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_1
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -1653.9639  (not concave)
Iteration 1:   log likelihood =  -1636.857  (not concave)
Iteration 2:   log likelihood = -1620.4655  (not concave)
Iteration 3:   log likelihood = -751.71258  (not concave)
Iteration 4:   log likelihood = -612.02721  (not concave)
Iteration 5:   log likelihood = -472.58637  (not concave)
Iteration 6:   log likelihood = -447.83119
Iteration 7:   log likelihood = -426.73142
Iteration 8:   log likelihood = -424.24206
Iteration 9:   log likelihood = -423.88791
Iteration 10:  log likelihood = -423.88667
Iteration 11:  log likelihood = -423.88667

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -423.88667
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |
fp() |   .0978945   .0043123    22.70   0.000     .0894426    .1063464
M1[id] |          1          .        .       .            .           .
_cons |   .5806025   .0650906     8.92   0.000     .4530272    .7081778
sd(resid.) |   .4913553   .0085914                      .4748017     .508486
-------------+----------------------------------------------------------------
logp:        |
fp() |   .0129505   .0007125    18.18   0.000      .011554    .0143469
M2[id] |          1          .        .       .            .           .
_cons |   2.368306   .0050532   468.68   0.000     2.358401     2.37821
sd(resid.) |   .0844013   .0014805                      .0815488    .0873536
-------------+----------------------------------------------------------------
stime:       |
trt |   .0347763   .1802431     0.19   0.847    -.3184936    .3880462
EV[] |    .910846   .1311311     6.95   0.000     .6538339    1.167858
EV[] |   9.262327   1.987964     4.66   0.000     5.365988    13.15867
_cons |  -26.49357   4.758907    -5.57   0.000    -35.82086   -17.16629
log(gamma) |   .0977752   .0776863     1.26   0.208    -.0544871    .2500375
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   1.109887   .0471348                      1.021245    1.206223
sd(M2) |   .0739872   .0040438                      .0664713    .0823529
corr(M2,M1) |   .7080859   .0376831                      .6262417    .7744934
------------------------------------------------------------------------------



where _cmp_3_2_1 and _cmp_3_3_1 are now our log hazard ratios for a one-unit increase in the biomarkers, at time $t$, again showing higher values increase the mortality rate.

Due to the flexibility of the complex linear predictor, we can further extend this framework, by investigating whether there is an interaction between the two biomarkers, and the risk of event. Our linear predictor is as follows,

$$\log (\lambda) = X \beta_{0} + \mu_{1}(t) \alpha_{1} + \mu_{2}(t) \alpha_{2} + \mu_{1}(t) \mu_{2}(t) \alpha_{2}$$

. merlin  (logb  fp(time,pow(1)) M1[id]@1,                       ///
>                family(gaussian) timevar(time))                 ///
>         (logp  fp(time,pow(1)) M2[id]@1,                       ///
>                family(gaussian) timevar(time))                 ///
>         (stime trt EV[logb] EV[logp] EV[logb]#EV[logp],        ///
>                family(weibull, failure(died)) timevar(stime))  ///
>          , covariance(unstructured) restartvalues(M2 0.1)
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_1
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -1653.9639  (not concave)
Iteration 1:   log likelihood = -1636.8638  (not concave)
Iteration 2:   log likelihood = -1620.4978  (not concave)
Iteration 3:   log likelihood = -680.83551  (not concave)
Iteration 4:   log likelihood =  -606.4221  (not concave)
Iteration 5:   log likelihood = -505.04613
Iteration 6:   log likelihood = -449.37862  (not concave)
Iteration 7:   log likelihood = -436.08313
Iteration 8:   log likelihood =  -425.1957
Iteration 9:   log likelihood = -422.01442
Iteration 10:  log likelihood = -421.84849
Iteration 11:  log likelihood = -421.77433
Iteration 12:  log likelihood = -421.75105
Iteration 13:  log likelihood = -421.75097
Iteration 14:  log likelihood = -421.75097

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -421.75097
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |
fp() |   .0979839   .0043123    22.72   0.000      .089532    .1064359
M1[id] |          1          .        .       .            .           .
_cons |   .5797618   .0649987     8.92   0.000     .4523667    .7071569
sd(resid.) |   .4913342   .0085901                       .474783    .5084623
-------------+----------------------------------------------------------------
logp:        |
fp() |   .0129836   .0007088    18.32   0.000     .0115943    .0143728
M2[id] |          1          .        .       .            .           .
_cons |    2.36813   .0050434   469.55   0.000     2.358246    2.378015
sd(resid.) |   .0843536   .0014779                      .0815061    .0873006
-------------+----------------------------------------------------------------
stime:       |
trt |   .0145318   .1775198     0.08   0.935    -.3334005    .3624642
EV[] |   10.09101   4.555123     2.22   0.027     1.163134    19.01889
EV[] |   15.27737   3.705786     4.12   0.000     8.014164    22.54058
EV[]#EV[] |  -3.739887   1.851981    -2.02   0.043    -7.369703    -.110072
_cons |  -41.07509    8.99078    -4.57   0.000    -58.69669   -23.45348
log(gamma) |   .0608195   .0804273     0.76   0.450    -.0968151    .2184542
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   1.108243   .0470686                      1.019726    1.204444
sd(M2) |   .0738394   .0040427                      .0663261    .0822037
corr(M2,M1) |   .7041582   .0381703                      .6212948    .7714465
------------------------------------------------------------------------------



There’s a lot of ways to extend this further, with more complex trajectory functions, random effects design matrices, further linking between outcomes, non-continuous longitudinal outcomes…the list goes on. I’ll add more to this later.

This example used merlin version 2.0.2.