# Multivariate joint longitudinal-survival models

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. We think merlin fills a few gaps.

In this post, we’ll take a look at the extension to modelling multiple continuous longitudinal outcomes, jointly with survival. For simplicity, I’ll concentrate on an example with two continuous biomarkers, where we 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)      //

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.451  (not concave)
Iteration 3:   log likelihood = -703.30831  (not concave)
Iteration 4:   log likelihood =  -595.3835
Iteration 5:   log likelihood = -453.10824
Iteration 6:   log likelihood = -428.34984
Iteration 7:   log likelihood = -426.98677
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
------------------------------------------------------------------------------
| Coefficient  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    -.3315195    .3745559
M1[id] |   .8489245   .1324586     6.41   0.000     .5893105    1.108539
M2[id] |   9.413099   2.061308     4.57   0.000     5.373009    13.45319
_cons |  -4.096369   .2912855   -14.06   0.000    -4.667278    -3.52546
log(gamma) |   .5006391   .0723945     6.92   0.000     .3587485    .6425298
-------------+----------------------------------------------------------------
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 coefficients on M1[id] and M2[id] 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.47338
Iteration 3:   log likelihood = -443.81745
Iteration 4:   log likelihood = -443.11497
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
------------------------------------------------------------------------------
| Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |
time |    .096978   .0043078    22.51   0.000     .0885349    .1054211
M1[id] |          1          .        .       .            .           .
_cons |    .330117   .0513636     6.43   0.000     .2294462    .4307879
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    .3696654
M1[id] |   .8744751   .1286834     6.80   0.000     .6222603     1.12669
M2[id] |   8.316883   1.890573     4.40   0.000     4.611427    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    .8754232
sd(M2) |   .0540816   .0032575                      .0480596    .0608583
corr(M2,M1) |   .7173809   .0464357                      .6136248    .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 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(t)) = 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(t)) = X \beta_{0} + \mu_{1}(t) \alpha_{1} + \mu_{2}(t) \alpha_{2}$

We fit our model with,

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

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.4631  (not concave)
Iteration 3:   log likelihood = -729.09743  (not concave)
Iteration 4:   log likelihood = -600.43545  (not concave)
Iteration 5:   log likelihood = -467.76905  (not concave)
Iteration 6:   log likelihood =  -445.6846
Iteration 7:   log likelihood = -431.65116
Iteration 8:   log likelihood = -424.61221
Iteration 9:   log likelihood = -424.14093
Iteration 10:  log likelihood = -423.88744
Iteration 11:  log likelihood = -423.88667
Iteration 12:  log likelihood = -423.88667

Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -423.88667
------------------------------------------------------------------------------
| Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |
time |   .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:        |
time |   .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    .180243     0.19   0.847    -.3184936    .3880461
EV[] |   .9108469   .1311328     6.95   0.000     .6538313    1.167862
EV[] |   9.262304   1.988052     4.66   0.000     5.365795    13.15881
_cons |  -26.49352   4.759115    -5.57   0.000    -35.82121   -17.16583
log(gamma) |   .0977754   .0776864     1.26   0.208    -.0544871    .2500379
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   1.109887   .0471348                      1.021245    1.206223
sd(M2) |   .0739872   .0040438                      .0664713    .0823529
corr(M2,M1) |   .7080859   .0376832                      .6262417    .7744934
------------------------------------------------------------------------------

where the coefficients on EV[logb] and EV[logp] are 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(t)) = X \beta_{0} + \mu_{1}(t) \alpha_{1} + \mu_{2}(t) \alpha_{2} + \mu_{1}(t) \mu_{2}(t) \alpha_{2}$

. merlin  (logb  time M1[id]@1,           ///
>                family(gaussian)         ///
>                timevar(time))           ///
>         (logp  time 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)     //

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.4626  (not concave)
Iteration 3:   log likelihood = -707.63326  (not concave)
Iteration 4:   log likelihood = -645.09357  (not concave)
Iteration 5:   log likelihood = -602.03532
Iteration 6:   log likelihood = -455.90324
Iteration 7:   log likelihood = -430.76596
Iteration 8:   log likelihood = -423.47881
Iteration 9:   log likelihood = -421.75788
Iteration 10:  log likelihood = -421.75097
Iteration 11:  log likelihood = -421.75097

Mixed effects regression model                           Number of obs = 1,945
Log likelihood = -421.75097
------------------------------------------------------------------------------
| Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
logb:        |
time |    .097984   .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:        |
time |   .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 |   .0145312   .1775196     0.08   0.935    -.3334008    .3624633
EV[] |    10.0914   4.554614     2.22   0.027     1.164524    19.01828
EV[] |   15.27767   3.704862     4.12   0.000     8.016271    22.53906
EV[]#EV[] |  -3.740047   1.851765    -2.02   0.043     -7.36944   -.1106548
_cons |   -41.0758   8.988554    -4.57   0.000    -58.69305   -23.45856
log(gamma) |   .0608178   .0804258     0.76   0.450    -.0968138    .2184495
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   1.108243   .0470686                      1.019726    1.204444
sd(M2) |   .0738394   .0040427                      .0663261    .0822037
corr(M2,M1) |   .7041581   .0381703                      .6212946    .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.

1. Ying Shang says: