key: cord-0209370-wnimqeek
authors: Konstantinoudis, Garyfallos; G'omez-Rubio, Virgilio; Cameletti, Michela; Pirani, Monica; Baio, Gianluca; Blangiardo, Marta
title: A framework for estimating and visualising excess mortality during the COVID-19 pandemic
date: 2022-01-17
journal: nan
DOI: nan
sha: e0ad2c820c7bcb732890de4d46d6a204140ee429
doc_id: 209370
cord_uid: wnimqeek

COVID-19 related deaths underestimate the pandemic burden on mortality because they suffer from completeness and accuracy issues. Excess mortality is a popular alternative, as it compares observed with expected deaths based on the assumption that the pandemic did not occur. Expected deaths had the pandemic not occurred depend on population trends, temperature, and spatio-temporal patterns. In addition to this, high geographical resolution is required to examine within country trends and the effectiveness of the different public health policies. In this tutorial, we propose a framework using R to estimate and visualise excess mortality at high geographical resolution. We show a case study estimating excess deaths during 2020 in Italy. The proposed framework is fast to implement and allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. This makes it particularly powerful and appealing for online monitoring of the pandemic burden and timely policy making.

Estimating the burden of the COVID-19 pandemic on mortality is an important challenge (Weinberger et al., 2020) . COVID-19 related deaths are subject to testing capacity and changes in definition and reporting, raising accuracy and completeness considerations (Aburto et al., 2021; Konstantinoudis et al., 2021) . In addition, COVID-19 related deaths give an incomplete picture of the burden of the COVID-19 pandemic on mortality, as they do no account for the indirect pandemic effects due to, for instance, disruption to health services (Kaczorowski and Del Grande, 2021) . Alternatively, excess mortality has been extensively used to evaluate the impact of the COVID-19 pandemic on mortality (Rossen et al., 2020; Islam et al., 2021; Kontis et al., 2020; Konstantinoudis et al., 2021; Verbeeck et al., 2021) .

Excess mortality is estimated by comparing the observed number of deaths in a particular time period with the expected number of deaths under the counterfactual scenario that the pandemic did not occur. Typically this counterfactual scenario is calculated using a comparison period, for instance several previous years (https://www.euromomo.eu/). Calculating accurately the expected deaths requires accounting for factors such as population trends, seasonality, temperature, public holidays and spatio-temporal dependencies. While accounting for the above-mentioned factors, most studies to date have estimated excess mortality at national level (e.g., Rossen et al., 2020; Weinberger et al., 2020) and a few have looked across different countries reporting important geographical discrepancies (Islam et al., 2021; Kontis et al., 2020 Kontis et al., , 2021 .

While national level estimates of excess mortality give valuable insights into the total number of excess deaths, they do not allow the evaluation of within country geographical disparities. Such information is essential to understand the country's transmission patterns and effectiveness of local policies and measures to contain the pandemic (Kontopantelis et al., 2021) . Temporal trends in excess can substantially differ across regions of the same country, which makes national-based comparisons even more challenging (Blangiardo et al., 2020) . Therefore, understanding the burden of the COVID-19 pandemic on mortality requires higher spatial resolution and models that account for spatial, temporal and spatio-temporal dependencies.

When working at high spatio-temporal resolution, data are generally sparse, leading to excess mortality estimates that are highly variable. This is aggravated by the fact that they are subject to spatial and temporal correlation, making it essential to use statistical methods that account for these dependencies in order to provide robust and accurate estimates. The disease mapping framework, which is commonly employed in epidemiology to study the spatio-temporal variation of diseases (Waller et al., 1997; Moraga, 2018) , can be exploited to estimate excess mortality at subnational and weekly level. The Bayesian approach is naturally suited in this context, as it is able to model complex dependency structures, as well as to incorporate uncertainty in the data and modelling process. In addition, while fully propagating the uncertainty, it allows us to summarise the results at any desired spatio-temporal aggregation (using for instance coarser geographical units suitable for policy implementation). This in combination with fast approximate methods to inference, such as the Integrated Laplace Approximation (INLA, Rue et al. (2009)) , make this framework particularly powerful and appealing for monitoring of the pandemic burden and timely policy making.

In this tutorial, we describe how to run a study on excess mortality at high spatial and temporal resolution using a Bayesian approach and R. This tutorial provides a detailed implementation of the approach followed previously in 5 European regions (Konstantinoudis et al., 2021) . Figure 1 illustrates graphically the workflow followed in this paper together with the main r-packages used. Source code for replicating the data wrangling (r-files starting with 01, 02 or 03), analysis (r-files starting with 04) and post-processing steps (r-files starting with 05 or 06 and the App folder) and data files are available from GitHub at https://github.com/gkonstantinoudis/TutorialExcess. This tutorial is structured as follows: we first describe the modelling framework and present the case study in Italy. We then show how to run and evaluate the model, and extract and visualise the results. Finally, we present an R-shiny app which makes the results effectively and easily disseminated. 

We propose a Bayesian hierarchical model to quantify spatio-temporal excess mortality under extreme events such as the COVID-19 pandemic, stratified by specific age-sex population groups. To do so, we first describe the statistical model for predicting the number of deaths from all-causes based on historical data, in the counterfactual scenario in which the pandemic did not take place. Then, we show how to estimate the magnitude of excess deaths over a specific period of time, with associated uncertainty, by comparing the predicted versus the actual number of deaths.

Let y jtsk and P jtsk be the number of all-cause deaths and the population at risk for the j-th time point (e.g., week or month) of the t-th year preceding the pandemic, the s-th spatial unit, and k-th age-sex group. Also let x j and z jts denote the adjustment covariates, respectively public holidays (i.e., x j = 1 if week j contains a public holiday and 0 otherwise) and temperature. We assume that the historical observed number of deaths, conditional on the risk r jtsk , follows a Poisson distribution, with a log-linear model on the risk: y jtsk |r jtsk ∼Poisson(µ jtsk = r jtsk P jtsk ),

Here, β 0t is the year specific intercept given by β 0t = β 0 + t , where β 0 is the global intercept and t is an unstructured random effect representing the deviation of each year from the global intercept, which is modelled using independent and identically (iid) distributed Gaussian prior distribution with zero-mean and variance equal to τ −1 . The parameter β is an unknown regression coefficient associated to the public holidays covariate x j . The effect of temperature is allowed to be non-linear f (·) by specifying a second-order random walk (RW2) model:

where τ −1 z is the variance. Terms b s and w j are spatial and temporal random effects, respectively. We specify the spatial random effect term, b s , with a convolution prior (Besag et al., 1991) , and the temporal random effect term, w j with a non-stationary in time prior. In detail, we model b using a reparameterisation of the popular Besag-York-Mollié prior, which is a convolution of an intrinsic conditional autoregressive (CAR) model and an iid Gaussian model. Let u s be the spatially structured component defined by an intrinsic CAR (Besag, 1974) 

whereū is the local mean and ∂ s and m s are respectively the set and the number of neighbours of area s, τ −1 u the conditional variance, and v s the unstructured component with prior v s iid ∼ Normal(0, τ −1 v ). We model b s as follows (Besag et al., 1991; Riebler et al., 2016; Konstantinoudis et al., 2020) :

where u s and v s are standardised versions of u s and v s to have variance equal to 1 (Simpson et al., 2017) . The term 0 ≤ φ ≤ 1 is a mixing parameter, which measures the proportion of the marginal variance explained by the structured effect. Finally, we assign to the temporal random effect term, w j , a Gaussian random walk model of order 1 (RW1), specified as:

The Bayesian representation of the above model is completed once we select priors for the fixed effects β 0 and β and the hyperparameters: τ , τ b , τ w , and φ. For the fixed effects we selected minimally informative Normal distributions whereas, we specified "penalising complexity" (PC) priors (Simpson et al., 2017) for the hyperparameters. PC priors are defined by penalising deviations from a "base" model (e.g., specified in terms of a specific value) and have the effect of regularising inference, while not implying too strong information. Technically, PC priors imply an exponential distribution on a function of the Kullback-Leibler divergence between the base model and an alternative model in which the relevant parameter is unrestricted. This translates to a suitable "minimally informative", regularising prior on the natural scale of the parameter.

In order to quantify the weekly excess mortality at sub-national level for specific age-sex population groups, we need to predict the number of deaths had COVID-19 not occurred. In Bayesian analysis, this can be performed by drawing random samples from the posterior predictive distribution (that is, the distribution of unobserved values conditional on the observed values from previous years). Specifically, letting θ θ θ be the model parameters, D be the observed data, and y jt * sk be the count of deaths that we want to predict, we have:

Operationally, we first generate random samples from the joint posterior marginal of the fitted linear predictor specified in equations (1) at the highest spatial resolution available (NUTS3 regions; Nomenclature of Territorial Units for Statistics 3 regions, https://ec.europa.eu/eurostat/web/nuts/ background/). Successively, we use these to be in turn the mean parameter of a Poisson distribution, to obtain the predicted number of deaths, which represents the baseline number of deaths assuming the pandemic did not take place.

Finally, to estimate the magnitude of excess deaths, the predicted number of deaths is compared against the actual observed number of deaths, allowing the computation of the relative change in the mortality (i.e., relative to what we could expect if the pandemic did not occur). This is obtained by (i) subtracting the predicted number of deaths from the observed number of deaths in each time point t and spatial unit s (number of excess deaths or NED), and (ii) dividing NED by the predicted number of deaths for each sample and multiplying by 100 (% relative excess mortality or REM).

Bayesian inference for the model is computed using Integrated Nested Laplace Approximation (INLA; Rue et al. (2009) , which performs approximate Bayesian inference on the class of latent Gaussian models (Rue and Held, 2005) . Unlike simulation based Markov chain Monte Carlo method, INLA is a deterministic algorithm, which employs analytical approximations and efficient numerical integration schemes to provide accurate approximations of the posterior distributions in short computing times. The INLA software is provided through the R package INLA, which can be downloaded from https: //www.r-inla.org/download-install. 

We retrieved all-cause mortality data during 2015-2020 in Italy from the Italian National Institute of Statistics (https://www.istat.it/). Data were available weekly (ISO week), by age (5-year age groups), sex and NUTS3 regions. As the COVID-19 mortality rates increase with age, we aggregated mortality counts based on the following age groups: <40, 40-59, 60-69, 70-79 and 80 years and older (Davies et al., 2021) .

Population data in Italy during 2015-2020 were retrieved from the Italian National Institute of Statistics. The data represent the population in Italy on January 1st of every year and are available by age (5-year age groups), sex and NUTS3 regions. To retrieve weekly population, we performed linear interpolation by the selected age groups (<40, 40-59, 60-69, 70-79 and 80+), sex and NUTS3 regions using populations at January 1st of the current and the next year. Population counts on January 1st 2021, which takes COVID-19 deaths in 2020 into consideration, were available at the time of analysis. Our goal was, however, to predict mortality for 2020 had the pandemic not occurred. Thus we performed an additional linear interpolation by age, sex and NUTS3 regions to predict the population at January 1st 2021, using the years 2015-2020 (Figure 2, panel A) . Object pop is a tibble containing the NUTS3 region ID (NUTS318CD), age group (ageg), sex (sex), year (year) and population counts ( We can use the following code (based on tidyverse and "piping" principles) to calculate the population of the 1st of January 2021 by NUTS3 regions, sex and age: pop %>% group_by(NUTS3, sex, age) %>% summarise(pop = as.vector(coef(lm(pop~year)) %*% c(1, 2021))) %>% mutate(year = 2021) -> pop2021

Once we obtained the year 2021 we performed an additional linear interpolation to calculate weekly number of population as shown on Figure 2 , panel B.

We used covariates related with ambient temperature and national holidays to help the model predictions. Data on air-temperature during 2015-2020 in Italy at 2m above the surface of land were retrieved from the ERA5 reanalysis data set of the Copernicus climate change program (Hersbach et al., 2020) . The geographical resolution of the ERA5 estimates is 0.25 • × 0.25 • (panel A of Figure 3 ). We calculated the weekly mean by the centroids of the 0.25 • × 0.25 • grid (panel B of Figure 3 ) and then averaged the weekly temperature over the ERA5 centroids that overlay with the NUTS3 regions (panels B and C of Figure 3 ).

The modelling process in INLA consists of three main steps: (1) the selection of priors; (2) definition of the model "formula" (which sets out the expression for the generalised linear predictor); and (3) the call to the main function inla, which computes the estimates.

In particular, we constructed the PC priors for σ = √ 1/τ ,

1/τ w based on statement that it is unlikely to have a relative risks higher than exp(2) based solely on spatial, yearly and seasonal variation, Figure 4 , panel A. For the mixing parameter φ, we set Pr(φ < 0.5) = 0.5 reflecting our lack of knowledge about whether overdispersion or strong spatial autocorrelation should dominate the field b, Figure 4, Figure 4 , panel B.

These assumptions can be encoded using the following code: # Defines the priors hyper.bym <-list( theta1 = list( PCprior , c(1, 0.01)), theta2 = list( PCprior , c(0.5, 0.5)) ) hyper.iid <-list(theta = list(prior="pc.prec", param=c(1, 0.01))) # Defines the model "formula" After fitting the model, we take 1000 samples from the (approximated) posterior distribution of the linear predictor and we use each drawn sample as the mean of a Poisson distribution to retrieve the predicted mortality counts: post.samples <-inla.posterior.sample(n = 1000, result = inla.mod) predlist <-do.call(cbind, lapply(post.samples, function(X) exp(X$latent[startsWith(rownames(X$latent), "Pred")]))) pois.samples <-apply(predlist, 2, function(Z) rpois(n = length(Z), lambda = Z)) This allows us to estimate the entire predictive posterior distribution of the mortality counts incorporating both the sampling and the linear predictor uncertainty.

To examine model validity, we performed a cross validation leaving out one historical year at a time and predicting the weekly number of deaths by NUTS3 regions for the year left out. For each stratum, we calculated the correlation between observed and fitted and a coverage probability, i.e. the probability that the observed death fall into the 95% credible interval (95% CrI) of the predicted.

We overall found good predictive ability of our model. 

The object pois.samples.list contains 1000 samples from the posterior predictive distribution 2, i.e. 1000 samples of the expected number of deaths by age, sex, NUTS3 regions and week, had the pandemic not occurred. We can access the different age-sex groups as follows:

names(pois.samples.list) # "F_less40" "F_40_59" "F_60_69" "F_70_79" "F_80plus" # "M_less40" "M_40_59" "M_60_69" "M_70_79" "M_80plus"

where F stands for females and M for males across the different age groups. To get an idea about the structure of the data, we can use the head ( We can also calculate median and 95% CrI expected number of deaths for this specific age-sex group over the year by NUTS3 region:

pois.samples.list$F_60_69 %>% select(starts_with("V"), "ID_space") %>% group_by(ID_space) %>% summarise_all(sum) %>% rowwise(ID_space) %>% mutate(median = median(c_across(V1:V1000)), LL = quantile(c_across(V1:V1000), probs= 0.025), UL = quantile(c_across(V1:V1000), probs= 0.975)) %>% select ( 

The above results can be combined in different ways using the functions get2020data() and get2020weeklydata() to calculate excess mortality. The function get2020data() aggregates over the entire country, NUTS2 regions, sex, age and time resulting in the object d, whereas the function get2020weeklydata() over the entire country, NUTS2 regions, sex and age but not over time resulting in the object d_week. names(d) # "province" "region" "country" names(d$province) # "none" "age" "sex" "agesex" names(d$province$age) # "40<" "40-59" "60-69" "70-79" "80+" names(d$province$sex) # "F" "M" names(d$province$agesex) # "F40<" "F40-59" "F60-69" "F70-79" "F80+" "M40<" "M40-59" "M60-69" "M70-79" "M80+" Notice that the object d$province$none is a simple feature collection, making mapping it straightforward.

In Figure 5 (plots 1A, 2A and 3A) we show the median posterior of REM for total age and sex at the national, NUTS2 and NUTS3 regional level. For these plots we used the object d with the selection "none" and plot the median REM (median.REM), for example for 3A: prov <-"Foggia" d$province$none %>% filter(DEN_UTS == prov) %>% select(geometry) %>% ggplot() + geom_sf(data = d$province$none, aes(fill = median.REM.cat)) + geom_sf(fill = NA, col = col.highlight, size = .8) + scale_fill_manual(values=colors, name = "", drop=FALSE) + theme_light() + ggtitle(paste0("3A. NUTS3 regions: ", prov)) Overall, the REM in Italy during 2020 was between 5-10%, Figure 5 , panel 1A. When the higher geographical resolution is assessed, it is revealed that north and in particular Lombardia is the region hit the worst, with the REM exceeding 20%, Figure 5 , panels 1B and 1C. Figure 6 shows a measure of uncertainty of the REM, now for the different age groups and both sexes (selection "age"). The probability of a positive excess (exceedance.REM) in older people and in the north of Italy is larger than 0.80, Figure 6 . Panels 1B and 1C of Figure 5 show the median temporal nationwide excess together with 95% CrI by sex after aggregating the different age groups (using d_week and the "sex" selection). We observe a clear first pandemic wave during March and May and a second one during mid October and December in 2020. During the first pandemic wave, there were weeks when the median relative excess death reached almost 100% in males, Figure 5 . Panels 2B, 2C, 3B and 3C of Figure 5 show the median spatio-temporal excess together with 95% CrI by sex after aggregating over the different age groups. Panels 2B and 2C highlight the region of Puglia, where during the first wave of the pandemic experienced a positive, similar with the nationwide excess. When we increase the spatial resolution in panels 3B and 3C, we highlight the region of Foggia, where there was insufficient evidence of a positive excess during the first wave, but strong during the second.

To be able to effectively examine and communicate the different aggregation levels of the output of our modelling framework, we have also developed a Shiny Web-Application (WebApp), Figure 7 . The WebApp provides spatial, temporal and spatio-temporal analysis tabs, and within each tab there are plots and summary statistics for the level of aggregation selected from the drop-down menu. Users can select across different variables (relative excess mortality or number of excess deaths), statistics (median or posterior probability), sex (males, females or both), age group (40<, 40-59, 60-69, 70-79, 80> and all) and different geographical level (national, NUTS2 or NUTS3 regions). Summary statistics for each area are available and they are displayed in a pop-up window, which is activated by clicking on the area of interest. In addition, graphical pop-ups are provided to show the weekly estimates for each area with the leafpop R-package (Appelhans and Detsch, 2021) , in the spatio-temporal analysis tab. The WebApp that we have developed is hosted at http://atlasmortalidad.uclm.es/italyexcess.

This tutorial provides a detailed implementation of the framework followed previously to calculate excess mortality during the COVID-19 pandemic in 5 European regions (Konstantinoudis et al., 2021) . We have proposed a Bayesian framework for estimating excess mortality and shown how to use R and INLA to retrieve fast and accurate estimates of the excess mortality. The proposed framework also allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. We have given a practical example of how to use the proposed framework modelling the excess mortality during the 2020 COVID-19 pandemic in Italy at small area level. We also developed a Webapp to effectively communicate the results. This framework can be extended to monitor mortality from other extreme events for instance natural hazards such as hurricanes (Acosta and Irizarry, 2020) . All the above make the proposed framework particularly powerful, generalisable and appealing for online monitoring of the pandemic burden and timely policy making. Potential extensions include different ways of modelling the younger age groups to increase the predictive ability, for instance using a zero-inflated Poisson. 

V.G.R. conceived the study. M.B. supervised the study. G.K. developed the initial study protocol and discussed it with M.B., M.C., M.P. and G.B.. G.K. developed the statistical model, prepared the population and covariate data and led the acquisition of mortality data. M.C. and V.G.R. validated and modified accordingly the code. G.K. ran the analysis. G.K., V.G.R., M.B., M.C. and M.P wrote the initial draft and all the authors contributed in modifying the paper and critically interpreting the results. V.G.R. developed the Shiny app. All authors read and approved the final version for publication.

A flexible statistical framework for estimating excess mortality. medRxiv

leafpop: Include Tables, Images and Graphs in Leaflet Pop-Ups

Spatial interactions and the statistical analysis of lattice systems (with discussion)

Bayesian image restoration, with two applications in spatial statistics

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

Community factors and excess mortality in first wave of the COVID-19 pandemic in England

The ERA5 global reanalysis

Excess deaths associated with COVID-19 pandemic in 2020: age and sex disaggregated time series analysis in 29 high income countries

Beyond the tip of the iceberg: direct and indirect effects of COVID-19

Discrete versus continuous domain models for disease mapping

Regional excess mortality during the 2020 COVID-19 pandemic: a study of five european countries. medRxiv

Magnitude, demographics and dynamics of the effect of the first wave of the COVID-19 pandemic on all-cause mortality in 21 industrialized countries

Lessons learned and lessons missed: impact of the coronavirus disease 2019 (covid-19) pandemic on all-cause mortality in 40 industrialised countries prior to mass vaccination

Excess mortality in england and wales during the first wave of the covid-19 pandemic

Small Area Disease Risk Estimation and Visualization Using R

An intuitive Bayesian spatial model for disease mapping that accounts for scaling

Excess deaths associated with COVID-19, by age and race and ethnicity-United States

Gaussian Markov random fields: theory and applications

Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations

Penalising model component complexity: A principled, practical approach to constructing priors

A linear mixed model to estimate covid-19-induced excess mortality

Hierarchical spatio-temporal mapping of disease rates

Estimation of excess deaths associated with the COVID-19 pandemic in the United States

All authors acknowledge infrastructure support for the Department of Epidemiology and Biostatistics provided by the NIHR Imperial Biomedical Research Centre (BRC). The authors also acknowledge infrastructure support for the domain of the shiny app from the University of Castilla-La Mancha. We thank Univesidad de Castilla-La Mancha for hosting the server on which the Shiny App is running.

Province stands for the NUTS3 regions (the resolution we used to fit the models), region for the NUTS2 (coarser than NUTS3, appropriate for policy making) and country for the nationwide aggregation. Within these aggregations users can select the option "none" being the total aggregation by age and sex, "age" by sex, "sex" by age and "agesex" refers to no age and sex aggregation. The objects d and d_week have similar structure and contain summary statistics for REM and NED and posterior probabilities of a positive REM or NED: 

The authors declare no competing interests.

The study is about secondary, aggregate anonymised data so no additional ethical permission is required. 
key: cord-0829925-zd59ip30
authors: Konstantinoudis, G.; Gómez-Rubio, V.; Cameletti, M.; Pirani, M.; Baio, G.; Blangiardo, M.
title: A framework for estimating and visualising excess mortality during the COVID-19 pandemic
date: 2022-01-17
journal: ArXiv
DOI: nan
sha: e0ad2c820c7bcb732890de4d46d6a204140ee429
doc_id: 829925
cord_uid: zd59ip30

COVID-19 related deaths underestimate the pandemic burden on mortality because they suffer from completeness and accuracy issues. Excess mortality is a popular alternative, as it compares observed with expected deaths based on the assumption that the pandemic did not occur. Expected deaths had the pandemic not occurred depend on population trends, temperature, and spatio-temporal patterns. In addition to this, high geographical resolution is required to examine within country trends and the effectiveness of the different public health policies. In this tutorial, we propose a framework using R to estimate and visualise excess mortality at high geographical resolution. We show a case study estimating excess deaths during 2020 in Italy. The proposed framework is fast to implement and allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. This makes it particularly powerful and appealing for online monitoring of the pandemic burden and timely policy making.

Estimating the burden of the COVID-19 pandemic on mortality is an important challenge (Weinberger et al., 2020) . COVID-19 related deaths are subject to testing capacity and changes in definition and reporting, raising accuracy and completeness considerations (Aburto et al., 2021; Konstantinoudis et al., 2021) . In addition, COVID-19 related deaths give an incomplete picture of the burden of the COVID-19 pandemic on mortality, as they do no account for the indirect pandemic effects due to, for instance, disruption to health services (Kaczorowski and Del Grande, 2021) . Alternatively, excess mortality has been extensively used to evaluate the impact of the COVID-19 pandemic on mortality (Rossen et al., 2020; Islam et al., 2021; Kontis et al., 2020; Konstantinoudis et al., 2021; Verbeeck et al., 2021) .

Excess mortality is estimated by comparing the observed number of deaths in a particular time period with the expected number of deaths under the counterfactual scenario that the pandemic did not occur. Typically this counterfactual scenario is calculated using a comparison period, for instance several previous years (https://www.euromomo.eu/). Calculating accurately the expected deaths requires accounting for factors such as population trends, seasonality, temperature, public holidays and spatio-temporal dependencies. While accounting for the above-mentioned factors, most studies to date have estimated excess mortality at national level (e.g., Rossen et al., 2020; Weinberger et al., 2020) and a few have looked across different countries reporting important geographical discrepancies (Islam et al., 2021; Kontis et al., 2020 Kontis et al., , 2021 .

While national level estimates of excess mortality give valuable insights into the total number of excess deaths, they do not allow the evaluation of within country geographical disparities. Such information is essential to understand the country's transmission patterns and effectiveness of local policies and measures to contain the pandemic (Kontopantelis et al., 2021) . Temporal trends in excess can substantially differ across regions of the same country, which makes national-based comparisons even more challenging (Blangiardo et al., 2020) . Therefore, understanding the burden of the COVID-19 pandemic on mortality requires higher spatial resolution and models that account for spatial, temporal and spatio-temporal dependencies.

When working at high spatio-temporal resolution, data are generally sparse, leading to excess mortality estimates that are highly variable. This is aggravated by the fact that they are subject to spatial and temporal correlation, making it essential to use statistical methods that account for these dependencies in order to provide robust and accurate estimates. The disease mapping framework, which is commonly employed in epidemiology to study the spatio-temporal variation of diseases (Waller et al., 1997; Moraga, 2018) , can be exploited to estimate excess mortality at subnational and weekly level. The Bayesian approach is naturally suited in this context, as it is able to model complex dependency structures, as well as to incorporate uncertainty in the data and modelling process. In addition, while fully propagating the uncertainty, it allows us to summarise the results at any desired spatio-temporal aggregation (using for instance coarser geographical units suitable for policy implementation). This in combination with fast approximate methods to inference, such as the Integrated Laplace Approximation (INLA, Rue et al. (2009)) , make this framework particularly powerful and appealing for monitoring of the pandemic burden and timely policy making.

In this tutorial, we describe how to run a study on excess mortality at high spatial and temporal resolution using a Bayesian approach and R. This tutorial provides a detailed implementation of the approach followed previously in 5 European regions (Konstantinoudis et al., 2021) . Figure 1 illustrates graphically the workflow followed in this paper together with the main r-packages used. Source code for replicating the data wrangling (r-files starting with 01, 02 or 03), analysis (r-files starting with 04) and post-processing steps (r-files starting with 05 or 06 and the App folder) and data files are available from GitHub at https://github.com/gkonstantinoudis/TutorialExcess. This tutorial is structured as follows: we first describe the modelling framework and present the case study in Italy. We then show how to run and evaluate the model, and extract and visualise the results. Finally, we present an R-shiny app which makes the results effectively and easily disseminated. 

We propose a Bayesian hierarchical model to quantify spatio-temporal excess mortality under extreme events such as the COVID-19 pandemic, stratified by specific age-sex population groups. To do so, we first describe the statistical model for predicting the number of deaths from all-causes based on historical data, in the counterfactual scenario in which the pandemic did not take place. Then, we show how to estimate the magnitude of excess deaths over a specific period of time, with associated uncertainty, by comparing the predicted versus the actual number of deaths.

Let y jtsk and P jtsk be the number of all-cause deaths and the population at risk for the j-th time point (e.g., week or month) of the t-th year preceding the pandemic, the s-th spatial unit, and k-th age-sex group. Also let x j and z jts denote the adjustment covariates, respectively public holidays (i.e., x j = 1 if week j contains a public holiday and 0 otherwise) and temperature. We assume that the historical observed number of deaths, conditional on the risk r jtsk , follows a Poisson distribution, with a log-linear model on the risk: y jtsk |r jtsk ∼Poisson(µ jtsk = r jtsk P jtsk ),

Here, β 0t is the year specific intercept given by β 0t = β 0 + t , where β 0 is the global intercept and t is an unstructured random effect representing the deviation of each year from the global intercept, which is modelled using independent and identically (iid) distributed Gaussian prior distribution with zero-mean and variance equal to τ −1 . The parameter β is an unknown regression coefficient associated to the public holidays covariate x j . The effect of temperature is allowed to be non-linear f (·) by specifying a second-order random walk (RW2) model:

where τ −1 z is the variance. Terms b s and w j are spatial and temporal random effects, respectively. We specify the spatial random effect term, b s , with a convolution prior (Besag et al., 1991) , and the temporal random effect term, w j with a non-stationary in time prior. In detail, we model b using a reparameterisation of the popular Besag-York-Mollié prior, which is a convolution of an intrinsic conditional autoregressive (CAR) model and an iid Gaussian model. Let u s be the spatially structured component defined by an intrinsic CAR (Besag, 1974) 

whereū is the local mean and ∂ s and m s are respectively the set and the number of neighbours of area s, τ −1 u the conditional variance, and v s the unstructured component with prior v s iid ∼ Normal(0, τ −1 v ). We model b s as follows (Besag et al., 1991; Riebler et al., 2016; Konstantinoudis et al., 2020) :

where u s and v s are standardised versions of u s and v s to have variance equal to 1 (Simpson et al., 2017) . The term 0 ≤ φ ≤ 1 is a mixing parameter, which measures the proportion of the marginal variance explained by the structured effect. Finally, we assign to the temporal random effect term, w j , a Gaussian random walk model of order 1 (RW1), specified as:

The Bayesian representation of the above model is completed once we select priors for the fixed effects β 0 and β and the hyperparameters: τ , τ b , τ w , and φ. For the fixed effects we selected minimally informative Normal distributions whereas, we specified "penalising complexity" (PC) priors (Simpson et al., 2017) for the hyperparameters. PC priors are defined by penalising deviations from a "base" model (e.g., specified in terms of a specific value) and have the effect of regularising inference, while not implying too strong information. Technically, PC priors imply an exponential distribution on a function of the Kullback-Leibler divergence between the base model and an alternative model in which the relevant parameter is unrestricted. This translates to a suitable "minimally informative", regularising prior on the natural scale of the parameter.

In order to quantify the weekly excess mortality at sub-national level for specific age-sex population groups, we need to predict the number of deaths had COVID-19 not occurred. In Bayesian analysis, this can be performed by drawing random samples from the posterior predictive distribution (that is, the distribution of unobserved values conditional on the observed values from previous years). Specifically, letting θ θ θ be the model parameters, D be the observed data, and y jt * sk be the count of deaths that we want to predict, we have:

Operationally, we first generate random samples from the joint posterior marginal of the fitted linear predictor specified in equations (1) at the highest spatial resolution available (NUTS3 regions; Nomenclature of Territorial Units for Statistics 3 regions, https://ec.europa.eu/eurostat/web/nuts/ background/). Successively, we use these to be in turn the mean parameter of a Poisson distribution, to obtain the predicted number of deaths, which represents the baseline number of deaths assuming the pandemic did not take place.

Finally, to estimate the magnitude of excess deaths, the predicted number of deaths is compared against the actual observed number of deaths, allowing the computation of the relative change in the mortality (i.e., relative to what we could expect if the pandemic did not occur). This is obtained by (i) subtracting the predicted number of deaths from the observed number of deaths in each time point t and spatial unit s (number of excess deaths or NED), and (ii) dividing NED by the predicted number of deaths for each sample and multiplying by 100 (% relative excess mortality or REM).

Bayesian inference for the model is computed using Integrated Nested Laplace Approximation (INLA; Rue et al. (2009) , which performs approximate Bayesian inference on the class of latent Gaussian models (Rue and Held, 2005) . Unlike simulation based Markov chain Monte Carlo method, INLA is a deterministic algorithm, which employs analytical approximations and efficient numerical integration schemes to provide accurate approximations of the posterior distributions in short computing times. The INLA software is provided through the R package INLA, which can be downloaded from https: //www.r-inla.org/download-install. 

We retrieved all-cause mortality data during 2015-2020 in Italy from the Italian National Institute of Statistics (https://www.istat.it/). Data were available weekly (ISO week), by age (5-year age groups), sex and NUTS3 regions. As the COVID-19 mortality rates increase with age, we aggregated mortality counts based on the following age groups: <40, 40-59, 60-69, 70-79 and 80 years and older (Davies et al., 2021) .

Population data in Italy during 2015-2020 were retrieved from the Italian National Institute of Statistics. The data represent the population in Italy on January 1st of every year and are available by age (5-year age groups), sex and NUTS3 regions. To retrieve weekly population, we performed linear interpolation by the selected age groups (<40, 40-59, 60-69, 70-79 and 80+), sex and NUTS3 regions using populations at January 1st of the current and the next year. Population counts on January 1st 2021, which takes COVID-19 deaths in 2020 into consideration, were available at the time of analysis. Our goal was, however, to predict mortality for 2020 had the pandemic not occurred. Thus we performed an additional linear interpolation by age, sex and NUTS3 regions to predict the population at January 1st 2021, using the years 2015-2020 (Figure 2, panel A) . Object pop is a tibble containing the NUTS3 region ID (NUTS318CD), age group (ageg), sex (sex), year (year) and population counts ( We can use the following code (based on tidyverse and "piping" principles) to calculate the population of the 1st of January 2021 by NUTS3 regions, sex and age: pop %>% group_by(NUTS3, sex, age) %>% summarise(pop = as.vector(coef(lm(pop~year)) %*% c(1, 2021))) %>% mutate(year = 2021) -> pop2021

Once we obtained the year 2021 we performed an additional linear interpolation to calculate weekly number of population as shown on Figure 2 , panel B.

We used covariates related with ambient temperature and national holidays to help the model predictions. Data on air-temperature during 2015-2020 in Italy at 2m above the surface of land were retrieved from the ERA5 reanalysis data set of the Copernicus climate change program (Hersbach et al., 2020) . The geographical resolution of the ERA5 estimates is 0.25 • × 0.25 • (panel A of Figure 3 ). We calculated the weekly mean by the centroids of the 0.25 • × 0.25 • grid (panel B of Figure 3 ) and then averaged the weekly temperature over the ERA5 centroids that overlay with the NUTS3 regions (panels B and C of Figure 3 ).

The modelling process in INLA consists of three main steps: (1) the selection of priors; (2) definition of the model "formula" (which sets out the expression for the generalised linear predictor); and (3) the call to the main function inla, which computes the estimates.

In particular, we constructed the PC priors for σ = √ 1/τ ,

1/τ w based on statement that it is unlikely to have a relative risks higher than exp(2) based solely on spatial, yearly and seasonal variation, Figure 4 , panel A. For the mixing parameter φ, we set Pr(φ < 0.5) = 0.5 reflecting our lack of knowledge about whether overdispersion or strong spatial autocorrelation should dominate the field b, Figure 4, Figure 4 , panel B.

These assumptions can be encoded using the following code: # Defines the priors hyper.bym <-list( theta1 = list( PCprior , c(1, 0.01)), theta2 = list( PCprior , c(0.5, 0.5)) ) hyper.iid <-list(theta = list(prior="pc.prec", param=c(1, 0.01))) # Defines the model "formula" After fitting the model, we take 1000 samples from the (approximated) posterior distribution of the linear predictor and we use each drawn sample as the mean of a Poisson distribution to retrieve the predicted mortality counts: post.samples <-inla.posterior.sample(n = 1000, result = inla.mod) predlist <-do.call(cbind, lapply(post.samples, function(X) exp(X$latent[startsWith(rownames(X$latent), "Pred")]))) pois.samples <-apply(predlist, 2, function(Z) rpois(n = length(Z), lambda = Z)) This allows us to estimate the entire predictive posterior distribution of the mortality counts incorporating both the sampling and the linear predictor uncertainty.

To examine model validity, we performed a cross validation leaving out one historical year at a time and predicting the weekly number of deaths by NUTS3 regions for the year left out. For each stratum, we calculated the correlation between observed and fitted and a coverage probability, i.e. the probability that the observed death fall into the 95% credible interval (95% CrI) of the predicted.

We overall found good predictive ability of our model. 

The object pois.samples.list contains 1000 samples from the posterior predictive distribution 2, i.e. 1000 samples of the expected number of deaths by age, sex, NUTS3 regions and week, had the pandemic not occurred. We can access the different age-sex groups as follows:

names(pois.samples.list) # "F_less40" "F_40_59" "F_60_69" "F_70_79" "F_80plus" # "M_less40" "M_40_59" "M_60_69" "M_70_79" "M_80plus"

where F stands for females and M for males across the different age groups. To get an idea about the structure of the data, we can use the head ( We can also calculate median and 95% CrI expected number of deaths for this specific age-sex group over the year by NUTS3 region:

pois.samples.list$F_60_69 %>% select(starts_with("V"), "ID_space") %>% group_by(ID_space) %>% summarise_all(sum) %>% rowwise(ID_space) %>% mutate(median = median(c_across(V1:V1000)), LL = quantile(c_across(V1:V1000), probs= 0.025), UL = quantile(c_across(V1:V1000), probs= 0.975)) %>% select ( 

The above results can be combined in different ways using the functions get2020data() and get2020weeklydata() to calculate excess mortality. The function get2020data() aggregates over the entire country, NUTS2 regions, sex, age and time resulting in the object d, whereas the function get2020weeklydata() over the entire country, NUTS2 regions, sex and age but not over time resulting in the object d_week. names(d) # "province" "region" "country" names(d$province) # "none" "age" "sex" "agesex" names(d$province$age) # "40<" "40-59" "60-69" "70-79" "80+" names(d$province$sex) # "F" "M" names(d$province$agesex) # "F40<" "F40-59" "F60-69" "F70-79" "F80+" "M40<" "M40-59" "M60-69" "M70-79" "M80+" Notice that the object d$province$none is a simple feature collection, making mapping it straightforward.

In Figure 5 (plots 1A, 2A and 3A) we show the median posterior of REM for total age and sex at the national, NUTS2 and NUTS3 regional level. For these plots we used the object d with the selection "none" and plot the median REM (median.REM), for example for 3A: prov <-"Foggia" d$province$none %>% filter(DEN_UTS == prov) %>% select(geometry) %>% ggplot() + geom_sf(data = d$province$none, aes(fill = median.REM.cat)) + geom_sf(fill = NA, col = col.highlight, size = .8) + scale_fill_manual(values=colors, name = "", drop=FALSE) + theme_light() + ggtitle(paste0("3A. NUTS3 regions: ", prov)) Overall, the REM in Italy during 2020 was between 5-10%, Figure 5 , panel 1A. When the higher geographical resolution is assessed, it is revealed that north and in particular Lombardia is the region hit the worst, with the REM exceeding 20%, Figure 5 , panels 1B and 1C. Figure 6 shows a measure of uncertainty of the REM, now for the different age groups and both sexes (selection "age"). The probability of a positive excess (exceedance.REM) in older people and in the north of Italy is larger than 0.80, Figure 6 . Panels 1B and 1C of Figure 5 show the median temporal nationwide excess together with 95% CrI by sex after aggregating the different age groups (using d_week and the "sex" selection). We observe a clear first pandemic wave during March and May and a second one during mid October and December in 2020. During the first pandemic wave, there were weeks when the median relative excess death reached almost 100% in males, Figure 5 . Panels 2B, 2C, 3B and 3C of Figure 5 show the median spatio-temporal excess together with 95% CrI by sex after aggregating over the different age groups. Panels 2B and 2C highlight the region of Puglia, where during the first wave of the pandemic experienced a positive, similar with the nationwide excess. When we increase the spatial resolution in panels 3B and 3C, we highlight the region of Foggia, where there was insufficient evidence of a positive excess during the first wave, but strong during the second.

To be able to effectively examine and communicate the different aggregation levels of the output of our modelling framework, we have also developed a Shiny Web-Application (WebApp), Figure 7 . The WebApp provides spatial, temporal and spatio-temporal analysis tabs, and within each tab there are plots and summary statistics for the level of aggregation selected from the drop-down menu. Users can select across different variables (relative excess mortality or number of excess deaths), statistics (median or posterior probability), sex (males, females or both), age group (40<, 40-59, 60-69, 70-79, 80> and all) and different geographical level (national, NUTS2 or NUTS3 regions). Summary statistics for each area are available and they are displayed in a pop-up window, which is activated by clicking on the area of interest. In addition, graphical pop-ups are provided to show the weekly estimates for each area with the leafpop R-package (Appelhans and Detsch, 2021) , in the spatio-temporal analysis tab. The WebApp that we have developed is hosted at http://atlasmortalidad.uclm.es/italyexcess.

This tutorial provides a detailed implementation of the framework followed previously to calculate excess mortality during the COVID-19 pandemic in 5 European regions (Konstantinoudis et al., 2021) . We have proposed a Bayesian framework for estimating excess mortality and shown how to use R and INLA to retrieve fast and accurate estimates of the excess mortality. The proposed framework also allows combining different models and presenting the results in any age, sex, spatial and temporal aggregation desired. We have given a practical example of how to use the proposed framework modelling the excess mortality during the 2020 COVID-19 pandemic in Italy at small area level. We also developed a Webapp to effectively communicate the results. This framework can be extended to monitor mortality from other extreme events for instance natural hazards such as hurricanes (Acosta and Irizarry, 2020) . All the above make the proposed framework particularly powerful, generalisable and appealing for online monitoring of the pandemic burden and timely policy making. Potential extensions include different ways of modelling the younger age groups to increase the predictive ability, for instance using a zero-inflated Poisson. 

V.G.R. conceived the study. M.B. supervised the study. G.K. developed the initial study protocol and discussed it with M.B., M.C., M.P. and G.B.. G.K. developed the statistical model, prepared the population and covariate data and led the acquisition of mortality data. M.C. and V.G.R. validated and modified accordingly the code. G.K. ran the analysis. G.K., V.G.R., M.B., M.C. and M.P wrote the initial draft and all the authors contributed in modifying the paper and critically interpreting the results. V.G.R. developed the Shiny app. All authors read and approved the final version for publication.

A flexible statistical framework for estimating excess mortality. medRxiv

leafpop: Include Tables, Images and Graphs in Leaflet Pop-Ups

Spatial interactions and the statistical analysis of lattice systems (with discussion)

Bayesian image restoration, with two applications in spatial statistics

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

Community factors and excess mortality in first wave of the COVID-19 pandemic in England

The ERA5 global reanalysis

Excess deaths associated with COVID-19 pandemic in 2020: age and sex disaggregated time series analysis in 29 high income countries

Beyond the tip of the iceberg: direct and indirect effects of COVID-19

Discrete versus continuous domain models for disease mapping

Regional excess mortality during the 2020 COVID-19 pandemic: a study of five european countries. medRxiv

Magnitude, demographics and dynamics of the effect of the first wave of the COVID-19 pandemic on all-cause mortality in 21 industrialized countries

Lessons learned and lessons missed: impact of the coronavirus disease 2019 (covid-19) pandemic on all-cause mortality in 40 industrialised countries prior to mass vaccination

Excess mortality in england and wales during the first wave of the covid-19 pandemic

Small Area Disease Risk Estimation and Visualization Using R

An intuitive Bayesian spatial model for disease mapping that accounts for scaling

Excess deaths associated with COVID-19, by age and race and ethnicity-United States

Gaussian Markov random fields: theory and applications

Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations

Penalising model component complexity: A principled, practical approach to constructing priors

A linear mixed model to estimate covid-19-induced excess mortality

Hierarchical spatio-temporal mapping of disease rates

Estimation of excess deaths associated with the COVID-19 pandemic in the United States

All authors acknowledge infrastructure support for the Department of Epidemiology and Biostatistics provided by the NIHR Imperial Biomedical Research Centre (BRC). The authors also acknowledge infrastructure support for the domain of the shiny app from the University of Castilla-La Mancha. We thank Univesidad de Castilla-La Mancha for hosting the server on which the Shiny App is running.

Province stands for the NUTS3 regions (the resolution we used to fit the models), region for the NUTS2 (coarser than NUTS3, appropriate for policy making) and country for the nationwide aggregation. Within these aggregations users can select the option "none" being the total aggregation by age and sex, "age" by sex, "sex" by age and "agesex" refers to no age and sex aggregation. The objects d and d_week have similar structure and contain summary statistics for REM and NED and posterior probabilities of a positive REM or NED: 

The authors declare no competing interests.

The study is about secondary, aggregate anonymised data so no additional ethical permission is required.