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 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 FEV1 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 FEV1 trajectory and PE onset. Jointly modeling BMIp, WFA or HFA as shown in respective models (ii) – (iv) did not impact association between FEV1 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 FEV1. There were nonlinear associations, as reflected by spline coefficients, between each of genotype, sex, SES, MRSA, Pa and CFRD and FEV1 decline. In the PE event submodel, patients who were male, used pancreatic enzymes, had lower SES, infection with MRSA and Pa and were diagnosed with CFRD had greater risk of PE onset. The multivariate joint models, each of which included a longitudinal submodel for FEV1, had similar association results for the clinical/demographic covariates and shared parameters regarding FEV1.
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.
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 FEV1, 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 FEV1 and BMIp. Her first recorded FEV1 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 FEV1 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.