Magnitude and time-course of excess mortality during COVID-19 outbreak: population-based empirical evidence from highly impacted provinces in northern Italy.

BACKGROUND: The real impact of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) on overall mortality remains uncertain as surveillance reports have attributed a limited number of deaths to novel coronavirus disease 2019 (COVID-19) during the outbreak. The aim of this study was to assess the excess mortality during the COVID-19 outbreak in highly impacted areas of northern Italy.
METHODS: We analysed data on deaths that occurred in the first 4 months of 2020 provided by the health protection agencies (HPAs) of Bergamo and Brescia (Lombardy), building a time-series of daily number of deaths and predicting the daily standardised mortality ratio (SMR) and cumulative number of excess deaths through a Poisson generalised additive model of the observed counts in 2020, using 2019 data as a reference.
RESULTS: We estimated that there were 5740 (95% credible set (CS) 5552-5936) excess deaths in the HPA of Bergamo and 3703 (95% CS 3535-3877) in Brescia, corresponding to a 2.55-fold (95% CS 2.50-2.61) and 1.93 (95% CS 1.89-1.98) increase in the number of deaths. The excess death wave started a few days later in Brescia, but the daily estimated SMR peaked at the end of March in both HPAs, roughly 2 weeks after the introduction of lockdown measures, with significantly higher estimates in Bergamo (9.4, 95% CI 9.1-9.7).
CONCLUSION: Excess mortality was significantly higher than that officially attributed to COVID-19, disclosing its hidden burden likely due to indirect effects on the health system. Time-series analyses highlighted the impact of lockdown restrictions, with a lower excess mortality in the HPA where there was a smaller delay between the epidemic outbreak and their enforcement.

Introduction

As the novel coronavirus disease 2019 (COVID-19) has continued to spread internationally since the beginning of the outbreak in December 2019, one of the most relevant public health concerns has been its mortality rate [1, 2].

Italy is among the countries with the highest number of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infections, reaching >240 000 confirmed cases and around 35 000 deaths in late June 2020 [3]. According to the epidemiological bulletin of June 25, 2020 of the Italian National Institute of Health, patients who died of COVID-19 had a median age of 82 years, while the median age of subjects who tested positive for SARS-CoV-2 infection was 62 years [4, 5]. Surveillance data also suggested that the principal risk factors for mortality in COVID-19 patients were older age, male sex and concomitant conditions [4, 6–8].

Current estimates of the case fatality ratio (CFR) vary considerably from one country to another, likely due to differences in age distribution and health status of the population. There is also general consensus that CFR variability may be explained by inaccurate estimates of the number of people who are infected with SARS-CoV-2. Asymptomatic cases of COVID-19, patients with mild symptoms or individuals who have been misdiagnosed could be omitted from the denominator, leading to its underestimation and thus overestimation of the CFR [1, 9]. As a consequence, CFRs do not currently provide reliable figures to assess the real impact of COVID-19 on a country-based mortality burden [10-12].

Moreover, early research assessed that the mortality burden during the first phase of the spread of COVID-19 in Italy is underestimated, because epidemiological surveillance reports have been attributing only a limited number of all deaths to COVID-19 [1, 2, 13, 14], suggesting that the SARS-CoV-2 epidemic had a bigger impact on overall mortality by direct and indirect effects [1, 2]. Direct effects might include a significant proportion of people who would eventually die from the disease before the diagnosis was made; indirect effects might involve the patient's comorbidity status, as well as the disproportionate hospital overload and severe shortage of healthcare resources in the most affected areas of the country [2, 13]. Analyses on mortality data published by the Italian Institute of Statistics (ISTAT) found increased standardised mortality ratios across different areas of the country in the first weeks of the COVID-19 epidemic. This excess mortality was found to be spread across the northern Italian regions [2]. Lombardy registered the biggest COVID-19 outbreak in Italy (and thus in Europe), accounting for >94 000 cases [3], and some of its provinces – including Bergamo and Brescia – were hit hardest in terms of the death toll during the peak epidemic weeks [13, 14].

Scarce evidence is available on excess overall mortality during the outbreak using population-based direct empirical observations and very little has been published on the time-course of this excess.

Based on these considerations, we decided to conduct a retrospective analysis of mortality data from a densely populated vast geographical area of northern Italy in the first 4 months of 2020, with the aim of describing the magnitude and time-course of the entire cycle of the excess mortality associated with the COVID-19 outbreak and investigating the role of lockdown measures in the mitigation of the disease impact on the general mortality.

Methods

Study population

The Lombardy region, the epicentre of the Italian SARS-CoV-2 outbreak, is served by eight local health protection agencies (HPAs), which cover mutually exclusive areas of the region.

We focused on all inhabitants of the HPAs of Bergamo and Brescia, where the outbreak magnitude has been particularly intense [13, 14].

Study design and data sources

Italy has a tax-funded, universal coverage National Health System that is organised at three levels: national, regional and local. In Lombardy, the local HPAs are primarily responsible for collecting healthcare information (e.g. inpatient and outpatient treatments), including the mortality data of all inhabitants under their jurisdiction [15].

We carried out a longitudinal retrospective time-series study on the overall mortality traced in these standardised healthcare regional administrative databases.

To accomplish the study objectives, we retrieved data on all-cause deaths that occurred between 1 January and 30 April in the years 2020 and 2019, the latter to be used as a reference period. For each deceased subject we linked the unique identification code in the mortality database with demographic information (i.e. age at death and sex).

We chose 2019 as a reference, because we observed that the seasonality of the flu outbreak, which is one of the drivers of variations in mortality during winter, was similar to that observed in 2020 up to the beginning of the COVID outbreak [16]; during previous years, peaks occurred in different periods. Therefore, we hypothesise that 2019 approximates the expected situation in 2020.

Statistical analysis

We compared the age and sex distribution among the deceased in 2019 and 2020 using Chi-squared tests and the number of daily deaths using the Wilcoxon rank-sum test.

We compiled the time-series of the number of observed daily deaths during January to April in 2019 and 2020 for each HPA. Based on these time-series, and considering that the population of each HPA remained stable from 2019 to 2020 (HPA Brescia: from 1 089 602 to 1 097 648; HPA Bergamo: from 1 160 374 to 1 163 243), we computed the expected daily number of deaths for 2020 as the average daily number of deaths that occurred in 2019.

In order to assess the magnitude and time-course of the mortality variation during the SARS-CoV-2 outbreak, we built a Poisson generalised additive model, with the daily death count during 2020 as the outcome, and the HPA, a non-linear function of the date and the number of expected deaths as the predictors.

In detail, the model was:where is the observed number of deaths during the i-th day in j-th HPA, is the daily expected number of deaths in the j-th HPA, and it is held fixed for each i-th day, is a categorical variable indicating the HPA where the number of deaths was measured in the i-th day, and is a thin plate regression spline of the calendar date [17]. The term indicates that a separate spline is fitted for each HPA, allowing for the estimation of a different epidemic time-course in the two areas. The use of the constant as an offset allows direct modelling of the daily observed/expected ratio of mortality, that is the standardised mortality ratio (SMR), instead of the daily counts, predicting both the estimated number of daily deaths and the estimated SMR [18, 19].

Therefore, based on the model results, we first predicted the SMR for each day of 2020 together with its 95% confidence interval. Thereafter, using the method reported by Rivera et al. [20], we estimated the cumulative number of excess deaths, as compared to 2019, and their correlated 95% credible sets (CS), which similarly to confidence interval synthesise the precision of the estimates but based on Bayesian methods. Finally, we expressed the excess in term of an n-fold increase in the cumulative number of deaths, using 2019 cumulative number of deaths as the reference.

We replicated the whole analysis stratifying by sex and by age, categorised as: <60 years, 60–69 years, 70–79 years, 80–89 years and 90+ years.

All analyses were carried out with statistical software SAS version 9.4 (SAS Institute, Cary, NC, USA) and R version 4.0.0 (R Project for Statistical Computing, www.R-project.org).

Results

During the first 4 months of 2020, 17 099 deaths were registered among the residents of the two HPAs of Bergamo and Brescia, as compared to 7592 in 2019 (table 1). In both provinces, the crude number of deaths increased by >90%. Overall, the daily average number of deaths rose significantly from 63 (sd: 10.7) to 141 (sd: 129.8), with a steeper wave of the cumulative time-series analysis in the area of Bergamo (figure 1).

TABLE 1

Characteristics of the deceased over the period 1 January to 30 April 2020, stratified by health protection agency (HPA); comparison between 2020 and 2019

	HPA Bergamo		HPA Brescia		Total
	2019	2020	2019	2020	2019	2020
Deaths n	3662	9433	3930	7666	7592	17 099
Sex
Female	1920 (52.4)	4407 (46.7)	2082 (53.0)	3823 (49.9)	4002 (52.7)	8230 (48.1)
Male	1742 (47.6)	5026 (53.3)	1848 (47.0)	3843 (50.1)	3590 (47.3)	8869 (51.9)
p-value versus 2019		<0.001		0.002		<0.001
Age-class years
<60	281 (7.7)	440 (4.7)	242 (6.2)	357 (4.7)	523 (6.9)	797 (4.7)
60–69	321 (8.8)	873 (9.3)	313 (8.0)	634 (8.3)	634 (8.4)	1507 (8.8)
70–79	741 (20.2)	2267 (24.0)	786 (20.0)	1713 (22.3)	1527 (20.1)	3980 (23.3)
80–89	1399 (38.2)	3809 (40.4)	1543 (39.3)	3019 (39.4)	2942 (38.8)	6828 (39.9)
90+	920 (25.1)	2044 (21.7)	1046 (26.6)	1943 (25.3)	1966 (25.9)	3987 (23.3)
p-value versus 2019		<0.001		<0.001		<0.001
Daily number of deaths
Mean±sd	30.5±6.39	78.0±82.75	32.8±6.94	63.4±49.49	63.3±10.65	141.3±129.80
p-value versus 2019		<0.001		<0.001		<0.001

Data are presented as n (%), unless otherwise stated.

FIGURE 1

Comparison of 2020 and 2019 cumulative daily counts of death from 1 January to 30 April, stratified by health protection agency (HPA).

The trends in the estimated daily SMR confirm the substantial differences in the two HPAs for both sexes and for all age-classes older than 60 years (figure 2). The overall daily estimated SMR peaked at 9.4 (95% CI 9.1–9.7) in Bergamo province and at 5.6 (95% CI 5.4–5.8) in Brescia. The highest SMRs were observed among the elderly, peaking in the age-class 70–79: 12.7 (95% CI 11.9–13.5) in Bergamo; 7.3 (95% CI 6.8–7.9) in Brescia.

FIGURE 2

Predicted standardised mortality ratio (SMR) from 1 January to 30 April 2020. Results from the overall analysis and from analyses stratified by age or sex, by health protection agency (HPA).

Overall, during the SARS-CoV-2 outbreak we estimated 5740 (95% CS 5552–5936) excess deaths in the area of Bergamo and 3703 (95% CS 3535–3877) in Brescia, corresponding to a 2.55-fold (95% CS 2.50–2.61) and 1.93-fold (95% CS 1.89–1.98) increase in the number of cumulative deaths (table 2; figure 3).

TABLE 2

Estimated cumulative excess deaths and relative increase as compared to 2019, over the period 1 January to 30 April 2020, stratified by health protection agency (HPA), age or sex

Disease	HPA Bergamo		HPA Brescia
	Estimated excess number of deaths (95% CS)	Estimated increase (n-fold) in the number of deaths (95% CS)	Estimated excess number of deaths (95% CS)	Estimated increase (n-fold) in the number of deaths (95% CS)
Total	5740 (5552–5936)	2.55 (2.50–2.61)	3703 (3535–3877)	1.93 (1.89–1.98)
Sex
Female	2471 (2346–2601)	2.28 (2.21–2.34)	1724 (1609–1851)	1.82 (1.77–1.88)
Male	3269 (3134–3414)	2.86 (2.78–2.94)	1980 (1862–2103)	2.06 (2.00–2.13)
Age-class years
<60	157 (120–203)	1.55 (1.42–1.72)	113 (80–155)	1.46 (1.33–1.64)
60–69	549 (496–612)	2.70 (2.53–2.89)	318 (274–374)	2.01 (1.87–2.18)
70–79	1520 (1429–1618)	3.03 (2.91–3.17)	920 (845–1009)	2.16 (2.07–2.27)
80–89	2398 (2283–2525)	2.70 (2.62–2.79)	1463 (1361–1577)	1.94 (1.87–2.01)
90+	1116 (1033–1212)	2.20 (2.11–2.31)	888 (807–981)	1.84 (1.76–1.93)

The summation of the estimated number of deaths within sex or age categories might not equal the total estimated excess, since stratified estimates are derived from separate models.

FIGURE 3

Estimated cumulative excess deaths and percentage increase compared to 2019, over the period 1 January to 30 April 2020, stratified by health protection agency (HPA) and age or sex.

Restricting the analyses for the months of March and April, the estimated number of excess deaths was 5719 (95% CS 5556–5892) in Bergamo and 3820 (95% CS 3676–3967) in Brescia, representing virtually the whole excess mortality of the entire 2020 analysed period (table 3). This is confirmed by crude rates (supplementary table S1): observed rates over the months of January and February during 2020 were comparable to those of 2019, if not slightly lower. The differential increase across sexes and age-classes was similar in the two provinces, but higher in Bergamo: in the age-class 70–79, the number of cumulative deaths increased by 5.05-fold (95% CS 4.83–5.28), compared to 2.18-fold (95% CS 1.96–2.44) among the youngest.

TABLE 3

Estimated cumulative excess deaths and relative increase as compared to 2019, over the period 1 March to 30 April 2020, stratified by health protection agency (HPA), age or sex

Disease	HPA Bergamo		HPA Brescia
	Estimated excess number of deaths (95% CS)	Estimated increase (n-fold) in the number of deaths (95% CS)	Estimated excess number of deaths (95% CS)	Estimated increase (n-fold) in the number of deaths (95% CS)
Total	5719 (5556–5892)	4.07 (3.98–4.17)	3820 (3676–3967)	2.91 (2.84–2.99)
Sex
Female	2422 (2314–2536)	3.48 (3.37–3.60)	1787 (1686–1894)	2.69 (2.59–2.79)
Male	3296 (3176–3426)	4.72 (4.59–4.87)	2033 (1930–2140)	3.16 (3.05–3.28)
Age-class
<60	168 (137–206)	2.18 (1.96–2.44)	108 (81–139)	1.87 (1.66–2.13)
60–69	538 (489–593)	4.30 (3.99–4.63)	334 (294–380)	3.10 (2.85–3.39)
70–79	1525 (1442–1613)	5.05 (4.83–5.28)	982 (913–1058)	3.46 (3.29–3.65)
80–89	2377 (2271–2487)	4.34 (4.19–4.50)	1530 (1439–1625)	2.95 (2.83–3.07)
90+	1108 (1034–1187)	3.37 (3.21–3.54)	866 (796–941)	2.63 (2.50–2.77)

The sum of the estimated number of deaths within sex or age categories might not equal the total estimated excess, since stratified estimates are derived from separate models.

Observing the time-series of the daily number of deaths (figure 4), the excess death wave started earlier in the area of Bergamo, roughly 1 week after the detection of the first COVID-19 case in Lombardy (February 21, 2020), and a few days later in the area of Brescia. In both areas, the excess death wave started before the introduction of lockdown measures (March 7, 2020), and the peaks were observed between March 20 and March 25, 2020, roughly 2 weeks after the introduction of lockdown measures, returning to pre-epidemic levels by the end of April.

FIGURE 4

Time-series of daily number of deaths, stratified by health protection agency (HPA). a) 21 February: first confirmed case in Lombardy; b) 22 February: lockdown and quarantine for 11 municipalities; c) 25 February: schools closed and crowd restrictions in six northern regions; d) 4 March: schools closed and crowd restrictions extended to the whole country, e) 7 March: lockdown and quarantine for Lombardy and 14 provinces, f) 11 March: lockdown and quarantine extended to the whole country.

Discussion

This retrospective observational study, based upon the complete data of two Italian HPAs where the outbreak was particularly intense, evaluated the temporal trend of all-cause mortality during the COVID-19 outbreak, describing its complete cycle using real world direct empirical observations.

During the 4-month period of study, we estimated 9443 excess deaths in the areas covered by the HPAs of Bergamo and Brescia, corresponding to a 2.2-fold increase in number of excess deaths as compared to those recorded in the same period of 2019, with a trend that showed a rapid rise starting 10 days after the first case of COVID-19 in the region. These results exceed by far the number of deaths officially attributed to the disease in the two HPAs, which registered 2973 and 2295 deaths in Bergamo and Brescia, respectively, as of 30 April 2020.

It is unlikely that the large increase in number of deaths is solely attributable to conditions other than COVID-19, which affected mortality directly and potentially also indirectly, e.g. through healthcare facilities being overwhelmed and delayed care for time-dependent conditions, such as stroke, myocardial infarction, etc. [2, 21]. Indeed, we observed that all excess deaths occurred in March and April, after the beginning of the disease outbreak, in line with preliminary findings from March only that the overall mortality did not increase in Italian provinces with low rates of infection [13].

Even though the two HPAs are adjacent, the spike in the daily number of deaths started a few days earlier in Bergamo than in Brescia, ∼2 weeks before the introduction of lockdown restrictions. These few days may have been crucial: although the waves had similar shapes, the province of Bergamo showed a significantly and substantially higher increase in mortality (2.55-fold) compared with Brescia (1.93-fold). These findings suggest that the national and local restrictions may have had a massive impact, resulting in a less steep wave in the time-series of the daily number of deaths in Brescia, where virus circulation was delayed. Moreover, in both HPAs we observed significant excess in mortality for 8 weeks, with a peak 4 weeks after the start of the epidemic wave, and an overlapping descending phase. These trends indicated that the length of the excess death period may not be dependent only on the circulation of infection, which is correlated to overall SMR, suggesting that other factors might have influenced the overall mortality in both areas. One hypothesis can be identified in the emergency measures taken during the epidemic phase, such as activation of Special Continuity Assistance Units, a primary care medical home service dedicated to COVID-19 patients, or the construction of field hospitals, which lessened the overload experienced by healthcare facilities in both areas [22, 23].

Confirming previous findings, we observed an increased risk of death among males and the elderly [6–8, 24]. Overall, the 70–79 years age group was the most affected in terms of cumulative deaths, although in COVID-19 patients the highest risk of death was observed among those aged 80 years or more [4]. The discrepancy might be due to a lower infection rate among the oldest, the mobility and social contacts of whom are frequently reduced. Noteworthy, we found a weaker excess mortality in nonagenarians than in younger elderly groups, possibly explained by the healthy survivor effect as well as a less clear-cut direct and indirect impact of COVID-19 in those age-classes in which the underlying functional status is a more significant predictor of all-cause mortality than acute life-threatening conditions [25]. Moreover, demographic characteristics of the population may have contributed to the shape of the curves, with the distribution of the elderly Italian population posing a challenge in reducing the impact of mortality due to the epidemic.

The strengths of the study include the use of administrative databases within a universal coverage system on an unselected population of residents with complete knowledge of vital status, within an area where the outbreak has been particularly strong, and the virus circulated long before the enforcement of lockdown measures. Also, the analysis did not focus only on deaths classified as due to COVID-19, thereby avoiding underestimating the number of deaths at the beginning of the epidemic, when only a selected sample of patients underwent COVID-19 testing. Indeed, our estimates provide inferential estimates of excess deaths and SMRs, based on flexible models that could be easily reproducible in other geographical areas. This is the first study describing the complete cycle of the outbreak, from the first cases to the baseline risk of death using regression. This also allows the comparison of the excess death waves between two geographical areas, providing an initial description of how the timing of lockdown measures can ultimately affect overall mortality. Finally, this research expressed robust data for evaluating the effectiveness of the epidemic response and for providing measures to inform policymakers on ending restrictive lockdown.

Some limitations should also be considered when interpreting the results. First, while accuracy is preserved, the precision of the estimated excess deaths is affected by the average number of daily events, decreasing in the subgroups that have low daily frequencies of death (e.g., the youngest age-class). Second, the reference level is computed based on 2019 only: using a longer time-window could define the reference level better. However, trends in mortality within these areas proved to be very stable during the last few years [26].

In conclusion, we documented a significant increase in the overall mortality during the first months of 2020, particularly March and April, indicating that the COVID-19 outbreak had a substantially larger impact than official estimates. Time-series analyses suggest that the national and local restrictions had a massive effect, providing a considerable reduction of COVID-19 burden. Furthermore, this study may serve as a model for country-based estimations of the overall (direct and indirect) impact of the COVID-19 epidemic on population mortality.

Please note: supplementary material is not edited by the Editorial Office, and is uploaded as it has been supplied by the author.

TABLE S1 Crude rates over the period January 1st to April 30th 2020 in the two HPAs, overall and stratified by sex and/or age class 00458-2020.tableS1