Integrating Latent Classes in the Bayesian Shared Parameter Joint Model of Longitudinal and Survival Outcomes
Eleni-Rosalina Andrinopoulou Kazem Nasserinejad Rhonda Szczesniak Dimitris Rizopoulos


	Integrating Latent Classes in the Bayesian Shared Parameter Joint Model of Longitudinal and Survival Outcomes
	
	

	
	Eleni-Rosalina Andrinopoulou^1, Kazem Nasserinejad^2, Rhonda Szczesniak^3,4 and Dimitris Rizopoulos^1
	
	
	

	
	
		1. Department of Biostatistics, Erasmus MC, Rotterdam, The Netherlands

		2. Department of Hematology, Erasmus MC, Rotterdam, The Netherlands

		3. Division of Biostatistics & Epidemiology and Division of Pulmonary Medicine, Cincinnati Children’s Hospital Medical Center, Cincinnati, United States

		4. Department of Pediatrics, University of Cincinnati, Cincinnati, United States
	
	

	
	Corresponding author: Eleni-Rosalina Andrinopoulou, Department of Biostatistics, Erasmus MC, PO Box 2040, 3000 CA Rotterdam, The Netherlands

		email: e.andrinopoulou@erasmusmc.nl, Tel: +31/10/7043731, Fax: +31/10/7043014
	
	
		Cystic fibrosis is a chronic lung disease which requires frequent patient monitoring to maintain lung function over time and minimize onset of acute respiratory events known as pulmonary exacerbations. It is important to characterize the association between key biomarkers such as FEV_1 and time-to first exacerbation. Progression of the disease is heterogeneous, yielding different sub-groups in the population exhibiting distinct longitudinal profiles. It is desirable to categorize these unobserved sub-groups according to their distinctive trajectories. Accounting for these latent classes, in other words heterogeneity, will lead to improved estimates of association arising from the joint longitudinal-survival model. 
		
		The joint model of longitudinal and survival data constitutes a popular framework to analyze longitudinal and survival outcomes simultaneously. Two paradigms within this framework are the shared parameter joint models and the joint latent class models. The former paradigm allows one to quantify the strength of the association between the longitudinal and survival outcomes but does not allow for latent sub-populations. The latter paradigm explicitly postulates the existence of sub-populations but does not directly quantify the strength of the association.  
		
		We propose to integrate latent classes in the shared parameter joint model in a fully Bayesian approach, which allows us to investigate the association between FEV_1 and time-to first exacerbation within each latent class. We, furthermore, focus on the selection of the optimal number of latent classes.  
	
	KEY WORDS: Cystic fibrosis, joint model, longitudinal outcome, survival outcome, latent class model
	
	

§ INTRODUCTION

	
	Cystic fibrosis (CF) is a lethal genetic disorder that primarily affects the lungs. The clinical course of CF is marked by progressive loss of lung function and typically results in respiratory failure. Forced expiratory volume in 1 second (hereafter, FEV_1) is the most important clinical indicator in monitoring lung function decline in patients with CF. Patients during follow-up might experience acute respiratory events referred to as pulmonary exacerbations. It is, therefore, of clinical interest to characterize the association between the longitudinal outcome FEV_1 and time-to first exacerbation. The motivation for our research comes from the US CF Foundation Patient Registry that consists of patients that were monitored from 2003 until 2015. In particular, we examined a subset of the Registry which consists of 1016 patients. These patients were six years and older and were observed with a median number of follow-up visits equal to six (with a range of 1-93 visits). The average age at baseline is 15 years (with a range of 6-21).  
	
	Several authors have studied the evolution of lung function over time, as summarized in a recent review <cit.>, however, to our knowledge little work has been done regarding the association of the lung function such as FEV_1 with time-to-event outcomes. In particular, joint modeling of longitudinal FEV_1 and survival outcomes in CF was introduced several years ago <cit.>, but has not been further used in CF epidemiology due to the computational burden of this approach. Furthermore, it is well recognized that different unobserved sub-groups of the biomarker FEV_1 exhibit different longitudinal profiles <cit.>. Patients can be categorized in several sub-groups (latent classes) with different trajectories. It is, therefore, of high clinical interest to measure the strength of association between FEV_1 with the risk of first exacerbation accounting for the latent trajectories. 
	
	The joint model of longitudinal and survival data constitutes a popular framework to analyze longitudinal and survival outcomes jointly <cit.>. In particular, two paradigms within this framework are the shared parameter joint models and the joint latent class models. The former paradigm links the longitudinal and the survival process via the random effects <cit.>, which does not allow for latent classes. The latter paradigm <cit.>, which associates the two processes through latent classes, explicitly postulates the existence of sub-populations but does not directly quantify the strength of the association. 
	
	The aim of the paper is twofold. Firstly, to model the relationship between FEV_1 and time-to first exacerbation. For this purpose, we propose a Bayesian shared parameter joint model that integrates latent classes inherent in this heterogeneous population. This model will assess the strength of association between the two outcomes while allowing for latent classes. Secondly, to address a problem that arises in latent class models, which is the selection of the optimal number of classes. Several approaches have been proposed in the literature both in frequentist and Bayesian frameworks, including among others the use of information criterion, Bayes factors and reversible jump MCMC. These approaches are computationally intensive and can require the fit of several models with different numbers of classes, which can be time-consuming. To overcome this problem, we will implement the method of <cit.> to our joint model. This method is a pragmatic extension of <cit.> criterion that showed that when we overfit a mixture model by assuming more latent classes than present in the data, the superfluous latent classes will asymptotically become empty if the Dirichlet prior on the class proportions is sufficiently uninformative. <cit.> performed an extensive simulation study to further investigate this approach and used it as a criterion also in longitudinal studies for obtaining the optimal number of classes by simply excluding latent classes that are negligible in proportion.  
	
	
	

§ JOINT MODEL ESTIMATION

	

 §.§ Longitudinal submodel

	
	To account for the fact that the population is heterogeneous and consists of G possible unobserved sub-groups, we postulate a latent class mixed-effects model <cit.>. We let y_i denote the longitudinal response vector for the ith patient (i = 1, …, n) obtained at different time points t_ij>0, (j = 1, …, n_i). In particular, we have
	
    y_i(t | v_i = g) = η_ig(t) + ϵ_i(t) = x^⊤_i(t)β_g + z^⊤_i(t)b_ig + ϵ_i(t),

	where v_i = g  (g = 1, …, G) presents the latent class indicator, x_i(t) denotes the design vector for the fixed effects regression coefficients β_g and z_i(t) the design vector for the random effects b_ig. Moreover, ϵ_i(t) ∼ N(0, σ_y^2). For the corresponding random effects, we assume a multivariate normal distribution, namely
	
    b_ig∼ N(0,Σ_bg),

	where N denotes the normal distribution and Σ_bg is the variance-covariance matrix of the random effects. An individual has a probability π_ig = P(v_i = g) of belonging to latent class g. Using a multinomial distribution we obtain the class of each individual as,
	
    v_i ∼ Multinomial(π_ig).

	According to the specification of the latent class mixed-effects submodel (<ref>), both fixed and random effects are class-specific, whereas the measurement error ϵ_i(t) is not.
	
	
	

 §.§ Survival submodel

	We let T_i^* denote the true failure time for the i-th individual, and C_i the censoring time. Moreover, T_i=min(T_i^*,C_i) denotes the observed failure time and δ_i = {0, 1} is the event indicator where zero corresponds to censoring.
	We postulate a joint model for the relationship between the survival and the longitudinal outcome. Specifically, we have
	
    h_i(t | v_i = g)=h_0g(t)exp [γ_g^⊤w_i+ α_g η_ig(t)],

	where w_i is a vector of baseline covariates with a corresponding vector of regression coefficients γ_g and h_0g(t) is the baseline hazard. Specifically, the B-splines baseline hazard function is assumed log h_0g(t)=γ_h_0g,0+∑_q=1^Qγ_h_0g,qB_q(t,ν),
	where B_q(t,ν) denotes the q-th basis function of a B-spline with knots ν_1,…, ν_Q and γ_h_0g the vector of spline coefficients. The knots are placed at equally spaced percentiles of the observed event times. Furthermore, α_g denotes the association parameter for the gth class. According to the specification of the survival submodel (<ref>) the baseline covariates, the baseline hazard and the association parameter are class-specific parameters. The proposed model goes beyond the standard joint model and joint latent class model where a single or no association parameter is assumed and provides a class-specific association. This is a more realistic assumption for the motivating data set since it is clinically expected that the risk of the first exacerbation will be higher when the rate of FEV_1 decline is faster. Accounting for these latent classes will lead to improved estimates of association arising from the joint model. 
	
	
	

§ BAYESIAN ESTIMATION

	
	We employ a Bayesian approach where inference is based on the posterior distribution of parameters in the model. We use Markov chain Monte Carlo (MCMC) methods to estimate the parameters of the proposed model. The likelihood of the model is derived under the assumption that the longitudinal and survival processes are independent given the random effects. Moreover, the longitudinal responses of each subject are assumed independent given the random effects <cit.>. The likelihood contribution for the i-th patient is written as
	
    [                                                             p(y_i,T_i,δ_i | v_i = g, θ, b_ig) =;                                                                                                 ∑_g=1^Gπ_ig{∏_j=1^n_i[ p(y_ij| v_i = g, θ_y, b_ig) ] p{T_i,δ_i| v_i = g, η_ig(T_i),θ_s, b_ig}}, ]

	
	where θ = (θ_s^⊤, θ_y^⊤, π_ig)^⊤ with θ_y = (β_g, σ_y, Σ_bg) and θ_s = (γ_g, α_g, γ_h_0g). 
	
	The likelihood contribution of the longitudinal outcome takes the form
	
    p(y_ij| v_i = g, θ_y, b_ig)=(2πσ_y)^-1/2exp[  -(y_ij -x_ij^⊤β_g - z_ij^⊤b_ig)^2/2σ_y^2] .
 
	
	The likelihood contribution of the survival model is given by
	
    p{T_i,δ_i| v_i = g, η_ig(T_i),θ_s , b_ig}= 
          exp[γ_h_0g,0 + ∑_q=1^Qγ_h_0g,qB_q(T_i,ν)+ γ_g^⊤w_i+ η_ig(T_i)  α_g] ^I(δ_i=1)×
          exp{-exp(γ_g^⊤w_i)∫_0^T_iexp[γ_h_0g,0 + ∑_q=1^Q γ_h_0g,qB_q(s,ν)+ η_ig(s) α_g,  ]ds  }.
 
	
	The posterior distribution is written as
	
    [ p(θ, b_g|y,T,δ) = 
    	∏_i=1^n  p(y_i,T_i,δ_i | v_i = g, θ, b_ig)p(b_ig| v_i = g, θ_y)p(θ). ]

	where
	
    p(b_ig| v_i = g, θ_y) = [2π(Σ_bg)]^-1/2exp(- b_ig^⊤Σ_bg^-1b_ig/2) ,

	and p(θ) denotes the prior distributions.
	
	A commonly used prior in mixture models for the class probability is a Dirichlet distribution. In particular,
	
    π_ig = P(v_i = g) ∼ Dirichlet(a).

	Small values of a = {a_1… a_G} correspond to a less informative prior and a flat prior distribution is obtained when each a_g is equal to 1. The selection of a is an important task and will be discussed in Section <ref>. Standard priors can be assumed for the rest of the parameters. In particular, for the coefficients of the longitudinal fixed effects, the survival covariates and the baseline hazard, normal priors can be taken. For the variance-covariance matrix of the random effects we can assume an inverse Wishart prior, while for the precision parameter of the longitudinal outcome we can assume an gamma prior. 
	
	
	
	

 §.§ Selection of Number of Classes

	An important task in latent class models is to identify the optimal number of classes. Several approaches have been previously proposed for choosing the optimal number of classes in both frequentist and Bayesian settings. Common examples are the Bayesian information criterion (BIC) <cit.>, deviance information criterion (DIC) <cit.> and other Bayesian approaches such as Bayes factor and reversible jump MCMC algorithm <cit.>. A drawback of the aforementioned approaches is that they are computationally intensive and some require the fit of models assuming different numbers of classes, which might be time-consuming for complex models such as the joint models of longitudinal and survival outcomes. 
	
	An interesting alternative was proposed by <cit.>, where they proved that in overfitted mixture models (with more latent classes than present in the data), the superfluous latent classes will asymptomatically become empty if the Dirichlet prior on the class proportion is sufficiently uninformative. Recently, <cit.> used this approach and proposed a latent class selection procedure for longitudinal models. An overfitted mixture model converged to the true mixture by assigning a small portion of individuals to empty classes, if the parameters of the Dirichlet prior a are smaller than d/2, where d is the number of class-specific parameters. Furthermore, uninformative priors for the rest of the parameters are required. The steps are described as follows: 
	
	
		
  * First, a latent class model with a large enough number of latent classes is fitted. 
		
  * Then, the number of non-empty classes at each iteration is calculated as:
		
    g_k,opt = G - ∑_g=1^G I(n_k,g/n≤ψ),

		where G is the total number of classes, k represents the iteration, n_k,g is the number of patients in class g at iteration k, n is the total number of patients and ψ is a predifined value. 
		
  * After obtaining the non-empty classes per iteration, the posterior mode of the non-empty classes is calculated. 
		
  * Finally, the model with the optimal number of classes which are the non-empty classes is refitted. 
	
	
	Advantages of this approach are that it is easy to implement even in such complex models and it is not influenced by the label switching problem since we observe the non-empty classes at each iteration. The only time that we need to correct for label switching is when we fit the final model with the optimal number of classes. Furthermore, this approach requires us to fit the model only two times, (namely one with the high number of classes and one with the optimal number of classes) instead of assuming all possible number of classes, therefore decreasing computational burden. It has been shown through extensive simulations in the longitudinal setting that this method performs better than alternative model selection criteria such as BIC and DIC <cit.>.
	
	
	

§ ANALYSIS OF THE CF DATA

	
		
	
	
	

	In this section we present the analysis of the motivating data set introduced in Section <ref>. Our primary focus is to investigate the association between FEV_1 and time-to first exacerbation by taking into account that we have sub-groups with different evolution over time for FEV_1. The first step is to obtain the optimal number of classes that can explain the heterogeneity of the population. From the literature, it is known that two or three classes are observed for the evolution of FEV_1 outcome <cit.>. Therefore, for the selection process, we fitted a joint model assuming six classes. For the longitudinal outcome, we assumed a linear mixed-effects submodel including natural cubic splines for time (modeled as age, in years) with two internal knots at 13.76 and 17.62 years (corresponding to 33.3% and 66.67% of the observed follow-up times) in both the fixed and random effects parts. The DIC criterion and subject-specific plots (with observed and predicted values) were used to investigate the need of non-linear evolution over time in a mixed-effects linear model. Furthermore, we corrected for some baseline characteristics. These variables, together with descriptive statistics, are presented in Table <ref>. In Figure <ref> the FEV_1 evolutions of 25 randomly selected patients with more than two repeated measurements are presented.

	
	
	
	
Specifically, the model takes the form,
	
    [                                                                               y_i(t)                            = η_ig(t) + ϵ_i(t) = β_0g + ∑_ω =1^3
    	β_ω gns(Age_i, ω) +;                                                                                               β_4gGender +  ∑_ω =5^7 β_ω gF508 + β_8gHispanic + β_9gWhite + β_10gSESlow +;                                                                                              β_11gMRSA + β_12gMSSA + β_13gPa + β_14gaspergillus +  ∑_ω =15^16β_ω  gCFRD +;                                                                                           β_17gPancEnzymes + β_18gnumVisityr + ∑_ω =19^21β_ω  gns(Age_i, ω - 18):Gender +;                                                                                      ∑_ω  =22^30β_ω  gns(Age_i,ω - 21):F508 +  ∑_ω =31^33β_ω  gns(Age_i, ω - 30):SESlow +;                                                                                                                                     ∑_ω =1^3 b_ω  gns(Age_i, ω) + ϵ_i(t). ]

	
	
	To investigate the association between FEV_1 and time-to first exacerbation, we postulated the proposed joint latent class model:
	
    h_i(t,θ_s)=h_0g(t)exp [γ_g Gender_i + α_g η_ig(t) ].

	For the baseline hazard we assumed a quadratic B-splines basis with 8 equi-distance internal knots ranging from zero until 19.25 years. 
	
	In the Dirichlet distribution for the prior of the class probability, following the recommendation in <cit.>, we assumed a smaller than d/2 (where d is the number of class-specific parameters). To ensure that we have the same scale for the coefficients of the covariates in order to easier select uninformative priors, we standardized the FEV_1 outcome and the continuous variables (age and numVisityr). Relatively uninformative priors were selected for the parameters in the model. These priors are as follows: 
	
		
  * β_g∼ N(0, 1000),
		
  * γ_g ∼ N(0, 1000),
		
  * γ_h_0g,q∼ N(0, 1000),
		
  * α_g ∼ N(0, 100),
		
  * σ_y^2 ∼ GA^-1(0.01,0.01)
		
  * Σ_bg∼ W^-1(M, df),
	
	where GA^-1 denotes the inverse gamma distribution and W^-1 denotes the inverse Wishart distribution with M=diag(0.01) being the scale matrix and df the degrees of freedom which is set as the total number of the random effects. For the variance of the association parameter no large variance was required to ensure that we have a uninformative prior since, with the standard joint model we obtained an association parameter smaller than 0.1. The selection of the variances of these priors was investigated with simulations. We ran the MCMC using a single chain with 300,000 iterations, 250,000 burn-in and 10 thinning. The results indicate the presence of three or four classes, assuming that a class is empty if it contains 10 to 15% of the patients (10%≤ψ≤15%). Since it is established in the literature that two or three classes are present in such populations, we decided to continue with three classes <cit.>.
		
	We reran the model assuming three classes and the normal scale of the continuous covariate age and FEV_1 outcome (we standardized only the numVisityr variable). We ran the MCMCs with a single chain for 500,000 iterations, with a burn-in of 450,000 and thinning of 10 and we fixed the label switching problem. Convergence was monitored by trace plots. Table A1 in the Appendix shows the mean and standard deviation of age (at baseline), FEV_1 (at baseline) and number of visits (at last follow-up) per class, while Table A2 shows the percentage of the categorical variables (at baseline) pes class. In Figure <ref> we illustrate the evolution of the longitudinal outcome in each class assuming patients who are F508del homozygotes, non-Hispanic, White, without low SES, are not infections with MRSA, MSSA or aspergillus, do not use pancreatic enzyme, do not have pseudomonas aeruginosa, have normal CFRD and had five visits within the prior year (which is the mean value of all observations). In particular, the upper plots represent female patients while the lower plots represent male patients. We obtain a faster progression in class one for both females and males. Patients in class two have a stable evolution in the middle of the follow-up and patients in class three are stable throughout the follow-up period. In addition, patients in class one and two start from a higher FEV_1 compared to patients in class three. In Figure <ref> we illustrate the evolution of the longitudinal outcome in each class assuming patients who are F508del homozygotes, Hispanic, White, have low SES, are infections with MRSA, MSSA, aspergillus, use pancreatic enzymes, have pseudomonas aeruginosa, have impaired CFRD and had eight visits within the prior year. Again, the upper plots represent female patients while the lower plots represent male patients. We obtain that patients in these classes start from a lower FEV_1 value compared Figure <ref>. In addition, we observe a faster progression in class one for both female and male patients. The mean and the credible interval of the MCMC samples of the association parameters per class are presented in Figure <ref>. We obtain a weak association between FEV_1 and time-to first exacerbation for the second and third class, while a stronger negative association for class one. 
	
	
	
	
	
	
	
	
	
	


§ SIMULATIONS

	
	We performed a series of simulations to investigate the proposed class selection method on the joint modeling framework. 
	
	

 §.§ Design

	
	We assumed N_1 = 350, N_2 = 525 and N_3 = 1050 patients with maximum number of repeated measurements equal to ten. To simulate the continuous longitudinal outcome, we used the following linear mixed-effects model per data set. In particular,
	
    y_i(t) = η_i(t) + ϵ_i(t) = β_0 +  β_1male_i + β_2 t + b_0i + b_1it+ ϵ_i(t),

	where ϵ_i ∼ N(0,σ_y^2) and b_i=(b_0i, b_1i) ∼ N_2(0,Σ_b). For simplicity, we adopted a linear effect of time for both the fixed and the random part, and corrected for a binary variable (male_i). Time t was simulated from a uniform distribution between zero and 19.5.
	For the survival part, we assumed the following model:
	
    [                            h_i(t) =h_0(t)exp{γ_^⊤Age_i + αη_i(t) }. ]

	The baseline risk was simulated from a Weibull distribution h_0g(t)=ξ t^ξ -1. For the simulation of the censoring times, an exponential censoring distribution was chosen so that the censoring rate was between 40% and 60%. Age was simulated from a normal distribution with mean 45 and standard deviation 15.7.
	Under this setting we simulated three different data sets that have different parameters for the fixed effects in the longitudinal submodel, the baseline covariates and baseline hazard in the survival submodel, the variance-covariance matrix of the random effects and the association parameter (more details are presented in Table <ref>). Figure <ref> illustrates the evolution of the longitudinal outcome per group from the simulation parameters for each one of the three data sets.
	
	
	
	
	
	
	
	
	
	

 §.§ Analysis

	
	In order to investigate the proposed class selection approach, we applied the model in three different Scenarios. For Scenario I we combined all three data sets assuming N_1 = 350 individuals in each of them, for Scenario II we combined the first two data sets with N_2 = 525 and finally, for Scenario III we used only the first data set with N_3 = 1050. We fitted the proposed joint model assuming six classes where all the parameters were class-specific except for the measurement error in the mixed-effects model. These include the fixed effects from the longitudinal submodel (3 parameters), the baseline covariates from the survival model (2 parameters), the baseline hazard (5 parameters) from the survival submodel, the variance-covariance matrix of the random effects (3 parameters) and the association parameter (1 parameter). To simplify the simulations, for the baseline hazard we assumed quadratic B-splines basis with 4 equally distance internal knots. We assumed a_g = 6.9 which is smaller than the total number of the class-specific parameters divided by two. The same priors were used as in the application and we ran the MCMC using a single chain with 50,000 iterations, 25,000 burn-in and 10 thinning. 
	
	We performed 150 simulations per Scenario. We compared our proposed method with the joint latent class model using function  from the  package in  developed by <cit.>, where we used the BIC as criterion. In particular, we assumed that the covariates from the fixed effects in the mixed-effects model, the variance-covariance of the random effects, the baseline covariate and baseline hazard in the survival model are class-specific. We, furthermore, assumed a cubic M-splines baseline risk function. 
	
	

 §.§ Results

	
	The results from the different Scenarios assuming different cut off percentage ψ indicating when a class is defined as empty are illustrated in Table <ref>. In particular, we present the percentage of true number of classes and the mode of the number of classes.
	
	For Scenario I, we obtain the highest percentage when assuming ψ to be between 12-15%. In particular, we obtain around 55% of the time the correct number of classes and a mode equal to the correct number of classes (three). On the other hand, the BIC in the joint latent class model selects only 20% of the time the correct number of classes. Furthermore, this method seems to underestimate the true number of classes (mode equal to one). 
	
	For Scenario II, we obtain the highest percentage when ψ is between 8-15%. In particular, we obtain around 50% of the time the correct number of classes and a mode equal to the correct number of classes (two). On the other hand, the BIC in the joint latent class model selects only 12% of the time the correct number of classes. Similar to Scenario I, this method seems to underestimate the true number of classes (mode equal to one). 
	
	Finally, for Scenario III, the BIC in the joint latent class model seems to perform better than the proposed approach where it always selects the correct number of classes (one). This is not surprising, since the BIC always underestimated the true number of classes in the previous Scenarios. Using the proposed approach and assuming that the ψ is equal to 15%, we obtain 43% of the time the correct number of class and a mode equal to two. In this Scenario some convergence problems were detected. When recalculating the % and mode including only the simulations that were converged we obtain similar percentages for the correct number of class and a mode equal to 1 when ψ is 15%.
	
	
	
	
	
	
	

§ DISCUSSION

	In this paper we proposed a shared parameter joint model incorporating latent classes. Applying it to CF data, this model accounted for patient heterogeneity inherent in the progression of FEV_1. Compared to previously proposed joint latent class models <cit.> we obtained the strength of the association between FEV_1 and time-to first exacerbation per group of patients. Finally, we focused on the selection of the optimal number of classes and used an overfitted mixture model (high number of classes) to obtain the non-empty classes.
	
	A limitation of this approach is that it requires an intensive computational effort. In particular, for the class selection, where a model with a high number of classes is required, the number of parameters increases drastically. This, in combination with the high number of observations in the CF application increases the computational time that is required. Considering the difficulty of this model, it is almost impossible to obtain the optimal number of classes with other Bayesian criterion. Implementing the proposed criterion is straightforward; however, due to the complexity of the model it is computationally expensive to fit a model with a larger number of classes, e.g., 10. This could also explain the fact the a higher percentage for the predefined number ψ was required in order to obtain the non-empty number of classes. It was shown in the simulation analysis that the BIC always underestimated the true number of parameters and it, therefore, performed better when the true number of classes was one. Even though, in that Scenario the proposed method did not work perfectly, it seems to be better than other criteria and easier to perform.
	
	Although there is a large database available in the US Registry, we used only a subset in order to make it feasible to run the proposed model. This subset has particular characteristics and it cannot be generalized to all patients in the Registry. Therefore, the presented results do not reflect the diversity of the whole database.
	
	Possible extensions would be to include more covariates also in the survival submodel in order to take into account extra information regarding the patients. Furthermore, using the proposed model for obtaining future FEV_1 measurement and time-to first exacerbation probabilities, could lead to more efficient treatment prioritization and clinical management for patients with CF.
	
	
	
	
	abbrvnat 
	
	
	
	

§ APPENDIX




		
	
	
	
	
	
	
	
	
