# 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(t) = h_{0}(t) \exp(X_{1ij}\mathbf{\beta_{1}} + \alpha m_{i}(t))$

with our longitudinal submodel, $y_{}(t) = m_{i}(t) + \epsilon_{ij}$

where $m_{i}(t) = X_{2ij}(t)\mathbf{\beta}_{2} + Z_{ij}(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_{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,

. //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(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,

. //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

At @RedDoorAnalytics, we develop methodology and software for efficient modelling of biomarkers, measured repeatedly over time, jointly with survival outcomes, which are being increasingly used in cancer settings. We have also developed methods and software for general non-Markov multi-state survival analysis, allowing for the development of more plausible natural history models, where patient history can […]

Videos

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

Hierarchical time-to-event data is common across various research domains. In the medical field, for instance, patients are often nested within hospitals and regions, while in education, students are nested within schools. In these settings, the outcome is typically measured at the individual level, with covariates recorded at any level of the hierarchy. This hierarchical structure […]

Statistical Primers

### What are competing risks?

Competing risks In survival analysis, competing risks refer to the situation when an individual is at risk of experiencing an event that precludes the event under study to occur. Competing risks commonly occur in studies of cause-specific mortality, as all other causes of death than the one under study might happen before the individuals “have […]