Flexible semiparametric joint modeling: an application to estimate individual lung function decline and risk of pulmonary exacerbations in cystic fibrosis

Background Epidemiologic surveillance of lung function is key to clinical care of individuals with cystic fibrosis, but lung function decline is nonlinear and often impacted by acute respiratory events known as pulmonary exacerbations. Statistical models are needed to simultaneously estimate lung function decline while providing risk estimates for the onset of pulmonary exacerbations, in order to identify relevant predictors of declining lung function and understand how these associations could be used to predict the onset of pulmonary exacerbations. Methods Using longitudinal lung function (FEV1) measurements and time-to-event data on pulmonary exacerbations from individuals in the United States Cystic Fibrosis Registry, we implemented a flexible semiparametric joint model consisting of a mixed-effects submodel with regression splines to fit repeated FEV1 measurements and a time-to-event submodel for possibly censored data on pulmonary exacerbations. We contrasted this approach with methods currently used in epidemiological studies and highlight clinical implications. Results The semiparametric joint model had the best fit of all models examined based on deviance information criterion. Higher starting FEV1 implied more rapid lung function decline in both separate and joint models; however, individualized risk estimates for pulmonary exacerbation differed depending upon model type. Based on shared parameter estimates from the joint model, which accounts for the nonlinear FEV1 trajectory, patients with more positive rates of change were less likely to experience a pulmonary exacerbation (HR per one standard deviation increase in FEV1 rate of change = 0.566, 95% CI 0.516–0.619), and having higher absolute FEV1 also corresponded to lower risk of having a pulmonary exacerbation (HR per one standard deviation increase in FEV1 = 0.856, 95% CI 0.781–0.937). At the population level, both submodels indicated significant effects of birth cohort, socioeconomic status and respiratory infections on FEV1 decline, as well as significant effects of gender, socioeconomic status and birth cohort on pulmonary exacerbation risk. Conclusions Through a flexible joint-modeling approach, we provide a means to simultaneously estimate lung function trajectories and the risk of pulmonary exacerbations for individual patients; we demonstrate how this approach offers additional insights into the clinical course of cystic fibrosis that were not possible using conventional approaches. Electronic supplementary material The online version of this article (10.1186/s12982-017-0067-1) contains supplementary material, which is available to authorized users.


Background
Maintaining pulmonary function is essential for survival in individuals with cystic fibrosis (CF), a lung disease that currently affects nearly 70,000 individuals worldwide [1]. Lung function is measured primarily as forced expiratory volume in 1 s of percent predicted (hereafter, FEV 1 ), which is the primary marker of disease severity in individuals with CF. FEV 1 declines in a nonlinear fashion over age and exhibits substantial heterogeneity both between patients and within an individual patient over time (Fig. 1). Most individuals with CF have decreased FEV 1 over age, but initial FEV 1 and rate of decline vary between patients. Linear mixed modeling is an established approach to estimate age-related FEV 1 progression. Historically, CF epidemiologic studies have not accounted for nonlinearity in the FEV 1 trajectory; however, recent approaches have included piecewise polynomials, either in the form of regression splines or change-point models, to estimate decline in FEV 1 [2,3]. These approaches have shown that FEV 1 decline is variable with maximal loss occurring in adolescence and early adulthood.
Meanwhile, the clinical course of CF is often marked by the occurrence of an acute respiratory event known as a pulmonary exacerbation, which intensifies the severity of this chronic disease. Increased pulmonary symptoms and decreased lung function and weight are common clinical indicators of this event [4]. A study of the United States Cystic Fibrosis Foundation Patient Registry showed that, after having a pulmonary exacerbation, 25% of patients' lung function levels did not recover to baseline (prepulmonary exacerbation) levels [5]; this study identified lower socioeconomic status (use of Medicaid insurance), malnutrition, having respiratory infections and pathogens, and being female as risk factors. Another study identified pulmonary exacerbations as a risk factor for more rapid lung function decline [6].
Although much information has been gleaned, these and other epidemiologic studies in CF have modeled the longitudinal and time-to-event processes of lung function and pulmonary exacerbations separately. Efficiency and additional information about the disease course may be gained by modeling them simultaneously. Joint modeling of longitudinal and event processes have enjoyed a renaissance in recent years among biostatisticians and epidemiologists, largely due to software developments and modern applications demonstrating their clinical utility for the purposes of monitoring disease progression and making predictions [7]. Developments in joint models have also been motivated by data on CF disease progression, focusing on simultaneously fitting a longitudinal submodel of lung function and a survival submodel. For example, Schluchter et al. [8] developed a joint longitudinal-survival model to create prognostic indicators of CF disease severity, such as predicted age at death and other empirical Bayes estimates of parameters in the longitudinal submodel, such as slope, which could be used to estimate rate of decline in FEV 1 adjusted for survivor bias. They combined a linear mixed effects model for longitudinal FEV 1 with a Gaussian model for age at death and applied it to data taken from a local CF center. They later extended this model to account for left truncation, as follow-up does not necessarily begin at age zero for CF patients [9]. Such joint models can be classified as shared parameter approaches, meaning that FEV 1 and the event of interest are assumed to be conditionally independent given a set of random effects.
In this study, we utilize the aforementioned advantages of joint modeling in the shared-parameter framework and combine them with flexible semiparametric regression splines to estimate FEV 1 trajectories in the longitudinal submodel. In this context, joint modeling would account for informative dropout due to pulmonary exacerbation events; not accounting for these events could bias estimates of the FEV 1 trajectories. Our objectives are to identify relevant predictors of lung function trajectories and pulmonary exacerbation events, and to understand how associations found from a joint model could be used to predict the onset of pulmonary exacerbations for individual CF patients.

Study design and cohort
This longitudinal cohort study consisted of clinical encounter and hospitalization data from the U.S. Cystic Fibrosis Foundation Patient Registry (CFFPR) between January 1, 2003, and December 31, 2011. This registry has been tracking CF demographic and clinical variables for more than 40 years, and has been thoroughly described elsewhere [10]. Because the majority of relevant predictors of lung function decline and pulmonary exacerbation events were consistently documented beginning in 2003, we considered data available from this year onward. Patients younger than 6 years of age were excluded due to the potential for unreliable data from pulmonary function testing. Data from those older than 45 years of age were excluded due to the possibility that these individuals have milder phenotypes that are not representative of typical CF disease progression.
The longitudinal outcome of interest, FEV 1 % predicted (hereafter, FEV 1 ), was obtained at each clinical encounter; % predicted values were calculated using Wang and Hankinson reference equations [11,12]. The event outcome of interest, onset of first pulmonary exacerbation, was subject to censoring. A pulmonary exacerbation was considered to have occurred if documented in the CFFPR as warranting treatment with intravenous antibiotics. Analyses were restricted to patients with at least ten measurements of FEV 1 . Longitudinal data were included up to the time of the first pulmonary exacerbation on record and censored thereafter. Baseline was defined as the time point at which the first FEV 1 measurement was available. Covariates of interest were selected from the literature on FEV 1 [13,14] and pulmonary exacerbation [5,6] as separate models, and included age, gender, initial FEV 1 measure and baseline lower socioeconomic status (lower SES, defined as having only state/federal or no insurance; recorded as 1, and 0 otherwise); a birth cohort covariate was used to adjust for potential left truncation bias; time-varying covariates included CF-related diabetes (CFRD, with or without fasting hyperglycemia), positive cultures for methicillin-resistant staphylococcus aureus (MRSA), Burkholderia cepacia (B. cepacia), and Pseudomonas aeruginosa (Pa).

Joint model formulation and estimation
Each joint model consists of two linked submodels, a mixed model to fit longitudinal FEV 1 and a Weibull model to fit pulmonary exacerbation data.

Longitudinal submodel for FEV 1
Both the separate modeling approach and the joint modeling approach require specification of a model for the longitudinal measure of FEV 1 . In this subsection we describe this model before linking it to the time-to-event model in subsequent sections. Suppose there are N patients in the CFFPR indexed by i = 1, 2, . . . , N. Let FEV 1 ij represent lung function measured for the ith patient at the jth time, represented as age ij (corresponding to patient age, in years), j = 1, . . . , n i . The longitudinal submodel consists of a linear mixed model with random effects: where the population-level mean response f age ij is modeled as a continuous function of time, which combines the time-related fixed effects terms and regression splines 15 ; m i age ij = x ′ 1i age ij α 1 refers to the subject-level fixed effects. The vector x 1i age ij represents static and time-varying covariates defined above, and the vector α 1 contains their corresponding regression coefficients. The expression W 1i age ij = z ′ 1i age ij U i incorporates random effects that describe how subjectspecific true FEV 1 levels deviate from their expected behavior, where the vector U i corresponds to subjectspecific random slopes and intercepts. We used the usual Laird-Ware form 8 , W 1i age ij = U 0i + U 1i age ij , corresponding to z 1i age ij = 1, age ij ′ ; this specification allows individuals to have varying baseline FEV 1 measurements and different rates of change in FEV 1 over time. Lastly, the term ε ij ∼ N 0, σ 2 ε denotes zero-mean Gaussian measurement error.
We considered two nonlinear representations of f age ij , which we hereafter refer to as semiparametric and cubic, and examined both in the joint modeling framework. First, we used penalized splines with a cubic truncated power basis in order to provide smooth estimates of the longitudinal course of FEV 1 measurements. Second, we considered only global cubic polynomials to represent population-level FEV 1 decline by excluding the basis functions in f age ij . Details of the two formulations are available in the "Appendix 1". Taking the first derivative with respect to age in Eq. (1) yields overall and subject-specific estimates of the rate of change in FEV 1 , and has been used previously to model FEV 1 decline [2].

Time-to-event submodel for pulmonary exacerbation
Let T i denote the possibly censored survival time to a pulmonary exacerbation event for the ith patient. In a Weibull model, we assume that the survival time follows a Weibull distribution, that is The hazard at time t for the ith patient is which monotonically increases with time if k > 1, decreases if k < 1, and reduces to the exponential hazard and remains constant if k = 1. The vectors x 2i (t) and α 2 represent possibly time-dependent covariates and the corresponding regression coefficients. Covariates x 2i (t) need not have elements in common with x 1i age in the longitudinal model as shown in Eq. (1). Notice that W 2i could have a time-dependent form [15], but we do not consider it here since this level of complexity is not required for our data. Similar to the form of W 1i age in Eq. (1), W 2i = θ 0 U 0i + θ 1 U 1i corresponds to patientspecific random intercepts and slopes, but has distinct regression parameter coefficients, θ 0 and θ 1 . If none of the covariates vary over time, W 2i reduces to 0 in the absence of random effects. The main idea of this approach is to connect the longitudinal and survival processes with a latent bivariate Gaussian process. We can add a frailty term in

Submodel links
Pulmonary exacerbation event times were associated with longitudinal FEV 1 measurements through stochastic dependence between W 1i and W 2i from Eqs. (1)-(2) by assuming: The subject-specific random intercept and slope, depicting varying initial values and rates of FEV 1 decline for each patient after accounting for the nonlinear FEV 1 trajectory, are contained in the vector U i , which is shared between the longitudinal and time-to-event submodels.
Here U i is assumed to follow a bivariate normal distribution N (0, ), where has diagonal entries σ 2 0 , σ 2 1 and off-diagonal entries σ 01 . The joint model with (4) allows both the random intercept U 0i and slope U 1i , involved in (3), to affect the risk of the event. Thus, deviations of the patient-specific FEV 1 trajectories from the population-level FEV 1 curve enter the pulmonary exacerbation model in the form of random intercepts and slopes. If this type of association exists between the longitudinal FEV 1 and pulmonary exacerbation event processes, then inference from the joint model should be less biased and more efficient, compared to modeling the processes separately [16]. Given that absolute FEV 1 (intercept) and rate of change in FEV 1 (slope) are on different scales, we used the estimated SDs of the random effects to express each corresponding HR as per SD change in the respective covariate. These quantities were obtained from the joint model by estimating exp {θ 0 * σ 0 } and exp {θ 1 * σ 1 }, where θ 0 and θ 1 are the association parameters estimated ordinarily from the HRs exp {θ 0 } and exp {θ 1 }, and σ 0 and σ 1 are the respective SDs of intercept U 0i and slope U 1i as defined previously.

MCMC sampling procedure
Conventional likelihood-based estimation of the parameters requires integration of the two submodels over the distribution of random effects. As the number of random effects in the model increase, the dimension of integration in the joint model likelihood increases; as a result, parameter estimation typically involves specialized numerical algorithms and becomes computationally burdensome [17]. Alternatively, we employed Markov-Chain Monte-Carlo (MCMC) via Gibbs sampling in WinBUGS [18] to simulate data from the respective posterior distributions under each model. The highest posterior density (HPD) and accompanying 95% credible interval (CI) were used to estimate each parameter of interest. Specification of priors and sampling procedures are detailed in the "Appendix 2".

Model comparisons
We calculate the deviance information criterion (DIC) to compare the performance of the separate and joint models [19]. The DIC measure balances the fit of a model to the data with its complexity. The components of DIC for the two submodels, denoted as DIC 1 and DIC 2 , were also provided to evaluate their relative contributions to the total DIC score. Formulas are provided in the "Appendix 3". A smaller value of DIC indicates the preferred model.

Research ethics
The Internal Review Board of Cincinnati Children's Hospital Medical Center (Cincinnati, OH, USA) approved the study.

Longitudinal submodel for FEV 1
Both the separate and joint modeling of longitudinal FEV 1 indicated that decline was nonlinear over age ( Table 1, Semiparametric longitudinal submodel). In both models, having higher FEV 1 at entry corresponded to higher FEV 1 over time; lower SES, diagnosis with CFRD and the presence of infections were associated with lower FEV 1 over age. The younger birth cohorts appeared to have similar FEV 1 progression, while being Table 1 Posterior estimates for lung function decline and pulmonary exacerbation onset based on separate and joint models * In the Bayesian sense, a 95% CI that excludes zero indicates statistical significance at the 0.05 level. Parameters b 1 − b 6 are regression coefficients of the cubic truncated power functions defined in "Appendix 1" in an older birth cohort was associated with higher FEV 1 .
Although results were consistent between the two models in Table 1, estimates from the joint model tended to be lower for some of the covariate effects. This could be attributable to incorporating the effect of pulmonary exacerbation through the shared intercept and slope terms in the joint modeling. Both models had similar variance component estimates and indicated substantial heterogeneity in the FEV 1 response, in terms of measurement error, as well as between and within subjects. We examined the smoothed posterior estimates of individual FEV 1 obtained from the longitudinal submodel. The estimated curves of FEV 1 for adolescents and young adults who were observed with pulmonary exacerbation showed more rapid decline than those who had not experienced a pulmonary exacerbation (Fig. 2).

Time-to-event submodel for pulmonary exacerbation
Both Weibull event models indicated that the hazard of pulmonary exacerbation significantly increased over age. Neither model suggested that age at entry was a significant factor (zero was an element of each 95% CI). The separate model indicates a negative association between age and risk of pulmonary exacerbation and the joint model suggests that the association was positive; however, these associations were not statistically significant. Both models imply that being male, having higher FEV 1 at entry, and belonging to one of the earlier birth cohorts were associated with decreased risk of pulmonary exacerbation; however, being born between 1989 and 1994 corresponded to a decreased risk of pulmonary exacerbation, compared to the youngest birth cohort (those born after 1994). Lower SES was associated with an increase in the hazard of pulmonary exacerbation.

Submodel links
Parameter estimates from the joint model indicate that lung function trajectory is associated with risk of having a pulmonary exacerbation. Formally, the intercept and slope parameters corresponding to the shared random effects in the joint model (Table 1, θ 0 and θ 1 estimates, respectively) imply that higher values along the FEV 1 trajectory correspond to a decreased hazard of having a pulmonary exacerbation (HR = exp θ 0 = 0.991 for every 1% predicted increase in FEV 1 , 95% CI 0.986-0.996; HR per one SD increase in FEV 1 = 0.856, 95% CI 0.781-0.937) and positive rates of change or improvements in the FEV 1 trajectory also correspond to decreased hazard (HR = exp θ 1 = 0.702 for every increase of one percentage point (1% predicted) in the rate of change per year in FEV 1 , 95% CI 0.663-0.742; HR per one SD increase in FEV 1 rate of change = 0.566, 95% CI 0.516-0.619). This is clinically reasonable, since a higher level of FEV 1 represents better lung function; patients with FEV 1 measurements that are low or with more rapid decline would be expected to have a higher hazard of pulmonary exacerbation.

Model comparisons
Obvious differences between separate and joint modeling can be graphically observed. We investigated patients with observed FEV 1 trajectories but unknown pulmonary exacerbation event times. We selected two such patients, Patient A and Patient B and compared their estimated posterior median pulmonary exacerbation event time distributions using an established approach [20]. Although neither patient had an observed pulmonary exacerbation, they were censored at ages 16.8 and 17, respectively. Patient A's FEV 1 trajectory began relatively high and remained steady, but Patient B had a trajectory that started relatively low and declined over time (Fig. 3a, b). We compared the median time to pulmonary exacerbation for the two selected patients using posterior distributions from the separate pulmonary exacerbation event model and the two joint models. The predicted median time to pulmonary exacerbation for Patient A was younger under the separate model, compared to each joint model (Fig. 3c). The semiparametric joint model with cubic splines indicated that pulmonary exacerbation onset would occur at an older age than estimated for the joint model with only cubic polynomials and no splines; however, the distributions for the joint model parameters had substantial overlap. The two joint models were in closer agreement regarding time to pulmonary exacerbation for Patient B, but the estimate from the separate pulmonary exacerbation model projected a later event time (Fig. 3d). This is due to joint model's accounting of the correlation between the longitudinal and survival data. The semiparametric mixed effects joint model had the best fit of all models considered based on total DIC (Table 2), followed by the cubic polynomial model without splines, then the separate modeling of the longitudinal and event processes. The association between submodels through the common random intercepts and slopes in each joint model yielded a slight decrease in DIC 1 for the longitudinal submodel and a substantial decrease in DIC 2 for the survival submodel, compared to separate modeling.
In comparing the two joint models, which only differed based on inclusion of splines in the longitudinal submodel, posterior estimates of age-related FEV 1 progression are quite similar until adolescence (Fig. 4a). Based on rates of change, the semiparametric model estimates more rapid decline into early adulthood (Fig. 4b).
Although the cubic regression closely resembles the spline-based estimates, it shows an upward trend in FEV 1 at the later range of age. Patients attained their most rapid decline at 17.6 (6.1) years; median (IQR) estimates   Figure 4b indicates the cubic regression provided poorer fit and unrealistic estimates of the trend in FEV 1 decline at the later range of age. The improvement afforded by the semiparametric model was likely due to the localized splines, which captured more variability in FEV 1 than cubic polynomial terms alone, and in particular allowed for a more reasonable fit at later ages. As a sensitivity analysis, we compared the joint model to a "two-stage" version of the joint model in "Appendix 5". The two-stage model yielded similar estimates, and detected the same significant effects in the random intercept and slope as the joint model. Based on DIC, the joint model had a better fit to the data than the two-stage model, while the two-stage model was superior to the separate modeling. It demonstrates that the joint model not only allows uncertainty in the random effects to carry through to the event model, but also allows the estimation of random effects to depend on both FEV 1 decline and the occurrence of pulmonary exacerbation.

Discussion
In this paper, we described a flexible joint modeling technique aimed at analyzing long sequences of longitudinal and time-to-event data, and used it to simultaneously characterize the nonlinear progression of FEV 1 and assess risk of pulmonary exacerbation events for individual CF patients. We demonstrated how this approach could be used to inform patient management regarding rapid decline in lung function and assessment of pulmonary exacerbation risk over time. On the population level, we identified clinical and demographic risk factors associated with more rapid decline in FEV 1 and earlier onset of pulmonary exacerbation, which could be used to target subpopulations at increased risk of rapid decline or pulmonary exacerbation. Translating these individualized results into clinical care is important because, once a pulmonary exacerbation occurs, it is possible that patients will not recover to their previous FEV 1 levels [5].
In addition to individualized estimates of pulmonary exacerbation risk that account for patient-specific longitudinal trajectories and sources of variation in FEV 1 , this study examined the effect of established risk factors for lung function decline and onset of pulmonary exacerbation. We found significant effects of birth cohort, socioeconomic status and infections with MRSA, Pa and B. cepacia on FEV 1 decline, as well as significant effects of gender, socioeconomic status and birth cohort on pulmonary exacerbation risk. These findings corroborate what has been shown in separate modeling of FEV 1 and pulmonary exacerbation onset [2,5,6]; thus, findings from joint and separate modeling of pulmonary exacerbation appear consistent on the levels of population and subpopulations stratified by risk factors. In addition, this study is the first to our knowledge that has implemented a joint modeling approach to examine the relationship between FEV 1 decline and pulmonary exacerbation events. Through this approach, further insights are gained on an individual patient, the level at which joint and separate modeling results appear to diverge.
Our findings have a number of implications for the epidemiologic study of CF lung disease. Longitudinal modeling of FEV 1 , joint or separate, shows that there is an inverse relationship between starting FEV 1 and rate of decline (Table 1), which corresponds to a correlation of − 0.60 between initial FEV 1 at age 6 years (the age at which patients typically begin performing reliable pulmonary function testing) and slope estimates in our selected model. This association echoes the "ceiling effect" discovered in a previous study of pediatric CF registry data [6]. Our model setup allows us to estimate the population-and subject-specific progression of lung function decline, with the latter being important for personalized medicine. We demonstrated through a specific example that a patient with a less severe FEV 1 trajectory is predicted to be free from a pulmonary exacerbation much longer than a patient with a more severe trajectory. It is likely that these differences among the separate and joint models regarding median pulmonary exacerbation time are attributable to the joint model's correctly accounting for the correlation between the longitudinal and survival data. In terms of model fit, including the time-to-event submodel for pulmonary exacerbation appeared to have a greater impact than including the splines to accommodate nonlinear FEV 1 progression, although both features improved fit (Table 2). Furthermore, having the localized differences fitted from the FEV 1 trajectory has been shown to offer more insights into rate of decline [2,3], particularly in examining the first derivative estimates from our study (Fig. 4). Higher-order random effects could be included to estimate more complex changes over time. Adding a frailty term to Eq. (2) did not improve the total DIC (results not shown). Similar findings have been reported regarding the fits of models with and without the frailty term, all other portions of the model being the same [17].
This study has some limitations that should be considered, which are related to survivor and delayed entry biases, reference equations for lung function and lack of a clear definition of pulmonary exacerbation. Our joint modeling approach corresponds to a missing not at random mechanism as described on pg. 90 by Rizopoulos [21]; the pulmonary exacerbation event process, which is modeled through the random effects, corresponds to the attrition mechanism in a drop-out model. There is a clear birth cohort effect on results from both separate and joint models, which has been noted in other CF registry studies [2,14,22] and could be reflective of advancements in care that were largely unavailable to older patients or survivor bias. Our inclusion of birth cohort in the modeling is only a partial adjustment for left truncation bias, as the correlation between age at registry entry and birth cohort will dictate the extent to which this approach combats left truncation bias. Current WinBUGS implementation is not flexible enough to specify the likelihood for each patient as being conditional on his or her entry time into the registry, which would be necessary to model the delayed entry; we refer to Crowther and colleagues' Eq. (5) [23]. Although there are alternative R packages, such as 'rstan' [24], for flexible likelihood formulations to include delayed entry, estimation can be slow with a large number of patients. Faster computing algorithms would be needed to practically implement this model in large patient registry studies. FEV 1 trajectories, which are modeled based on % predicted values over age, may differ according to the type of spirometry reference equation applied to the raw FEV 1 data, which is expressed in liters. In terms of fitting the FEV 1 trajectory, CFFPR data for patients aged 8-17 years taken from 2013 was analyzed in another study; Wang and Hankinson equations yielded higher median FEV 1 values, compared to values obtained using the Global Lung Initiative equations [25]. Monitoring changes in FEV 1 , by contrast, appear to be less susceptible to this effect. There is no standard definition for what constitutes a pulmonary exacerbation; however, CFFPR studies typically infer its occurrence based on intravenous treatment with antibiotics documented in the registry [10]. This definition overlooks less severe exacerbations that impact lung function but do not warrant intravenous antibiotics. Sensitivity to this definition could be assessed with additional data from an individual CF center, which would likely have documentation on mild pulmonary exacerbation events. These events periodically occur throughout the course of CF; thus, extensions to the joint modeling approach presented here could be used to examine risk of recurrence (Additional files 1, 2).
In conclusion, we have utilized novel statistical modeling of data from a national patient registry to provide more realistic estimates of the FEV 1 trajectory and individualized assessment of pulmonary exacerbation risk in patients with CF. Through the Bayesian approach implemented here using existing software, posterior predictive distributions of this model could be used to aid clinicians in estimating risk of pulmonary exacerbation and rate of lung function decline for individual patients. This approach could be extended to a multivariate joint model [26], in which temporal associations of the evolutions of infections and other characteristics are assessed in conjunction with lung function and exacerbations in CF patients. Furthermore, this joint modeling approach could be used to characterize lung function decline in other diseases and disorders, and to identify subgroups of individuals who may benefit from novel therapeutics being tested in clinical trials.
In the joint modeling literature, various two-stage methods have been proposed. Here we compare the predictive power of a joint model with that of a simple (naïve) two-stage model. The naïve two-stage model consists of: (Stage I) longitudinal model for FEV 1 ; and (Stage II) Weibull model for pulmonary exacerbation with random effects from Stage I as covariates. To be consistent with the joint model (Model III), we use the same semiparametric mixed effects modeling for FEV 1 . The posterior means of subject-specific random effects (i.e., random intercept Ũ 0i and slope Ũ 1i ) from Stage I are included as covariates in the Weibull model in Stage II. The parameters of the two-stage model are given the same priors as the joint model, and are estimated using the fully Bayesian-MCMC (WinBUGS) approach. The results of Table 3 Computing time for each model in WinBUGS * The data set contains 7672 individuals, and the total number of longitudinal observations (FEV 1 measurements) is 36,051. For each model, we run 40,000 iterations, discard the first 20,000 as a warm-up phase and thin by 10