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.4547  (not concave)
Iteration 3:   log likelihood =  -702.3799  (not concave)
Iteration 4:   log likelihood =  -628.0059  (not concave)
Iteration 5:   log likelihood = -581.70293  
Iteration 6:   log likelihood =  -457.1158  
Iteration 7:   log likelihood = -427.20592  
Iteration 8:   log likelihood = -426.98114  
Iteration 9:   log likelihood = -426.98094  
Iteration 10:  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 |    .021518   .1801246     0.12   0.905    -.3315197    .3745558
      M1[id] |   .8489245   .1324586     6.41   0.000     .5893105    1.108538
      M2[id] |   9.413098   2.061308     4.57   0.000     5.373009    13.45319
       _cons |  -4.096368   .2912855   -14.06   0.000    -4.667277   -3.525459
  log(gamma) |   .5006391   .0723945     6.92   0.000     .3587484    .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.53991  
Iteration 3:   log likelihood = -443.82154  
Iteration 4:   log likelihood = -443.11626  
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   .0464356                      .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.964  (not concave)
Iteration 1:   log likelihood =  -1636.857  (not concave)
Iteration 2:   log likelihood = -1620.4807  (not concave)
Iteration 3:   log likelihood = -690.03146  (not concave)
Iteration 4:   log likelihood = -581.00553  
Iteration 5:   log likelihood = -459.81099  (not concave)
Iteration 6:   log likelihood = -437.64469  
Iteration 7:   log likelihood = -425.50314  (not concave)
Iteration 8:   log likelihood = -424.21738  
Iteration 9:   log likelihood =  -423.8876  
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    .180243     0.19   0.847    -.3184935    .3880462
        EV[] |   .9108463   .1311305     6.95   0.000     .6538352    1.167857
        EV[] |    9.26232   1.987933     4.66   0.000     5.366043     13.1586
       _cons |  -26.49356   4.758831    -5.57   0.000    -35.82069   -17.16642
  log(gamma) |   .0977752   .0776863     1.26   0.208     -.054487    .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.964  (not concave)
Iteration 1:   log likelihood = -1636.8638  (not concave)
Iteration 2:   log likelihood = -1620.4578  (not concave)
Iteration 3:   log likelihood = -714.04208  (not concave)
Iteration 4:   log likelihood = -610.50966  (not concave)
Iteration 5:   log likelihood = -466.59374  
Iteration 6:   log likelihood = -436.75564  
Iteration 7:   log likelihood = -430.20305  (not concave)
Iteration 8:   log likelihood = -429.55238  
Iteration 9:   log likelihood = -423.87796  
Iteration 10:  log likelihood = -421.94622  
Iteration 11:  log likelihood = -421.75216  
Iteration 12:  log likelihood = -421.75097  
Iteration 13:  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 |   .0145316   .1775196     0.08   0.935    -.3334005    .3624637
        EV[] |   10.09113   4.554526     2.22   0.027      1.16442    19.01784
        EV[] |   15.27746   3.704811     4.12   0.000     8.016167    22.53876
   EV[]#EV[] |  -3.739935    1.85173    -2.02   0.043    -7.369259   -.1106113
       _cons |  -41.07531   8.988426    -4.57   0.000     -58.6923   -23.45832
  log(gamma) |   .0608188   .0804258     0.76   0.450    -.0968128    .2184504
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   1.108243   .0470686                      1.019726    1.204444
      sd(M2) |   .0738394   .0040427                      .0663261    .0822037
 corr(M2,M1) |   .7041582   .0381703                      .6212947    .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.

Related

comments powered by Disqus