Revised estimates of influenza-associated excess mortality, United States, 1995 through 2005

Background Excess mortality due to seasonal influenza is thought to be substantial. However, influenza may often not be recognized as cause of death. Imputation methods are therefore required to assess the public health impact of influenza. The purpose of this study was to obtain estimates of monthly excess mortality due to influenza that are based on an epidemiologically meaningful model. Methods and Results U.S. monthly all-cause mortality, 1995 through 2005, was hierarchically modeled as Poisson variable with a mean that linearly depends both on seasonal covariates and on influenza-certified mortality. It also allowed for overdispersion to account for extra variation that is not captured by the Poisson error. The coefficient associated with influenza-certified mortality was interpreted as ratio of total influenza mortality to influenza-certified mortality. Separate models were fitted for four age categories (<18, 18–49, 50–64, 65+). Bayesian parameter estimation was performed using Markov Chain Monte Carlo methods. For the eleven year study period, a total of 260,814 (95% CI: 201,011–290,556) deaths was attributed to influenza, corresponding to an annual average of 23,710, or 0.91% of all deaths. Conclusion Annual estimates for influenza mortality were highly variable from year to year, but they were systematically lower than previously published estimates. The excellent fit of our model with the data suggest validity of our estimates.

Recent U.S. estimates of average annual excess mortality due to seasonal influenza exceed 30,000 [11,14].
To estimate excess mortality due to influenza, two fundamental approaches have previously been used. The most popular one is based on Serfling's seasonal regression method [18] and has resulted in numerous estimates of excess mortality due influenza [3][4][5][6][7][8]12,14,15,19]. This periodical regression approach is based on parametric estimation of a sinusoidal "baseline" function that repre-sents mortality in absence of influenza. The difference between the baseline function and the observed numbers of deaths is then interpreted as the number of excess deaths due to influenza. Typically, the baseline mortality function is fitted to weekly or monthly mortality rates or numbers during non-influenza months, using two or more Fourier terms [18]. This approach is intuitively appealing as it captures the strong seasonal periodicity of mortality. However, the particular choice of a parametric baseline function lacks epidemiological justification: Why should the baseline function be sinusoidal rather than of any other periodic form? Depending on the shape of the "true" baseline function, under-or overestimation of excess mortality due to influenza might result. If, for example, the true baseline function is "higher" (i.e. the definite integral of the true function is larger) than the assumed sinusoidal function, then overestimation would result and vice versa. Another, potentially more important shortcoming of the periodical regression approach lies in the fact that seasonally correlated causes of mortality, including influenza, are not controlled for, which might lead to confounded estimates of excess mortality.
To avoid this difficulty, one could gauge all-cause mortality with some independent measure of influenza transmission (or mortality). Following this rationale, Thompson et al. [11] estimated excess mortality due to both influenza and respiratory syncytial virus (RSV). They used a generalized linear model (GLM) with a Poisson distribution and a logarithmic link function to model the weekly number of deaths. They also used two Fourier terms in their model, but, in addition, used indicators of influenza and RSV transmission. These indicator variables were defined by the proportions of specimens testing positive for influenza A(H1N1), influenza A(H3N2), influenza B and RSV. Several potential shortcomings of this methods are, however, apparent. First, this model also makes a priori assumptions about the baseline mortality function-in this case an exponentiated sinusoidal function. Although this might conceivably be true, there is lit- tle empirical evidence to support this assumption. Second, the multiplicative form of the model implies that excess mortality, given a certain amount of influenza activity, depends on the current level of all-cause mortality. Again, this does not appear to be a well-founded assumption. Finally, the proportion of test positive specimens is likely to be a poor measure of excess mortality. While a high proportion of test positive specimens is compatible with high levels of influenza transmission (and excess mortality), this is not necessarily true. The model, however, implies that five hundred influenza positives, obtained from a thousand tests, are associated with less excess mortality than two influenza positives, obtained from three tests. This appears to be an unrealistic assumption. The seasonally changing frequency of influenza testing [20] is, at least partly, due to the seasonally changing incidence of other agents causing influenza-like illness (ILL).
Alternatively, one could postulate that mortality directly attributed to influenza (influenza-certified mortality) represents a certain proportion of all mortality attributable to influenza. This assumption implies that the coefficient associated with influenza-certified mortality represents the ratio of total influenza mortality to influenza certified mortality [17,21]. Here we use a method for the estimation of influenza excess mortality which is similar to the one recently presented by Schanzer and colleagues [17]: we adopt the proportionality assumption and avoid specific parametric assumptions about the baseline function. In addition, and deviating from the Schanzer model, we allow for random variability of influenza-certified mortality by adopting a hierarchical modeling approach. We present the resulting estimates of U.S. excess mortality due to influenza for the years 1995 to 2005. We compare these to estimates obtained from a Thompson-like model [11], as well as to previously published estimates of influenzaassociated excess mortality.

Data
We used Multiple Cause-of-Death Data for the years 1995 to 2005 (Multiple Cause-of-Death Microdata, 1995-2005, National Center for Health Statistics, Hyattsville, Maryland). This dataset is in the public domain and can be electronically downloaded from the web site of the National Bureau of Economic Research http:www.nber.org/data/vital-statistics-mortalitdatmulti ple-cause-of-death.html. We defined deaths as influenzacertified if influenza was given as underlying cause of death. The corresponding diagnostic code for ICD-9 (1995 to 1998) was 487 and and for ICD-10 (1999 onwards) the code range was J10-J12. Influenza years were defined as lasting from July 1 of one year to June 30 of the following year. We defined four age categories: < 18, 18-49, 50-64 and 65+. Observations with missing age (N = 4,490) were not included in this analysis.

Statistical model
The epidemiological model on which our analyses are based implies that monthly all-cause mortality is the sum of "baseline mortality", i.e. mortality that is independent of influenza, and mortality that is a direct or indirect result of influenza. Based on the pronounced seasonal periodicity of all-cause mortality we assume that baseline mortality is a function of calender month. In addition, we allow for a linear and/or non-linear trend in all-cause mortality that takes into account demographic or other population level changes resulting in linear/non-linear changes in baseline mortality over time. Finally, we accommodate extra variability of baseline mortality that is not accounted for by calender month and trend. This epidemiological model corresponds to the following statistical model, hence referred to as "current model": where Y i is the observed all-cause mortality count during index months i = 1, ..., 132. The variable Y i which represents a number and not a rate, is assumed to follow a Poisson distribution with a mean parameter θ i . The Poisson mean parameter θ i has an identity link and is distributed as Normal with mean μ i and variance τ. This parametrization for the Poisson mean θ i allows for overdispersion. In the implementation, θ i is restricted to positive values to ensure the positivity of the generated samples. The model for μ i has two parts. The first part concerns mortality due to non-influenza related causes (baseline mortality) which includes a random intercept for calendar month m i (m i = 1, ..., 12) that models the seasonal background mortality, and also includes linear, quadratic and cubic effects for temporal changes due to health, demographic or socioeconomic factors. The variable t i (t i = 0, ..., 10) indicates the calendar year; t i = 0 corresponds to the year 1995. The regression coefficients β 1 , β 2 and β 3 measure these changes. The second part of the model for μ i concerns mortality due to influenza. The symbol γ i is the λ m i Poisson parameter from the second level of hierarchy for the observed influenza-certified mortality, X i . The parameter ϕ measures the effect of influenza-certified mortality on all cause mortality assuming that all other effects are fixed. This is the parameter of interest. It can also be interpreted as the ratio of total influenza mortality to influenza-certified mortality. Thus, the total excess influenza mortality for index month i, , is given by To estimate excess mortality due to influenza, is calculated using expression 5, with posterior estimates of γ i and ϕ. As total influenza mortality cannot be lower than influenza-certified mortality, the minimum value for the range in the prior distribution for ϕ was set to one (see additional file 1).
To assess the performance of the current model, we also analyzed the data with a modified form of the model proposed by Thompson et al. [11]. The modified model has the following form: where β 0 is an intercept, β 1 and β 2 are defined as above, α 1 and α 2 represent the parameters associated with the Fourier terms and λ is the natural logarithm of the rate ratio associated with influenza-certified mortality. In contrast to Thompson et al. we used monthly, rather than weekly data and used observed influenza-certified mortality, rather than proportion of positive influenza tests, as indicator for total influenza mortality. For this Thompson- The parameters of the TL model could easily be estimated using a GLM procedure in any standard statistical software package. However, to allow for direct comparison of the model fit we used the same estimation procedure as for the current model. The fit of the two age-specific models was compared using the deviance information criterion (DIC) [23]. DIC penalizes the model goodness-of-fit for additional complexity. The complexity is measured by the effective number of parameters.
In order to quantify the variance explained by the fitted models, we used a Bayesian version of the classical Rsquared [24]. i.e.
where E(·) and V(·) are the operators for the posterior mean and empirical variance, respectively and e i = Y i -μ i .
The empirical variance of e is computed for each iteration.
Separate analyses were performed for the four age categories, because of substantial differences in seasonality of all-cause mortality: While seasonal periodicity is quite obvious in the oldest category, it distinctly decreases with age and becomes inapparent in the youngest category ( Figure 2).

Results
A total of 26,262,147 deaths over the period of eleven years were included in this analysis. Of these deaths, 12,922 (0.05%) were influenza-certified. The current model (1)(2)(3)(4) fit the data very well (Figure 3a) and explained most of the variance of all-cause mortality    of the age group-specific TL models combined was worse than the corresponding fit of the current model ( Figures  3a and 3b). Accordingly, the variance explained by the TL model was also lower than for our model (BRSQ for age categories < 18, 18-49, 50-64 and 65+: 0.61, 0.31, 0.84 and 0.83, respectively). Yet, the pattern of predicted total influenza deaths was virtually identical to the one predicted by the current model ( Figure 4).

Discussion
Our  [14] estimated that 25,071 deaths were attributable to influenza in ages 65+ alone. Thompson et al. [11] estimated the number of excess deaths during that influenza year at 36,280-more than three times our estimate. The obvious question arises, which of these estimates are closest to the true excess mortality? As pointed out above, the method of Simonsen et al. [14] is problematic for two reasons. First, it does not account for temporal correlation between baseline mortality and influenza excess mortality. The resulting estimates of influenza excess mortality may therefore be confounded. Second, their model makes a priori assumptions about the parametric shape of the baseline function; these assumptions may or may not be true. They should, in any event, be validated. The Thompson model [11], which superficially resembles a hybrid between the Simonsen model and the model proposed by Schanzer et al. [17] (or the current model), addresses the issue of temporal confounding by controlling for the proportion of influenza test positives. As pointed out in the Background section, the use of that specific variable to control for influenza mortality may not be appropriate. We compared estimates from the TL model with estimates from the current model. The TL model is based on the Thompson model, but influenza-certified mortality is substituted for proportion positives. Although the resulting estimates were about a sixth higher than our estimates, the seasonal pattern was highly consistent with the pattern seen with the current model. This consistency implies relative robustness of excess deaths estimates to the choice of a specific baseline function. The vast difference between our and Thompson's estimates [11] can therefore not be explained by differences in model structure, nor in the way the baseline function is modeled. They may rather be due to the use of proportion of specimens testing positive to control for influenza mortality.
Schanzer et al. [17], like us, used a Poisson model with linear (rather than logarithmic) link function, to analyze weekly mortality data from Canada. Modeling weekly mortality has the advantage of giving higher temporal resolution to the analysis. On the other hand, deaths associated with, but not attributed to influenza may occur with some delay and may thus be partially decoupled from influenza-certified mortality. However, Schanzer and colleagues did not find an obvious lag between weekly influenza-certified mortality and mortality due to other causes. Future studies will be needed to determine what level of temporal aggregation results in the best estimates.
To take into account random variability in influenza-certified mortality, we used a hierarchical model. While the point estimates for ϕ (corresponding to β 3 in [17]) obtained from a GLM are very similar to the ones obtained from the hierarchical model (21.35 and 21.16, respectively, for 65+), the confidence limits are much   . This may even be more pronounced for weekly data, where numbers of influenza-certified deaths are often quite small. To the extent that our hierarchical model takes into account random variability of influenza-certified deaths and thus leads to wider confidence limits around the resulting excess mortality estimates, it is more conservative than non-hierarchical GLM models.
The potentially most serious shortcoming of our approach to influenza excess mortality estimation relates to the possibility that influenza-certified mortality is a poor indicator of total influenza mortality. Although a death certificate diagnosis of influenza will likely only be given under strong suspicion of that cause, the diagnosis will rarely have been laboratory-confirmed and therefore is likely of low specificity. As this indicator explains mor-tality in excess of the seasonal baseline very well one could speculate that a better indicator of influenza mortality would lead to even lower estimates of influenza excess mortality. A preliminary comparison of the number of influenza A(H3N2) isolates with influenza-certified death for two influenza seasons (2002/03 and 2003/04) revealed a remarkably close correspondence of the two indicators ( Figure 5). Could it be that the number of influenza isolates informs the death certificate diagnosis? The cross-validation of various indicators for influenza mortality will be an important target of future research. The very high proportion of all-cause mortality explained by our model makes it appear quite unlikely that a substantially better model of all-cause mortality can be constructed. We therefore believe that our estimates of excess mortality are better than previous estimates, which are invariably larger than ours. Future studies of excess mortality due to influenza should also consider the possibility that the relationship between the chosen influenza indicator and total influenza excess mortality might change from season to season or even over the season. Studies examining the public health impact of influenza should address some additional issues. First, mortality may not be a good measure of the burden a disease inflicts upon a population. It may be that most deaths that are triggered by an influenza infection occur in people on the verge of dying from other causes, so that the time of death is advanced only by a short period of time. Clearly, the burden of disease would be much lighter in that case than if most people were to die from influenza prematurely by many years. A more adequate measure of disease burden therefore may be disability-adjusted life years (DALYs) [25]. By estimating a more refined age distribution of those who died from influenza, calculations of DALYs, or at least of potential years of life lost (PYLL), should be relatively straightforward. As persons with underlying illnesses are particularly vulnerable to fatal influenza infection, and also might have reduced life expectancy at the time of death, resulting PYLL estimates might be inflated. Second, by modeling monthly mortality independent of mortality during previous periods we ignore, like others, the demographic process that may lead to reduced mortality after epidemic mortality: Not only will lower mortality result from the depletion of the population at risk (smaller denominator), but influenza is also likely to disproportionately affect the frailest individuals, thus leaving the remaining population less frail and thus less susceptible. The first part of this problem could be remedied by modeling the mortality rates rather than numbers of deaths. The second one could be addressed by making specific assumptions on the frailty distribution in the population and on the association between that frailty characteristic and relative mortality risk.

Conclusion
Previous estimates of excess mortality due to influenza may be biased and inflated. We propose a two-level hierarchical Poisson model where the baseline mortality varies with time. The goodness-of-fit statistic indicates that this model fits the data very well, explaining well above 90% of the observed variation of all-cause mortality during the eleven years study period. The resulting estimates are therefore likely of high validity. Future attempts to quantify the public health burden of influenza should also explore demographic approaches that take into account life expectancy.