Estimation of the mortality burden associated with influenza relies on statistical association rather than on direct measurement of laboratory-confirmed deaths, both because it is expensive and clinically unnecessary to test all patients with respiratory symptoms, and because many sequelae of influenza that may lead to mortality (eg, bacterial pneumonia) may be diagnosed only after influenza virus has become undetectable in the host.^{1} This article describes a new statistical method for estimating influenza-related mortality.

This inference method relies on incidence proxies (introduced by Goldstein and colleagues)^{2} for the three major influenza strains (A/H3N2, A/H1N1, and B). These influenza incidence proxies are different from those used in several previous studies^{1}^{,}^{3}^{,}^{4} estimating influenza-associated mortality. The weekly incidence proxies are defined in this study as the proportion of physician visits for influenza-like illness multiplied by the proportion of respiratory specimen testing positive for particular influenza subtypes, weighted by regional populations in the United States. Under certain assumptions, this proxy should be approximately linearly related to the weekly incidence of subtype-specific influenza cases in the population.^{2} The influenza-like illness measure accounts for the changing numbers of sick individuals from whom tests are sampled, whereas the viral positivity indicator accounts for variation in the proportion of influenza-like illness, that is, in fact, influenza-related.

Those incidence proxies can be used to define a linear relation between influenza incidence and excess mortality. More precisely, weekly all-cause or cause-specific mortality data are regressed linearly against the weekly influenza incidence proxies, temporal trends, and a periodic annual baseline. We do not assume a trigonometric model to account for seasonality in mortality baselines,^{1}^{,}^{3}^{,}^{5}^{,}^{6} but instead model them using periodic cubic splines to accommodate their (a priori unknown) shape.^{7} Moreover, we allow the baseline seasonal pattern of mortality to change following the introduction of pneumococcal conjugate vaccine in the United States in 2000, which led to a reduction in invasive pneumococcal disease.^{8–11} In particular, the holiday season spikes for pneumococcal disease in older adults are reduced in the postvaccination era,^{12} which may affect the estimates of the contribution of influenza to various underlying causes of death. In addition, given the prominent role played by influenza A/H3N2 in influenza-associated mortality,^{13} we allow the impact of the A/H3N2 incidence proxy to vary over time, and to take into account the emergence of the A/H3N2 Fujian new antigenic variant strain in the United States in 2003–2004.

Besides studying influenza’s contribution to all-cause mortality and mortality underlying respiratory and circulatory causes, we examine the potential contribution of influenza to mortality for several underlying causes: diabetes, cancer, Alzheimer disease, renal disease, and chronic liver disease. The latter underlying causes were not included in the recent revision of the U.S. Centers for Disease Control (CDC) estimation procedures for flu-associated mortality.^{14}

## METHODS

### Mortality Data

Weekly mortality data were extracted from the National Center for Health Statistics database between 1997 and 2007––the most recent year for which U.S. vital statistics were available at the time of this writing. The mortality data are stratified by underlying causes, with each mortality record having a unique underlying cause.^{15} The mortality coding system changed in 1999 owing to the transition to a new revision of the International Classification of Diseases in the United States (ICD-9 to ICD-10). Conversion factors for mortality counts adjusting for the above change were adopted from the report by Anderson and colleagues.^{16}

In addition to all-cause mortality, we investigated the potential contribution of influenza to several specific underlying causes (see eTable 1 [https://links.lww.com/EDE/A611] for a list of causes and specific codes). Those causes were chosen because of their prominence in mortality data and the relatively high prevalence of the associated conditions among adults hospitalized for influenza.^{17} In addition, the same analysis was performed for deaths owing to septicemia and unintentional injuries, which we used as controls that are independent of influenza infection.

For a specific underlying cause (eg, myocardial infarction) or a collection of underlying causes (eg, all circulatory deaths or all-cause mortality), *c*, we denote by *M*_{c}(*t*) the reported mortality for cause(s) *c* in week *t*.

### Incidence Proxies and Their Estimation From Data

Each influenza season is defined as the time between week 40 of one calendar year and week 20 of the next calendar year, with 10 seasons from 1997–1998 to 2006–2007.^{18} We adopt the definition of the weekly incidence proxy for each of the three influenza subtypes from the article by Goldstein and colleagues.^{2} Specifically, for the influenza season, the weekly incidence proxy is defined as the proportion of influenza-like illness among all outpatient visits to sentinel physicians, multiplied by the proportion of respiratory viral isolates tested that were positive for a particular strain, weighted by the regional populations for the ten regions defined by the U.S. Department of Health and Human Services.^{18} For the noninfluenza season, the weekly incidence proxy for each strain is set to zero. We denote by *H*3(*t*),*H*1(*t*),*B*(*t*) the incidence proxies of influenza A/H3N2, A/H1N1, and B on week *t* in the data as defined above.

### Statistical Approach

As an initial step, we inspected the feasibility of relating mortality to the influenza incidence proxies by comparing excess pneumonia and influenza mortality above a Serfling seasonal baseline (thought to be a specific indicator of influenza-related mortality) to influenza incidence proxies shifted forward by 2 weeks.^{19} These analyses, detailed in the eAppendix (https://links.lww.com/EDE/A611), revealed good correlation between excess pneumonia and influenza mortality and the A/H3N2 incidence proxy (eFigure 1, https://links.lww.com/EDE/A611), in line with the prominent contribution of influenza A/H3N2 to mortality.^{13} In addition, a “bump” in excess mortality was observed around the first week of January in the first 4 years of our study period and disappeared in more recent years, coincidental with the introduction of the pneumococcal conjugate vaccination in 2000. Further, the ratio of excess mortality to A/H3N2 incidence was bigger before than after the 2003–2004 season, when the Fujian variant of A/H3N2 appeared.

Based on these preliminary analyses, the A/H3N2 incidence proxy was split as

Here, *H*3^{1}(*t*) is the incidence proxy of A/H3N2 before the appearance of the Fujian strain in the 2003–2004 season (defined to be zero afterward), and *H*3^{2}(*t*) is the incidence proxy of A/H3N2 starting with the 2003–2004 season. Second, we split

Here, Base_{1}(*t*) is periodic (period 52 weeks) before the 2000–2001 season and zero after that (starting week 23 of 2001), representing the baseline before the introduction of pneumococcal conjugate vaccine in the United States. Base_{2}(*t*) is zero before week 23 of 2001 and periodic thereafter, representing the baseline in the presence of the pneumococcal vaccine.^{8–12} Given that 2 years in the CDC data (1997 and 2003)^{18} have 53 epidemiologic weeks in order to make the baselines periodic, we discarded mortality data for the week in mid-July for the subsequent years (1998 and 2004), and specified the 1997–1998 and 2003–2004 influenza seasons to end on calendar week 19 rather than 20 of both 1998 and 2004.

Our basic model for relating mortality for a collection of causes *c* to influenza circulation is

Here, the temporal trend is modeled by a polynomial in the calendar year. Both Base_{1}(*t*) and Base_{2}(*t*) are assumed to be elements of the space *SB*_{4}-periodic cubic splines on the interval [0,52] (with a period of 52 weeks), with a knot every 4 weeks. The noise in our primary inference process is assumed to be white noise (ordinary least squares) inference. Finally, the operator *S* (applied to each incidence proxy) is the forward-shifting operator by the time from incidence to mortality. The distribution of the time from incidence to death is unknown and may be specific to a particular cause of death; the specific lag between incidence and mortality will also depend on details of the reporting system for influenza-like illness and viral activity. Analysis of the mean-square accuracy of the ordinary least squares fits (*R*^{2}) showed that the best fits for all causes occur when incidence is shifted forward between 1 and 2 weeks. Correspondingly, if *S*_{i} is the forward shift operator by *i* weeks, we use

and the parameter 0≤*f*≤1 is determined for each cause by maximizing the *R*^{2} of the fit. The components of the model described by equations (3) and (4) are summarized in Table 1.

In the full model given by equation (3), the regression coefficient for influenza A/H1N1 was close to zero for all categories of death, and negative for most categories studied (including all-cause mortality). We report the results of the model with influenza A/H1N1 excluded, and present the results with A/H1N1 incidence included in the inference process in eAppendix ( https://links.lww.com/EDE/A611).

### Autocorrelated Noise and the Bootstrapped Confidence Bounds

Examination of the residuals in the model given by equation (3) reveals autocorrelation for each of the studied underlying causes of death. However, even in the presence of autocorrelation––which might have a complicated, nonstationary (eg, seasonal) structure––the ordinary least squares estimates are generally unbiased.^{20} The confidence bounds for these estimates are bootstrapped by fitting the residuals to an autoregressive moving-average model––in fact, examination of the autocorrelation, as well as the partial autocorrelation structure of the residuals, suggest that modeling the residuals as an autoregressive AR(1) process is adequate. Thebootstrap works as follows in two steps: Let *M*_{c}(*t*) be the mortality data,

be the ordinary least squares estimates, *e*(*t*) = *M*_{c}(*t*) − *M*_{c}(*t*) be the residuals.

- Step 1: (Bootstrap estimate of the autoregressive parameter). Let
- be the initial estimate of the autoregressive parameter for the AR(1) noise. This estimate is generally biased downward.
^{21}The bias is estimated through the following bootstrap procedure: - Define
- where
*V**is drawn (with replacement) from {*V*(*j*)}. Define - .
*M**_{c}(*t*) is then regressed against the covariates, and an estimate ρ* of the autoregressive parameter is generated as the correlation between the residuals in the above linear regression. The bias is estimated as- ,
^{21}where - is the mean of the bootstrapped sample (ρ*). The bootstrapped estimate of the autoregressive parameter is
- .
- Step 2: (Bootstrap estimate of the regression parameters). Repeat Step 1 with
- replaced by ρ
_{boot}to obtain a sample of the regression parameters (β*).^{22}The latter sample is used to obtain confidence bounds on those parameters (Table 1) as well as on the contribution of influenza to mortality during different time periods.

### Testing the Presence of Harvesting: Association Between Influenza-Associated Mortality and Mortality in the Following Summer

It has been suggested that excess influenza-related mortality is associated with a decrease in mortality in the following summer.^{23}^{,}^{24} To investigate this further, we examined the correlation between our annual estimates of influenza-associated all-cause mortality (between week 41 of calendar year *X* and week 22 of year *X*+1 and

- (1) Summertime excess mortality (defined as observed mortality minus [baseline + trend]). Summer was defined as weeks 23–40 of year
*X*+1 - (2) Overall annual excess mortality (between week 41 of year
*X*and week 40 of year*X*+1 - (3) Wintertime annual excess mortality [(2)−(1)].

## RESULTS

### All-Cause Mortality

In this section, we estimate the contribution of influenza to all-cause mortality. The forward-shifting parameter ƒ in equation (4) is estimated to be 0.445, suggesting a balanced contribution of influenza proxies in weeks *t*-1 and *t*-2 to mortality in week *t* (*R*^{2} = 0.9613). The estimated regression coefficients (multiplied by 1,000) related to influenza subtype-specific incidence proxies are β^{2}_{H3} = 15.06 (95% CI = 12.89–17.24), β^{2}_{H3} = 10.84 (8.94–12.88), and β_{B} = 14.07 (6.72–21.45). The coefficient for A/H3N2 before 2003–2004 is almost 40% higher than that starting in 2003–2004.

The fit of the all-cause mortality model was excellent, and the model predictions were temporally consistent with the observed values, as illustrated by Figure 1.

Figure 2 plots the weekly mortality baselines Base_{1}(*t*) and Base_{2}(*t*) before and after 2001. As expected from preliminary analysis (eAppendix, https://links.lww.com/EDE/A611), the main change is the reduction of the January “bump” in baseline for the years after the introduction of the pneumococcal conjugate (PCV7) vaccine as well as a slight shift of estimated baseline mortality from fall and winter toward early spring.

Annual estimates of the incidence proxies for the three influenza strains are given in the article by Goldstein and colleagues.^{2} One can combine these incidence estimates with the estimated influenza regression coefficients to generate annual estimates (and CIs) for each strain’s contribution to all-cause mortality. The latter estimates are summarized in Table 2, and they highlight the important variability in influenza mortality burden across years, ranging from 4.96 (winter 2000–2001) to 20.49 (winter 1998–1999) per 100,000 population.

The estimates in Table 2 refer to a model without influenza A/H1N1; for comparison, a model including A/H1N1 is presented in eAppendix (https://links.lww.com/EDE/A611) section S3.

### Specific Mortality Causes

In this section, we repeat the same analysis for several specific underlying mortality causes. Table 3 gives estimates of the influenza-specific model parameters and average influenza-associated annual mortality rate during the study period.

All causes of death that were potentially associated with influenza based on their prominence in mortality statistics or influenza-related hospitalization records (including respiratory, cardiac, cancer, diabetes, renal, chronic liver, and central nervous system diseases) have substantial positive influenza estimates (Table 3). For control causes of death, coefficients were generally small and sometimes negative.

A few disjoint causes within each of the broader categories explained nearly the whole estimate for those categories. Specifically, circulatory disease, respiratory disease, cancer, diabetes, Alzheimer disease, renal disease, and chronic liver disease added up to about 84% of the all-cause estimate for influenza-associated mortality; heart disease and cerebrovascular disease estimates added up to about 97% of all circulatory flu-associated deaths; and pneumonia and influenza and lower respiratory disease estimates added up to about 96% of all respiratory flu-associated deaths.

Model fit to weekly data was temporally consistent for individual causes of death (Fig. 3), and *R*^{2} was above 83% for 13 of the 14 causes of death associated with influenza (Table 3). Chronic liver disease had the lowest model *R*^{2} (51%) and the smallest influenza contribution of all studied death categories.

Similar to the patterns seen in estimated baselines for all-cause mortality, the holiday season bump was observed for a number of underlying causes before the year 2000, especially among respiratory and circulatory causes, and nearly disappeared in recent years (Fig. 4).

### Harvesting: Influenza-Associated Mortality Versus Mortality in the Following Summer

We studied the correlation between annual influenza-associated mortality and excess mortality during different time periods, based on all-cause mortality estimates. The results for the correlation above are presented in Table 4, the relation of which is also plotted in Figure 5.

Wintertime influenza mortality was positively correlated with annual excess mortality and wintertime excess mortality: overall mortality was higher in years and in winters with higher influenza-associated mortality. Mortality in the summer after high influenza-associated winter mortality tended to be lower. However, the relationship was weak, and the wintertime flu-associated mortality was an order of magnitude larger than the estimated reduction in mortality the following summer. Overall, these findings give little evidence of substantial “harvesting” by influenza-associated mortality of persons who would have died within the year anyway.

## DISCUSSION

Although the evaluation of mortality attributable to influenza is of great interest, several aspects of the existing estimation methodologies merit reconsideration. Some of these aspects are the usage of the percent-positive proxy for the incidence of influenza, an exponential (rather than linear) relation assumed between influenza incidence and mortality,^{1}^{,}^{3}^{,}^{7}^{,}^{25} and the use of a trigonometric baseline for the annual non–flu-associated mortality.^{1}^{,}^{3}^{,}^{6} Moreover, we hypothesized that it might be important to consider two epidemiologic phenomena that occurred during the study period: the introduction of pneumococcal conjugate vaccination^{8}^{,}^{9}^{,}^{11}^{,}^{12} and the appearance of the antigenically novel Fujian A/H3N2 strain during the 2003–2004 season in the United States.

Our estimation procedure for the contribution of influenza to several underlying causes of death (including all-cause mortality) addresses the above aspects. We use weekly incidence proxies for the major influenza subtypes A/H3N2, A/H1N1, and B that combine information on the incidence of influenza-like illness and laboratory testing,^{2} and relate those incidence proxies linearly to influenza-associated mortality. In addition to respiratory and circulatory causes, we examine the potential contribution of influenza to mortality for several underlying causes frequently present in adults hospitalized for influenza^{17}: diabetes,^{26} cancer, Alzheimer disease, renal disease, and chronic liver disease. We note that some of the mortalities for other underlying causes not considered here may also be influenza-associated.

We estimate that 69% of flu-associated mortality in the United States is attributable to respiratory and circulatory causes. The model also shows a substantial contribution of influenza to mortality from cancer, diabetes, renal disease, central nervous system diseases, and chronic liver disease; the latter underlying causes account for an additional 16% of influenza-associated mortality. This stratification of influenza’s contribution to mortality is different from the recent CDC estimation,^{27} which limits influenza mortality to underlying respiratory and circulatory deaths. Moreover, we estimate the average annual influenza-associated mortality rate to be 8.2 per 100,000 individuals between 1997 and 2007 in the United States for respiratory and circulatory causes, whereas the CDC estimate^{27} is 11.4––closer to our estimate of 11.9 for the influenza contribution to all-cause mortality. The difference might stem largely from the fact that the CDC estimates^{27} use data for the whole period from 1976 to 2007 in one inference framework. If one restricts the Poisson regression estimation procedure^{3} to 1997–2007, the average annual rate of all-cause influenza-associated deaths is 13.5 per 100,000 individuals, and the average rate of annual respiratory and circulatory deaths is 8.7 per 100,000 individuals (Section S6 of the eAppendix, https://links.lww.com/EDE/A611), fairly close to our estimates.

The model indicates a change in the annual mortality baselines following the introduction of the pediatric pneumococcal vaccine (PCV7) in 2000. The timing of the change, and the fact that most of the underlying causes here are documented risk factors for pneumococcal disease, suggests that the introduction of PCV7 is a plausible explanation for this change. Between-country comparisons could take advantage of variations in the introduction of PCV7 worldwide to assess the consistency of this effect.

The model also suggests a substantial decline in the number of fatalities per 100,000 population associated with each unit increase in our incidence proxy for influenza A/H3N2 viruses circulating after the emergence of the Fujian antigenic variant in 2003. The latter difference may be attributable in part to a potentially lower impact of flu on mortality after the introduction of PCV7 (which also covers the period of the Fujian strain circulation), and to the fact that an antigenically novel strain affects children disproportionately compared with a strain that has previously circulated, lowering the overall case–fatality ratio. Given those considerations, one should be cautious in attributing a biological significance to the observed difference in the case–fatality ratios for the pre-Fujian vs. post-Fujian A/H3N2 strains.

The contribution of the influenza A/H1N1 subtype to each mortality cause was consistently small or negative for all causes of death. For this reason, we excluded A/H1N1 from the main analysis. As shown in the eAppendix ( https://links.lww.com/EDE/A611), inclusion of A/H1N1 has a very minor impact on other model parameters and resulting influenza-burden estimates. Our results are in line with the previous U.S. mortality studies that have shown that A/H1N1 consistently causes the lowest mortality burden of all three influenza subtypes.^{1}^{,}^{6} In addition, our study period was relatively short, and did not include any A/H1N1-dominant season. Further, our study ended before the emergence of the A/H1N1 pandemic virus in 2009, which might have contributed to more mortalities than the seasonal influenza A/H1N1 viruses. The mortality from this pandemic cannot yet be compared with that of previous influenza seasons owing to lack of recent U.S. vital statistics data.^{28}

The CIs for the regression coefficients for influenza A/H3N2 are narrower than those for influenza B. For all respiratory causes, pneumonia and influenza, and chronic liver disease, the contribution of influenza B is uncertain. Moreover, influenza B’s contribution to chronic lower respiratory diseases is estimated to be higher than its contribution to all respiratory causes, further emphasizing the lack of reliability in the latter estimate. Our results suggest that mortality attributable to influenza A/H3N2 is 3.8-fold higher than that attributable to influenza B during the study period. This is consistent with several earlier studies showing that A/H3N2-dominant seasons have two- to four-fold higher mortality than seasons dominated by other subtypes, including influenza B.^{1}^{,}^{2}^{,}^{6}^{,}^{27}

The basic measure for goodness of fit in our model was *R*^{2} for the linear regression. This measure dictated some of the choices made in the model, such as the distribution of time from incidence to mortality for the various causes of death and the selection of polynomial coefficients for the temporal trend. Additionally, *R*^{2} was the principal measure of the model’s accuracy in the sensitivity analysis described in eAppendix (https://links.lww.com/EDE/A611). This analysis revealed, in particular, that use of proxies for influenza incidence based only on data from tests of respiratory specimens results in lower *R*^{2} values in the linear model and higher estimates of influenza-associated mortality.

The residuals in the ordinary least squares model used for inference were found to be autocorrelated. Maximum likelihood inference using autoregressive moving-average models for the noise produced somewhat smaller estimates for the regression coefficients for each subtype and each underlying cause compared with the ordinary least squares estimates. This suggests that some bias is possible if the autoregressive moving-average structure for the noise is assumed in the model inference for our sample (538 observations and 33 regression coefficients).^{29} Correspondingly, (the unbiased) ordinary least squares inference was used for our estimates of the contribution of influenza to mortality, and bootstrapping based on the autoregressive AR(1) model for the noise^{21}^{,}^{22} was used to estimate the confidence bounds.

Although the relationship between influenza and excess wintertime mortality is well-established, it has also been suggested that excess wintertime mortality reduces the heat effect on mortality in the following summer.^{23}^{,}^{24} We have found that flu-associated mortality is substantially larger than summertime excess mortality in the following summers, suggesting that if there is a harvesting effect, it is relatively small. Although influenza-associated mortality and summertime mortality were moderately negatively correlated, the correlation was not significant, whereas the correlation between annual flu-associated mortality and overall excess mortality was positive and significant.

Overall, our methodology renders temporally consistent estimates of flu-attributable mortality for a number of underlying causes, some of which have received little attention previously. Those estimates are “almost” additive in the sense that separate estimates for disjoint causes that make up most fatalities in a certain collection of causes nearly add up to the contribution of influenza for that whole collection. The linear relation between influenza incidence and mortality improves the mechanistic plausibility of the traditional regression models using a Poisson link between influenza incidence and deaths. Our model allows for the assessment of the role of related epidemiologic phenomena, such as changes in mortality baselines after the introduction of PCV and changes in case–fatality ratios associated with the appearance of the novel Fujian A/H3N2 strain. It suggests that, although some influenza deaths may occur in persons who would have died anyway within the coming year, this “harvesting” effect is a modest proportion of total influenza-attributable mortality. Our estimation method of the baseline for noninfluenza-associated mortality may be helpful in assessing the influenza-associated mortality during pandemic periods, such as the 2009 H1N1 pandemic. Finally, the linear relation between influenza incidence and excess mortality can be used in a prediction framework^{2} to explore the potential ramifications of an evolving influenza season, particularly one dominated by influenza A/H3N2.

## ACKNOWLEDGMENTS

*We thank Dan Weinberger and Lone Simonsen for their helpful discussions.*