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.