Timedependent ROC curve analysis in medical research: current methods and applications
 Adina Najwa Kamarudin^{1}Email author,
 Trevor Cox^{1} and
 Ruwanthi KolamunnageDona^{1}
DOI: 10.1186/s1287401703326
© The Author(s). 2017
Received: 27 October 2016
Accepted: 28 March 2017
Published: 7 April 2017
Abstract
Background
ROC (receiver operating characteristic) curve analysis is well established for assessing how well a marker is capable of discriminating between individuals who experience disease onset and individuals who do not. The classical (standard) approach of ROC curve analysis considers event (disease) status and marker value for an individual as fixed over time, however in practice, both the disease status and marker value change over time. Individuals who are diseasefree earlier may develop the disease later due to longer study followup, and also their marker value may change from baseline during followup. Thus, an ROC curve as a function of time is more appropriate. However, many researchers still use the standard ROC curve approach to determine the marker capability ignoring the time dependency of the disease status or the marker.
Methods
We comprehensively review currently proposed methodologies of timedependent ROC curves which use single or longitudinal marker measurements, aiming to provide clarity in each methodology, identify software tools to carry out such analysis in practice and illustrate several applications of the methodology. We have also extended some methods to incorporate a longitudinal marker and illustrated the methodologies using a sequential dataset from the Mayo Clinic trial in primary biliary cirrhosis (PBC) of the liver.
Results
From our methodological review, we have identified 18 estimation methods of timedependent ROC curve analyses for censored event times and three other methods can only deal with noncensored event times. Despite the considerable numbers of estimation methods, applications of the methodology in clinical studies are still lacking.
Conclusions
The value of timedependent ROC curve methods has been reestablished. We have illustrated the methods in practice using currently available software and made some recommendations for future research.
Keywords
ROC curve Timedependent AUC Biomarker evaluation Eventtime Longitudinal data SoftwareBackground
In a screening process, an appropriate marker is used to provide information on the individual risk of disease onset. Information and signalling of future disease identification may be given by a single continuous measurement marker or a score. A single measurement could be any clinical measure such as cell percentage in the synthesis phase to detect breast cancer [1], CD4 cell counts to detect AIDS [2] or HIV1 RNA to detect HIV [3]. A score from a regression of potential factors or some other model to detect disease can also be used as a marker. Chambless and Diao [4] used the score from a logistic regression model, including several traditional and newer risk factors, to detect Coronary Heart Disease (CHD). Lambert and Chevret [5] used the prognostic score of four covariates (age, platelet count, prothrombin time, and serum alphafetoprotein level) to predict compensated cirrhosis patients’ survival and also used a score of three baseline characteristics (age, white blood cell and performance status) to predict eventfree survival (EFS) in acute leukaemia patients. Moreover, some studies used a published score as a marker in which the score considers the most important mortality predictors of a certain disease. For example, the Framingham risk score is used for cardiovascular patients [6] and the Karnofsky score is used for lung cancer patients [7].
The decision from a diagnostic test is often based on whether the marker value exceeds a threshold value, in which case the diagnosis for the individual is “diseased” and “nondiseased” otherwise. There is a possibility that the diagnostic test gives a positive result for a nondiseased individual or a negative result for a diseased individual. The sensitivity is defined as the probability of a diseased individual being predicted as having the disease (truepositive) and the specificity as the probability of a nondiseased individual being predicted as not having the disease (truenegative). These probabilities change as the threshold value for the marker changes and the value or range of threshold values chosen depends on the tradeoff that is acceptable between failing to detect disease and falsely identifying disease with the test [8]. In relation to this, the receiver operating characteristic (ROC) curve is a tool that simply describes the range of tradeoffs achieved by a diagnostic test. ROC curve analysis is extensively used in biomedical studies for evaluating the diagnostic accuracy of a continuous marker. It is a graphical display which plots sensitivity estimates (probability of a true positive) against one minus specificity (probability of a false positive) of a marker for all possible threshold values. The performance of a marker is evaluated by the area under the ROC curve (AUC) in which a higher AUC value indicates a better marker performance. The AUC is also equal to the probability of a diseased individual having a higher marker value than a healthy individual [8]. It is usually assumed that a higher marker value is more indicative of disease [8, 9] and we assume this for the rest of this article.
Recent research has incorporated time dependency in the sensitivity and specificity in disease (event)time data for individuals instead of using the standard ROC curve approach. Such methods are proven more effective; however, these methods are still underused in medical research. Once the timedependent setting has been applied, the disease status is observed at each time point which yields different values of sensitivity and specificity throughout the study.
with \( \left[1 Sp\left( c, t\right)\right]=\frac{\partial \left[1 Sp\left( c, t\right)\right]}{\partial c} dc. \)
The AUC is equal to the probability that the diagnostic test results from a randomly selected pair of diseased and nondiseased individuals are correctly ordered [10, 11].
Cumulative sensitivity and dynamic specificity (C/D)
At each time point t, each individual is classified as a case or control. A case is defined as any individual experiencing the event between baseline t = 0 and time t (individual A, B or E in Fig. 1a) and a control as an individual remaining eventfree at time t (individual C, D or E in Fig. 1a). The cases and controls are changing over time and each individual may play the role of control at the earlier time (when the event time is greater than the target time, i.e. T _{ i } > t) but then contributes as a case for later times (when the event time is less than or equal to the target time, i.e.T _{ i } ≤ t).
It is more appropriate to apply the C/D definitions when there is a specific time of interest that is used to discriminate between individuals experiencing the event and those eventfree prior to the specific time. This type of discrimination has more clinical relevance than the other definitions (I/D and I/S) and hence C/D definition has commonly been used by clinical applications [5, 13]. However, since some individuals may contribute as controls at an earlier time and then contribute as cases later, this definition uses redundant information in separating cases and controls [5].
Incident sensitivity and dynamic specificity (I/D)
A case for I/D definition is defined as an individual with an event at time t (individual A in Fig. 1a) while the control is an eventfree individual at time t. (individual C, D or F in Fig. 1a). In this definition, there are individuals neither a control nor case (when the event time is less than the target time, i.e. T _{ i } < t, individual B or E in Fig. 1a). Each individual who had an event may play the role of control at the earlier time (when the event time is greater than target time, i.e. T _{ i } > t) but then contributes as a case at the later incident time (when the event time is the same as the target time, i.e.T _{ i } = t).
Incident sensitivity and static specificity (I/S)
As illustrated in Fig. 1a, the controls are static and do not change (individuals D and F), and each individual only contributes once as a case or as eventfree individual within the fixed followup (0, t*).
The above definitions facilitate the use of standard regression approaches for characterizing sensitivity and specificity because the time prior to event T _{ i } − s _{ ik } can simply be used as a covariate.
Blanche et al. [13] have reviewed methodologies of timedependent ROC curve analysis under the C/D definition only; however, in this article, we have undertaken a comprehensive review of the current estimation methods under each definition and also identify additional methods, aiming to provide clarity for each methodology. We illustrate how each method is implemented on a timevarying disease status or over a time course of a longitudinal marker using a sequential dataset from Mayo Clinic trial in primary biliary cirrhosis (PBC) of the liver. We identify the software that implements current methods in practice and the need for further methodologies.
Benefits of timedependent ROC curve analysis
In the standard ROC curve analysis, the individual’s disease status is defined once a marker value is measured and it is assumed to be fixed for the whole study period. The study period is usually take a long followup and during this, the individual without disease earlier may have the disease. In contrast, the disease status of an individual is observed and updated at each time point in timedependent ROC curve analysis. With additional information of time of disease onset for each individual, a ROC curve can be constructed at several time points and the marker’s predictive ability can be compared. Thus, timedependent ROC curve is an efficient tool in measuring the performance of a candidate marker given the true disease status of individuals at certain time points. In general, a baseline marker value is used for computing the predictive ability but it can become weaker as the target time gets further from the baseline.
In longitudinal studies, the marker is measured several times within a fixed followup. If a marker measurement has ability to signify a pending change in the clinical status of a patient, then a timedependent ROC curve on a timevarying marker can be used to guide key medical decisions.
Challenges of timedependent ROC curve analysis
The most common problem is censoring, in which some individuals may be lost during the followup period. If the censored individuals are ignored, the estimation of the sensitivity and specificity may be biased as the information from the individual before censoring may contribute to the estimation. In a timedependent ROC curve analysis, the sensitivity and specificity are defined at each time point, where not all individuals are equally informative, and their contributions differing according to the aims and definitions used. A longitudinal biomarker brings an additional challenge to take account of the marker measurements at a number of visits for each individual. In the I/S definition, not all marker values are used but only the most recent, which is assumed more reliable for predicting the disease status [6, 14]. Other timedependent ROC curve approaches currently proposed for a longitudinal marker either assume noncensored eventtimes or ignore censored individual records.
Methods
We have used MEDLINE (Ovid), Scopus and the internet to search for relevant papers for our review. We restricted our search to English language published papers between years 1995 to August 2016 to ensure all methodology papers of timedependent ROC curves analysis were included. A total of 332 papers were found and 24 of these discussed timedependent ROC curve methodology. The remaining 308 papers included only an application of standard or timedependent ROC curves. For each methodology paper, the following details were extracted: definition of sensitivity and specificity (whether C/D, I/D, I/S or other), estimation method, type of estimation (nonparametric, semiparametric or parametric), limitations and availability of software. Only 16 methodology papers are within the scope of this review, and out of the 16 methodology papers, 10 (63%) discussed methodologies along the lines of the C/D definition. Three papers (19%) proposed methodologies based on the I/D definition, only one paper (6%) proposed methodology based on the I/S and another two papers (12%) proposed other methodologies for longitudinal marker measurements. Full details of the review are available as Additional file 1.
Summary of current methods for each definition
Definition and marker time  Sensitivity and specificity  Estimation method and R software (when available)  Pros/Cons  

C/D t = 0  \( \begin{array}{l} S{e}^C\left( c, t\right)= P\left({X}_i> c\Big{T}_i\le t\right)\\ {} S{p}^D\left( c, t\right)= P\left({X}_i\le c\Big{T}_i> t\right)\end{array} \)  CD1 survivalROC  Pro: Easy Cons: a) Produce nonmonotone sensitivity and specificity b) Not robust to markerdependent censoring  Pro: Clinically relevant since many clinical experiments aim to discriminate individuals with disease prior to specific time and healthy individual beyond that time Con: Use redundant information in separating cases and controls 
CD2 survivalROC  Pros: a) Produce monotone sensitivity and specificity b) Allow censoring to depend on marker Con: Does involve smoothing parameter  
CD3 Programme code  Pro: Does not involve any smoothing parameter Cons: a) Does involve recursive computation b) Produce nonmonotone specificity  
CD4 survAUC  Pros: a) Produce monotone sensitivity and specificity if the score is produced from a survival model b) Allow censoring to depend on marker Con: Not invariant to an increasing transformation of the marker  
CD5 timeROC  Pros: a) Produce monotone sensitivity and specificity b) More robust than CD1 and CD3 Con: Does not robust to marker dependent censoring  
CD6 timeROC  Pros: a) Produce monotone sensitivity and specificity b) Robust to markerdependent censoring c) More less biased than CD2 when censoring strongly depends on marker  
CD8, VL Cox survival  Pro: Straightforward to implement  
VL Aalen timereg  
VL KM prodlim  
C/D A longitudinal time point  \( \begin{array}{l} S{e}^C\left( c, t\right)= P\left({Y}_{i k}> c\Big{T}_i{s}_{i k}\le t\right)\\ {} S{p}^D\left( c,{t}^{*}\right)= P\left({Y}_{i k}\le c\Big{T}_i{s}_{i k}> t\right)\end{array} \)  AD4 (ECD2)  Pros: a) Produce monotone sensitivity and specificity b) Allow censoring to depend on marker Con: Does involve smoothing parameter  Pro: Use the most recent marker value prior to prediction time Con: Just use a marker value at a particular time instead of using all serial of marker value 
I/D t = 0  \( \begin{array}{l} S{e}^I\left( c, t\right)= P\left({X}_i> c\Big{T}_i= t\right)\\ {} S{p}^D\left( c, t\right)= P\left({X}_i\le c\Big{T}_i> t\right)\end{array} \)  ID1 risksetROC  Pros: Produce consistent sensitivity and specificity if the control set is unbiased  Pro: Allow timeaveraged summaries that directly relate to a familiar concordance measures such as Kendall’s tau or cindex Con: Require an exact time of interest which often just a few individual has an event at a particular point 
ID2  Pro: Potentially more robust than ID1 Con: Computationally intensive  
ID3 Programme code  Pros: a) Easier especially when involve a large number of marker b) Understandable since it is a “regressiontype” model  
I/S t = 0  \( \begin{array}{l} S{e}^I\left( c, t\right)= P\left({Y}_i> c\Big{T}_i= t\right)\\ {} S{p}^S\left( c,{t}^{*}\right)= P\left({Y}_i\le c\Big{T}_i>{t}^{*}\right)\end{array} \)  None  Pro: Allow separation of longterm survivors from healthy individual within a fixed followup Con: Require an exact time of interest which often just a few individual has an event at a particular point  
I/S A longitudinal time point  \( \begin{array}{l} S{e}^I\left( c, t\right)= P\left({Y}_{i k}> c\Big{T}_i{s}_{i k}= t\right)\\ {} S{p}^S\left( c,{t}^{*}\right)= P\left({Y}_{i k}\le c\Big{T}_i{s}_{i k}>{t}^{*}\right)\end{array} \)  IS1  Pro: Provides unbiased estimates of model parameters of sensitivity and specificity Con: Computationally intensive since involve spline functions  Pro: Use the most recent marker value prior to prediction time Con: Just use a most recent of marker value instead of all marker values 
IS2  Pro: Use the most recent marker value prior to prediction time Con: not a natural companion to hazard models  
Other All longitudinal time points  ROC(t, p) = S[a _{0}(T _{ ik }) + a _{1}(T _{ ik })S ^{− 1}(p)]  AD1 Programme code  Pro: Use all marker value along visit times in the estimation of ROC curve Con: Do not incorporate censored outcomes 
Naïve estimator of timedependent ROC curve analysis
Many studies have used an empirical estimator as a basis for comparison with other estimation methods. This estimator only considers observed events and, the sensitivity and specificity are calculated by the observed proportions of truepositives and truenegatives respectively.
Cumulative sensitivity and dynamic specificity (C/D)
Ten estimation methods have been proposed under C/D definition, and these are discussed in CD1 – CD8 below. CD8 describes three estimation methods.
(CD1) KaplanMeier estimator of Heagerty et al. [1]
However, this estimator yields nonmonotone sensitivity and specificity, and not bounded in [0, 1]. This problem is illustrated by the authors using a hypothetical dataset, and is due to the quadrant probability estimator \( \widehat{P}\left({X}_i> c,{T}_i> t\right)=\widehat{S}\left( t\Big{X}_i> c\right)\left(1{\widehat{F}}_X(c)\right) \) not necessarily producing a valid bivariate distribution as the redistribution to the right of the probability mass associated with censored observations will change as the conditioning set (X > c) changes. Another problem is that it is not robust to markerdependent censoring since the conditional KaplanMeier estimator, Ŝ(tX _{ i } > c), assumes the censoring process does not depend on the marker.
(CD2) nearest neighbour estimator of Heagerty et al. [1]
(CD3) KaplanMeier like estimator of Chambless and Diao [4]
An advantage of this method is the sensitivity is monotone and bounded in [0, 1]. A nice property of this nonparametric estimator is that it does not involve any smoothing parameter, unlike CD2. Chambless and Diao [4] have compared CD3 with the cstatistic gained from the logistic regression model of baseline values in a simulation study and apparently it shows little bias. In order to compute variances and confidence intervals of this estimator, Chambless and Diao [4] suggested using bootstrap resampling.
(CD4) alternative estimator of Chambless and Diao [4]
(CD5) Inverse probability of censoring weighting
(CD6) Conditional IPCW
(CD7) Weighted AUC (t)
(CD8) Viallon and Latouche [20] Estimators
VL Cox is very similar to the estimator proposed by Heagerty and Zheng [12] that will be introduced in method ID1 but it does not involve the computation of the bivariate expectation [20].
VL estimators are straightforward to implement since they just plugin the estimates of the conditional absolute risk function into the timedependent AUC estimator. This plugin nature allows their theoretical properties to follow the other established estimators of the conditional absolute risk function. Moreover, it is advisable to use CD8 compared to CD2 in the situations where the independence assumption between censoring time C, and the pair (T, Z) might be violated [20].
Incident sensitivity and dynamic specificity (I/D)
There are three estimation methods proposed under the I/D definition, these are discussed in ID1 – ID3 below.
Specific notation: Let R _{ i }(t) = I(Z _{ i } ≥ t) denote the atrisk indicator. Let ℛ_{ i }(t) = (i : R _{ i }(t) = 1) denote the individuals that are in the riskset at time t, which ℛ _{ t } ^{1} = (i; T _{ i } = t), are individuals with the event (cases) and ℛ _{ t } ^{0} = (i; T _{ i } > t) are individuals without the event (controls). Let n _{ t } = ℛ _{ t } ^{0}  be the size of the control set at time t and d _{ t } = ℛ _{ t } ^{1}  the size of case set at time t. Note that the riskset at time t can be represented as ℛ_{ t } = (ℛ _{ t } ^{1} ∪ ℛ _{ t } ^{0} ).
(ID1) Cox regression
 (i)
Fit a Cox model λ _{0}(t)exp(X _{ i } γ) where \( \gamma \) is the proportional hazard regression parameter. In order to relax the proportionality assumption, use a regressionsmoothing method to estimate the timevarying coefficient \( \widehat{\gamma}(t) \) and use it to estimate sensitivity in (ii) instead of y.
 (ii)The sensitivity can be evaluated using \( \widehat{\gamma}(t) \) from (i) as follows$$ \widehat{Se}\left( c, t\right)={\displaystyle {\sum}_i} I\left({X}_i> c\right){\pi}_k\left(\widehat{\gamma}(t), t\right). $$
 (iii)The specificity can be estimated empirically as follow$$ \widehat{Sp}\left( c, t\right)=1{\displaystyle {\sum}_k} I\left({X}_k> c\right)\frac{{\mathrm{\mathcal{R}}}_i^0(t)}{n_t}. $$
Heagerty and Zheng [12] suggested using flexible semiparametric methods such as locally weighted maximum partial likelihood (MPL) by Cai and Sun [25] as the regressionsmoothing method in (i), and simple local linear smoothing of the scaled Schoenfeld residuals [26] for reducing the bias [12].
The sensitivity is consistent for both the proportional and nonproportional hazards models whenever a consistent estimator of \( \widehat{\gamma}(t) \) is used [27]. Since the specificity is an empirical distribution function calculated over the control set, it is consistent provided the control set represents an unbiased sample [12]. It is suggested that the computation of standard errors and confidence intervals is carried out using the nonparametric bootstrap based on resampling of observations (X _{ i }, Z _{ i }, δ _{ i }) [12].
(ID2) weighted mean rank
where \( {K}_{h_n} \) is a standardized kernel function such that \( {\displaystyle {\sum}_j}{K}_{h_n}\left( t{t}_j\right)=1 \). Eq. (1) is a smoothed version of Eq. (2) and it is based on local Ustatistics summaries. SahaChaudhuri and Heagerty [28] suggested integrated mean square error (IMSE) as a potential method to select an optimal bandwidth.
Under certain conditions, SahaChaudhuri and Heagerty [28] showed that WMR(t) follows a normal distribution. It is suggested that this variance estimator for inference can be used in practice since it is simple and does not require resampling methods. Moreover, SahaChaudhuri and Heagerty [28] also provided the details of large sample properties of this estimator, and then the construction of a confidence interval for WMR(t) using the asymptotic properties is straightforward. Although it is desirable to obtain the simultaneous confidence bands for the function WMR(t), the theory may not be applicable in this case since the limiting process may not possess an independent increment structure. Instead, a simulation of a Gaussian process while keeping the estimates of ID2 fixed is needed to approximate the distribution of the Gaussian process and to estimate the quantiles. ID2 also has the advantage to be potentially robust since the relative bias remains significantly lower than for the ID1estimator.
(ID3) fractional polynomial
As the ID2 method is computationally intensive, especially in the selection of the bandwidth, Shen et al. [29] proposed a semiparametric timedependent AUC estimator which is easier and more applicable when comparing and screening a large number of candidate markers. The suggested model used fractional polynomials [30], the parameters of which are estimated by using a pseudo partiallikelihood function.
Maximizing this pseudo partiallikelihood yields parameter estimates \( \widehat{\beta} \). Then the timedependent AUC estimate is obtained from Eq. (3) as a smooth function of time t and β. In practice, the integrated AUC is always of interest for the I/D definition and it can be defined by \( {\displaystyle {\int}_0^{\tau}}\omega \left( t;\tau \right) A U C\left( t;\widehat{\beta}\right) d t \). When the weight function ω(t; τ) is invariant to time, the integrated AUC can be viewed as the global average of the AUC curve [29]. One major advantage of this estimator compared to ID2 is that the proposed method estimates the entire curve as a function of t and β while ID2 just uses a pointwise approach to estimate AUC. Further, this method is understandable and it is easier to make inference since it is a “regressiontype” method, with covariates being functions of time. In estimating the integrated AUC, the ID3 method is more convenient since it uses an analytical expression while ID2 computation is more complex since the kernelbased estimation procedure has to be repeated N times, and also the selection of bandwidths has to be considered. However, SahaChaudhuri and Heagerty [28] decreased the computational burden by calculating the integrated AUC as an average of AUC(t) at 10 time points, which can lead to approximation errors.
Incident sensitivity and static specificity (I/S)
There is only one estimation method proposed under the I/S definition found from the methodological review and one extended method which are discussed below.
(IS1) Marginal regression modelling approach
where g _{ D } and \( {g}_{\overline{D}} \) are specified inverse link functions, h _{0} and c _{0} are baseline functions of the threshold c that are completely unspecified. These nonparametric baseline functions of c represent the shape and location of the sensitivity and specificity functions while the parameters β _{0} and b _{0} quantify the covariate effects on them and η α _{0} and \( \xi \) α _{0} are the time effects. The dependence on time for sensitivity is through the parametric functions η α _{0}(t, s) = α ' _{0} η(t, s) and ξ α _{0}(s _{ ik }) = α ' _{0} ξ(s) where η and ξ are vectors of polynomial or spline basis functions.
and the corresponding score equation is solved to estimate the nonparametric baseline functions, Η _{0}(c). Further, θ _{0} is estimated by solving the integration of the corresponding score equation. Cai et al. [6] also proposed an approach that ignores censored observations.
Simulation studies [6] showed that the above method provides reasonably unbiased estimates of model parameters of sensitivity and specificity. The approach which includes the censored observations is always more precise than the one that excludes the censored observations.
(IS2) extended Cox regression
The main difference between I/D and I/S definitions is related to the controls. The controls in I/D are changing based on the target time whereas in I/S, controls are static survivors beyond a fixed time. This difference has motivated us to extend the Cox Regression method (ID1) to incorporate a longitudinally repeated marker using the I/S definition. A marker value at a particular visit time s is considered. Thus, we have changed the definition of the riskset as those individuals beyond target time by including those beyond a fixed followup. However, as I/S is not based on classification of the riskset at time t like I/D, this extended method cannot be said as a natural companion to hazard models. We have also extended the current software of ID1 (see Section for Software below) by redefining the riskset according to the I/S definition. The extended software can also be used with the baseline value of the marker.
Additional methods for longitudinal outcomes
Three estimation methods have been proposed for a longitudinal marker in addition to those described above under I/S definition, although some do not incorporate censoring. These estimation methods are discussed below. An extension of the C/D definition for a longitudinally repeated marker is suggested as a fourth method.
Specific notation: Let \( n={n}_D + {n}_{\overline{D}} \) denote the total number of individuals which is the summation of the where n _{ D } is the total number of cases and \( {n}_{\overline{D}} \) is the total number of controls. Let U _{ ik } ^{ T } = vec(T _{ i }, s _{ ik }) = U _{ D } denote the vector of covariates associated with U _{ ik }. The total number of longitudinally repeated marker values for cases is \( {N}_D={\displaystyle {\sum}_i^{n_D}{K}_i} \). The time prior to an event is defined as the time lag between the measurement time and the event time: T _{ ik } = T _{ i } − s _{ ik } as above. Similarly for controls, let Y _{ jl } be the biomarker value obtained from individual j at the l ^{ th } visit time s _{ jl } with \( j={n}_D+1, \dots,\ {n}_D+{n}_{\overline{D}} \) and l = 1, …, L _{ j }. Let \( {\boldsymbol{U}}_{jl}^T= v e c\left({s}_{jl}\right)={\boldsymbol{U}}_{\overline{D}} \) denote the vector of covariates associated with Y _{ jl }. The total number of longitudinally repeated marker values for controls is \( {N}_{\overline{D}}={\displaystyle {\sum}_j^{n_{\overline{D}}}}{L}_j \). Thus, the total number longitudinally repeated marker values in study is \( N={N}_D+{N}_{\overline{D}} \).
(AD1) Linear mixedeffect regression model
where \( {\boldsymbol{U}}_{\boldsymbol{D}}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \boldsymbol{s}\hfill & \hfill \begin{array}{cc}\hfill \boldsymbol{t}\hfill & \hfill \boldsymbol{s}\boldsymbol{t}\hfill \end{array}\hfill \end{array}\right],{\boldsymbol{\beta}}^{\boldsymbol{D}}={\left[\begin{array}{ccc}\hfill {\widehat{\boldsymbol{\beta}}}_0\hfill & \hfill {\widehat{\boldsymbol{\beta}}}_1\hfill & \hfill \begin{array}{cc}\hfill {\widehat{\boldsymbol{\beta}}}_2\hfill & \hfill {\widehat{\boldsymbol{\beta}}}_3\hfill \end{array}\hfill \end{array}\right]}^{\boldsymbol{T}},{\boldsymbol{U}}_{\overline{\boldsymbol{D}}} = \left[\begin{array}{cc}\hfill 1\hfill & \hfill \boldsymbol{s}\hfill \end{array}\right] \) and \( {\boldsymbol{\beta}}^{\overline{\boldsymbol{D}}}={\left[\begin{array}{cc}\hfill {\widehat{\boldsymbol{\beta}}}_0\hfill & \hfill {\widehat{\boldsymbol{\beta}}}_1\hfill \end{array}\right]}^{\boldsymbol{T}}. \)
(AD2) Model of ROC as a function of time prior to disease
 (i)
Construct a dataset of {(Y _{ ik }, Y _{ jl }), D = I(Y _{ ik } ≥ Y _{ jl })}
 (ii)
Calculate the quantile p in the control population (control observations in each pair as defined in step 1 above). It can be estimated by the empirical cumulative distribution function in the control sample.
 (iii)
The indicator I(.) in step 1 is estimated conditional on p in step 2. Thus, the ROC(p) is estimated by fitting a generalized linear model to data I(.), where the family is binomial, the link is probit and the covariates are Φ^{1}(p) and T _{ ik }.
There are a few advantages of this method compared to the first method in which the range of setting of this method is much broader [34]; the range of models that allowed for the ROC curve is broader; the model can include the interactions between p and U; the distributions of the test result in cases and controls do not need to be derived from the same family. Indeed, no assumptions are made regarding the distribution of marker for controls but only on the relationship between cases and controls which made through the ROC curve model.
(AD3) Semiparametric regression quantile estimation
where μ _{ D } σ _{ D }, \( {\mu}_{\overline{D}} \) and \( {\sigma}_{\overline{D}} \) are the location and scale functions. Instead of using a quasilikelihood method to estimate μ _{ D } ,σ _{ D }, \( {\mu}_{\overline{D}} \) and \( {\sigma}_{\overline{D}} \) [35], Zheng and Heagerty [9] used regression splines. In order to estimate the conditional baseline distribution function F and G, Zheng and Heagerty [9] proposed using an empirical distribution function of the standardized residuals if the baseline functions are independent of covariates, and to consider the symmetrized nearest neighbour (SNN) estimator [36] if the baseline functions are smooth functions of covariates. Thus, this semiparametric estimation method gives greater flexibility than the parametric method [32] by allowing separate model choices for each of the key distributional aspects.
(AD4) Cumulative/Dynamic definition extending for a longitudinal marker
Zheng and Heagerty [14] proposed a generalization of CD1 by Heagerty et al. [1] for longitudinal marker measurements. The key idea was the same as for the IS2 method in which the most recent marker is used to discriminate between cases prior to time t from controls after time t. Contrasted with CD1, it is no longer just the baseline marker or prognostic information that will be used but the updated information. The proposed sensitivity and specificity take the same form of CD1. In order to estimate the distribution function \( {\widehat{F}}_Y(c) \)(see CD1), Zheng and Heagerty [14] used the semiparametric regression quantile method for longitudinal data [35]. For the bivariate survival function, S(c, t), and the marginal survival function, S(t), Zheng and Heagerty [14] used a partly conditional hazard model as proposed by Zheng and Heagerty [37].
Motivated by the above methodology, we have extended CD2 to incorporate the most recent marker value instead of baseline marker value. CD2 is chosen rather than CD1 because CD1 produces nonmonotone sensitivity or specificity. The sensitivity and specificity are defined the same as CD2. This extended CD2 (denoted as ECD2) is assumed to have all the advantages of CD2 with an extra advantage of using the most recent marker value which is more reliable in depicting current status of an individual.
Software
The current software for computing the timedependent ROC curves are available as R packages. These are briefly described below.
survivalROC
The survivalROC [38] package estimates CD1and CD2. The R documentation includes worked examples using the builtin dataset called mayo (Primary Biliary Cirrhosis (PBC) dataset from Mayo Clinic). The estimators can be chosen by the type of method “KM” or “NNE” in the function syntax.
survAUC
The package [39] provides a variety of functions to estimate timedependent true/false positive rates and AUC for censored data. The AUC.cd can be used to calculate CD4 and it is restricted to Cox regression. The estimates obtained from this function are valid as long as the Cox model is specified correctly. The values returned by this function are AUC, integrated AUC and times at which the AUC are evaluated.
timeROC
The package [40] provides the functions to compute confidence intervals of AUC and tests for comparing AUC of two markers measured on the same individuals. Both CD5 and CD6 estimators can be computed by this package. It is also capable of allowing for competing risks event times.
survival, timereg and prodlim
The Basehaz function in the “survival” package [41] in R is used to obtain the VL Cox estimates which uses the baseline hazard under a Cox model. The Aalen function in the “timereg” package [42] can be used to estimate the conditional absolute risk under VL Aalen; it returns estimated coefficients β _{0} and β _{1}. The VL KM estimator can be computed using the “prodlim” package [43]. For the selection of the smoothing parameter l _{ n }, a direct plugin method can be used by setting l _{ n } to 0.25 n ^{− 1/5}.
risksetROC
This risksetROC package [44] estimates the timedependent ROC curves under I/D definition and produces accuracy measures for censored data under proportional or nonproportional hazard assumption of ID1 estimator.
Results
Examples of applications
Among the three definitions for sensitivity and specificity, C/D has been the most commonly applied in clinical papers (69/308, 22%). The I/D definitions have been applied in 14 papers (4.6%) while none was found for the I/S definitions. The detail on the review strategy is presented as a CONSORT diagram (Additional file 1: Fig. S1A) with a brief description of the process and the discussion about the remaining papers. Since the publication by Heagerty and Zheng [12] who introduced the three definitions, the number of clinical papers that used an I/D methodology has been increased (Additional file 1: Fig. S1B). Lung, breast and liver cancer are the most common areas for the application of C/D and I/D (Additional file 1: Fig. S1C). Some of the applications of C/D and I/D from cancer are described below.
Lu et al. [45] aimed to determine a robust prognostic marker for tumour recurrence as 30% of Stage I nonsmall cell lung cancer (NSCLC) patients will experience the tumour recurrence after therapy. They used timedependent ROC curve analysis to assess the predictive ability of gene expression signatures. The recurrencerelated genes were identified by performing a Cox proportional hazards analysis. A 51gene expression signature was validated as highly predictive for recurrence in Stage I NSCLC with AUC values greater than 0.85 from baseline up to 100 months of followup. The highest AUC values have been seen after 60 months to 100 months of followup with AUC(t) = 0.90, implying the 51gene expression signature is a better marker in discriminating between Stage 1 NSCLC patients who will experience tumour recurrence up to 60 months and patients who will not experience tumour recurrence beyond 60 months of followup. Lu et al. [45] concluded that this gene expression signature has important prognostic and therapeutic implications for the future management of these patients.
Tse et al. [46] has developed a prognostic risk prediction model for silicosis among workers exposed to silica in China using a Cox regression analysis to screen the potential predictors. The score from this model was then developed as a unique score system which includes 6 covariates: age at entry, mean concentration of respirable silica, net years of dust explore, smoking illiteracy and number of jobs. This score system was regarded as accurate in discriminating the workers with silicosis and healthy workers up to 600 months of followup since the AUC values are more than 0.80. These AUC values seems to decrease from baseline AUC(t = 0) = 0.96 to the end of followup AUC(t = 600) = 0.83 which indicates the discrimination potential of baseline score is diminished across study followup. This study provides scientific guidance to the clinicians to identify highrisk workers.
Yue et al. [47], [48] have used pretreatment 18 FFDGPET/CT imaging and combinatorial biomarkers respectively to stratify the risk of TNBC (Triplenegative breast cancer) patients. TNBC is considered as a high risk disease and normally associated with poor survival. A stratification of prognosis of this disease can help in identifying the patients with good prognosis for less aggressive therapy. The eventtime outcome of the studies was defined as the time to recurrence from TNBC disease. Timedependent ROC curve was used to assess the prognostic value of the biomarkers, EFGR and CK5/6 at different cutoff points and the optimal cutoff was obtained based on the AUC values. The cutoff values were estimated by maximizing both sensitivity and specificity of the eventtime outcome. The optimal values of 15% with AUC = 0.675 and 50% with AUC = 0.611 for EFGR and CK5/6 were found respectively. AUC values obtained were used as a basis of a decision rule. By using the optimal cutoff value, the patients were stratified into two different risk level groups which helps in selecting the appropriate treatment strategies for patients.
Desmedt et al. [49] have studied the performance of the gene expression index (GGI) in predicting relapses in postmenopausal women who were treated with tamoxifen (T) or letrozole (L) within the BIG 198 trial. The predictive ability of GGI was estimated using timedependent AUC and has been plotted as a function of time to characterize temporal changes in accuracy of the GGI marker. They have calculated AUC(t = 24) = 0.73 which implies that 73% of the patients who relapse at 24^{th} month have greater GGI score than patients who relapse after 24^{th} month. Further, AUC at = 27) was found to be the highest which indicates that the maximal discrimination occurs near the median followup time.
George et al. (2012) aimed to determine the predictive ability of lesions texture along with traditional features in order to detect the early tumour response. Texture features are important in detecting the progression of tumour among cancer patients, e.g. s (18)Ffluorodeoxyglucose (FDG) followed with positron emission tomography (PET) estimates. The eventtime outcome was defined as the time to tumour progression, which is the distance between subspaces from baseline scan and followup scan. Timedependent ROC curve is used to obtain the predictive ability of the weighted subspacesubspace distance from the baseline and the followup scan as a marker for predicting early tumour response. In a study of 15 patients who had metastatic colorectal cancer, the followup scan was taken at the first week after the first dose of the treatment. As a result, a concordance summary of 0.68 is found from the predictive model using weighted subspacesubspace distance metrics. This result helps as an added value in using textural information for therapy response evaluation.
Illustrative application
We have used sequential data from a randomized placebocontrolled trial of the drug Dpenicillamine (DPCA) for the treatment of primary biliary cirrhosis (PBC) conducted at the Mayo Clinic between 1974 and 1984 [50] in order to illustrate the performance of the current methods in estimating the timedependent ROC curves. The eventtime outcome of this study is the time to death due to PBC liver disease. The original clinical protocol for the study specified visits at 6 months, 1 year, and annually thereafter. We use a model score estimated from the Cox model containing five covariates: log(bilirubin), albumin, log(prothrombin time), edema and age as the marker [12].
Estimated timedependent AUC for Year 1, Year 5 and Year 10
Definitions  Marker time  Method  AUC (SD)  

Year 1  Year 5  Year 10  
C/D  t = 0  Naïve  0.846 (0.023)  0.885 (0.022)  0.883 (0.030) 
CD1  0.922 (0.041)  0.921 (0.021)  0.878 (0.027)  
CD2  0.895 (0.056)  0.897 (0.024)  0.869 (0.028)  
CD3  0.922 (0.042)  0.917 (0.020)  0.898 (0.031)  
CD5  0.922 (0.042)  0.915 (0.021)  0.866 (0.028)  
CD6  0.922 (0.038)  0.915 (0.020)  0.870 (0.030)  
C/D  Last value prior to:  ECD2  
Year 1  0.926 (0.039)  0.918 (0.019)  0.871 (0.027)  
Year 5    0.911 (0.019)  0.910 (0.021)  
Year 10      0.899 (0.022)  
I/D  t = 0  ID1  0.845 (0.010)  0.791 (0.028)  0.692 (0.024) 
ID3  0.893 (0.048)  0.757 (0.041)  0.716 (0.143)  
I/S  t = 0  IS2  0.939 (0.025)  0.836 (0.028)  0.698 (0.034) 
I/S  Last value prior to:  IS2  
Year 1  0.968 (0.003)  0.872 (0.024)  0.698 (0.043)  
Year 5    0.957 (0.003)  0.698 (0.031)  
Year 10      0.768 (0.038) 
Parameter estimates for linear mixed effect model
Case  Control  

Fixed Effect  \( {\widehat{\beta}}_0(SE) \)  1.139(8.865 × 10^{2})  0.569 (0.043) 
\( {\widehat{\beta}}_1(SE) \)  4.813 × 10^{4}(4.419 × 10^{5})  2.906 × 10^{4} (2.502 × 10^{5})  
\( {\widehat{\beta}}_2(SE) \)  2.283 × 10^{4}(5.696 × 10^{5})  
\( {\widehat{\beta}}_3(SE) \)  1.083 × 10^{7}(1.605 × 10^{8})  
Random Effect  \( {\widehat{\sigma}}_{int} \)  0.593  0.550 
\( {\widehat{\sigma}}_{VT} \)  3.448 × 10^{4}  2.615 × 10^{4}  
\( {\widehat{\rho}}_{int, VT} \)  0.378  0.209  
\( {\widehat{\sigma}}_{Res} \)  0.293  0.220 
Discussion and conclusions
Although C/D is the most commonly being applied, if a researcher has a specific time point of interest in order to distinguish between individuals with an event and individuals without event at that time point, I/D or I/S is more appropriate. Since I/S requires a fixed followup to observe the clinical outcome of interest, it can be applied in long followup studies with longitudinally measured markers. C/D and I/D are usually used for a single biomarker value while I/S can include a longitudinal biomarker. As the disease status of an individual may change during followup, the biomarker values may also change, and hence, the most recent marker value may be best related to the current disease status of an individual. Thus, usage of the most recent marker value prior to a target prediction time t is acceptable as we discussed using the extended methods.
The optimal cutoff is determined by choosing the highest AUC value in which describes the marker has the largest separation between cases and controls. In general, the cutoff (also called as threshold) is chosen based on the availability of the healthcare resources and the level of disease severity.
None of the methods discussed earlier used a complete history of longitudinal marker conditional on an eventtime. The approach of considering a more complete record of each individual when estimating the ROC summaries over time can be more appropriate. A joint modelling framework in an attempt to estimate the timedependent ROC curve is recommended since it considers the association between longitudinal marker and the corresponding eventtime processes. Further, it is also suggested to assume the event times to be parametrically distributed which is then be easier to estimate the survival function if a researcher is attempting to extend for measurement error.
Abbreviations
 18 FFDGPET/CT:

Fluorine18fluorodeoxyglucose positron emission tomography/computed tomography
 AIDS:

Acquired immune deficiency syndrome
 AUC:

Area under ROC curve
 C/D:

Cumulative/Dynamic
 CD4:

T–lymphocyte cell bearing CD4 receptor
 DPCA:

Drug Dpenicillamine
 FDG:

Ffluorodeoxyglucose
 GGI:

Gene expression index
 HIV:

Human immunodeficiency virus
 I/D:

Incident/Dynamic
 I/S:

Incident/Static
 NSCLC:

NonSmall Cell Lung cancer
 OCT:

Optical coherence tomography
 PBC:

Primary biliary cirrhosis
 PET:

Positron emission tomography
 RNFL:

Retinal nerve fiber layer
 ROC:

Receiver operating characteristic
 TNBC:

Triplenegative breast cancer
Declarations
Acknowledgements
ANK is funded by Majlis Amanah Rakyat (MARA), Malaysian Government PhD studentship.
Funding
Not applicable
Availability of data and materials
The PBC dataset analysed during the current study is publically available from http://www.mayo.edu/research/documents/pbcseqhtml/doc10027141
Authors’ contributions
ANK performed the review and analysis of data, and prepared the manuscript for publication. RKD and TC jointly conceived the initial idea, and provided comments and assistance throughout the study on both the extracted data and drafts of the manuscript. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable
Ethics approval and consent to participate
Not applicable
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis 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
 Heagerty PJ, Lumley T, Pepe MS. Timedependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.View ArticlePubMedGoogle Scholar
 Hung H, Chiang CT. Estimation methods for timedependent AUC models with survival data. Can J Stat Revue Can Stat. 2010;38(1):8–26.Google Scholar
 Song X, Zhou XH, Ma S. Nonparametric receiver operating characteristicbased evaluation for survival outcomes. Stat Med. 2012;31(23):2660–75.View ArticlePubMedPubMed CentralGoogle Scholar
 Chambless LE, Diao G. Estimation of timedependent area under the ROC curve for longterm risk prediction. Stat Med. 2006;25(20):3474–86.View ArticlePubMedGoogle Scholar
 Lambert J, Chevret S. Summary measure of discrimination in survival models based on cumulative/dynamic timedependent ROC curves. Stat Methods In Med Res. 2014;25(5):2088–102.
 Cai T, Pepe MS, Lumley T, Zheng Y, Jenny NJ. The sensitivity and specificity of markers for event times. Biostatistics. 2006;7(2):182–97.View ArticlePubMedGoogle Scholar
 Kalbfleisch JD, Prentice RL. The statistical analysis of failure time data, vol. 360. New Jersey: Wiley; 2011. https://books.google.co.uk/books?hl=en&lr=&id=BR4Kqa1MIMC&oi=fnd&pg=PR7&dq=The+statistical+analysis+of+failure+time+data&ots=xCsg6MQU7_&sig=gf4iHW8SkymUSR4RSzNCZsbTGdY#v=onepage&q=The%20statistical%20analysis%20of%20failure.
 Pepe MS. The statistical evaluation of medical tests for classification and prediction. USA: Oxford University Press; 2003.Google Scholar
 Zheng Y, Heagerty PJ. Semiparametric estimation of timedependent ROC curves for longitudinal marker data. Biostatistics. 2004;5(4):615–32.View ArticlePubMedGoogle Scholar
 Bamber D. The area above the ordinal dominance graph and the area below the receiver operating characteristic graph. J Math Psychol. 1975;12(4):387–415.View ArticleGoogle Scholar
 Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982;143(1):29–36.View ArticlePubMedGoogle Scholar
 Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005;61(1):92–105.View ArticlePubMedGoogle Scholar
 Blanche P, Dartigues JF, JacqminGadda H. Review and comparison of ROC curve estimators for a timedependent outcome with markerdependent censoring. Biom J. 2013;55(5):687–704.View ArticlePubMedGoogle Scholar
 Zheng Y, Heagerty PJ. Prospective accuracy for longitudinal markers. Biometrics. 2007;63(2):332–41.View ArticlePubMedGoogle Scholar
 Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53(282):457–81.View ArticleGoogle Scholar
 Akritas MG. Nearest neighbor estimation of a bivariate distribution under random censoring. Ann Stat. 1994;1299–1327.
 Cai T, Gerds TA, Zheng Y, Chen J. Robust Prediction of t‐Year Survival with Data from Multiple Studies. Biometrics. 2011;67(2):436–44.View ArticlePubMedGoogle Scholar
 Hung H, Chiang CT. Optimal Composite Markers for TimeDependent Receiver Operating Characteristic Curves with Censored Survival Data. Scand J Stat. 2010;37(4):664–79.View ArticleGoogle Scholar
 Song X, Zhou XH. A semiparametric approach for the covariate specific ROC curve with survival outcome. Statistica Sinica. 2008;18(3):94765.
 Viallon V, Latouche A. Discrimination measures for survival outcomes: connection between the AUC and the predictiveness curve. Biom J. 2011;53(2):217–36.View ArticlePubMedGoogle Scholar
 Uno H, Cai TX, Tian L, Wei LJ. Evaluating prediction rules for tyear survivors with censored regression models. J Am Stat Assoc. 2007;102(478):527–37.
 Royston P, Parmar MK. The use of restricted mean survival time to estimate the treatment effect in randomized clinical trials when the proportional hazards assumption is in doubt. Stat Med. 2011;30(19):2409–21.View ArticlePubMedGoogle Scholar
 Cox DR. lRegression Models and Life Tables. mJ R Stat Soc Ser B. 1972;34(2):187–220.Google Scholar
 Aalen OO. A linear regression model for the analysis of life times. Stat Med. 1989;8(8):907–25.View ArticlePubMedGoogle Scholar
 Cai Z, Sun Y. Local Linear Estimation for Time‐Dependent Coefficients in Cox's Regression Models. Scand J Stat. 2003;30(1):93–111.View ArticleGoogle Scholar
 Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994;81(3):515–26.View ArticleGoogle Scholar
 Xu R, O'Quigley J. Proportional hazards estimate of the conditional survival function. J R Stat Soc Ser B (Stat Methodol). 2000;62(4):667–80.View ArticleGoogle Scholar
 SahaChaudhuri P, Heagerty PJ. Nonparametric estimation of a timedependent predictive accuracy curve. Biostatistics. 2013;14(1):42–59.View ArticlePubMedGoogle Scholar
 Shen W, Ning J, Yuan Y. A direct method to evaluate the time‐dependent predictive accuracy for biomarkers. Biometrics. 2015;71(2):439–49.View ArticlePubMedPubMed CentralGoogle Scholar
 Royston P, Altman DG. Regression Using Fractional Polynomials of Continuous Covariates  Parsimonious Parametric Modeling. Appl StatJ Roy St C. 1994;43(3):42967.
 Leisenring W, Pepe MS, Longton G. A marginal regression modelling framework for evaluating medical diagnostic tests. Stat Med. 1997;16(11):1263–81.View ArticlePubMedGoogle Scholar
 Etzioni R, Pepe M, Longton G, Hu C, Goodman G. Incorporating the time dimension in receiver operating characteristic curves: a case study of prostate cancer. Med Decis Mak. 1999;19(3):242–51.View ArticleGoogle Scholar
 Tosteson ANA, Begg CB. A general regression methodology for ROC curve estimation. Med Decis Mak. 1988;8(3):204–15.View ArticleGoogle Scholar
 Pepe MS. Three approaches to regression analysis of receiver operating characteristic curves for continuous test results. Biometrics. 1998; 54(1):12435.
 Heagerty PJ, Pepe MS. Semiparametric estimation of regression quantiles with application to standardizing weight for height and age in US children. J R Stat Soc: Ser C: Appl Stat. 1999;48(4):533–51.View ArticleGoogle Scholar
 Yang SS. Linear combination of concomitants of order statistics with application to testing and estimation. Ann Inst Stat Math. 1981;33(1):463–70.View ArticleGoogle Scholar
 Zheng Y, Heagerty PJ. Partly conditional survival models for longitudinal data. Biometrics. 2005;61(2):379–91.View ArticlePubMedGoogle Scholar
 Heagerty PJ, SahaChaudhuri P, SahaChaudhuri MP. Package ‘survivalROC’. 2013.Google Scholar
 Potapov S, Adler W, Schmid M: survAUC: Estimators of Prediction Accuracy for TimetoEvent Data. R package version 1.05. In.; 2012.
 Blanche P. TimeROC: Timedependent ROC curve and AUC for censored survival data. R package version 02, URL https://cran.rproject.org/web/packages/timeROC/timeROC.pdf.
 Therneau TM, Lumley T. Package ‘survival’. In.: Verze; 2015
 Scheike T. Timereg Package. In.: R Package Version; 2009.
 Gerds TA, Rcpp I, Rcpp L, Gerds MTA. Package ‘prodlim’. 2015.Google Scholar
 Heagerty PJ, SahaChaudhuri P, SahaChaudhuri MP. Package ‘risksetROC’. 2012.Google Scholar
 Lu Y, Wang L, Liu P, Yang P, You M. Geneexpression signature predicts postoperative recurrence in stage I nonsmall cell lung cancer patients. PLoS One. 2012;7(1):e30880.View ArticlePubMedPubMed CentralGoogle Scholar
 Tse LA, Dai JC, Chen MH, Liu YW, Zhang H, Wong TW, Leung CC, Kromhout H, Meijer E, Liu S et al. Prediction models and risk assessment for silicosis using a retrospective cohort study among workers exposed to silica in China. Scientific Reports. 2015;5.
 Yue Y, Cui X, Bose S, Audeh W, Zhang X, Fraass B. Stratifying triplenegative breast cancer prognosis using 18 FFDGPET/CT imaging. Breast Cancer Res Treat. 2015;153(3):607–16.View ArticlePubMedPubMed CentralGoogle Scholar
 Yue Y, Astvatsaturyan K, Cui X, Zhang X, Fraass B, Bose S. Stratification of Prognosis of TripleNegative Breast Cancer Patients Using Combinatorial Biomarkers. PLoS One. 2016;11(3):e0149661.View ArticlePubMedPubMed CentralGoogle Scholar
 Desmedt C, GiobbieHurder A, Neven P, Paridaens R, Christiaens MR, Smeets A, Lallemand F, HaibeKains B, Viale G, Gelber RD. The Gene expression Grade Index: a potential predictor of relapse for endocrinetreated breast cancer patients in the BIG 1–98 trial. BMC Med Genet. 2009;2(1):1.Google Scholar
 Fleming TR, Harrington DP. Counting processes and survival analysis, vol. 169. New Jersey: Wiley; 2011. doi:10.1002/9781118150672.