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.