Study design and cohort

We conducted a retrospective longitudinal cohort study using the CFFPR, which is a national patient registry that collects demographic and clinical data on individuals with CF who are patients at care centers across the U.S [14]. The study included data from 2003 to 2018. To construct the analysis cohort, we considered patients with a valid CF diagnosis (e.g., blood test, sweat test, genetic test) who carry the G551D mutation and received ivacaftor any time after January 1, 2012, which represents the era of widespread U.S. Food and Drug Administration (FDA) approval for this therapy. Data observed when patients were younger than 6 years old was excluded, given the potential for unreliable pulmonary function test (PFT) results in very young patients with CF. Data observed after lung transplantation were censored.

Determining ivacaftor initiation. We omitted PFTs observed during the first 30 days after the CFFPR-recorded ivacaftor start date (Fig. 1). Our rationale was (i) we wanted to avoid estimating the initial increase in FEV1pp related to ivacaftor that has been previously reported in other analyses [2], and (ii) we sought to reduce the potential bias between the ivacaftor prescription date and the actual start date recorded in the CFFPR. To ensure valid estimation before and after ivacaftor initiation, we restricted the cohort to patients who had at least one PFT before initiation and at least two PFTs after, and for whom the earliest and latest post-initiation PFTs were separated by at least 6 months in time.

Fig. 1

Lung function trajectory and ivacaftor initiation shown for a male with cystic fibrosis in the analysis cohort. The outcome was measured as the percent predicted forced expiratory volume in 1 s (FEV1pp, y-axis) observed over time (x-axis). (A) before and (B) after ivacaftor initiation. His pre-ivacaftor FEV1 ranged from 68 to 118% predicted, while post-ivacaftor was 82 to 102% predicted. His overall FEV1pp increased, but he experienced FIES events (examples of events shown using red arrows in A and B) before and after ivacaftor initiation

FIES definition. We derived our criteria at each clinical encounter from the definition provided by the CF Learning Network [12], which was applied as follows:

1)

For baseline FEV1pp ≥ 50, if the current FEV1pp represents a 10% or more relative decline in lung function, compared to the baseline.

2)

For baseline FEV1pp < 50, if the current FEV1pp represents a 5% or more relative decline in lung function, compared to the baseline.

In this definition, baseline is the average of the two highest FEV1pp values in the past 12 months that were not recorded during intravenous antibiotic treatment. Expanded details on the FIES definition are provided as supplemental material (Part 1, Section I).

Statistical analysis

Outcome, covariates, and missing data. The study outcome was the FEV1pp value observed at each clinical encounter. Covariates included observed demographic and clinical characteristics that have been previously associated with accelerated FEV1pp decline: the time-varying variables were Medicaid insurance use, infection with methicillin-resistant Staphylococcus aureus (MRSA), infection with Pseudomonas aeruginosa (Pa), diagnosis of CF-related diabetes (CFRD), and numbers of acute pulmonary exacerbations and outpatient visits within the previous year; the non-time-varying variables were age and FEV1pp at study entry, birth cohort (defined based on year of birth), sex, and pancreatic insufficiency (defined as any reported use of pancreatic enzymes). All subsequently described prediction modeling assumed that outcome data were missing at random [15].

Prediction model setup. A previously described longitudinal model framework with nonstationary stochastic process to the prediction of FEV1pp and the use of clinically relevant target functions in CF was adapted for this study [10]. Particularly, the variance terms in the linear mixed effects model included a random intercept to account for between-patient variation, an integrated Brownian motion covariance function to account for within-patient variation over time and to allow us to create predictive probability distributions using a prior approach, and a residual measurement error [16]. We first set up a saturated model within this framework to examine various covariate effects. The time since study entry (in years) was used as the time variable. A change point term, which represented pre- and post-ivacaftor initiation periods, was included as a main effect, and its interaction with time was used to examine associations between ivacaftor response and absolute FEV1pp, as well as the difference in slopes between pre- and post-ivacaftor initiation periods. Covariates were considered as both main effects and interactions with the time variable in the saturated model. Results were scaled to the time since ivacaftor initiation (in years) for presentation purposes. Details of the model setup are presented as supplemental material (Part 2, Section I) and the residual diagnostics of the selected model are shown in the supplemental material (Part 1, Section II).

Covariate selection. Reduced forms of the previously described terms in the saturated model were examined with the Akaike and Bayesian information criteria (AIC and BIC, respectively) and the likelihood ratio test (LRT).

Predictive probabilities for FIES events. The target function derived from the above FIES definition was implemented as part of the model fitting in the R package “lmenssp” version 1.2 [17]. The formulas and code can be found from supplementary material Part 2, Sections II-III and Part I, Section V, respectively.

Validation. Two types of validation were performed that are relevant to clinical scenarios: (i) predictions for “new patients,” and (ii) forecasting for patients who were part of the model building but return for follow-up visits (i.e., updated predictions). To accomplish both types of predictions, the analysis cohort was randomly split into 80% for training and 20% for independent testing (Fig. 2). For type (i), we examined predictions within the testing cohort. For type (ii), we examined data that were held out for the last 6 months of follow-up in the training cohort. Both types of predictions were evaluated using 5-fold cross-validation.

Fig. 2

Dataflow for model fitting and validation. The overarching analysis cohort (the first level) was segmented into training and testing cohorts (the second level). The training cohort was further split into a cohort for model fitting and a masked cohort for forecast validation (the third level)

Evaluation of predictive performance. Metrics to evaluate predictive accuracy included the root-mean-square error (RMSE), the mean absolute error (MAE), the Brier score, and the area under the receiver operating characteristic curve (AUC). Formulas for the metrics are provided as supplemental material (Part 1, Section III). Lower values of the RMSE, MAE, and Brier score imply higher predictive accuracy, while lower AUC values indicate lower predictive accuracy. The 95% confidence interval (CI) for each AUC estimate was obtained through nonparametric bootstrapping with 1,000 replicates via the R package “boot” version 1.3–28.1 [18].