Nonlinear joint models for individual dynamic prediction of risk of death using Hamiltonian Monte Carlo: application to metastatic prostate cancer
 Solène Desmée^{1}Email author,
 France Mentré^{1},
 Christine VeyratFollet^{2},
 Bernard Sébastien^{3} and
 Jérémie Guedj^{1}
DOI: 10.1186/s1287401703829
© The Author(s) 2017
Received: 8 February 2017
Accepted: 30 June 2017
Published: 17 July 2017
Abstract
Background
Joint models of longitudinal and timetoevent data are increasingly used to perform individual dynamic prediction of a risk of event. However the difficulty to perform inference in nonlinear models and to calculate the distribution of individual parameters has long limited this approach to linear mixedeffect models for the longitudinal part. Here we use a Bayesian algorithm and a nonlinear joint model to calculate individual dynamic predictions. We apply this approach to predict the risk of death in metastatic castrationresistant prostate cancer (mCRPC) patients with frequent ProstateSpecific Antigen (PSA) measurements.
Methods
A joint model is built using a large population of 400 mCRPC patients where PSA kinetics is described by a biexponential function and the hazard function is a PSAdependent function. Using Hamiltonian Monte Carlo algorithm implemented in Stan software and the estimated population parameters in this population as priors, the a posteriori distribution of the hazard function is computed for a new patient knowing his PSA measurements until a given landmark time. Timedependent area under the ROC curve (AUC) and Brier score are derived to assess discrimination and calibration of the model predictions, first on 200 simulated patients and then on 196 real patients that are not included to build the model.
Results
Satisfying coverage probabilities of Monte Carlo prediction intervals are obtained for longitudinal and hazard functions. Individual dynamic predictions provide good predictive performances for landmark times larger than 12 months and horizon time of up to 18 months for both simulated and real data.
Conclusions
As nonlinear joint models can characterize the kinetics of biomarkers and their link with a timetoevent, this approach could be useful to improve patient’s followup and the early detection of most at risk patients.
Keywords
Calibration Discrimination Hamiltonian Monte Carlo Individual dynamic prediction Nonlinear joint modelBackground
One potential major application of mathematical modelling to personalized medicine is to provide dynamic prediction of a disease progression and occurrence of a severe clinical event. For that purpose an increasingly popular approach in the statistical community is to use joint models, which simultaneously handle longitudinal and timetoevent data by maximizing the joint likelihood of both processes to avoid the bias due to informative dropouts or measurements error [1–3]. In this approach, the longitudinal part is described by a mixedeffect model and the survival part is described by a parametric or semiparametric hazard function which depends on the true, unobserved, biomarker kinetics. After the parameters have been estimated, they can be used as priors to make “individual dynamic predictions” of the disease progression [4–6], i.e., prospective predictions that can be updated over the patient’s followup. Nevertheless the numerical difficulties have long limited the use of joint models, and hence of dynamic predictions, to linear models for the longitudinal processes. Although this approach can be made more flexible by using splines [3, 7], it does not handle models that are nonlinear in the parameters, i.e., nonlinear mixedeffect models (NLMEM), such as mechanistic models defined by differential equations.
Recently, we and others have shown that the SAEM algorithm, initially developed for inference in NLMEM, could be used to provide precise estimates of nonlinear joint models [8, 9]. However, even after the population parameters have been estimated, one still needs to characterize the entire a posteriori distribution of the individual parameters, which has no analytical form, in order to perform individual dynamic predictions. In that context, Bayesian inference using Markov Chain Monte Carlo (MCMC) appears naturally suited. Nevertheless, traditional MCMC are based on random walk which provides estimators with good properties of convergence, but in practice and especially in a highdimensional context, this asymptotic behavior is of limited use because of finite computational resources. In that case, the Halmitonian Monte Carlo (HMC), implemented in Stan [10], uses the geometry of the parameters space to generate effective and rapid exploration of this space, and stronger guarantees on the convergence [11–14].
Here we propose to use nonlinear joint model and HMC to characterize the a posteriori distribution of the individual survival probability. We apply in this nonlinear context novel tools that have been developed to assess timedependent discrimination and calibration metrics of dynamic models, such as the Area under the Receiver Operating Characteristic (ROC) curve (AUC) and the Brier score (BS) [5, 15, 16].
The approach is applied to a phase III clinical trial in metastatic castrationresistant prostate cancer (mCRPC) patients [17] where prostatespecific antigen (PSA) is frequently measured. The priors are obtained using a nonlinear joint model whose parameters are estimated on a training dataset of 400 patients using the SAEM algorithm implemented in Monolix. We show how dynamic predictions can be performed by characterizing the full a posteriori distribution of the risk of death for new individuals using 200 simulated patients and 196 mCRPC patients from a validation dataset. The timedependent metrics for discrimination and calibration are obtained for different landmark and horizon times and we discuss the potential applications of this approach.
Methods
General framework
Joint model with a nonlinear longitudinal biomarker
where b(t, ψ _{ i }) is the true and unobserved value of the biomarker at time point t, nonlinear in regards to the individual parameters ψ _{ i }, and e _{ ij } is the residual Gaussian error of mean 0 and variance σ ^{2}. ψ _{ i }=g(μ,η _{ i }) depends on the fixed effects μ identical for all patients and on the random effects η _{ i } specific for each individual, and the function g defines the transformations of the individual parameters. The random effects are assumed to be normally distributed with mean zero and variancecovariance matrix Ω (\(\eta _{i}\sim \mathcal {N}(0,\ \Omega)\)), and are assumed independent of the residual error e _{ ij }.
the survival function conditionally on the random effects.
In a NLMEM framework, the Stochastic Approximation ExpectationMaximization (SAEM) algorithm [18] implemented in Monolix (www.lixoft.eu) provides unbiased estimates for both longitudinal and survival parameters [8, 9]. As in other EM algorithms, this algorithm is an iterative process where each iteration is divided into a step where the complete likelihood conditional on observations is calculated (Estep), and a step where the complete likelihood is maximized (Mstep). In addition, in the SAEM algorithm, the Estep is divided into two parts: a simulation of individual parameters using a Markov Chain Monte Carlo (MCMC) algorithm (Sstep), and then a calculation of the expected likelihood using a stochastic approximation (Astep). Once parameters are estimated, the Fisher information matrix (FIM) can be stochastically approximated to obtain the relative standard errors (r.s.e).
Individual dynamic predictions
 1.Draw a realization \(\eta ^{(\ell)}_{i}\) from the posterior distribution of the random effects:$${} {{\begin{aligned} \left\{\eta_{i}  \ X_{i}>s,\mathcal{Y}_{i}(s) ;\ \theta\right\} \sim \left\{\!\prod\limits_{j=1}^{n_{i}(s)} p(y_{ij}\ \ \eta_{i}; \theta)\!\right\}\! S_{i}\!\left(s\  \ g(\mu, \eta_{i}); \theta\right)p(\eta_{i}; \theta) \end{aligned}}} $$(4)
where n _{ i }(s) is the number of longitudinal measurements of patient i available at the landmark time s.
 2.
Infer \(\psi ^{(\ell)}_{i} = g\left (\mu, \eta ^{(\ell)}_{i} \right)\)
 3.
Compute: \(y^{(\ell)}_{i}(u)=b\left (u,\psi ^{(\ell)}_{i}\right)\) with u>0
And: \(\pi ^{(\ell)}_{i}\left (s+t\mathrel {\left \right.\kern \nulldelimiterspace }s\right)=\frac {S_{i}\left (s\mathrel {\left \right.\kern \nulldelimiterspace }\psi ^{(\ell)}_{i} ; \theta \right)S_{i}\left (s+t\mathrel {\left \right.\kern \nulldelimiterspace }\psi ^{(\ell)}_{i} ; \theta \right)}{S_{i}\left (s\mathrel {\left \right.\kern \nulldelimiterspace }\psi ^{(\ell)}_{i} ;\theta \right)}\)
The difficulty is the first step since the posterior distribution of the random effects has no analytical solution when the model for the biomarker is nonlinear in regard to the parameters. We use Hamiltonian Monte Carlo (HMC) algorithm implemented in Stan software version 2.8 [10] and its R (version 3.1.3) interface. The a priori distribution of the random effects is assumed to be normal of mean zero and variancecovariance matrix Ω, estimated on the training dataset. The a posteriori distribution of the individual random effects defined by the Eq. (4) requires the integration of the hazard function in Eq. (3) of the survival. For that purpose, we use a GaussLegendre quadrature of order 8. Of note, and unlike what was proposed by Rizopoulos (2011) [4] the uncertainty in θ is neglected (see “Results”).
And 95% prediction intervals for the individual predictions are obtained using the 2.5% and 97.5% Monte Carlo sample percentiles.
Model discrimination and calibration
The scaled Brier score is obtained using the estimated Brier score: \(\widehat {sBS}(s,t)=\frac {\widehat {BS}_{KM}(s,t){\widehat {BS}(s,t)}}{{\widehat {BS}}_{KM}(s,t)}\) where \({\widehat {BS}}_{KM}(s,t)\) is the Brier Score obtained using the KaplanMeier estimate at s+t in the training dataset.
Application to risk of death in patients with metastatic prostate cancer
Illustration focuses on metastatic CastrationResistant Prostate Cancer (mCRPC) for which PSA is frequently measured and survival is the primary endpoint. First we develop a reference nonlinear joint model on a training dataset, second we simulate mCRPC patients to evaluate dynamic predictions when the model is known and last we apply this approach to real mCRPC patients from a validation dataset.
Building a reference nonlinear joint model
Real data come from mCRPC patients of the control arm of a phase III clinical trial [17] that included 598 men treated with the firstline reference chemotherapy (docetaxel in combination with prednisone). All PSA measurements at baseline (i.e., measured within 8 days before treatment initiation) and after treatment initiation are used. Two patients that have no PSA measurements are not included in the analysis. PSA was measured every 3 weeks during treatment and every 12 weeks after treatment cessation, and the date of death or last visit was obtained for all patients. Data are randomly split into two datasets as described in [26]: a training dataset containing N=400 patients to develop the reference nonlinear joint model and to estimate the population parameters θ, and a validation dataset containing the N ^{′}=196 remaining patients to provide individual dynamic predictions and assess the predictability of the model on real data.
In order to describe the kinetics of PSA, we use the biexponential model already described in [9]. In brief, this model assumes that PSA is produced by prostatic cells and that chemotherapy inhibits prostatic cells proliferation with a constant effectiveness ε until the time T _{ esc }. This leads to the following analytical solution for PSA kinetics:
where P S A _{0} (ng.mL ^{−1}) is the PSA value at treatment initiation, δ (day ^{−1}) is the rate of PSA elimination, r (day ^{−1}) the rate of prostatic cells proliferation in absence of treatment and d (day ^{−1}) the rate of prostatic cells elimination. ε is the constant treatment effect and T _{ esc } the time at which treatment has no longer an effect.
Because only 4 parameters can be identified from Eq. (6), we fix d to 0.046 day ^{−1}, corresponding to a halflife of tumor cells of 15 days, consistent with an apoptotic index of 5% in metastatic prostate cancer [27]. Moreover we fix δ to 0.23 day ^{−1}, corresponding to a PSA halflife in blood of about 3 days [28]. Finally PSA kinetics is defined by the vector of parameters ψ={r, P S A _{0}, ε, T _{ esc }}.
Here in the NLMEM (Eq. (1)), the observed biomarker y _{ ij } corresponds to the j ^{ t h } measurement of log(P S A+1) for the patient i at time t _{ ij } and b(t,ψ) is log(P S A(t,ψ)+1), consistent transformation with an additive residual error. Lognormal distributions for r, P S A _{0} and T _{ esc } (i.e. ψ _{ i }= log(μ)+η _{ i }) and logitnormal distribution for ε (i.e. ψ _{ i }=l o g i t(μ)+η _{ i } with \(logit(x)=log\left (\frac {x}{1x}\right)\) for 0<x<1) are assumed. The variancecovariance matrix of the random effects is diagonal with parameters: \(\Omega = diag\left (\omega ^{2}_{r}, \omega ^{2}_{PSA_{0}}, \omega ^{2}_{\varepsilon }, \omega ^{2}_{T_{esc}}\right)\).

No link: f(t, ψ _{ i })=0

Current PSA: f(t, ψ _{ i })= log(P S A(t, ψ _{ i })+1)

Current PSA slope: \(f\left (t,\ \psi _{i}\right)=\frac {d\log (PSA\left (t,\ \psi _{i}\right)+1)}{dt}\)

Area under PSA: \(f\left (t,\ \psi _{i}\right)=\int _{0}^{t} \log (PSA\left (u,\ \psi _{i}\right)+1)du\)
BIC and parameters estimates (r.s.e (%)) of PSA kinetics and survival in the N=400 patients for the 4 joint models
Models  No link  Current PSA  PSA slope  Area under PSA 

BIC  14350  14192  14291  14327 
r (day^{−1})  0.054 (1)  0.054 (1)  0.055 (1)  0.054 (1) 
P S A _{0}(ng.mL^{−1})  74.6 (8)  73.9 (8)  73.4 (8)  74.9 (8) 
ε  0.35 (5)  0.34 (5)  0.35 (5)  0.35 (5) 
T _{ esc } (day)  138 (4)  138 (4)  142 (4)  136 (4) 
λ (day)  885 (4)  3800 (9)  1500 (9)  1410 (13) 
k  1.52 (3)  1.19 (1)  1.33 (9)  1.15 (7) 
β    0.32 (4)  100 (10)  0.00025 (20) 
ω _{ r }  0.098 (5)  0.098 (4)  0.11 (5)  0.10 (5) 
\(\omega _{PSA_{0}}\)  1.57 (4)  1.57 (4)  1.55 (4)  1.56 (4) 
ω _{ ε }  1.35 (5)  1.34 (5)  1.22 (5)  1.36 (5) 
\(\omega _{T_{esc}}\)  0.68 (5)  0.64 (5)  0.63 (5)  0.66 (5) 
σ  0.38 (1)  0.38 (1)  0.38 (1)  0.38 (1) 
Simulation
PSA and timetodeath are simulated for N _{ sim }=200 patients using the vector θ estimated with the reference nonlinear joint model. PSA was measured every 3 weeks for 30 months or until the simulated timetodeath, if it occurs before. No other mechanism than death is considered for dropout. Dynamic individual predictions are performed using Stan and R (codes available as Additional file 2) as explained previously in the general framework with L=200 the number of Monte Carlo samples. For each landmark s∈{0,6,12,18} months and each horizon time t>2 months, we calculate the coverage probabilities of the 95% prediction intervals for both PSA and hazard rate, i.e., the proportion of simulated patients for whom the true value of interest (either simulated PSA value or the simulated risk of death) is contained in the corresponding Monte Carlo 95% prediction interval.
Then, using the package timeROC of R [25] we estimate timedependent AUC, BS and sBS for each landmark time s and horizon time t, using the KaplanMeier estimates computed in the N _{ sim }=200 simulated patients themselves to calculate B S _{ KM }(s,t).
Real data
Individual dynamic predictions are calculated following the same approach than in the simulation in the N ^{′}=196 mCRPC patients of the validation dataset using the reference nonlinear joint model and L=200. Likewise, timedependent AUC, BS and sBS are estimated using the N ^{′}=196 real patients. For B S _{ KM }(s,t) the KaplanMeier estimates in the training dataset are used.
Results
Reference nonlinear joint model
All PSA measurements and vital status in the N=400 patients from the treatment initiation to the end of followup are used to estimate the population parameters of the 4 proposed joint models. Overall 5 710 PSA measurements are used with median [minimum; maximum] number of measurements per patient of 13 [1 ; 55]. Regarding survival, 286 deaths occur (71.5%), leading to a median survival [KaplanMeier 95% confidence interval] of 656 days [598 ; 741].
All parameters are precisely estimated with relative standard errors (r.s.e) smaller than 9% for both fixed effects and variance components (Table 1), and thus the uncertainty on θ is neglected in the following. Details of the model evaluation (individual fits, residuals for longitudinal and survival processes and mean survival curves for both training and validation datasets) are provided in the Additional file 1.
Simulated data
Real data
In the validation dataset, a total of 2,720 PSA measurements are collected and the median [min ; max] number of measurements per patient is 13 [2 ; 57]. 145 deaths occur (74.0%), with a median survival time [KaplanMeier 95% confidence interval] of 598 days [547 ; 732].
Discussion
This work is the first one to our knowledge to perform individual dynamic predictions in nonlinear joint models. The following approach is used: i) the priors for the parameters are found by estimating population parameters in a large training dataset using the SAEM algorithm [8, 9], ii) the distribution of the individual parameters is found using the Hamiltonian Monte Carlo (HMC) algorithm and prediction interval for the risk of death is derived accordingly, iii) the predictive performances are assessed using timedependent discrimination and calibration metrics previously developed in a context of linear model.
Here we use HMC implemented in Stan to characterize the full a posteriori distribution of the individual random effects. Of note softwares for nonlinear mixedeffects models (R, SAS or more specifically Monolix or Nonmem in pharmacometrics) can also produce individual “posthoc” parameters, typically the mode (or the mean) of the conditional distribution of the random effects. Yet, in clinical practice, having only the most likely value of the prediction does not account for the uncertainty on the individual parameter estimates. In order to characterize the prediction interval, one frequent approach is to use asymptotic Gaussian approximations [29]. However this may not always be accurate, for instance when the data are limited and additionally it does not take into account the correlations between the parameters. We show by simulation that using HMC implemented in Stan provides good coverage probabilities, except for long followup where the prediction interval tended to be overoptimistic (s=18 months, see Fig. 1). Whether this is specific to this simulation framework or is a more general pattern will need to be verified. Likewise a formal comparison between HMC and traditional MCMC methods in context of individual dynamic prediction using nonlinear joint model could be of interest. In terms of model prediction assessment, the AUC and BS metrics are modelfree and thus can be applied to a nonlinear context using existing packages [25]. Here, while the AUC and the BS improve over the landmark time in the simulation study, they tend to stagnate in the real data. This is likely due to the fact that in the simulation the amount of data increases linearly with the landmark time (since we assume measurements every 3 weeks), while in the real data PSA measurements become less frequent over time in patients after the end of treatment.
Our model framework contains several limitations. First the training and the validation dataset come from the same clinical trial. Second dynamic predictions are performed neglecting the uncertainty on the population parameters θ. In our context this approximation is reasonable because the training dataset is large and the relative standard errors are small compared to the betweenpatient variability (Table 1). However in other contexts this may not be true, for instance in the first steps of adaptive schemes where each new individual is used to update the model prediction. In this case or in external validation, a full Bayesian approach, that can also be done with Stan [30], could be relevant. Further, the biological model, albeit nonlinear, remains very simplistic. For instance effect of covariates like age could be investigated on the longitudinal process. Moreover the model does not account for the mechanisms leading to resistance and then relapse to treatment that we identified previously [26]. Rather we assume that PSA kinetics and risk of death are not modified after treatment cessation and continue at the same pace than before. Moreover PSA kinetics only was assumed to drive the complex process leading to death. These simplifications may explain in part why the model is good at identifying patients at higher risk but does less well at predicting the timetodeath.
Conclusion
Beside the concrete application shown here, we believe that this approach can be exemplified to develop more biologically relevant models in various medical context. In that respect the recent release of Stan software for stiff ODEs will make possible to use more mechanistic joint models naturally integrating the correlation between several longitudinal processes. Thus, the development of nonlinear models that will accompany the collection of new biomarkers in routine [31] may be an important step towards a better prediction of the risk of death and improve the early identification of patients at greater risk.
Declarations
Acknowledgements
The authors thank MarieKarelle Rivière for her help in Stan and Hervé Le Nagard and François Cohen for the use of the computer cluster services hosted on the “Centre de biomodélisation UMR1137”.
Funding
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.
Availability of data and materials
Simulated data are provided in the Additional file 2. Researchers may access the real anonymous data by sending a formal request to the authors.
Ethical approval and consent to participate
Ethical approval from Institutional Review Board (IRB) and Independent Ethics Committee (IEC) was obtained for this clinical trial, which allowed the use of data recorded in this database. Each patient was informed that medical data can be use in medical research.
Authors’ contributions
SD, FM and JG designed the study and performed the analysis. All authors drafted the manuscript. All authors read and approved the final manuscript.
Consent for publication
Yes.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Tsiatis AA, Davidian M. Joint modeling of longitudinal and timetoevent data: an overview. Statistica Sinica. 2004; 14:809–34.Google Scholar
 Wu L, Liu W, Yi GY, Huang Y. Analysis of longitudinal and survival data: joint modeling, inference methods, and issues. J Probab Stat. 2011;2012;17.Google Scholar
 Rizopoulos D. Joint Models for Longitudinal and TimetoEvent Data: With Applications in R.Boca Raton: CRC Press; 2012.View ArticleGoogle Scholar
 Rizopoulos D. Dynamic predictions and prospective accuracy in joint models for longitudinal and timetoevent data. Biometrics. 2011; 67(3):819–29.View ArticlePubMedGoogle Scholar
 Rizopoulos D, Murawska M, Andrinopoulou ER, Molenberghs G, Takkenberg JJM, Lesaffre E. Dynamic predictions with timedependent covariates in survival analysis using joint modeling and landmarking. arXiv preprint arXiv:1306.6479. 2013. https://arxiv.org/abs/1306.6479.
 ProustLima C, Taylor JM. Development and validation of a dynamic prognostic tool for prostate cancer recurrence using repeated measures of posttreatment PSA: a joint modeling approach. Biostatistics. 2009; 10:535–49.View ArticlePubMedPubMed CentralGoogle Scholar
 Brown ER, Ibrahim JG, DeGruttola V. A Flexible BSpline Model for Multiple Longitudinal Biomarkers and Survival. Biometrics. 2005; 61(1):64–73.View ArticlePubMedGoogle Scholar
 Mbogning C, Bleakley K, Lavielle M. Joint modelling of longitudinal and repeated timetoevent data using nonlinear mixedeffects models and the stochastic approximation expectation maximization algorithm. J Stat Comput Simul. 2015; 85(8):1512–28.View ArticleGoogle Scholar
 Desmée S, Mentré F, VeyratFollet C, Guedj J. Nonlinear mixedeffect models for prostatespecific antigen kinetics and link with survival in the context of metastatic prostate cancer: a comparison by simulation of twostage and joint approaches. AAPS J. 2015; 17(3):691–9.View ArticlePubMedPubMed CentralGoogle Scholar
 Stan Development Team. Stan Modeling Language User’s Guide and Reference Manual, Version 2.8.0. 2015. http://mcstan.org/users/documentation/index.html.
 Neal RM, et al.MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo. 2011; 2:113–62.Google Scholar
 Hoffman MD, Gelman A. The NoUturn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J Mach Learn Res. 2014; 15(1):1593–623.Google Scholar
 Betancourt M, Girolami M. Hamiltonian Monte Carlo for hierarchical models. Current trends in Bayesian methodology with applications. 2015; 79:30.Google Scholar
 Betancourt M. A Conceptual Introduction to Hamiltonian Monte Carlo. arXiv preprint arXiv:1701.02434. 2017. https://arxiv.org/abs/1701.02434.
 Blanche P, ProustLima C, Loubère L, Berr C, Dartigues JF, JacqminGadda H. Quantifying and comparing dynamic predictive accuracy of joint models for longitudinal marker and timetoevent in presence of censoring and competing risks. Biometrics. 2015; 71(1):102–13.View ArticlePubMedGoogle Scholar
 Mauguen A, Rachet B, MathoulinPélissier S, Lawrence GM, Siesling S, MacGrogan G, et al.Validation of death prediction after breast cancer relapses using joint models. BMC Med Res Methodol. 2015; 15(1):1.View ArticleGoogle Scholar
 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 castrationresistant prostate cancer (VENICE): a phase 3, doubleblind randomised trial. Lancet Oncol. 2013; 14(8):760–8.View ArticlePubMedGoogle Scholar
 Delyon B, Lavielle M, Moulines E. Convergence of a stochastic approximation version of the EM algorithm. Ann Stat. 1999; 27:94–128.View ArticleGoogle Scholar
 Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005; 61(1):92–105.View ArticlePubMedGoogle Scholar
 Parast L, Cheng SC, Cai T. Landmark prediction of longterm survival incorporating shortterm event time information. J Am Stat Assoc. 2012; 107(500):1492–501.View ArticlePubMedPubMed CentralGoogle Scholar
 Schoop R, Graf E, Schumacher M. Quantifying the predictive performance of prognostic models for censored survival data with timedependent covariates. Biometrics. 2008; 64(2):603–10.View ArticlePubMedGoogle Scholar
 Graf E, Schmoor C, Sauerbrei W, Schumacher M. Assessment and comparison of prognostic classification schemes for survival data. Stat Med. 1999; 18(1718):2529–45.View ArticlePubMedGoogle Scholar
 Blanche P, Latouche A, Viallon V. Timedependent AUC with rightcensored data: a survey. In: Risk Assessment and Evaluation of Predictions. New York: Springer: 2013. p. 239–51.Google Scholar
 Gerds TA, Schumacher M. Consistent estimation of the expected brier score in general survival models with rightcensored event times. Biometrical J. 2006; 48(6):1029–40.View ArticleGoogle Scholar
 Blanche P. Package ’timeROC’: timedependent ROC curve and AUC for censored survival data. Vienna: R Foundation for Statistical Computing; 2013.Google Scholar
 Desmée S, Mentré F, VeyratFollet C, Sébastien B, Guedj J. Using the SAEM algorithm for mechanistic joint models characterizing the relationship between nonlinear PSA kinetics and survival in prostate cancer patients. Biometrics. 2016; 73(1):305–12.View ArticlePubMedGoogle Scholar
 Tu H, Jacobs SC, Borkowski A, Kyprianou N. Incidence of apoptosis and cell proliferation in prostate cancer: Relationship with TGF β1 and bcl2 expression. Int J Cancer. 1996; 69(5):357–63.View ArticlePubMedGoogle Scholar
 Polascik TJ, Oesterling JE, Partin AW. Prostate specific antigen: a decade of discoverywhat we have learned and where we are going. J Urol. 1999; 162(2):293–306.View ArticlePubMedGoogle Scholar
 Prague M, Commenges D, Drylewicz J, Thiébaut R. Treatment Monitoring of HIVinfected patients based on mechanistic models. Biometrics. 2012; 68(3):902–11.View ArticlePubMedGoogle Scholar
 Brilleman SL, Crowther MJ, May MT, Gompels M, Abrams KR. Joint longitudinal hurdle and timetoevent models: an application related to viral load and duration of the first treatment regimen in patients with HIV initiating therapy. Stat Med. 2016; 35(20):3583–94.View ArticlePubMedGoogle Scholar
 Wilbaux M, Tod M, De Bono J, Lorente D, Mateo J, Freyer G, et al.A joint model for the kinetics of CTC count and PSA concentration during treatment in metastatic castrationresistant prostate cancer. CPT: Pharmacometrics Syst Pharmacol. 2015; 4(5):277–85.Google Scholar