Abstract
Background
Multistate models have become increasingly useful to study the evolution of a patient’s state over time in intensive care units ICU (e.g. admission, infections, alive discharge or death in ICU). In addition, in criticallyill patients, data come from different ICUs, and because observations are clustered into groups (or units), the observed outcomes cannot be considered as independent. Thus a flexible multistate model with random effects is needed to obtain valid outcome estimates.
Methods
We show how a simple multistate frailty model can be used to study semicompeting risks while fully taking into account the clustering (in ICU) of the data and the longitudinal aspects of the data, including left truncation and right censoring. We suggest the use of independent frailty models or joint frailty models for the analysis of transition intensities. Two distinct models which differ in the definition of time t in the transition functions have been studied: semiMarkov models where the transitions depend on the waiting times and nonhomogenous Markov models where the transitions depend on the time since inclusion in the study. The parameters in the proposed multistate model may conveniently be computed using a semiparametric or parametric approach with an existing R package FrailtyPack for frailty models. The likelihood crossvalidation criterion is proposed to guide the choice of a better fitting model.
Results
We illustrate the use of our approach though the analysis of nosocomial infections (ventilatorassociated pneumonia infections: VAP) in ICU, with “alive discharge” and “death” in ICU as other endpoints. We show that the analysis of dependent survival data using a multistate model without frailty terms may underestimate the variance of regression coefficients specific to each group, leading to incorrect inferences. Some factors are wrongly significantly associated based on the model without frailty terms. This result is confirmed by a short simulation study. We also present individual predictions of VAP underlining the usefulness of dynamic prognostic tools that can take into account the clustering of observations.
Conclusions
The use of multistate frailty models allows the analysis of very complex data. Such models could help improve the estimation of the impact of proposed prognostic features on each transition in a multicentre study. We suggest a method and software that give accurate estimates and enables inference for any parameter or predictive quantity of interest.
Background
Multistate models have become increasingly useful to understand complicated biological processes and to evaluate the relations between different types of events. These methods have been developed to study simultaneously several competing causes of failure (e.g. competing risks of death) or to study the evolution of a patient’s state over time (e.g. admission in intensive care units (ICU), infections, alive discharge or death in ICU) and the focus is in the process of going from one state to another.
Furthermore, many studies include clustering of survival times. For instance, in critically ill patients, data come from different ICUs and because observations are clustered into groups (or units), the observed outcomes cannot be considered as independent. Thus a flexible multistate model with random effects is needed to obtain valid outcome estimates. Ignoring the existence of clustering may underestimate the variance of regression coefficients specific to each group, leading to incorrect inferences.
Ripatti et al. [1] proposed a threestate frailty model to model age at onset of dementia and death in Swedish twins. The intrapairs correlations and the other parameters were estimated using hierarchical bayesian model formulation and Gibbs sampling, both of which can be timeconsuming. Katsahian et al. [2,3] extended Fine and Gray’s [4] model to the case of clustered data, by including random effects in the subdistribution hazards. They first used the residual maximum likelihood then the penalized partial loglikelihood to estimate the parameters. However, the estimation approach does not directly yield estimators of the transition intensities, which often have a meaningful interpretation in epidemiological studies. Most of the time, the baseline intensity estimate is based on Breslow’s estimate leading to a piecewiseconstant baseline hazard function or unspecified baseline hazard function.
In this paper, we show how a simple multistate frailty model can be used to study semicompeting risks [5] while fully taking into account the clustering (in ICU) of the data and the longitudinal aspects of the data, including left truncation and right censoring. We suggest the use of independent frailty models or joint frailty models for the analysis of transition intensities. This approach is of interest for several reasons. First, it allows to deal with heterogeneity between ICUs. We do this by including clusterspecific random effects or frailties in the multistate model. Frailties represent the unmeasured covariates at the cluster level, which may affect the rate of occurrence of each of the events. Moreover, this approach allows us to evaluate different prognostic effects of covariates on each transition probability. Another interesting and perhaps underrated advantage of multistate models is the possibility to use them to predict clinical prognosis whereby a patient will be in a given health state at time u given a particular history at time t. This work extends previous research by dealing with clustered competing risks and by giving smoothed estimates of the transition rates. In addition the joint approach allows the analysis of two processes that evolve with time leading to more accurate estimates.
Two distinct approaches are often used in multistate models. They differ in the definition of time t in the transition functions. In the first approach, the transition probability between two states depends only on the waiting times, the clock is reset to zero every time a patient enters a new state and a semiMarkov model is used. In the second approach, the transition depends only on the time since inclusion in the study, and nonhomogenous Markov models are used. The proposed approach can deal with both situations and is illustrated in this article. The choice between one the two approaches depends on the clinical knowledge of the events of interest. If it is expected that the transition probability (for instance toward death) will not change as a function of the time since randomization or inclusion, the analysis can be based solely on the semiMarkov model and it can thus be studied how the transition probability evolves after an event has taken place (for instance after nosocomial infections). However, if it is difficult to choose clinically between the semiMarkov or the nonhomogeneous Markov approach, one can use statistical criteria [6].
As discussed in the section on the estimation procedure method, an important advantage of our proposed approach is that the parameters in the multistate model may conveniently be computed using a semiparametric or parametric approach with an existing R package FrailtyPack for frailty models.
The paper is organized as follows. First, the ICU data is briefly presented. The next section describes the statistical multistate frailty model for clustered data with estimation procedure. Then, the model is applied to the analysis of nosocomial infections (ventilatorassociated pneumonia infections) in ICUs, with “alive discharge” and “death” in ICU as other endpoints. Results from a simulation study are reported. Finally, a concluding discussion is presented.
Methods
Motivating example
Data Source
We conducted a prospective observational study using data from the multicenter Outcomerea database between November 1996 to April 2007. The database contains data from 16 French ICUs, among which data on admission features and diagnosis, daily disease severity, iatrogenic events, nosocomial infections, and vital status. Every year, the data of a subsample of at least 50 patients per ICU were entered in the database; patients had to be older than 16 years and to have stayed in ICU for more than 24 hours. To define this random subsample, each participating ICU selected either consecutive admissions to specific ICU beds throughout the year or consecutive admissions to all ICU beds over a single month.
Data collection
Database quality measures were taken such as the continuous training of investigators in each ICU or regular data quality checks (see [7]). A one day coding course was organized annually with the study investigators and research monitors. In all ICUs, as previously reported [8,9], VAP was suspected based on the development of persistent pulmonary infiltrates on chest radiographs combined with purulent tracheal secretions, and/or body temperature ≥ 38.5°C or ≤ 36.5°C, and/or peripheral blood leukocyte count ≥ 10 × 10^{9} / L(Giga/liter) or ≤ 4 × 10^{9} / L. The definite diagnosis of VAP required a positive culture result from a protected specimen brush (≥ 10^{3}cfu/ml), plugged telescopic catheter specimen (≥ 10^{3}cfu/ml), BAL fluid specimen (≥ 10^{4}cfu/ml), or quantitative endotracheal aspirate (≥ 10^{5}cfu/ml). Investigators systematically performed bacteriological sampling before changing antimicrobials.
Study population
We considered death in ICU and discharge to be absorbing state and VAP as a non absorbing state. Patients were included in the study if they had stayed in the ICU for at least 48 hours and had received mechanical ventilation (MV) within 48 hours after ICU admission. We obtained 2871 patients, corresponding to 37395 ICU days. The median MV duration was 7 days with inter quartile range (IQR = [413]).
The multistate approach and estimation
Multistate model
We consider the multistate model represented in Figure 1 corresponding to our motivating example. This discrete stochastic process (X(t))_{t ≥ 0} with state space {0,1,2,3} called disability model is often used to describe the occurrence of nosocomial infections in ICU [10]. Here we focus on the occurrence of VAP modelled via the transition 0 (admission in ICU without VAP) into state 1 (VAP infections). Discharge alive and death in ICU are modelled via transitions into states 3 and 2, respectively. A transition will be simply denoted hk. The set of transitions for the disability model will be denoted ( ). The distribution of this multistate process is characterized by the transition intensities:
where are the transition probabilities. The filtration stands for the history of the process just before time t (t^{−}).
Figure 1. The “disability Model”.
From the definition of the transition intensities, we focus on two kinds of models. First, we consider the nonhomogeneous Markov model implying that the transition intensities depend only on the current time t. In our application, the time t represent time since entry in ICU (and for all patients X(0) = 0). So for instance, when death occurs after VAP infections, t represents time since ICU entry, and not time since infections. We obtain:
depends on t but not on the entry time into state 1 (noted T_{1} in the sequel). Second, we define the semiMarkov model implying that the transitions intensities depend only on the time spent in the current state. In particular, we get in our application:
Intensity functions with a shared frailty term
Take the case of a study with Gindependent intensive care units (i = 1,…,G). Let be a vector of covariates (specific to the transition hk) for subject j (j = 1,…,n_{i}) from group i. In order to take into account the heterogeneity of the population due to the fact data come from different ICUs, we simply model the transition intensities by a proportional intensities model with frailty. For any , the transition intensity for the jth patient in the ith ICU with random ICU effects given a vector of covariates is defined by
where is the baseline transition intensity specific to the transition hk (with specific parameter θ_{hk} for parametric model) and the random ICU effects are also specific to the transition hk, independent and follow a Gamma distribution ( , , ). The variance γ_{hk} of the represents the heterogeneity between ICU of the overall underlying baseline risk for the transition hk.
Remark
this definition is correct for the nonhomogeneous Markov model. We get a similar definition for the semiMarkov model by replacing the current time t in the equation (4) by the time spent in the current state t − T_{1}(for the transition 12 and 13).
Intensity functions with a joint frailty term
In the model (4) we assume that the different times to transitions are independent. In some cases this assumption may be violated, for instance in our motivating example, the transition times to death with the VAP and the transition times to discharge with VAP may be dependent. This dependency should be accounted for in the joint modelling of these two survival endpoints. There can be many reasons to use joint models of two survival endpoints, including giving a general description of the data, correcting for bias in survival analysis due to dependent dropout or censoring, and improving efficiency of survival analysis due to the use of auxiliary information [11].
In this work we will thus also consider some transitions jointly, in a joint frailty model setting as follows:
The random effects ω_{j} (frailties) are assumed independent. Mainly for reasons of mathematical convenience, the frailty terms are often assumed to follow a gamma distribution. The gamma frailty density is adopted here with unit mean and variance η. This choice and other possibilities such as lognormal, positive stable distributions are discussed in several papers [12,13]. The common frailty parameter ω_{j} will take into account the heterogeneity of the data, associated with unobserved covariates.
In the traditional model, the assumption is that ζ = 0 in (5), that is α_{hk}(t) does not depend on , and thus the two intensity functions α_{hk}(t) and are not associated, conditional on covariates.When ζ = 1, the effect of the frailty is identical for the two transition times. When ζ > 0, the two transition times are positively associated; higher frailty will result for instance in a higher risk of discharge and a higher risk of death; while ζ < 0 implies a negative association. This means that unobserved individual factors produce higher transition rates to death and smaller transition rates to discharge, or inversely. We can think that for sicker patients, the mortality will be high but with a low discharge rate, conversely for healthy patients, the discharge rates will be high with a low death rate. The interpretation of the value of ζ only makes sense in case of heterogeneity, i.e. when the variance of the random effects is significantly different from zero. However, in this model we assume that there is no intracluster correlation anymore after having taken into account prognostic factors and after adjusting for a subject specific random effect term.
Estimation procedure
First, in our study, we consider that the process (X(t))_{t ≥ 0} is observed in continuous time, i.e., we know at each time t the state of the process for each subjects. Secondly, for the model with shared frailties terms, we specify that each transition intensity has its own set of parameters. In other words, the parameter θ_{hk}, β_{hk} and γ_{hk} are specific to the transition hk. Such an assumption is common when dealing with multistate models and [14] has shown that this assumption allows us to consider the problem of estimating (parametrically or semiparametrically) the transition intensities separately (that is transition by transition). Thus, each transition intensity can be evaluated by estimation methods used in survival analysis in the presence of rightcensored data only (for instance when semiMarkov models are used) or in the presence of rightcensored and lefttruncated data (for instance when nonhomogeneous Markov models are used, the transition intensities 12 and 13 are evaluated using delayed entry).
In this paper, we use two different approaches. First, we use a parametric model for the multistate model specifying each baseline transition intensity by a Weibull distribution ( ; θ_{hk} = (a_{hk}, b_{hk})). For each transition, we estimate the different parameter (ζ_{hk} = (θ_{hk}, β_{hk}, γ_{hk})) by maximizing the full marginal loglikelihood [15]. In this case, the full log likelihood for rightcensored data and lefttruncated data takes a simple form with an analytical solution for the integrals [15]. In practice, we use the FrailtyPack R package [16] to estimate all parameters in the different multistate model performed [17]. Secondly, we consider a semiparametric estimation by introducing a semiparametric penalized likelihood approach to estimate the different parameters β_{hk}, γ_{hk} and the baseline intensity function α_{0,hk}(t). This approach is more flexible than a too rigid parametric approach to estimate a smooth baseline intensity function. To obtain a smooth baseline intensity function, we penalize the likelihood by the squared norm of the second derivative of the intensity function. The penalized loglikelihood for the transition hk is thus defined as
where is l(β_{hk}, γ_{hk}, α_{0,hk}(·)) the full log likelihood, is the second derivate of the baseline intensity function, and κ_{hk} is a positive smoothing parameter that controls the tradeoff between the data fit and the smoothness of the functions. Maximization of (6) defines the maximum penalized likelihood estimators (MPnLEs). The estimator (MPnLEs) cannot be calculated explicitly but can be approximated on the basis of splines and the smoothing parameters κ_{hk} can be chosen by maximizing a likelihood crossvalidation criterion as in Joly et al. [18]. The maximum penalized likelihood method is also implemented in the FrailtyPack R package [16].
Remark
In the presence of intensity functions with a joint frailty term in the model, a maximum penalized likelihood estimation is also used [11] and implemented in the FrailtyPack R package.
Model choice
Several models with different approaches have been defined: a semiMarkovian model (noted ( )) or a nonhomogeneous Markov model (noted ( )) with a totally parametric approach or a semiparametric approach. For example, we define and representing respectively the estimators of ( ) and ( ) based on a sample of i.i.d. observations which come from the true unknown distribution P^{∗}. To assess the risk of each estimator we propose to use the Expected KullbakLeibler risk (EKL) [6] defined as:
where is the loglikelihood ratio and is a new i.i.d observation. Commenges et al. [6] proposed the likelihood crossvalidation (LCV) criterion to estimate this risk:
where represents the likelihood of the observation based on the estimator defined without this observation.
Finally, to choose between the two estimators and , we compute some differences of Expected KullbakLeibler risk (based on difference of LCV criterion). Commenges et al. [6] give an interpretation of the magnitude of these risks: values of 10^{−1}, 10^{−2}, 10^{−3}, 10^{−4} are respectively qualified “large”, “moderate”, “small” and “negligible”. The FrailtyPack R package provides the LCV criterion for survival analysis (for us a transition hk). As we can estimate each transition separately, it is straightforward to get the LCV criterion for a particular Multistate model. Concerning the selection of the covariates in each model, we use a manual selection procedure motivated by epidemiological/clinical knowledge and also based on statistical significance of hazard ratios.
Prediction
A posterior probability of any event can be easily derived from the multistate model and its parameters. This probability, which can be computed for a new subject using a given set of covariates at the current time, and given a postinclusion events, constitutes a dynamic tool of prediction. For instance, the aim may be to predict if VAP occurs between time t^{∗} and horizon time t^{∗} + h. The posterior probability of developing VAP on [t^{∗};t^{∗} + h] given the random effects and the covariates can be easily computed using the following expression with k = 1:
with,
The estimated posterior probabilities, can be obtained by substituting the maximum penalized likelihood estimates of parameters β_{0k}, γ_{0k}, α_{0,0k}, the empirical Bayes estimates for and the individual information for the covariates by equation (7).
Results
Application revisited
The methodology exposed in the previous section is now applied to the OUTCOMEREA database. Table 1 gives first a description of the number of patients for each transition. The study initially included 2438 patients from 16 hospitals (between 3 and 1027 patients by centre). Patients were mostly discharged alive (69%), but there are also individuals who died in ICU (16%) or developed a nosocomial infection, a VAP (15%). The database contains many covariates that are related with the states modelled. Based on our previous studies [8,9,19,20] and literature reviews [21], we selected timefixed covariates candidates as potential predictors of VAP, discharge and death in the ICU. A classical descending method was applied to select prognostic factors in each model (without frailties). This backward or descending elimination technique begins by calculating statistics for a model which includes all of the independent variables. Then the variables are deleted from the model one by one until all the variables remaining in the model produce statistics significant at the 0.10 level. At each step, the variable showing the smallest contribution to the model is deleted. In this application to assist in the variables interpretation, continuous covariates “age” and “Simplified Acute Physiology Score” (SAPS) were respectively transformed as categorial variables by median or quartiles. The different models were estimated using the FrailtyPack package developed by [16]. The Likelihood cross validation criterion is used to compare the different models. According to this criterion, the best is a semiMarkov model estimated by a semiparametric approach (penalized likelihood estimator) with LCV = 4.044. Note that the difference of LCV between the semiparametric estimator from the semiMarkov model and from the nonhomogeneous Markov model (LCV = 4.055) is equal to 0.011 qualified as “moderate”. We have also compared a parametric approach versus a semiparametric approach: for the nonhomogenous model we obtained a difference of LCV equal to 0.407 qualified as “Large” between the parametric estimation (LCV = 4.462) and the semiparametric estimation (LCV = 4.055); for the semiMarkov model we obtained a difference equal to 0.416 between the parametric estimation (LCV = 4.460) and the semiparametric estimation (LCV = 4.044).
Table 1. Number of patients from the OUTCOMEREA database present in each transition
Table 2 presents the semiparametric estimation results of the semiMarkov model. A complete description of these selected covariates in this model is provided elsewhere ( [7]).
Table 2. Result of the semimarkov model with penalized likelihood estimation incorporating frailty terms: HR (exp(β)) estimates and corresponding confidences intervals for the different transition intensities
We also fitted a semiMarkov model without taking into account the intracentre correlation (results not shown). Using this model some factors were wrongly significantly associated. For instance, the effects of coma (Relative Risk= 1.39 (95%CI 1.071.79)), shock (Relative Risk= 1.28 (95%CI 1.011.62)) and the presence of an enteral nutritional therapy (Relative Risk= 1.24 (95%CI 1.001.53)) were incorrectly observed as significantly associated with the risk of VAP using the multistate model without frailty term.
The variance of the frailty is significantly different from 0 for the transitions 01, 03 and not significantly different from 0 for the other transitions (Table 2). This means a significant heterogeneity of the transition rates (01, 03) between centres, but we did not observe a betweencentre heterogeneity in the risk of death among patients with or without a VAP, nor for the risk of discharge with VAP. This is also depicted with the posterior frailty mean estimation for each centre in Figure 2. The number of subjects recruited greatly varied between centres (between 3 and 1027 patients by centre) therefore we also conducted an analysis excluding the two smallest hospitals (with 3 and 7 patients in each). The findings were very similar in terms of parameter’ estimates, and variance components.
Figure 2. Random centrespecific effect as estimated by posterior distribution for the different transitions. The clusters have been ordered (ascending order) by the number of subjects in the cluster.
We previously evaluated the heterogeneity between centres, but considering a different random effect for each transition of the multistate model. We also fitted a joint frailty model for the two transitions 12 and 13, with a shared subjectspecific random effect for the two transitions [11]. This approach allows us to simultaneously evaluate the prognostic effects of covariates on the two survival endpoints, discharge or death with a VAP. This joint frailty model accounts for the dependency between the two outcomes, and corrects for bias in the analysis due to dependency. Results are exposed at the bottom of Table 2. When comparing the joint model for transitions 12 and 13 to the reduced shared frailty models for the same transitions, covariates effects are similar, while the hazard ratio of SAPSII is slightly greater in the joint frailty model. This illustrates that ignoring the dependence between time to death and time to discharge may result in biases in the reduced shared frailty model compared to the joint model.
The variance of the joint random effect ( , onesided Wald test = 0.77/0.13 = 5.92 > 1.64) is significantly different from 0 but the power coefficient ζ is not significantly different from 0 (twosided Wald test = 0.16/0.29 = 0.55 < 1.96). This shows that times of deaths and discharge are correlated, but with a slightly nonsignicant negative association ( ).
To illustrate the use of posterior probabilities of VAP between times t∗ and t∗ + h given the information collected until time t∗, we show in Figure 3 the predicted probabilities for a patient in state 0 of developing VAP between times t∗ and t∗ + 3 days for patients of the same profile but belonging to each of the 16 centres. We observe a large difference of VAP probability according to the centre. These individual predictions underline the usefulness of dynamic prognostic tools that can take into account the clustering of observations.
Figure 3. Predicted probabilities (see formula (7)) of a patient in state 0 developing VAP betweent∗andt∗ + 3days for a patient from the 16 different centers. Predictions given for a men, aged > 62, SAPSII<33, with medicine admission, no chronic diseases, no diabetes, no ARDS, no Trauma, with shock, no acute respiratory failure, no antimicrobials, with inotropes, with enteral and parenteral nutrition.
A Simulation study
Simulation details
In this section, we present a short simulation study, which in particular reveals the importance of taking into account the heterogeneity of the population for estimating the different intensity transition. We have considered the semiMarkov multistate model (with four states) described in Figure 1. Each transition intensity is modelled by a proportional intensities model with a shared frailty term as defined in (4): the baseline intensity corresponds to the one of a Weibull distribution, the random effects follow a Gamma distribution and two covariates are used. These same two covariates are used for all the transitions: a binary variable modelled by a Binomial distribution with probability 0.5 (with specific effect ) and a continuous variable (with specific effect ) modelled by a Normal distribution (μ = 2 and σ^{2}=0.1). The continuous variable is specific to the cluster; i.e. the same value for each subject belonging to the same cluster. There were either G=20, 30, 50 or 100 clusters with n_{i} = 75 or 150 subjects in each cluster. Two different values (0.15 and 0.30) of the variance of the frailty γ_{hk} are used for each transition. There were 500 simulated data sets for each case. The complete description of the simulation parameters is provided in Table 3.
Table 3. Simulation parameters of the semimarkov model
Briefly, a simulation of a semiMarkov model consists of: (i) for each subject generate 3 times T_{01},T_{02} and T_{03} (representing respectively the occurrence of VAP, Death or Discharge) according to the intensities transitions defined in (4); (ii) if T_{01} ≠ min(T_{01},T_{02},T_{03}) then the subject dies at time T_{02} (if T_{02} = min(T_{01},T_{02},T_{03})) or transits into the state Discharge (if T_{03} = min(T_{01},T_{02},T_{03})) and then back to (i) for a new subject; (iii) if T_{01} = min(T_{01},T_{02},T_{03}) then the subject contracts VAP at time T_{01}and then generates 2 new times (T_{12},T_{13}) representing respectively the occurrence of Death or Discharge for a subject with VAP; (iv) if T_{12} = min(T_{12},T_{13}) then the subject die at time T_{01} + T_{12} else the subject transits in the state Discharge at time T_{01} + T_{13}; and then back to (i) for a new subject.
For simulation run, we estimated two parametric semimarkov models: 1) with specific frailty term in each transition and 2) without frailty term. For the two models, we computed the mean, the empirical standard errors (SEs), i.e. the SEs of estimates and the mean of the estimated SEs for , and .
Simulation results
The results of the simulation studies using parametric estimation are summarized in Table 4 and 5. To shorten the presentation of the results, we only present the case with G = 30 and 100 clusters. In general for the multistate model incorporating frailty terms, the fixed effect and the variance γ_{hk}of the frailties are well estimated. These bias decrease with the increasing sample size. The variability of estimates of β_{hk}, measured by the empirical SEs or the mean of the estimated SEs, are similar. We observe the same result for the variability of estimates of γ_{hk}. However, if we omit the frailties term in the estimation of the semiMarkov multistate model, the estimators of the fixed effect β_{hk} are biased and the variability of estimates of β_{hk}, measured by the mean of the estimated SEs, is lower than the variability measured by the empirical SEs. These results are clearly highlighted for the variable specific (representing by the effect β_{2}) to the cluster (i.e. the same value for each subject belonging to the same cluster).
Discussion
We have described a multistate model with frailty terms to account for heterogeneity between clusters on each transition. Such models appear promising in the setting of competing risk analyses using clustered data (i.e., multicentre clinical trials, metaanalysis). Lack of software is a potential obstacle. We propose here a tractable model, semiMarkov as well as nonhomogenous Markov, with semiparametric or parametric estimates. The model can be readily derived with the R package FrailtyPack a simple and free approach, which does not require any timeconsuming calculation. We also proposed a measure of models selection which evaluates the relative goodness of fit among a collection of models. To give an example, we provided a R code for simulating a data set and to analyze it (see Additional file 1).
Additional file 1. We provide a R code for simulating a data set and to analyze it: file name “supplementarymaterialRcode.R”.
Format: R Size: 11KB Download file
Vital status and time of death or time of discharge in ICU are known exactly. However, there may be more complex schemes with intervalcensored times to events, i.e., the event occurs in a known time interval L,R. The semiMarkov or nonhomogenous Markov models discussed in this paper do not allow the direct treatment of these intervalcensored data. It would be interesting to extend the multistate frailty model in the case of intervalcensoring [22]. This would lead to numerical integration in the estimation process due to the lack of a closedform solution of the multiple integrals, and this can be very timeconsuming especially when the number of states rises. Also for large numbers of states, it is clear that substantial datasets are required for frailty variance estimation.
The proposed approach can also be used to predict probabilities of future events, given a patient’s history, covariates, and random effets, using parameter estimates and the estimates of corresponding baseline hazards and survival functions. Open research questions include prediction assessment with timedependent prognostic factors. The aim would be to develop an updating mechanism which would allow dynamic updating of the predictions for a given patient in case of important changes in biomarkers.
A recent article discussed the identifiability and the (im)possibilities of frailties in multistate models [23] but without considering covariates. They also compared predictive accuracy of different multistate models with or without frailties using a kfold cross validated partial loglikelihood [24]. They obtained that frailty models showed the best predictive performance in the comparison.
VAP represents an important and challenging example in which multistate frailty models should be used. This nosocomial infection is very frequent in ICU and is associated with an increase in ICU mortality, length of stay and cost. Many risk factors have been described in the past, and some new preventive interventions have been tested with often conflicting results [25]. This result could be due to the absence of control of appropriate confounding factors, to the heterogeneity of effects according to ICU subpopulations or to discrepancies in the diagnostic procedure. Indeed the diagnosis is based on the physician’s behaviour when clinical, biological and radiological signs of pneumonia occur. Even when diagnostic criteria are fixed in randomized control trials, the reported incidences vary from nine to 70%. Even if all fixed covariates are taken into account, and definitions carefully followed, the incidence densities still vary from 9.7 to 26.1 per 1000 mechanicalventilation days [26]. The use of frailty terms in multistate models might then be important to unmask residual sources of variability (centre effect) when looking for risk factors of VAP. It may also be useful to compare incidences of VAP between hospitals if used as a quality indicator. It may also be important to explain the huge variability in the estimation of the attributable mortality of the disease.
Since the multistate model under consideration contains clusterspecific random effects, the definition of the predictions is not straightforward. The proposed posterior prediction probabilities may be used to predict survival functions of subjects from existing clusters. In this cluster focus, the random effects are themselves of interest, and they are parameters to be estimated. In contrast, when the interest is on predicting survival for patients from new clusters, a marginal approach is better and corresponds to a population focus rather than a cluster focus. In this population focus, each conditional survival function is replaced by its marginal version. The marginal or observable survival function for the shared gamma frailty model (4) is with the cumulative transition intensity function specific to the transition hk.
Conclusions
The use of multistate frailty models allows the simple analysis of very complex data. Such models could help improve the estimation of the impact of proposed prognostic features on each transition in a multicentre study. We have suggested a method and software that gives accurate estimates and enables inference for any parameter or predictive quantity of interest.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
BL and VR developed the methodology, performed the simulation and the analysis on the dataset as well as wrote the manuscript. JFT collected and interpreted the dataset as well as wrote the description and the interpretation on the Application section. All authors read and approved the final manuscript.
Acknowledgements
This work was supported by the ANR grant 2010 PRSP 006 01 for the MOBIDYQ project (Dynamical Biostatistical models). We would like to thank the members of the Outcomerea Study Group for sharing their database.
References

Ripatti S, Gatz M, Pedersen N, Palmgren J: Threestate frailty model for age at onset of dementia and death in Swedish twins.
Genet Epidemiol 2003, 24(2):139149. PubMed Abstract  Publisher Full Text

Katsahian S, Boudreau C: Estimating and testing for center effects in competing risks.
Stat Med 2011, 30(13):16081617. PubMed Abstract  Publisher Full Text

Katsahian S, RescheRigon M, Chevret S, Porcher R: Analysing multicentre competing risks data with a mixed proportional hazards model for the subdistribution.
Stat Med 2006, 25(24):42674278. PubMed Abstract  Publisher Full Text

Fine J, Gray R: A proportional hazards model for the subdistribution of a competing risk.
J Am Stat Assoc 1999, 94:496509. Publisher Full Text

Commenges D, Joly P, GegoutPetit A, Liquet B: Choice between semiparametric estimators of Markov and nonMarkov multistate models from coarsened observations.
Scand J Stat 2007, 34:3352. Publisher Full Text

NguileMakao M, Zahar J, Français A, Tabah A, GarrousteOrgeas M, Allaouchiche B, GoldgranToledano D, Azoulay E, Adrie C, Jamali S, et al.: Attributable mortality of ventilatorassociated pneumonia: respective impact of main characteristics at ICU admission and VAP onset using conditional logistic regression and multistate models.
Intensive Care Med 2010, 36(5):781789. PubMed Abstract  Publisher Full Text

Clech́ C, Timsit J, Lassence A, Azoulay E, Alberti C, GarrousteOrgeas M, Mourvilier B, Troche G, Tafflet M, Tuil O, et al.: Efficacy of adequate early antibiotic therapy in ventilatorassociated pneumonia: influence of disease severity.
Intensive Care Med 2004, 30(7):13271333. PubMed Abstract  Publisher Full Text

Moine P, Timsit J, De Lassence A, Troche G, Fosse J, Alberti C, Cohen Y: Mortality associated with lateonset pneumonia in the intensive care unit: Results of a multi center cohort study.
Intensive Care Med 2002, 28(2):154163. PubMed Abstract  Publisher Full Text

Beyersmann J, Wolkewitz M, Allignol A, Grambauer N, Schumacher M: Application of multistate models in hospital epidemiology: Advances and challenges.
Biometrical J 2011, 53(2):332350. Publisher Full Text

Rondeau V, MathoulinPelissier S, JacqminGadda H, Brouste V, Soubeyran P: Joint frailty models for recurring events and death using maximum penalized likelihood estimation: application on cancer events.
Biostatistics 2007, 8(4):708721. PubMed Abstract  Publisher Full Text

Hougaard P: Frailty models for survival data.
Lifetime Data Anal 1995, 1(3):255273. PubMed Abstract  Publisher Full Text

Pickles A, Crouchley R: A comparison of frailty models for multivariate survival data.
Stat Med 1995, 14(13):14471461. PubMed Abstract  Publisher Full Text

Hougaard P: Multistate Models: A Review.
Lifetime Data Anal 1999, 5:239264. PubMed Abstract  Publisher Full Text

Rondeau V, Commenges D, Joly P: Maximum Penalized Likelihood Estimation in a GammaFrailty Model.
Lifetime Data Anal 2003, 9:139153. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rondeau V, Mazroui Y, Diakite A, Gonzalez JR: FRAILTYPACK.
An R package for General frailty models using a semiparametrical penalized likelihood estimation or a parametrical estimation 2011.
[http://CRAN.Rproject.org/package=frailtypack webcite]. [R package version 2.220]

Rondeau V, Mazroui Y, Gonzalez J: FRAILTYPACK: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation.

Joly P, Commenges D, Letenneur L: A Penalized Likelihood Approach for Arbitrarily Censored and Truncated Data: Application to AgeSpecific Incidence of Dementia.
Biometrics 1998, 54:185194. PubMed Abstract  Publisher Full Text

Timsit JF, Zahar J, Chevret S: Attributable mortality of ventilatorassociated pneumonia.
Curr Opin Crit Care 2011, 17(5):464471. PubMed Abstract  Publisher Full Text

Bekaert M, Timsit JF, Vansteelandt S, Depuydt P, Vesin A, GarrousteOrgeas M, Decruyenaere J, Clec’h C, Azoulay E, Benoit D: OUTCOMEREA study group: Attributable Mortality of Ventilator Associated Pneumonia: A Reappraisal Using Causal Analysis.
Am J Respir Crit Care Med 2011, 184:11331139. PubMed Abstract  Publisher Full Text

Chastre J, Fagon J: Ventilatorassociated pneumonia.
Am J Respir Crit Care Med 2002, 165(7):867. PubMed Abstract  Publisher Full Text

Joly P, Commenges D, Helmer C, Letenneur L: A penalized likelihood approach for an illness–death model with intervalcensored data: application to agespecific incidence of dementia.
Biostatistics 2002, 3(3):433443. PubMed Abstract  Publisher Full Text

Putter H, van Houwelingen H: Frailties in multistate models: Are they identifiable? Do we need them?

Bøvelstad H, Nygård S, Størvold H, Aldrin M, Borgan Ø, Frigessi A, Lingjærde O: Predicting survival from microarray data–a comparative study.
Bioinformatics 2007, 23(16):20802087. PubMed Abstract  Publisher Full Text

Melsen W, Rovers M, Koeman M, Bonten M: Estimating the attributable mortality of ventilatorassociated pneumonia from randomized prevention studies.
Crit Care Med 2011, 39(12):273642. PubMed Abstract  Publisher Full Text

Zahar J, NguileMakao M, Français A, Schwebel C, GarrousteOrgeas M, GoldgranToledano D, Azoulay E, Thuong M, Jamali S, Cohen Y, de Lassence A, Timsit J: Predicting the risk of documented ventilatorassociated pneumonia for benchmarking: construction and validation of a score.
Prepublication history
The prepublication history for this paper can be accessed here: