Multivariate joint longitudinal-survival models
Extending joint models to handle multiple continuous longitudinal outcomes modelled jointly with a survival outcome.
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 biomarkers 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. Rather than analyse a real dataset, I’ll simulate one, so we know the true values the models should recover.
The model
I model the survival outcome using a Weibull model, and each biomarker with a normal model, so we have:
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.
The data
We’ll simulate 400 patients, each randomised to treatment (trt) with probability 0.5.
. clear
. set seed 3172
. set obs 400
Number of observations (_N) was 0, now 400.
. gen id = _n
. gen trt = runiform() > 0.5Each patient has a random intercept for each biomarker, and , drawn from a bivariate normal distribution with standard deviations 1 and 0.5 and a correlation of 0.5. It’s worth checking what we actually drew:
. matrix Sigma = (1, 0.25 \ 0.25, 0.25) // SDs 1 and 0.5, correlation 0.5
. drawnorm b1 b2, cov(Sigma)
. correlate b1 b2
(obs=400)
| b1 b2
-------------+------------------
b1 | 1.0000
b2 | 0.4183 1.0000The true trajectories are and , so both biomarkers start at zero on average and rise slowly. We simulate survival from a Weibull model with and , a log hazard ratio of −0.5 for treatment, and the current values of both biomarkers in the linear predictor, with log hazard ratios and . That’s the current value model we fit further down. avalon, the merlin family’s simulator, takes the hazard as a function of time with avalon user and does the numerical integration and root-finding for us, and anyone still alive at 10 years is censored.
. avalon user stime died, /// new variables
> hazard(0.03:*1.3:*{t}:^(1.3-1) /// Weibull baseline
> :* exp(0.5 :* (b1 :+ 0.1:*{t}) /// alpha1 x current value 1
> :+ 1.0 :* (b2 :+ 0.05:*{t}))) /// alpha2 x current value 2
> covariates(trt -0.5) /// treatment effect
> maxtime(10) // admin. censoring
note: this hazard is close to singular at the lower limit.
Gauss-Kronrod puts the 15-point quadrature error at 2.3e-06; the run
used nodes() points and is better than that, but in the same regime. Raisi
> ng
nodes() helps only slowly here -- see r(qerror).avalon adds a note that the hazard is close to singular at the lower limit, and that more quadrature nodes help only slowly. It is a cautious check, and here a harmless one: rerunning with nodes(200) instead of the default gives the same simulated times, to within 1e-12. Both biomarkers are then measured at entry and annually until death or censoring, at the same visits, with measurement error standard deviations of 0.5 and 0.25, and we keep a single row of survival data per patient.
. expand 10
(3,600 observations created)
. bys id : gen time = _n - 1
. drop if time > stime
(1,026 observations deleted)
. gen y1 = b1 + 0.1*time + rnormal(0, 0.5)
. gen y2 = b2 + 0.05*time + rnormal(0, 0.25). bys id (time) : replace stime = . if _n > 1
(2574 real changes made, 2574 to missing)
. bys id (time) : replace died = . if _n > 1
(2,574 real changes made, 2,574 to missing)Let’s have a look at a particular patient’s data:
. list id y1 y2 time trt stime died if id==15, noobs
+-----------------------------------------------------------+
| id y1 y2 time trt stime died |
|-----------------------------------------------------------|
| 15 -1.407958 .8162562 0 1 3.0939363 1 |
| 15 -1.106535 .4909965 1 1 . . |
| 15 -.0829626 .385632 2 1 . . |
| 15 -.6855976 .1765713 3 1 . . |
+-----------------------------------------------------------+where y1 and y2 are our two simulated biomarkers, time is the time at which the biomarkers were measured, trt is treatment group (0 for placebo, 1 for treated), stime is the observed survival time (death or censoring), and died is the event indicator. Patient 15 had four repeated measures of both biomarkers, measured at the same time points, recorded in the variable time. They were on treatment and died after 3.09 years.
Fitting the model
We fit our model with merlin as follows,
. merlin (y1 time M1[id]@1 ///
> , family(gaussian)) ///
> (y2 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 = -6394.0019
Iteration 1: Log likelihood = -4616.4127
Iteration 2: Log likelihood = -4389.3962
Iteration 3: Log likelihood = -4353.7652
Iteration 4: Log likelihood = -4353.7096
Iteration 5: Log likelihood = -4353.7096
Mixed effects regression model Number of obs = 2,974
Log likelihood = -4353.7096
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
y1: |
time | .1080563 .0036577 29.54 0.000 .1008874 .1152253
M1[id] | 1 . . . . .
_cons | .0589698 .0524317 1.12 0.261 -.0437945 .1617341
sd(resid.) | .5096345 .0070989 .495909 .5237398
-------------+----------------------------------------------------------------
y2: |
time | .0480728 .0018057 26.62 0.000 .0445338 .0516119
M2[id] | 1 . . . . .
_cons | .030193 .0262409 1.15 0.250 -.0212382 .0816242
sd(resid.) | .2516678 .0035072 .2448868 .2586366
-------------+----------------------------------------------------------------
stime: |
trt | -.5767498 .1321037 -4.37 0.000 -.8356682 -.3178314
M1[id] | .5391087 .0787494 6.85 0.000 .3847627 .6934547
M2[id] | .8215445 .1464382 5.61 0.000 .5345308 1.108558
_cons | -3.443443 .2098124 -16.41 0.000 -3.854668 -3.032218
log(gamma) | .4517125 .05714 7.91 0.000 .3397201 .5637048
-------------+----------------------------------------------------------------
id: |
sd(M1) | .9965954 .0372234 .9262454 1.072289
sd(M2) | .499473 .0186409 .4642418 .537378
corr(M2,M1) | .4086983 .0441872 .3185883 .4914875
------------------------------------------------------------------------------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] are the log hazard ratios for a one-unit increase in the subject-specific deviations from the mean intercepts for y1 and y2, respectively. Both show that higher values increase the mortality rate.
We simulated from the current value of each biomarker, not its random intercept, but with a fixed slope the two differ only by a function of time that’s the same for everyone: . So the model we simulated from is a random intercept model with the same log hazard ratios, and , and a baseline hazard multiplied by . A Weibull can only approximate that baseline, but the estimates are close to the truth: 0.539 and 0.822 against 0.5 and 1, and −0.577 for treatment (−0.5). The longitudinal submodels recover their slopes, 0.108 and 0.048 (0.1 and 0.05), their residual standard deviations, 0.510 and 0.252 (0.5 and 0.25), and the random intercept standard deviations, 0.997 and 0.499 (1 and 0.5). The estimated correlation, 0.409, is about two standard errors below the 0.5 we simulated from, but it’s close to the 0.418 that our 400 draws happened to have, as correlate showed above.
A sensitivity analysis with t-distributed random effects
We can conduct a sensitivity analysis for the above model by using 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)
. merlin (y1 time M1[id]@1 ///
> , family(gaussian)) ///
> (y2 time M2[id]@1 ///
> , family(gaussian)) ///
> (stime trt M1[id] M2[id] ///
> , family(weibull ///
> , failure(died))) ///
> , covariance(unstructured) ///
> redistribution(t) df(3) ///
> from(startvals) //
Fitting full model:
Iteration 0: Log likelihood = -4393.1646
Iteration 1: Log likelihood = -4374.995
Iteration 2: Log likelihood = -4374.93
Iteration 3: Log likelihood = -4374.93
Mixed effects regression model Number of obs = 2,974
Log likelihood = -4374.93
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
y1: |
time | .1080629 .0036572 29.55 0.000 .1008949 .1152309
M1[id] | 1 . . . . .
_cons | .0606781 .0519693 1.17 0.243 -.0411798 .162536
sd(resid.) | .50984 .007107 .4960991 .5239615
-------------+----------------------------------------------------------------
y2: |
time | .0479611 .0018056 26.56 0.000 .0444222 .0515
M2[id] | 1 . . . . .
_cons | .0221629 .0255875 0.87 0.386 -.0279878 .0723135
sd(resid.) | .2517612 .0035109 .2449731 .2587373
-------------+----------------------------------------------------------------
stime: |
trt | -.5767442 .1320987 -4.37 0.000 -.8356528 -.3178356
M1[id] | .5404562 .0784645 6.89 0.000 .3866687 .6942437
M2[id] | .8060864 .1453929 5.54 0.000 .5211215 1.091051
_cons | -3.450372 .2101319 -16.42 0.000 -3.862223 -3.038521
log(gamma) | .4524762 .0571365 7.92 0.000 .3404907 .5644617
-------------+----------------------------------------------------------------
id: |
sd(M1) | .8029791 .0389018 .730241 .8829626
sd(M2) | .396539 .019387 .360305 .4364169
corr(M2,M1) | .4413838 .0505132 .3372526 .5348449
------------------------------------------------------------------------------where we specify redistribution(t) and choose the degrees of freedom, here 3, with df(). merlin integrates over the t-distributed random effects with its default adaptive Gauss–Hermite quadrature, with the t density carried in the weights, so nothing more needs specifying, and the model converges in three iterations. Its log-likelihood, −4374.9, is about 21 lower than the −4353.7 of the model with normally distributed random effects, which has the same number of parameters, so the heavier-tailed model isn’t needed here. That is as it should be: we simulated normally distributed random effects.
Alternative association structures
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. Outcomes 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 use timevar() to specify the variable which holds the time of measurements in each submodel; it doesn’t have to be the same between outcomes, and here it isn’t (time and stime). 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, for the Weibull scale , which now varies over time, is,
With continuous outcomes, and identity link functions, we have,
We fit our model with,
. merlin (y1 time M1[id]@1, ///
> family(gaussian) ///
> timevar(time)) ///
> (y2 time M2[id]@1, ///
> family(gaussian) ///
> timevar(time)) ///
> (stime trt ///
> EV[y1] EV[y2] ///
> , family(weibull ///
> , failure(died)) ///
> timevar(stime)) ///
> , covariance(unstructured) ///
> technique(bfgs) //
Fitting fixed effects model:
Fitting full model:
Iteration 0: Log likelihood = -6394.0019
Iteration 1: Log likelihood = -6162.6139 (backed up)
Iteration 2: Log likelihood = -5789.1623 (backed up)
Iteration 3: Log likelihood = -5446.6342 (backed up)
Iteration 4: Log likelihood = -5437.3814 (backed up)
Iteration 5: Log likelihood = -5406.0928 (backed up)
Iteration 6: Log likelihood = -5398.6417 (backed up)
Iteration 7: Log likelihood = -5372.536 (backed up)
Iteration 8: Log likelihood = -5364.657 (backed up)
Iteration 9: Log likelihood = -5344.7101 (backed up)
Iteration 10: Log likelihood = -5296.3359 (backed up)
Iteration 11: Log likelihood = -5149.9826 (backed up)
Iteration 12: Log likelihood = -4922.5186 (backed up)
Iteration 13: Log likelihood = -4457.4156
Iteration 14: Log likelihood = -4440.7317
Iteration 15: Log likelihood = -4416.5674
Iteration 16: Log likelihood = -4383.2556
Iteration 17: Log likelihood = -4353.1663
Iteration 18: Log likelihood = -4351.427
Iteration 19: Log likelihood = -4351.1968
Iteration 20: Log likelihood = -4351.1852
Iteration 21: Log likelihood = -4351.1833
Iteration 22: Log likelihood = -4351.1831
Iteration 23: Log likelihood = -4351.1831
Mixed effects regression model Number of obs = 2,974
Log likelihood = -4351.1831
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
y1: |
time | .1082094 .0036569 29.59 0.000 .101042 .1153767
M1[id] | 1 . . . . .
_cons | .0587273 .052435 1.12 0.263 -.0440434 .161498
sd(resid.) | .5096451 .0070994 .4959188 .5237514
-------------+----------------------------------------------------------------
y2: |
time | .0481309 .0018057 26.66 0.000 .0445918 .0516699
M2[id] | 1 . . . . .
_cons | .0300896 .0262418 1.15 0.252 -.0213433 .0815226
sd(resid.) | .2516724 .0035074 .2448911 .2586416
-------------+----------------------------------------------------------------
stime: |
trt | -.5834369 .1322134 -4.41 0.000 -.8425704 -.3243034
EV[] | .5542788 .0791606 7.00 0.000 .3991269 .7094307
EV[] | .844021 .1468829 5.75 0.000 .5561357 1.131906
_cons | -3.384123 .20124 -16.82 0.000 -3.778547 -2.9897
log(gamma) | .2405593 .0626121 3.84 0.000 .1178417 .3632768
-------------+----------------------------------------------------------------
id: |
sd(M1) | .9966968 .0372265 .9263407 1.072396
sd(M2) | .4995017 .0186419 .4642686 .5374087
corr(M2,M1) | .4087873 .0441834 .3186843 .4915686
------------------------------------------------------------------------------I’ve added technique(bfgs) to this model and the next. With these data, merlin’s default Newton–Raphson optimiser, starting from its default starting values, stalls in a region where the log-likelihood isn’t concave, and gives up altogether on the model with the interaction below. The BFGS algorithm, which builds up its own approximation to the curvature as it goes, converges on both. If a joint model won’t converge, trying another optimiser, or better starting values, is a good first step.
The coefficients on EV[y1] and EV[y2], the first and second EV[] rows of the output, are our log hazard ratios for a one-unit increase in the biomarkers, at time t, again showing higher values increase the mortality rate. This is the model we simulated from, and the estimates are 0.554 (95% CI: 0.399, 0.709) and 0.844 (0.556, 1.132), against true values of 0.5 and 1, with −0.583 for treatment (−0.5). As we’d expect from the argument above, they’re close to the random intercept model’s, and the log-likelihood improves a little, from −4353.7 to −4351.2, because the baseline hazard is now the right shape. You can see where that shows: in the random intercept model the Weibull shape had to stretch to mimic the extra time trend, with a log(gamma) of 0.452, whereas here it’s 0.241, close to .
Due to the flexibility of the complex linear predictor, we can further extend this framework, by investigating whether the two biomarkers interact in their effect on the risk of event. Our linear predictor is as follows,
. merlin (y1 time M1[id]@1, ///
> family(gaussian) ///
> timevar(time)) ///
> (y2 time M2[id]@1, ///
> family(gaussian) ///
> timevar(time)) ///
> (stime trt ///
> EV[y1] EV[y2] ///
> EV[y1]#EV[y2] ///
> , family(weibull ///
> , failure(died)) ///
> timevar(stime)) ///
> , covariance(unstructured) ///
> technique(bfgs) //
Fitting fixed effects model:
Fitting full model:
Iteration 0: Log likelihood = -6394.0019
Iteration 1: Log likelihood = -6162.4783 (backed up)
Iteration 2: Log likelihood = -5789.4883 (backed up)
Iteration 3: Log likelihood = -5446.4986 (backed up)
Iteration 4: Log likelihood = -5436.4454 (backed up)
Iteration 5: Log likelihood = -5404.5247 (backed up)
Iteration 6: Log likelihood = -5397.1618 (backed up)
Iteration 7: Log likelihood = -5373.133 (backed up)
Iteration 8: Log likelihood = -5368.3959 (backed up)
Iteration 9: Log likelihood = -5343.0615 (backed up)
Iteration 10: Log likelihood = -5247.906 (backed up)
Iteration 11: Log likelihood = -4933.2608 (backed up)
Iteration 12: Log likelihood = -4830.2631 (backed up)
Iteration 13: Log likelihood = -4658.5529 (backed up)
Iteration 14: Log likelihood = -4505.3595
Iteration 15: Log likelihood = -4472.4248
Iteration 16: Log likelihood = -4469.1097
Iteration 17: Log likelihood = -4380.3735
Iteration 18: Log likelihood = -4361.67
Iteration 19: Log likelihood = -4353.4678
Iteration 20: Log likelihood = -4351.0097
Iteration 21: Log likelihood = -4350.7813
Iteration 22: Log likelihood = -4350.7362
Iteration 23: Log likelihood = -4350.7241
Iteration 24: Log likelihood = -4350.7225
Iteration 25: Log likelihood = -4350.7222
Iteration 26: Log likelihood = -4350.7222
Mixed effects regression model Number of obs = 2,974
Log likelihood = -4350.7222
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
y1: |
time | .10824 .0036565 29.60 0.000 .1010734 .1154065
M1[id] | 1 . . . . .
_cons | .0587146 .0524282 1.12 0.263 -.0440428 .161472
sd(resid.) | .5096338 .007099 .4959082 .5237393
-------------+----------------------------------------------------------------
y2: |
time | .0481771 .0018061 26.68 0.000 .0446372 .0517169
M2[id] | 1 . . . . .
_cons | .0300815 .0262408 1.15 0.252 -.0213495 .0815125
sd(resid.) | .2516653 .0035072 .2448844 .258634
-------------+----------------------------------------------------------------
stime: |
trt | -.5739741 .1325181 -4.33 0.000 -.8337047 -.3142434
EV[] | .6122112 .100526 6.09 0.000 .4151838 .8092386
EV[] | .9890418 .2119836 4.67 0.000 .5735616 1.404522
EV[]#EV[] | -.1415821 .1490245 -0.95 0.342 -.4336648 .1505006
_cons | -3.415263 .2057831 -16.60 0.000 -3.81859 -3.011935
log(gamma) | .2348909 .063083 3.72 0.000 .1112505 .3585313
-------------+----------------------------------------------------------------
id: |
sd(M1) | .9965619 .0372219 .9262147 1.072252
sd(M2) | .4994846 .0186415 .4642522 .5373909
corr(M2,M1) | .4074803 .0442711 .3172088 .490434
------------------------------------------------------------------------------We simulated no interaction, and the model agrees: the interaction is estimated as −0.142 (95% CI: −0.434, 0.151), comfortably compatible with zero, and the log-likelihood barely moves, from −4351.2 to −4350.7. With the interaction in the model, the main effects, 0.612 and 0.989, are the log hazard ratios for each biomarker when the other is zero, which is where both start on average.
There are lots 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.
MethodsJoint longitudinal–survival models Longitudinal data analysis