Study population
The ETHOS study (NCT02465567) evaluated the efficacy and safety of the inhaled corticosteroid (ICS)/long-acting muscarinic antagonist/long-acting β2-agonist fixed-dose combination budesonide/glycopyrronium/formoterol fumarate dihydrate metered dose inhaler (BGF MDI; referred to hereafter as BGF) at two budesonide dose levels. Patients with moderate-to-very severe COPD were randomised 1:1:1:1 to BGF 320/18/9.6 μg, BGF 160/18/9.6 μg, glycopyrronium/formoterol fumarate dihydrate MDI (GFF) 18/9.6 μg or budesonide/formoterol fumarate dihydrate MDI (BFF) 320/9.6 μg [25]. All treatments were administered over 52 weeks as two actuations, twice-daily, using an Aerosphere™ inhaler. Triple therapy with BGF 320/18/9.6 μg was found to reduce the risk of moderate or severe exacerbations vs dual therapy with BFF and GFF, and mortality vs GFF [25]. All methods were carried out in accordance with relevant guidelines and regulations, see the ‘Ethics approval and consent to participate’ section.
Inclusion/exclusion criteria
Inclusion and exclusion criteria have been published previously [25]. Briefly, eligible patients were 40 to 80 years of age and had symptomatic COPD (defined as a score of ≥ 10 on the COPD Assessment Test) with a post-bronchodilator ratio of FEV1/forced vital capacity < 0.7, with a post-bronchodilator FEV1 of 25% to 65% predicted; and had a documented history of ≥ 1 moderate or severe COPD exacerbations (if their FEV1 was < 50% predicted) or ≥ 2 moderate or ≥ 1 severe COPD exacerbations (if their FEV1 was ≥ 50% predicted) in the year prior to screening. Patients were excluded from covariate-adjusted analyses if they had any missing baseline covariate data.
Statistical analysis
This post hoc analysis of the ETHOS study included on-treatment events from the modified intent-to-treat population (mITT; all patients who underwent randomisation, received a study treatment and had post-randomisation data obtained before discontinuation of treatment).
As an initial investigation into the association between exacerbations and mortality, Kaplan–Meier estimates for the proportion of patients dying following exacerbation events were obtained. Survival times were censored at the time of transition to another exacerbation event or time of treatment discontinuation. A hazard ratio (HR) for the risk of death following a severe exacerbation was estimated using a Cox proportional hazards model with an indicator of a severe exacerbation as a time-dependent covariate. Initial analyses made no adjustment for baseline patient characteristics.
For each patient, times of exacerbation events and death were observed exactly or right-censored (i.e., the data points were known to be above a certain value, but it was unknown by how much). Exacerbations were considered distinct events if there were > 7 days between the end of one event and the start of the next. This produced multi-state survival data, where there is a series of event times for each patient, corresponding to times of transition to the next state. At any point in time, a patient occupied one of five well-defined states based upon the events cumulatively experienced while on treatment in the ETHOS study:
-
State 1: No exacerbation since study start
-
State 2: 1 moderate exacerbation (without severe exacerbation) since study start
-
State 3: ≥ 2 moderate exacerbations (without severe exacerbation) since study start
-
State 4: ≥ 1 severe exacerbations since study start
-
State 5: Death
The last observed state of a patient was either death or censored.
To formally characterise the association between exacerbations and mortality, multi-state modelling based on the Markov assumption (where transition to future states depends only on the current state at time t and not previously occupied states) was used [24, 26]. Semi-parametric and fully parametric multi-state models were employed. Each multi-state model estimated transition intensities or transition hazards, representing the instantaneous risk of a patient moving from one state to another in continuous time [24, 27, 28].
All patients were in the ‘no exacerbation’ state at randomisation. From this state, the model only allowed a patient to transition to a worsened health state. Transitions did not necessarily have to be to the next adjacent state (i.e. an individual could transition from no exacerbations to death without an exacerbation); however, transitions directly from no exacerbations to ≥ 2 moderate exacerbations were not permitted (Fig. 1). Additional details can be found in the supplementary material.
For the semi-parametric model, a standard Cox proportional hazards regression was fitted to the multi-state data. This specified a different baseline hazard for each transition by including each transition type as separate strata. The baseline hazard was estimated non-parametrically, using the Breslow estimator.
Standard parametric models (exponential, Weibull, Gompertz, log-logistic, log-normal, gamma and generalised gamma) were fitted to the multi-state data and compared using visual assessments and Akaike information criterion (AIC). Times from baseline to transition to next state were modelled, assuming a distinct rate parameter for each transition type. Visual assessments of model fit included comparisons of the Aalen-Johansen estimates [29] of state occupancy probabilities and the Nelson-Aalen estimates [30] of the cumulative hazard (i.e. non-parametric estimates) with estimates based on the fitted parametric models. For each parametric model, state occupancy probabilities were obtained and plotted for various time points up to 5 years, which was beyond the typical duration of an outcome trial in COPD [31,32,33]. The best fitting parametric model was that which fitted well to the observed data and gave clinically plausible state occupancy probability estimates up to 5 years.
To assess whether risk of death increases as patients accumulate exacerbations, estimated regression coefficients and corresponding standard errors from the best fitting parametric model were used to obtain hazard ratios with 95% confidence intervals (CIs), comparing the risk of death following 1 moderate exacerbation, ≥ 2 moderate exacerbations and ≥ 1 severe exacerbations with the risk of death following no exacerbation. Similarly, hazard ratios were obtained to compare the risk of ≥ 1 severe exacerbations from states of exacerbation with the risk from no exacerbation.
The transition probabilities at 4 weeks from study entry were estimated based on the best fitting parametric multi-state model, without adjustment for covariates. To estimate probabilities of occupying each state at a fixed point in time t, we computed a transition probability matrix, where the (r,s) entry is the probability of occupying state s at time t, given the state at time 0 is r. Since we used Markov models, this was obtained by solving the Kolmogorov differential equations [34].
To explore the influence of baseline patient characteristics on each of the transition intensities, the multi-state models were adjusted for pre-specified covariates, including treatment, number of exacerbations in previous year, eosinophil count (log-transformed), ICS use at screening, FEV1% predicted, sex and age. The pairwise interaction between eosinophil count (log-transformed) and treatment was included in the final covariate-adjusted multi-state model, based on the likelihood ratio test with a significance level < 0.05. The covariate-adjusted models assumed that hazards (transition intensities) were proportional between covariate values/patient groups. Covariate effects were allowed to vary for each transition, with the exception of transitions from exacerbation states (States 2–4) to death, where common effects were assumed due to insufficient data. The transition to death was, for the same reason, assumed to be independent of treatment and sex.
For eight patient profiles defined by FEV1 (40% vs 55%), exacerbation at baseline (1 vs 2) and age (55 vs 70 years), transition probabilities at numerous time points within and beyond the study period were calculated based on the best fitting parametric model, adjusted for covariates. These transition probabilities represented the proportion of patients incurring moderate and/or severe exacerbations or dying over up to 5 years.
Since results may be biased if the Markov assumption is violated, as a simple test, we included time spent in State 1 (no exacerbation) as a covariate in the best fitting parametric model for transitions out of exacerbation states (i.e. transitions out of States 2–4). The estimated hazard ratios were all around the null value 1 (Table S1), supporting the Markov assumption, therefore the results presented are based on the fitted Markov models.
Data were analysed using R version 3.6.3 [35]. The “survival” [36] and “mstate” [37] R packages were used for non-parametric and semi-parametric approaches. The “flexsurv” [38] R package was used to fit fully parametric Markov models via maximum likelihood estimation and predict state occupancy probabilities.