Performance of joint modelling of time-to-event data with time-dependent predictors: an assessment based on transition to psychosis data

Joint modelling has emerged to be a potential tool to analyse data with a time-to-event outcome and longitudinal measurements collected over a series of time points. Joint modelling involves the simultaneous modelling of the two components, namely the time-to-event component and the longitudinal component. The main challenges of joint modelling are the mathematical and computational complexity. Recent advances in joint modelling have seen the emergence of several software packages which have implemented some of the computational requirements to run joint models. These packages have opened the door for more routine use of joint modelling. Through simulations and real data based on transition to psychosis research, we compared joint model analysis of time-to-event outcome with the conventional Cox regression analysis. We also compared a number of packages for fitting joint models. Our results suggest that joint modelling do have advantages over conventional analysis despite its potential complexity. Our results also suggest that the results of analyses may depend on how the methodology is implemented.


103
h i (t) = h 0 (t)exp( T ω i + αm i (t)), i = 1, 2, …, n, t > 0. (1) 104 105 Here h 0 (t) denotes the baseline hazard rate, ω i is a vector of baseline predictors (e.g. treatment 106 indicator, gender, age, etc.) and  is the corresponding vector of regression coefficients. The 107 time dependent predictor is represented by m i (t) with α being the corresponding coefficient 108 vector. A commonly used model for m i (t) is the linear mixed-effects model. Specifically, m i (t) is 109 given by:  The JM package is very versatile and allows many variations to the fitting of joint models. 135 Firstly, it allows the baseline hazard to be unspecified, to take the form of the hazard 136 corresponding to the Weibull distribution for the event times or to be approximated by (user-137 controlled) piecewise-constant functions or splines. For ordinary Cox regression, the baseline 138 hazard is usually left unspecified. This is of course a well-known advantage of Cox regression. 139 This advantage avoids the restriction resulting from specifying a certain form for the baseline 140 hazard and at the same time still can offer valid statistical inference through the use of partial 141 likelihood. However, in the context of joint modelling, this advantage no longer holds because a 142 completely unspecified baseline hazard will generally lead to underestimation of the standard 185 (iii) a 0 and a 1 were the fixed effects with given values.
186 (iv) b 0i and b 1i were the random effects generated from a bivariate normal distribution with mean 187 0 and a given covariance matrix. 188 (v)  i (t) was the random error generated from a normal distribution with mean 0 and a given 189 variance.
190 The time-to-event data was generated as follows: 191 (i) The hazard rate, h i (t), for subject i at time t (t = 0, 1, 2, . . ., 364), was computed as follows: 204 (iv)If T i  C i , the survival status for subject i was taken to be 1 and the time to event occurrence 205 was taken to be T i . Otherwise, the survival status was taken to be 0 and the censoring time 206 was taken to be C i .

207
To complete the generation of the simulated data, the data collection of the time-dependent 208 predictor was taken to occur at regular time points, specifically at day 0 (i.e. baseline) and then at 209 30-day intervals thereafter. Also, the data of the time-dependent predictor was taken to be 210 unavailable after the event time or censoring time, whichever was applicable. Therefore, for each 211 subject, non-missing data for the time-dependent predictor were taken to be those at days 0, 30, 212 60 and so on up to the measurement occasion prior to the event time or censoring time. Any 213 post-event or post-censoring data were not used.

214
The parameters associated with the simulations were given the following values: The fixed-effect intercept, a 0 , was given the value 40. splines. This is an alternative semi-parametric approach following the same rationale as using 294 a piecewise-constant function.

295
The JM package offers two options for numerical integration: the standard Gauss-Hermite 296 rule and the pseudo-adaptive Gauss-Hermite rule. It has been shown that the latter can be more  Manuscript to be reviewed 329 be around 50 (as illustrated in the boxplots of Figures 1a and 1c). Table 2 shows these results for 330 the estimation of  1 in each set of simulation. For the confidence intervals, it is expected that 331 good performance should correspond to a coverage of approximately 95%, say 90% or more. For 332 unbiasedness, as mentioned above, the percentage of estimates less than the true parameter value 333 should be approximately 50 for good performance, say between 40 and 60. The shaded entries in 334 Table 2 are those scenarios which did not perform well. It can be seen that joineR and stjm 335 tended to show better results. Note also that, for a small number of the simulated datasets, the 336 estimates were not available when JM was used with a piecewise-constant baseline hazard due to 337 convergence problems.

338
The results in Table 2  355  1 is negative, this suggests that these two analysis methods tended to underestimate or even 356 reverse the value of  1 . Conversely, for the two JM analyses, a substantial number of datasets 357 had the percentage above 60 with some way above 60. This suggests that these two analysis 358 methods tended to overestimate the value of  1 .
359 Figure 3 shows the corresponding results for the estimation of. It can be seen from Figure   360 3a that all analysis methods had 90% or more for the confidence interval coverage. Figure 3b 361 shows that, except for a few occasions, the percentage of estimates less than the true parameter 362 value were all between 40 and 60 for all the analysis methods. These results suggest that the 363 performance of the different analysis methods were all good for the estimation of the group 364 effect.

365
Recall that all of the 32 sets of simulations were for a zero group effect. Simulations were 366 repeated on four of the simulation sets for a non-zero group effect. These four sets were sets 3, 7, 367 19 and 23 in Table 2. These sets were chosen because they showed the worst results for the two 368 Cox analyses and the two JM analyses. The reason that a non-zero group effect was not applied 369 to all the 32 simulation sets was that it was very time consuming to run the simulations. The non-370 zero group effect was taken to be -0.5. Table 3 shows the results of these simulations, which are 371 very similar to the corresponding results in Table 2. 373 As a further assessment of the joint modelling methodology and the various software packages, 374 the same analyses described above were applied to a set of real data collected in a study of risk For the estimation of , while there were also substantial differences between estimates, the 405 standard errors were similar. It is again of interest to note that there was considerable variation 406 among the p-values, although all of them were non-significant at the 0.05 level.