key: cord-0999466-owdrzv2z
authors: Blangiardo, M.; Cameletti, M.; Pirani, M.; Corsetti, G.; Battagli, M.; Baio, G.
title: Estimating weekly excess mortality at subnational level in Italy during the COVID-19 pandemic
date: 2020-06-09
journal: nan
DOI: 10.1101/2020.06.08.20125211
sha: fa28c65c049398eba8e4713618c6a24fa282748c
doc_id: 999466
cord_uid: owdrzv2z

Background Excess mortality from all-cause has been estimated at national level for different countries, to provide a picture of the total burden of the COVID-19 pandemic. Nevertheless, there have been no attempts at modelling it at high spatial resolution, needed to understand geographical differences in the mortality patterns, to evaluate temporal lags and to plan for future waves of the pandemic. Methods This is the first subnational study on excess mortality during the COVID-19 pandemic in Italy, the third most-hit country. We considered municipality level and estimated all-cause mortality weekly trends based on the first four months of 2016 -- 2019. We specified a Bayesian hierarchical model allowing for spatial heterogeneity as well as for non-linear smooth spatio-temporal terms. We predicted the weekly mortality rates at municipality level for 2020 based on the modelled spatio-temporal trends (i.e.~in the absence of the pandemic) and estimated the excess mortality and the uncertainty surrounding it. Results We found strong evidence of excess mortality for Northern Italy, with higher mortality rates than expected from the end of February in Lombardia, with total excess deaths of 23,946 (23,013 -- 24,786), and the beginning of March for North East and North West with total excess deaths of 8,033 (7,061 -- 9,044) and 1,588 (404 -- 2,700) respectively. We found marked geographical differences, with percent excess of up to 88.9% (81.9% -- 95.2%) at the peak of the pandemic, in the city of Bergamo (Lombardia).

Excess mortality from all-cause has been estimated at national level for different countries, to provide a picture of the total burden of the COVID-19 pandemic. Nevertheless, there have been no attempts at modelling it at high spatial resolution, needed to understand geographical differences in the mortality patterns, to evaluate temporal lags and to plan for future waves of the pandemic.

This is the first subnational study on excess mortality during the COVID-19 pandemic in Italy, the third most-hit country. We considered municipality level and estimated all-cause mortality weekly trends based on the first four months of 2016 -2019. We specified a Bayesian hierarchical model allowing for spatial heterogeneity as well as for non-linear smooth spatio-temporal terms. We predicted the weekly mortality rates at municipality level for 2020 based on the modelled spatio-temporal trends (i.e. in the absence of the pandemic) and estimated the excess mortality and the uncertainty surrounding it.

We found strong evidence of excess mortality for Northern Italy, with higher mortality rates than expected from the end of February in Lombardia, with total excess deaths of 23,946 (23,013 -24,786) , and the beginning of March for North East and North West with total excess deaths of 8,033 (7,061 -9,044) and 1,588 (404 -2,700) respectively. We found marked geographical differences, with percent excess of up to 88.9% (81.9% -95.2%) at the peak of the pandemic, in the city of Bergamo (Lombardia).

We provide policy makers with a flexible modelling framework for estimating trends of excess mortality at subnational level. This can be used to highlight geographical differences and temporal lags. The model will be used as a real-time surveillance tool, to track mortality rates at area level informing about the end of the first wave, or the start of a second one.

MBl was partially covered by National Institute of Health research grant, number R01HD092580-01A1.

We searched PubMed, medRxiv, bioRxiv, arXiv, and Wellcome Open Research using the terms COVID-19, coronavirus and mortality or excess mortality and similar terms up to 31st of May. We found only two studies which estimated excess mortality at regional level, in England and Wales and in Italy, while the other studies analysed mortality at national level or focused only on specific cities. In order to understand the geographical differences of mortality during the COVID-19 pandemic, there is the urgent need to study mortality at subnational level.

We present the first spatio-temporal analysis of all-cause excess mortality across Italy, which is one of countries most dramatically hit by COVID-19 pandemic. We modelled the weekly counts of deaths for all-causes at municipality level over the first four months of 2016 -2019, then we predicted the mortality for the corresponding months in 2020 accounting for spatial and temporal variations and according to the historical trends estimated by the model (i.e. as if the COVID-19 pandemic did not occur). The comparison between the observed and predicted deaths in 2020 allowed us to estimate the weekly excess of mortality at municipal level. We identified considerable geographical variation in the local trends of the excess mortality. Our analysis also makes use of fully probabilistic and fully reproducible methods, which are able to flexibly account for latent spatial factors, long-and short-term trends and for local fluctuations in the temperature-mortality relationship. These can be used to timely monitor mortality and can be easily extended to analyse routinely collected mortality data in other countries.

The impact of the COVID-19 outbreak on population health is made by the direct and the indirect effects of the disease, as a result of the social and economic measures implemented to contain it and the resulting pressure on health-care provision.

While mortality should be the least controversial outcome to measure, its analysis is in fact complicated by the lack of real time cause specific data; a potential additional issue concerns the quality of coding on the death certificates, particularly in the earlier stages of the pandemic. In addition, there are important differences in the recording systems, both across and within countries (e.g. in the UK, England and Wales have consistently reported data on daily deaths based on different time of recording 1 ).

In this context, estimating excess deaths for all causes with respect to past year trends has been used in several countries as an effective way to evaluate the total burden of the COVID-19 pandemic 2-5 , including direct COVID-19-related as well as indirect effects (e.g. people not being able to access health-care). At the same time, looking at all-cause mortality is not affected by mis-coding on the death certificates 6 .

The recent contributions to estimate excess mortality have focused on nationwide data 7 and have proved accurate, if only presenting global pictures of the total burden of the first wave of the infection. Nevertheless, in order to understand the dynamics of the pandemic, there is the need to analyse data at subnational level, accounting for geographical differences due to the infectious nature of the disease, as well as the differences in the population characteristics and health system provision. Additionally, time trends can vary substantially, rendering comparisons even more complex. Recently, important contributions have emerged calling for health authorities of each country to release mortality data at high spatial and temporal resolution in order to evaluate the pattern of the disease and to compare it with the expected trend from previous years 8 . Up to date, to the best of our knowledge, only two papers have analysed mortality at regional level in England and Wales 9 and in Italy 10 , while there are no comprehensive studies that have looked at the impact of COVID-19 at subnational level.

In this paper, we present the first comprehensive study of excess mortality at subnational level in Italy, one of the countries with the largest number of deaths for COVID-19 (as of 4 June 2020, 33,601 confirmed fatalities, making it the third most affected country in the world, in absolute terms); we consider municipalities as geographical unit and we model weekly number of deaths for all causes. A report by the Italian National Institute of Health 11 showed the characteristics of the epidemics, both in terms of cases and mortality, across Italian regions. Conversely, our objective is to provide the community with a scalable method to obtain high resolution temporal trends of excess deaths across municipalities, in order to highlight similarities and differences in space and time.

We use 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 model weekly mortality trends over the years 2016 -2019, while accounting for air temperature. We use these trends to predict the weeklyand municipality-based mortality over the period 1 January to 29 April 2020. Then, we compare the observed mortality for this period with the model-based predicted mortality. 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 calculate 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 apply the age-specific rates to the year-and age-specific population in each municipality to estimate the expected number of cases as denominator.

We used official weekly mortality data, available from the Italian Institute of Statistics (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 sensors 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.

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 models have been run in 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.

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.

The funder of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report. All the authors had full access to all the data in the study and the corresponding author had final responsibility for the decision to submit for publication.

For males, the total number of deaths in the first four months of 2020 was 122,129, progressing from 6,278 in the week commencing on 1 January to 5,748 in the week commencing on the 22 April; for females the total number of deaths was 129,044, going from 6,633 to 6,614 during the same observation period. Nationally, the highest weekly death toll was of 11,301 in the week of 18 March and of 10,612 in the week of 25 March for males and females, respectively.

We present the trend in mortality rates split by five macroareas: the North West (Piemonte, Valle d'Aosta, Liguria); Lombardia; the North East (Veneto, Trentino-Alto Adige, Friuli Venezia Giulia, Emilia-Romagna); the Centre (Lazio, Marche, Toscana and Umbria); and the South and major Islands (Abruzzo, Basilicata, Calabria, Campania, Molise, Puglia, Sardegna and Sicilia). These are consistent with the Eurostat classification, although we report Lombardia (which would normally be included in the North West) separately, as it was the most affected region (Figures 1 and 2) . For both sexes, the model based on 2016 -2019 predicts a slightly decreasing trend for the entire period, with the blue curve and the gray ribbon indicating respectively the mean and 95% interval estimate for the posterior predictive distribution of the number of deaths.

Lombardia is the first region to show a sharp increase from the estimated trend, starting with males in in the week of the 26 of February (observed rate = 22.6 per 100,000 vs upper limit of the 95% posterior interval = 22.1 per 100,000). Mortality rates for males in North West and the North East regions start to deviate substantially from the estimated trend a week after Lombardia (starting 4 March), while the Centre shows a small deviation from the week starting with 11 March and the South does not show substantial differences from the model estimates of the time trend. For females there is a lagged effect of the pandemic, with mortality rates deviating from the expected trend on 4 March in Lombardia and on the 11 March in North East and North West. The Centre and South do not show any evidence of excess death. It is worth noting that at the end of our follow up, while sharply falling down, mortality rates are at the level of (or above) those normally observed at the beginning of January.

The extent of geographical differences is clear if looking at the trend in the percent excess mortality by province (Figures 3 and 4 for males and females, respectively). In the maps, darker shades of colour indicate larger percent excess mortality in a given week (we only present the progression from the week commencing on 26 February and for the provinces in the North West, Lombardia and North East; maps for all the 107 provinces are presented in the Supplementary material). There is a clear increasing darkening of the hot-spots at the heart of Lombardia (starting with Lodi and moving to the provinces of Cremona, Bergamo and Brescia), spreading mostly to the West into Piemonte, South-West into parts of Emilia-Romagna and mainly North-East into Trentino-Alto Adige. Interestingly, despite recording the first COVID-19 death in Italy, the region of Veneto (in the North East) is affected to a much lower extent, in comparison to neighbouring areas of Lombardia.

Looking at specific municipalities the full extent of the geographical variation is even more striking: for males ( Figure 5 ) the city of Piacenza is the first to reach the peak of the pandemic on the week starting on 11 March, when the observed mortality rates reach 177 per 100,000 residents, corresponding to an excess 5 mortality of 86.9%. Bergamo and Pesaro (in the region of Marche, in Central Italy) peak one week later, with mortality rates of 182 and 134 per 100,000, corresponding to an excess mortality of 88.9% and 84.2%. Brescia shows a lagged peak, occurring over the weeks of 18 and 25 March and slightly lower mortality rates of 100 per 100,000, with a corresponding excess of 80.6%. Meanwhile, Milano, (the second largest city in Italy) is characterised by a more diluted peak and much lower rates (51 per 100,000) with an excess of 63.2%. Finally, Verona (located in Veneto and bordering the Lombard province of Brescia) shows rates barely above the expected interval in the peak weeks. This effect is possibly associated with the comprehensive measures and testing system implemented by the Veneto region to contain the epidemic 21 . Females show slightly lower rates and a seem to peak more or less at the same time, on the week of 18 March. The only exceptions are possibly represented by Milano, which still presents a diluted peak and Verona, which shows no clear evidence of a peak (Figure 6 ). The plot of the weekly excess mortality for the same municipalities is presented in the Supplementary material.

We combined the output from our analysis with official data on confirmed deaths for COVID-19 in Italy in the March-April period. Table 1 shows that if we crudely remove the confirmed COVID-19 deaths given by the estimates from the posterior predictive distributions of the total excess deaths obtained from our model, Lombardia is still left, on average, with a staggering 10,197 (9,264 -11,037) excess deaths. These may be attributable to the fact that some individuals who would normally access the health care provision were somehow prevented or delayed in doing so during the past few months, resulting in deaths that would have occurred at a lower pace in normal circumstances. Similarly, but to a lower extent, the North West and North East show indirect effect of the pandemic on mortality, with an excess, after roughly discounting the confirmed COVID-19 fatality of 2,572 (1,772 -3,297) and 2,047 (1,075 -3,058) respectively. The Centre shows results that span across 0 in the 95% posterior interval (indicating weak evidence of any difference in comparison to the estimated trends), while the South of Italy shows a negative excess after accounting for COVID-19 (with on average 1,056 deaths lower than estimated from the model).

In this paper we presented the first fully probabilistic analysis of the excess mortality at weekly subnational level in Italy, one of countries with the highest number of deaths during the COVID-19 epidemic. We used data from 1 January to 28 April 2016 -2019 to estimate the mortality rates for each sex, adjusting for the population age structure. We then predicted the mortality rates for 2020 under the estimated trend from past years, under the alternative scenario that the pandemic did not take place. Substantively, 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 the previous years, albeit we see strong geographical differences. While an early comparative study on COVID-19 related mortality suggested that the region Lombardia was not heavily impacted by the epidemic wave 22 , in accordance with more recent studies 23, 24 we found that Lombardia is the most affected region, for both sexes; nevertheless, we are able to detect marked differences at municipality level.

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 the mortality. This is consistent with typical modelling in disease mapping studies 25 and surveillance studies 26 . Additionally, we are framed in the same perspective as forecasting studies 27,28 and we take advantage of 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 which are only varying in space as 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.

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 returns to mortality rates in line with the previous years and the South even is characterised by lower death rates than estimated, 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, mortality rates return to or above levels normally observed at the beginning of January, by the end April. 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, the aggressive programme of mass-testing implemented in Veneto is potentially a driving factor in the much lower excess mortality experience by the region, in comparison to neighbouring areas such as Lombardia and, albeit to a lower extent, Trentino-Alto Adige.

It emerged that deaths attributed to COVID-19 may be significantly undercounted 29 . To quantify the under-reporting for Italy, we compared the estimated excess mortality with the official mortality statistics reported by ISTAT, which distinguish COVID-19 registered deaths from the others. 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, especially, the South of Italy show some signs of reduced net mortality, possibly due to the effects of milder flu season.

The model presented 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 are 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 account for less than 1% of the total deaths for the considered period and our results show that the excess mortality more than compensate that. Additionally, there have been several reports which have showed improvement in air quality during the lockdown in different countries 30,31 , leading to a consistent reduction of the pollution-related mortality 32, 33 . 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.

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 makes our modelling framework particularly valuable in the current pandemic, as it would allow, as new data come in, 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 perspective real-time surveillance, as allows to effectively follow each municipality in time and flag unusual behaviour as soon as they happen, aiming at identifying potential new public health threats and contain them.

MBl, MP, MC and GB contributed equally to the design, literature search, statistical analyses, data interpretation and writing. GC and MBa were responsible for data collection and validation. All the authors have commented and approved the final manuscript. 

Why counting coronavirus deaths is not an exact science. The Guardian

Excess mortality from COVID-19. Weeekly excess death rates by age and sex for Sweden

Baseline characteristics and outcomes of 1591 patients infected with SARS-CoV-2 admitted to ICUs of the Lombardy Region

Mortality impacts of the coronavirus disease

outbreak by sex and age: rapid mortality surveillance system

Disease mapping. Encyclopedia of

ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate

Bayesian disease mapping for public health

A comparison of Bayesian spatial models for disease mapping

Bayesian image restoration, with two applications in spatial statistics

Projections of temperature-related excess mortality under climate change scenarios. The Lancet Planetary Health

Mortality risk attributable to high and low ambient temperature: a multicountry observational study

In one Italian town, we showed mass testing could eradicate the coronavirus. The Guardian

The spread of COVID-19 in six western metropolitan regions: a false myth on the excess of mortality in Lombardy and the defense of the city of Milan

COVID-19 deaths in Lombardy, Italy: data in context. The Lancet Public Health

COVID-19 and Italy: what next? The Lancet 2020

A space-time multivariate Bayesian model to analyse road traffic accidents by severity

A Bayesian mixture modeling approach for public health surveillance

The future of life expectancy and life expectancy inequalities in England and Wales: Bayesian spatiotemporal forecasting

Small area forecasts of cause-specific mortality: application of a Bayesian hierarchical model to US vital registration data

Assessment of deaths from COVID-19 and from seasonal influenza

Impact of coronavirus outbreak on NO 2 pollution assessed using TROPOMI and OMI observations

COVID-19 pandemic and environmental pollution: A blessing in disguise?

000 air pollution-related deaths avoided in Europe as coal, oil consumption plummet

Air pollution reduction and mortality benefit during the COVID-19 outbreak in China. The Lancet Planetary Health

Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Mortality rate for all causes (x100,000)

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 

Predicted mortality rate x100,000 (mean and 95% interval) Observed mortality rate x100,000 Week Deaths for all causes (x100,000)