Estimating COVID-19-induced excess mortality in Lombardy, Italy.

We compare the expected all-cause mortality with the observed one for different age classes during the pandemic in Lombardy, which was the epicenter of the epidemic in Italy. The first case in Italy was found in Lombardy in early 2020, and the first wave was mainly centered in Lombardy. The other three waves, in Autumn 2020, March 2021 and Summer 2021 are also characterized by a high number of cases in absolute terms. A generalized linear mixed model is introduced to model weekly mortality from 2011 to 2019, taking into account seasonal patterns and year-specific trends. Based on the 2019 year-specific conditional best linear unbiased predictions, a significant excess of mortality is estimated in 2020, leading to approximately 35000 more deaths than expected, mainly arising during the first wave. In 2021, instead, the excess mortality is not significantly different from zero, for the 85+ and 15-64 age classes, and significant reductions with respect to the 2020 estimated excess mortality are estimated for other age classes.

Introduction

Official publicly available data mainly report incidence indicators. Incidence indicators measure the number of individuals with a particular condition, related with the epidemic, recorded during a given period. They can be referred to different time periods; in particular, in the Italian Civil Protection Department dataset (https://github.com/pcm-dpc/COVID-19), daily incidence counts are available for the following indicatorsThe number of positives available from official reporting, however, underestimates the true number of cases since there exists a vast proportion of asymptomatic or mildly symptomatic patients, among all infected individuals, who are not detected [1-4]. Complex methods are then required to provide reasonable forecasts of the epidemic evolution [5-8] and avoid unreliable predictions, which make people worry (on this point, see e.g. [9, 10]). Thus, the attention has recently moved to the analysis of mortality, as the gold standard measure of the impact of COVID-19 worldwide [11, 12]. This is mainly because mortality data suffer less of underreporting, measurement errors and other data issues, i.e. are in general more reliable. We would like to contribute to this growing literature by introducing a generalized linear mixed model to estimate COVID-19-induced excess mortality in Lombardy, Italy.

positives, which are classified into

hospitalized (either in regular wards or in ICU)

isolated-at-home

tested

deceased

recovered/discharged

Italy was dramatically hit by the COVID-19 pandemic, with 4.8 million cases and 132000 deaths as of October 31st, 2021. Lombardy, the most populated and industrialized Italian region, was the epicenter of the outbreak, in particular during the first wave in early 2020. It experienced one of highest case-fatality rates with nearly 900,000 cases (about one fifth of those registered in Italy) and nearly 34,000 deaths out of 10 million inhabitants (for a discussion see e.g. [13]). The impact of COVID-19 on mortality rates is rather clear and can be easily depicted by a simple data visualization. Quantifying the excess deaths induced by COVID-19 is more difficult and may reveal both direct and indirect effects of the pandemic on mortality. Statistical methods have been recently introduced ranging from simple averaging methods [14] to structured models [15-17].

Here we introduce a model to estimate excess mortality in 2020 and 2021 in Lombardy stratified by various age classes. The proposed generalized linear mixed model, in the spirit of [18], is able to capture year-specific trends, seasonality and overdispersion in the number of all-cause deaths and is based on historical data collected on a weekly basis from 2011 to 2019.

Data and methods

Open-source daily all-cause Italian mortality data, stratified by Region and age classes, from 2011 to 2021 are available from the Italian National Institute of Statistics (ISTAT, https://www.istat.it/it/archivio/240401). These data were downloaded on October 25th 2021 and temporally aggregated at the week level. The weeks number 53, present in some years, were dropped from the analysis. Data on historical all-cause mortality are displayed in Fig. 1. Seasonality is clearly present, and is shown for all the considered years. Peaks are rather evident at the major COVID-19 waves, with a maximum of 7729 weekly deaths at the end of March 2020, during the first pandemic wave. Similar behaviors, on different scales, are observed for age-specific data, not reported here for the sake of brevity. We report age-specific plots in the Supplementary Material.

Fig. 1

Observed total number of deaths per week (solid black line) and model's estimated values (thick solid blue line). The dashed (red) line shows the expected weekly deaths under the historical model for years 2020 and 2021

To better appreciate age-class peculiarities, different models are fitted to observed deaths by age classes. Five age classes are considered here, namely (0–14), (15–64), (65–74), (75–84), 85+. This is very useful to detect if and how COVID-19 has an impact on mortality at different age classes of the population and to plan interventions to protect at-risk sub-populations.

The crucial aspect in estimating excess mortality is the definition of a reliable benchmark mortality model, i.e. a model providing good estimation of the expected mortality under pre-epidemic conditions. We model the weekly mortality at week and years . The count data process on weekly deaths is modeled with a Negative Binomial distribution, i.e.withCount data may be distributed as Negative Binomial if the rate, at which events occur, is heterogeneous, and consequently the counts are characterised by overdispersion compared to the Poisson (as typically happens in our deaths time series). The Poisson distribution is nested within the Negative Binomial, in the sense that if no overdispersion/heterogeneity is present, the Negative Binomial distribution converges to the Poisson distribution.

In count data analysis, the interest is usually focused upon the parameter vector which is modelled, in a regression context, by defining a generalized linear mixed model with log link for the analysed response, i.e.Cyclical patterns, with-year fluctuations, are captured by Fourier series, whose term I is defined by using model selection criteria; to account for year-specific mortality baselines and to accommodate, at least partially, auto-correlation of the weekly counts, random effects are introduced. Model’s fitting to the data is shown in Fig. 1 (see the thick solid blue line). The weekly predictions of mortality data for years 2020 and 2021 are based on the 2019 year-specific conditional best linear unbiased predictions of the generalized linear mixed model. Excess mortality is obtained by subtracting the expected deaths based on the pre-pandemic period from the registered all cause deaths in the pandemic period. Finally, the weekly estimates of excess mortality are summed over the year (or part of it, for 2021).

Results and discussion

Results are shown in Table 1. Prediction intervals are found through parametric bootstrap; if zero is included in the intervals, no difference with the expected number of deaths is estimated, i.e. there is no excess mortality. The first interesting result is that no excess mortality is estimated in Lombardy neither in 2020 nor in the first 34 weeks of 2021 for the youngest sub-population aged (0–14); indeed, a significant reduction in the number of deaths is estimated for this specific sub-population in 2021. This reinforces the idea that children did not experience excess mortality due to COVID-19 infections. For all the other classes, significant excesses are estimated in 2020, leading to approximately 35,000 more deaths than expected. The largest excess of deaths, in relative terms, is in the classes 65–74 and 75–84, with an excess of approximately 28%. A 25% excess is estimated in the oldest class 85+. For the working age class 15–64, instead, the estimated excess mortality is approximately 18%. There is, thus, a clear age effect on the risk of death induced by COVID-19. As a by-product of the analysis, we are able to detect weeks where the estimated excess mortality is relevant. In Lombardy, a large part of the excess of mortality is estimated during the first wave in March–April 2020, though similar behaviors, with smaller effects, are estimated in Autumn 2020.

Table 1

Estimated excess mortality in Lombardy, with 95% prediction intervals in brackets

Age class	Observed deaths2020	Excess mortality2020
0–14	275	− 33 (− 96; 18)
15–64	12260	2019 (1575; 2440)
65–74	18150	5090 (4512; 5639)
75–84	41275	11568 (10604; 12487)
85+	63260	15920 (12407; 19226)
	Observed	Excess	Observed	Excess
	Deaths	Mortality	Deaths	Mortality
	Weeks 1 to 34	Weeks 1 to 34	Weeks 1 to 34	Weeks 1 to 34
	2020	2020	2021	2021
0–14	182	− 20 (− 61; 13)	146	− 56 (− 96; -24)
15–64	8275	1525 (1232; 1802)	7043	293 (− 1; 572)
65–74	12384	3765 (3384; 4126)	10045	1425 (1047; 1789)
75–84	28253	8594 (7953; 9197)	21315	1655 (1032; 2264)
85+	43177	11864 (9521; 14093)	32004	691 (− 1573; 2966)

According to our analysis, the population 65+ in Lombardy has the highest excess mortality induced by the COVID-19 epidemic. This justifies the priority assigned to the elderly during the vaccination process, as they showed the highest mortality both in absolute and relative terms. As we do not have data for the entire 2021, we focus on the first 34 weeks (up to the end of August), for which we have complete data. Excess mortality is not significantly different from zero for the 85+ and 15–64 age classes, and the estimated excess mortality shows significant reductions with respect to the 2020 estimated values for other age classes. The 65–74 and 75–84 classes are the ones more at risk, with 14.2% and 7.8% of excess mortality estimated; this is somehow not surprising because these sub-populations were the latest to be vaccinated among the elderly (for the relevant data, please, refer to https://github.com/italia/covid19-opendata-vaccini). The positive effect on excess mortality reduction of the vaccine can be better appreciated if we compare excess mortality in the first 34 weeks of 2020 and 2021. A strong reduction in the excess mortality induced by COVID-19 is estimated, even in age classes showing a significant excess.We further discuss results for all age classes except the (0–14) on, which has no excess mortality estimated, focusing on specific weeks where excess mortality arose. In 2020, we estimated an excess mortality on all other classes starting from mid-March to mid-May and, then, from the beginning of November to the end of the year, with peaks in weeks 45–48. Our model does not estimate any excess mortality at the beginning of 2021 for all age classes but the (75–84) one, whose estimated excess mortality on weeks 2 to 4 is statistically significant. The excess mortality in 2021 is mainly concentrated in some specific weeks, from mid-March to mid-May, with small time-scale differences across (65–74) and (75–84) age classes. In 2021, instead, the vaccine limits the rise of complications possibly leading to death. This result is even more robust as non-pharmaceutical interventions, as lockdowns, are definitely less restrictive during 2021 than during 2020.

Below is the link to the electronic supplementary material.

Supplementary file1 (PDF 11 KB)

Supplementary file2 (PDF 11 KB)

Supplementary file3 (PDF 11 KB)

Supplementary file4 (PDF 381 KB)

Supplementary file5 (PDF 10 KB)