The determinants of COVID-19 case fatality rate (CFR) in the Italian regions and provinces: An analysis of environmental, demographic, and healthcare factors.

The Italian government has been one of the most responsive to COVID-2019 emergency, through the adoption of quick and increasingly stringent measures to contain the outbreak. Despite this, Italy has suffered a huge human and social cost, especially in Lombardy. The aim of this paper is dual: i) first, to investigate the reasons of the case fatality rate (CFR) differences across Italian 20 regions and 107 provinces, using a multivariate OLS regression approach; and ii) second, to build a "taxonomy" of provinces with similar mortality risk of COVID-19, by using the Ward's hierarchical agglomerative clustering method. I considered health system metrics, environmental pollution, climatic conditions, demographic variables, and three ad hoc indexes that represent the health system saturation. The results showed that overall health care efficiency, physician density, and average temperature helped to reduce the CFR. By the contrary, population aged 70 and above, car and firm density, air pollutants concentrations (NO2, O3, PM10, and PM2.5), relative average humidity, COVID-19 prevalence, and all three indexes of health system saturation were positively associated with the CFR. Population density, social vertical integration, and altitude were not statistically significant. In particular, the risk of dying increases with age, as 90 years old and above had a three-fold greater risk than the 80-to-89 years old and four-fold greater risk than 70-to-79 years old. Moreover, the cluster analysis showed that the highest mortality risk was concentrated in the north of the country, while the lowest risk was associated with southern provinces. Finally, since prevalence and health system saturation indexes played the most important role in explaining the CFR variability, a significant part of the latter may have been caused by the massive stress of the Italian health system.

Introduction

The novel coronavirus disease (COVID-19) is a severe acute respiratory syndrome detected for the first time in December 2019 in Wuhan, Hubei province, China. On March 11, 2020, the World Health Organization (WHO, 2020) declared that COVID-19 could be characterized as a pandemic. As of July 26, 2020, according to Worldometer (2020), the virus has spread across 213 countries and territories, affecting over 16,4 million people and causing more than 650 thousand deaths. Italy was one of the countries hit the worst by the pandemic, with almost 250 thousand confirmed cases and over 35 thousand deaths at the time of writing. Despite the widely recognized excellence of the Italian health system (World Health Organization, 2010; GBD, 2017; Bloomberg, 2019), the country has paid a very high price, with one of the highest case fatality rate (CFR) in the world (https://ourworldindata.org/grapher/coronavirus-cfr). However, at the peak of the epidemic, mortality was affected by a large spatial heterogeneity; in fact, the northern regions were characterized on average by a significantly higher CFR than southern regions (Fig. 1 ).

Fig. 1

The case fatality rate (CFR) at the peak of the epidemic in the Italian provinces.

The aim of this work is to contribute to the existing literature by investigating the main reasons and determinants of COVID-19 CFR at the peak of the virus outbreak in 20 Italian regions and 107 Italian provinces. To reach this goal, I considered several variables: health system metrics, air pollutants concentrations, climatic conditions, demographic factors, and three ad hoc indexes that represent the health system saturation.

Specifically, I chose the peak of mild and severe cases of COVID-19, which was approximately reached on April 3–4, 2020 (Fig. 2 ), because it can be considered as the moment of maximum health care system saturation. This is an important aspect, because at the early outbreak stages the number of infected usually growths exponentially and health system cannot systematize its response. So, the health system saturation may assume greater importance.1

Fig. 2

The number of people recovered from COVID-19 in the period February 23, 2020–July 23, 2020.

Relevant literature

In the last months, a plenty and increasing body of literature focused its attention on the environmental, meteorological, demographic, and social factors that may affect COVID-19 mortality (Bayer and Kuhn, 2020; Brandt et al., 2020; Comunian et al., 2020; Du et al., 2020; Ma et al., 2020; Pansini and Fornacca, 2020; Sannigrahi et al., 2020; X. Wu et al., 2020; Verity et al., 2020; Zhu et al., 2020). What follows is a summary of the main findings on the relationship between COVID-19 CFR and climatic conditions, demographic variables, and air pollutants concentrations, respectively.

Ma et al. (2020) used a generalized additive model (GAM) to study the relationship between the meteorological factors and the daily deaths of COVID-19 in Wuhan from January 20, 2020 to February 29, 2020. They found a positive association between the daily deaths of COVID-19 and the diurnal temperature range (DTR), and a negative relationship between the former and the relative humidity and temperature. Similarly, Y. Wu et al. (2020) used a log-linear GAM to analyze the effect of temperature and humidity on daily new deaths of COVID-19 in 166 countries, as of March 27, 2020. Temperature and relative humidity were found to be both significantly and negatively associated with daily new deaths. Rahman et al. (2020) performed partial correlation analysis and linear mixed effect modeling to analyze the effect of temperature on COVID-19 mortality risk in 149 countries. They showed that higher temperatures were negatively associated with mortality in high-income countries, while extreme temperature may increase mortality risk in low-and middle-income countries.

Sannigrahi et al. (2020) used spatial regression models to analyze the spatial association between the key demographic variables and COVID-19 deaths across 31 European countries. They found that the incidence of the population aged 80 and above on overall casualties caused by COVID-19 was highly significant. Verity et al. (2020) used a model-based approach to estimate the case fatality ratio associated with age groups in mainland China. They found that the population aged 80 and above had the highest case fatality ratio (13.4%). Similarly, Du et al. (2020), by implementing a univariate and multivariate logistic regression analysis, showed that advanced age is a significant risk factor for COVID-19 mortality. Moreover, Ioannidis et al. (2020), by using official data from 14 countries and 13 US states as of June 17, 2020, estimated that age risk gradients were highly significant. In particular, people aged less than 65 had 16–100-fold lower risk of COVID-19 deaths than older people.

In early March 2020, Bayer and Kuhn (2020) advanced the hypothesis that Italian higher vertical social integration may arise the COVID-19 fatality rate. They used a sample of 24 countries with at least 200 COVID-19 confirmed cases and found a positive correlation between the share of the population aged 30–49 living with parents and COVID-19 death rate. However, this relation was strongly criticized by Belloc et al. (2020), which showed how the sign of the correlation turned negative when considering the variation within Italian regions.

Regarding the effect of air pollutants concentrations on COVID-19 related deaths, Conticini et al. (2020), by investigating the relevant literature, concluded that prolonged exposure to air pollutants may lead to chronic respiratory conditions, even in healthy and young people. X. Wu et al. (2020) analyzed COVID-19 death counts for 3087 counties in the USA, covering 98% of the population, by using a negative binomial mixed model. They found that a positive and significant association between PM2.5 and COVID-19 mortality rates. In particular, a 1 unit increase in PM2.5 is related to an 8% increase in the COVID-19 fatality rate. Pansini and Fornacca (2020), using Kendall's tau and Pearson correlation coefficient, found a positive and significant association between COVID-19 mortality and several air pollutants (CO, NO2, PM10, PM2.5) in China and the USA.

Ogen (2020) used spatial analysis to examine the relationship between long-term exposure to NO2 and COVID-19 mortality in 66 administrative regions in France, Germany, Italy, and Spain. He showed that 78% of the total COVID-19 deaths were located in north Italy and central Spain, i.e. the regions with the highest level of NO2. Bianconi et al. (2020), using multiple linear regression models, showed that mean annual exposure to PM10 and PM2.5 was significantly and positively associated with COVID-19 death rate in Italian 20 regions. Yao et al. (2020), by using spatial analysis and multivariate linear regression for China as of April 12, 2020, found that every 10 μg/m3 increase in PM10 and PM2.5 concentrations is associated with a 0.24% and 0.26% increase in the COVID-19 mortality rate, respectively. Brandt et al. (2020) also stressed that air pollution levels are strongly associated with densely populated urban areas in the USA. Hamidi et al. (2020) investigated both direct and indirect impacts of population density on death rates in 913 USA metropolitan countries by using structural equation modeling. They found that larger metropolitan areas and counties with higher population density were significantly associated with higher mortality rates.

Finally, several studies have shown that the presence of at least one comorbidity, such as hypertension, diabetes, cardiovascular disease, and chronic lung disease, negatively affects outcomes of patients hospitalized with COVID-19 (Guan et al., 2020; Istat-ISS, 2020; Jordan et al., 2020; Wang et al., 2020; Wu and McGoogan, 2020).

Material

At regional level, I used the following 17 explanatory variables (Table 1 ): an overall index (IPS) of the Italian health system performance in 2017–2018, the public health expenditure per capita in the period 2015–2017, the total specialist doctors and general practitioners per 1000 inhabitants in the period 2016–2018, the total ordinary hospital beds per 1000 inhabitants in the period 2016–2018, an index of car and firm density in 2015–2017, the electric power consumption (kWh per capita) in the period 2016–2018,2 the proportion of population aged 70 and over in 2019, the proportion of population aged 80 and over in 2019, the proportion of population aged 90 and over in 2019, the vertical social integration proxied by the share of adults aged 18–34 living with their parents in 2019, the average relative humidity levels registered during March 2020, the average historical diurnal temperature range (DTR) in March, the average historical temperature in March, the average prevalence of COVID-19 on April 3 and 4, 2020, the ratio between the COVID-19 prevalence and the ordinary hospital beds, an two ad hoc indexes that indicate the saturation of ordinary hospital beds (OB) and critical care beds (CCB) at the peak of the outbreak, respectively. The dependent variable is represented by the average case fatality rate of COVID-19 on April 3 and 4, 2020.3

Table 1

Definitions of all variables used for OLS regional analysis.

Variables	Definitions	Sources
Dependent variables
CFR	The average case fatality rate for COVID-19 in each region, obtained by dividing the average confirmed deaths by the average confirmed cases on April 3 and 4, 2020.	Italian Ministry of Healtha
Independent variables
IPS	A synthetic index of the Italian health system performance in the period 2017–2018, which includes eight different parameters.b	Demoskopika, 2018, Demoskopika, 2019
Health expenditure	The average public health expenditure per capita for each region, in the period 2015–2017.	I.Stat (database)c
Physicians	The average total specialist doctors and general practitioners (per 1000 inhabitants) for each region, in the period 2016–2018.	I.Stat (database)
Hospital beds	The average ordinary hospital beds (per 1000 inhabitants) for each region, in 2016–2018.	Italian Ministry of Health
Cars & Firms	A synthetic index of car and firm (> 250 employees) density for each region, in 2015–2017.d	I.Stat (database)
kWh per capita	The average electric power consumption in kilowatt-hours (kWh) per capita for each region, in the period 2016–2018.	Terna (2019)
Ages 70+	The proportion of population aged 70 and over for each region, in 2019.	I.Stat (database)
Ages 80+	The proportion of population aged 80 and over for each region, in 2019.	I.Stat (database)
Ages 90+	The proportion of population aged 90 and over for each region, in 2019.	I.Stat (database)
Vertical integration	The share of unmarried young adults aged 18–34 living with at least one parent for each region, in 2019.	I.Stat (database)
Humidity	The average relative humidity levels registered during March 2020, for each region.e	www.il meteo.itf
DTR	The historical diurnal temperature range in March, for each region.	Mipaaf (2019a)
Temperature	The historical average temperature in March, for each region.	Mipaaf (2019a)
Prevalence	The average ratio between the people who have been tested positive for COVID-19 and the overall population of each region on April 3 and 4, 2020.	I.Stat (database), Italian Ministry of Health
Preval./Beds	The ratio between the average COVID-19 prevalence on April 3 and 4, 2020, and the average number of ordinary hospital beds in 2016–2018, for each region.	Italian Ministry of Health
CCB saturation	The ratio between the average people who have been recovered from COVID-19 in intensive care on 3 and 4 April 2020, and the average number of critical care beds (CCB) in the period 2016–2018, for each region.	Italian Ministry of Health
OB saturation	The ratio between the average people who have been recovered from COVID-19 with mild symptoms on 3 and 4 April 2020, and the average number of ordinary hospital beds in the period 2016–2018 for each region.	Italian Ministry of Health

Data are available at URL: www.salute.gov.it.

The parameters used are the following: patient satisfaction, active patient mobility, passive patient mobility, legal fees for disputes, operating result, life expectancy, equality in health treatment, and economic hardship. In particular, each parameter is standardized, with mean = 100 and standard deviation = 10, and the final synthetic index is obtained by calculating the simple average of them.

Data are available at URL: http://dati.istat.it/.

The number of cars refers to those recorded in the Pubblico registro automobilistico (Public vehicle register). The number of largest firms (> 250 employees) refers to those that operate in the following sectors: (1) mining and minerals from quarries and mines; (2) manufacturing activities; (3) supply of electricity, gas, vapors, and air conditioning; and (4) supply sewerage, waste management and remediation activities. The index is compiled according to the following analytical method: i) first, I standardized the data according to surface area (cars and firms for 100 sq. km.); ii) than, the respective outputs are switched to fixed-base indexes (with mean = 100); iii) finally, I computed the simple arithmetic mean of the latter.

The average values have been calculated by dividing the data coming from 62 different official weather stations, managed by the Italian Air force and located in the main Italian provinces.

This is one of the most trusted Italian weather forecast website. https://www.ilmeteo.it/business/assets//images/aboutUs/pdf/Google%2009_04_2020.pdf.

At the province level, I used the following 18 explanatory variables (Table 2 ): the average general practitioners per 1000 inhabitants in the period 2015–2017, the historical average temperature in March, the historical average diurnal temperature range (DTR) in March, an ordinal index of the population structure (rural-intermediate-urban) in 2013, the population density (people per sq. km) in 2019, the proportion of population aged 70 and over in 2019, the proportion of population aged 70–79 in 2019, the proportion of population aged 80–89 in 2019, the proportion of population aged 90 and over in 2019, the average concentrations (in μg/m3) of particulate matter less than 10 μm in diameter (PM10) in 2017–2018, the number of days in which PM10 concentrations exceeded the legal limit of 50 μg/m3,4 the average concentrations (in μg/m3) of particulate matter less than 2.5 μm in diameter (PM2.5) in the period 2017–2018, the average concentrations of nitrogen dioxide (NO2) expressed in μg/m3 in 2017–2018, the number of days in which ozone (O3) concentrations exceeded the limit of 120 μg/m3,5 the average altitude of the capital city of each province, the prevalence of COVID-19 on March 31, 2020, and the ratio between the COVID-19 prevalence on March 31 and the average number of ordinary hospital beds in the period 2016–2018. The dependent variable is represented by the case fatality rate of COVID-19 on March 31, 2020.6 Further details are provided in Table 1, Table 2.7

Table 2

Definitions of all variables used for OLS province analysis.

Variables	Definitions	Sources
Dependent variables
Death rate	The average case fatality rate for COVID-19 in each province, obtained by dividing the confirmed deaths by the number of confirmed cases, on 31 March 2020.	Istat-ISS (2020)
Independent variables
General practitioners	The average general practitioners for each province, in 2019.	Il Sole 24 Ore (2019)
Temperature	The historical average temperature in March, for each province.	Mipaaf (2019b)
DTR	The historical average diurnal temperature range in March, for each province.	Mipaaf (2019b)
Urbanization	An ordinal index that ranks population of each province by urban-rural structure: predominantly rural (1), intermediate (2), and predominantly urban (3).	European Commission (2013)
Density	The number of human inhabitants per square kilometer (sq. km.) of land area for each province, in 2019.	I.Stat (database)
Ages 70+	The proportion of population aged 70 and over for each province, in 2019.	I.Stat (database)
Ages 70–79	The proportion of population aged 70–79 for each province, in 2019.	I.Stat (database)
Ages 80–89	The proportion of population aged 80–89 for each province, in 2019.	I.Stat (database)
Ages 90+	The proportion of population aged 90 and over for each province, in 2019.	I.Stat (database)
PM10	The average concentrations of particulate matter less than 10 μm in diameter, expressed in μg/m3 for each province, in the period 2017–2018.	Istat (2019)
PM10 (>50)	The average number of days in which PM10 exceeded the limit of 50 μg/m3 for each province, in the period 2017–2018.	Istat (2019)
PM2.5	The average concentrations of particulate matter less than 2.5 μm in diameter, expressed in μg/m3 for each province, in the period 2017–2018.	Istat (2019)
NO2	The average concentrations of nitrogen dioxide, expressed in μg/m3 for each province, in the period 2017–2018.	Istat (2019)
O3	The average number of days in which ozone exceeded the limit of 120 μg/m3 for each province, in the period 2017–2018.	Istat (2019)
Altitude	The average altitude of the capital city of each province.	Istat (2019)
Prevalence	The ratio between the people who have been tested positive for COVID-19 on March 31, 2020, and the total population of each province in 2019	I.Stat (database), Italian Ministry of Health
OB saturation	The ratio between the COVID-19 prevalence on March 31, 2020, and the average number of ordinary hospital beds in 2016–2018, for each province.	Italian Ministry of Health

Methods

First, I used a multivariate cross-sectional OLS (ordinary least squares) approach to identify the main determinants at regional at province level, and then I applied the Ward's hierarchical agglomerative clustering method (Ward Jr, 1963) to build a “taxonomy” of provinces with similar mortality risk of COVID-19. The rationale behind the choice of cross-sectional regression methodology instead of a panel approach is as follows: i) comprehensive daily data on March 2020 are not currently available for air pollutants concentrations and climatic variables; and ii) demographic and health system variables change very little or not at all in the short term.

The OLS equation estimated at the regional and provincial level was as follows:where i represented the regions/provinces, β 0 was the intercept, X1 to Xn were the independent variables for each region/province, and εi was the random error.

Then, I conducted a cluster analysis on selected significant variables obtained from multivariate OLS on province data. The procedure encompassed the following 4 sequential steps. First, since variables are not on the same scale, I calculated the standardized data matrix by using the following formula:where x was the value of the variable in the original dataset, μ was the arithmetic mean of the original variable, and σ was the standard deviation of the latter. Second, I computed the Euclidean distance. Given two points X and Y in d dimensional space, the Euclidean distance between X and Y was equal to:

Then, I applied Ward's hierarchical clustering method, which allowed to obtain clusters of provinces with features as similar as possible by minimizing the total within-cluster variance. Specifically, at each step, the pair of clusters that are characterized by minimum between-cluster distance were merged. Therefore, Ward's method merging cost formula between two clusters, p and q, was given by:

From this, it was obtained:where m was the center of cluster j, and n was the overall number of points included in cluster j.

Results and discussion

OLS at the regional level

In Table 3a, Table 3b ), I presented the OLS estimations at the regional level. All the 12 estimated OLS models were statistically significant; in fact, the Fisher-Snedecor distribution assumed values far higher than the tabulated critical values at the 1% level of significance. In particular, the r-square showed that models were able to explain from 0.44% to 0.88% of CFR variability. Furthermore, since Breusch and Pagan (1979) and Shapiro and Wilk (1965) tests allowed to accept the null hypothesis of homoscedasticity and normality of residuals, models seemed well specified. However, due to the small sample, I preferred to adopt a conservative approach, by applying the HC2 correction proposed by MacKinnon and White (1985).8 It is important to stress that OLS cross-sectional analysis is very sensitive to the presence of outliers, which can cause misspecification issues (Mur and Lauridsen, 2007), especially in such small samples (Wooldridge, 2015, p. 334). Therefore, I investigated the presence of highly influential points. The analysis revealed that none of the leverage points (h) are beyond the cutoff value (h > 2n/k) proposed by Belsley et al. (1980). Finally, I also identified the presence of multicollinearity. The variance inflation factors (VIF) ranged from 1.46 to 8.289 and were less than 10, i.e. the rule of thumb suggested by the relevant literature (Belsley, 1982; Hair Jr. et al., 1995). Therefore, I concluded that the independent variables just suffered from weak linear dependency and there were no serious multicollinearity issues.

Table 3a

(Models 1–6). OLS regression at the regional level between CFR and environmental, demographic, and healthcare factors.

Variables	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6
Constant	16.4788[13.9432]	21.4019[12.8697]	27.8448*[12.7974]	−12.1666[28.5844]	−11.7866[12.8946]	−14.8309[13.3257]
IPS	−0.4788*[0.2303]	−0.5451**[0.2221]	−0.5731**[0.2153]	−0.5758**[0.2034]	−0.4511**[0.1402]	−0.4808***[0.1468]
Health exp.	0.0129**[0.0059]	0.0157**[0.0056]	0.0144**[0.0058]	0.0237**[0.0096]	0.0146*[0.0074]	0.0168**[0.0067]
Physicians	−3.2988**[1.2894]	−3.1879**[1.1264]	−3.1683**[1.136]	−2.3358[1.8958]	−3.0878**[1.3113]	−2.5227*[1.1406]
H. Beds	3.8867[2.8595]	2.7984[2.8415]	3.2841[3.1279]	−1.0369[5.1024]	−1.4253[1.7597]	−2.1594[2.191]
Car & Firm	6.2449***[1.7936]	7.0056***[1.7842]	6.9565***[1.93]	8.5611***[1.9709]	4.5281**[1.4903]	5.4669***[1.4142]
Kilowatt	−0.0007[0.0006]	−0.0005[0.0006]	−0.0005[0.0007]	−0.0015*[0.0007]	−0.0009[0.0007]	−0.0012*[0.0006]
Aged 70+	0.7821**[0.3183]			1.0916**[0.4199]	0.7973***[0.2343]	0.8982***[0.234]
Aged 80+		1.5**[0.5131]
Aged 90+			6.3807**[2.5141]
Humidity				0.449**[0.1886]	0.1413[0.1535]	0.2615**[0.0808]
DTR				0.3657[1.5456]	2.679***[0.6612]	2.3164***[0.6588]
Temperature				−1.1144*[0.5752]	0.5302[0.4919]
Prevalence					26.5164***[4.3121]	22.0411***[5.3084]
Breusch-P. (p)	0.2454	0.5101	0.3985	0.7505	0.4038	0.6391
Shapiro-W. (p	0.9892	0.9969	0.9713	0.7997	0.0115	0.1289
Influential (h)	0.13–0.78	0.15–0.74	0.17–0.69	0.23–0.82	0.25–0.87	0.2–0.87
VIF	1.56–4.44	1.46–4.45	1.49–4.88	2.49–5.43	2.65–8.28	2.04–5.79
F-statistic	5.85***	5.28***	5.5***	5.18***	22.13***	22.55***
Observations	20	20	20	20	20	20
Adjusted R2	0.4425	0.4753	0.4531	0.4899	0.8102	0.8127

Notes: h, leverage; p, p-value. Standard errors (in brackets) are based on HC2 method developed by MacKinnon and White (1985). Significance level: p-value <0.01***; p-value <0.05**: p-value <0.1*.

Table 3b

(Models 7–12). OLS regression at the regional level between CFR and environmental, demographic, and healthcare factors.

Variables	Model 7	Model 8	Model 9	Model 10	Model 11	Model 12
Constant	−9.9732[12.8992]	−3.5934[12.9892]	−12.2744[13.7815]	168.707**[65.3856]	−4.8115[20.9305]	8.2251[12.6372]
IPS	−0.5065***[0.1497]	−0.5475***[0.1504]	−0.4995***[0.1506]	−0.671***[0.1663]	−0.6399***[0.1731]	−0.4174**[0.1348]
Health exp.	0.0151**[0.0056]	0.0155**[0.0056]	0.0163**[0.0069]	0.012**[0.0051]	0.0156[0.0099]	0.0062[0.0074]
Physicians	−2.4044*[1.2974]	−1.8737[1.2409]	−2.4275*[1.1621]	−1.8232[1.0917]	−1.838[1.3513]	−0.8314[0.9748]
H. Beds			−0.7423[2.3357]	−1.2938[1.4252]	1.3549[3.2707]	−1.0043[2.4216]
Car & Firms	5.2252***[1.391]	5.7172***[1.4174]	5.395***[1.4313]	7.2545***[1.2172]	6.4978***[1.4789]	0.6217[1.6553]
Kilowatt	−0.0012[0.0007]	−0.0009[0.0007]	−0.0012*[0.0006]	−0.001*[0.0004]	−0.0007[0.0006]	0.0000[0.0004]
Aged 70+	0.7931**[0.263]		0.7943**[0.2501]		0.8673**[0.3488]	0.2744[0.3046]
Aged 80+		1.3558**[0.4649]		2.1616***[0.6187]
Vertical Integ.	−0.008[0.1366]	0.0254[0.1343]
Humidity	0.2455**[0.0956]	0.2108**[0.0927]	0.2536**[0.084]	−4.6525**[1.8055]	0.2824**[0.1213]	−0.0854[0.0711]
Humidity^2				0.036**[0.0134]
DTR	1.9813**[0.7199]	1.6247**[0.6522]	2.0091**[0.6946]	2.5458***[0.6437]	1.18[1.0982]	2.3931**[0.8386]
Preval./beds	7.2081***[2.0728]	7.121***[2.0874]	7.3835***[1.8403]	6.9567***[1.1719]
ICB saturation					4.6762**[1.5899]
OB saturation						45.4811***[7.9362]
Breusch-P. (p)	0.5386	0.6034	0.6502	0.5497	0.8678	0.645
Shapiro-W (p)	0.087	0.8456	0.0951	0.9033	0.1169	0.6978
Influential (h)	0.2–0.85	0.21–0.85	0.21–0.87	0.23–0.92	0.2–0.81	0.2–0.81
VIF	2.32–5.89	2.08–5.91	1.95–5.43	–	2.31–5.48	2.03–7.25
F-statistic	18.96***	25.07***	21.32***	20.24***	9.57***	54.14***
Observations	20	20	20	20	20	20
Adjusted R2	0.7965	0.799	0.7982	0.8782	0.7024	0.8576

Notes: h, leverage; p, p-value. Standard errors (in brackets) are based on HC2 method developed by MacKinnon and White (1985). Significance level: p-value <0.01***; p-value <0.05**: p-value <0.1*.

The results showed that the IPS index and physicians were significantly and negatively associated with COVID-19 CFR.10 By the contrary, health expenditure, car and firm density, people aged 70 and above, humidity, DTR, prevalence, the ratio of prevalence/ordinary beds, and saturation indexes of ordinary and critical care beds were significantly and positively correlated with the CFR. Hospital beds, kWh per capita, and social vertical integration showed no statistical significance. Even if the sign of health expenditure is unexpected, the IPS index seemed to confirm the importance of health system effectiveness in all its dimensions more than only health expenditure. The regression coefficients also gave interesting information. The change in 1 unit of physician per 1000 inhabitants and IPS index was correlated on average with a change of −2.87% and − 0.53% in the CFR, respectively.11 Consistent with Du et al. (2020), Ioannidis et al. (2020), Sannigrahi et al. (2020), the risk of dying increased progressively with the increase in the elderly population. In fact, the regression coefficients of the population aged 90 and above were 4 times larger than that for the population aged 80 and above, and 8 times larger than that for the population aged 70 and above.

Moreover, each 1 unit increase in humidity and DTR was positively associated on average with a change of 0.28%12 and 2.22% in the CFR. The results for humidity are consistent with a recent study by Bianconi et al. (2020) on the Italian case, while DTR outcomes are in line with Ma et al. (2020). Prevalence and health system saturation indexes were highly significant in all the models. In particular, saturation of ordinary beds had the largest coefficient, followed by COVID-19 prevalence, the ratio prevalence/beds, and the saturation of intensive care beds. They played the most important role in explaining CFR variability; in fact, the average r-square of models (5–12) with prevalence and saturation indexes was 0.81, i.e. almost twice the models (1–4) without them (r = 0.46). In particular, the change in 0.1 unit of CCB and OB saturation was correlated with a change of 0.47% and 4.54% in the CFR. Therefore, a significant part of the CFR may be caused by the massive stress of the Italian health system.

Finally, the correlation matrix in Fig. S1 indicated that the average temperature and humidity were significantly and negatively correlated to the COVID-19 prevalence, with a Pearson's r of −0.8 and −0.46, respectively. By the contrary, the kWh per capita and car and firm density showed a good positive correlation with the latter, with a Pearson's r of 0.66 and 0.52, respectively.

OLS at the provincial level

In Table 4a, Table 4b , I presented the OLS estimations at the provincial level. All 11 estimated OLS models were statistically significant; in fact, the Fisher-Snedecor distribution assumed values far higher than the tabulated critical values at 1% level of significance. In particular, r-square showed that models were able to explain from 0.27% to 0.45% of CFR variability.

Table 4a

(Models 1–6). OLS regression at the provincial level between CFR and environmental, demographic, and healthcare factors.

Variables	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6
Constant	21.3371*[10.8017]	22.645**[11.0182]	21.4665**[10.4233]	28.7773***[10.0096]	27.8601**[11.075]	28.3833**[11.7957]
G. P.	−10.0768**[4.2187]	−9.4475**[4.1747]	−10.7185**[4.298]	−10.8571**[4.2279]	−9.9708**[4.0876]	−12.8379***[4.2603]
Temperature	−0.5506**[0.242]	−0.5966**[0.236]	−0.5292**[0.2431]	−0.5296**[0.2449]	−0.4515*[0.2455]	−0.4038[0.255]
DTR	−2.5289**[1.0323]	−2.6147**[1.0539]	−2.3998**[1.0138]	−2.6982**[1.0269]	−2.516**[1.0833]	−2.5241**[1.049]
PM10 (μg/m3)	0.2569***[0.0939]	0.2529***[0.0952]	0.256***[0.0938]	0.252***[0.0931]
Urbanization	0.8979[0.9731]	0.7571[0.9805]	1.0869[0.97]	1.0815[0.9602]
Pop. Density	−0.0004[0.0012]	−0.0005[0.0012]	−0.0004[0.0012]	−0.0003[0.0012]
Aged 70+	0.7596**[0.3106]				0.6823**[0.3157]	0.6955**[0.3149]
Aged 70–79		1.2866**[0.6163]
Aged 80–89			1.9984***[0.7511]
Aged 90+				5.5753**[2.2503]
PM10 (> 50)					0.0739***[0.0257]
NO2 (μg/m3)						0.1268**[0.055]
Breusch-P p.	0.0334	0.05	0.043	0.0385	0.0049	0.0083
F-statistic	6.37***	5.9***	6.33***	6.95***	9.81***	9.0808***
VIF	1.11–1.56	1.09–1.56	1.12–1.55	1.18–1.57	1.06–1.43	1.07–1.52
Observations	107	107	107	107	107	107
Adjusted R2	0.2944	0.2775	0.306	0.2903	0.301	0.266

Notes: p, p-value. Standard errors (in brackets) are based on HC2 method developed by MacKinnon and White (1985). Significance level: p-value <0.01***; p-value <0.05**: p-value <0.1*.

Table 4b

(Models 7–12). OLS regression at the provincial level between CFR and environmental, demographic, and healthcare factors.

Variables	Model 7	Model 8	Model 9	Model 10	Model 11	Model 12
Constant	31.4907***[9.6085]	41.079***[11.7581]	14.7429[11.7796]	18.1456[11.3803]	15.5861[11.579]	18.8624[11.2617]
G. P.	−12.4582***[4.642]	−13.0142***[4.2388]	−7.0346**[3.4178]	−7.6595**[3.5788]	−6.9201*[3.5378]	−7.578**[3.7071]
Temperature	−0.22[0.262]	−0.4403*[0.2471]	−0.1162[0.2952]	−0.129[0.2905]	−0.2588[0.2778]	−0.2667[0.2745]
DTR	−2.7569***[1.0324]	−4.0316***[1.2007]	−1.7712*[1.0379]	−1.7872*[1.0266]	−1.8775*[1.055]	−1.8874*[1.0434]
PM10 (μg/m3)			0.0476[0.0476]	0.0454[0.09]	0.0961[0.0989]	0.0939[0.0989]
Aged 70+	0.5006[0.3969]	0.6215*[0.3324]	0.64**[0.3028]		0.6583**[0.3031]
Aged 80+				1.178**[0.5477]		1.2287**[0.5514]
O3 (> 120)	0.1084***[0.0293]
PM2.5 (μg/m3)		0.3127***[0.1112]
Altitude			−0.0002[0.0021]	−0.0005[0.0021]	−0.0002[0.0022]	−0.0005[0.0022]
Prevalence			16.7325***[3.7856]	16.5201***[3.842]
OB saturation					4.4258***[1.0094]	4.3711***[1.0332]
Breusch-P p.	0.0003	0.0045	0.0000	0.0000	0.0003	0.0004
F-statistic	14.22***	10.48***	11.67***	11.77***	14.24***	14.14***
VIF	1.07–1.54	1.05–1.33	1.13–2.52	1.14–2.5	1.13–2.44	1.14–2.42
Observations	88	92	107	107	107	107
Adjusted R2	0.4241	0.3531	0.4488	0.4492	0.4142	0.416

Notes: p, p-value. Standard errors (in brackets) are based on HC2 method developed by MacKinnon and White (1985). Significance level: p-value <0.01***; p-value <0.05**: p-value <0.1*.

Since the Breusch-Pagan (1979) test revealed heteroscedasticity issues, I applied the HC2 correction. As stated by Ghasemi and Zahediasl (2012), the violation of normality assumption should not be a major problem in a sample with enough observations (n > 40), such as in this case. Furthermore, the variance inflation factors (VIF) ranged from 1.05 to 2.56 and are much less than 10; thus, there were no multicollinearity issues.

The OLS models showed that general practitioners, average temperature, and DTR were significantly and negatively associated with the CFR.13 By the contrary, all the considered air pollutants (PM10, PM2.5, NO2, and O3), people aged 70 and above, prevalence, and saturation index for ordinary beds were significantly and positively correlated with the CFR. Urbanization, population density, and altitude showed no statistical significance. The change in 1 degree Celsius of temperature and DTR was correlated with an average change of −0.53% and −2.44% in the CFR, respectively. If the results for the temperature are consistent with other studies (Bianconi et al., 2020; Rahman et al., 2020; Y. Wu et al., 2020), the sign of DTR is unexpected and in contrast with previous findings. Therefore, the relative impact is ambiguous and unclear.

Among air pollutants, PM2.5 and PM10 had the largest impact on COVID-19 CFR, followed by NO2 and O3. The change in 10 μg/m3 of PM2.5 and PM10 was associated with a change of 3.13% and 2.54% in the CFR, respectively. The change in 10 μg/m3 of NO2 was “only” correlated with a change of 1.27% in the CFR. While, each day in which O3 exceeded the limit of 120 μg/m3, the CFR have increased of 0.11%. The inclusion of O3 is very significant, since it allowed r-square to increase to 0.42. The positive and significant relationship between air pollutants and CFR is consistent with the relevant literature (Bianconi et al., 2020; Ogen, 2020; Pansini and Fornacca, 2020; X. Wu et al., 2020; Yao et al., 2020). As found at the regional level, the CFR increased progressively with population age. Specifically, 90 years old and above had a three-fold greater risk than the 80–to–89 years old and four-fold greater risk than 70–to–79 years old. This result is consistent with Livingston and Bucher, 2020 and Onder et al. (2020), and assumes great importance, since the susceptibility to the COVID-19 is constant across all age groups (Istat-ISS, 2020). The prevalence and saturation indexes were highly and significantly correlated with the CFR in all the considered models (9–12). Every 0.1 unit increase in COVID-19 prevalence and saturation of hospital ordinary beds was associated with an average change of 1.66% and 0.44% in the CFR, respectively. Moreover, the inclusion of prevalence and saturation in the models allowed an increase in the average r-square from 0.31 to 0.43. This is consistent with the previous findings.

The correlation matrix (Fig. S2) gave other interesting information, especially regarding air pollutants. Air pollutants concentrations were highly and positively correlated with each other, with a Pearson's r that ranged from 0.48 to 0.84. Moreover, consistent with Brandt et al. (2020), population density was moderately and positively correlated with air pollutants concentrations. Finally, as found in other studies (Bianconi et al., 2020; Pansini and Fornacca, 2020; Setti et al., 2020; Zhu et al., 2020), COVID-19 prevalence and air pollutant showed a good positive correlation, with a Pearson's r ranged from 0.36 to 0.65.

Cluster analysis

In Fig. 3 , I presented the dendrogram obtained using Ward's method. The optimal number of clusters was identified by using two different methods: i) the EM (Expectation-Maximization) algorithm for Gaussian finite mixture model, proposed by Fraley et al. (2012) in the package ‘Mclust’ (R environment); and iii) the package ‘NbClust’ (R environment) that proposes 30 different indexes.14 In the first case, according to BIC (Bayesian Information Criterion) score, the best model was VVI (varying volume, varying shape, and equal orientation) with 3 clusters. The majority of indices (7 out of 30) included in ‘NbClust’, also proposed 3 clusters. Therefore, by cutting the dendrogram at an approximately height of 13, I obtained three clusters with an increasing risk of mortality (Table 5 ). Each cluster was identified with a grey dotted rectangle (Fig. 3) and represented graphically through a map (Fig. 4 ). The cluster with the highest mortality risk (risk = 3) was composed by 25 provinces of northern Italy, 1 province of central Italy (Pesaro and Urbino), and 2 provinces of southern Italy (Frosinone and Avellino). The cluster with the medium mortality risk (risk = 2) was composed by 22 provinces of northern Italy, 18 provinces of central Italy, and 6 provinces of southern Italy (Chieti, L'Aquila, Oristano, Pescara, South Sardinia, and Teramo). Finally, the cluster with the lowest mortality risk (risk = 1) was composed by 2 provinces of central Italy (Latina and Rome), and 31 provinces of southern Italy.15 Therefore, the highest mortality risk was concentrated in the north of the country, while the lowest risk was associated with southern provinces. Specifically, cluster 3 had an 8.62% higher CFR than cluster 1, and 5.3% higher CFR than cluster 2. Cluster 3 had a temperature of 2.58 degrees lower and a number of general practitioners per 1000 inhabitants of 0.24 lower than that for cluster 1. Moreover, cluster 3 had a PM10 of 11.23 μg/m3 higher and a proportion of population aged 70 and above of 1.49% higher than that for cluster 1. Finally, the hospital bed saturation in cluster 3 was more than double than that for cluster 2, and 8 times greater than that for cluster 1. The results are consistent with Grasselli et al. (2020), according to which Lombardy's intensive care units already had an 85% to 90% occupancy ahead of the outbreak.16

Fig. 3

Dendrogram of provinces obtained using Ward's method.

Table 5

The clusters obtained by cutting dendrogram at an approximately height of 13.

Variables	Cluster (CL1)(Low risk)	Cluster (CL2)(Medium risk)	Cluster (CL3)(High risk)	CL3 - CL1
CFR	4.5758	7.8961	13.1957	8.6199
G. practitioners	1.0322	0.939	0.7947	−0.2375
Temperature	10.2758	7.1516	7.6946	−2.5812
Aged 70+	15.9355	19.3297	17.4285	1.493
PM10 (μg/m3)	21.9211	21.9833	33.1464	11.2253
OB saturation	0.1393	0.5019	1.161	1.0217
Numerosity	33	46	28	–

Fig. 4

Map of the 107 Italian provinces divided into three increasing clusters of risk.

Conclusions

To the best of my knowledge, this is one of the first studies to investigate the relationship between a wide set of heterogeneous factors and COVID-19 mortality. The OLS analysis showed that environmental, demographic, and healthcare factors played an important role in explaining the CFR variability. In particular, population aging, air pollutants concentrations (NO2, O3, PM10, and PM2.5), relative humidity, COVID-19 prevalence, and critical care and ordinary beds saturation were positively correlated with the CFR. By the contrary, overall health care efficiency (IPS), physician density, and average temperature were negatively associated with CFR. Specifically, the inclusion of the COVID-19 prevalence and saturation indexes of ordinary and critical care beds explained up to 86% of the CFR variability. Therefore, a significant part of the CFR variability may have been caused by the massive stress of the Italian health system. The results are robust across several model specifications. Moreover, cluster analysis showed that the highest mortality risk was concentrated in northern Italy, while the lowest risk was associated with southern provinces.

However, this study also has some limitations that can be summarized as follows: i) first, a significant part of the patients died in hospital presented at least one comorbidity ahead of COVID-19 infection (ISS, 2020); ii) then, as pointed out in other studies (Spalt et al., 2016; Wang et al., 2019), the utilization of air pollution implies unavoidable measurement errors since most people usually stay indoors; iii) third, climatic variables, such as average temperature and DTR, refer just to the average historical values; and iv) finally, the data on prevalence are probably underestimated due to the limited number of swabs carried out in the early stages of the epidemic.

Finally, the study seemed to stress the importance of implementing quick and rational lockdown measures, of making patients comfortable, of implementing an action plan to discourage car use and decrease firm's pollution, and of buying ad hoc health care facilities, medical equipment, and devices to adequately tackle similar and unforeseeable emergencies.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

CRediT authorship contribution statement

Gaetano Perone: Writing - original draft, Writing - review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.