Alexandra Lavalley-Morelle(1), France Mentré(1,2), Jimmy Mullaert(1,3), Emmanuelle Comets(1,4)
(1) Université Paris Cité, INSERM, IAME, F-75018 Paris, France. (2) Department of Epidemiology, Biostatistics and Clinical Research, AP-HP, Bichat-Claude Bernard University Hospital, F-75018 Paris, France. (3) Université Paris-Saclay, UVSQ, Institut Curie, Cancer et Génome, 92210 Saint-Cloud, France. (4) Université de Rennes, Inserm, EHESP, Irset - UMRS 1085, F-35000 Rennes, France.
Introduction: Joint modeling of longitudinal and time-to-event data has gained attention over recent years with extensive developments [1,2]. When the longitudinal outcome is described by a nonlinear mixed-effects model (NLMEM) and the survival outcome by a competing risk model, few developments have been proposed, all of which rely on pharmacometrics software such as Monolix and NONMEM [3,4]. On the other hand, general-purpose software such as R lack tools for nonlinear joint model estimation, particularly in the context of competing risks. Spreading the use of these models could be particularly beneficial for the scientific community. Therefore, the objective of this work is to extend the code from the saemix package (version 3.1 on CRAN) [5] to fit parametric joint models with full user control over the model function. We also aim to evaluate the performances of the estimation in a simulation study and apply it to a real dataset previously analyzed using Monolix software.
Methods: We used the saemix package, designed to fit nonlinear mixed-effects models (NLMEM) through the Stochastic Approximation Expectation Maximization (SAEM) algorithm [6], and extended the main functions to joint model estimation. Any parametric model and link function between longitudinal and survival parts can be fully written by the user. We also implemented a recently developed stochastic algorithm requiring only first derivatives [7] to compute standard errors (SE) of parameter estimates. A simulation study was proposed to assess (i) the relative bias (RB) and relative root mean square errors (RRMSE) of the estimated parameters, (ii) the accuracy of the estimated SE and (iii) the adequacy of 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 time-to-event or a competing risk model. We considered a natural link setting where the (predicted) longitudinal values are directly related to the survival process. For each joint model, we simulated 200 datasets of 100 patients. We assumed a rich design and parameters were chosen to obtain about 50% of events in single event models, and about 45% for each of the two events in competing risks models. Finally, a real case study in patients hospitalized for SARS-COV-2 infection recently published was considered as an illustration [8]. In this study, authors used Monolix software to estimate a multivariate joint model including (1) three linear and nonlinear mixed-effects models to describe biomarker dynamics (neutrohpils, C-reactive protein (CRP) and arterial pH) and (2) a competing risk model to describe the risk of in-hospital death and discharge from hospital. We used our saemix extension to fit this model in R.
Results: The new extension of saemix is available on Github [9]. Parameters were precisely and accurately estimated with low bias and uncertainty in all simulation scenarios. The link coefficients had respectively RB and RMMSE less than 7% and 25% for all joint models considered. Moreover, the type 1 error of the Wald test on the link coefficients was close to 5% for each scenario. The empirical SE of parameters obtained over all simulations were very close to those computed with the stochastic algorithm. For complex joint models (involving 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 (RB = 20% for one coefficient). In addition, some estimates of random effect variances had higher uncertainty (RRMSE greater than 50%) and their SE were moderately underestimated. Finally, in the real case study, we found similar results as in the previous analysis using Monolix software: an increase of the neutrophils/CRP and a decrease of the pH were significantly associated with a higher risk of death and lower risk of discharge.
Conclusions: saemix is a flexible open-source package and we adapted it to fit complex parametric joint models. This work extends the use of nonlinear joint models mainly available in specialized pharmacometrics software such as Monolix or NONMEM and makes therefore saemix the lone R package supporting estimation of nonlinear joint models with competing risks. Code and examples to help users get started are freely available on Github.
References:
[1] Desmée et al. Nonlinear Mixed-effect Models for Prostate-specific Antigen Kinetics and link with survival in the context of Metastatic Prostate Cancer, The AAPS Journal, 2015
[2] Kerioui et al. Bayesian inference using Hamiltonian Monte-Carlo algorithm for nonlinear joint modeling in the context of cancer immunotherapy, Statistics in Medicine, 2020
[3] Lavalley-Morelle et al. Joint modeling under competing risks: Application to survival prediction in patients admitted in Intensive Care Unit for sepsis with daily Sequential Organ Failure Assessment score assessment, CPT: Pharmacometrics & Systems Pharmacology, 2022
[4] Krishnan et al. Multistate model for Pharmacometric analyses of overall survival in HER2-negative breast cancer patients treated with docetaxel. CPT: Pharmacometrics & Systems Pharmacology, 2021
[5] Comets et al. Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm, Journal of Statistical Software, 2017
[6] Delyon et al. Convergence of a stochastic approximation version of EM algorithm, The Annals of Statistics, 1997
[7] Delattre and Kuhn, Estimating Fisher Information Matrix in Latent Variable Models based on the Score Function, https://doi.org/10.48550/arXiv.1909.06094
[8] Lavalley-Morelle et al. Multivariate joint model under competing risks to predict death of hospitalized patients for SARS-CoV-2 infection, Biometrical Journal, 2023
[9] https://github.com/saemixdevelopment/saemixextension/tree/master/joint
Reference: PAGE 32 (2024) Abstr 11060 [www.page-meeting.org/?abstract=11060]
Poster: Methodology - New Modelling Approaches