Multivariate joint modeling to identify markers of growth and lung function decline that predict cystic fibrosis pulmonary exacerbation onset



Attenuated decreases in lung function can signal the onset of acute respiratory events known as pulmonary exacerbations (PEs) in children and adolescents with cystic fibrosis (CF). Univariate joint modeling facilitates dynamic risk prediction of PE onset and accounts for measurement error of the lung function marker. However, CF is a multi-system disease and the extent to which simultaneously modeling growth and nutrition markers improves PE predictive accuracy is unknown. Furthermore, it is unclear which routinely collected clinical indicators of growth and nutrition in early life predict PE onset in CF.


Using a longitudinal cohort of 17,100 patients aged 6–20 years (US Cystic Fibrosis Foundation Patient Registry; 2003–2015), we fit a univariate joint model of lung-function decline and PE onset and contrasted its predictive performance with a class of multivariate joint models that included combinations of growth markers as additional submodels. Outcomes were longitudinal lung function (forced expiratory volume in 1 s of % predicted), percentiles of body mass index, weight-for-age and height-for-age and PE onset. Relevant demographic/clinical covariates were included in submodels. We implemented a univariate joint model of lung function and time-to-PE and four multivariate joint models including growth outcomes.


All five joint models showed that declining lung function corresponded to slightly increased risk of PE onset (hazard ratio from univariate joint model: 0.97, P < 0.0001), and all had reasonable predictive accuracy (cross-validated area under the receiver-operator characteristic curve > 0.70). None of the growth markers alongside lung function as outcomes in multivariate joint modeling appeared to have an association with hazard of PE. Jointly modeling only lung function and PE onset yielded the most accurate (area under the receiver-operator characteristic curve = 0.75) and precise (narrowest interquartile range) predictions. Dynamic predictions were accurate across forecast horizons (0.5, 1 and 2 years) and precision improved with age.


Including growth markers via multivariate joint models did not yield gains in prediction performance, compared to a univariate joint model with lung function. Individualized dynamic predictions from joint modeling could enhance physician monitoring of CF disease progression by providing PE risk assessment over a patient’s clinical course.



Background
Cystic fibrosis (CF) is a chronic lung disease in which death commonly results from respiratory failure (1). A cycle of prolonged and acute drops in lung function marks the clinical course of CF (2), most notably during adolescence and early adulthood (3) (4) (5). Forced expiratory volume in 1 s of % predicted (hereafter, FEV 1 ), is the primary surrogate of lung function and strongest predictor of mortality in the CF population (6). FEV 1 remains an important outcome in clinical trials and the most relevant clinical indicator in monitoring lung function decline (7). It is also of clinical importance to identify an acute respiratory event referred to as a pulmonary exacerbation (PE), which frequently corresponds to sharp decreases in lung function. Diagnosis of a PE depends on various factors corresponding to lung function and nutritional status, often resulting in a patient being hospitalized and administered intravenous antibiotics (8). Further, patients who experience a PE often fail to recover to their corresponding pre-event (i.e., "baseline") levels of FEV 1 (9). For these reasons, it has been desirable in epidemiologic studies to identify risk factors for having a PE.
Although maintaining lung function and minimizing risk of PEs are essential to survival, CF is a multisystem disease in which malnutrition and poor growth have been shown to adversely impact pulmonary function (10). Body mass index percentile (BMIp) is the primary nutritional indicator in CF (11); guidelines published by the CF Foundation recommend maintaining BMIp > 50 th percentile (12).
A recent CF registry analysis found that weight-for-age (WFA) and height-for-age (HFA) percentiles more accurately identify CF patients of poor nutritional status and stunting, respectively, compared to BMI percentiles (13).
These and other studies have collectively demonstrated that multiple causes of disease progression are typically observed from indicators of pulmonary and nutritional decline, yet it remains unclear how the underlying longitudinal processes collectively evolve over time. Joint modeling is a popular statistical approach used to link longitudinal and time-to-event processes and is usually implemented through shared parameters (14). The idea behind these models is that the longitudinal and event processes share common random effects, inducing correlation between the longitudinal and event 4 processes. Specifically, we can construct a mixed-effects model to describe the evolution over time for each longitudinal outcome and use these estimated evolutions as a time-dependent covariate in a time-to-event model.
Novel applications of so-called univariate joint models (refers to a single longitudinal marker being modeled simultaneously with an event outcome) in CF epidemiologic studies have yielded new insights and predictions regarding how changes in FEV 1 are associated with survival (i.e., time to death/lung transplant) (15) (16) (17). Although these applications elucidate CF disease progression into its end stage, little is known about the utility of joint models in assessing onset of PEs; these events initialize earlier in life but could be attributable to downstream impacts. Furthermore, to our knowledge, the association of one or more longitudinal outcomes with onset of PE has not yet been addressed. Univariate joint modeling potentially ignores the fact that different markers of interest, specifically those related to growth and nutrition status, could also have an influence on PE in combination.
Due to computational advancements, such as new software packages (18), joint models have enjoyed a renaissance in the last decade (19). In this paper, we focus on multiple longitudinal outcomes and their potential to predict PE onset in children with CF. Specifically, we propose a multivariate joint model of FEV 1 and markers of nutritional status with time-to-PE. The objective of our study was to construct these models to estimate how commonly collected markers of these longitudinal processes evolve through the clinical course of CF, and to compare results from these models regarding predictive accuracy. We utilize data from the CFF patient registry and provide implementation of the joint models as supplemental material.

Cohort description
The CF cohort for this empirical study consists of patients with longitudinal data recorded from childhood until early adulthood (ages: 6-20 years) in the US Cystic Fibrosis Foundation Patient Registry (CFFPR) between January 1, 2003, and December 31, 2015. Patients younger than 6 years of age were excluded due to potentially unreliable pulmonary function testing. To focus the study on the 5 most modern era of CF care, we considered available data from 2003 onward. A detailed description of this registry and its contents is provided elsewhere (20).

Outcomes
To build the multivariate joint model for this cohort, we focused on established outcomes of lung function and nutrition. Outcomes on nutrition and growth, available from the CFFPR and computed based on formulas from the Centers for Disease Control and Prevention, included the aforementioned percentiles for BMIp, WFA and HFA. Observed FEV 1 was expressed as percentage of predicted using established reference equations for age, sex and race (21,22). Occurrence of first PE, the event outcome, was subject to censoring and considered to have occurred if documented in the CFFPR as warranting treatment with intravenous antibiotics in the hospital.

Covariates
Covariates from the aforementioned literature were considered, as well as findings from analyses specified in the most recent CFFPR report. These included genotype (F508del homozygous, heterozygous or neither/unknown), sex, Hispanic ethnicity, socioeconomic status (SES, having only state/federal or no insurance was recorded as 1, and 0 otherwise); time-varying covariates included use of pancreatic enzymes (corresponding to pancreatic insufficiency), CF-related diabetes (CFRD, with or without fasting hyperglycemia), positive cultures for methicillin-resistant staphylococcus aureus (MRSA) and Pseudomonas aeruginosa (Pa).. Baseline was defined as the time at which all longitudinal outcomes were first recorded during the eligibility period. A birth cohort variable was used to account for potential delayed entry into the registry and changes due to advancements in CF care and therapeutics approvals. Another timevarying covariate was used to account for irregular sampling bias, potentially induced by patients having a varying number of clinical encounters over time; for a given patient and encounter, this variable was the total number of encounters within the prior year. Data acquired after lung transplant were excluded.

Ethical aspects
The Institutional Review Board at Cincinnati Children's Hospital Medical Center approved the study.
Upon local IRB approval of this study, we submitted a data use proposal to the CFF Patient Registry Committee, which undertakes an established peer-review process of every request. We refined our request according to Committee feedback and, once they gave final approval, we received nonidentifiable and non-traceable data through an encrypted data delivery Citrix ShareFile system.

Descriptive analysis
Descriptive statistics, including median (Q1-Q3) for continuous variables and n (%) for categorical variables, were used to summarize cohort characteristics and extent of follow-up available for each patient.

Joint model setup and notation
Prior to multivariate joint modeling, we performed a univariate joint model that included only FEV 1 and time-to-first PE. To model the association between combinations of longitudinal FEV 1 and other outcomes together with time-to-first PE, we examined multivariate joint models. We describe the notation using a particular multivariate joint model with two longitudinal outcomes, which can be generalized to the case of three outcomes (e.g. FEV 1 , WFA and HFA). Assuming n subjects, we let y 1i (t) and y 2i (t) denote the follow-up measurements (e.g., FEV 1 and BMIp, respectively) for patient i (i = 1,…, n) at a specific time t. To describe the subject-specific evolution over time of the longitudinal outcomes, we rely on mixed-effects models. In particular we postulate where β1 and β2 denote the vector with the regression coefficients of the row vector of the design matrices for the fixed effects x1iTt and x2iTt and, z1iTt and z2iTt denote the row vectors of the design matrices for the random effects b1i and b2i. We included natural cubic splines to accommodate nonlinear progression of different covariate levels over age with respect to each longitudinal outcome, using an approach similar to our previous work in modeling nonlinear age-related FEV 1 progression 7 (5) (23). Furthermore, to build the correlation between FEV 1 and BMIp, we assume a multivariate normal distribution for the random effects bi = {b1i, b2i}.
For the PE part, we let Ti denote the observed failure time for patient i, taken as Ti = min(Ti*, Ci) with Ti* indicating the true failure time (occurrence of PE) that the ith individual experiences the event and Ci the censored time. To model the risks of the event, we postulate the proportional hazard model: where ωiT denotes row vectors of the design matrix of the baseline covariates, γ is the corresponding regression coefficients vector, and a1 and a2 are the coefficients that link the longitudinal and event parts. Particularly, they denote the strength of the association between FEV 1 and BMIp with onset of PE. To better understand the connection between the event and longitudinal parts in our model, we explain the meaning of the a1 and a2 parameters. Specifically, for a one-unit increase in the underlying value of FEV 1 for patient i, the hazard ratio is exp (a1) assuming that the baseline covariates and BMIp remain the same. Similarly, exp (a2) is the hazard ratio when the underlying value of BMIp for the i -th patient is increased by one unit assuming that the baseline covariates and FEV 1 remain the same.

Implementation
We employed Markov-Chain Monte-Carlo (MCMC) via Gibbs sampling using the 'JMbayes' package (Version 0.8-82) in R (Version 3.5.3., R Foundation for Statistical Computing, Vienna Austria) to obtain the parameters from the respective posterior distributions under the multivariate joint model (18). The highest posterior density (HPD) and accompanying standard errors were used to estimate each parameter of interest. Due to the large patient sample size that led to memory problems, we randomly divided the dataset into three parts (each part had similar percentage of PE at onset). Each of the datasets was fitted separately, and for the final results we combined the MCMC samples as described previously (24). Code implementation is provided as supplemental material. In total, we fit five joint models wherein PE onset was the event: i) only FEV 1 ; others included ii) BMIp; iii) HFA; iv) WFA; v) WFA and HFA. The multivariate joint models in (ii)-(iv) were implemented to evaluate the 8 impact of different measures of nutritional status on predicting PE onset. Due to the large number of parameters in scenario (v), we assumed a less flexible evolution for FEV 1 and BMIp over time. In particular, we postulate natural quadratic splines to estimate this evolution instead of natural cubic splines assumed in the other scenarios.

Evaluating predictive performance
We calculated the area under the receiver-operator characteristic curve (AUC) to evaluate predictive accuracy of each model with respect to PE events using ten-fold cross validation. A large value of AUC indicates the preferred model. To obtain correct estimates of the AUC we performed a 5-fold cross validation procedure, wherein each time we selected 700 individuals from the large data set and we split it in five subsets. Each time we fitted the model in four of the subsets, we calculated the AUC in the subset that was excluded. The calculation of the AUC was performed at 12 and 16 years of age with prediction windows of 0.5, 1, 2 years, in order to mimic clinically meaningful age strata and prediction windows. This procedure was repeated 100 times. Appropriate diagnostics for joint models were examined. The code of the analysis and cross validation are shown in the Supplement.

Results
Analysis cohort There were 16,455 patients contributing 245,513 observations to the analysis cohort (Table 1).
Patients who had a PE during follow-up tended to be female and had the F508del homozygous genotype. They had slightly higher proportions of infections with MRSA and Pa than patients who did not have a PE during follow up. By contrast, patients who were PE-free during follow-up had a higher reported use of pancreatic enzymes and had higher levels of lung function and better nutrition status at baseline. Extent of follow-up also differed according to PE onset. Representative trends of lung function decline according to PE status are shown in Figure 2 and reflect summary results from Table   1. Median (95% CI) age of PE onset was 19.5 (19.3-19.6) years ( Figure 2).

Joint model estimates
Each Gibbs sampler ran for 900 iterations, and the first 600 iterations were discarded for burn-in.
Shared parameter estimates, which represent associations between each longitudinal process and the 9 PE event, are shown for all five fitted joint models ( Table 2). The log-hazard estimates of the shared parameters, which were introduced in Section 3.2, indicate that FEV 1 was consistently and negatively associated with PE onset. Corresponding hazard ratio estimates ranged from 0.96 to 0.97, suggesting a small association between changes in the FEV 1 trajectory and PE onset. Jointly modeling BMIp, WFA or HFA as shown in respective models (ii)-(iv) did not impact association between FEV 1 and PE onset, and their associations with PE onset did not reach statistical significance. The univariate joint model (i) in Table S1 indicated that being Hispanic and having lower SES, MRSA and CFRD and using pancreatic enzymes corresponded to worse overall FEV 1   The multivariate joint model (ii) in Table S2, which included a longitudinal submodel for BMIp, showed that being born into an older birth cohort, taking pancreatic enzymes and having more frequent clinic visits corresponded to higher BMIp. There were nonlinear associations between each of genotype, SES, birth cohort and Pa with BMIp trajectory. In Table S3, using height in multivariate joint model (iii), we found that having a F508del homozygous genotype and belonging to an older birth cohort were positively associated with overall HFA. By contrast to the BMIp submodel, having more frequent clinic visits corresponded to lower HFA. Similar nonlinear associations were observed in the BMIp and HFA submodels. Results for the WFA submodel in (iv) were also similar ( Table S4). The multivariate joint model in (v), compared to models (iii) and (iv), had slight changes in associations shown in the WFA submodel; however, the HFA submodels were similar.

Predictive performance
The results of the predictive performance of the joint models presented in Table 2 are summarized in Figure 3. In particular, the AUC of each model is presented assuming a different prediction interval (0.5, 1 and 2 years), which starts at adolescence or early adulthood (ages 12 and 16 years, respectively). The median AUC values ranged from 0.7 to 0.8. Moreover, the univariate joint model, assuming only FEV 1 ,provided the highest AUC value and was most robust across age strata and prediction windows. Precision of the predictions was assessed by examining the interquartile ranges of the boxplots in Figure 3. Precision was lower for predictions done at adolescence, regardless of the interval, compared to precision estimated during young adulthood. Precision was highest under the univariate joint model for predicting PE risk at age 16 years out to 2 years. In terms of overall precision, the multivariate joint models performed similarly well.
Dynamic predictions for an individual CF F508del homozygous patient (born 1990-1994) are shown in Figure 4 for the multivariate joint model of longitudinal FEV 1 and BMIp. Her first recorded FEV 1 and BMIp were 98.5% predicted and at the 48.2 percentile, respectively. She began taking pancreatic enzymes and had a positive culture for Pa around 10 years old. Her monitored outcomes are depicted across four clinic visits from ages 8.2 to 14.3 years, along with probability of not having a PE over follow up. Her long-term risk of PE appears higher during an earlier visit (e.g., PE-free probability is 0.6), compared to later visits that are informed by dynamic predictions from the multivariate joint models (e.g., PE-free probability increased to 0.8). Commensurate with these results are the changes in her FEV 1 and BMIp trajectories, which imply minimal rate of decline and improving nutritional status, respectively. She did not experience a PE event during follow up. Her projected PE-free probabilities over follow up using a multivariate joint model replacing BMIp with WFA and HFA were similar ( Figure 5); however, precision with which PE-free probability could be estimated was decreased using this combination of markers rather than BMIp.

Discussion
Past epidemiologic studies of PE onset have largely focused on explanation of risk factors rather than identifying predictors of this event. Clinical consensus on PE symptoms is lacking. Acute drops in lung function and changes in nutritional status have been recognized as key determinants (25); however, relative contributions of these markers to prediction of PE onset have not been studied. Further, markers of growth/nutritional status (WFA, HFA and BMIp) and lung function (FEV 1 ) are endogenous variables. In this epidemiologic study, we assessed predictive performance of these markers using a series of joint models. Our novel application confirmed the reliance of changes in the FEV 1 trajectory on predicting PE events but revealed the clinical utility that these novel models could have for monitoring multiple indicators of disease progression for an individual patient.
We found that the univariate joint model, which utilized FEV 1 and PE event submodels without accounting for growth/nutritional status, produced the most accurate and precise predictions of all models considered. In addition, these predictions were largely robust across a range of age and prediction window strata. A potential reason for this finding is that the diagnosis of PE events may rely more heavily upon acute changes in lung function than changes in nutritional status. Changes in growth are relatively slow over time, compared to precipitous drops in lung function, which could be responsible for the lack of improved predictive performance when adding WFA, HFA or BMIp. Of the three markers available from the CFFPR to study growth/nutritional status, it is plausible that BMIp would be most likely to improve prediction of PE onset (after including FEV 1 ); however, it appears that including BMIp as a submodel introduces additional noise into the joint modeling of these processes, thereby making predictions less accurate and more variable (Figure 3). In addition, associations between clinical/demographic covariates and each outcome were largely consistent across FEV 1 submodels. This finding further suggests that WFA, HFA and BMIp do not play a large role in clinical diagnosis of PE.
The multivariate joint models afforded the opportunity to obtain dynamic predictions at the individual patient level ( Figure 4); however, using multiple nutrition markers may decrease precision of PE predictions ( Figure 5). Using these dynamic predictions revealed the heterogeneous nature of CF disease progression but also highlighted the benefit of updated predictions across clinic visits. The ability to simultaneously monitor multiple markers of disease progression, which may serve as a clinical decision support tool, provides an advantage of the multivariate joint model over the univariate joint model and the more conventional Cox model. A limitation of this study is that the models cannot be used to draw conclusions about recurrent PE events, since the PE event submodel estimates time to first PE. Although this focus simplifies analyses, time-to-first PE during follow up tended to occur in early adulthood (Figure 2). It is possible that a portion of patients experienced a PE prior to follow-up available in the data (e.g., patients born in older birth cohorts). In addition, the advent of newer therapies and/or quality improvement initiatives likely increase the overall trajectories of outcomes studied. We accounted for these potential sources of bias by including levels of birth cohort through linear and nonlinear terms in each of the longitudinal submodels (Tables S1-S5). Consequently, including higher-order associations through splines required additional degrees of freedom. Future work is needed on covariate selection strategies in the context of multivariate joint models (26) (27). Fitting these more complex models to a large number of patients followed over long periods of time creates a computational challenge. We addressed this issue by partitioning into distinct sets of patients, which was the only feasible manner to estimate the model parameters; however, the bootstrap approach that we employed does not account for all the variability of the outcomes. Lastly, the multivariate joint models presented here allow for different patterns of missing data among the longitudinal outcomes but assume that the missingness is conditional on observed data only (i.e., missing at random), specifically occurrence of PE. It is possible that missingness could be related to patient severity. A recent study demonstrates that linear mixed effects modeling, which is the framework for the longitudinal submodels used in this application, can alleviate this potential source of bias (28).

Conclusion
The univariate joint modeling of the CFFPR demonstrates the reliance of PE diagnosis on FEV 1 , and multivariate joint modeling of FEV 1 and other outcomes related to nutritional status and growth can be used for routine medical monitoring of an individual patient and dynamic assessment of PE risk over time. All individuals diagnosed with CF and associated disorders who receive care at a center accredited by the CFF who have provided informed consent are considered eligible to have their data collected in the CFFPR. Per CFF requirements to acquire data in the CFFPR, each center obtained approval from their local institutional review board and written informed consent and assent from individuals and respective legal guardians as required. The analysis cohort derived from the CFFPR and utilized in this particular study was therefore retrospective in nature and de-identified prior to receiving data from the CFF.

Consent for publication
Not applicable.      Multivariate joint model predictions with other growth markers. Same female CF patient as in Figure 4, but simultaneously modeling lung function (FEV1, upper-left y-axis), growth as weight-for-age (WFA, middle-left y-axis) height-for-age (HFA, lower-left y-axis) percentiles and probability of not having a pulmonary exacerbation (PE-free probability, repeated along the right y-axis). Her probability of remaining PE-free is similar to Figure 4 but estimates are less precise based on 95% CI over the four clinical encounters.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.