Extending the code in the open-source saemix package to fit joint models of longitudinal and time-to-event data

Background and objective: Joint modeling of longitudinal and time-to-event data has gained attention over recent years with extensive developments including nonlinear models for longitudinal outcomes and flexible time-to-event models for survival outcomes, possibly involving competing risks. However, in popular software such as R, the function used to describe the biomarker dynamic is mainly linear in the parameters, and the survival submodel relies on pre-implemented functions (exponential, Weibull, ...). The objective of this work is to extend the code from the saemix package (version 3.1 on CRAN) to fit parametric joint models where longitudinal submodels are not necessary linear in their parameters, with full user control over the model function.

Methods: We used the saemix package, designed to fit nonlinear mixed-effects models (NLMEM) through the Stochastic Approximation Expectation Maximization (SAEM) algorithm, and extended the main functions to joint model estimation. To compute standard errors (SE) of parameter estimates, we implemented a recently developed stochastic algorithm. A simulation study was proposed to assess (i) the performances of parameter estimation, (ii) the SE computation and (iii) the type I error when testing independence between the two submodels. Four joint models were considered in the simulation study, combining a linear or nonlinear mixed-effects model for the longitudinal submodel, with a single terminal event or a competing risk model.

Results: For all simulation scenarios, parameters were precisely and accurately estimated with low bias and uncertainty. For complex joint models (with NLMEM), increasing the number of chains of the algorithm was necessary to reduce bias, but earlier censoring in the competing risk scenario still challenged the estimation. The empirical SE of parameters obtained over all simulations were very close to those computed with the stochastic algorithm. For more complex joint models (involving NLMEM), some estimates of random effects variances had higher uncertainty and their SE were moderately under-estimated. Finally, type I error was controlled for each joint model.

Conclusions: saemix is a flexible open-source package and we adapted it to fit complex parametric joint models that may not be estimated using standard tools. Code and examples to help users get started are freely available on Github.