An exploratory study to assess patterns of influenza- and pneumonia-related mortality among the Italian elderly

ABSTRACT

Older adults are at disproportionately high risk of severe influenza-related outcomes and represent the main target of the annual influenza vaccination. The protective effect of seasonal influenza vaccination on the observed mortality indicators is controversial. In this ecological study, spatiotemporal patterns of pneumonia- and influenza-related mortality registered in the Italian elderly over seven (2011–2017) consecutive seasons were explored and the epidemiological association between the observed local pneumonia- and influenza-related mortality and influenza vaccination campaign features were modeled by using both fixed- and random-effects panel regression models. The descriptive spatiotemporal analysis showed a clear North–South gradient, where northern regions tended to report more pneumonia- and influenza-related deaths. After adjustment for potential confounders, it was found that each 1% increase in influenza vaccination coverage rate would be associated (P < .001) with a 1.6–1.9% decrease in pneumonia- and influenza-related mortality. Moreover, each 1% increase in the use of MF59®-adjuvanted trivalent influenza vaccine would be associated (P < .05) with a further 0.4% decrease in pneumonia- and influenza-related mortality. This study supports the increase in annual influenza vaccination in Italy and suggests that a higher level of use of the adjuvanted influenza vaccine in the elderly may be beneficial.

KEYWORDS:

• Influenza
• vaccine
• adjuvanted influenza vaccine
• elderly
• ecological study
• spatial analysis
• Italy

Introduction

Worldwide, influenza is one of the leading infectious disease in terms of both incidence and mortality rates.^1,2 Seasonal influenza vaccination represents the most effective public health intervention able to reduce the burden of disease.^2,3 Indeed, the World Health Organization’s (WHO) most recent position paper² has listed several priority targets for annual influenza vaccination: pregnant women, children aged 6 months to 5 years, the elderly, subjects with specific chronic conditions, healthcare workers, and international travelers. Among these, the elderly is probably the most recognized target group; indeed, as per the European Center for Disease Prevention and Control (ECDC),⁴ all European Union (EU) Member States recommend seasonal influenza vaccination for older adults.

Influenza vaccine (IV)-induced immunogenicity and/or protection is often poor in the elderly as a result of immunosenescence.⁵ In order to circumvent this unmet need, alternative IV formulations have been developed. The first worldwide available IV specifically developed for the elderly was that formulated as a standard-dose egg-based subunit trivalent IV including MF59® (Seqirus UK Ltd.) adjuvant (adjuvanted trivalent influenza vaccine; aTIV).⁶ Italy was the first country to adopt aTIV in 1997,⁶ where it was still available during 2020/21 influenza season. Other historically or currently commercialized IVs may also address immunosenescence, including (i) virosomal;⁷ (ii) intradermal⁸ and (iii) high-dose IVs.⁹

The rationale for this study was primarily driven by a gap in the understanding of association between influenza vaccination coverage (IVC) rates and influenza-associated mortality. For instance, a previous Italian ecological study¹⁰ did not find any meaningful association between the IVC rate and influenza excess mortality over time. However, the study by Rizzo et al.¹⁰ did not distinguish between different types of available IVs at that time. Indeed, in the paper by Rizzo et al.,¹⁰ dating back to 2006, it has been stated that “In Italy, more immunogenic vaccine with novel adjuvants has been introduced since 1997 … but it is too early to evaluate their population impact.” On the other hand, a more recent study¹¹ conducted in the Province of Treviso (northeastern Italy) found that the risk of all-cause death was significantly lower (by 33–39%) in the vaccinated elderly (as compared with unvaccinated subjects) in three consecutive seasons (2014/15–2016/17). Of note, aTIV was the most frequently administered IV in Treviso.¹¹

In the context of Italian fiscal federalism single regions, the autonomous provinces of South Tyrol and Trento (henceforth referred to as “regions”) are granted a certain level of freedom to achieve their own public health goals.¹² Regarding influenza immunization, each year the Italian Ministry of Health issues a circular on the prevention and control of influenza;¹³ each region may then fully adopt the national recommendations or provide its own circular/recommendations.^12,14

Diversity in the adopted policies may result in both (i) the so-called “jeopardization”^12,14 of IVC rates [up to double difference reported in the observed 2019/20 season IVC rates among older adults aged ≥65 years¹⁵ and (ii) different patterns of use for the available types of IV.^12,14

The primary aim of this study was to explore spatiotemporal patterns of pneumonia- and influenza (P&I)-related mortality observed in Italian older adults aged 65 years or above. The second goal was to investigate the epidemiological association between the observed local P&I mortality among subjects aged ≥65 years and IVC rates and IV policy patterns.

Materials and methods

Overall study design

This is a typical ecological study: we investigated officially registered P&I-related mortality in the elderly (defined here as subjects aged ≥65 years) at the level (i.e., unit) of Italian provinces (N = 110) and/or regions (N = 21) over seven consecutive post-pandemic seasons (2010/11–2016/17). In other words, we analyzed population groups and not single individuals. Both exploratory and analytical approaches were considered.

Readers interested in both strengths and limitations of the ecological study design are invited to read the paper by Morgenstern;¹⁶ moreover, the limitations specific to this exploratory study will be discussed later in the manuscript.

In this paper, we considered only the post-pandemic period (i.e., starting from season 2010/11). This choice was based on the fact that the pandemic 2009 A/H1N1 (A/H1N1pdm09) virus completely replaced the so-called seasonal A/H1N1 (A/H1N1s) that circulated before 2009.^17,18 Data for 2018 onwards were not considered since no officially reported P&I-related mortality estimates were available at the time of data extraction (as of December 2020).^19,20

Data sources

Most data came from the official Italian data flows publicly available from the Italian Ministry of Health 15, National Institute of Health,¹⁸ National Institute of Statistics,^19,20 and National Institute for Environmental Protection and Research.²¹ Data regarding quotas for different types of IV were provided by Seqirus S.r.l., Italy (company database of regional demands for individual types of IV, i.e., data on tender allotments). The variables considered and corresponding data sources are reported in Supplementary Material, Table S1.

Study outcome

The study outcome was the country-/province-/region- and year-specific estimate of P&I-related mortality in older adults aged ≥65 years as per the European Shortlist for Causes of Death (N = 65 causes) compatible with the three most recent International Classification of Diseases (ICD) versions for influenza (ICD-8: 470–474; ICD-9: 487; ICD-10: J10–J11) and pneumonia (ICD-8: 480–486; ICD-9: 480–486; ICD-10: J12–J18) codes.^20,22 For this reason, we extracted the readily available dataset²⁰ on the P&I mortality rate (per 10,000 inhabitants) until the last available year of 2017 for the whole country, regions, and provinces.

Spatiotemporal analysis of pneumonia- and influenza-related mortality

Depending on data availability,²⁰ the spatiotemporal analysis could be conducted at the level of single provinces (N = 110). For this reason, first we visually explored the observed province-specific P&I mortality rates separately by year. This was done by plotting choropleth maps.

Moran’s I global spatial autocorrelation coefficients²³ were then computed in order to measure the overall clustering pattern of the observed mortality rates. The interpretation of Moran’s I is similar to that of Pearson’s r correlation coefficient: positive statistically significant I values indicate geographic patterns of spatial clustering, negative significant I estimates show clustering of dissimilar values, while non-significant values at α < 0.05 indicate complete spatial randomness. Considering that Italy has the two Islands of Sicily and Sardinia, for Moran’s I statistics the k nearest neighbor spatial weights matrix was used.²⁴ As a “rule-of-thumb”²⁵ we set the value of k as the square root of the total number of observations (N = 110); this means the k-value used was 10.

Providing that all and year-specific global I coefficients were statistically significant, we then further investigated the local indicators of spatial association (LISA).²⁶ Choropleth maps were created to visualize the four types of clusters/outliers, namely hot-hot (hotspots), i.e., observations that signified provinces with a higher than average mortality rate surrounded by provinces with a higher than average mortality rate, while cold–cold (coldspots) signified provinces with a lower than average mortality rate surrounded by provinces with a lower than average mortality rate. Low–high and high–low outcomes represented outliers: these were provinces with low/high average mortality rates surrounded by provinces with high/low average mortality rates, respectively.

The spatiotemporal analysis was performed in R stats packages, version 2.15.2.²⁷

Spatiotemporal analysis of pneumonia- and influenza-related mortality

The independent variables of interest were regional IVC rates and the proportion of aTIV doses to the total number of IV doses put into tender allotments. During the study period, IVC was recommended and fully reimbursed for all subjects aged ≥65 years, people ≥6 months affected by certain health conditions and some other categories. The Italian Ministry of Health routinely report region- and season-specific IVC rates for both the general population and older adults aged ≥65 years.¹⁵ Data on province-specific IVC rates are not publicly available. Therefore, the unit of this analytical part of the analysis was a region (N = 21). The primary predictor of interest was IVC in older adults aged ≥65 years. However, a higher IVC rate in younger age groups may exercise some protective effect on the elderly owing to the phenomenon of herd protection. Indeed, some studies underlined the important role of children and adolescents in spreading influenza virus in their households^28–30 and therefore to their grandparents. A model by Fumanelli et al.³¹ has suggested a significant social interaction between Italian elderly and younger individuals. For this reason and in order to account for the possible effects of herd protection, we also included a variable of IVC in subjects aged <65 years.

Another independent variable of interest was the proportion of potential aTIV users to the total number of IV doses put into tender allotments, and the data from the Seqirus Italy tender department. According to the latest Italian official recommendations,¹³ aTIV may be used only for people aged ≥65 years. Therefore, we hypothesized that the higher local use of aTIV may be associated with better health-related outcomes among the Italian elderly.

To establish an association (or lack of association) between the region- and year-specific P&I mortality rates in the elderly and predictors of interest (i.e., IVC rates and share of aTIV doses) panel regression analysis was undertaken. Briefly, the panel considered 21 spatial units (i.e., regions) followed over seven consecutive post-pandemic years and therefore consisted of 21 × 7 = 147 observations. However, an important assumption has to be highlighted here. P&I mortality data are routinely reported for the whole calendar year,^19,20 while the IV campaign usually starts in mid-October/November and almost all IV doses are administered by the end of December³². In Italy, according to the Italian National Institute of Health, most laboratory-confirmed influenza deaths occur between January and March,³² and influenza-like illness (ILI) peaks were usually reached in late January or February (Supplementary Material, Table S2).¹⁸ Moreover, considering the time lag of 2–6 weeks between IV administration and the peak of the vaccine-induced immune response,³³ it is more likely that IV administered in autumn/winter of a year t-1 will mainly exercise its effect (if any) on bacterial influenza-related complications (that are the most frequent^35,36 and require some time to be developed) leading to death in the first months of the following year t.

Both the fixed-effects (FE) and random-effects (RE) methods were applied. The FE approach may be useful in the context of causal inference: while standard regression techniques provide biased estimates of causal effects in case there are unobserved confounders, FE regression may provide unbiased estimates in this situation.^32,34 In other words, in our models, region-level FEs were included to absorb unobserved region-level heterogeneity in the observed P&I mortality rates not explained by other covariates in the model.³⁷ Indeed, unobserved effects are typical in ecological and social research.³⁸ By contrast, the RE approach assumes that region-specific effects are not correlated with independent variables.³⁹ In any case, the Hausman’s specification test was applied⁴⁰ to formally differentiate between FE and RE models; the null hypothesis of this test is that the RE model estimates are consistent and efficient.

The following socioeconomic, environmental, and virological variables were selected as potential confounders: public health expenditure per capita (€), population density (inhabitants per km²), average winter temperature, and the predominant influenza virus (sub)type(s). The reasons for inclusion of these variables are described below.

Public health expenditure per capita represents a proxy measure of regional welfare and is commonly used in health-related econometric studies.^41–44 Indeed, this parameter varies substantially among the Italian regions^20,44 and has been found to be a significant predictor of regional measles, mumps, and rubella (MMR) vaccination uptake in Italy.⁴⁴

As per environmental factors, we selected two potential confounders, namely: population density and mean winter temperature regimens. The empirical idea for the former variable was that a higher population density would be associated with a higher virus transmission.^37,45 In fact, the population density in Italy is highly non-homogeneous.²⁰ Second, single Italian regions lay in different climatological areas with highly different daily temperature paradigms; this fact could have direct implications on the influenza-related outcomes since the so-called “cold waves” usually interfere with the mortality rate.⁴⁶ Moreover, Lytras et al.⁴⁷ have concluded that in Greece the winter excess mortality rates attributable to cold temperatures were substantially higher than those attributable to influenza. Therefore, we proxied the cold waves in a given year and region as an average minimum temperature observed in the winter period. In our analysis, the winter period started at week 40 of the previous year and ended at week 20 of the next year, as per the FluMOMO model.⁴⁸

Finally, circulation patterns of influenza virus (sub)types (A/H1N1pdm09, A/H3N2 and B) may determine the magnitude of influenza-related outcomes. For instance, in Italy the predominance of the A/H3N2 subtype was associated with significantly higher excess mortality in the elderly.¹⁰ The predominance of a single virus (sub)type over other (sub)types was a priori set to 50% of the total national detections. This assumption was however, formally proved by performing a single-proportion z-test. Moreover, the adopted classification rule was compared with the previously published Italian studies,^49,50 meeting full agreement. Otherwise [i.e., when the most prevalent virus (sub)type was detected in <50% cases], the overall virological picture was dubbed as co-circulation (Supplementary Table S2).

As recommended,^37,51 in all panel regression models performed, the continuous variables (i.e., P&I mortality rates, public health expenditure per capita, population density, and average winter temperature) that were not percentages were transformed using natural logarithms (loge). The regression coefficients are therefore interpreted as elasticities. For instance, the model coefficient for IVC rate should be interpreted as the percent change in P&I mortality rate associated with a 1% change in coverage.^37,51

The following panel model specification was considered:

$\begin{array}{rl}& \mathit{lo}{g}_e{\left(\mathit{PI}\mathrm{\_}\mathrm{mort}\mathrm{\_}65+\right)}_{\mathit{i},\mathit{t}}=b_1{\left(\mathrm{IVC}\mathrm{\_}65+\right)}_{\mathit{i},\mathit{t}}\\ & +b_2{\left(\mathrm{IVC}\mathrm{\_}65\right)}_{\mathit{i},\mathit{t}}+b_3{\left(\mathrm{aTIV}\right)}_{\mathit{i},\mathit{t}}\\ & +b_4\mathrm{loge}{\left(\mathrm{PHexp}\right)}_{\mathit{i},\mathit{t}}+b_5{\left[\mathrm{loge}\left(\mathrm{Dens}\right)\right]}_{\mathit{i},\mathit{t}}\\ & +b_6{\left[\mathrm{loge}\left(\mathrm{Temp}\right)\right]}_{\mathit{i},\mathit{t}}+b_7{\left(\mathrm{Virus}\right)}_{\mathit{i},\mathit{t}}\\ & +\mathit{\alpha }\left(\mathit{i}\right)+{\epsilon }_{\mathit{i},\mathit{t}},\end{array}$

for i = 1 … 21 and t = 2011 … 2017, where “P&I_mort_65+” is P&I mortality rate in subjects aged ≥65 years; bs are regression coefficients; α is the unobserved time-invariant regional effect (in FE model) or constant intercept (in RE model); i is a region; t is a year; ε is the error term; “IVC_65+” is IVC in subjects aged ≥65 years; “IVC_<65” is IVC in subjects aged <65 years “PHexp” is public health expenditure per capita; “Dens” is population density; “Temp” is average low winter temperature; “Virus” is a dummy variable indicating the predominant virus (sub)type.

Taking into account a high probability of heteroscedasticity and/or autocorrelation, all models considered also the Arellano’s heteroscedasticity-autocorrelation (HAC) robust standard errors (SEs). We performed the model diagnostics by applying the Breusch–Godfrey test for panel models to detect serial correlation for the errors and Pesaran cross-sectional dependence (CD) and Breusch–Pagan Lagrange multiplier tests for CD in the constructed panel models.^52,53

All the modeling was made in R stats packages.²⁷

Results

Exploratory spatiotemporal analysis of pneumonia- and influenza-related mortality in the Italian older adults aged ≥65 years

Over seven years (2011–2017), a total of 71,876 P&I-related deaths were reported in older adults aged ≥65 years. There was some variability in terms of P&I-related mortality rates observed between years. The highest rates were observed in years 2017 (10.05 per 10,000) and 2015 (8.78 per 10,000). In the remaining years the mortality rate was lower and 6.7 < 8 per 10,000 (Figure 1).

Figure 1. Choropleth maps of pneumonia- and influenza-related mortality rates (per 10,000) in older adults aged ≥65 years, by year.

*Data from four provinces of Sardinia Region (Olbia-Tempio, Ogliastra, Medio-Campidano, Carbonia-Iglesias) were not available for year 2017.

We then explored choropleth charts by mapping province-specific P&I-related death rates. A clear north–south gradient (especially in 2011, 2015, and 2017) was evident: compared with central and southern provinces, those located in the Northern Italy displayed higher P&I mortality rates (Figure 1). As shown by Moran’s Is, a significant (P < .001) clustered pattern of the observed mortality rates took place in all years with the I–value ranging from 0.28 to 0.36.

The LISA analysis (Supplementary Material, Figure S1) confirmed the north-south gradient: most hotspots and coldspots were located among northern and southern provinces, respectively. The few outliers detected (mainly cold–hot) were located in Sicily.

Association between pneumonia- and influenza-related mortality and influenza vaccination patterns

The total panel was composed of 147 observations and was balanced (i.e., no single space-time observations were missing). Table 1 reports principal descriptive statistics of the continuous variables of interest. During the study period, an average IVC in the elderly was 55.0%. Two significant drops in IVC were observed: the first occurred in 2012 (from 62.7% to 54.2%), the second in 2014 (from 55.4% to 48.6%). During the study period, only one region reached the recommended target of 75% (Umbria in the 2010/11 season). The proportion of aTIV use was highly non-homogeneous with a range of 0–76% (Table 1).

Table 1. Summary statistics of the continuous independent variables considered

Description	Mean	SD	Median	Min	Max	Ref
IVC in subjects aged ≥65 years, %	54.3	8.2	54.2	33.9	75.2	^Citation15
IVC in subjects aged <65 years, %	4.9	2.2	4.4	1.6	14.5	^Citation15,Citation20
aTIV regional allotments to total vaccine doses, %	28.3	18.0	28.6	0	76.3	Seqirus data
Public health expenditure per capita, €	1,885	148	1,854	1,652	2,283	^Citation20
Population density, inhabitants per km^2	178.0	110.5	162.7	38.8	429.1	^Citation20
Average low winter temperature, °C	6.9	2.7	6.8	1.0	12.9	^Citation21

aTIV: MF59®-adjuvanted trivalent influenza vaccine; IVC: influenza vaccination coverage.

For what concerns the seasonal dummy variables, seasons 2011/12, 2013/14 and 2016/17 were dominated by A/H3N2, season 2010/11 by A/H1N1pdm09, seasons 2012/13 and 2015/16 by B type, while the remaining 2014/15 season was ascribed by a co-circulation of A/H1N1pdm09 and A/H3N2 (Supplementary Table S2).

Results of FE and RE panel regression models are reported in Table 2. In the FE model, both IVC rate in the elderly, proportion of aTIV use, and average winter temperature were negatively associated with the observed P&I mortality rate in the Italian elderly. In particular, both FE and RE models predicted that each 1% increase in IVC rate in the elderly would be associated (P < .001) with a 1.6–1.9% decrease in P&I mortality. Analogously, each 1% increase in aTIV use would be associated (P < .05) with a 0.4% decrease in P&I mortality. By contrast, the co-circulation of type A virus subtypes was a significant positive predictor. No statistically significant association was observed for other independent variables. The output of the RE model was similar to that of the FE model. However, the Hausman’s test suggested (P < .001) that the FE model should be retained. The model diagnostics justified the use of both HAC robust standard errors (Table 2).

Table 2. Panel regression analysis on the association between pneumonia- and influenza-related mortality and potential predictors

Variable	Fixed-effects model				Random-effects model
	b	SE (P)	HAC SE (P)	b		SE (P)	HAC SE (P)
IVC in subjects aged ≥65 years	−0.019	0.004 (<0.001)***	0.004 (<0.001)***	−0.016		0.004 (<0.001)***	0.003 (<0.001)***
IVC in subjects aged <65 years	0.006	0.012 (0.61)	0.012 (0.63)	−0.006		0.012 (0.64)	0.011 (0.61)
Proportion of aTIV	−0.004	0.001 (0.013)*	0.001 (0.003)**	−0.004		0.001 (0.014)*	0.001 (0.006)**
Public health expenditure per capita (€ 1.000)	−0.714	0.656 (0.28)	0.869 (0.41)	0.287		0.574 (0.62)	0.906 (0.75)
Population density	−2.215	1.592 (0.17)	1.940 (0.26)	0.145		0.103 (0.16)	0.109 (0.19)
Average low winter temperature	−0.299	0.097 (0.003**)	0.072 (<0.001***)	−0.341		0.084 (<0.001***)	0.071 (<0.001***)
Predominance of A/H1N1pdm	Ref				Ref
Predominance of A/H3N2	0.036	0.043 (0.40)	0.029 (0.21)	0.053		0.045 (0.24)	0.028 (0.062∙)
Predominance of B	−0.059	0.055 (0.29)	0.040 (0.14)	−0.033		0.057 (0.56)	0.039 (0.40)
Co-circulation A/H1N1pdm09 and A/H3N2	0.108	0.068 (0.11)	0.049 (0.030*)	0.153		0.071 (0.033*)	0.047 (0.001**)
R^2, %	41.6				38.4
Hausman test, P	<0.001***
Breusch-Pagan test, P	<0.001***				<0.001***
Pesaran test, P	<0.001***				<0.001***
Breusch-Godfrey/Wooldridge test, P	<0.001***				<0.001***

P< .10; *P< .05; **P< 0.01; ***P< .001; aTIV: MF59-adjuvanted trivalent influenza vaccine; HAC SE: heteroscedasticity and autocorrelation robust standard errors; IVC: influenza vaccination coverage.

Discussion

This study confirms that annual influenza vaccination in the elderly reduces overall P&I mortality rates and therefore that increased immunization rates are desirable. From the ecological and policy-making perspectives this study has also confirmed the usefulness of aTIV in preventing P&I-related mortality in the elderly Italian population. This is in line with a recently proposed concept of the appropriate use of IVs in Italy.⁵⁴

The protective effect of IV on (excess) mortality is still controversial. For instance, Rizzo et al.¹⁰ and Simonsen et al.⁵⁵ have not documented any meaningful temporal association between IVC rate and winter excess mortality in Italy and the United States, respectively. By contrast, available meta-analyses of observational studies^56,57 suggest a significant reduction in mortality among vaccinated individuals. Contrary to the previous Italian time-trend study by Rizzo et al.,¹⁰ we were able to demonstrate a protective effect of IVs on P&I mortality. The reasons for this discrepancy are likely to be multiple. First, two different and non-overlapping time periods with both different circulating viruses and available IVs were assessed. Second, two different proxy outcomes to quantify influenza-related mortality were used. Third, in the present study both time and space were incorporated in the analysis; this may provide additional benefits in countries like Italy with its “jeopardized” pattern of IVCs.⁵⁸

Our second main finding was that regions using a higher proportion of aTIV showed significantly lower P&I mortality in the elderly independent of IVC, virus circulation pattern, and other potential confounders. In the elderly population aTIV has systematically been shown to be both more immunogenic⁵⁹ and effective⁶⁰ than standard-dose non-adjuvanted IVs. In particular, meta-analysis by Nicolay et al.⁵⁹ has shown that the use of aTIV was associated with significantly higher seroconversion rates and geometric mean titers against both vaccine antigen strains and heterologous strains, independently from the (sub)type analyzed. Analogously, the systematic review by Domnich et al.⁶⁰ concluded that the available observational studies had usually displayed a greater effectiveness of aTIV against various influenza-related outcomes, as compared with non-adjuvanted counterparts. In this regard the ecological study design may confirm findings coming from primary experimental or observational research and may therefore be useful in the decision-making process.

The main strength of this study lies in the methodology adopted. While a limited number of time-space observations does not allow the model to be adjusted for “anything we want to” as well as the fact that some potentially useful data may be not available, the FE panel models provide unbiased estimates in situations when unobserved confounders are present (as in the case of this study).³⁸

Apart from the well-known general limitations of ecological studies (ecological fallacy in primis),¹⁶ specific limitations apply to this study that must be considered. First, owing to data availability, the final regression models could be performed at the regional level only. A more space-detailed (e.g., provincial or local health unit level) evaluation is warranted. Second, the data on market share of single IVs come from regional tender allotments and therefore may not exactly correspond to effective IV administration. However, we believe that this limitation had a limited impact on the results since it is unlikely that the wastage/non-utilization rate differs among single IVs. Third, although we have no reason to believe that the automatic coding of causes of death may differ between regions, we could not completely rule out the between-region differences in reporting quality. However, this limitation may have only a small impact on the study results since the Italian National Institute of Statistics performs quality checks on a regular basis. Moreover, the FE panel model adopted may absorb this eventual heterogeneity.

To conclude, our analysis supports the increase in annual influenza vaccination in Italy and suggests that a higher IV uptake in the Italian elderly population would be beneficial. The use of aTIV in older adults is advised to reduce the burden of seasonal influenza disease.

Institutional review board statement

Ethical review and approval were waived for this study, due to the fact that this is based on the publicly available aggregated information sources. No single subject data were available.

Author contributions

Conceptualization, E.F. and A.D.; methodology, E.F. and A.D.; formal analysis, A.D. and A.S.; investigation, E.F. and A.D.; data curation, E.F. A.D. and A.S.; writing—original draft preparation, E.F. and A.D.; writing—review and editing, A.O. and G.I; supervision, A.O. and G.I.; project administration, A.D. All authors have read and approved the submitted manuscript.

Disclosure statement

E.F. is PhD student at Siena University whose program was fully funded by Seqirus, a pharmaceutical company who manufacture and commercialize influenza vaccines. At the time of submission, E.F. became a Seqirus permanent employee. A.D. was a permanent employee of Seqirus at the time study conception and realization. A.S. was remunerated by Seqirus for his statistical analysis. A.O and G.I. declare no conflicts of interest regarding this publication.

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials. In the supplementary materials, Table S1 summarizes all the data sources (and related url) used in the analysis.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website at https://doi.org/10.1080/21645515.2021.2005381.

Additional information

Funding

This research was funded by Seqirus, a pharmaceutical company who manufacture and commercialize influenza vaccines. The APC was funded by Seqirus.