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.

Multivariate joint model under competing risks to predict death of hospitalized patients for SARS-CoV-2 infection.

During the coronavirus disease 2019 (COVID-19) pandemic, several clinical prognostic scores have been proposed and evaluated in hospitalized patients, relying on variables available at admission. However, capturing data collected from the longitudinal follow-up of patients during hospitalization may improve prediction accuracy of a clinical outcome. To answer this question, 327 patients diagnosed with COVID-19 and hospitalized in an academic French hospital between January and July 2020 are included in the analysis. Up to 59 biomarkers were measured from the patient admission to the time to death or discharge from hospital. We consider a joint model with multiple linear or nonlinear mixed-effects models for biomarkers evolution, and a competing risks model involving subdistribution hazard functions for the risks of death and discharge. The links are modeled by shared random effects, and the selection of the biomarkers is mainly based on the significance of the link between the longitudinal and survival parts. Three biomarkers are retained: the blood neutrophil counts, the arterial pH, and the C-reactive protein. The predictive performances of the model are evaluated with the time-dependent area under the curve (AUC) for different landmark and horizon times, and compared with those obtained from a baseline model that considers only information available at admission. The joint modeling approach helps to improve predictions when sufficient information is available. For landmark 6 days and horizon of 30 days, we obtain AUC [95% CI] 0.73 [0.65, 0.81] and 0.81 [0.73, 0.89] for the baseline and joint model, respectively (p = 0.04). Statistical inference is validated through a simulation study.

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

Joint models of longitudinal process and time-to-event data have recently gained attention, notably to provide individualized dynamic predictions. In the presence of competing risks, models published mostly involve cause-specific hazard functions jointly estimated with a linear or generalized linear model. Here we propose to extend the modeling to full parametric joint estimation of a nonlinear mixed-effects model and a subdistribution hazard model. We apply this approach on 6046 patients admitted in intensive care unit (ICU) for sepsis with daily Sequential Organ Failure Assessment (SOFA) score measurements. The joint model is built on a randomly selected training set of two thirds of patients and links the current predicted SOFA measurement to the instantaneous risks of ICU death and discharge from ICU, both adjusted on the patient age. Stochastic Approximation Expectation Maximization algorithm in Monolix is used for estimation. SOFA evolution is significantly associated with both risks: 0.37, 95% confidence interval (CI) = [0.35, 0.39] for the risk of death and -0.38, 95% CI = [-0.39, -0.36] for the risk of discharge. A simulation study, inspired from the real data, shows the good estimation properties of the parameters. We assess on the validation set the added value of modeling the longitudinal SOFA follow-up for the prediction of death compared with a model that includes only SOFA at baseline. Time-dependent receiver operating characteristic area under the curve and Brier scores show that when enough longitudinal individual information is available, joint modeling provides better predictions. The methodology can easily be applied to other clinical applications because of the general form of the model.