Longitudinal FEV₁ and Exacerbation Risk in COPD: Quantifying the Association Using Joint Modelling

Introduction

Exacerbations of chronic obstructive pulmonary disease (COPD) are episodes of respiratory symptom worsening requiring additional therapy (eg oral steroids, antibiotics) and/or hospitalization. Prevention of exacerbations is an important goal in COPD treatment and the exacerbation event is recommended by regulatory agencies as an endpoint for assessing efficacy in clinical studies.¹ Effects of drug treatment on exacerbation risk are commonly analyzed based on event rates, eg using negative binominal regression, or time-to-first-event, eg using Cox regression. The relatively low frequency of exacerbations in COPD patients means that long and large studies are needed to get desirable precision in treatment effect estimates. Consequently, exacerbations are usually not studied until Phase III, while there is a need to predict efficacy and dosing regimen of novel treatments in earlier phases of development.

Lung function, measured as the change from baseline in forced expiratory volume in one second (FEV₁), is another frequently used primary endpoint to assess treatment effects in COPD in Phase III, as well as Phase II dose-finding trials of new treatments. FEV₁ is usually measured at baseline and then repeatedly during a study (ie longitudinally), and through analysis of repeated measurement data, using linear or nonlinear mixed-effects (LME, NLME) models, increased precision in treatment effects can be obtained.

It is well established that there is correlation between FEV₁ and exacerbation risk; an improved lung function is related to lower risk of exacerbation.^2–8 The association between improvements in FEV₁ and reduced exacerbation risk was quantified using meta-analysis approaches,^4,5 and based on individual patient data.^6–9 These results indicate that there may be value in co-analyzing longitudinal FEV₁ and exacerbations in clinical trials; ie capturing treatment effects on both endpoints, and the inter-dependencies between endpoints, in the same model could potentially increase statistical power, and enable assessment of exacerbation efficacy in shorter and smaller trials.

The concept of joint modelling of longitudinal biomarker data with time-to-event data provides an approach to adequately handle the influence of an endogenous time-dependent covariate, as it allows the longitudinal biomarker to affect the hazard in a time-dependent manner, while accounting for the measurement error, and can provide more efficient and less biased estimates of treatment effects. This method has received increasing attention, eg in oncology, focusing on association between quality of life¹⁰ or tumor size dynamics¹¹ and survival, and treatment of HIV linking longitudinal CD4+ count and survival (eg¹²). To our knowledge, this method has not been applied in COPD.

The aim of the present work was to evaluate joint modelling as a tool for analyzing treatment effects in COPD clinical trials by (1) quantifying the association between longitudinal improvements in FEV₁ and the risk of exacerbation, (2) comparing a joint model of the two endpoints to a Cox proportional hazards model of only exacerbations, and (3) assessing the consistency in parameter estimates across several clinical studies and treatments. To this end, we used a large set of patient-level data from three Phase III clinical studies in moderate-to-severe COPD, evaluating the efficacy of fixed-dose combinations of a bronchodilator (long-acting beta-agonist, LABA) and an anti-inflammatory (inhaled corticosteroid, ICS), two compound classes included in the standard of care treatment for COPD.

Methods

Data Description

Data from three clinical COPD studies, evaluating fixed-dose combinations of formoterol (LABA) and budesonide (ICS) versus the mono-components and placebo, were used in the analysis.^13–15 A subset of the original study data was available for analysis, including patients who had provided informed consent for data re-use, as summarized in Table 1. Full details of the original studies have been published elsewhere.^13–15 Data from the largest study, Study A (NCT00206167, 1740 patients),¹³ was used for model development and qualification. Two additional studies: Study B (NCT00206154; 1604 patients)¹⁴ and Study C (NCT00419744; 1042 patients)¹⁵ were used for external model validation and subsequent re-estimation of model parameters. All three studies were performed across multiple geographical regions and the only country represented in all three studies was the USA. Summaries of baseline characteristics, per study, are shown in Table 2. The studies were selected to provide adequate numbers of exacerbations and included longitudinally measured FEV₁, to enable appropriate analysis of the association between these clinical endpoints.

Table 1 Study Data Details: Duration, Treatment Arms, FEV1 Measurement Schedule

Table 2 Summary of Baseline Characteristics, by Study

Pre-dose FEV₁ was measured repeatedly over time, and the time to first exacerbation was captured in the studies. Exacerbations were defined as in the original study protocols, ie hospitalization and/or oral steroid treatment due to worsening of COPD. The percentage of patients with at least one exacerbation during the study was 34%, 25% and 44% for study A, B and C, respectively.

Analysis

Cox Proportional Hazards Model for Time-to-First Exacerbation

A conventional Cox proportional hazards model for time-to-first exacerbation was used as a reference model to estimate exacerbation treatment effects in Study A. Categorical treatment arm, country and baseline FEV₁ were included as covariates.

A Joint Model of FEV₁ and Exacerbations

A joint model (referred to as the base model) of exacerbation risk and longitudinal pre-dose FEV₁ was constructed. The model consisted of two sub-models, the Cox proportional hazards model for time-to-first exacerbation (above) and an LME model for longitudinal FEV₁, linked with an association parameter describing how the estimated individual response in FEV₁ affects the exacerbation hazard (Figure 1). For the Cox proportional hazards sub-model, the baseline hazard was defined as a piecewise constant function.

Figure 1 Overview of the joint models of individual longitudinal FEV1 and exacerbation hazard, indicating the covariates, and longitudinal data, included in the base model and the model for prediction of FEV1 and exacerbation hazard. Prognostic factors refer to additional baseline covariates identified in the covariate search.

The LME sub-model used longitudinal change from baseline FEV₁ (∆FEV₁) as the independent variable and was built using natural splines. Different placements of the spline knots were investigated, for both the fixed and random effects of the model. Treatment was initially included as a categorical covariate, with one fixed-effect spline estimated for each treatment arm.

This base joint model differed from the reference Cox model (above) only in that the effect of longitudinal ∆FEV₁ on the exacerbation hazard was included. Parameter estimates of the model, when applied to Study A, were compared to those of the reference Cox model. For additional details, see the Supplementary Material.

Updated Joint Model for Prediction

After evaluation of the base joint model, it was updated to allow for prediction of outcome in other studies with budesonide and/or formoterol treatments (Figure 1). The categorical treatment arm covariate was removed from the exacerbation hazard and was substituted by separate components to reflect the bronchodilator and anti-inflammatory effects.

Formoterol (bronchodilator) was used as a covariate in the longitudinal FEV₁ model and was thereby assumed to only affect the exacerbation hazard via effects on FEV₁. Absolute FEV₁ was modelled as the independent variable in this model and thus replaced baseline and ∆FEV₁. Budesonide (anti-inflammatory) was used as a dose-dependent covariate in the Cox proportional hazards sub-model.

An extensive covariate search was further performed. The following baseline covariates were tested on the exacerbation hazard: country, exacerbation history, breathlessness score, sleep score, sputum score, cough score, eosinophils, sex, age, race, FVC, FEV₁/FVC and season when the treatment started. Up to five covariate combinations were tested. The final set of covariates was selected based on the Bayesian Information Criterion (BIC).

Model Validation and Parameter Re-Estimation

The updated model was validated by predicting the outcome of Studies B and C, using a 1-month cutoff for longitudinal FEV₁ data, thus including at least one FEV₁ measurement post-baseline (see the Supplementary Material for details). Only patients from the USA were included in the validation procedure, since it was the only country represented in all 3 studies (see Table 2 for details). The model was also re-estimated on Studies B and C to investigate the consistency of parameter estimates across studies.

Lastly, the consistency of the association parameter linking FEV₁ to the exacerbation hazard across treatment mechanisms was assessed. This was done by re-estimating the association for each treatment arm (vs reference), and for the reference arm alone, in each study separately (for details see the Supplementary Material).

Software

All models were estimated using the JM package version 1.4–7¹⁶ in R version 3.2.4.¹⁷ The joint models were estimated and validated using a joint conditional formulation of the likelihood function, incorporating the impact of FEV₁ and hazard function, as described in.¹⁶

Results

Comparison of the Cox Proportional Hazards Model and the Base Joint Model

A significant association between longitudinal ∆FEV₁ and the risk of exacerbation was estimated (p<0.0001) in the joint model for Study A (Table 3). The estimated association constant for ∆FEV₁ was of similar magnitude to that of the baseline FEV₁ covariate effect. The hazard ratios (HRs) for each of the treatment arm comparisons estimated with the Cox proportional hazards model and the base joint model are also shown in Table 3 (full model outputs can be found in Supplementary Material - Tables S1 and S2). The estimated HRs are higher for the joint model vs the Cox model (by approximately 5 percentage points) since part of the treatment effect is captured via the longitudinal ∆FEV₁ model. The effect of ∆FEV₁, however, is small in relation to the overall treatment exacerbation risk reduction (30–35% for treatment arms including budesonide). No improvement in precision was seen in the estimated total treatment effect of the joint model (ie the combined effect of longitudinal FEV₁, via the association, and the treatment coefficients in the Cox model) (data not shown). For details on the longitudinal sub-model see the Supplementary Material.

Table 3 Estimated Hazard Ratiosa for Exacerbation Risk (Cox Model and Base Joint Model, Study A)

The Updated Joint Model for Prediction

The updated joint model for prediction (Figure 1), where the categorical treatment effect was substituted by separate components for the bronchodilator (included in the longitudinal sub-model) and anti-inflammatory effects (included in the hazard), resulted in an adequate fit to the data. The estimated HR per 160 µg budesonide was 0.85 (95% CI: 0.76–0.93). The covariate search identified the following baseline covariates: country, breathlessness score (HR 0.84 (95% CI: 0.74–0.95) per decrease by 1 unit) and the number of exacerbations in the previous year (HR 0.88 (95% CI: 0.83–0.94) per decrease by 1 unit).

The variability in individual FEV₁ profiles was large, as can be seen in Figure 2A. However, the longitudinal behavior of FEV₁ was adequately described (Figure 2A and B). The model also adequately described time-to-event data for all four treatment arms in Study A (Figure 2C). Model parameter estimates and additional goodness-of-fit figures can be found in the Supplementary Material (Table S3 and Figure S1).

Figure 2 Diagnostics of the joint model for prediction fitted to Study A. Mean of individual predictions of FEV1 (A) and change from baseline FEV1 (B) with 95% CI (solid colored lines and shaded areas), superimposed with individual profiles (thin colored lines) and means of the original data (stars) vs time. (C) Predicted mean exacerbation-free probabilities (solid line) with interquartile range (shaded area) and Kaplan–Meier estimates (dashed line).

Model Validation and Parameter Re-Estimation

The updated joint model developed based on Study A was assessed in terms of its predictive ability for Studies B and C. As shown in Figure 3, the model successfully predicted 6- and 12-month exacerbation outcomes (for patients from the USA) based on a data cutoff at the second visit (1 month), including the two additional treatment arms in Study B which were not present in Study A: budesonide 320 μg and the combination of budesonide 320 μg and formoterol 9 μg using separate inhalation devices (Figure 3A).

Figure 3 Prediction of exacerbation outcomes, for patients in the USA, in Study B (A) and Study C (B), using the joint model for prediction developed on Study A, and a 1-month data cut-off. Predicted exacerbation-free probabilities (means and 95% CI, blue) vs Kaplan–Meier estimates (dashed line). The shaded area in (A) denotes additional treatment arms in Study B which was not included in Study A.

When re-estimating model parameters on Studies B and C, the current longitudinal value of FEV₁ was found to be significantly associated with the risk of exacerbation also in these studies (p<0.0001). Moreover, parameter estimates were consistent across all three studies; an increase of 100 mL in FEV₁ decreased the risk of exacerbation by 8–9% (Figure 4). Full model outputs can be found in the Supplementary Material (Tables S4 and S5, Figures S2 and S3). Parameter estimates were also similar across treatment arms (Figure 5).

Figure 4 Parameter estimates of the joint model for prediction with 95% CI per study. Histograms show the instantaneous exacerbation risk change with respect to the relative parameter change.

Figure 5 Estimated hazard ratio per 100 mL improvement in FEV1 per study arm for each study; for the reference arm alone (filled circles) and each treatment arm vs reference (open circles). Error bars represent 95% CI. Abbreviations: Form, formoterol; Bud, budesonide.

In addition, a consistent dose-dependent effect of budesonide was estimated; the exacerbation risk decreased by 10–15% per 160 µg of budesonide. A patient’s baseline breathlessness score and previous history of exacerbations influenced the risk of exacerbation, even though the point estimates of the effects varied approximately 2-fold across studies, as shown in Figure 4.

Discussion

The primary aim of this work was to evaluate joint modelling as a tool for analyzing treatment effects in COPD clinical trials. This was achieved by applying joint models to a large set of patient-level data on longitudinal FEV₁ and exacerbations, two of the most established endpoints of drug effect in COPD, measured in three phase III studies of budesonide/formoterol. The developed joint models, which included an LME model for longitudinal pre-dose FEV₁ and a Cox proportional hazards model for time-to-first exacerbation, estimated statistically significant association between the two endpoints (p<0.0001). Notably, the estimate of the association parameter was consistent across studies and treatments, including placebo, confirming that longitudinal FEV₁ contains information about exacerbation risk.

The size of the association parameter implied an approximate 10% reduction in instantaneous exacerbation risk for a 100 mL improvement in pre-dose FEV₁. When considering the limited treatment effects on FEV₁ in the COPD studies used in our analysis, on average around 50–90 mL, the average exacerbation risk reduction associated with improvements in FEV₁ is in the range of 4–7%. The magnitude of the effect is thus relatively small and only a minor part of the exacerbation risk reduction is accounted for by longitudinal changes in FEV₁. This can be partly attributed to variability in FEV₁ measurements,⁶ however, it also emphasizes that other factors than lung function improvement are important for exacerbation risk reduction. This is particularly true for anti-inflammatory drugs like ICS, which have limited effects on FEV₁ in COPD patients but significantly reduces exacerbation risk.

We found the association between lung function and exacerbation risk to be consistent across three studies of an ICS/LABA combination treatment. Furthermore, we found the FEV₁-exacerbation association parameter to be similar across all treatment arms (budesonide, formoterol, budesonide/formoterol and placebo) when estimated separately. Although the studies used in our analysis are of similar design, in similar populations, this suggests that the effect on exacerbation risk, for a specific change in FEV₁, is unrelated to the mechanism of action. This is also supported by the meta-analyses done by Zider et al⁴ and Ribbing et al,⁵ who identified similar slopes of the relationship for anti-inflammatory and bronchodilator compounds.

The fact that the FEV₁-exacerbation association is of a similar magnitude to that previously reported based on other methods^7,8 indicates that the joint modelling methodology works and has the potential to be used also with other variables and potentially extended to multivariate models (eg¹⁸). It also suggests a strong prior probability can be used for the association between FEV₁ and exacerbation risk when applying this type of joint model in a Bayesian setting.¹⁹ However, given the modest effect of longitudinal FEV₁ on exacerbation risk reduction and large variability in FEV₁ data, we see no immediate impact on trial design based on these results. In our updated joint model for prediction, we estimated a significant dose-dependent effect of budesonide (in addition to the longitudinal lung function effect) in the exacerbation hazard model. Ideally, this should be replaced by (longitudinal) biomarkers capturing anti-inflammatory effects. Breathlessness score, at baseline, was found to be a covariate strongly affecting the risk of exacerbations. There may therefore be value in accounting for longitudinally measured breathlessness score in a multivariate joint model. Inclusion of multiple longitudinal variables, such as symptom scores, quality of life assessments (patient-reported outcomes), and biomarkers related to inflammation has the potential to maximize the information from the data and improve statistical inference and predictive ability.

An additional consideration is the link between treatment response in FEV₁ and risk of early drop-out. It has been shown that treatment failure of a bronchodilator (ie lower FEV₁ response) increased the risk of early drop-out.²⁰ Not accounting for such informative drop-out could lead to both bias and imprecision of parameters,^21–23 hence extending the modelling to competing events could be beneficial.

Conclusions

Joint modelling can be used to co-analyze longitudinal FEV₁ and exacerbation data in COPD clinical trials. An approximate 10% reduction in exacerbation risk was estimated per 100 mL improvement in FEV₁, consistent across three trials and multiple treatment arms, confirming that treatment effects on FEV₁ contain information about exacerbation risk reduction. However, due to the relatively small contribution of improved lung function to the overall risk reduction, no major impact on exacerbation trial design can be expected based on FEV₁ alone. Further exploration with other longitudinal endpoints should be considered to evaluate the use of joint modelling in analyzing and predicting outcome of COPD clinical trials.