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_{1}^{2})$$
$$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.
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)]\alpha_{1} + E[y_{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) \times \mu_{2}(t)\alpha_{3}$$
. 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.
Latest Resources
Videos
State-of-the-art statistical models for modern HTA
Videos