Mortality attributable to seasonal influenza in Greece, 2013 to 2017: variation by type/subtype and age, and a possible harvesting effect

Introduction Estimating the contribution of influenza to excess mortality in the population presents substantial methodological challenges. Aim In a modelling study we combined environmental, epidemiological and laboratory surveillance data to estimate influenza-attributable mortality in Greece, over four seasons (2013/14 to 2016/17), specifically addressing the lag dimension and the confounding effect of temperature. Methods Associations of influenza type/subtype-specific incidence proxies and of daily mean temperature with mortality were estimated with a distributed-lag nonlinear model with 30 days of maximum lag, separately by age group (all ages, 15–64 and ≥ 65 years old). Total and weekly deaths attributable to influenza and cold temperatures were calculated. Results Overall influenza-attributable mortality was 23.6 deaths per 100,000 population per year (95% confidence interval (CI): 17.8 to 29.2), and varied greatly between seasons, by influenza type/subtype and by age group, with the vast majority occurring in persons aged ≥ 65 years. Most deaths were attributable to A(H3N2), followed by influenza B. During periods of A(H1N1)pdm09 circulation, weekly attributable mortality to this subtype among people ≥ 65 years old increased rapidly at first, but then fell to zero and even negative, suggesting a mortality displacement (harvesting) effect. Mortality attributable to cold temperatures was much higher than that attributable to influenza. Conclusions Studies of influenza-attributable mortality need to consider distributed-lag effects, stratify by age group and adjust both for circulating influenza virus types/subtypes and daily mean temperatures, in order to produce reliable estimates. Our approach addresses these issues, is readily applicable in the context of influenza surveillance, and can be useful for other countries.


Introduction
Seasonal influenza is a considerable public health problem worldwide. A 2018 modelling study estimated that 291,000-646,000 annual deaths globally are attributable to influenza, most of those in older adults [1]. The impact of influenza on mortality can vary widely each year and across different populations, due to factors such as the types/subtypes of circulating viruses, the age structure of the population, its level of health and access to healthcare, and its degree of immunity due to vaccination or previous influenza epidemics [2,3].
There is intense interest in obtaining accurate and timely estimates of influenza-attributable mortality across all age groups, in order to inform public health policy and evaluate the effectiveness of preventive interventions. Numerous studies have attempted to assess influenza-attributable mortality, mostly in highincome countries, and with heterogeneous results [4]. Differences in study methodology are important contributors to this heterogeneity, besides actual hostand pathogen-related factors [4].
The variety in the employed methods is a reflection of the many challenges involved in accurately estimating the part of population mortality that is attributable to influenza. First, the types/subtypes of circulating influenza viruses need to be taken into account in any statistical model, as not all types/subtypes may have the same effects on all age groups [5]. This highlights the need to have solid laboratory data to inform an influenza activity proxy, as models without such a proxy tend to produce higher mortality estimates [4]. Second, estimation needs to adjust for the confounding effect of ambient temperature, not only at the extremes but also at milder cold levels, which may account for substantially more temperature-related deaths than extreme cold [6]. Finally, the association between influenza circulation and mortality is characterised by a lag that cannot be assumed to be fixed; not all people die at the same time after falling ill with influenza, therefore the effect of a given level of influenza activity on mortality is likely to be distributed over a period of time. Such distributed-lag associations with mortality have been previously described for temperature [6] and air pollution [7], but are also likely for influenza as well. To our knowledge though, very few studies to date have attempted to account for the lag dimension in the association between influenza and mortality [8][9][10][11], despite the availability of powerful relevant methods [12].
The objective of the current study was to examine the mortality attributable to influenza over four winter seasons in Greece, and across different age groups, while addressing the methodological issues outlined above. For this purpose we combined epidemiological, virological and mortality data routinely collected as part of public health surveillance, with temperature data covering the entire country. A secondary objective was to compare the mortality attributable to cold temperatures with that attributable to influenza, in order to assess which one is the primary driver of excess winter deaths [13,14].

Data sources
The study encompassed the entire population of Greece, a high-income country with a temperate climate and a population of slightly less than 11 million. Individual notifications of all-cause deaths from International Organization for Standardization (ISO) week 22/2013 (May 2013) to week 41/2017 (October 2017) were available and obtained in December 2017 from the Hellenic Ministry of the Interior. In this dataset 99.3% of all deaths were registered within 1 week of their occurrence, thus there is practically no bias due to reporting delay. Deaths coded as violent (0.5% of the total) were excluded. We summed the notifications by age to create daily counts of deaths for three groups: all ages, 65 years or older, and 15-64 years.
In the 0-14 years age group the data were too sparse (median: 2 deaths/day; range 0-19) and were not analysed further. Population denominators by age group were also obtained from the Hellenic Statistical Authority (http://www.statistics.gr/).
Average daily temperatures for the same period recorded at seven weather stations in Greece (Supplement S1 Figure 1) were downloaded from the website of the National Centers for Environmental Information (NCEI) of the National Oceanic and Atmospheric Administration (NOAA) [15]. A populationweighted daily mean over all stations was used as the countrywide temperature of that day.
Weekly influenza-like illness (ILI) rates for the same period were obtained from the Hellenic Centre for Disease Control and Prevention (HCDCP). These rates are recorded by a sentinel surveillance system operated by the HCDCP, based on a countrywide geographically representative network of primary care physicians who voluntarily report every week their total number of patient consultations and the number of patients with ILI. The ILI rate is calculated using a mixed-effects Poisson model (Supplement S1 Figure  7). As primary care physicians in Greece do not have a defined catchment area, the ILI rate denominator is total consultations rather than population; this ILI rate is linearly related to ILI incidence, assuming a constant average consultation rate per person [16]. Aggregate laboratory data were also obtained by the two National Influenza Reference Laboratories operating in Greece (Hellenic Pasteur Institute, Athens, and Microbiology Laboratory, Aristotle University of Thessaloniki School of Medicine, Thessaloniki), namely the weekly number of all respiratory specimens tested by real-time PCR (RT-PCR) for influenza, and the number of specimens positive for subtype A(H1N1)pdm09 and A(H3N2) viruses respectively, as well as those positive for type B influenza virus. These two laboratories perform the bulk of influenza RT-PCR testing for surveillance purposes in Greece.
To create incidence proxies for type/subtype-specific influenza activity, we multiplied the weekly ILI rates by the proportion of specimens testing positive to each of the three influenza virus types/subtypes, as previously suggested in the literature [17]. Under certain assumptions, these proxies are approximately linearly related to the weekly type/subtype-specific incidence of influenza in the population [18]. Each weekly timeseries of influenza incidence proxies was converted to a daily time-series by assigning the weekly values to all days of each week and applying a cubic smoothing spline with one degree of freedom per week, with the smoothing parameter chosen using the Generalised Cross-Validation (GCV) criterion. This method introduces minimal smoothing to the daily time series, essentially only removing the jags between adjacent weeks (Supplement S1 Figure 2).

Statistical analysis
The association between daily death counts for each age group, average daily temperature and the three influenza incidence proxies was modelled using distributed-lag nonlinear models (DLNM) of quasi-Poisson family and with log link [12]. In this type of models, both the 'conventional' exposure-response and the additional lag-response association are simultaneously assessed using a pair of functions known as the crossbasis functions, which transform the exposure variable to a cross-basis matrix. The regression coefficients for this matrix, expressed as a relative risk (RR), then describe the relationship of the exposure to the outcome across both dimensions [12]. The general shape of both the exposure-response and the lag-response association needs to be described by specifying appropriate cross-basis functions, which can be nonlinear. The maximum lag also needs to be specified.
The exposure-response curve for temperature was modelled with a quadratic B-spline with three internal knots placed at the 10th, 75th and 90th percentile of the temperatures distribution, similar to previous studies [6]. For the three influenza incidence proxies a linear relationship with mortality was specified; this implies a fixed case-fatality ratio per influenza type/ subtype, i.e. that mortality is proportional to the incidence of infection with a particular virus type/subtype. We modelled the lag-response association for all four variables using a natural cubic spline with three internal knots placed at equally spaced values on the log scale [6]. The maximum lag was set to 30 days, in order to capture the delayed effects of both influenza and cold temperatures, but also in order to assess potential mortality displacement (harvesting) over this period, i.e. detect deaths that would occur anyway and were only brought forward in time by the exposure, rather than being 'genuine' excess deaths that would not have occurred without the exposure [19]. The assumptions regarding the lag-response shape and the maximum lag were tested in a sensitivity analysis (Supplement S1 Figures 3-4).
Besides the four cross-basis matrices for temperature and influenza, we included in the regression models a linear term for trend, an indicator variable for day of the week, and a periodic cubic B-spline with three equidistant knots for day of the year (1-366) to model seasonality; the latter allows a more flexibly-shaped seasonal pattern of mortality compared with a classical Serfling-like sinusoidal term. One such model was fit on the daily death counts for every age group (all ages, 15-64 years, 65 years or older).
For every day in the series, the number of deaths attributable to each influenza type/subtype and to cold temperatures was calculated using a previously described method for DLNMs [20]. Specifically, the attributable risk at each day was treated as being caused by multiple lagged effects of the exposure in the previous days up to the maximum lag back ('backward' attributable risk). The total attributable number of deaths is then given by summing the contributions of each day. Empirical 95% confidence intervals (95% CI) are calculated using Monte Carlo simulations, assuming a multivariate normal distribution of the cross-basis regression coefficients [20]. We used 50,000 replicates for that purpose, to ensure precision. We calculated influenza-and temperature-attributable mortality for each of the four winter seasons available in our dataset (defined as week 40 of each year to week 20 of next year, from 2013/14 to 2016/17), and also for every week in the study period to create an epidemic curve. In this context, 'cold temperatures' means lower than the minimum mortality temperature, which is the point where mortality is lowest; however, our model can estimate the mortality attributable to any temperature range, be it mild cold or extreme cold.
All calculations were performed with the R software environment, version 3.4.3 [21], and with package 'dlnm' for R, version 2.3.2 [22].

Ethical statement
This study used only aggregate and non-identifiable surveillance data, therefore no ethical approval was necessary.

Results
From The overall association (across all lags up to the maximum of 30 days) of temperature with mortality, as provided by the fitted DLNMs, is illustrated in Figure  2. It was found to be slightly attenuated for those aged 15-64 years, although the difference was nonsignificant. The temperature of minimum mortality was 25 °C, for both age groups. With respect to the three influenza types/subtypes, the overall RR is a linear function of the respective incidence proxies and it is illustrated in Table 1    In comparison, cold temperatures (lower than the minimum mortality temperature) were associated with 74.7 excess deaths per 100,000 population (95% CI: 35.3 to 111.7), and accounted for the larger part of the wintertime seasonal excess mortality ( Figure 3). Most of the influenza-attributable mortality was accounted for by subtype A(H3N2) affecting older people, with type B also making a substantial contribution. In contrast, even in seasons with high A(H1N1)pdm09 activity, influenza A(H1N1)pdm09 was not associated with a net excess of deaths among people 65 years and older (Table 2).
However, when plotting weekly influenza-attributable mortality by type/subtype for the entire study period (Figure 4), an interesting pattern emerged: A(H1N1) pdm09 activity appeared to be associated with a substantial excess of influenza-attributable deaths for several weeks, but towards the end of the seasonal epidemic wave the number of associated deaths fell to zero and even became negative (for season 2015/16). Even though confidence bands were wide, this is suggestive of a possible mortality displacement (harvesting) effect for influenza A(H1N1)pdm09. Stratifying weekly A(H1N1)pdm09-attributable mortality by age group shows this pattern only among those 65 years or older, and not among people aged 15-64 years Similar results were obtained in sensitivity analyses under alternative model specifications (Supplement S1 Figures 3-4). A longer maximum lag appeared to increase the number of deaths attributable to influenza A(H3N2) and augment the observed mortality displacement for A(H1N1)pdm09, while a shorter maximum lag appeared to attenuate it. Interestingly, a model based on weekly data generated attributable mortality estimates similar to those of the main model. Finally, a model without the cross-basis matrix for temperature resulted in twice as many deaths attributable to influenza (all types) compared with the main model (4,552 vs 2,559 per year -Supplement S1 Figures 3 and 5).

Discussion
Our results confirm that the mortality attributable to influenza varies greatly between seasons, and indicate that the types/subtypes of circulating viruses are a major determinant of this variation, with A(H3N2) virus being associated with many more deaths compared with A(H1N1)pdm09 and influenza B virus. As such, it is essential to individually account for each influenza virus type/subtype when modelling their association with mortality, using detailed laboratory data, rather than consider all influenza together in a single incidence proxy. We have also confirmed that the majority of influenza-attributable deaths occur among persons 65 years and older, and shown that each influenza virus type/subtype has different effects on mortality among different age groups, with A(H3N2) affecting older adults, and A(H1N1)pdm09 younger persons. Therefore, in order to fully explore the effects of influenza on population mortality, stratification by age group should always be performed. Although it is difficult to compare studies due to different methodologies, our results from Greece are overall very similar to analogous estimates from the United States, with respect to both the total mortality burden attributable to influenza and the relative contribution of each influenza type/subtype [17]. They are also within the range reported in a recent systematic review of studies examining influenza-attributable mortality [4] ( Table 3).
The addition of the lag dimension in the association between influenza and all-cause mortality offers a novel and important insight. To our knowledge this is the first time a short-term harvesting effect is suggested for influenza, for a particular type/sybtype and in a particular age group. Such an effect could not have been detected by previous studies that did not employ distributed-lag modelling, and only considered fixed lags selected according to goodness of model fit [4]. One prior study from Sweden did observe that high winter mortality, including deaths coded as influenza, led to fewer heat-related deaths over the following summer [23]. Our study, however, is the first to directly associate with influenza A(H1N1)pdm09 an excess of deaths among people 65 years or older, over several weeks, followed by a period with no excess or even a possible mortality deficit for several more weeks within the same season. Due to low precision, this finding will need to be confirmed in further studies applying similar distributed-lag models in bigger datasets from larger populations (and in finer age groups, such as 65-74 years and 75 years and above). If confirmed, it would suggest that among the elderly influenza A(H1N1)pdm09 affects those with the most vulnerable health, many of whom would have died in the short-tomedium term from other causes independent of influenza. In contrast, A(H3N2) and type B appear to be associated with a 'genuine' excess of influenza-related deaths, among older people who would not otherwise have died in the same year. This has important implications for correctly assessing the effect of circulating influenza virus types/subtypes on population mortality, as well as the potential benefit of targeted vaccination campaigns.
There is an ongoing debate about whether the wintertime seasonal mortality excess is primarily driven by seasonal influenza [13,24] or cold temperatures [14,25]. In our study we found that mortality attributable to cold temperatures was several times higher

Figure 3
Weekly observed all-cause deaths, deaths attributable to influenza and to cold temperatures, Greece, May 2013-October 2017 (n = 518,688 deaths) Week than that attributable to influenza, with temperature explaining most of the within-year seasonal variation in mortality (Supplement S1 Figure 5). The strong confounding effect of temperature on the association between influenza and mortality is highlighted by our sensitivity analysis, which found a doubling of influenza-attributable deaths when not adjusting for temperature. As a result, to minimise residual confounding we specified the same maximum lag and lag-response specification for both influenza and temperature, and used daily counts of death and daily temperatures rather than weekly averages. On the other hand, long-term temperature trends are themselves highly seasonal (Figure 1), therefore it has been argued that distributed-lag models with long lags may capture part of this seasonality and produce larger effects of cold temperatures due to collinearity [24]. However, in our analyses we did not detect substantial collinearity (Variance Inflation Factors for all coefficients of interest was under 1.7), while in sensitivity analyses models with shorter or longer lag produced similar attributable mortality estimates for both influenza and temperature (Supplement S1 Figure 3). Although it is true (and our results also indicate) that the magnitude of the overall winter excess mortality does correlate with the dominant influenza type/subtype for each season [13], influenza-attributable deaths appear to be a minority among all winter deaths. The majority is attributable to cold temperatures, independent of influenza, and that part is relatively constant each season.

Figure 4
Weekly numbers of estimated influenza-attributable deaths for all ages, by influenza type/subtype (panels A, B, C) during four winter seasons in Greece, weeks 40/2013-20/2017

A. Influenza A(H1N1)pdm09
Week number/year Week number/year Week number/year Weekly number of influenza-attributable deaths Blue bars represent weekly numbers of influenza-attributable deaths; the grey areas are 95% confidence bands.
Our proposed approach uses routine epidemiological and laboratory surveillance data, provides detailed weekly estimates of influenza-attributable mortality, and can be readily implemented with open-source software. As a result, it is particularly useful and applicable in an influenza surveillance context [26]. Further research and evaluation in settings other than Greece is warranted in this regard, as well as a comparison with other methods to estimate influenza-attributable mortality that are presently used in Europe [27,28].
Besides the advanced DLNM statistical methodology, a strength of the current study is the integration of comprehensive epidemiological, laboratory and weather data in a single detailed analytical model, offering more reliable estimates of influenza-attributable mortality. Thus we could account for individual influenza types/subtypes, closely adjust for daily mean temperatures, flexibly model seasonality using periodic splines, and allow for distributed-lag effects. We also used daily death counts as the outcome, rather than the usual weekly counts, in order to minimise residual confounding.
On the other hand, the study has important limitations.
With only four seasons of available mortality data on a relatively small population, statistical power to detect associations was suboptimal particularly for the younger age groups. Especially in the highly interesting 0-4 year group, deaths were so sparse that no informative conclusions could be drawn. In addition, we had only one season with significant influenza B circulation (2014/15, entirely of B/Yamagata lineage), which might not necessarily be representative of the actual effects on mortality of that influenza type, or of the B/ Victoria lineage. We could only study all-cause deaths, as no comprehensive coded cause-specific mortality data are available in Greece; for example, it would have been very useful to examine cardiac or 'pneumonia and influenza' (P and I) deaths, but unfortunately this was not possible. We also had no available laboratory surveillance data on respiratory syncytial virus (RSV) infection; RSV is known to account for a substantial part of winter excess deaths [10,11], and might co-circulate with influenza virus [29], thereby confounding its association with mortality. Finally, we did not specifically consider the effect of seasonal influenza vaccination in the population; however, vaccination coverage in Greece is reportedly very low even among high risk groups [30], and is also unlikely to vary significantly between seasons.
In conclusion, we have analysed the mortality burden attributable to influenza in Greece, and shown it to be Coloured bars (in pink or orange) represent weekly numbers of influenza-attributable deaths; the grey areas: are 95% confidence bands.
substantial especially among older people, as well as highly variable depending on the circulating influenza types/subtypes; therefore stratification by age and adjustment for influenza types/subtypes and ambient temperature are essential to obtain accurate estimates of influenza-attributable mortality. We have also shown how the DLNM framework can be used to explore the lag dimension, and identify potential mortality displacement in the association between influenza and mortality. More studies with similar methodology are warranted, to confirm these results, further investigate the observed harvesting effect, and evaluate the utility of this approach in the context of routine influenza surveillance; our methods are readily applicable, and we hope will find wide use in other countries as well. A final takeaway from this study: if the majority of influenza-attributable deaths occur in the age group (65 years and older) that has the lowest immune response to the seasonal influenza vaccine [31], and by the subtype (A(H3N2)) that the vaccine usually offers the least protection against [3], then improved and more effective influenza vaccines are urgently needed to protect the population, with a view to the ultimate goal of a 'universal' influenza vaccine [32].