Multistate model for studying an intermediate event using timedependent covariates: application to breast cancer
 Carolina MeierHirmer^{1}Email author and
 Martin Schumacher^{2}
DOI: 10.1186/147122881380
© MeierHirmer and Schumacher; licensee BioMed Central Ltd. 2013
Received: 6 January 2013
Accepted: 8 June 2013
Published: 20 June 2013
Abstract
Background
The aim of this article is to propose several methods that allow to investigate how and whether the shape of the hazard ratio after an intermediate event depends on the waiting time to occurrence of this event and/or the sojourn time in this state.
Methods
A simple multistate model, the illnessdeath model, is used as a framework to investigate the occurrence of this intermediate event. Several approaches are shown and their advantages and disadvantages are discussed. All these approaches are based on Cox regression. As different timescales are used, these models go beyond Markov models. Different estimation methods for the transition hazards are presented. Additionally, timevarying covariates are included into the model using an approach based on fractional polynomials. The different methods of this article are then applied to a dataset consisting of four studies conducted by the German Breast Cancer Study Group (GBSG). The occurrence of the first isolated locoregional recurrence (ILRR) is studied. The results contribute to the debate on the role of the ILRR with respect to the course of the breast cancer disease and the resulting prognosis.
Results
We have investigated different modelling strategies for the transition hazard after ILRR or in general after an intermediate event. Including timedependent structures altered the resulting hazard functions considerably and it was shown that this timedependent structure has to be taken into account in the case of our breast cancer dataset. The results indicate that an early recurrence increases the risk of death. A late ILRR increases the hazard function much less and after the successful removal of the second tumour the risk of death is almost the same as before the recurrence. With respect to distant disease, the appearance of the ILRR only slightly increases the risk of death if the recurrence was treated successfully.
Conclusions
It is important to realize that there are several modelling strategies for the intermediate event and that each of these strategies has restrictions and may lead to different results. Especially in the medical literature considering breast cancer development, the timedependency is often neglected in the statistical analyses. We show that the timevarying variables cannot be neglected in the case of ILRR and that fractional polynomials are a useful tool for finding the functional form of these timevarying variables.
Background
In this article the disease course of primary, nonmetastatic breast cancer patients is analysed taking first isolated locoregional recurrence (ILRR) of the tumour into account as an intermediate event. In this context, usually separate analyses are carried out for each endpoint and also for the intermediate events. As pointed out by [1] “these separate analyses are not completely satisfying, since they fail to reveal the relations between different types of events”. In this article the illnessdeath model is used. This is a type of multistate model nowadays widely used in medical research [2] for describing chronic diseases and intermediate events, particularly in oncology [1, 3, 4].
The presented approach focuses on the analysis of the hazard rate after the first recurrence (intermediate event) in order to understand the influence of the time of occurrence of this event on the further disease development. As endpoints we consider either “distant metastases or death” or “death”. Different methods for modelling and for estimating the hazard function after recurrence are provided. These methods are very general and can be applied to other problems if the course of the disease to be modelled includes an intermediate event. The methods are applied to a dataset containing data from four studies conducted by the German Breast Cancer Study Group (GBSG) see e.g.[5].
In the medical literature, the role of the ILRR with respect to the course of the breast cancer disease and the resulting prognosis is not clearly established. Different hypotheses are under debate: On the one hand, the ILRR could be seen as an independent event, which does not alter prognosis after successful removal, or, on the other hand, the ILRR could be a first sign of disease progression [6–8] and even the cause of distant metastases [9]. The impact of time to recurrence on the death hazard after an ILRR was demonstrated in [10].
The models proposed in this article are more general than Markov models for the following reasons:

The transition hazards between different disease states are allowed to depend on previous state occupation times. This is necessary for the analysis of the impact of the time to recurrence on the future disease course.

After the occurrence of an intermediate event, there are two different time scales. The first one is the time since entry into the study or since the diagnosis of the primary tumour, which is usually the time axis used before the intermediate event. The second possible time scale is the time since the occurrence of the intermediate event, which can also play an important role. When constructing the statistical model for a particular disease, it has to be be taken care of which time scale and which previous state occupation times are important for the further development of the disease.

As pointed out by [11], it is important to use the information contained in covariates changing over time. This information is a potential predictor that should not be neglected or substituted by timeconstant covariates.
Concerning the isolated local recurrence of breast cancer, the dependence of the hazard function after the ILRR on the sojourn time is analysed. This means that the effect of the time since occurrence of the ILRR on the further disease development is examined. As stated above, the models proposed in this article are constructed in order to allow time dependence. In medical literature, this timedepending structure of the ILRR is often neglected. In some models the ILRR is even analysed as a “standard” prognostic factor, i.e. like a covariate measured at the time of diagnosis of the primary tumour.
In the next section of this article, the illnessdeath model is described in detail. Based on a Cox regression model, it is shown how different approaches determine the shape of the transition hazard function. A method based on fractional polynomials is proposed for the analysis of the functional form of the timevarying covariates. In Section Application, the presented methods are used to analyse data from the German Breast Cancer Study Group (GBSG). Additionally to the timevarying covariates, also the effects of “standard” prognostic factors are considered. The results are then analysed and discussed.
Methods
When using common regression models for censored data, there are three possible approaches for the estimation of the transition hazards: separate models for every transition (approach S), timedependent covariates (approach J) [12], and the stratified model [13]. There are also mixture models, e.g. one proposed in [2]. In the following paragraphs, we explain how these methods estimate the transition hazards and which assumptions are (indirectly) made.
which we denote only by d and t  d in the following. We are now interested in the functional form with which these variables enter the hazard function. We therefore transform (t  d) and d by a function f: f(d,t  d). At this stage, the shape of the function f is not yet determined. In paragraph Shape of the timevarying covariate some possible choices are proposed for doing this.
Problems can arise when sojourn times and in general random timedependent covariates are included into multistate models. For example, cumulative incidence functions cannot longer be estimated in a straightforward way. Such situations were studied by [15] and [11]. The latter paper points out that “studying [...] sojourn times as timedependent covariates may be useful for testing model assumptions and for investigating their effect on the survival” and treats this kind of problem in section 4.1 [11]. The model mentioned by these authors is equivalent to our straight line model of Section Shape of the timevarying covariate.
The following two paragraphs discuss the joint transition hazards approach J and the separate transition hazards approach S.
Joint approach
Possible relations between transitions 1 → 3 and 2 → 3 if approach J is used
J.I  λ _{23}(t,d)=λ _{13}(t)  The hazard rate remains the same 
J.II  λ _{23}(t,d)=c·λ _{13}(t)  The hazard rate changes by 
the factor c (constant)  
J.III.1  λ _{23}(t,d)=c(d)·λ _{13}(t)  c depends on d, the 
time of ILRR  
J.III.2  λ _{23}(t,d)=c(t  d)·λ _{13}(t)  c depends on t  d, the 
time since ILRR  
J.IV  λ _{23}(t,d)=c(d,t  d)·λ _{13}(t)  c depends on d as well as 
on t  d the time since ILRR 
Models J.III (1 and 2) and model J.IV violate the Markov assumption. The estimated parameters of the timevarying covariates allow for testing directly the difference between the hazards before and after the intermediate event. The models are nested and it is therefore possible to use directly a criterion (AIC for example) to assess the goodness of fit. This constitutes an advantage of this approach.
Separate approach
The timeinvariant covariates are allowed to have different effects on all hazards. The timedependency of the hazard 2 →3 on the time of recurrence d and the time since recurrence t  d is modelled via the function f. Approach J is a special case of model S.
Approach S is a special case of the stratified model [13]. The stratified model assumes different baseline hazards for all transitions. The covariates can be chosen to have identical effects for all transitions or can be stratified by transition whereas the baseline hazards remain different for all three transitions. As the main interest of this article is the modelling of the timedependent covariates and not the effect of the “standard” covariates, we use approach S as described above.
Shape of the timevarying covariate
Fractional polynomials (FP) are used for the modelling of the function f. Fractional polynomials were introduced in [16] and are explained in detail in [17]. In order to carry out model selection using FPs, a sequential selection procedure was proposed in [18]. Our approach is similar to the approach of [19], where time  covariate interaction was investigated using FPs; for a comprehensive overview see [20].
By contrast, in the situation investigated here, fractional polynomials are directly applied on the timevarying covariates d and t  d. The resulting FP gives the function f which describes the shape of the impact of the timedependent covariates on the hazard function λ _{23}. Using approach J, it also determines the relation between function λ _{13} and λ _{23} as c(d,t  d)= exp (β ^{ T } f(d,t  d)).
The procedure is the following: Suppose that X is a timedependent covariate. We are interested in finding the function f(.) that describes the impact of the timedependent covariate on the hazard function. The starting point is a straight line model f(X) = β _{1} X. In some cases, this is already an adequate description of the relationship, but other models are analysed in order to improve the fit. The extension of the straight line model are power functions β _{1} X ^{ p }. The values of p are chosen from the set S = {2,1,0.5,0,0.5,1,2,3} where X ^{0} denotes log(X). The resulting functions are called oneterm fractional polynomials or FP1 functions. The variable that counts the number of FPterms m is set to 1. In this article also FP2 functions are used: $f\left(X\right)\phantom{\rule{.1em}{0ex}}=\phantom{\rule{.1em}{0ex}}{\beta}_{1}{X}^{{p}_{1}}\phantom{\rule{.1em}{0ex}}+\phantom{\rule{.1em}{0ex}}{\beta}_{2}{X}^{{p}_{2}}$, with ${p}_{1},{p}_{2}\phantom{\rule{.1em}{0ex}}\in \phantom{\rule{.1em}{0ex}}S$. In this case the number of FPterms increases m = 2. If p _{1} = p _{2}, the function ${\beta}_{1}{X}^{{p}_{1}}+{\beta}_{2}{X}^{{p}_{1}}log\left(X\right)$ is used.
In order to find the best model, the 8 different FP1 functions and the 36 different FP2 functions are fitted and compared using the AIC criterion [21].
If there are several timedependent covariates, linear combinations of the FP of every covariate are used. A model selection procedure treating several covariates is described in [18].
Application
In this section, the methods described above are applied to a real dataset. In the first paragraph, a description of the dataset is given. Then it is shown how the analysis of the data was carried out.
Data
Between 1983 and 1989, four studies including patients with primary, histologically proven, nonmetastatic breast cancer were carried out by the German Breast Cancer Study Group (GBSG). 2746 patients from 118 clinical institutions entered the studies. The treatment design of the studies is shown in [5], Table 1. The following covariates were registered at primary diagnosis: patient’s age, menopausal status, number of positive axillary lymph nodes, tumour location, tumour size, histologic tumour grade, oestrogen, and progesterone receptor status. Patients were examined at regularly scheduled followup visits. The study was carried out with a very smallmeshed followup scheme in particular in the first years after primary surgery in order to detect any kind of recurrence at the earliest time possible [10]. In the first two years, for example, the examinations took place every 3 month. More details are available in [10]. Followup was continued until death. Censoring took place at the end of the study or if followup was not possible.
We reduced the data to 2390 patients in order to avoid a missing value problem. This number does not match the number of patients used in [10] as we also excluded patients with survival time equal to zero (N = 15). These patients are all censored individuals, lost to followup after primary surgery. If ILRR was detected at death we defined the recurrence to have taken place one month earlier. If ILRR and distant metastases were detected at the same time, the time of recurrence was distributed uniformly on an interval of three month before manifestation of the distant disease.
In order to simplify the application of the proposed methods, the timeinvariant covariates were categorized, using predefined cutpoints. Only patient’s age, number of positive lymph nodes, tumour size, tumour grade and progesterone receptor status were considered. The other covariates were eliminated within the first model selection using backward elimination at a level α=0.01.
Implementation
Two illnessdeath models are analysed. In the fist one, the overall survival is estimated, the absorbing state is therefore “death”. Patients with distant disease but without ILRR are censored. In the second model, the endpoint is “distant metastases or death”.
We start with approach J which forces the hazard ratio λ _{23}/λ _{13} to have the shape of function c. Thereafter, it is investigated if this assumption has to be relaxed by calculating model S with t as underlying time scale.
Then, the time of recurrence d is included into the analysis. Using stepfunctions to model this timedependence, it was already shown by [10] that d has a linear impact on the overall survival after recurrence. This is why d will be included in the model without transformation. It is not known if the effect of time d is also linear with respect to distant disease free survival but we will assume this first.
Concerning the time after ILRR, there is a hypothesis that the risk of death is very high just after the recurrence and diminishes with time. But the ILRR could also be an indicator of an accelerated disease development and therefore increase the hazard rate enduringly.
Based on empirical knowledge, a possible shape of the time after recurrence could be exp ((t  d)). The more the corresponding parameter estimate is positive, the more serious is the occurrence of an ILRR but the faster is the reduction of the risk to a normal level. In this case, the ILRR could be considered as a new and independent disease that does not alter the hazard rate on the long term. If the parameter estimate is negative, the risk increases up to a saturation level. This corresponds to the hypothesis that the ILRR is a sign of a deterioration of the patient’s health.
which is a special case of relation J.IV. By omitting successively β _{2}, β _{1} and β _{0}, the model is reduced to relation J.III (1 or 2), J.II and J.I, respectively. These models are then compared using the AIC criterion.
In a further step the intuitive choice of the function exp ((t  d)) is justified by comparing the resulting hazard function to the hazard function of a timedependent FP.
The analysis is carried out a second time for the distant disease free survival.
The resulting hazard ratios for approach J and S are plotted and compared graphically.
The analysis was done in R[22]. The Cox model was estimated using the procedure cph of the rms package written by Frank Harrell [23] and on the coxph procedure of the survival package written by Terry Therneau [24]. It is necessary to use the countingprocess notation for survival data in order to account for timedependent covariates and lefttruncation. The modification of the data and the implementation of the programs are described in [13] and [25]. Cumulative baseline hazards were smoothed and first derivatives were calculated by the smoothing function smooth.spline of R. Other methods for deriving the baseline hazard functions like kernel estimators were also applied, but are not shown here. The FP fit was done manually. At present, only software for fitting timeinvariant FPs [26] or time  covariate interaction FPs [19, 20] exists.
Results
The following two paragraphs contain the results of the analysis respectively for overall and distant disease free survival.
Results for overall survival
In order to know if the function exp ((t  d)), based only on expert knowledge and describing the dependence of the hazard function on the time since recurrence (t  d), accords with reality, FPs were applied on this timevarying covariate. All FP1 and FP2 models were tested and compared using the AIC. The FP1 model (p _{1}=1) is the best one and is denoted by J.FP. As can be seen in Figure 4, the J.FP model has a slightly lower AIC but the shape of log (c(t,t  d)) remains nearly the same.
Results for distant disease free survival
The same analysis was done for the outcome variable “distant metastases or death”. The parameter estimate of the time to recurrence d was not significantly different from zero in all approaches (α=0.05). This covariate is therefore omitted. Then the functional form of the covariate t  d is analysed using model J.III.2.
In the next step, approach S is applied to the data.
The functional form of t  d was determined by the FPapproach and the FP m = 2,p _{1} = log,p _{2} = log2 was the best model for both approaches, J and S. The shapes of the hazard ratios are similar for both approaches. The resulting curve of approach J corresponds approximately to the curve with d = 2.5 years of approach S, indicating that approach J takes into account a mean value of time of recurrence d as this covariate was omitted from the model. Approach S shows this dependence on d via the baseline hazards.
Discussion
The often referenced Stanford heart transplant data provided a starting point in illustrating the problem of the comparison of survival curves before and after an intermediate event [27]. Since, disease courses which should be modelled including intermediate events were identified in many medical fields. The course of bone marrow transplantation [28] or the benefit of lung transplantation [29, 30] are examples for this development. In the latter two publications, approach J was applied including timevarying effects.
In this article, a multistate model is used for the analysis of breast cancer data including an intermediate event (the ILRR). The study presented here is exactly in the scope of the articles [31] and [32]. These articles deal with illnessdeath models, where the intermediate event is intervalcensored and where there are patients with unknown status at the time of death or followup. In our case, ILRR and metastases events are intervalcensored as the followup took place periodically. The four studies had a median followup time between 7 and 9 years, going up to 10 years for the overall survival [10]. The patients are examined at intervals, the ILRR event is the time of the first follow up by which ILRR has occurred. The sojourn time without recurrent event could in this way be overestimated and therefore the lifetime after the recurrent event underestimated. But the followup periods are very short in comparison to the study period (cf. Section Data). We therefore chose to not take into account the intervalcensoring. The bias related to an analysis without intervalcensoring should in that case be comparatively small. There was a high number of patients dying from other causes than breastcancer (67 patients) and autopsy was not compulsory. This could suggest that ILRR or distant disease were present at death but the status was unknown by the study organisers. It has to be taken into account, however, that one third of the patients were older than 60 years at the beginning of the study. So there was a high probability for these patients to die from other causes. An “unknown status” was therefore not incorporated into the analysis that again should not be associated with a large bias.
This article demonstrates that the shape of the hazard ratio in the “illnessdeath” model depends on the inclusion or exclusion of time depending covariates. Comparisons between model J.I up to model J.IV with “death” as outcome variable confirm that the ILRR should be modelled as a timedependent intermediate event. The results allow the determination of the impact of the ILRR on the course of the disease. Comparing the results of model J.I and J.II it is obvious that the ILRR is a timedependent phenomenon and affects the hazard rate. Entering the study, the patients are under risk of ILRR or death. After having experienced an ILRR they are only at risk of death. Therefore a model has to be chosen where the assignment of the observations to different risks is not fixed in advance but dynamic in time. This is ignored in many publications when the patients are grouped retrospectively by recurrence as in [9]. Models that ignore this timedependency cannot reflect reality and are prone to the socalled timedependent bias, see e.g.[4]. The real effect of the intermediate event cannot be assessed in this way and may be underestimated.
The comparison of model J.II and J.III.1 indicates that the time of recurrence affects overall survival. The later the recurrence appears, the lower is the risk of death. The comparison between model J.III.1 and J.IV implies that the hazard function depends also on time since ILRR. With respect to overall survival (cf. Figure 5), approach J restricts too much the shape of the hazard ratio. The relation between the hazards in approach J does not model the difference between the hazard rates satisfyingly and is therefore not the appropriate model. Approach S shows that an early ILRR leads to a sudden increase of the risk of death which decreases relatively fast after the recurrence, but that a late ILRR (after 2 years) only increases the hazard by a factor of about 2.5 which does not longer depend on time since recurrence t  d. From a medical point of view, it is interesting to know whether the appearance of an ILRR is due to the first disease and expresses the declining health of the patient, or if the ILRR is independent from the first disease. The results show that the ILRR is not completely independent from the first disease. Indeed, the risk of death remains on a higher level after an ILRR (ratio >1). But there is also an independent part: approach S indicates that an early recurrence leads to an additional risk, but this risk diminishes after successful removal of the second tumour.
With respect to distant disease free survival, approach J models the shape of the ratio correctly, but the dependency on time of recurrence d could not be analysed. In fact, in all J models this variable was not significant. The model S however shows this dependency on d via the baseline hazard. This proves that the proportionality assumption between the hazard functions λ _{13} and λ _{23} in model J does not match the data. Whenever the ILRR appears, the risk of distant disease has a peak thereafter. This peak is more pronounced if the ILRR occurs early (Figure 7). After two years, the risk has diminished. In this case, the ILRR can be seen as an independent event. If the recurrence is treated with success, the risk of a distant disease decreases and attains the previous level.
The results also show that fractional polynomials are useful to determinate the functional form of timevarying covariates. Even if the results of the J.IV model using a predefined exponential function are almost equivalent in the case of overall survival, the results show that the use of FPs is the best choice for distant disease free survival.
Conclusions
Summing up, we have investigated different modelling strategies for the transition hazard after ILRR or in general after an intermediate event. An example is given in which including timedependent structures alters the resulting hazard function considerably. It is important to realize that there are several modelling strategies and that each of these strategies has certain restrictions and may lead to different results.
It was also shown that fractional polynomials are a useful tool for finding the functional form of the timevarying variables.
Diagnostic methods have changed over the last 20 years. The aim of the paper was to show that the timevarying effect of the ILRR has not to be neglected. The methods developed for doing this are still valid nowadays. The aim of this article was not to reanalyse the four studies of the German BreastCancer Study Group.
Declarations
Acknowledgements
We like to thank Claudia Schmoor for her help to use the dataset.
Authors’ Affiliations
References
 Putter H, van der Hage, de Bock GH, Elgalta R, van de Velde: Estimation and predicition in a multistate model for breast cancer. Biom J. 2006, 48: 366380. 10.1002/bimj.200510218.View ArticlePubMedGoogle Scholar
 Andersen PK, Keiding N: Multistate models for event history analysis. Stat Methods Med Res. 2002, 11: 91115. 10.1191/0962280202SM276ra.View ArticlePubMedGoogle Scholar
 Broët P, Scholl SM, Fourquet A, De Rycke Y, Pouillart P, Mosseri V, Asselain B, de la, Rochefordière A: Analyzing prognostic factors in breast cancer using a multistate model. Breast Cancer Res Treat. 1999, 54: 8389. 10.1023/A:1006197524405.View ArticlePubMedGoogle Scholar
 Beyersmann J, Wolkewitz M, Allignol A, Grambauer N, Schumacher M: Application of multistate models in hospital epidemiology: Advances and challenges. Biom J. 2011, 53: 332350. 10.1002/bimj.201000146.View ArticlePubMedGoogle Scholar
 Schmoor C, Olschewski M, Sauerbrei W, Schumacher M: Longterm followup of patients in four prospective studies of the german breast cancer study group (GBSG): a summary of key results. Onkologie. 2002, 25: 143150. 10.1159/000055224.PubMedGoogle Scholar
 Fisher B, Anderson S, Fisher ER, Redmond CK, Wickerham DL, Wolmak N, Mamounas EP, Deutsch M, Margolese RG: Significance of ipsilateral breast tumor recurrence after lumpectomy. The Lancet. 1991, 338: 327331. 10.1016/01406736(91)904755.View ArticleGoogle Scholar
 Veronesi U, Marubini E, Del Vecchio M, Manzari A, Andreola S, Greco M, Luini A, Merson M, Saccozzi R, Rilke F, Salvadori B: Local recurrences and distant metastases after conservative breast cancer treatments: partly independent events. J Natl Canc nst. 1995, 87: 1927. 10.1093/jnci/87.1.19.View ArticleGoogle Scholar
 Arriagada R, Rutqvist LE, Mattsson A, Kramar A, Rotstein S: Adequate locoregional treatment for early breast cancer may prevent secondary dissemination. J Clin Oncol. 1995, 13: 28692878.PubMedGoogle Scholar
 Dunst J, Steil B, Furch S, Fach A, Lautenschläger C, Diestelhorst A, Lampe D, Kölbl H, Richter C: Prognostic significance of local recurrence in breast cancer after postmastectomy radiotherapy. Strahlenther Onkol. 2001, 10: 504510.View ArticleGoogle Scholar
 Schmoor C, Sauerbrei W, Bastert G, Schumacher M: Role of isolated locoregional recurrence of breast cancer: Results of four prospective studies. J Clin Oncol. 2000, 18: 16961708.PubMedGoogle Scholar
 Cortese G, Andersen PK: Competing risks and timedependent covariates. Biom J. 2009, 51: 138158.Google Scholar
 Klein J P Keiding, Copelan E: Plotting summary predicitions in multistate survival models: probabilities of relapse and death in remission for bone marrow transplantation patients. Stat Med. 1993, 12: 23152332. 10.1002/sim.4780122408.View ArticlePubMedGoogle Scholar
 Therneau TM, Grambsch PM: Modelling Survival Data  Extending the Cox Model. New York: Springer, 2000 Statistics for Biology and Health
 Cox DR: Regression models and lifetables (with discussion). J R Stat Soc B. 1972, 34: 187220.Google Scholar
 Beyersmann J, Schumacher M: Timedependent covariates in the proportional subdistribution hazards model for competing risks. Biostatistics. 2008, 9: 765776. 10.1093/biostatistics/kxn009.View ArticlePubMedGoogle Scholar
 Royston P, Altman DG: Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling (with discussion). Appl Stat. 1994, 43: 429467. 10.2307/2986270.View ArticleGoogle Scholar
 Royston P, Sauerbrei W: Multivariable ModelBuilding: A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables. Chichester: Wiley, 2008 Wiley Series in Probability and Statistics
 Sauerbrei W, Royston P: Building multivariable prognostic and diagnostic models: transformations of the predictors by using fractional polynomials. J R Stat Soc A. 1999, 162: 7194. 10.1111/1467985X.00122. [Corrigendum (2002) 165 399–400]View ArticleGoogle Scholar
 Berger U, Schäfer J, Ulm K: Dynamic Cox modelling based on fractional polynomials: timevariations in gastric cancer prognosis. Stat Med. 2003, 22: 11631180. 10.1002/sim.1411.View ArticlePubMedGoogle Scholar
 Buchholz A, Sauerbrei W: Comparison of procedures to assess nonlinear and timevarying effects in multivariable models for survival data. Biom J. 2011, 53: 308331. 10.1002/bimj.201000159.View ArticlePubMedGoogle Scholar
 Akaike H: A new look at the statistical model identification. IEEE Trans Automatic Control. 1974, 19: 716723. 10.1109/TAC.1974.1100705.View ArticleGoogle Scholar
 R project: The R project for Statistical Computing. [http://www.rproject.org]
 Harrell FE: Package ’rms’. Rproject 2013. 2009, [http://www.rproject.org]Google Scholar
 Therneau TM: Package ‘survival’: Survival Analysis, Including Penalised Likelihood. Rproject. 2011, [http://www.rproject.org]Google Scholar
 Harrell FE: Regression Modeling Strategies. New York: Springer, 2001 Springer Series in Statistics
 Sauerbrei W, MeierHirmer C, Benner A, Royston P: Multivariable regression model building by using fractional polynomials: Description of SAS, STATA and R programs. Comput Stat Data Anal. 2006, 50: 34643485. 10.1016/j.csda.2005.07.015.View ArticleGoogle Scholar
 Kalbfleisch JD, Prentice RL: The Statistical Analysis of Failure Time Data. 1980, New York: Wiley, Wiley Series in Probability and Mathematical StatisticsGoogle Scholar
 Klein JP, Shu Y: Multistate models for bone marrow transplantation studies. Stat Methods Med Res. 2002, 11: 117139. 10.1191/0962280202sm277ra.View ArticlePubMedGoogle Scholar
 Charman SA, Sharples LD, McNeil KD, Wallwork J: Assessment of survival benefit after lung transplantation by patient diagnosis. J Heart Lung Transplant. 2002, 21: 226232. 10.1016/S10532498(01)003527.View ArticlePubMedGoogle Scholar
 De Meester J, Smits JMA, Persijn GG, Haverich A: Listing for lung transplantation: Life expectancy and transplant effect, stratified by type of endstage lung disease, the Eurotransplant experience. J Heart Lung Transplant. 2001, 20: 518524. 10.1016/S10532498(01)002418.View ArticlePubMedGoogle Scholar
 Frydman H, Szarek M: Nonparametric estimation in a Markov ‘illnessdeath’ process from interval censored observations with missing intermediate transition status. Biometrics. 2009, 65: 143151. 10.1111/j.15410420.2008.01056.x.View ArticlePubMedGoogle Scholar
 Frydman H, Szarek M: Estimation of overall survival in an ‘illnessdeath’ model with application to the vertical transmission of HIV1. Stat Med. 2010, 29: 20452054. 10.1002/sim.3949.View ArticlePubMedGoogle Scholar
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14712288/13/80/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.