Estimating weekly excess mortality at sub-national level in Italy during the COVID-19 pandemic

Study design

We used a spatio-temporal disease mapping approach to evaluate the excess mortality in Italy at municipality level, while detecting and predicting its evolution on a weekly basis. First, we modelled weekly mortality trends over the years 2016 – 2019, accounting for air temperature. We used these trends to predict the weekly- and municipality-specific mortality over the period 1 January to 28 April 2020. Then, we compared the observed mortality for this period with the model-based predictions. This allows us to quantify the excess, defined as the number of deaths from all causes relative to what it would have been expected in the absence of the COVID-19 outbreak, based on the model.

We performed separate analyses for males and females, as previous studies have found differences in mortality between sexes. [12,13] The main results are presented for the total population; we adjust for the age structure at municipality level and across the study period through an internal standardisation. [14] We calculated age-specific rates (0-14, 15-24, 25-34, 35-44, 45-54, 55-64, 65-74, 75+) across the country for the period under study in 2016 – 2019 as these represent the all-cause mortality in years without the pandemic. We then applied the age-specific rates to the year-and age-specific population in each municipality to estimate the expected number of cases as denominator.

Data

We used official weekly mortality data, available from the Italian Institute of Statistics (ISTAT; source: https://www.istat.it/en/archivio/240106). The data included counts by age, sex and municipality of residence of the cases. All-cause mortality cases were recorded for the 7,904 municipalities that comprise the whole Italian territory. Administratively, the Italian municipalities are nested into 107 provinces, which are themselves grouped into 20 regions.

We used the weekly deaths in 2020, available from the same source, as a comparison with the values predicted from our model. Currently the 2020 data are only available on 7,251 municipalities, which cover 91.7% of the Italian territory and 93.4% of the Italian population; hence, we limit the comparison to this subset.

To capture the seasonal variability in death counts, we also included air temperature data (at 2 meters height), obtained from Copernicus ERA5 global weather and climate reanalysis dataset, [15] which combines a weather model with observational data from satellites and ground monitors to produce global hourly data at a 30 km grid spatial resolution. Operationally, we processed air temperature data within the Google Earth Engine cloud-based platform to obtain weekly series of air temperature for each Italian municipality.

Municipality weekly trends for 2016 – 2019

In line with the approach typically considered for disease mapping, [16] we modelled the total number of deaths at the municipality level, separately for males and females, for each week and year using a Poisson distribution with age-adjusted expected number of cases as offset. We specified a Bayesian hierarchical model on the log mortality relative risk, which is a common approach to overcome the high variability in the estimates driven by the small numbers of cases, due to the high spatio-temporal resolution. [17] In order to capture changes in space, we assumed that the log relative risk for each municipality is smoothed locally through a conditional autoregressive structure, [18] which models similarity across neighboring municipalities. To capture the long-term temporal trend, we considered a year-specific random intercept and accounted for the potential non-linearity in the mortality relative risks through a weekly smooth effect, where each week is assumed to depend on the value of the previous one (order one random walk). As the temporal trend might be different in space, we allowed a different random walk for each province. Finally, as there is extensive evidence of a non-linear effect of temperature on mortality, [19,20] we included a smoothed effect on weekly temperature at municipality level, through a second order random walk, so that each week is influenced by the previous two.

The model has been estimated using Integrated Nested Laplace Approximation and the software R-INLA (www.r-inla.org). A full specification of the model structure and the underlying distributional assumptions and priors are presented in details in the Supplementary material.

Weekly death prediction for 2020

We used the output from the model for 2016 – 2019 to predict the number of deaths expected for each week of our follow up period. Specifically, for each municipality and week, we sampled 1000 values from the posterior predictive distribution, where the log relative risk is the linear combination of the corresponding spatial, temporal and temperature effects estimated given the data from 2016 – 2019. To report results at the province (regional) level, we summed each draw across the municipalities belonging to the same province (region). This simulation-based approach to the post-processing of the model output allows to fully characterise the uncertainty in the mortality rates and to propagate it in a principled way through to the predicted mortality counts.

We present plots of observed and predicted trends in mortality rates as well as maps and plots of percentage excess mortality, estimated as the difference between the observed number of deaths and the predicted number under current trends (based on past data and on the assumption that the pandemic did not occur), divided by the observed number.

Patients and Public Involvement

It was not appropriate or possible to involve patients or the public in the design, or conduct, or reporting, or dissemination plans of our research

Statement of principal findings

In this paper we presented the first fully probabilistic analysis of the excess mortality at weekly sub- national level in Italy, one of countries with the highest number of deaths during the COVID-19 pandemic. We found that 2020 starts with lower rates than in previous years [13], but from the beginning of March 2020 there is evidence of an excess mortality compared with the trend in 2016 – 2019. However, our model suggests large geographical differences in the excess mortality across Italy; while it is now widely established that Lombardia is the most affected region (and one of the hardest hit, globally), we estimated large heterogeneity even with the closest neighbours. For instance, the city of Verona in Veneto, one of the regions in the North East, which borders Brescia (one of the hot-spots in Lombardia) was associated with mortality rates essentially in line with the historical trend. Similarly, Pesaro (in the central region of Marche) is associated with excess mortality rates that exceed those of Milan (the capital of Lombardia), while belonging to a region that has by and large not been severely affected. We compared our estimated excess mortality with the official statistics reported by ISTAT, which distinguish COVID-19 registered deaths from other causes. While it is impossible, nor is it our objective, to draw causal conclusions from the combination of these data, our analyses suggest a potential large increase in mortality, after discounting the confirmed COVID-19 deaths, particularly in Lombardia, but also in the other regions of Northern Italy. On the other hand, the Centre and the South of Italy do not show evidence of increased net mortality.

Strengths and weaknesses of the study

The specification of the model at municipality level, coupled with the inclusion of weekly non-linear terms, which are allowed to vary for each province, ensures that we can detect heterogeneity across space, while still retaining the flexibility to model the temporal pattern of mortality. This is consistent with typical modelling in disease mapping [22] and surveillance studies [23]. Additionally, we are framed in the same perspective as forecasting studies [24,25] and we take advantage of the Bayesian nature of our model to characterise the full uncertainty over the estimate of the rates to predict the weekly trends in 2020. At the same time, as our analysis estimates temporal trends for each municipality with the aim of predicting 2020 under the alternative scenario of the absence of a pandemic, we do not need to explicitly include covariates only varying in space. This is because their effect is captured by the municipality effect, which estimates the spatial heterogeneity in the rates. Nevertheless, we included weekly mean air temperature for each municipality to adjust the trends in mortality, which is crucial as 2020 was generally warmer than the previous years.

Our model has some limitations. At the time of writing, cause-specific mortality for the year 2020 was not available. Thus, while we were able to estimate the non-COVID-19-related mortality, we cannot disentangle the specific causes behind it. At the same time, we were not able to account for the potential reduction in the number of deaths decreasing during lockdown, due for example to a decrease of traffic accidents and injuries at work. This would lead to a conservative estimate of the excess mortality in our model. However, these causes generally account for less than 1% of the total deaths. Our results show that the excess mortality more than offsets that. Additionally, there have been several reports which have showed improvement in air quality during the lockdown in different countries, [26,27] leading to a consistent reduction of the pollution-related mortality. [28,29] However, the lack of real-time air pollution data covering the entire study region meant that we could not consider this aspect in the analyses.

Strengths and weaknesses in relation to other studies, discussing important differences in results

While an early comparative study on COVID-19-related mortality suggested that the region Lombardia was not heavily impacted by the epidemic wave, [30] in accordance with more recent studies, [31,32] we found that it was the most affected region, for both sexes. Other studies have tried to assess the impact of the pandemic using a lower geographical resolution, e.g. by considering national or regional data [9,10]. Conversely, our model is specified at the finest administrative level in Italy, which allows us to detect a very large heterogeneity both in space and time, in terms of the spread and the resulting impact of the pandemic on the health-care services. This would not have been picked up by an analysis at regional level.

Meaning of the study: possible explanations and implications for clinicians and policymakers

From the analysis of our results, it appears that the lockdown measures implemented by the Italian government have managed, over the course of several weeks, to contain the excess mortality, as evidenced by the steep decrease during the month of April. However, this is characterised by large spatial heterogeneity: while the Centre and the South return to mortality rates in line with the previous years, the regions of the North-West, North-East and, especially, Lombardia are still showing excess mortality by the end of our observation period. In the most affected region, males mortality rates by the end April return to levels normally observed at the beginning of January, while females rates remain still higher than expected. This suggests perhaps the necessity of staggered relaxation of the lockdown, in order to prevent a second wave and continue to limit the burden on the health-care system. As suggested elsewhere, [21] the aggressive programme of mass-testing implemented in Veneto is potentially a driving factor in the much lower excess mortality experienced by the region, in comparison to neighbouring areas such as Lombardia and, albeit to a lower extent, Trentino-Alto Adige.

The first confirmed COVID-19 death in Italy was on 21 of February and from the model there is evidence of a deviation from the expected temporal trend since the last week of February, with the national lockdown implemented on 9 of March. This suggests that our model is able to track the changes in the historical dynamics, thus making our framework particularly valuable to aid policy-making: as new data become available, the model can be used to identify when the rates return to the normal range, but also if and where they start deviating again due to the emergence of a potential second wave requiring for instance social restrictions targeted in space. Additionally, the modelling framework is a useful tool for prospective real-time surveillance, as it allows to effectively follow each municipality in time and flag unusual behaviours as soon as they happen, aiming at identifying potential new public health threats and contain them.

Unanswered questions and future research

In addition to providing estimates for the specific Italian case, the proposed modelling framework is highly scalable and adaptable. A natural extension of our work consists in the analysis of cause-specific mortality, when the data become available: this will lead to a better understanding of the trends in COVID- 19-related mortality, with respect to other competing mortality risks. Moreover, it would be possible to replicate almost without changes the model for other countries, once data become available at a relevant spatial resolution (e.g. local authorities, or even middle-super output areas in the UK).