- Research article
- Open access
- Published:

# Investigating hospital heterogeneity with a multi-state frailty model: application to nosocomial pneumonia disease in intensive care units

*BMC Medical Research Methodology*
**volume 12**, Article number: 79 (2012)

## 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 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 multi-state model with random effects is needed to obtain valid outcome estimates.

### Methods

We show how a simple multi-state frailty model can be used to study semi-competing 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: semi-Markov 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 multi-state model may conveniently be computed using a semi-parametric or parametric approach with an existing R package FrailtyPack for frailty models. The likelihood cross-validation 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 (ventilator-associated 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 multi-state 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 multi-centre 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 multi-state 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 three-state frailty model to model age at onset of dementia and death in Swedish twins. The intra-pairs correlations and the other parameters were estimated using hierarchical bayesian model formulation and Gibbs sampling, both of which can be time-consuming. 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 log-likelihood 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 piecewise-constant baseline hazard function or unspecified baseline hazard function.

In this paper, we show how a simple multi-state frailty model can be used to study semi-competing 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 cluster-specific random effects or frailties in the multi-state 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 multi-state 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 multi-state 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 semi-Markov 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 semi-Markov 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 semi-Markov 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 multi-state model may conveniently be computed using a semi-parametric 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 (ventilator-associated 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 multi-center **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 = [4-13]).

### The multi-state approach and estimation

#### Multi-state model

We consider the multi-state 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 \mathcal{E} (\mathcal{E}\phantom{\rule{0.3em}{0ex}}:=\{01,02,03,12,13\}). The distribution of this multi-state process is characterized by the transition intensities:

where {P}_{\mathrm{hk}}(t,t+\mathrm{\Delta t})=P\left(X\right(t+\mathrm{\Delta t})=k|X\left(t\right)=h;{\mathcal{F}}_{{t}^{-}}) are the transition probabilities. The filtration {\mathcal{F}}_{{t}^{-}} stands for the history of the process just before time *t* (*t*
^{−}).

From the definition of the transition intensities, we focus on two kinds of models. First, we consider the *non-homogeneous 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 *semi-Markov 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 *G*-independent intensive care units (*i* = 1,…,*G*). Let {\mathbf{Z}}_{\mathrm{hk}}^{\mathrm{ij}} 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 \mathrm{hk}\in \mathcal{E}, the transition intensity for the *j*th patient in the *i*th ICU with random ICU effects {u}_{\mathrm{hk}}^{i} given a vector of covariates {\mathbf{Z}}_{\mathrm{hk}}^{\mathrm{ij}} is defined by

where \left(\right.separators="">\n \n \n \n \alpha \n \n \n 0\n ,\n hk\n \n \n \n \n \theta \n \n \n hk\n \n \n \n \n (\n t\n )\n \n is the baseline transition intensity specific to the transition *hk* (with specific parameter *θ*
_{
hk
} for parametric model) and the random ICU effects {u}_{\mathrm{hk}}^{i} are also specific to the transition *hk*, independent and follow a Gamma distribution ({u}_{\mathrm{hk}}^{i}\sim \Gamma \left(\frac{1}{{\gamma}_{\mathrm{hk}}},\frac{1}{{\gamma}_{\mathrm{hk}}}\right), E\left[{u}_{\mathrm{hk}}^{i}\right]=1, \mathrm{Var}\left[{u}_{\mathrm{hk}}^{i}\right]={\gamma}_{\mathrm{hk}}). The variance *γ*
_{
hk
} of the {u}_{\mathrm{hk}}^{i} represents the heterogeneity between ICU of the overall underlying baseline risk for the transition *hk*.

##### Remark

this definition is correct for the *non-homogeneous Markov model*. We get a similar definition for the *semi-Markov 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 log-normal, 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 {\alpha}_{{\mathrm{hk}}^{\prime}}\left(t\right), and thus the two intensity functions *α*
_{
hk
}(*t*) and {\alpha}_{{\mathrm{hk}}^{\prime}}\left(t\right) 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 intra-cluster 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 multi-state models and [14] has shown that this assumption allows us to consider the problem of estimating (parametrically or semi-parametrically) 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 right-censored data only (for instance when *semi-Markov models* are used) or in the presence of right-censored and left-truncated data (for instance when *non-homogeneous 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 multi-state model specifying each baseline transition intensity by a Weibull distribution (\left(\right.separators="">\n \n \n \n \alpha \n \n \n 0\n ,\n hk\n \n \n \n \n \theta \n \n \n hk\n \n \n \n \n \n (\n t\n )\n \n =\n \n (\n \n \n a\n \n \n hk\n \n \n \n /\n \n \n \n b\n \n \n hk\n \n \n )\n \n \n \n (\n t\n /\n \n \n b\n \n \n hk\n \n \n )\n \n \n (\n \n \n a\n \n \n hk\n \u2212\n 1\n \n \n )\n \n \n \n; *θ*
_{
hk
} = (*a*
_{
hk
}, *b*
_{
hk
})). For each transition, we estimate the different parameter (*ζ*
_{
hk
} = (*θ*
_{
hk
}, *β*
_{
hk
}, *γ*
_{
hk
})) by maximizing the full marginal log-likelihood [15]. In this case, the full log likelihood for right-censored data and left-truncated 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 multi-state model performed [17]. Secondly, we consider a semi-parametric estimation by introducing a semi-parametric 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 log-likelihood for the transition *hk* is thus defined as

where is *l*(*β*
_{
hk
}, *γ*
_{
hk
}, *α*
_{0,hk
}(·)) the full log likelihood, {\alpha}_{0,\mathrm{hk}}^{\mathrm{\prime \prime}}\left(t\right) is the second derivate of the baseline intensity function, and *κ*
_{
hk
} is a positive smoothing parameter that controls the trade-off 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 cross-validation 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 *semi-Markovian* model (noted ({P}_{1}^{\zeta})) or a *non-homogeneous Markov model* (noted ({P}_{2}^{\zeta})) with a totally parametric approach or a semi-parametric approach. For example, we define {P}_{1}^{\widehat{\zeta}} and {P}_{2}^{\widehat{\zeta}} representing respectively the estimators of ({P}_{1}^{\zeta}) and ({P}_{2}^{\zeta}) based on a sample of i.i.d. observations ({\mathcal{O}}_{1},\dots ,{\mathcal{O}}_{n}) which come from the true unknown distribution *P*
^{∗}. To assess the risk of each estimator we propose to use the Expected Kullbak-Leibler risk (EKL) [6] defined as:

where {L}_{{\mathcal{O}}_{n+1}}^{{P}_{1}^{\widehat{\zeta}}/{P}^{\ast}}=\mathrm{log}{\frac{d{P}_{1}^{\widehat{\zeta}}}{d{P}^{\ast}}}_{|{\mathcal{O}}_{n+1}} is the log-likelihood ratio and {\mathcal{O}}_{n+1} is a new i.i.d observation. Commenges et al. [6] proposed the likelihood cross-validation (LCV) criterion to estimate this risk:

where {\mathcal{\mathcal{L}}}_{{\mathcal{O}}_{i}}^{{P}_{1}^{\widehat{\zeta}-i}} represents the likelihood of the observation {\mathcal{O}}_{i} based on the estimator {P}_{1}^{\widehat{\zeta}-i} defined without this observation.

Finally, to choose between the two estimators {P}_{1}^{\widehat{\zeta}} and {P}_{2}^{\widehat{\zeta}}, we compute some differences of Expected Kullbak-Leibler 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 Multi-state 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 multi-state 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 post-inclusion 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 {U}_{i}=({u}_{01}^{i},{u}_{02}^{i},{u}_{03}^{i}) and the covariates {Z}_{i}=({Z}_{01}^{i},{Z}_{02}^{i},{Z}_{03}^{i}) can be easily computed using the following expression with *k* = 1:

with,

The estimated posterior probabilities, {\widehat{\Pi}}_{0k}^{i}({t}^{\ast},{t}^{\ast}+h|{U}_{i},{Z}_{i})\phantom{\rule{2.77695pt}{0ex}}i=1,\dots ,N can be obtained by substituting the maximum penalized likelihood estimates of parameters *β*
_{0k
}, *γ*
_{0k
}, *α*
_{0,0k
}, the empirical Bayes estimates for {u}_{0k}^{i} and the individual information for the covariates {Z}_{0k}^{i} 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 time-fixed 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 semi-Markov model estimated by a semi-parametric approach (penalized likelihood estimator) with LCV = 4.044. Note that the difference of LCV between the semi-parametric estimator from the semi-Markov model and from the non-homogeneous Markov model (LCV = 4.055) is equal to 0.011 qualified as “moderate”. We have also compared a parametric approach versus a semi-parametric approach: for the non-homogenous model we obtained a difference of LCV equal to 0.407 qualified as “Large” between the parametric estimation (LCV = 4.462) and the semi-parametric estimation (LCV = 4.055); for the semi-Markov model we obtained a difference equal to 0.416 between the parametric estimation (LCV = 4.460) and the semi-parametric estimation (LCV = 4.044).

Table 2 presents the semi-parametric estimation results of the semi-Markov model. A complete description of these selected covariates in this model is provided elsewhere ([7]).

We also fitted a semi-Markov model without taking into account the intra-centre 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.07-1.79)), shock (Relative Risk= 1.28 (95%CI 1.01-1.62)) and the presence of an enteral nutritional therapy (Relative Risk= 1.24 (95%CI 1.00-1.53)) were incorrectly observed as significantly associated with the risk of VAP using the multi-state 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 between-centre 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.

We previously evaluated the heterogeneity between centres, but considering a different random effect for each transition of the multi-state model. We also fitted a joint frailty model for the two transitions 12 and 13, with a shared subject-specific 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 (\widehat{\eta}=0.77, one-sided 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 (two-sided Wald test = 0.16/0.29 = 0.55 < 1.96). This shows that times of deaths and discharge are correlated, but with a slightly non-signicant negative association (\widehat{\zeta}=-0.16).

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.

### 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 semi-Markov multi-state 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 {\beta}_{\mathrm{hk}}^{1}) and a continuous variable (with specific effect {\beta}_{\mathrm{hk}}^{2}) 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.

Briefly, a simulation of a semi-Markov 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 semi-markov 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 {\widehat{\beta}}_{\mathrm{hk}}=({\widehat{\beta}}_{\mathrm{hk}}^{1},{\widehat{\beta}}_{\mathrm{hk}}^{2}), and {\widehat{\gamma}}_{\mathrm{hk}}.

#### 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 multi-state model incorporating frailty terms, the fixed effect {\beta}_{\mathrm{hk}}=({\beta}_{\mathrm{hk}}^{1},{\beta}_{\mathrm{hk}}^{2}) 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 semi-Markov multi-state 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., multi-centre clinical trials, meta-analysis). Lack of software is a potential obstacle. We propose here a tractable model, semi-Markov as well as non-homogenous Markov, with semi-parametric or parametric estimates. The model can be readily derived with the R package FrailtyPack a simple and free approach, which does not require any time-consuming 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).

Vital status and time of death or time of discharge in ICU are known exactly. However, there may be more complex schemes with interval-censored times to events, i.e., the event occurs in a known time interval *L,R*. The semi-Markov or non-homogenous Markov models discussed in this paper do not allow the direct treatment of these interval-censored data. It would be interesting to extend the multistate frailty model in the case of interval-censoring [22]. This would lead to numerical integration in the estimation process due to the lack of a closed-form solution of the multiple integrals, and this can be very time-consuming 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 time-dependent 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 multi-state models [23] but without considering covariates. They also compared predictive accuracy of different multi-state models with or without frailties using a k-fold cross validated partial log-likelihood [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 mechanical-ventilation 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 multi-state model under consideration contains cluster-specific 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 {u}_{\mathrm{hk}}^{i} 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 S\left(t\right|{Z}_{\mathrm{hk}}^{\mathrm{ij}})={E}_{u}[S\left(t\right|{Z}_{\mathrm{hk}}^{\mathrm{ij}},{u}_{\mathrm{hk}}^{i}\left)\right]=
\frac{1}{{(1+{\gamma}_{\mathrm{hk}}{\Lambda}_{\mathrm{hk}}^{\mathrm{ij}}\left(t\right))}^{1/{\gamma}_{\mathrm{hk}}}} with {\Lambda}_{\mathrm{hk}}^{\mathrm{ij}}\left(t\right) 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 multi-centre study. We have suggested a method and software that gives accurate estimates and enables inference for any parameter or predictive quantity of interest.

## References

Ripatti S, Gatz M, Pedersen N, Palmgren J: Three-state frailty model for age at onset of dementia and death in Swedish twins. Genet Epidemiol. 2003, 24 (2): 139-149. 10.1002/gepi.10209.

Katsahian S, Boudreau C: Estimating and testing for center effects in competing risks. Stat Med. 2011, 30 (13): 1608-1617. 10.1002/sim.4132.

Katsahian S, Resche-Rigon M, Chevret S, Porcher R: Analysing multicentre competing risks data with a mixed proportional hazards model for the subdistribution. Stat Med. 2006, 25 (24): 4267-4278. 10.1002/sim.2684.

Fine J, Gray R: A proportional hazards model for the subdistribution of a competing risk. J Am Stat Assoc. 1999, 94: 496-509. 10.1080/01621459.1999.10474144.

Fine JP, Jiang H, Chappell R: On semi-competing risks data. Biometrika. 2001, 88: 907-920.

Commenges D, Joly P, Gegout-Petit A, Liquet B: Choice between semi-parametric estimators of Markov and non-Markov multi-state models from coarsened observations. Scand J Stat. 2007, 34: 33-52. 10.1111/j.1467-9469.2006.00536.x.

Nguile-Makao M, Zahar J, Français A, Tabah A, Garrouste-Orgeas M, Allaouchiche B, Goldgran-Toledano D, Azoulay E, Adrie C, Jamali S, et al: Attributable mortality of ventilator-associated pneumonia: respective impact of main characteristics at ICU admission and VAP onset using conditional logistic regression and multi-state models. Intensive Care Med. 2010, 36 (5): 781-789. 10.1007/s00134-010-1824-6.

Clech́ C, Timsit J, Lassence A, Azoulay E, Alberti C, Garrouste-Orgeas M, Mourvilier B, Troche G, Tafflet M, Tuil O, et al: Efficacy of adequate early antibiotic therapy in ventilator-associated pneumonia: influence of disease severity. Intensive Care Med. 2004, 30 (7): 1327-1333. 10.1007/s00134-004-2292-7.

Moine P, Timsit J, De Lassence A, Troche G, Fosse J, Alberti C, Cohen Y: Mortality associated with late-onset pneumonia in the intensive care unit: Results of a multi- center cohort study. Intensive Care Med. 2002, 28 (2): 154-163. 10.1007/s00134-001-1172-7.

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): 332-350. 10.1002/bimj.201000146.

Rondeau V, Mathoulin-Pelissier S, Jacqmin-Gadda 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): 708-721.

Hougaard P: Frailty models for survival data. Lifetime Data Anal. 1995, 1 (3): 255-273. 10.1007/BF00985760.

Pickles A, Crouchley R: A comparison of frailty models for multivariate survival data. Stat Med. 1995, 14 (13): 1447-1461. 10.1002/sim.4780141305.

Hougaard P: Multi-state Models: A Review. Lifetime Data Anal. 1999, 5: 239-264. 10.1023/A:1009672031531.

Rondeau V, Commenges D, Joly P: Maximum Penalized Likelihood Estimation in a Gamma-Frailty Model. Lifetime Data Anal. 2003, 9: 139-153. 10.1023/A:1022978802021.

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

Rondeau V, Mazroui Y, Gonzalez J: FRAILTYPACK: An R package for the analysis of correlated survival data with frailty models using penalized likelihood estimation. J Stat Software. 2012, 47 (4): 1-24.

Joly P, Commenges D, Letenneur L: A Penalized Likelihood Approach for Arbitrarily Censored and Truncated Data: Application to Age-Specific Incidence of Dementia. Biometrics. 1998, 54: 185-194. 10.2307/2534006.

Timsit JF, Zahar J, Chevret S: Attributable mortality of ventilator-associated pneumonia. Curr Opin Crit Care. 2011, 17 (5): 464-471. 10.1097/MCC.0b013e32834a5ae9.

Bekaert M, Timsit JF, Vansteelandt S, Depuydt P, Vesin A, Garrouste-Orgeas 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: 1133-1139. 10.1164/rccm.201105-0867OC.

Chastre J, Fagon J: Ventilator-associated pneumonia. Am J Respir Crit Care Med. 2002, 165 (7): 867-

Joly P, Commenges D, Helmer C, Letenneur L: A penalized likelihood approach for an illness–death model with interval-censored data: application to age-specific incidence of dementia. Biostatistics. 2002, 3 (3): 433-443. 10.1093/biostatistics/3.3.433.

Putter H, van Houwelingen H: Frailties in multi-state models: Are they identifiable? Do we need them?. Stat Methods Med Res. 2011, 23: 1-19.

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): 2080-2087. 10.1093/bioinformatics/btm305.

Melsen W, Rovers M, Koeman M, Bonten M: Estimating the attributable mortality of ventilator-associated pneumonia from randomized prevention studies. Crit Care Med. 2011, 39 (12): 2736-42.

Zahar J, Nguile-Makao M, Français A, Schwebel C, Garrouste-Orgeas M, Goldgran-Toledano D, Azoulay E, Thuong M, Jamali S, Cohen Y, de Lassence A, Timsit J: Predicting the risk of documented ventilator-associated pneumonia for benchmarking: construction and validation of a score. Crit Care Med. 2009, 9: 2545-51.

### Pre-publication history

The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2288/12/79/prepub

## 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.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### 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.

## Electronic supplementary material

### 12874_2012_840_MOESM1_ESM.r

Additional file 1: We provide a R code for simulating a data set and to analyze it: file name “supplementary-material-R-code.R”. (R 11 KB)

## Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

## Rights and permissions

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.

## About this article

### Cite this article

Liquet, B., Timsit, JF. & Rondeau, V. Investigating hospital heterogeneity with a multi-state frailty model: application to nosocomial pneumonia disease in intensive care units.
*BMC Med Res Methodol* **12**, 79 (2012). https://doi.org/10.1186/1471-2288-12-79

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/1471-2288-12-79