Joint modelling for nonlinear longitudinal PSA kinetics and survival data in metastatic prostate cancer patients
Solène Desmée (1,2), France Mentré (1,2), Christine Veyrat-Follet (3), Bernard Sébastien (4) and Jérémie Guedj (1,2)
(1) INSERM, IAME, UMR 1137, F-75018, Paris, France (2) Université Paris Diderot, IAME, UMR 1137, Sorbonne Paris Cité, F-75018, Paris, France (3) Drug disposition, Disposition Safety and Animal Research Department, Sanofi, Bridgewater, USA (4) Biostatistics and Programming, Sanofi, Chilly-Mazarin, France
Introduction
Since its introduction, serum prostate specific antigen (PSA) is massively used in the screening, the management and the response to treatment of prostate cancer. However there is no general consensus on the relevance of PSA for prediction of disease evolution and survival, in particular in patients with advanced cancer. The lack of consensus might be partly related to the use of some specific aspects of the PSA trajectory, such as doubling time or nadir value during treatment, and not of all the kinetic information available, in the evaluation of the PSA predictive ability. A more sophisticated approach is to use a joint model, which characterizes simultaneously the entire kinetics of a biomarker and its impact on a time-to-event [1,2].
Objectives
The objectives of this PhD thesis were the following:
1. To compare by simulation the precision of parameters obtained using joint and two-stage approaches when PSA kinetics is described by a nonlinear mixed-effect model (NLMEM) using the SAEM algorithm implemented in Monolix
2. To characterize on real data the relationship between PSA kinetics and survival in metastatic castration-resistant prostate cancer (mCRPC) patients treated by docetaxel
3. To evaluate the predictability on survival of Bayesian individual dynamic predictions of PSA kinetics
Methods
Data came from one arm of a phase 3 clinical trial [3]. In this study 596 mCRPC men were treated by docetaxel, the reference chemotherapy, and PSA had to be measured every 3 weeks during treatment, then every 12 weeks until the end of study, when vital status was collected. For the sake of internal validation, the dataset was randomly split into a training sample and a validation sample of 400 and 196 patients, respectively.
1. We performed a clinical trial simulation to compare joint and two-stage approaches. PSA kinetics was described by a biexponential model with parameter values inspired from the real data. 100 samples with 500 patients were simulated using the R software with PSA measurements every 3 weeks for 2 years. The survival model was a Weibull proportional hazard model with increasingly high levels of association between the current PSA kinetics and survival. Using the SAEM algorithm, parameters estimation using joint approach was compared in terms of bias, type 1 error and power to two-stage approach which consists in fitting the PSA data using a NLMEM then inserting the obtained individual longitudinal parameters into a survival model [4].
2. To analyze the training sample, we developed a mechanistic joint model. The PSA kinetic model relied on 3 ordinary differential equations (ODEs) and assumed that PSA is produced by two types of cells, namely treatment-sensitive cells and -resistant cells, on which docetaxel has no effect [5]. Several survival models relying on functions of ODE outputs, such as PSA kinetics or number of resistant cells, were compared using the Bayesian Information Criterion (BIC). Model evaluation was based on individual weighted residuals (IWRES) and on Cox-Snell and Martingale residuals for longitudinal and survival parts, respectively. Model prediction was assessed on the validation sample.
3. We focused on the possibility to predict survival in a new patient of the validation sample based on his individual PSA kinetics and assuming that the model and the population parameters are known and used as priors. Thus the PSA observations of the new patient were assumed to be available until a landmark time s (if the patient is still alive) and we aimed to predict the conditional probability of survival up to the prediction horizon s+t with t>0. To take into account the individual estimation uncertainty [6], for each patient of the validation sample and each s∈{0,3,6,12,18} months, we drew 200 Monte-Carlo samples of individual parameter using the STAN software [7] and computed 200 PSA trajectories and survival functions. This procedure provided individual dynamic predictions of PSA evolution and survival with 95% prediction intervals using the 2.5% and 97.5% simulated percentiles. Model discrimination and calibration were assessed using time-dependent area under the ROC curve (AUC), Brier Score (BS) and scaled BS (sBS) [8–10] along with Monte-Carlo confidence intervals.
Results
1. As described in detail in [11], we found that estimation with the correct joint models provided small biases regardless of the association between PSA and survival. Inversely the two-stage approach led to increasingly high level of bias when the association increased. In particular there was a systematic underestimation of the PSA effect on survival. Type 1 error, i.e., the probability to detect an effect of PSA on survival when there is none, was equal to 4 and 12% for joint and two-stage approaches, respectively. Power was 100% for the two methods.
2. The mechanistic nature of the model allowed us to consider other markers for survival that are not observed. Thus a joint model relying on the non-observed number of both resistant and sensitive cells led to the lowest BIC. No misspecification was revealed by residuals. The relevance of this model, instead of a model relying on the sole PSA, was reinforced by the fact that it could be used to correctly predict the survival curve of the validation sample using only PSA measurements.
3. This approach allowed to plot for each new patient individual dynamic predictions of PSA evolution and survival with prediction intervals which got wider when the horizon t increased, and shrank when s increased. Using the AUC, BS and sBS, we showed that joint modelling, provided that PSA kinetics is observed for at least 6 months (s>6), could identify the most-at-risk patients and precisely predict their survival in the next 12 months (i.e., t<12).
Conclusion
Our work investigated the use of mechanistic nonlinear joint models to predict the survival in patients with mCRPC treated by docetaxel. SAEM algorithm implemented in Monolix was shown by simulation to provide precise estimates for nonlinear joint models. When used on real data, a model accounting for the kinetics of both docetaxel-resistant and -sensitive cells provided a better fit to the data than a model relying on the sole PSA kinetics. Lastly we showed how this approach could be used in personalized medicine to prospectively predict patient’s survival using STAN to get the full conditional distribution and using new metrics for evaluation of predictability. Time-dependent discrimination and calibration metrics allowed to define predictive capacities of the joint model according to landmark time and prediction horizon. This work opens the way for the use of more complex and physiological joint models that could incorporate other relevant biomarkers to improve treatment evaluation and predictions in prostate cancer.
Acknowledgements:
The authors would like to thank the Drug Disposition Department, Sanofi, Paris, which supported Solène Desmée by a research grant during this work.
References:
[1] Tsiatis AA, Davidian M. Joint modeling of longitudinal and time-to-event data: an overview. Stat. Sin. 2004;14:809–34.
[2] Rizopoulos D. Joint Models for Longitudinal and Time-to-Event Data: With Applications in R. CRC Press; 2012.
[3] Tannock IF, Fizazi K, Ivanov S, Karlsson CT, Fléchon A, Skoneczna I, et al. Aflibercept versus placebo in combination with docetaxel and prednisone for treatment of men with metastatic castration-resistant prostate cancer (VENICE): a phase 3, double-blind randomised trial. Lancet Oncol. 2013;14:760–8.
[4] Ibrahim JG, Chu H, Chen LM. Basic Concepts and Methods for Joint Models of Longitudinal and Survival Data. J. Clin. Oncol. 2010;28:2796–801.
[5] Seruga B, Ocana A, Tannock IF. Drug resistance in metastatic castration-resistant prostate cancer. Nat. Rev. Clin. Oncol. 2011;8:12–23.
[6] Rizopoulos D. Dynamic Predictions and Prospective Accuracy in Joint Models for Longitudinal and Time-to-Event Data. Biometrics. 2011;67:819–29.
[7] Stan Development Team. Stan Modeling Language Users Guide and Reference Manual, Version 2.8.0. 2015.
[8] Graf E, Schmoor C, Sauerbrei W, Schumacher M. Assessment and comparison of prognostic classification schemes for survival data. Stat. Med. 1999;18:2529–45.
[9] Steyerberg EW, Vickers AJ, Cook NR, Gerds T, Gonen M, Obuchowski N, et al. Assessing the performance of prediction models: a framework for some traditional and novel measures. Epidemiol. Camb. Mass. 2010;21:128–38.
[10] Blanche P, Proust-Lima C, Loubère L, Berr C, Dartigues J-F, Jacqmin-Gadda H. Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and time-to-event in presence of censoring and competing risks. Biometrics. 2015;71:102–13.
[11] Desmée S, Mentré F, Veyrat-Follet C, Guedj J. Nonlinear Mixed-Effect Models for Prostate-Specific Antigen Kinetics and Link with Survival in the Context of Metastatic Prostate Cancer: a Comparison by Simulation of Two-Stage and Joint Approaches. AAPS J. 2015;17:691–9.
