The Association Between Rate and Severity of Exacerbations in Chronic Obstructive Pulmonary Disease: An Application of a Joint Frailty-Logistic Model.

Exacerbations are a hallmark of chronic obstructive pulmonary disease (COPD). Evidence suggests the presence of substantial between-individual variability (heterogeneity) in exacerbation rates. The question of whether individuals vary in their tendency towards experiencing severe (versus mild) exacerbations, or whether there is an association between exacerbation rate and severity, has not yet been studied. We used data from the MACRO Study, a 1-year randomized trial of the use of azithromycin for prevention of COPD exacerbations (United States and Canada, 2006-2010; n = 1,107, mean age = 65.2 years, 59.1% male). A parametric frailty model was combined with a logistic regression model, with bivariate random effects capturing heterogeneity in rate and severity. The average rate of exacerbation was 1.53 episodes/year, with 95% of subjects having a model-estimated rate of 0.47-4.22 episodes/year. The overall ratio of severe exacerbations to total exacerbations was 0.22, with 95% of subjects having a model-estimated ratio of 0.04-0.60. We did not confirm an association between exacerbation rate and severity (P = 0.099). A unified model, implemented in standard software, could estimate joint heterogeneity in COPD exacerbation rate and severity and can have applications in similar contexts where inference on event time and intensity is considered. We provide SAS code (SAS Institute, Inc., Cary, North Carolina) and a simulated data set to facilitate further uses of this method.

RCT Entities:   Population Interventions Outcomes

Chronic obstructive pulmonary disease (COPD) is a common chronic disease of the airways characterized by the progressive loss of lung function (1). COPD is the third leading cause of years of life lost in the United States (2), and evidence suggests that its burden in comparison with other common chronic diseases is increasing (3). The main symptoms of COPD include shortness of breath, cough, and sputum (4). Periods of intensified disease activity (typically short in duration), referred to as exacerbations, are a prominent feature of COPD and a major cause of morbidity, mortality, and economic burden (4). COPD exacerbation is the number 1 cause of medical hospitalizations in Canada (5).

The average annual rate of exacerbations in well-characterized COPD cohorts is between 0.85 and 1.30 episodes (6–8), but individual patients vary significantly in their propensity to experience exacerbations, with the “frequent exacerbator” trait considered to be relatively stable over time (7). Exacerbations also vary in severity, with current consensus reports categorizing them as mild (symptom-based), moderate (requiring a physician visit and/or treatment with oral corticosteroids), or severe (requiring hospitalization) (9). Approximately 20% of COPD exacerbations are severe (10). Whether individuals vary in their tendency to experience more severe exacerbations versus milder exacerbations largely remains unknown.

Prevention, early detection, and appropriate treatment of exacerbations are a focal point of attention in COPD care and research. In the presence of significant between-individual variability (heterogeneity) in the rate and severity of exacerbations, it is likely that the benefit-risk profile of prevention and early detection strategies for exacerbations will be a function of individuals’ background exacerbation risk and their propensity to experience severe exacerbations. Use of the ubiquitous proportional hazards model, with time to first exacerbation as the outcome, is a common mode of inference in contemporary clinical trials of COPD. While it is robust in estimating treatment effect in randomized controlled trials, this analytical method falls short of providing other features, such as background rate of exacerbations or the shape of the incidence function, to enable predictions about the rate and (absolute or relative) duration of time to future events for a given patient. As mentioned by Cox et al. (11), making such informative predictions has been hindered by the widespread use of semiparametric proportional hazards models.

Parametric frailty models are a class of survival models that enable the analysis of recurrent events, can explicitly model heterogeneity across individuals through random-effect terms, and can quantify the incidence rate and the shape of the incidence function, thus providing all of the quantities necessary for specification of the complete natural history of recurrent events (12). While there are dedicated computer routines to fit certain groups of such models, general nonlinear mixed-model optimizers, such as PROC NLMIXED in SAS (SAS Institute, Inc., Cary, North Carolina), give the analyst the necessary power to fully specify the likelihood function for statistical inference in such models (see chapter 15 in SAS for Mixed Models (13)). The added benefit of such generic optimizers is that the model can be modified to better address the specific question at hand (14).

Our overall objective in the present study was to quantify jointly the heterogeneity in the rate and severity of COPD exacerbations, with the specific objective of testing whether there is an association between exacerbation rate and severity. Our hypothesis was that persons who frequently experience exacerbation also tend to experience a higher proportion of severe exacerbations. This hypothesis was based on the assumption that being more vulnerable to the causes and drivers of COPD exacerbation (inflammation and/or infection (15)) would increase individuals’ risk of experiencing more frequent exacerbations as well as more severe events. A modification of parametric shared frailty models enabled us to conduct simultaneous modeling of the rate and severity of exacerbations to address these objectives.

METHODS

The data

This study was approved by the University of British Columbia's Human Ethics Board. The analysis was based on the data from the MACRO Study (16), a multicenter North American clinical trial designed to determine whether azithromycin, an antibiotic, was effective in reducing COPD exacerbations. Details on the study and its findings are presented elsewhere (16). In brief, the study randomized 1,142 patients to placebo or azithromycin (250 mg/day). Eligible patients included those aged 40 years or older with a clinical diagnosis of COPD who were either using continuous supplemental oxygen therapy or had received systemic glucocorticoids, had visited an emergency room, or had been hospitalized for COPD exacerbation within the previous year. The primary endpoint was time to first exacerbation, categorized into mild (requiring symptomatic treatment), moderate (requiring oral corticosteroids), severe (requiring inpatient care), or very severe (requiring intubation). The time of exacerbation was considered to be the date on which treatment was initiated. Follow-up time was 385 days.

Statistical analysis: joint modeling of rate and severity of exacerbations

We specified a joint survival-logistic model to simultaneously estimate the association between clinical features and the rate and severity of exacerbations. The framework used is generally similar to one described by Berridge and Whitehead (17), who investigated the efficacy of a treatment for headache in terms of both time to recurrence of headaches and their severity. The authors fitted a nonparametric proportional hazards model for time to first event and separately fitted an ordinal regression model for severity categories. In the present work, we used a joint parametric recurrent-event and logistic regression model to enable full quantification of exacerbation incidence and severity and their correlation. Compared with the alternative framework of modeling exacerbations of different severities as separate types of events, this framework is more conceptually aligned with the conventional wisdom that exacerbations share a common pathophysiology (15). In addition, while enabling us to model exacerbation severity, this model retained the ability to allow direct inferences on the occurrence of any exacerbations, which is the outcome of interest in many COPD studies (18). Finally, this fully parametric model enabled us to specify the likelihood (see Web Appendix 1, available at http://aje.oxfordjournals.org/) and to use generic algorithms that are developed for fitting nonlinear mixed models (13).

The rate component

The basis of our model was a parametric random-effects (frailty) accelerated failure time model. For the ith individual, the instantaneous exacerbation rate at time t (hazard) is with being the vector of observed, time-fixed characteristics (covariates), β the vector of regression coefficients, and z an unobserved normally distributed zero-mean random-effect term that is specific to each person and captures between-individual variability in exacerbation rate over and beyond the variability due to observed characteristics. λ0 is the baseline hazard function. The association between covariates and the event of interest in accelerated failure time models is modeled in terms of “speeding up” or “slowing down” the time to the event (19). Accelerated failure time models are popular models for predicting time to events and event rates. For example, they were the basis of the original Framingham Study risk prediction equations for cardiovascular disease (20).

Time zero in this analysis was the randomization date. Because our interest was in estimating the cumulative incidence of exacerbation from this time onward, a “total time” model for recurrent events was used (21). The total time model enables quantification of complete event incidence by specifying the hazard function from a given point forward, as opposed to the “gap time” analysis, which resets the hazard function after each event (21). The total time approach does not directly model any association between the occurrence of an exacerbation and the risk or severity of future ones. Instead, the model assumes that individuals have intrinsic properties determining the relationship between their exacerbation rate and severity, enabling statistical testing of the study's specific hypothesis.

For a parametric survival model, a baseline hazard (or survival) function needs to be specified (11). Common choices include Weibull (which includes the constant-hazard exponential survival model as a special case), lognormal, and log-logistic survival functions (11). We evaluated all such functions and used the Akaike Information Criterion as an objective measure to find the best-fitting function (22). We also visually inspected the predicted cumulative number of exacerbations versus the observed cumulative number of exacerbations to check the goodness-of-fit of the chosen survival function. With the Weibull survival function, the accelerated failure time model is also a proportional hazards model (23, p. 45), and thus regression coefficients can be presented in terms of both the acceleration factor and the more familiar hazard ratio (Web Appendix 1) (19). Due to the small number of deaths, and in line with the primary analysis of the MACRO Study (16), we did not adjust for the competing risk of death in the main analysis but investigated its effect in a sensitivity analysis.

The severity component

The second component was a random-effects logit (logistic) model for the severity of exacerbations. While this could be a multilevel ordinal logit model capturing the 4 levels of exacerbation severity, we felt that the most natural distinction lies between severe-to-very-severe exacerbations and mild-to-moderate exacerbations, as the former type requires expensive inpatient care, while the latter can be managed in a physician's office or by the patient at home. Thus, we used a binary categorization contrasting severe/very severe (referred to collectively as severe) exacerbations with mild/moderate (referred to collectively as mild) exacerbations. This model predicts the probability that the jth exacerbation in the ith individual (denoted by Y) will be severe (coded as 1, compared with 0 for mild exacerbation): with being covariates (all time-fixed), β′ the regression coefficients, and a normally distributed zero-mean random-effect term that is specific to each person and captures between-individual variability in the proportion of severe exacerbations (out of the total number). Regression coefficients from this component can be expressed in terms of odds ratios associating the covariate with the probability of an exacerbation's being severe. Persons with at least 1 exacerbation contribute to the severity component.

The same set of covariates was used for both the survival and logistic components and consisted of age at baseline, sex, treatment group assignment, smoking status at baseline (current smoker vs. ex-/never smoker), whether the patient had been hospitalized for a COPD exacerbation during the 12-month period before the study, whether the patient had received home oxygen therapy during the 12-month period before the study, degree of airflow obstruction (measured in terms of forced expiratory volume in 1 second at baseline; the higher the forced expiratory volume in 1 second, the lower the degree of obstruction), and COPD-specific health status (captured by total score on the St. George's Respiratory Questionnaire at baseline; the higher the score, the lower the health status) (24).

To examine the extent of heterogeneity, we calculated the coefficient of variation, defined as the ratio of the standard deviation to the mean, for model-estimated individualized exacerbation rate as well as the proportion of severe exacerbations. We also determined the lower and upper bounds of these quantities that contained 95% of the sample. Individualized rates and proportions were calculated using the maximum-likelihood estimates of the fixed-effect and empirical Bayes estimates of the random-effect terms (13).

This framework was also used to test the hypothesis that across individuals, there is a relationship between the rate of exacerbation and the severity of exacerbation. The 2 random-effect terms are governed by 3 parameters (2 variance parameters and 1 covariance parameter). A positive (negative) covariance indicates that persons with a higher rate of exacerbation tend to have a higher (lower) risk of their exacerbations being severe.

We performed sensitivity analyses to evaluate the robustness of the results. These included treating death as an event to evaluate its impact as a competing risk (25), as well as fitting the rate and severity components separately. SAS, version 9.4, was used for statistical analysis, with PROC NLMIXED used for maximum likelihood estimation (13). All tests of statistical significance were 2-sided at a significance level of 0.05. The SAS code and a simulated data set are provided in Web Appendix 2 and on our website (http://resp.med.ubc.ca).

RESULTS

Of the 1,142 persons recruited, 1,117 MACRO participants had at least 1 follow-up visit and were included in the original analysis (16). Of these, 10 had missing or invalid data for outcome or independent variables and were excluded from this analysis. The final sample constituted 1,107 individuals (mean age = 65.2 years; 59.1% male) who were followed for a mean of 343.7 days (0.94 years) and experienced 1,601 exacerbations, 348 (21.7%) of them severe. The rate of exacerbation was 1.53 episodes per patient-year for any exacerbation and 0.33 for severe exacerbations. The baseline characteristics of the final sample and their follow-up statistics are provided in Table 1. During the course of the study, 21 and 18 persons died in the placebo and treatment groups, respectively.

Table 1.

Baseline Characteristics and Follow-up Statistics for 1,107 Chronic Obstructive Pulmonary Disease Patients, According to Treatment Status, MACRO Study, 2006–2010[a]

Variable	Placebo Group (n = 552)			Treatment Group (n = 555)			Total (n = 1,107)
	No. of Persons	%	Mean (SD)	No. of Persons	%	Mean (SD)	No. of Persons	%	Mean (SD)
Baseline Characteristics
Male sex	327	59.2		327	58.9		654	59.1
Age, years			65.5 (8.5)			64.8 (8.7)			65.2 (8.6)
Current smoker	127	23.0		117	21.1		244	22.0
Oxygen therapy during prior year	323	58.5		332	59.8		655	59.1
Hospitalization during prior year	281	50.9		276	49.7		557	50.3
FEV1, L			1.12 (0.52)			1.10 (0.50)			1.11 (0.51)
SGRQ score[b]			51.7 (16.2)			52.8 (16.3)			50.6 (16.4)
Follow-up Statistics
Follow-up time, years			0.94 (0.18)			0.94 (0.18)			0.94 (0.18)
No. and rate of COPD exacerbations[c]
Total	878	1.69		723	1.39		1,601	1.53
Mild/moderate	681	1.31		572	1.10		1,253	1.20
Severe/very severe	197	0.38		151	0.29		348	0.33
Exacerbation frequency[d]
0	172	31.2		238	42.9		410	37.0
1	166	30.1		144	25.9		310	28.0
2	85	15.4		72	13.0		157	14.2
3	56	10.4		48	8.7		104	9.4
4	30	5.4		19	3.4		49	4.4
5	19	3.4		14	2.5		33	3.0
≥6	24	4.3		20	3.6		44	4.0

Abbreviations: COPD, chronic obstructive pulmonary disease; FEV1, forced expiratory volume in 1 second; SD, standard deviation; SGRQ, St. George's Respiratory Questionnaire.

a Data are presented as mean (SD) for continuous variables and number of subjects (% of column total) for dichotomous variables, except where noted.

b Between 0 and 100, with a higher score indicating worse status.

c Numbers in the % column show the annual rate of exacerbations (episodes/patient).

d Number of participants with the specified number of COPD exacerbations during follow-up.

Among the 3 different survival functions examined, the best fit was from the Weibull model, both in terms of the Akaike Information Criterion (Web Table 1) and in terms of the observed versus expected cumulative number of exacerbations during follow-up (Figure 1). As such, the Weibull model was chosen, and the results are reported in terms of both time acceleration and hazard ratio. The match in the observed and expected cumulative number of severe exacerbations also indicated that the logistic regression model could explain the distribution of severity in the sample. The addition of random-effect terms for both rate and severity further improved model fit, indicating the presence of unexplained heterogeneity in both the rate and severity of exacerbations (Web Table 1). The estimated shape parameter of the Weibull distribution was 0.967 (95% confidence interval (CI): 0.920, 1.013). Because the Weibull model with a unity shape parameter corresponds to the exponential model, the incidence function of exacerbations was compatible with a constant rate over the study period.

Figure 1.

Observed number of chronic obstructive pulmonary disease (COPD) exacerbations versus predicted number of exacerbations (A) and number of severe exacerbations (B) during follow-up time, MACRO Study, 2006–2010. The lower slope of the curve on the right side of the x-axis is due to attrition occurring towards the end of the study.

Table 2 provides the maximum likelihood estimates of model parameters. Treatment with azithromycin reduced the rate of exacerbation by 23% (hazard ratio (HR) = 0.77, 95% CI: 0.67, 0.89; P < 0.001) but had no effect on the severity of exacerbations (odds ratio (OR) = 0.93, 95% CI: 0.62, 1.38; P = 0.703). Among covariates, male sex was associated with a lower rate of exacerbation (HR = 0.84, P = 0.019) but a higher likelihood of exacerbations being severe (OR = 1.84, P = 0.005). A history of COPD-related hospital admission in the previous year was associated with a higher rate of exacerbation (HR = 1.48, P < 0.001) as well as increased severity of exacerbations (OR = 2.93, P < 0.001). Other factors that were significantly associated with exacerbation rate included being a smoker at baseline (HR = 0.82, P = 0.036), baseline forced expiratory volume in 1 second (per 1-L increase, HR = 0.82; P = 0.014), and COPD-related health status (per 10-unit increase in St. George's Respiratory Questionnaire score, HR = 1.12; P < 0.001). Among these variables, only the latter was associated with severity (per 10-unit increase, OR = 1.15; P = 0.030).

Table 2.

Maximum-Likelihood Estimates of the Parameters in a Joint Rate-Severity Model of Chronic Obstructive Pulmonary Disease Exacerbations, MACRO Study, 2006–2010

Variable	Rate Component					Severity Component
	HR	95% CI	β for AFT[a]	95% CI	P Value[b]	OR	95% CI	P Value[b]
Treatment group (vs. placebo)	0.77	0.67, 0.89	1.30	1.13, 1.50	<0.001	0.93	0.62, 1.38	0.703
Male sex (vs. female)	0.84	0.72, 0.97	1.20	1.03, 1.40	0.019	1.84	1.20, 2.82	0.005
Age at baseline (per 10-year increase)	0.96	0.88, 1.05	1.05	0.96, 1.14	0.332	0.97	0.76, 1.24	0.799
Current smoker[c] at baseline	0.82	0.68, 0.99	1.23	1.01, 1.48	0.036	1.29	0.76, 2.16	0.343
Oxygen therapy[c] during year before baseline	1.05	0.90, 1.23	0.95	0.81, 1.11	0.517	1.26	0.80, 1.98	0.320
Hospitalization[c] during year before baseline	1.48	1.28, 1.70	0.67	0.58, 0.77	<0.001	2.93	1.92, 4.47	<0.001
Baseline FEV1 (per 1-L increase)	0.82	0.69, 0.96	1.23	1.04, 1.46	0.014	0.74	0.46, 1.21	0.234
SGRQ score[d] (per 10-unit increase)	1.12	1.07, 1.17	0.89	0.85, 0.93	<0.001	1.15	1.01, 1.31	0.030

Abbreviations: AFT, accelerated failure time; CI, confidence interval; FEV1, forced expiratory volume in 1 second; HR, hazard ratio; OR, odds ratio; SGRQ, St. George's Respiratory Questionnaire.

a The regression coefficient y for covariate x means that persons with x = 1 experience as many exacerbations in y time units as persons with x = 0 experience in 1 time unit (unlike the case with the HR, the higher the acceleration factor the lower the event rate).

b P values were considered significant at the 0.05 level. All P values and confidence limits were computed from the final Hessian matrix based on a t distribution with default degrees of freedom (number of subjects minus number of random effects) in SAS NLMIXED.

c Binary variable.

d Between 0 and 100, with a higher score indicating worse status.

There was substantial unexplained variability (heterogeneity) in both the rate and severity of exacerbations (variance for random-effect terms: 0.55 (95% CI: 0.41, 0.68; P < 0.001) for rate and 2.25 (95% CI: 1.24, 3.25; P < 0.001) for severity). The correlation coefficient for the random-effect terms for rate and severity was −0.18 (95% CI: −0.40, 0.03; P = 0.099); thus, the null hypothesis that there is no association between exacerbation rate and severity could not be rejected at the significance level of 0.05.

Figure 2 provides a histogram of the distribution of the individualized rate of exacerbations and the individualized proportion of severe exacerbations (out of the total number). The coefficient of variation for the individualized rate of exacerbation was 0.69. Ninety-five percent of the subjects had an individualized rate between 0.47 episodes per year and 4.22 episodes per year. The coefficient of variation for the individualized proportion of severe exacerbations was 0.81, with 95% of the sample having an individualized proportion between 0.04 and 0.60.

Figure 2.

Individualized rate of chronic obstructive pulmonary disease (COPD) exacerbation (episodes/year) (A) and individualized ratio of severe exacerbations to total exacerbations (B), MACRO Study, 2006–2010. Values are model-estimated for the rate component and for the severity component (see text and Web Appendix 1 for full description). Values for random-effect terms (z and ) were obtained using the empirical Bayes method, implemented with the PREDICT statement of SAS PROC NLMIXED.

Results of sensitivity analyses are provided in Table 3. In general, the major findings of the study remained robust against alternative assumptions and design features.

Table 3.

Results From Sensitivity Analyses of Rate-Severity Models of Chronic Obstructive Pulmonary Disease Exacerbations, MACRO Study, 2006–2010

Scenario	Effect of Treatment on Rate		Effect of Treatment on Severity		Correlation Coefficient for Rate and Severity			CV for Exacerbation Rate	CV for the Probability of Exacerbation Being Severe
	HR	95% CI	OR	95% CI	ρ	95% CI	P Value[a]
Death defined as an event[b]	0.77	0.68, 0.89	0.96	0.65, 1.41	−0.20	−0.42, 0.00	0.054	0.68	0.79
Marginal models[c]	0.77	0.67, 0.89	0.92	0.62, 1.36		N/A		0.68	0.85

Abbreviations: CI, confidence interval; CV, coefficient of variation; HR, hazard ratio; N/A, not applicable; OR, odds ratio.

a The P value for the correlation coefficient was calculated using the delta method (default option in the ESTIMATE statement of SAS PROC NLMIXED).

b Severity coding: 0 = mild or moderate exacerbation; 1 = severe or very severe exacerbation or death.

c Two separate random-effect models were fitted for the rate and severity components.

DISCUSSION

Using a joint frailty-logistic model, we quantified the extent of heterogeneity in the rate and severity of COPD exacerbations and quantified the association between several clinical features and exacerbation rate and severity. The frailty (rate) component of the model, which modeled the cumulative incidence of exacerbations during follow-up, produced an effect estimate for treatment with azithromycin (HR = 0.77) that was slightly less favorable than the effect reported in the original study (HR = 0.73), which modeled duration of time to the first exacerbation (16). The results of the logistic (severity) component did not support the presence of an association between treatment and severity of exacerbations. We documented substantial levels of heterogeneity in both exacerbation rate and severity over and beyond variability due to clinical characteristics and treatment group assignment. The joint survival-logistic model could be fitted in standard statistical software and can have potential applications in many other contexts where time to recurrent events and intensity of events are both of interest.

Our a priori hypothesis of a positive association between exacerbation rate and severity was not confirmed. The weight of evidence, while not achieving statistical significance, was towards a negative association. A negative association could have several possible explanations. A plausible mechanism is a subject-specific threshold effect. Persons who seek care only when their symptoms (e.g., dyspnea) are severe will have a low overall rate of exacerbation but experience a high proportion of severe exacerbations. Another explanation is reverse causation, in which severe exacerbation might result in a change (e.g., intensification) of care, which may reduce the rate of subsequent exacerbations. A third explanation is that severe exacerbations and mild/moderate exacerbations are different (but related) phenotypes with different etiologies and molecular drivers.

To our knowledge, ours is the first study to have used a frailty model to quantify the extent of heterogeneity in the COPD exacerbation rate, as well as the first one to have investigated the presence and extent of heterogeneity in exacerbation severity. Previously, Hurst et al. (26) evaluated 2,138 COPD patients followed for up to 3 years and concluded that individuals differed in their rate of exacerbations. However, the investigators only partially quantified the degree of heterogeneity in rate, by categorizing individuals into those who had 0, 1, and 2 or more exacerbations in the first year of follow-up. The occurrence of a previous exacerbation was the strongest predictor of future exacerbations in a multinomial regression model for the categories of exacerbation frequency. Our approach in modeling exacerbation rate was different in that the variation in rate that was unexplained by observable characteristics was explicitly modeled through a random-effect term. In this framework, the previous rate of exacerbation, while not a covariate in the regression model, carries information about the value of the unobserved random-effect term and thus becomes associated with future exacerbations. In addition, through a second random-effect term, we captured heterogeneity in the severity of exacerbations and its relationship with exacerbation rate. The observed heterogeneity in exacerbation severity implies that a patient's previous history of exacerbation severity (which provides information about the value of the unobserved random-effect term for severity) can be informative in determining the severity of future exacerbations. This framework enabled us to fully quantify the spectrum of heterogeneity (Figure 2) and has the capacity to make “individualized” predictions of the rate, timing, and severity of future exacerbations that incorporate the uncertainty in the unobserved value of the random-effect terms in a given patient.

The limitations of the present study should be acknowledged. The confidence interval for the correlation coefficient between exacerbation rate and severity contained relatively strong negative values. Studies with larger sample sizes might be able to provide more precise estimates of the association between rate and severity. Only patients with a previous history of COPD exacerbation were enrolled in the MACRO Study (16). This was to ensure that there would be enough exacerbations during follow-up time to provide sufficient statistical power at a reasonable sample size. This renders the MACRO cohort relatively homogeneous in terms of exacerbation rate, and it is expected that heterogeneity will be even higher in the general COPD population. However, it can be argued that the utility of a risk-stratification tool for future exacerbations is the greatest among COPD patients who have experienced exacerbation at least once in the past, thus justifying the use of estimates based on this sample. Another limitation of the present study was that we could not model all potentially relevant clinical characteristics because they were not captured in MACRO. For instance, previous studies have found respiratory symptoms (e.g., wheezing), history of gastroesophageal reflux, and blood biomarkers (white blood cell count, fibrinogen) to be associated with exacerbation rate (26). It is likely that some component of the heterogeneity in the rate and severity of exacerbations may be explained by inclusion of these and other covariates. Finally, death is a competing risk for COPD exacerbation. Fewer than 4% of participants died during the study, and the sensitivity analysis that included death as an event produced similar results.

In addition to severity, the duration of exacerbations is another determinant of their burden. Exacerbations can last for a few days to a few weeks, during which time an individual is not at risk for another exacerbation (18). There are extensions of survival models that can take the duration of events into account. For example, Xue and Brookmeyer (27) have proposed bivariate random-effect models that can jointly model between-episode and within-episode intervals. Yan and Fine (28) employed a “temporal process regression” framework to model marginal means of the cumulative number of exacerbations, number of days with exacerbations, and proportion of follow-up time spent in an exacerbation episode among patients with cystic fibrosis (29). Extension of the present method to include exacerbation duration, or extension the above-mentioned methods to include exacerbation severity, may result in more comprehensive models that capture 3 fundamental aspects of exacerbations: their occurrence, duration, and severity. This could be the topic of further research.

In addition to variable exacerbation duration, some evidence suggests that COPD patients are more likely to experience another exacerbation during the period immediately after a previous one (30). As such, the hazard function is a complex function of time in the short period after the occurrence of exacerbations. These temporal changes in hazard do not, per se, result in biased estimates from a model for cumulative incidence. However, they may affect the accuracy of model-based predictions if the patient has recently experienced an exacerbation. A gap time analysis could potentially provide more flexibility in modeling postexacerbation hazard as well as possible causal associations between exacerbations (21), but such a model is more suitable for predictions related to the next exacerbation than for characterizing exacerbation history. In addition, the follow-up time, starting from randomization, is a left-truncated interval (with unknown duration) in this type of analysis, causing additional challenges in proper gap-time model specification. It can be argued that in a total time model, the randomization time is an arbitrary milestone for defining time zero. While this is generally a valid concern, the compatibility of the fitted model with a constant hazard model makes this argument less relevant in this context.

Bearing these limitations in mind, the results of this study can have important implications. The findings indicate that COPD patients fall on a 2-dimensional intrinsic rate-severity spectrum, with significant variability in both dimensions. The future burden of exacerbations in a patient can be predicted by knowing where the patient belongs on this spectrum. This can be partially determined through observed characteristics and the patient's previous history of exacerbations. The significant unexplained variability in both rate and severity also means that attempts in finding other predictors (e.g., biomarkers) can potentially significantly refine our abilities to risk-stratify patients. Developing and validating prediction models for such risk stratification could be the focus of future research.

Despite the devastating burden of COPD, many features of the natural history of this disease are not well understood. COPD exacerbations (especially severe ones) are watershed moments in the course of COPD but are poorly characterized (31). In this analysis, we demonstrated substantial between-individual variability not only in the rate of exacerbations but also in the severity of exacerbations. From a methodological perspective, frailty models for recurrent-event data such as COPD exacerbations are underutilized. Compared with conventional survival analysis methods, they can more optimally use the available data and can provide important insights about the extent and determinants of heterogeneity. On a broader scale, more nuanced analyses of COPD outcomes aimed at improving our understanding of the natural history of COPD and enabling risk prediction are required.