# Joint longitudinal-survival models with time-dependent effects (non-proportional hazards)

In this post we’ll focus on how to model time-dependent effects (non-proportional hazards), specifically within a joint longitudinal-survival model.

Now joint models are becoming commonplace in medical research, but as always, the fundamentals still matter, and indeed are often ignored. We’re going to look at how to account for time-dependency in both baseline covariates in the survival submodel, and allowing for the association (between the longitudinal and survival models) to vary over time. In my experience, these two assumptions are almost never considered – part of the reason is a lack of software allowing easy application and testing. Let’s remedy that.

First, I recommend you take a read of this introduction to joint models, which defines the standard joint model. Let’s define it here again anyway. Our proportional hazards submodel is, assuming the current value association structure, is

$$h_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{1} \alpha m(t))$$

with our longitudinal submodel,

$$y(t) = m_{i}(t) + \epsilon_{ij}$$

where

$$m_{i}(t) = X_{2i}(t)\beta_{2} + Z_{i}(t)b_{i}$$

and we have our normally distributed random effects,

$$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_{1i}(t)\), 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_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{1}(t) \alpha m(t))$$

where if we assume \(X_{1i}\) 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_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{1} + X_{1i} \times\beta_{2} \times \log(t) + \alpha m(t))$$

We can fit this model with `merlin`

as follows,

```
. //load example dataset
. use https://www.mjcrowther.co.uk/data/jm_example.dta, clear
(Example dataset for joint modelling, Michael Crowther 2011)
. //estimate our joint model
. merlin (stime /// surv. time
> trt /// treatment
> trt#fp(stime,pow(0)) /// treatment x log(time)
> EV[logb] /// current val.
> , family(weibull, failure(died)) /// surv. dist.
> timevar(stime)) /// timevar
> (logb /// biomarker
> time /// fixed linear time
> time#M2[id]@1 /// random linear trend
> M1[id]@1 /// random intercept
> , family(gaussian) /// distribution
> timevar(time)) // timevar
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.3991
Iteration 1: log likelihood = -2305.8791
Iteration 2: log likelihood = -2276.4815
Iteration 3: log likelihood = -1932.9473
Iteration 4: log likelihood = -1932.1891
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
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | .1236873 .2526511 0.49 0.624 -.3714998 .6188744
trt#fp() | -.0929626 .1662177 -0.56 0.576 -.4187434 .2328182
EV[] | 1.261099 .0955683 13.20 0.000 1.073788 1.448409
_cons | -4.50861 .3321885 -13.57 0.000 -5.159687 -3.857532
log(gamma) | .0631754 .1101581 0.57 0.566 -.1527305 .2790813
-------------+----------------------------------------------------------------
logb: |
time | .1690605 .013109 12.90 0.000 .1433674 .1947537
time#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .4993858 .0595809 8.38 0.000 .3826093 .6161623
sd(resid.) | .345686 .006619 .3329534 .3589055
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.026386 .0433006 .9449323 1.114861
sd(M2) | .1805098 .0123637 .1578334 .2064441
------------------------------------------------------------------------------
```

where `trt#fp(stime,pow(0))`

forms an interaction between our binary treatment variable, `trt`

, and a fractional polynomial of time, with power 0, which gives us the log of time, as desired.

We can also relax the proportional hazards assumption on our association parameter, as follows,

$$h_{i}(t) = h_{0}(t) \exp (X_{1i}\beta_{1} + \alpha_{1} m(t) + \alpha_{2} \times \log(t) \times m_{i}(t))$$

which can be fitted using,

```
. //estimate our joint model
. merlin (stime /// surv. time
> trt /// treatment
> EV[logb] /// current val.
> EV[logb]#fp(stime,pow(0)) /// curr. val. x log(time)
> , family(weibull, failure(died)) /// surv. dist.
> timevar(stime)) /// timevar
> (logb /// biomarker
> time /// fixed linear time
> time#M2[id]@1 /// random linear trend
> M1[id]@1 /// random intercept
> , family(gaussian) /// distribution
> timevar(time)) // timevar
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_1
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.7168
Iteration 1: log likelihood = -2571.7802 (not concave)
Iteration 2: log likelihood = -1989.3788
Iteration 3: log likelihood = -1934.5506
Iteration 4: log likelihood = -1929.9274
Iteration 5: log likelihood = -1929.8677
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
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0760779 .1805076 0.42 0.673 -.2777105 .4298663
EV[] | 1.518904 .1637653 9.27 0.000 1.19793 1.839878
EV[]#fp() | -.2095974 .0977579 -2.14 0.032 -.4011993 -.0179954
_cons | -5.362482 .5439993 -9.86 0.000 -6.428701 -4.296263
log(gamma) | .3840707 .1604339 2.39 0.017 .069626 .6985155
-------------+----------------------------------------------------------------
logb: |
time | .1698027 .0132098 12.85 0.000 .1439119 .1956936
time#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .500317 .0595577 8.40 0.000 .3835861 .6170479
sd(resid.) | .3454356 .0066109 .3327185 .3586388
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.025622 .0432787 .9442108 1.114054
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,

```
. //estimate our joint model
. merlin (stime /// surv. time
> trt /// treatment
> trt#fp(stime,pow(0)) /// treatment x log(time)
> EV[logb] /// current val.
> EV[logb]#fp(stime,pow(0)) /// curr. val. x log(time)
> , family(weibull, failure(died)) /// surv. dist.
> timevar(stime)) /// timevar
> (logb /// biomarker
> time /// fixed linear time
> time#M2[id]@1 /// random linear trend
> M1[id]@1 /// random intercept
> , family(gaussian) /// distribution
> timevar(time)) // timevar
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
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.3991
Iteration 1: log likelihood = -2513.1771 (not concave)
Iteration 2: log likelihood = -2179.3925 (not concave)
Iteration 3: log likelihood = -1959.5643 (not concave)
Iteration 4: log likelihood = -1937.9958
Iteration 5: log likelihood = -1930.2868
Iteration 6: log likelihood = -1929.5501
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
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | .2345152 .2600759 0.90 0.367 -.2752241 .7442545
trt#fp() | -.1396865 .1661589 -0.84 0.401 -.4653521 .185979
EV[] | 1.544048 .167597 9.21 0.000 1.215564 1.872532
EV[]#fp() | -.2205428 .0990945 -2.23 0.026 -.4147644 -.0263211
_cons | -5.544988 .5856658 -9.47 0.000 -6.692872 -4.397104
log(gamma) | .443782 .1675403 2.65 0.008 .1154091 .7721548
-------------+----------------------------------------------------------------
logb: |
time | .1698975 .0132107 12.86 0.000 .1440049 .19579
time#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 a higher degree fractional polynomial,

```
. merlin (stime /// surv. time
> trt /// treatment
> EV[logb] /// current val.
> EV[logb]#fp(stime,pow(1 2)) /// curr. val. x FP2
> , family(weibull, failure(died)) /// surv. dist.
> timevar(stime)) /// timevar
> (logb /// biomarker
> time /// fixed linear time
> time#M2[id]@1 /// random linear trend
> M1[id]@1 /// random intercept
> , family(gaussian) /// distribution
> timevar(time)) // timevar
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_2
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -3147.7168
Iteration 1: log likelihood = -2579.8758 (not concave)
Iteration 2: log likelihood = -2047.692
Iteration 3: log likelihood = -1966.9529
Iteration 4: log likelihood = -1931.3611
Iteration 5: log likelihood = -1930.0721
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
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0678671 .1803206 0.38 0.707 -.2855547 .421289
EV[] | 1.489871 .2001476 7.44 0.000 1.097589 1.882153
EV[]#fp():1 | -.0563393 .0658008 -0.86 0.392 -.1853065 .0726279
EV[]#fp():2 | .0012085 .0044444 0.27 0.786 -.0075024 .0099194
_cons | -4.948364 .4383387 -11.29 0.000 -5.807492 -4.089236
log(gamma) | .2286896 .1430908 1.60 0.110 -.0517633 .5091425
-------------+----------------------------------------------------------------
logb: |
time | .1689365 .0131413 12.86 0.000 .14318 .194693
time#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 /// surv. time
> trt /// treatment
> EV[logb] /// current val.
> EV[logb]#rcs(stime, df(3) log event) /// curr. val. x rcs
> , family(weibull, failure(died)) /// surv. dist.
> timevar(stime)) /// timevar
> (logb /// biomarker
> time /// fixed linear time
> time#M2[id]@1 /// random linear trend
> M1[id]@1 /// random intercept
> , family(gaussian) /// distribution
> timevar(time)) // timevar
variables created for model 1, component 3: _cmp_1_3_1 to _cmp_1_3_3
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 = -1960.2108
Iteration 3: log likelihood = -1931.6713
Iteration 4: log likelihood = -1929.4128
Iteration 5: log likelihood = -1929.3164
Iteration 6: log likelihood = -1929.316
Iteration 7: log likelihood = -1929.316
Mixed effects regression model Number of obs = 1,945
Log likelihood = -1929.316
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | .0860604 .1812674 0.47 0.635 -.2692172 .441338
EV[] | 1.574606 .2154759 7.31 0.000 1.152281 1.996931
EV[]#rcs():1 | -.1524371 .1325329 -1.15 0.250 -.4121968 .1073227
EV[]#rcs():2 | -.0205143 .1182018 -0.17 0.862 -.2521856 .211157
EV[]#rcs():3 | .0474099 .1849679 0.26 0.798 -.3151205 .4099402
_cons | -5.395894 .5640025 -9.57 0.000 -6.501318 -4.290469
log(gamma) | .3872794 .1669175 2.32 0.020 .0601271 .7144317
-------------+----------------------------------------------------------------
logb: |
time | .1695192 .0132021 12.84 0.000 .1436436 .1953947
time#M2[id] | 1 . . . . .
M1[id] | 1 . . . . .
_cons | .5005547 .0596044 8.40 0.000 .3837322 .6173771
sd(resid.) | .3454686 .0066128 .3327479 .3586755
-------------+----------------------------------------------------------------
id: |
sd(M1) | 1.026383 .0433163 .9449014 1.114892
sd(M2) | .1812489 .0124586 .1584039 .2073886
------------------------------------------------------------------------------
```

where `rcs(stime, df(3) log event)`

creates a restricted cubic spline of log time, with 3 degrees of freedom, with internal knots based on centiles of the observed event distribution. Phew.

As always, there’s a lot more `merlin`

can do, but hopefully this serves as a taster to modelling non-proportional hazards within a joint model framework.

Videos

### State-of-the-art statistical models for modern HTA

Videos

### Multilevel (hierarchical) survival models: Estimation, prediction, interpretation

Statistical Primers