INTRODUCTION
The burden of hospitalizations and deaths caused by influenza is generally estimated using ecological trend studies (e.g. [Reference Alling, Blackwelder and Stuart-Harris1–Reference Thompson10]). These studies use population-level rates of outcomes that are not specific to influenza, such as all-cause mortality. Outcome rates during time periods when influenza did not circulate are used to predict rates due to non-influenza causes during periods of influenza circulation. During influenza seasons, the difference between observed rates (which are due both to influenza and to non-influenza causes) and predicted rates (due to non-influenza causes only) is then attributed to influenza. Accurate estimates of the burden of influenza from these studies depend on accurate estimates of the predicted rates of non-influenza outcomes.
A key challenge to estimating non-influenza rates in these ecological studies is confounding by season [Reference Jackson11]. Influenza tends to circulate in winter, which is when many other seasonal causes of morbidity and mortality also peak. Ecological study designs must account for this seasonal confounding; a variety of approaches have been used for this (e.g. [Reference Choi and Thacker2, Reference Foppa and Hossain4, Reference Serfling7, Reference Simonsen8, Reference Thompson10, Reference McBean, Babish and Warren12]). Importantly, while the assumptions that underlie these methods have been described [Reference Jackson11], few data exist on whether these methods accurately account for seasonal confounding. Because seasonal influenza viruses circulate every winter in temperate regions, winter rates of outcomes in the absence of influenza are usually not observable. So it is unknown whether ecological studies successfully estimate the winter incidence of outcomes due to causes other than influenza.
The emergence and global circulation of the 2009 influenza A(H1N1) pandemic virus (2009 H1N1pdm) provided a unique opportunity to test the accuracy of ecological study designs. The 2009 H1N1pdm virus circulated early relative to the typical influenza season in the United States, with peak circulation occurring in September/October, and with influenza circulation essentially absent by mid-December 2009 [13]. Therefore, observed morbidity and mortality during the winter of 2009–2010 represents winter baselines of these events occurring in the absence of influenza. We used this winter to test the accuracy of common ecological methods used to estimate the winter incidence of outcomes due to causes other than influenza.
METHODS
Study population
We conducted this study within the Vaccine Safety Datalink (VSD) Project, a collaboration between the Centers for Disease Control and Prevention (CDC), America's Health Insurance Plans, and ten geographically diverse healthcare systems (‘sites’) [Reference Baggs14]. In combination, the VSD contains data on site enrolment, healthcare utilization, and mortality for about 3% of the US population [Reference Baggs14]. Our study used data from seven VSD sites with complete enrolment, demographic, hospitalization, and mortality data for the study period (1 September 1997 to 31 August 2010): Kaiser Permanente of Northern California (NCK; Oakland, CA); Kaiser Permanente of Colorado (KPC; Denver, CO); Health Partners Research Foundation (HPM; Minneapolis, MN); Marshfield Clinic Research Foundation (MFC; Marshfield, WI); Kaiser Permanente Northwest (NWK; Portland, OR); Kaiser Permanente of Southern California (SCK; Los Angeles, CA); and Group Health Cooperative (GHC; Seattle, WA).
We defined our study cohort as all seniors enrolled in one of the study sites between 1 September 2002 (KPC) or 1 September 1997 (other sites) and 31 August 2010. Seniors began contributing person-time to the study following their first year of continuous enrolment, and continued to contribute person-time until the earliest of death, disenrolment, or the study end date of 31 August 2010. We restricted our study population to seniors (adults aged ⩾65 years), for two reasons. First, older adults were less susceptible to 2009 H1N1pdm than were other age groups, due to cross-protective antibodies from A(H1N1) influenza strains that circulated prior to 1957 [Reference Hancock15]. Thus, seniors would be unlikely to experience delayed health outcomes in winter that could have resulted from influenza infections in autumn. Second, seniors are at high risk of influenza complications and are the most common group for whom influenza burden is estimated (e.g. [Reference McBean, Babish and Warren12, Reference Brinkhof16, Reference Nguyen-Van-Tam17]).
Health outcomes
The primary health outcomes studied were all-cause mortality; hospitalizations due to pneumonia or influenza (PI); and hospitalizations for respiratory or cardiovascular (RC) conditions. We also studied hospitalizations for acute myocardial infarction (AMI) as a secondary outcome. We identified all onset dates of these health outcomes (which aggregate influenza-attributed and non-influenza-attributed events) in study population members during the entire follow-up period. We determined dates of death from administrative records at the participating sites, which combine death data from multiple sources (including hospital discharge records, state mortality records, and enrolment databases).
Hospitalization outcomes were defined by International Classification of Diseases, version 9, Clinical Modification (ICD-9-CM) codes assigned to inpatient visits: codes 480–487 (PI hospitalizations), codes 390–519 (RC hospitalizations), and code 410 (AMI hospitalizations). These codes were chosen for consistency with previous studies [Reference Thompson9, Reference Nichol18].
Influenza circulation data
We used data on positive influenza tests from the US World Health Organization and National Respiratory and Enteric Virus Surveillance System (WHO/NREVSS) collaborating laboratories. Publically available WHO/NREVSS surveillance data are stratified by each of the ten Department of Health and Human Services regions of the United States. For each site, we calculated the percent of specimens testing positive for influenza, overall and by type and subtype [A(H1N1), 2009 H1N1pdm, A(H3N2), and B], from the region in which the site is located as a measure of influenza circulation. For each year and region we defined influenza seasons as the consecutive weeks with at least 10% of isolates testing positive for influenza [Reference Thompson9].
Analysis
Individual-level data on person-time and counts of health outcomes were aggregated by age group (65–69, 70–74, 75–79, 80–84, ⩾85 years), sex, site, and study week. The aggregated weekly data at each site were merged with weekly influenza data from the corresponding region. We defined the prediction period as the weeks of 14 December 2009 to 1 March 2010, as >95% of influenza infections for 2009–2010 occurred prior to 14 December 2009 [13]. The remaining weeks were the baseline period used to fit the models.
For each of the health outcomes we fit three different statistical models to the data. We first estimated each model's parameters using the baseline period. We then used the baseline model to predict the expected incidence during the prediction period. Because the 2009–2010 winter was effectively influenza-free, the predicted rates are predictions of the outcome rates due to causes other than influenza.
The first statistical model was a cyclic regression model, which was first introduced by Serfling in 1963 [Reference Serfling7]. In this approach, data from weeks when influenza circulated are removed from the time series. A cyclic regression model, using sine and cosine terms to represent seasonality, is then fit to the remaining data. Non-influenza incidence during periods when influenza circulates is interpolated from the model parameters, and differences between the observed and predicted influenza season incidence are attributed to influenza. We fit the following Poisson regression model to the data from weeks when influenza did not circulate, modelling the weekly count of events as a function of calendar time:
$${Y_t} = {\alpha }\,{\rm exp}\left[ \matrix{{{\beta }_{0j}} + {{\beta }_1}\left( t \right) + {{\beta }_2}\left( {{t^2}} \right) + {{\beta }_3}\left( {{\rm sin}\left( {2\pi t/k} \right)} \right) \hfill \cr +\; {{\beta }_4}\left( {{\rm cos}\left( {2\pi t/k} \right)} \right) + {{\beta }_5}\left( {male} \right) + {{\bi \beta }_{\bf 6}}({\bf age}) \hfill} \right],$$
where Y t is the number of events during week t, k is the period of the time series (k = 52·177 for weekly data), β 0j is the site-specific intercept (i.e. random intercept model), β 6 is a vector of parameters for the age strata, and α is the offset terms for log of weekly person-time. Predicted incidence in the total population was calculated by computing a weighted average of the stratum-specific predicted incidences. Because model-based standard errors do not account for the autocorrelation of the data, we calculated 95% confidence limits using seasonal block bootstrapping [Reference Politis19, Reference Weinberger20].
The second statistical model we used was the ‘virological’ regression model, which uses the Serfling model as a foundation and adds data on influenza circulation. This model has been the standard approach used by the CDC for estimating the burden of influenza from ecological studies since 2003 [Reference Thompson9, Reference Thompson10, 21, Reference Zhou22]. In the virological regression model, all time points during the baseline period are used, including weeks from both influenza and non-influenza seasons. Parameters are included for percent of tests positive for each influenza type and subtype. For consistency with the standard use of these models [Reference Thompson10, 21, Reference Zhou22], we did not include lagged effects of influenza. Incidence of non-influenza outcomes during influenza season is then predicted from the model parameters, setting the influenza covariates to zero. As with the Serfling model, we calculated 95% confidence limits using a seasonal block bootstrap.
The third statistical model we tested was an autoregressive, integrated moving average (ARIMA) time-series model [Reference Box and Jenkins23]. Use of these models for predicting the burden of influenza has been described in detail elsewhere [Reference Choi and Thacker2, Reference Jackson11]. In brief, an ARIMA model assumes that the incidence rate at time t (Y t ) is influenced by a ‘random shock’ (α t ) to the population at time t. The random shock is the cumulative effect of all factors affecting incidence, such as weather, temperature, pathogens, and air pollution. The effect of the random shock may persist for several time periods, so incidence at time t may depend on prior random shocks, α t–p for some values of p. In addition, Y t may depend on prior incidence, Y t–q for some values of q. Conceptually, this could be due to depletion of susceptibles during the early stages of an epidemic, leaving fewer susceptibles in later stages. Thus, an ARIMA model has the form:
$${Y_t} = {\o_1}\left( {{{\alpha} _t}} \right) + {{\o}_2}\left( {{{\alpha} _t}_{-1}} \right) + \ldots + {{\o}_{p}}\left( {{{\alpha} _{t-p}}} \right) + {\theta _1}\left( {{Y_{t-}}_1} \right) + \ldots + \theta {q}\left( {{Y_{t-q}}} \right).$
ARIMA models may also use differencing to remove trend or drift in the time series. In differencing, the dependent variable is the differenced time series Z t :
$${Z_{t}} = Y_t-{Y_{t-d}},$
where d represents the level, or order, of differencing. ARIMA models can include seasonal lags in α t and in Y t to account for cyclic trends in incidence. ARIMA models are described by the orders of autoregressive, differencing, and moving average terms; for example, a (p,d,q)(1,0,0) model refers to a model with a first-order autoregressive term and no differencing or moving average terms; upper-case (P,D,Q) are used to refer to seasonal terms. Finally, ARIMA models can include other time series as independent variables; these models are sometimes referred to as ‘ARIMAX’ models. In our study, we used ARIMAX methods to model outcome incidence, where we included weekly counts of positive influenza tests as independent variables.
Unlike the cyclic regression models, where the model covariates may be chosen a priori, ARIMA and ARIMAX models are built empirically. The modeller identifies values for p, q, and d that are specific to the time series being modelled. For each outcome and site, we first fit an ARIMA model to the weekly observed data during the entire baseline period. We then added weekly percents of tests positive for influenza (by type and subtype) as predictors if they were significantly associated with the outcome after fitting the initial ARIMA model and if their inclusion in the model did not decrease the fit of the ARIMA model to the data. We used the resulting ARIMAX model to estimate (forecast) weekly outcome incidence rates during the prediction period.
Because ARIMA models are fit to a single outcome time series, an ARIMAX model cannot be adjusted for age or sex. Stratifying by age would have required fitting separate models for each age stratum, a fivefold increase in the number of models to fit. Instead, we fit an ARIMAX model to each outcome time series at each site, aggregated across all age/sex groups. To test whether stratifying by age might improve the accuracy of the ARIMAX model, we fit separate ARIMAX models to each of the five age groups for the PI outcome in one site (NCK), and compared the overall predicted incidence from the unstratified model with the combined predicted incidence from the stratified models. Forecasts from the unstratified model differed from the age-stratified models by <0·5 cases/10 000 person-years (data not shown).
Accuracy endpoints
The study endpoints to assess accuracy of the statistical methods were the errors between the observed and predicted rates of the health outcomes during the prediction period. We compared the predicted incidence of each outcome during the prediction period with the observed incidence in each site. We quantified the prediction accuracy of the statistical models by calculating the difference between the observed and predicted incidence rates. We calculated prediction error as both absolute differences and relative differences as a percentage of the observed incidence. We used US census data to standardize prediction error to event counts in the US population and compared annual predicted deaths in the United States from the virological regression model to recent CDC estimates for adults aged ⩾65 years based on the same model [21] during 1997–1998 to 2006–2007, years for which predictions were available in the present study and the CDC estimates.
This study was approved by the Institutional Review Boards of all participating sites. Analyses were conducted using SAS v. 9.2 (SAS Institute, USA) and Stata version 12 (StataCorp, USA).
RESULTS
We observed a total of 10 947 081 person-years of follow-up time during the study period, of which 31% was in seniors aged 65–69 years, 26% was in seniors aged 70–74 years, 20% was in seniors aged 75–79 years, 13% was in seniors aged 80–84 years, and 10% was in seniors aged ⩾85 years (Table 1). Our study population experienced 408 437 deaths; 149 630 PI hospitalizations; 1 507 965 RC hospitalizations; and 102 810 AMI hospitalizations during the study period. Incidence of all health outcomes fluctuated seasonally (Fig. 1).
Fig. 1. Observed weekly incidence rates per 10 000 person-years, for (a) deaths; (b) pneumonia/influenza (PI) hospitalizations; (c) respiratory circulatory (RC) hospitalizations; (d) acute myocardial infarction (MI) hospitalizations. Grey bars indicate prediction periods.
Table 1. Distribution of person-time and outcomes by site, age, sex, and influenza year
NCK, Kaiser Permanente of Northern California (Oakland, CA); KPC, Kaiser Permanente of Colorado (Denver, CO); MFC, Marshfield Clinic Research Foundation (Marshfield, WI); HPM, Health Partners Research Foundation (Minneapolis, MN); NWK, Kaiser Permanente Northwest (Portland, OR); SCK, Kaiser Permanente of Southern California (Los Angeles, CA); GHC, Group Health Cooperative (Seattle, WA).
* Influenza years run from 1 September to 31 August.
In virological regression models, influenza A(H3N2) and B were significantly associated with deaths and with PI and RC hospitalizations. Influenza A(H1N1) was only significantly associated with RC hospitalizations, while influenza A(H1N1)pdm was only associated with PI hospitalizations. The exact form of the final ARIMAX models varied by health outcome and by site, but the most typical model included first-order differencing and a first-order moving average as well as a seasonal moving average term [i.e. a (p,d,q)(P,D,Q)(0,1,1)(0,0,1) model]. One or two influenza parameters [typically A(H3N2) or B] were usually included in the final models for death and for PI and RC hospitalizations. Influenza parameters were never significantly associated with AMI hospitalizations in either virological or ARIMA models. Yearly predicted US influenza-related deaths from the virological regression model during 1997–1998 to 2006–2007 ranged from 8531 to 36 972 and were well correlated (R 2 = 0·74) with CDC's estimated deaths during the same time period, which ranged from 10 800 to 43 727.
During the prediction period, Serfling and virological regression models were more accurate that the ARIMAX models for all health outcomes except AMI hospitalizations (Table 2). When averaged across all seven sites, all three statistical methods predicted non-influenza mortality rates during the winter of 2009–2010 with reasonable accuracy (Table 2). All three methods underestimated non-influenza mortality slightly, by 5% (ARIMAX), 2% (Serfling), and 2% (virological). Accuracy for the PI hospitalization outcomes was worse for the ARIMAX model, with 10% under-estimation. Accuracy was high for the RC hospitalization outcome for all three methods. However, confidence intervals were wide for all outcomes and all methods and spanned 0% prediction error. For example, the Serfling method predicted winter mortality incidence was 2% lower than the observed incidence, but the confidence interval on this prediction error was −21% to 17%.
Table 2. Observed and predicted health outcome rates per 10 000 person-years and prediction errors, with 95% confidence intervals, using three statistical methods, during the influenza-free winter of 2009–2010
PI, Pneumonia/influenza; RC, respiratory/circulatory, AMI, acute myocardial infarction.
We age- and sex-standardized the virological prediction errors to the US population. In a typical 12-week influenza season, a 2% underestimate of non-influenza mortality corresponded to overestimating deaths due to influenza by 9694 deaths. The 4% underestimate of PI hospitalizations corresponded to overestimating influenza-related PI hospitalizations by 5381.
DISCUSSION
This study provides new insights into the accuracy of methods used in ecological studies of the burden of influenza. In a winter with negligible influenza circulation, we found that all three of the methods predicted winter incidence of the health outcomes with fairly high accuracy. This finding was unexpected. The accuracy of the Serfling method in particular is surprising. This cyclic regression makes the assumption that the winter peak in non-influenza incidence is exactly equal to the summer trough in duration and amplitude [Reference Jackson11, Reference Gilca24]. Seasonal changes in mortality and hospitalizations are likely influenced by numerous factors that vary over time in a complex way, such as temperature, weather, air pollution, hours of daylight, and circulation of other pathogens. A priori, we did not expect the simplistic cyclic regression function to accurately account for seasonal changes in mortality and in hospitalizations. Our findings suggest that changes in seasonal incidence of mortality and hospitalizations are primarily driven by factors that vary with the same timing and amplitude from year to year, such as hours of sunlight per day (e.g. [Reference Kriszbacher25]).
We also found that the cyclic regression models, with or without viral circulation data, performed better than the ARIMAX model in predicting winter non-influenza incidence. This is also somewhat surprising, as ARIMAX models were developed specifically for handling some of the unique features of time-series data, such as trends over time, seasonal fluctuations, and autocorrelations [Reference Box and Jenkins23]. We expected that fitting an ARIMAX model to each individual time series would result in better prediction than applying identical cyclic regression models to each time series. The fact that cyclic regression models appear to predict non-influenza incidence as well as ARIMAX models is also evidence that the seasonal variation in mortality and hospitalizations can be well modelled by a simple cyclic regression function.
Despite the high accuracy of the Serfling and virological regression methods, our results suggest that caution is needed in using these methods to estimate the burden of influenza. A 2% underestimate in non-influenza mortality corresponds to overestimating US influenza-related deaths by over 9694 deaths per year. A recent study using virological regression estimated that influenza causes an average of 21 098 deaths per year [21]. Thus, our results suggest that nearly half of this estimate could be attributable to prediction error rather than influenza. Because influenza only accounts for a small proportion of all winter deaths, even a small error in estimating non-influenza deaths can lead to large errors in deaths attributed to influenza. By contrast, Thompson et al. [Reference Thompson9] used virological regression to estimate that influenza causes an annual average of 66 373 PI hospitalization in US adults aged ⩾65 years. In our study, the virological method overestimated PI hospitalizations by 5300, which implies that prediction error only accounts for 8% of the PI hospitalizations attributed to influenza.
A potential limitation of this study is that, in contrast to typical years, influenza viruses were circulating intensely during autumn 2009. It is possible that seniors who would typically have had influenza-related hospitalizations or deaths during winter were affected instead in autumn. This effect in turn might cause models to overestimate 2009–2010 winter mortality on the basis of atypically high mortality in autumn 2009. However, we think this is unlikely. Influenza viruses circulating in autumn 2009 were almost entirely 2009 H1N1pdm, which caused comparatively little morbidity and mortality in seniors [13] due to cross-protective antibodies from influenza A(H1N1) strains that circulated before 1957 [Reference Hancock15]. Thus, the impact of the autumn pandemic wave on deaths in seniors was modest at best, and should not substantively affect model predictions. A second limitation is the VSD population used for this study represents only 3% of the total population of seniors in the United States. However, at the time we began this study in 2011, the mortality and hospitalization data that are commonly used for US burden of influenza studies were not yet available for the winter of 2009–2010. Extrapolating from our study population to the entire US for 1997–1998 to 2006–2007 gave similar mortality estimates as a recent study based on the entire US population [21], which increases our confidence that our sample is representative of the United States as a whole. The smaller sample size may lead to wider confidence intervals than are found in similar studies that use data for the entire US population, although similar studies have often ignored autocorrelation when estimating standard errors [Reference Foppa and Hossain4, 21, Reference Zhou22] or have not reported standard errors [Reference Thompson9, Reference Thompson10]. We cannot rule out chance as an explanation for our findings that the methods tend to underestimate non-influenza morbidity and mortality.
Estimating the burden of morbidity and mortality caused by influenza remains a matter of public health importance. Our study suggests that ecological estimates of non-influenza outcome rates may be sufficiently accurate when applied to health outcomes where the contribution of influenza is large. By contrast, these models are probably not sufficiently accurate for outcomes where the contribution of influenza is low.
ACKNOWLEDGEMENTS
This work was supported by a subcontract with America's Health Insurance Plans (200–2002–00732) from the Centers for Disease Control and Prevention.
DECLARATION OF INTEREST
None.
