Skip to main content

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



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.


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.


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.


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.


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, FEV1), which is the primary marker of disease severity in individuals with CF. FEV1 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 FEV1 over age, but initial FEV1 and rate of decline vary between patients. Linear mixed modeling is an established approach to estimate age-related FEV1 progression. Historically, CF epidemiologic studies have not accounted for nonlinearity in the FEV1 trajectory; however, recent approaches have included piecewise polynomials, either in the form of regression splines or change-point models, to estimate decline in FEV1 [2, 3]. These approaches have shown that FEV1 decline is variable with maximal loss occurring in adolescence and early adulthood.

Fig. 1
figure 1

FEV1 versus age (in years) for 100 randomly selected patients with cystic fibrosis from the U.S. CFFPR, 2003–2011. Points have been connected over age for the 50 patients who were observed with pulmonary exacerbation (in blue)

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 (pre-pulmonary 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 FEV1 adjusted for survivor bias. They combined a linear mixed effects model for longitudinal FEV1 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 FEV1 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 FEV1 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 FEV1 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, FEV1 % predicted (hereafter, FEV1), 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 FEV1. 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 FEV1 measurement was available. Covariates of interest were selected from the literature on FEV1 [13, 14] and pulmonary exacerbation [5, 6] as separate models, and included age, gender, initial FEV1 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 FEV1 and a Weibull model to fit pulmonary exacerbation data.

Longitudinal submodel for FEV1

Both the separate modeling approach and the joint modeling approach require specification of a model for the longitudinal measure of FEV1. 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, \ldots ,N\). Let \(FEV1_{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, \ldots ,n_{i}\). The longitudinal submodel consists of a linear mixed model with random effects:

$$FEV1_{ij} = f\left( {age_{ij} } \right) + m_{i} \left( {age_{ij} } \right) + W_{1i} \left( {age_{ij} } \right) + \varepsilon_{ij} ,$$

where the population-level mean response \(f\left( {age_{ij} } \right)\) is modeled as a continuous function of time, which combines the time-related fixed effects terms and regression splines15; \(m_{i} \left( {age_{ij} } \right) = {\mathbf{x}}_{1i}^{{\prime }} \left( {age_{ij} } \right)\varvec{\alpha}_{1}\) refers to the subject-level fixed effects. The vector \({\mathbf{x}}_{1i} \left( {age_{ij} } \right)\) represents static and time-varying covariates defined above, and the vector \(\varvec{\alpha}_{1}\) contains their corresponding regression coefficients. The expression \(W_{1i} \left( {age_{ij} } \right) = {\mathbf{z}}_{1i}^{{\prime }} \left( {age_{ij} } \right){\mathbf{U}}_{i}\) incorporates random effects that describe how subject-specific true FEV1 levels deviate from their expected behavior, where the vector \({\mathbf{U}}_{i}\) corresponds to subject-specific random slopes and intercepts. We used the usual Laird-Ware form8, \(W_{1i} \left( {age_{ij} } \right) = U_{0i} + U_{1i} age_{ij}\), corresponding to \({\mathbf{z}}_{1i} \left( {age_{ij} } \right) = \left( {1,age_{ij} } \right)^{{\prime }}\); this specification allows individuals to have varying baseline FEV1 measurements and different rates of change in FEV1 over time. Lastly, the term \(\varepsilon_{ij} \sim N\left( {0,\sigma_{\varepsilon }^{2} } \right)\) denotes zero-mean Gaussian measurement error.

We considered two nonlinear representations of \(f\left( {age_{ij} } \right)\), 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 FEV1 measurements. Second, we considered only global cubic polynomials to represent population-level FEV1 decline by excluding the basis functions in \(f\left( {age_{ij} } \right)\). 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 FEV1, and has been used previously to model FEV1 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 \(t_{i} \sim {\text{Weibull}}\left( {k, \mu_{i} \left( t \right)} \right)\), where \(\mu_{i} \left( t \right) = { \exp }\left\{ {{\mathbf{x}}_{2i}^{{\prime }} \left( t \right)\varvec{\alpha}_{2} + W_{2i} } \right\}\) and \(k > 0.\) The hazard at time \(t\) for the ith patient is

$$h_{i} (t) = kt^{(k - 1)} \exp \left( {{\mathbf{x}}_{2i}^{{\prime }} (t)\varvec{\alpha}_{2} + W_{2i} } \right),$$

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 \({\mathbf{x}}_{2i} \left( t \right)\) and \(\varvec{\alpha}_{2}\) represent possibly time-dependent covariates and the corresponding regression coefficients. Covariates \({\mathbf{x}}_{2i} \left( t \right)\) need not have elements in common with \({\mathbf{x}}_{1i} \left( {age} \right)\) 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} \left( {age} \right)\) in Eq. (1), \(W_{2i} = \theta_{0} U_{0i} + \theta_{1} U_{1i}\) corresponds to patient-specific random intercepts and slopes, but has distinct regression parameter coefficients, \(\theta_{0}\) and \(\theta_{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 \(W_{2}\), i.e., \(W_{2i} = \theta_{0} U_{0i} + \theta_{1} U_{1i} + U_{3i}\), where \(\left( {U_{0i} ,U_{1i} } \right)^{{\prime }}\) follows a bivariate Gaussian distribution, while \(U_{3i} \sim N\left( {0,\sigma_{3}^{2} } \right)\), independent of \(\left( {U_{0i} ,U_{1i} } \right)^{{\prime }} .\)

Submodel links

Pulmonary exacerbation event times were associated with longitudinal FEV1 measurements through stochastic dependence between \(W_{1i}\) and \(W_{2i}\) from Eqs. (1)–(2) by assuming:

$$W_{1i} \left( {\text{age}} \right) = U_{0i} + U_{1i} \,{\text{age}};$$
$$W_{2i} \left( t \right) = W_{2i} = \theta_{0} U_{0i} + \theta_{1} U_{1i} .$$

The subject-specific random intercept and slope, depicting varying initial values and rates of FEV1 decline for each patient after accounting for the nonlinear FEV1 trajectory, are contained in the vector \({\mathbf{U}}_{i}\), which is shared between the longitudinal and time-to-event submodels. Here \({\mathbf{U}}_{i}\) is assumed to follow a bivariate normal distribution \(N\left( {0,{\varvec{\Sigma}}} \right),\) where \({\varvec{\Sigma}}\) has diagonal entries \(\sigma_{0}^{2}\), \(\sigma_{1}^{2}\) and off-diagonal entries \(\sigma_{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 FEV1 trajectories from the population-level FEV1 curve enter the pulmonary exacerbation model in the form of random intercepts and slopes. If this type of association exists between the longitudinal FEV1 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 FEV1 (intercept) and rate of change in FEV1 (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 \left\{ {\theta_{0} *\sigma_{0} } \right\}\) and \(\exp \left\{ {\theta_{1} *\sigma_{1} } \right\}\), where \(\theta_{0}\) and \(\theta_{1}\) are the association parameters estimated ordinarily from the HRs \(\exp \left\{ {\theta_{0} } \right\}\) and \(\exp \left\{ {\theta_{1} } \right\}\), and \(\sigma_{0}\) and \(\sigma_{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 DIC1 and DIC2, 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.


Study cohort

There were 7672 individuals who met inclusion criteria, yielding a total of 136,051 FEV1 measurements. Median (IQR) follow-up was 5.8 (4.0–7.7) years; 3349 (43.7%) of the patients experienced a pulmonary exacerbation. Median age at entry was 10.8 (6.9–16.6) years; 4298 (56%) of the patients were male; 2690 (35%) had lower SES (only state/federal or no insurance). In total, 1036 (13.5%), 1176 (15.3%), 2144 (27.9%) and 3316 (43.2%) of the patients were born before 1981, between 1981 and 1988, between 1989 and 1994, and after 1994, respectively. Few patients had infection with B. cepacia (2.5%), but 1589 (20.7%) had Pa infection at baseline. Mean (SD) FEV1 at entry was 92.3 (19.6) % predicted.

Longitudinal submodel for FEV1

Both the separate and joint modeling of longitudinal FEV1 indicated that decline was nonlinear over age (Table 1, Semiparametric longitudinal submodel). In both models, having higher FEV1 at entry corresponded to higher FEV1 over time; lower SES, diagnosis with CFRD and the presence of infections were associated with lower FEV1 over age. The younger birth cohorts appeared to have similar FEV1 progression, while being in an older birth cohort was associated with higher FEV1. 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 FEV1 response, in terms of measurement error, as well as between and within subjects.

Table 1 Posterior estimates for lung function decline and pulmonary exacerbation onset based on separate and joint models

We examined the smoothed posterior estimates of individual FEV1 obtained from the longitudinal submodel. The estimated curves of FEV1 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).

Fig. 2
figure 2

Smoothed posterior estimates of individual FEV1 for the 100 randomly selected patients presented in Fig. 1. Red lines are the smoothed estimates for individuals who were observed with pulmonary exacerbation, while black lines are the smoothed estimates for individuals who were not observed with pulmonary exacerbation

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 FEV1 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, \(\theta_{0}\) and \(\theta_{1}\) estimates, respectively) imply that higher values along the FEV1 trajectory correspond to a decreased hazard of having a pulmonary exacerbation (HR = \(\exp \left( {\hat{\theta }_{0} } \right) =\) 0.991 for every 1% predicted increase in FEV1, 95% CI 0.986–0.996; HR per one SD increase in FEV1 = 0.856, 95% CI 0.781–0.937) and positive rates of change or improvements in the FEV1 trajectory also correspond to decreased hazard (HR = \(\exp \left( {\hat{\theta }_{1} } \right) =\) 0.702 for every increase of one percentage point (1% predicted) in the rate of change per year in FEV1, 95% CI 0.663–0.742; HR per one SD increase in FEV1 rate of change = 0.566, 95% CI 0.516–0.619). This is clinically reasonable, since a higher level of FEV1 represents better lung function; patients with FEV1 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 FEV1 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 FEV1 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.

Fig. 3
figure 3

Observed data and estimated FEV1 over time (a), the rate of decline in FEV1 over time (b) and estimated posteriors of the median time to pulmonary exacerbation for two patients (c) and (d)

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 DIC1 for the longitudinal submodel and a substantial decrease in DIC2 for the survival submodel, compared to separate modeling.

Table 2 Model comparison results

In comparing the two joint models, which only differed based on inclusion of splines in the longitudinal submodel, posterior estimates of age-related FEV1 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 FEV1 at the later range of age. Patients attained their most rapid decline at 17.6 (6.1) years; median (IQR) estimates are 17.5 (IQR 13.7–18) years. The rate of decline continues at a slower pace until patients approximately reach 30 years of age. Figure 4b indicates the cubic regression provided poorer fit and unrealistic estimates of the trend in FEV1 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 FEV1 than cubic polynomial terms alone, and in particular allowed for a more reasonable fit at later ages.

Fig. 4
figure 4

Posterior estimates obtained from joint models with semiparametric (solid line) and cubic (dashed line) submodels, and separate model with semiparametric (dash-dotted line) submodel of decline (a), and the rate of decline (b) in FEV1 over age (in years) for the overall population

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 FEV1 decline and the occurrence of pulmonary exacerbation.


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 FEV1 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 FEV1 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 FEV1 levels [5].

In addition to individualized estimates of pulmonary exacerbation risk that account for patient-specific longitudinal trajectories and sources of variation in FEV1, 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 FEV1 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 FEV1 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 FEV1 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 FEV1, joint or separate, shows that there is an inverse relationship between starting FEV1 and rate of decline (Table 1), which corresponds to a correlation of − 0.60 between initial FEV1 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 FEV1 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 FEV1 progression, although both features improved fit (Table 2). Furthermore, having the localized differences fitted from the FEV1 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. FEV1 trajectories, which are modeled based on % predicted values over age, may differ according to the type of spirometry reference equation applied to the raw FEV1 data, which is expressed in liters. In terms of fitting the FEV1 trajectory, CFFPR data for patients aged 8–17 years taken from 2013 was analyzed in another study; Wang and Hankinson equations yielded higher median FEV1 values, compared to values obtained using the Global Lung Initiative equations [25]. Monitoring changes in FEV1, 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 FEV1 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.


B. cepacia :

Burkholderia cepacia


Cystic fibrosis


Cystic fibrosis foundation patient registry


Cystic fibrosis related diabetes


Methicillin-resistant staphylococcus aureus

Pa :

Pseudomonas aeruginosa


Socioeconomic status


  1. Farrell PM. The prevalence of cystic fibrosis in the European Union. J Cyst Fibros. 2008;7(5):450–3.

    Article  PubMed  Google Scholar 

  2. Szczesniak RD, et al. A semiparametric approach to estimate rapid lung function decline in cystic fibrosis. Ann Epidemiol. 2013;23(12):771–7.

    Article  PubMed  Google Scholar 

  3. Moss A, et al. A comparison of change point models with application to longitudinal lung function measurements in children with cystic fibrosis. Stat Med. 2016;35(12):2058–73.

    Article  PubMed  Google Scholar 

  4. Ferkol T, Rosenfeld M, Milla CE. Cystic fibrosis pulmonary exacerbations. J Pediatr. 2006;148(2):259–64.

    Article  PubMed  Google Scholar 

  5. Sanders DB, et al. Failure to recover to baseline pulmonary function after cystic fibrosis pulmonary exacerbation. Am J Respir Crit Care Med. 2010;182(5):627–32.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Konstan MW, et al. Risk factors for rate of decline in forced expiratory volume in one second in children and adolescents with cystic fibrosis. J Pediatr. 2007;151(2):134–9.

    Article  PubMed  Google Scholar 

  7. Asar O, et al. Joint modelling of repeated measurement and time-to-event data: an introductory tutorial. Int J Epidemiol. 2015;44(1):334–44.

    Article  PubMed  Google Scholar 

  8. Schluchter MD, Konstan MW, Davis PB. Jointly modelling the relationship between survival and pulmonary function in cystic fibrosis patients. Stat Med. 2002;21(9):1271–87.

    Article  PubMed  Google Scholar 

  9. Piccorelli AV, Schluchter MD. Jointly modeling the relationship between longitudinal and survival data subject to left truncation with applications to cystic fibrosis. Stat Med. 2012;31(29):3931–45.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Knapp, E.A., et al., The Cystic Fibrosis Foundation Patient Registry. Design and Methods of a National Observational Disease Registry. Ann Am Thorac Soc, 2016;13(7):1173–79.

  11. Wang X, et al. Pulmonary function between 6 and 18 years of age. Pediatr Pulmonol. 1993;15(2):75–88.

    Article  PubMed  Google Scholar 

  12. Hankinson JL, Odencrantz JR, Fedan KB. Spirometric reference values from a sample of the general US population. Am J Respir Crit Care Med. 1999;159(1):179–87.

    Article  CAS  PubMed  Google Scholar 

  13. Corey M, et al. Longitudinal analysis of pulmonary function decline in patients with cystic fibrosis. J Pediatr. 1997;131(6):809–14.

    Article  CAS  PubMed  Google Scholar 

  14. Taylor-Robinson D, et al. Understanding the natural progression in %FEV1 decline in patients with cystic fibrosis: a longitudinal study. Thorax. 2012;67(10):860–6.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Henderson R, Diggle P, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics. 2000;1(4):465–80.

    Article  CAS  PubMed  Google Scholar 

  16. Tsiatis AA, Davidian M. Joint modeling of longitudinal and time-to-event data: an overview. Stat Sin. 2004;14:809–34.

    Google Scholar 

  17. Fitzmaurice G, Davidian M, Verbeke G, Molenberghs G. Longitudinal data analysis. In: Molenberghs G, Fitzmaurice G, editors. Incomplete data: introduction and overview. Boca Raton: Taylor & Francis Group; 2009. p. 398–400.

    Google Scholar 

  18. Lunn DJT, Best N, Spiegelhalter D. WinBUGS: a Bayesian modelling framework—concepts, structure, and extensibility. Stat Comput. 2000;10:325–37.

    Article  Google Scholar 

  19. Spiegelhalter DJ, Best NGC, Van der Linde A. Bayesian measures of model complexity and fit (with discussion). J R Stat Soc Ser B. 2002;64(4):583–616.

    Article  Google Scholar 

  20. Guo XC, Carlin BP. Separate and joint modeling of longitudinal and event time data using standard computer packages. Am Stat. 2004;58(1):16–24.

    Article  Google Scholar 

  21. Rizopoulos D. Joint models for longitudinal and time-to-event data: with applications in R, vol. xiv., Chapman & Hall/CRC biostatistics seriesBoca Raton: CRC Press; 2012. p. 261.

    Book  Google Scholar 

  22. VanDevanter DR, Pasta DJ, Konstan MW. Improvements in lung function and height among cohorts of 6-year-olds with cystic fibrosis from 1994 to 2012. J Pediatr. 2014;165(6):1091e2–1097e2.

    Article  Google Scholar 

  23. Crowther MJ, et al. Joint modelling of longitudinal and survival data: incorporating delayed entry and an assessment of model misspecification. Stat Med. 2016;35(7):1193–209.

    Article  PubMed  Google Scholar 

  24. Stan Development Team. Rstan: the R interface to Stan, version 2.16.2; 2017. Accessed 9 Nov 2017.

  25. Foundation CF, Cystic fibrosis foundation patient registry. In: 2013 annual report to the center directors. Cystic Fibrosis Foundation: Bethesda; 2014.

  26. Albert PS, Shih JH. An approach for jointly modeling multivariate longitudinal measurements and discrete time-to-event data. Ann Appl Stat. 2010;4(3):1517–32.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Ngo L, Wand M. Smoothing with mixed model software. J Stat Softw. 2004;9:1–54.

    Google Scholar 

  28. Lunn DJ, Thomas A, Best N, Spiegelhalter D. WinBUGS—a Bayesian modelling framework: concepts, structure, and extensibility. Stat Comput. 2000;10:325–37.

    Article  Google Scholar 

  29. Spiegelhalter D, Thomas A, Best N, Gilks W. WinBUGS user manual.

Download references

Authors’ contributions

DL, RDS: Developed the study concept, analysis plan, results summaries and drafted the initial manuscript. DL: Implemented all model analyses. RK: Critically reviewed the manuscript, provided statistical input on the evaluation of the models and interpretation of findings, and contributed to the writing in all sections. JPC: Critically reviewed the manuscript, provided interpretation of the statistical model results into the clinical context, and contributed to the writing in all sections. All authors read and approved the final manuscript.


The authors would like to thank the Cystic Fibrosis Foundation for the use of CF Foundation Patient Registry data to conduct this study. Additionally, we would like to thank the patients, care providers, and clinic coordinators at CF centers throughout the United States for their contributions to the CF Foundation Patient Registry.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The dataset analyzed for this study is not publicly available. Data is available upon request through the Cystic Fibrosis Foundation Patient Registry Comparative Effectiveness Research Committee. Additional information may be available by contacting Restrictions on access to data are to ensure patient privacy for all persons in the Cystic Fibrosis Foundation Patient Registry. Implementation materials in the form of program files are available as additional files accompanying this work.

Consent for publication

Not applicable.

Ethics approval and consent to participate

This retrospective study of the United States Cystic Fibrosis Foundation CF Patient Registry was approved by the Institutional Review Board of the Cincinnati Children’s Hospital Medical Center (IRB # 2015-4515). The Cystic Fibrosis Foundation Patient Registry Committee approved the request for data access and dispensed the data (Request # PRR088).

Financial disclosure statement

The authors have no financial relationships relevant to this article to disclose.


This study was funded by grants K25 HL125954 from the National Heart, Lung and Blood Institute of the National Institutes of Health and R457-CR11 from the Cystic Fibrosis Foundation Research and Development Program.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Rhonda D. Szczesniak.

Additional files


Appendix 1: Regression splines used in the longitudinal submodel

The function \(f\left( t \right)\) is defined as a cubic truncated power basis function for the semiparametric representation:

$$f\left( t \right) = \beta_{0} + \beta_{1} t + \beta_{2} t^{2} + \beta_{3} t^{3} + \mathop \sum \limits_{l = 1}^{L} b_{l} \left( {t - \kappa_{l} } \right)_{ + }^{3} ,$$

where \(\beta_{0} , \beta_{1} , \beta_{2}\) and \(\beta_{3}\) are respectively the coefficients for intercept, linear, quadratic and cubic terms; the term \(\sum\nolimits_{l = 1}^{L} {b_{l} \left( {t - \kappa_{l} } \right)_{ + }^{3} }\) represents basis functions for the cubic splines evaluated at time \(s\), and \(\left( {t - \kappa_{l} } \right)_{ + } = \hbox{max} \left( {t - \kappa_{l} , 0} \right)\). We selected the number and location of knots (\(\kappa_{l} , l = 1, \ldots ,L\)) for the corresponding spline basis functions using the quantile-based approach described in Ngo and Wand [27]. Please see “Appendix 2” for more discussion on knot specification. Six knots used across the [6, 45] age (in years) interval were \(\varvec{\kappa}= \left( {10.78, 15.31, 19.89, 24.96, 30.49, 37.34} \right)\).

The second formulation of the longitudinal submodel defines \(f\left( t \right)\) without the basis functions, which we refer to as the cubic representation:

$$f\left( t \right) = \beta_{0} + \beta_{1} t + \beta_{2} t^{2} + \beta_{3} t^{3}$$

Appendix 2: Knot specification

The main idea of choosing knots is to choose enough knots to resolve the essential structure in the underlying regression function, but to keep the number of knots relatively low considering the computational advantages. A reasonable default is to choose the knots to ensure that there are a fixed number of unique observations (e.g., 4–5) between each knot. However, for large data sets, this may lead to an excessive number of knots, so a maximum number of knots (e.g., 20–40) is recommended. A simple default choice of \(L\) that always works well is:

$$L = \hbox{max} \left\{ {5,\hbox{min} \left( {\frac{1}{4} \times {\text{number of unique}} t_{i}^{{\prime }} s, 35} \right)} \right\}.$$

A reasonable default rule for the knot location is:

$$K_{l} = \{ (l + 1)/(L + 2)\} {\text{th}}\,{\text{sample}}\,{\text{quantile}}\,{\text{of}}\,{\text{the}}\,{\text{unique's}},$$

for \(l = 1, \ldots ,L\). The default choice of knots can be generated in R using the function defined as below:

figure a

Appendix 3: Prior specification and sampling

We used a fully Bayesian version of the joint modeling approach implemented via Markov chain Monte Carlo (MCMC) method to fit the joint model. Although we focus on a particular form, the Bayesian framework allows for easier comparison among many forms of association between \(W_{1}\) and \(W_{2}\). The data analyses are conducted using MCMC through WinBUGS software [28]. R and WinBUGS codes implementing the proposed method are available in Additional files 1 and 2. We selected vague prior distributions in the Bayesian analysis. For each of the regression coefficients of the cubic penalized spline and the subject-level covariate effects \(\varvec{\beta}, \varvec{b}\) and \(\varvec{\alpha}_{1}\), we used a multivariate normal prior with mean vector \(0\) and precision matrix \(0.01 \cdot {\mathbf{I}}\), where \({\mathbf{I}}\) is the identity matrix; and an inverse gamma prior19 \(IG\left( {0.01,0.01} \right)\) for the measurement error variance \(\sigma_{\varepsilon }^{2}\). Similarly, we selected a vague normal prior in the survival submodel, \(\varvec{\alpha}_{1} \sim N\left( {0,\varvec{ }0.01 \cdot {\mathbf{I}}} \right)\). For the random effects \({\mathbf{U}}_{i}\) common to both submodels, we took a non-informative multivariate inverse-Wishart distribution with degrees of freedom equals to 2 (the number of random effect components). Finally, for the association parameters, we set \(\theta_{0} \sim N\left( {0,100} \right)\) and \(\theta_{1} \sim N\left( {0,100} \right)\).

Complete conditional distributions were used to approximate the posterior distributions for the parameters of interest. MCMC diagnostics to assess convergence were performed in accordance with standard techniques [29]. Early samples were discarded based on graphical inspection of burn-in across MCMC iterations. The 95% credible intervals (CIs) were obtained by taking the 2.5th and 97.5th percentiles of the posterior samples post burn-in. For each model, we conducted the MCMC simulation for 40,000 iterations. The first 20,000 draws were discarded as a burn-in phase. The effective iterations were thinned by storing every 10th iteration in order to reduce level of correlation between successive values of the sample. The posterior mean and 95% confidence interval (CI) were calculated with the resulting 2000 values. The computing time for each model is shown in Table 3.

Table 3 Computing time for each model in WinBUGS

Appendix 4: DIC formulation

Each deviance information criterion (DIC) for assessing model fit is calculated as \({\text{DIC}} = D + p_{D}\), with the posterior mean of deviance \(D\) penalized by the effective number of parameters \(p_{D}\) under the Bayesian framework. Since small \(D\) indicates good fit while small \(p_{D}\) indicate a parsimonious model, smaller values indicate better model fit, compared to larger values. We respectively write the components of DIC for the longitudinal and survival submodels as \({\text{DIC}}_{1}\) and \({\text{DIC}}_{2}\), which evaluate the contributions to the overall DIC, denoted as \({\text{DIC}}_{\text{total}}\) in Table 2. All DIC values are output directly from WinBUGS.

Appendix 5: Comparison between two-stage model and joint model

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 FEV1; 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 FEV1. The posterior means of subject-specific random effects (i.e., random intercept \(\tilde{U}_{0i}\) and slope \(\tilde{U}_{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 posterior estimates and model comparison for two-stage model and joint model (Model III) are given in Tables 4 and 5 .

Table 4 Posterior estimates under two-stage model and joint model for lung function decline and pulmonary exacerbation onset
Table 5 Model comparison results

The two-stage model yields similar estimates, and detects the same significant effects in the random intercept and slope as the joint model. However, DIC in Table A3 suggests that the joint model (DICtotal = 973,015) gives a much better fit of the data than the two-stage model (DICtotal = 973,294).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Li, D., Keogh, R., Clancy, J.P. et al. Flexible semiparametric joint modeling: an application to estimate individual lung function decline and risk of pulmonary exacerbations in cystic fibrosis. Emerg Themes Epidemiol 14, 13 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Bayesian analysis
  • Cystic fibrosis
  • Functional data analysis
  • Longitudinal studies
  • Mixed model analysis
  • Pulmonary decline
  • Pulmonary function
  • Joint modeling
  • Registry analyses
  • Spline regression