Influenza‐associated mortality determined from all‐cause mortality, Denmark 2010/11‐2016/17: The FluMOMO model

Background In temperate zones, all‐cause mortality exhibits a marked seasonality, and influenza represents a major cause of winter excess mortality. We present a statistical model, FluMOMO, which estimate influenza‐associated mortality from all‐cause mortality data and apply it to Danish data from 2010/11 to 2016/17. Methods We applied a multivariable time series model with all‐cause mortality as outcome, influenza activity and extreme temperatures as explanatory variables while adjusting for time trend and seasonality. Three indicators of weekly influenza activity (IA) were explored: percentage of consultations for influenza‐like illness (ILI) at primary health care, national percentage of influenza‐positive samples, and the product of ILI percentage and percentage of influenza‐positive specimens in a given week, that is, the Goldstein index. Results Independent of the choice of parameter to represent influenza activity, the estimated influenza‐associated mortality showed similar patterns with the Goldstein index being the most conservative. Over the 7 winter seasons, the median influenza‐associated mortality per 100 000 population was 17.6 (range: 0.0‐36.8), 14.1 (0.3‐31.6) and 8.3 (0.0‐25.0) for the 3 indicators, respectively, for all ages. Conclusion The FluMOMO model fitted the Danish data well and has the potential to estimate all‐cause influenza‐associated mortality in near real time and could be used as a standardised method in other countries. We recommend using the Goldstein index as the influenza activity indicator in the FluMOMO model. Further work is needed to improve the interpretation of the estimated effects.

related to Christmas and New Year holidays 1 . However, influenza is recognised as one of the main contributing factors of winter excess mortality.
A common method to assess winter excess mortality estimates, based on predefined periods with no or ignorable influenza activity, is by using the Serfling method 2,3 . The baseline for the Serfling model, which is usually fitted with a cyclic variation over the year, can be compared with the observed number of deaths, and the difference (the residual) presented as an estimate of excess mortality.
Since 2009, the network for European monitoring of excess mortality for public health action (EuroMOMO) has monitored weekly excess all-cause mortality 4 , using this method. Using this method, excess mortality may be an over-or underestimation of influenzaassociated mortality, depending on the impact of non-influenza morbidity on mortality during the winter.
To obtain specific estimates of influenza-associated mortality, parameters representing influenza activity can be included in a time series model, instead of excluding periods with influenza activity from the estimation of baseline mortality 3,[5][6][7] . In this approach, excess mortality associated with influenza is determined directly from influenza activity data and not as a residual, which will also be affected by non-influenza factors.
Within the EuroMOMO network, we developed a model to estimate influenza-associated mortality based on influenza activity in the population while controlling for extreme ambient temperatures, the FluMOMO model. This model can be used at the national level or at the European level to support risk assessment and health planning in connection with seasonal influenza epidemics and in a pandemic situation, and has the potential for timely, in-season use.
The objective of this study was to describe the statistical model, FluMOMO, and apply it to estimate influenza-associated mortality in Denmark for the seasons 2010/11 to 2016/17.

| MATERIAL S AND ME THODS
First, we describe the FluMOMO model and consider how the influenza activity and extreme ambient temperatures variables were defined. Secondly, we describe data sources, and finally analytical methods.

| The FluMOMO model
The model is a multivariable time series model with weekly number of deaths as dependent variable. We use a multiplicative Poisson regression with overdispersion and ISO-week as time unit. Influenza activity (IA) and extreme temperature (ET) are included as independent variables. The residual variation is post-regression corrected for skewness by applying a ⅔-power correction 8 A general formulation of the regression: where l is the link function (ie logarithm in a multiplicative model), X is the observed number of deaths and t is calendar week. Seasonyear (s) is defined as the period from week 27 to week 26 the following year, and p represents 4 parameters reflecting extreme temperature.
We estimate a separate effect of IA for each season and not for the entire period. The rationale for this relates to the fact that the type of circulating influenza virus and the match with the influenza vaccine will vary from one season to another.
Four parameters (ET t,1 to ET t,4 ) for ambient temperature were included, representing extreme cold in the winter (week 40 to week 20 the following year), where it may be associated with excess deaths, and extreme cold in the summer (week 21-39), where it may have a benign effect. Conversely, periods with extreme heat may be associated with excess mortality in the summer but a benign effect in the winter.
We include a lag-effect of the explanatory variables, that is, a carry-over effect of IA or ET from earlier weeks. Length of the lags is predefined as external parameters, and dIA and dET represent the delay in the formula.
The model estimates the baseline as well as the effect of IA and ET simultaneously, controlled for one another, that is, the baseline correlates to IA and ET.

| Expected number of deaths
The expected mean number of deaths including both IA and ET as explanatory variables will be The expected baseline will be Expected mean number of deaths if only IA or ET had an effect will be Hence, the expected mean excess number of deaths associated with IA t and ET t will be ( Figure 1): For post-estimation ⅔-power correction of residual variation and cumulated estimates see Appendix A.

| Influenza activity (IA)
A commonly used syndromic indicator of IA is the frequency of influenza-like illness (ILI), measured for example, through a sentinel network of general practitioners (GP). ILI can be expressed as a percentage of consultations with ILI (consultation%). If the GPs' catchment population is known, ILI can be expressed as an incidence.
Should information on ILI be unavailable, acute respiratory illness (ARI) may be used as an alternative indicator.
In Denmark, ILI is reported as the number of ILI consultations by age group, but only the total number of consultations. Hence, we do not have a denominator for each age group's consultation%.
However, assuming consultations have a constant distribution between age groups through a season, but vary between seasons, using the total number of consultations as the denominator in each age group, the consultation% in each age group will be proportional to the pattern of the correct consultation% over a specific season, thus providing valid estimates of the effect of IA, as the effect of IA is estimated separately for each season. Therefore, we used the total number of consultations as the denominator for all age groups.
Another indicator of IA is the percentage of influenza-positive samples (positive%). In Denmark, all influenza samples are registered on a personal level, that is, with known age for each person, in a national centralised database 9 . For this study, we considered more than one sample from the same person, on the same date, as one sample.
If at least one sample tested positive for influenza, we considered

| Extreme ambient temperatures
Danish weather stations register daily average temperature as well as daily minimum and maximum temperature. Averages over each of these daily temperatures, weighted by the populations in the 5 Danish regions, were used as the overall Danish daily temperature, and minimum and maximum temperature for that day. Weekly temperature as well as weekly minimum and maximum temperature was calculated as the weekly mean of each of the overall daily temperatures, respectively. Using each of these mean weekly temperatures, we estimated the expected weekly temperature, and expected minimum and maximum temperatures, in general linear models with a yearly seasonal variation, fitted as a sine wave.
Weeks with extreme temperatures were defined as weeks with a mean temperature above the expected maximum temperature or below the expected minimum temperature ( Figure 2). Nominal values were defined as how far the observed temperature was below or above the expected minimum or maximum.

| Deaths
Individual notifications of all-cause dates of deaths and dates of birth were obtained from the Danish Civil Registry.

| Influenza activity
The Danish sentinel surveillance system for influenza is where these samples have been registered since 2010.

| Ambient temperatures
Data on daily temperatures registered at Danish weather stations was downloaded from the National Oceanic and Atmospheric Administration Online Climate Data Directory 13 .
F I G U R E 1 Baseline, effect of influenza activity and extreme temperature (ET)

| Population
The size of the Danish population by age at the start of every quarter of the year was downloaded from Statistics Denmark 14 . The weekly size of each age group specific population was achieved by linear interpolation.

| Analyses
Analyses were performed separately for the age groups 0-4, 5-14, 15-64 and 65+ years of age, as well as for all ages.
A fatal outcome after infection with influenza virus will normally occur after a time-delay, where co-morbidity and complications may contribute. On this basis, we applied a 2 weeks lag for both IA and ET.
We included a trend and a half-yearly cycle if they contributed on a 5% level (P < .05), while yearly seasonality was included on a 10% level.
There may be other changes in mortality over calendar time than seasonality. Therefore, a trend is included in the baseline. Long-term changes may be modelled using a non-linear trend, for example, a spline. However, non-linear trends tend to be less stable at the ends, and if we intend to use the model for monitoring, interest will be in the ends, especially the right end, which informs the most recent impact. As the model potentially can be used for in-season and near timely monitoring, we have limited the model to use a linear trend.
When estimating over more than 5 seasons, the linear trend assumption did not always hold (likelihood-ratio test between a spline and a linear trend), and including <3 seasons seems to make the fit less     Table 1 and illustrated in Figure 5.

| D ISCUSS I ON
Infection with the influenza virus may often lead to exacerbation of underlying chronic conditions or to secondary bacterial infections 15 .
Therefore, influenza is often not recorded as the primary cause of death, and using cause-specific data as an outcome may lead to underestimation of influenza-associated mortality. Furthermore, in most countries, cause of death is available with large delays.
Therefore, timely estimation of influenza-associated mortality should be based on all-cause mortality, although this outcome may be associated with loss of some specificity. The present work started as an initiative within the EuroMOMO network, with the aim of obtaining a direct and timely estimation of mortality associated with influenza.
In the present study based on Danish data, the model generally captures the background seasonal mortality pattern well. However, Age: All ages Other respiratory pathogens are not directly included in the model, even though they may affect winter mortality. If they should be included in the model, care has to be taken to avoid collinearity, in particular when lags are included, as a lagged effect of one pathogen may be confounded by a direct impact of another pathogen 5,15 .
Additionally, if included in the model, potential interactions may complicate the model and blur the estimates. However, this should not prevent further investigation.
Another factor of importance in the seasonal variation in mortality is ambient temperature or more precisely deviation from expected temperature, for example, cold snaps in the winter or heat waves in the summer 17,18 . We have controlled influenza-associated mortality for temperature, but the pattern of the effect of temperature on F I G U R E 5 Influenza activity associated winter (week 40 to week 20 the following year) mortality rates using the FluMOMO model mortality has to be investigated further to be able to report valid estimates of this effect.
We estimated Danish influenza-associated mortality using 3 indicators of influenza activity (Table 1 and Figure 5). All 3 indicators show a similar pattern. ILI possibly overestimates influenza mortality as ILI may be caused by other respiratory pathogens. Using posi-tive% as indicator for IA will reduce the effect of other respiratory pathogens, but may not, as ILI, reflect the dynamic in the population.
A combination of the 2, the Goldstein index, will reflect both the dynamic of transmission and limit overestimation due to other pathogens 10,11 . Mortality associated with influenza based on the Goldstein index seems to be the most appropriate for Danish data, primarily because it is conservative. Using the Goldstein index, we found a median influenza-associated mortality rate of 8.2 per 100 000 population which is substantially lower-and probably more realistic-than a previous applied methodology that obtained a median mortality of season, although raised among children below 5 years of age 24  We also acknowledge the support provided by the World Health Organization, WHO, regional office for Europe, and the European Centre for Disease Prevention and Control, ECDC (Framework contract for service-ECDC/2016/041).

CO N FLI C T O F I NTE R E S T S
The authors declare that they have no competing interests.

APPENDIX A P OS T-E S TI M ATI O N ⅔-P OWE R CO R R EC TI O N O F R E-S I D UA L VA R I A N CE
Formulation of the ⅔-power variance corrected residual variance (RVar), that is, the uncertainty or confidence intervals of the estimated mean effects: Using the delta-method, the residual ⅔-power variances of ∆IA t or ∆ET t becomes where E(B t ) is the estimated baseline number of deaths and E(IA t ) is the expected number of deaths had the effect of ET t been zero.
Var(E(B t )) and Var(E(IA t )) are the regression variances. However, in a regression with a link function l, these are "only" known at the lscale. Using the link function log, as we do, the regression variance become where Var log (E(X t )) is the regression variance on the log-scale. Hence and likewise for ET t : Confidence intervals are then calculated as: where X t may be IA or ET, respectively.

CU M U L ATE D N U M B E R O F E XCE SS D E ATH S
The mean cumulated excess number of deaths over s weeks will be Assuming independence between weeks, that is, no residual autocorrelation, then the ⅔-power variance corrected residual variance of the cumulated excess deaths over s weeks will be Confidence intervals are then calculated as: where EB is baseline and X t may be IA or ET, respectively.