The standard joint model was described and illustrated with merlin here.
In this example we consider time-dependent effects in the survival submodel. Once again we have,
$$h(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\beta_{1}} + \alpha m_{i}(t))$$
and
$$y_{}(t) = m_{i}(t) + \epsilon_{ij}$$
where
$$m_{i}(t) = X_{2ij}(t)\mathbf{\beta}_{2} + Z_{ij}(t)b_{i}$$
and
$$b_{i} \sim N(0,\Sigma)$$
Given the survival submodel is a proportional hazards model, we are making the assumption that both our hazard ratios $\exp (\beta_{1})$ for the covariates $X_{1ij}$, and the hazard ratio for the effect of our biomarker, $\exp(\alpha)$, are constant across follow-up time. We can relax these assumptions to allow for non-proportional hazards i.e. by modelling time-dependent effects. First we consider the following model,
$$h(t) = h_{0}(t) \exp(X_{1ij}\beta_{1}(t) + \alpha m_{i}(t))$$
where if we assume $X_{1ij}$ is a binary treatment variable, we can allow its associated hazard ratio $\beta_{1}(t)$ to be time-dependent, by making it a function of $\log (t)$
$$h(t) = h_{0}(t) \exp(X_{1ij} \beta_{1} + X_{1ij} \times \beta_{2} \times \log (t) + \alpha m_{i}(t))$$
We can fit this model with merlin as follows
. use https://www.mjcrowther.co.uk/data/jm_example.dta, clear
(Example dataset for joint modelling, Michael Crowther 2011)
. merlin (stime trt trt#fp(stime,pow(0)) EV[logb], ///
> family(weibull, failure(died)) timevar(stime)) ///
> (logb fp(time,pow(1)) fp(time,pow(1))#M2[id]@1 M1[id]@1, ///
> family(gaussian) timevar(time)) ///
> , restartvalues(M2 0.1)
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1
variables created for model 2, component 2: _cmp_2_2_1 to _cmp_2_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.3991
Iteration 1: log likelihood = -2305.879
Iteration 2: log likelihood = -2288.7621
Iteration 3: log likelihood = -1934.0652
Iteration 4: log likelihood = -1932.1905
Iteration 5: log likelihood = -1932.1887
Iteration 6: log likelihood = -1932.1887
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1932.1887
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime: |
trt | .1236849 .2526514 0.49 0.624 -.3715027 .6188725
trt#fp() | -.0929609 .166218 -0.56 0.576 -.4187422 .2328205
EV[] | 1.261099 .0955683 13.20 0.000 1.073788 1.448409
_cons | -4.508606 .3321891 -13.57 0.000 -5.159685 -3.857527
log(gamma) | .063174 .1101584 0.57 0.566 -.1527326 .2790806
-------------+----------------------------------------------------------------
logb: |
fp() | .1690605 .013109 12.90 0.000 .1433674 .1947537
fp()#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .4993858 .0595809 8.38 0.000 .3826093 .6161622
sd(resid.) | .345686 .006619 .3329534 .3589055
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.026386 .0433006 .9449323 1.114861
sd(M2) | .1805098 .0123637 .1578334 .2064441
------------------------------------------------------------------------------
We can also relax the proportional hazards assumption on our association parameter, as follows,
$$h(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\beta_{1}} + \alpha_{1} m_{i}(t) + \alpha_{2} \times \log (t) \times m_{i}(t) )$$
which can be fitted using,
. merlin (stime trt EV[logb] EV[logb]#fp(stime,pow(0)) ///
> , family(weibull, failure(died)) timevar(stime)) ///
> (logb fp(time,pow(1)) fp(time,pow(1))#M2[id]@1 M1[id]@1 ///
> , family(gaussian) timevar(time))
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_1
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1
variables created for model 2, component 2: _cmp_2_2_1 to _cmp_2_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3328.8882 (not concave)
Iteration 1: log likelihood = -2385.0776
Iteration 2: log likelihood = -1962.0446
Iteration 3: log likelihood = -1933.2694
Iteration 4: log likelihood = -1929.9663
Iteration 5: log likelihood = -1929.8696
Iteration 6: log likelihood = -1929.8674
Iteration 7: log likelihood = -1929.8674
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1929.8674
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0760778 .1805074 0.42 0.673 -.2777102 .4298658
EV[] | 1.518905 .1637616 9.28 0.000 1.197938 1.839872
EV[]#fp() | -.2095983 .097755 -2.14 0.032 -.4011945 -.0180021
_cons | -5.362486 .5439833 -9.86 0.000 -6.428673 -4.296298
log(gamma) | .3840729 .1604283 2.39 0.017 .0696392 .6985066
-------------+----------------------------------------------------------------
logb: |
fp() | .1698027 .0132098 12.85 0.000 .1439119 .1956936
fp()#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .5003171 .0595577 8.40 0.000 .3835862 .6170479
sd(resid.) | .3454356 .0066109 .3327185 .3586388
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.025622 .0432787 .9442107 1.114053
sd(M2) | .1815505 .0124789 .1586683 .2077327
------------------------------------------------------------------------------
For completeness we can model non-proportional hazards in both our treatment effect and our association parameter as follows,
. merlin (stime trt trt#fp(stime,pow(0)) EV[logb] EV[logb]#fp(stime,pow(0)), ///
> family(weibull, failure(died)) timevar(stime)) ///
> (logb fp(time,pow(1)) fp(time,pow(1))#M2[id]@1 M1[id]@1, ///
> family(gaussian) timevar(time)) ///
> , restartvalues(M2 0.1)
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1
variables created for model 1, component 4: _cmp_1_4_1 to _cmp_1_4_1
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1
variables created for model 2, component 2: _cmp_2_2_1 to _cmp_2_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.3991
Iteration 1: log likelihood = -2513.1808 (not concave)
Iteration 2: log likelihood = -2179.3897 (not concave)
Iteration 3: log likelihood = -1959.3924 (not concave)
Iteration 4: log likelihood = -1937.9185
Iteration 5: log likelihood = -1930.2868
Iteration 6: log likelihood = -1929.5508
Iteration 7: log likelihood = -1929.5161
Iteration 8: log likelihood = -1929.516
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1929.516
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime: |
trt | .2345127 .2600756 0.90 0.367 -.275226 .7442515
trt#fp() | -.1396854 .1661587 -0.84 0.401 -.4653505 .1859798
EV[] | 1.544048 .1675965 9.21 0.000 1.215564 1.872531
EV[]#fp() | -.2205425 .0990941 -2.23 0.026 -.4147635 -.0263216
_cons | -5.544985 .5856636 -9.47 0.000 -6.692865 -4.397106
log(gamma) | .4437813 .1675397 2.65 0.008 .1154096 .772153
-------------+----------------------------------------------------------------
logb: |
fp() | .1698975 .0132107 12.86 0.000 .1440049 .19579
fp()#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .5003622 .0595612 8.40 0.000 .3836243 .6171001
sd(resid.) | .3454272 .00661 .3327118 .3586286
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.025705 .0432803 .9442896 1.114139
sd(M2) | .1815854 .0124773 .1587055 .2077637
------------------------------------------------------------------------------
In any of the models above, there is no restriction on what functional form we use to model each of the time-dependent
effects. If a more complex form is desired, then we can use the fp() or rcs() syntax to get more flexible time
functions. For example using fractional polynomials,
. merlin (stime trt EV[logb] EV[logb]#fp(stime,pow(1 2)) ///
> , family(weibull, failure(died)) timevar(stime)) ///
> (logb fp(time,pow(1)) fp(time,pow(1))#M2[id]@1 M1[id]@1 ///
> , family(gaussian) timevar(time)) ///
> , restartvalues(M2 0.1)
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_2
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1
variables created for model 2, component 2: _cmp_2_2_1 to _cmp_2_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.7168
Iteration 1: log likelihood = -2579.8759 (not concave)
Iteration 2: log likelihood = -2047.6717
Iteration 3: log likelihood = -1969.8766
Iteration 4: log likelihood = -1931.0458
Iteration 5: log likelihood = -1930.0696
Iteration 6: log likelihood = -1930.066
Iteration 7: log likelihood = -1930.066
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1930.066
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0678679 .1803206 0.38 0.707 -.285554 .4212899
EV[] | 1.489871 .2001501 7.44 0.000 1.097584 1.882158
EV[]#fp():1 | -.0563388 .0658018 -0.86 0.392 -.185308 .0726305
EV[]#fp():2 | .0012085 .0044445 0.27 0.786 -.0075025 .0099195
_cons | -4.948364 .4383433 -11.29 0.000 -5.807501 -4.089227
log(gamma) | .2286889 .1430927 1.60 0.110 -.0517676 .5091454
-------------+----------------------------------------------------------------
logb: |
fp() | .1689365 .0131413 12.86 0.000 .14318 .194693
fp()#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .5000151 .0596196 8.39 0.000 .3831628 .6168673
sd(resid.) | .3455656 .0066158 .3328391 .3587787
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.026716 .0433295 .9452087 1.115251
sd(M2) | .1805794 .0123842 .1578673 .206559
------------------------------------------------------------------------------
or restricted cubic splines,
. merlin (stime trt EV[logb] EV[logb]#rcs(stime, df(3) log event) ///
> , family(weibull, failure(died)) timevar(stime)) ///
> (logb fp(time,pow(1)) fp(time,pow(1))#M2[id]@1 M1[id]@1 ///
> , family(gaussian) timevar(time)) ///
> , restartvalues(M2 0.1)
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_3
variables created for model 2, component 1: _cmp_2_1_1 to _cmp_2_1_1
variables created for model 2, component 2: _cmp_2_2_1 to _cmp_2_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.7168 (not concave)
Iteration 1: log likelihood = -2274.4939
Iteration 2: log likelihood = -1957.0237
Iteration 3: log likelihood = -1931.4924
Iteration 4: log likelihood = -1929.4036
Iteration 5: log likelihood = -1929.3162
Iteration 6: log likelihood = -1929.316
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1929.316
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0860589 .1812669 0.47 0.635 -.2692177 .4413355
EV[] | 1.574589 .2154628 7.31 0.000 1.152289 1.996888
EV[]#rcs():1 | -.15243 .1325297 -1.15 0.250 -.4121834 .1073235
EV[]#rcs():2 | -.0205161 .1181801 -0.17 0.862 -.2521448 .2111126
EV[]#rcs():3 | .0474125 .1849339 0.26 0.798 -.3150512 .4098762
_cons | -5.395841 .5639922 -9.57 0.000 -6.501246 -4.290437
log(gamma) | .3872657 .1669159 2.32 0.020 .0601165 .714415
-------------+----------------------------------------------------------------
logb: |
fp() | .1695191 .0132021 12.84 0.000 .1436436 .1953947
fp()#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .5005547 .0596044 8.40 0.000 .3837323 .6173772
sd(resid.) | .3454686 .0066128 .3327479 .3586755
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.026383 .0433163 .9449014 1.114892
sd(M2) | .1812488 .0124586 .1584038 .2073885
------------------------------------------------------------------------------
I’ve also only shown a time-dependent effect under the current value parameterisation. Any of the other structures shown here can also be made time-dependent in a similar way.
This example used merlin version 2.0.2.