The regular seasonality of influenza in temperate countries is recognized, but regional differences in patterns of influenza-related mortality are poorly understood. Identifying patterns could improve epidemic prediction and prevention. The authors analyzed the monthly percentage of deaths attributable to pneumonia and influenza among people aged 65 or more years in the contiguous United States, 1968–1998. The local Moran's I test for spatial autocorrelation and correlograms assessing space-time synchrony within each influenza season were applied to detect and to characterize mortality patterns. Western US regions experienced epidemics of greater magnitude than did eastern regions. Positive spatial autocorrelation (two-sided p = 0.001) revealed the similarity in influenza mortality of neighboring states, with several western states forming a focus of high mortality. In transmission seasons dominated by virus subtype A(H3N2), mortality was correlated at a high and consistent level across the United States (mean correlation = 0.56, standard deviation = 0.134). However, when subtype A(H1N1) or type B dominated, the average synchrony was lower (mean correlation = 0.23, standard deviation = 0.058). These novel analyses suggest that causes of spatial heterogeneity (e.g., large-scale environmental drivers and population movement) have impacted influenza-associated mortality.
Influenza is highly seasonal in temperate regions, with most activity occurring during winter (1, 2). Despite this regularity, interannual differences in the onset, duration, and magnitude of each transmission period exist. There has been little systematic examination or explanation of regional influenza spatiotemporal patterns, which could enhance our understanding of influenza transmission and seasonality as well as build capacity for epidemic prediction and prevention.
A few empirical studies have explored spatiotemporal patterns of influenza. Viboud et al. (3) analyzed epidemic synchrony across Australia, France, and the United States from 1972 to 1997 and identified a high correlation in impact and timing between the latter two countries. Bonabeau et al. (4) modeled influenza spread within France, 1988–1995, finding that illness rapidly diffused across the country, with homogeneous global mixing occurring before epidemics could peak in local patches. Mugglin et al. (5) modeled patterns in Scotland during the 1989/1990 transmission season, identifying two clusters in the north and south of a high-population corridor. Paget et al. (6) found evidence for a west-to-east spread of peak influenza activity across 20 European countries in three of five recent winters (1999/2000–2003/2004). Sakai et al. (7) characterized influenza-like illness patterns across Japan from 1992 to 1999 by use of kriging and simple statistics, and they described nationwide epidemic patterns spreading in concentric circles from the west-central to the east.
The factors driving influenza epidemics in different parts of the world may be diverse, so observations from these European and Asian studies cannot be directly applied to understanding influenza epidemiology in the United States. In 1987, Lui and Kendal (8) briefly explored differences in pneumonia and influenza mortality across four broad US regions, but there has been no such effort recently within a large and ecologically heterogeneous country over a lengthy time series. The United States is rich with demographic and environmental diversity that could affect influenza patterns, yet it maintains a vital statistics reporting system that is relatively standardized across states, allowing for reliable regional comparisons. Studying many transmission seasons could allow identification of persistent patterns, which might be associated with more predictable (and potentially modifiable) risk factors than transient patterns.
Demonstration of region-specific patterns that cannot be explained by differences in demographic characteristics could suggest a relation between influenza mortality and environmental variables. Pronounced seasonal patterns and a simultaneous appearance of cases across widespread locales have motivated the search for climatic drivers of infectious diseases (9). Indeed, the La Niña phase of the El Niño Southern Oscillation has been linked to greater influenza-associated morbidity and mortality in the United States and France (10–12). If disease-associated environmental variables have an element of predictability (as does climate), then early warning systems for epidemic spread could be improved (13). Space-time pattern recognition in influenza activity is essential before associations with environmental variables can be defined. Accordingly, we analyzed time and space patterns of pneumonia and influenza mortality among elderly US residents with a test for spatial autocorrelation and a method for assessing space-time synchrony.
MATERIALS AND METHODS
During interpandemic periods, influenza mortality occurs primarily among the elderly (14). Most of these deaths reflect secondary bacterial lung infection (15), so pneumonia is usually included when analyzing epidemic patterns. Thus, we analyzed pneumonia and influenza mortality from the US Multiple Cause-of-Death files, 1968–1998 (16), among people aged 65 or more years. The following four variables were extracted from each record: month of death, age group (quinquennial age groupings beginning at age 65, for a total of seven categories), state of residence, and underlying cause of death (i.e., the condition initiating the sequence of events leading to death).
For coding the underlying cause of death, the International Classification of Diseases (ICD), Eighth Revision, revised for use in the United States (17), was used for 1968–1978, and the ICD, Ninth Revision (18), was used for the remaining years, 1979–1998. In the ICD, Eighth Revision, pneumonia and influenza were coded as 470–474 and 480–486, while the corresponding codes in the ICD, Ninth Revision, were 480–487.
To calculate death rates, we obtained annual US population estimates for those aged 65 or more years from 1970 to 1998 from the US Census (http://wonder.cdc.gov/census.html, accessed May 2004). Pneumonia and influenza deaths for each month were divided by the corresponding year's total population count, with the 1970 population applied to deaths occurring in years 1968–1970.
The burden of influenza is often studied by use of the pneumonia and influenza death rate or number of excess deaths (3, 19–25). Excess mortality measures are typically used to pinpoint when mortality exceeds a threshold above a seasonal baseline, and they are useful as an early quantitative index of an influenza epidemic. Our purpose was not surveillance, but rather a retrospective study of mortality trends and spatial patterns. Excess mortality calculations are sensitive to the methods and period used to estimate the expected baseline, and the selection of a threshold above which an epidemic is declared is essentially arbitrary (26). Furthermore, a living population at risk is required. These counts rely on accurate census enumeration, and undercounts may vary by age or geographic location, thereby introducing biases when comparing subgroups.
Another measure that avoids these potential limitations is the “pneumonia and influenza percentage,” defined at a monthly resolution as the pneumonia and influenza deaths/ (total deaths – deaths due to accidents, injuries, poisonings, or violence). We excluded accidents from analysis, as have other studies of mortality (27–30), as they are unrelated to influenza and peak in the summer and fall (31). To determine the burden of influenza-associated mortality each transmission season, we summed pneumonia and influenza deaths from December to March and then divided by the total deaths for these same 4 months. The pneumonia and influenza percentage was compared with both death rates and the number of excess all-cause deaths in persons aged 65 or more years each season (25).
The Centers for Disease Control and Prevention uses a measure similar to the pneumonia and influenza percentage in conducting routine influenza surveillance, fitting a baseline curve to a time series of the percentage of all deaths due to pneumonia and influenza in 122 US cities (32). The purpose of that effort is to determine whether an epidemic threshold has been exceeded in a particular week. The pneumonia and influenza percentage also has been used to develop other methods for detecting influenza activity and epidemics (33, 34).
A(H3N2) seasons are known to be associated with greater mortality than are other subtypes (35), so this was considered throughout analysis. Following the method of Reichert et al. (24), seasons were dichotomized as having been dominated by either subtype A(H3N2) or by subtype A(H1N1) or type B. We did not evaluate any effects of multiple circulating subtypes with interference and cross-subtype protection. The following spatiotemporal analyses were restricted to the 48 contiguous US states and Washington, DC (henceforth, “48 states”).
The pneumonia and influenza percentage was entered into Space-Time Intelligence System (STIS), version 1.03, software (TerraSeer, Inc., Crystal Lake, Illinois) at a state-level spatial resolution. STIS was used to calculate Moran's I statistic (36) to detect spatial autocorrelation in the subtype-specific pneumonia and influenza percentage for each month from December to March, averaged across 18 A(H3N2) seasons and 12 A(H1N1)/B seasons. A positive spatial autocorrelation value for Moran's I statistic would indicate that nearby states had similar pneumonia and influenza mortality levels; if levels were dissimilar, the statistic would be negative. To produce a reference distribution for evaluation of the test statistic, 999 Monte Carlo randomizations of the observations were generated.
After identifying spatial autocorrelation in pneumonia and influenza mortality in the United States overall, we performed a local Moran analysis, such that states with autocorrelation significant at α = 0.05 were identified on a map. The Simes correction to the p value accounted for the lack of independence of neighboring areas. There are four possibilities for local autocorrelation: 1) high-high (a state and its neighbors had high pneumonia and influenza mortality); 2) low-low (a state and its neighbors had low pneumonia and influenza mortality); 3) high-low (a state was a high outlier among states with low pneumonia and influenza mortality); and 4) low-high (a state was a low outlier among states with high pneumonia and influenza mortality). The local Moran analysis is useful for identifying specific clusters and spatial outliers.
Data from each state were aggregated into nine climate regions (henceforth “regions”) as defined by the US National Climatic Data Center (NCDC) by use of temperature and precipitation distributions for groups of contiguous states (37). To highlight temporal trends, we fit a LOESS smoothing curve (38) to the time series for each region by use of a first-degree local polynomial and a smoothing parameter of 0.3 in R, version 1.9.1, software (R Foundation for Statistical Computing, Vienna, Austria). The monthly values of this smoothing curve were subtracted from the absolute values of the pneumonia and influenza percentage, to detrend and center each of the nine regions around zero. This centering procedure allowed for regional comparisons of the magnitude of pneumonia and influenza mortality, while minimizing any effect of regional differences in background mortality outside the influenza season.
To determine whether variation in age structure distorted regional comparisons, we standardized the monthly pneumonia and influenza percentage of each region according to the age distribution of mortality within elderly residents of the contiguous United States in the year 1968. This permitted a comparison of age-adjusted pneumonia and influenza percentages across regions.
Space-time synchrony in mortality patterns
Measures of synchrony of the monthly pneumonia and influenza percentage across states during the influenza season were estimated by use of the Sncf function in the “ncf” package for R (O. N. Bjørnstad; software downloadable from http://asi23.ent.psu.edu/onb1/) (39–42). Correlograms representing the synchrony of the pneumonia and influenza percentage from December to March over distance among the 48 states' centroids were generated and analyzed by subtype. Subtype-specific associations between the average synchrony and the slope of the synchrony were calculated. The slope captured the degree to which synchrony in pneumonia and influenza mortality lessened with distance between states and was defined as the range between the 2.5 percent and 97.5 percent quantile estimates of average synchrony. A simple regression model was constructed to determine the relation between synchrony and subtype after adjustment for the magnitude of pneumonia and influenza mortality.
A(H3N2) seasons have a strong mortality signal and are of considerable public health interest. We further focused synchrony analysis within A(H3N2) seasons only. We consulted an antigenic map developed by Smith et al. (43) of the HA1 domain of A(H3N2) viruses collected from 1968 to 2002 to determine which transmission seasons occurred immediately following an antigenic cluster transition. We classified such seasons as experiencing “major” drift. Subsequent A(H3N2) seasons within any antigenic cluster were designated as having “minor” drift. Within A(H3N2) seasons, the space-time synchronies for seasons with major versus minor drift were compared.
Validation of pneumonia and influenza percentage and effects of subtype and age
Transmission season (December to March) pneumonia and influenza percentages (16), excess all-cause deaths for each season (25), and the dominant circulating influenza subtype (24) between 1968/1969 and 1997/1998 show complex variation in terms of interannual variability in epidemic magnitude and unpredictable cycling of subtypes (table 1). Nevertheless, the Pearson correlation coefficient between the pneumonia and influenza percentage and excess all-cause deaths over the 30 seasons was high (r = 0.83, two-sided p < 0.0001), indicating the similarity of these two measures. They retained their strong correlation even when stratified by dominant circulating subtype (18 A(H3N2) seasons: r = 0.78, p < 0.0001; 12 A(H1N1)/B seasons: r = 0.80, p < 0.002). The overall correlation between the pneumonia and influenza percentage and the pneumonia and influenza death rate was r = 0.98 (p < 0.0001).
December to March pneumonia and influenza percentage (%)
No. of excess all-cause deaths*
Dominant circulating influenza subtype†
December to March pneumonia and influenza percentage (%)
No. of excess all-cause deaths*
Dominant circulating influenza subtype†
The seasonal pneumonia and influenza percentage was greater during the 18 seasons when A(H3N2) was dominant (mean = 5.21 percent, standard deviation (SD) = 0.5 percent), compared with the 12 seasons when A(H1N1) or B was dominant (mean = 4.30 percent, SD = 0.7 percent). The nonparametric Wilcoxon two-sample test comparing these distributions was significant (two-sided p = 0.005). This significant difference in severity would be even more pronounced if the excess deaths, rather than total pneumonia and influenza mortality, related to each subtype were compared.
Older age groups had progressively higher pneumonia and influenza percentages compared with younger elderly (figure 1). Other than a general aging of the population nationwide, the age structure within the elderly did not vary in any notable, spatially dependent way. Age therefore would not distort interregion comparisons of the same influenza seasons.
Moran's I test for spatial autocorrelation
We tested whether each month's pneumonia and influenza percentage, averaged across seasons by dominant subtype, was spatially autocorrelated. Moran's I statistic for each of the eight month-subtype combinations ranged from 0.297 (p = 0.002, December; A(H1N1)/B) to as much as 0.461 (p = 0.001, January; A(H1N1)/B), showing that nearby states tended to have similar pneumonia and influenza mortality levels. We then tested for local autocorrelation to identify specific clusters of states with high or low pneumonia and influenza mortality (figure 2). On each of these eight maps, states were identified as “high-high” in the western regions and Great Plains. Thus, these states with a high pneumonia and influenza percentage were adjacent to other states with a high pneumonia and influenza percentage. A few eastern states were also identified as “low-low.”
Magnitude and trends of pneumonia and influenza mortality by climate region
The southwest, west north central, west, and northwest climate regions exhibited consistently higher pneumonia and influenza percentages than did the other five regions (figure 3). The southeast region experienced the lowest burden of pneumonia and influenza mortality. To explore temporal trends, we fit a smoothing curve to the monthly mortality for each of the nine NCDC regions (figure 3). Although all regions experienced some increase in smoothed pneumonia and influenza mortality over the time series, the western region had a remarkably steeper increase. This region includes California and Nevada, but it was an increase in pneumonia deaths in those aged 80 or more years in California that was responsible for that trend. After centering the time series for each region around zero by detrending, the four western regions still averaged a higher December to March pneumonia and influenza percentage (1.5 percent) versus the five eastern regions (1.2 percent). The regional differences in pneumonia and influenza percentage could therefore not be explained by spatial variation in background levels of pneumonia and influenza mortality.
The pneumonia and influenza percentage for each region was age adjusted to the 1968 distribution of deaths. The mean December to March pneumonia and influenza percentages for the four western regions (range: 4.6–5.3 percent) all remained higher than those for the five eastern regions (range: 4.0–4.4 percent). Explicit consideration of age structure did not change our conclusions about regional pneumonia and influenza mortality differences; our other spatial and temporal analyses did not stratify by elderly age group.
Synchrony of pneumonia and influenza mortality by dominant influenza subtype
The monthly pneumonia and influenza percentages for all 30 transmission seasons and each state, along with latitudinal and longitudinal coordinates of state centroids, were submitted to the nonparametric correlation function for spatiotemporal data. When A(H3N2) dominated, synchrony tended to be relatively high (average(ρij) = 0.56, SD = 0.13) and was remarkably consistent across the country during some seasons. When A(H1N1)/B dominated, the average synchrony was lower (average(ρij) = 0.23, SD = 0.13). For some seasons, states close together had higher synchrony than did those farther apart. To illustrate these contrasts, we present the spatial dependence of synchrony for two seasons with different dominant circulating subtypes (figure 4).
Each season's average synchrony was compared with the slope of synchrony over distance (figure 5, top). The slope for the 18 seasons when A(H3N2) dominated (mean = 0.21, SD = 0.058) was smaller than the slope for the 12 seasons when A(H1N1) or B dominated (mean = 0.24, SD = 0.045). For A(H3N2) seasons, as the average regional synchrony in a season decreased, the slope in synchrony across distance became steeper (r = −0.72, p = 0.0006). For A(H1N1) or B seasons, the correlation between the average regional synchrony and slope was not significant.
There is considerable overlap in the magnitude of mortality, as measured by the pneumonia and influenza percentage, between more severe A(H1N1)/B seasons and milder A(H3N2) seasons. However, we see markedly more separation between subtypes in the average synchrony (figure 5, bottom). This cannot be explained if the A(H1N1)/B lack of synchrony were simply due to lack of severity. After adjustment for the magnitude of mortality each season, the dominant circulating subtype is still a significant factor (β = 0.24, p = 0.0001) in a regression model predicting synchrony. We know that there is a true difference in the magnitude of mortality between A(H3N2) and A(H1N1)/B seasons (35), and using our index there is an even more prominent separation in the average synchrony associated with these subtypes. The two-sample Wilcoxon test statistic comparing the distributions of average regional synchrony by subtype was 86.5 (p < 0.0001).
Within the 18 A(H3N2) seasons only, 10 major drift events were identified (1968/1969, 1972/1973, 1975/1976, 1977/1978, 1980/1981, 1987/1988, 1989/1990, 1993/1994, 1996/1997, and 1997/1998; figure 1) (43). The other A(H3N2) seasons were, by definition, associated with minor drift events. The December to March average regional synchrony in major drift seasons (n = 10; mean = 0.63, SD = 0.095) was greater than that of minor drift seasons (n = 8; mean = 0.48, SD = 0.132). This difference was statistically significant (two-sample Wilcoxon test statistic = 50.0; p < 0.02). There was no clear association between synchrony and the antigenic distance between viral clusters.
This study presents novel descriptions of space-time patterns using a test for spatial autocorrelation and a spatiotemporal correlation function. Results suggest that climate regions experienced different burdens of seasonal mortality, neighboring states had similar mortality levels, and the degree of regional correlation in mortality patterns varied according to the dominant circulating subtype.
Others have confirmed a recent rise in influenza-associated mortality in the elderly, particularly in the 1990s, attributing the increase to population aging, A(H3N2) prominence, and a trend for influenza viruses to be detected for longer periods (23). We identified a dramatic increase in pneumonia and influenza mortality in California that would be interesting for further investigation. This trend would have been hidden to analyses conducted on mortality data for the country as a whole.
Temporal trends may also have been affected by changes in both the prevalence of underlying medical conditions that predispose to mortality from pneumonia and influenza and the proportion of mortality attributable to other causes. Partial protection against influenza A(H3) viruses acquired before 1892 may have impacted trends in the very elderly (44), while immunologic adaptation to A(H3N2) after the 1968 pandemic may explain the decline in influenza-related deaths among those aged less than 75 years (25). Pneumonia and influenza mortality is a reflection of infection as well as risk factors for death, including the presence of comorbidities and lack of access to care. Infection data are not available for the entire country during three decades. An analysis of infection/morbidity data, however, would be a potentially interesting complement to this study.
Using a spatial autocorrelation test, we found clustering in the pneumonia and influenza percentages of states. States sharing borders frequently had similar levels of pneumonia and influenza mortality, potentially reflecting the role of large-scale environmental drivers. Local fashions in how physicians complete death certificates, interstate population movement, the local mix of circulating influenza strains and protective immunity, and stochastic effects could also contribute to spatial autocorrelation. Differences in immunity/vaccination coverage, local population density, and international and domestic traffic may have further modified the spatiotemporal patterns of influenza epidemics (7) and should be integrated into future analyses.
Western states had higher transmission season pneumonia and influenza percentages than did eastern states, even after considering background levels of pneumonia and influenza mortality and age structure. Aggregation of pneumonia and influenza mortality data to the spatial scale of state or climate region and to the temporal resolution of month potentially resulted in the loss of important detail (5). Our analyses were not conducted on fine spatial units, but our identification of broad patterns motivates further research within each heterogeneous climate region. The nine NCDC climate regions were originally defined by Karl and Koscielny (45) from statewide averages of temperature and precipitation during 1895–1981. Whether other geographic divisions may better reflect environmental determinants of influenza epidemic patterns is unknown. This grouping may have failed to capture some spatial aspects potentially relevant to influenza transmission, but the identified associations suggest that these regions are not unreasonable and that regional climates may affect pneumonia and influenza mortality.
Furthermore, subtype-specific spatiotemporal synchrony in pneumonia and influenza mortality was found. Seasons when A(H3N2) dominated were more synchronous than when A(H1N1)/B dominated. Such subtype-specific synchrony is consistent with another recent study of influenza-like illness in Japan from 1992 to 1999, which determined that new antigenic variants of A(H3N2) were associated with a monotonous spreading pattern across the country, with peak activity covering the country within 3–5 weeks (7). During seasons without new antigenic variants of A(H3N2) in Japan, multitonous patterns were observed in smaller epidemics, with peak activity taking 12–15 weeks to reach across the country (7). The interpretation of spatiotemporal patterns of influenza in the 1990s in Japan may be complicated by the concurrent cessation of their childhood vaccination program (21).
Our finding of high values of average synchrony suggests a relatively efficient spread of A(H3N2). This ease of spread might be related to a brief time to maximum infectiousness, as the peak of clinical incidence for A(H3N2) is reached more quickly than that for other subtypes (46). The lower mean in synchrony when A(H1N1) or B dominated and steep gradients in synchrony over distance could represent multiple seeding or less efficient person-to-person spread, and they might reflect a difficulty in complete mortality signal detection during mild transmission seasons. Bjørnstad et al. (39) suggest that regional synchrony can be induced by diffusion or dispersal, community processes, common seasonality, nonlinear phase-locking, large-scale climatic forcing, and/or environmental shocks. The relation between synchrony and dominant viral subtype warrants further exploration, ideally with infection or morbidity data. Future analyses would be strengthened by incorporating regional influenza virus surveillance data.
We would expect that cumulative population immunity would be lower and that A(H3N2) would be more clearly the dominant subtype in years with major A(H3N2) drift compared with minor drift. Seasons with major drift events had significantly greater average synchrony in pneumonia and influenza mortality nationwide than did those seasons with minor drift. This analysis may be limited by the data Smith et al. (43) used to determine cluster transitions. They focused exclusively on the HA1 domain, while changes in neuraminidase and internal genes may also affect the severity of an influenza season. Nevertheless, low population immunity (using major drift events as a proxy) and mortality synchrony were positively associated, a novel finding.
Our study adds to others that have analyzed influenza patterns among the United States, France, and Australia (3), across Europe (6), and within France (4), Scotland (5), and Japan (7), emphasizing strong spatiotemporal patterns in influenza-related illness/mortality. Ecologic factors underlying these patterns remain unknown, but elucidating them will improve our capacity for early warning of influenza epidemics. Recognition and explanation of persistent patterns in historical influenza mortality, in conjunction with real-time monitoring of surveillance data, could improve control strategies. These may include targeted vaccination campaigns, antiviral medication stockpiling and distribution, and travel advisories.
Our findings of marked regional differences in pneumonia and influenza mortality, of similarity in the pneumonia and influenza percentages of nearby states, and of differences that cannot be explained by the age distribution of the elderly across states support the possible effect of the environment on observed pneumonia and influenza mortality patterns. Future studies should adapt our quantitative approach to examine the possible role of climatic factors in modifying regional pneumonia and influenza mortality (10–12, 47). Our results could also inform mathematical models of influenza spread. Such models have historically focused on transportation networks (48–53), with little consideration of differences in internal dynamics within a region. Subtype-specific effects on epidemic synchrony may also be important to future modeling efforts aimed at predicting and preventing the spread of influenza.
This project was supported in part by a grant from the National Oceanic and Atmospheric Administration's Joint Program on Climate Variability and Human Health, a consortium including the Environmental Protection Agency, the National Aeronautics and Space Administration, the National Science Foundation, and the Electric Power Research Institute (NA16GP2361).
The authors thank Dr. Ottar Bjørnstad for help with assessing spatiotemporal synchrony. Dr. Pierre Goovaerts and Dr. Dunrie Greiling assisted with STIS. Dr. Jim Koopman and Dr. Cécile Viboud offered useful insights into influenza mortality analysis. Dr. Hal Morgenstern and Justin Cohen constructively reviewed a previous manuscript version.
Conflict of interest: none declared.