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.

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, 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

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.

Methods

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 patients in the CFFPR indexed by . Let represent lung function measured for the ith patient at the jth time, represented as (corresponding to patient age, in years), . The longitudinal submodel consists of a linear mixed model with random effects:where the population-level mean response is modeled as a continuous function of time, which combines the time-related fixed effects terms and regression splines15; refers to the subject-level fixed effects. The vector represents static and time-varying covariates defined above, and the vector contains their corresponding regression coefficients. The expression incorporates random effects that describe how subject-specific true FEV1 levels deviate from their expected behavior, where the vector corresponds to subject-specific random slopes and intercepts. We used the usual Laird-Ware form8, , corresponding to ; this specification allows individuals to have varying baseline FEV1 measurements and different rates of change in FEV1 over time. Lastly, the term denotes zero-mean Gaussian measurement error.

We considered two nonlinear representations of , 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 . 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 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 , where and The hazard at time for the ith patient iswhich monotonically increases with time if , decreases if , and reduces to the exponential hazard and remains constant if . The vectors and represent possibly time-dependent covariates and the corresponding regression coefficients. Covariates need not have elements in common with in the longitudinal model as shown in Eq. (1). Notice that 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 in Eq. (1), corresponds to patient-specific random intercepts and slopes, but has distinct regression parameter coefficients, and . If none of the covariates vary over time, 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 , i.e., , where follows a bivariate Gaussian distribution, while , independent of

Submodel links

Pulmonary exacerbation event times were associated with longitudinal FEV1 measurements through stochastic dependence between and from Eqs. (1)–(2) by assuming: 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 , which is shared between the longitudinal and time-to-event submodels. Here is assumed to follow a bivariate normal distribution where has diagonal entries , and off-diagonal entries . The joint model with (4) allows both the random intercept and slope , 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 and , where and are the association parameters estimated ordinarily from the HRs and , and and are the respective SDs of intercept and slope 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.

Results

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

Parameter	Separate model (Model I)		Joint model (Model III)
	Posterior mean	95% HPD CI	Posterior mean	95% HPD CI
Semiparametric longitudinal submodel for FEV 1
Curve
β_0, intercept	18.57	[16.33, 20.68]	16.89	[14.63, 20.39]
β_1, age	2.846	[2.545, 3.143]	3.121	[2.74, 3.382]
β_2, age^2	− 0.1295	[− 0.1393, − 0.1197]	− 0.1367	[− 0.1481, − 0.124]
β_3, age^3	− 0.0012	[− 0.0014, − 0.0009]	− 0.0016	[− 0.0017, − 0.0014]
b_1	0.0058	[0.0044, 0.0077]	0.0073	[0.0058, 0.0084]
b_2	0.0083	[0.0076, 0.0089]	0.0103	[0.0095, 0.0111]
b_3	− 0.0143	[− 0.0187, − 0.0114]	− 0.0206	[− 0.0248, − 0.0167]
b_4	− 0.0023	[− 0.0056, 0.0027]	0.0029	[− 0.001, 0.0085]
b_5	0.0041	[− 0.0015, 0.0101]	0.0015	[− 0.0033, 0.007]
b_6	− 0.0011	[− 0.012, 0.0102]	0.0008	[− 0.0095, 0.0101]
Baseline FEV1	0.6828	[0.6722, 0.6933]	0.6834	[0.6676, 0.6945]
Male	− 0.1665	[− 0.6403, 0.3137]	− 0.1650	[− 0.6054, 0.2971]
Birth cohort
< 1981	11.37	[9.651, 12.97]	12.39	[10.89, 13.95]
1981–1988	3.766	[2.749, 4.839]	4.936	[3.995, 5.916]
1989–1994	− 0.1912	[− 0.698, 0.3584]	0.0608	[− 0.4059, 0.5908]
> 1994 (reference)
Lower SES	− 0.6318	[− 1.075, − 0.1686]	− 0.5984	[− 1.062, − 0.1119]
CFRD	− 0.4528	[− 0.765, − 0.1465]	− 0.3695	[− 0.6731, − 0.0633]
MRSA	− 0.8209	[− 1.033, − 0.6121]	− 0.8114	[− 1.016, − 0.6134]
B. cepacia	− 0.8318	[− 1.531, − 0.1367]	− 0.7705	[− 1.441, − 0.1022]
Pa	− 0.5661	[− 0.6963, − 0.4314]	− 0.5673	[− 0.6971, − 0.4373]
Sources of variation
σ^2 (measurement error)	59.29	[58.83, 59.76]	59.24	[58.77, 59.73]
σ02 (between patients, intercept)	324.9	[307.6, 342.2]	326.8	[311, 343.8]
σ12 (between patients, slope)	2.566	[2.441, 2.694]	2.602	[2.484, 2.727]
σ_01 (covariance, intercept and slope)	− 25.51	[− 26.82, − 24.16]	− 25.76	[− 27.11, − 24.49]
Weibull event submodel for pulmonary exacerbation
k	2.467	[2.395, 2.537]	2.564	[2.487, 2.634]
Intercept	− 3.705	[− 3.988, − 3.416]	− 3.568	[− 3.804, − 3.296]
θ_0 (random intercept)	–	–	− 0.0086	[− 0.0137, − 0.0036]
θ_1 (random slope)	–	–	− 0.3532	[− 0.4104, − 0.2979]
Baseline age	− 0.0021	[− 0.0187, 0.0147]	0.0023	[− 0.0131, 0.0176]
Baseline FEV1	− 0.0158	[− 0.0178, − 0.0138]	− 0.0198	[− 0.0218, − 0.0179]
Male	− 0.2973	[− 0.3618, − 0.2314]	− 0.2940	[− 0.3606, − 0.2267]
Birth cohort
< 1981	− 0.3785	[− 0.7538, − 0.0028]	− 0.5348	[− 0.9056, − 0.1726]
1981–1988	− 0.1805	[− 0.3947, 0.0182]	− 0.2348	[− 0.4405, − 0.0339]
1989–1994	0.1431	[0.031, 0.2544]	0.1261	[0.0178, 0.2335]
> 1994 (reference)
Lower SES	0.1173	[0.049, 0.1872]	0.1059	[0.0382, 0.1804]

* In the Bayesian sense, a 95% CI that excludes zero indicates statistical significance at the 0.05 level. Parameters are regression coefficients of the cubic truncated power functions defined in “Appendix 1”

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

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, and estimates, respectively) imply that higher values along the FEV1 trajectory correspond to a decreased hazard of having a pulmonary exacerbation (HR =  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 =  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

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

Model	Submodel 1	Submodel 2	DIC1	DIC2	D¯	p D	DICtotal
Separate model
I	Semiparametric mixed effects	Weibull	952,566	21,110.9	962,619	11,057.7	973,676
Joint model
II	Cubic mixed effects	Weibull	952,915	20,712.5	962,485	11,142.6	973,628
III	Semiparametric mixed effects	Weibull	952,409	20,605.4	961,923	11,091.4	973,015

Each of DIC1 and DIC2 is the sum of the posterior mean deviance () and the effective number of parameters (); DICtotal obtained by summing DIC1 and DIC2 from separate or submodels (see “Appendix 4” for details). Lower values indicate better model fit

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

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.

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

Additional files

Additional file 1. R file for running WinBUGS from R.

Additional file 2. WinBUGS file for analysis implementation.

Table 3

Computing time for each model in WinBUGS

Model	Submodel 1	Submodel 2	Computing time
Separate model
I	Semiparametric mixed effects	Weibull	9 h 57 min
Joint model
II	Cubic mixed effects	Weibull	8 h 45 min
III	Semiparametric mixed effects	Weibull	11 h 14 min

* The data set contains 7672 individuals, and the total number of longitudinal observations (FEV1 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, yielding a total of 2000 posterior samples

Table 4

Posterior estimates under two-stage model and joint model for lung function decline and pulmonary exacerbation onset

Parameter	Two-stage model		Joint model (Model III)
	Posterior mean	95% HPD CI	Posterior mean	95% HPD CI
Semiparametric longitudinal submodel for FEV 1
Curve
β_0, intercept	18.15	[15.54, 20.95]	16.89	[14.63, 20.39]
β_1, age	2.86	[2.477, 3.189]	3.121	[2.74, 3.382]
β_2, age2	− 0.1306	[− 0.1418, − 0.117]	− 0.1367	[− 0.1481, − 0.124]
β_3, age3	− 0.0011	[− 0.0013, − 0.0009]	− 0.0016	[− 0.0017, − 0.0014]
b_1	0.0051	[0.0039, 0.0061]	0.0073	[0.0058, 0.0084]
b_2	0.0108	[0.0097, 0.0119]	0.0103	[0.0095, 0.0111]
b_3	− 0.0177	[− 0.0226, − 0.0133]	− 0.0206	[− 0.0248, − 0.0167]
b_4	− 0.0001	[− 0.0044, 0.0065]	0.0029	[− 0.001, 0.0085]
b_5	0.0033	[− 0.0024, 0.0083]	0.0015	[− 0.0033, 0.007]
b_6	− 0.0009	[− 0.0114, 0.0095]	0.0008	[− 0.0095, 0.0101]
Baseline FEV1	0.686	[0.6728, 0.6978]	0.6834	[0.6676, 0.6945]
Male	− 0.15	[− 0.582, 0.2736]	− 0.1650	[− 0.6054, 0.2971]
Birth cohort
< 1981	11.43	[9.837, 13.29]	12.39	[10.89, 13.95]
1981–1988	4.0	[2.992, 4.978]	4.936	[3.995, 5.916]
1989–1994	− 0.1924	[− 0.7038, 0.3305]	0.0608	[− 0.4059, 0.5908]
> 1994 (reference)
Lower SES	− 0.623	[− 1.091, − 0.175]	− 0.5984	[− 1.062, − 0.1119]
CFRD	− 0.4613	[− 0.7895, − 0.1506]	− 0.3695	[− 0.6731, − 0.0633]
MRSA	− 0.816	[− 1.026, − 0.6022]	− 0.8114	[− 1.016, − 0.6134]
B. cepacia	− 0.7964	[− 1.465, − 0.1291]	− 0.7705	[− 1.441, − 0.1022]
Pa	− 0.5643	[− 0.6953, − 0.4345]	− 0.5673	[− 0.6971, − 0.4373]
Sources of variation
σ^2 (measurement error)	59.29	[58.83, 59.78]	59.24	[58.77, 59.73]
σ02 (between patients, intercept)	324.1	[307.6, 341.3]	326.8	[311, 343.8]
σ12 (between patients, slope)	2.564	[2.443, 2.69]	2.602	[2.484, 2.727]
σ_01 (covariance, intercept and slope)	− 25.46	[− 26.81, − 24.18]	− 25.76	[− 27.11, − 24.49]
Weibull event submodel for pulmonary exacerbation
k	2.524	[2.456, 2.601]	2.564	[2.487, 2.634]
intercept	− 3.509	[− 3.802, − 3.26]	− 3.568	[− 3.804, − 3.296]
θ_0 (random intercept)	− 0.0098	[− 0.0149, − 0.0048]	− 0.0086	[− 0.0137, − 0.0036]
θ_1 (random slope)	− 0.3273	[− 0.3825, − 0.2728]	− 0.3532	[− 0.4104, − 0.2979]
Baseline age	0.0002	[− 0.0136, 0.0161]	0.0023	[− 0.0131, 0.0176]
Baseline FEV1	− 0.0193	[− 0.0211, − 0.0175]	− 0.0198	[− 0.0218, − 0.0179]
Male	− 0.2919	[− 0.36, − 0.2194]	− 0.2940	[− 0.3606, − 0.2267]
Birth cohort
< 1981	− 0.5046	[− 0.8512, − 0.1807]	− 0.5348	[− 0.9056, − 0.1726]
1981–1988	− 0.2297	[− 0.4382, − 0.0415]	− 0.2348	[− 0.4405, − 0.0339]
1989–1994	0.1171	[0.0086, 0.2197]	0.1261	[0.0178, 0.2335]
>1994 (reference)
Lower SES	0.1091	[0.0369, 0.1799]	0.1059	[0.0382, 0.1804]

Table 5

Model comparison results

Model	Submodel 1	Submodel 2	DIC1	DIC2	D¯	p_D	DICtotal
Two-stage model	Semiparametric mixed effects	Weibull	952,553	20,740.6	951,177	11,063.7	973,294
Joint model (Model III)	Semiparametric mixed effects	Weibull	952,409	20,605.4	961,923	11,091.4	973,015