Pandemic Data Quality Modelling: A Bayesian Approach

When pandemics like COVID-19 spread around the world, the rapidly evolving situation compels officials and executives to take prompt decisions and adapt policies depending on the current state of the disease. In this context, it is crucial for policymakers to have always a firm grasp on what is the current state of the pandemic, and to envision how the number of infections and possible deaths is going to evolve over the next weeks. However, as in many other situations involving compulsory registration of sensitive data from multiple collectors, cases might be reported with errors, often with delays deferring an up-to-date view of the state of things. Errors in collecting new cases affect the overall mortality, resulting in excess deaths reported by official statistics only months later. In this paper, we provide tools for evaluating the quality of pandemic mortality data. We accomplish this through a Bayesian approach accounting for the excess mortality pandemics might bring with respect to the normal level of mortality in the population.


Introduction
Historically, human populations have always dealt with major outbreaks of infectious diseases and effective reaction to these outbreaks has always been sought.The numerous black plague outbreaks in Europe in medieval and post-medieval times, the cholera pandemic in the late nineteen-early twentieth century, and the Spanish flu in the early twentieth century are only a few examples in modern history.However, the COVID-19 pandemic, which brought the world close to a halt in 2020 and 2021 and has killed almost seven million people as of early 2023, has brought to light a new reality of a global pandemic never experienced before.Worldwide governments initially underscored its urgency and their healthcare systems were quickly overwhelmed in a tsunami-like fashion.Nearly all of the affected countries progressively implemented measures to slow down the spread of the virus, ranging from recommended social distancing to almost complete lockdowns of social and economic activity.These measures eventually proved to be effective, allowing numerous states to relax restrictions, in an attempt to gradually return to normality.At the same time, with the threat posed by the virus still looming, decision-makers were forced to strike a balance between epidemiological risk and allowance of socioeconomic activity.In this type of context, surveillance of the number of new infections, mortality monitoring and quantification of the effects of social distancing became increasingly important [10,20,35], particularly so on a regional level.Given the local nature of the phenomenon, such a regional view appears to be of crucial importance.One of the difficulties lies in the fact that exact numbers of new infections, reported deaths and recoveries from the disease are often available only with a certain probability to be later corrected or are reported with a delay of -sometimes -several days.Moreover, there are discrepancies among data reported at different levels (regional, provincial, etc.) sometimes crucially causing wrong responses to face the pandemic.
Another important aspect to keep in mind is that the deaths due to COVID-19 nationwide are underestimated by official data because deaths that have not been tested are not counted.In addition, the pandemic had an indirect effect on mortality, preventing timely treatment for other diseases or limiting preventive examinations that could have anticipated critical situations.For all these reasons, official data on COVID-19 deaths in some countries were fairly inadequate to measure the effect of the pandemic [11].In aid of this, the World Health Organization (WHO) suggests studying excess mortality for assessing the death burden (both direct and indirect) of COVID-19.Excess mortality refers to the number of deaths from all causes during the pandemic more than what we would have expected under "normal" conditions, which is a valid measure of the total effect of the COVID-19 pandemic [2].Worldwide studies on excess mortality are deepened as reported in [28] and [33].Multiple model-based mortality estimates have been proposed in the literature, [2], [5], [24], [25], and [8] at national, regional and also county level.Several studies have shown a critical increase in excess mortality for higher age groups and for males, see for instance [17].
Excluding COVID-19, the ongoing global disease outbreaks reported by the WHO between 7/4/2022 and 7/4/2023 with at least one death are those listed in Table 1.Most of the diseases are endemic, especially in Africa and Asia.Diseases first reported in the last 25 years are, apart from Sars-Cov2, Mers-Cov (first reported in 2012), the Nipha virus (first reported in 1999) and a new type of hepatitis in children (first reported in 2022).Of the 16 diseases in an outbreak phase, 3 were first detected in the 18th century or before, 3 in the 19th century, 8 in the 20th century and 2 in the first 23 years of the 21st century.for example, [9] for possible new dangerous and unknown viruses generated from the thawing of permafrost in Siberia), new intensive farming and cultivation techniques (see, for example, [27] on how high population density in farming single animal species may integrate with other factors to increase the risk of mutation, re-assortment, and the generation of new pathogens) and increased migration, urbanization and the number of wars at a global level (see, for example, [18]), and probably with characteristics similar to that of COVID-19 in terms of speed of spreading, contagion rates and impact on societies, the pandemic "data challenge" should be one of the most important global tasks we must focus on.It is not a matter of chance that during the COVID-19 pandemic, among western countries, those most hit had problems with their data collecting systems.It is true that the viral characteristics of COVID-19 do not allow a clear country comparison with respect to the efficiency of facing the pandemic.In fact, [23] found that European and American countries were less efficient than South Asian and African countries, but they also admitted that this was mainly due to demographic features of the populations, hidden mortality in African and South Asian countries and the virus hitting mostly the elderly and fragile people.However, other viruses with similar contagion speeds and rates but hitting youngsters rather than the elderly might be hard to face if data collection is not efficient.
In this paper, we consider the bias between excess mortality and the official Italian COVID-19 data in the first 2020 outbreak for evaluating data quality in a space-time context.To model this bias we opted to use a Bayesian framework where two different quality measures ought to be evaluated: (i) the share in the population dying because of a particular infectious disease without being officially reported, so large values of this measure represent a worse scenario, and (ii) the coverage of the epidemic by the health systems, which can be considered an adequate indicator of their quality and a proxy for the efficacy of the crisis response.Among the factors explaining the bias variability, the focus is on detecting the most important component among spatial, temporal and interaction components.
The paper is organized as follows.Section 2 is devoted to a description of the data used, Section 3 focuses on proposed metrics for quantifying bias in official death data.In Section 4 we present the spatial, temporal, and spatial-temporal Bayesian model for the proposed metrics.The results are displayed in Section 5. Finally, Section 6 concludes the paper with an overview of possible future work.

Data
In order to evaluate the data quality on epidemic mortality in general and on COVID-19 mortality in particular, we considered two data sources: (i) official data on the pandemic evolution considering the essential variables to be considered for classic SE(I)RD models (daily or weekly new infections -i.e.cases -, susceptible, exposed, recovered and deceased people).More often similar data collection starts ad hoc at the request of national and international organisms to face the epidemic as was the case for the COVID-19 pandemic.In the following, we refer to this data as official data.(ii) National or supranational statistical institute data on population mortality.Henceforth, we refer to this data as ISTAT data as our main application will be on the Italian case and ISTAT is the Italian national statistical agency.
Official data In Italy, official data about deaths related to COVID-19 was reported daily from the 24th of February 2020 to the 30th of October 2022, at the European Union NUTS-2 level (i.e.regions), by the Italian Ministry of Health.While information about new cases was also published at the more refined NUTS-3 level (i.e.provinces), the number of daily COVID-19 new deaths was not officially available at this level.Nevertheless, it was possible to reconstruct the time series of COVID-19 at a NUTS-3 level indirectly, using other official sources like regional authorities' daily bulletins on provincial new cases, hospitalization and deaths and other information sources [14].Bulletins were in general in a "pdf" format so we were able to scrape data from these documents and to retrieve the data of interest for the majority of the Italian provinces.
In March 2022, ISTAT published data from Istituto Superiore di Sanità (ISS -the main public health institute in Italy) about the monthly evolution of COVID-19 deaths in each province.Using the regional weekly trends, it was possible to reconstruct an estimate for the weekly provincial COVID-19 officially reported deaths for the provinces which were missing from the scraped data.This technique has also been applied for some provinces for which the data had been scraped but the discrepancy with this monthly report was significant for some provinces, especially those experiencing low levels of mortality.
Official data about deaths caused by COVID-19 is assumed to have been subject to delays, errors, inconsistencies between reporting protocols, etc.Moreover, it is sensible to assume that the data has a systematic underestimation and is therefore biased as an estimate for actual COVID-19-related deaths: especially in the first pandemic stages, testing and care were not available for all people infected with the disease, and the official data only reflected the share of people that went through at least a minimal contact with the health system [6,7,30].Of course, this bias is supposed to change over time, as well as space, and in general, it is assumed to be directly proportional to the severity of the epidemic.
Additionally, COVID-19 did not only cause deaths directly because of infection, but the burden on the health system caused by the epidemic prevented to cure other diseases and accidents, causing indirect deaths of other individuals as well.Although lockdown measures might have prevented some of these "traditional" causes of death, this aspect should nevertheless be taken into account when assessing the impact on mortality caused by COVID-19.
With respect to the period under consider ISTAT data Each year, ISTAT provides a daily record of deaths reported in each municipality of Italy1 .In this study, the ISTAT data is aggregated by province and weekly: this is because the corresponding COVID-19 official data has a strong seasonality over the weekdays.This detailed time series data can help estimate the actual deaths in 2020 caused by COVID-19, both directly and indirectly.A straightforward way to estimate this is through the concept of "excess mortality", as the excess between the 2020 overall mortality and the average in the previous few years.
In this work, it is proposed to use a 5-year window consisting of the period 2015-2019 to represent the stable mortality level.The excess mortality is then found by subtracting this stable level from the 2020 deaths data, in each province and week of the year.While the 2015-2019 average will be smoother, the 2020 time series is subject to a great amount of noise, so it is chosen to apply minimal smoothing, consisting of a 3-week moving average, to make it more stable.
Since no other relevant information is available, this excess mortality can be imputed to COVID-19 and be an estimate for the actual number of deaths caused by COVID-19, both directly and indirectly.Indeed, in many works the indirect contribution of the excess mortality due to COVID-19 is highlighted, see, for example, [12], [1], [26] and [34].Mixing this estimate for the actual number of deaths with the official data allows for creating measures for the under-reporting bias present in the official figures.Such metrics can be interpreted as a proxy for the quality of the health system response to the crisis and to assess how this evolved and changed over time and space.In the end, the spatial distribution of these metrics can be extremely useful for the policy-maker to identify hotspots from the health system network, with the best and worst response to the emergency.Finally, studying the local policies and procedures implemented in those highlighted areas can help in the definition of best practices for future emergencies and epidemics.

Proposed metrics
The aim of the metrics to be defined is to provide an estimate for the under-reporting mortality bias that the official data has been subject to.However, there is no unique way to define such bias.Here, two different metrics with different interpretations for the policy-maker are proposed.

Additive bias b A
Let D ij be the officially reported total number of deaths in province i and week j which exceeds the average of the previous 5 years, assuming that they can all be imputed to the COVID-19 emergency.Let Dij be an estimate for D ij .Let Y ij be the officially reported number of COVID-19-related deaths in province i and week j.Finally, let P OP i be the average population in province i along the considered period.The additive bias is built as the difference between the actual mortality D ij /P OP i and the official mortality Y ij /P OP i .This bias is defined as "additive" because it must be added to the official mortality to get the unbiased value: Although, at least theoretically, we have that D ij ≥ Y ij , bA ij can take in practice negative values because Dij can take values lower than Y ij and even be negative by design.Moreover, even assuming that Dij correctly estimates D ij , errors and delays potentially occurred in the reporting of Y ij .Negative values of bA ij can be removed and/or treated as 0.
In terms of interpretation, b(A) defines the share in the population that died because of COVID-19 without being officially reported, so large values represent a negative scenario.Its trend over time and space (but not its magnitude) is a rough proxy for the part of the pandemic that was concealed and undetected by the public administration.Also, it can be roughly interpreted as the risk for an individual in a population of being infected and die of COVID-19 without having had the possibility to access the healthcare system because of capacity limits: this is of course dangerous for the population, since access to the healthcare system reduces the consequent risk of dying because of the disease.However, it is a personal risk because it does not consider the severity of the outbreak in each point in time and space, but it only assesses the remaining part of the pandemic that was unfortunately missed out by the competent authorities.
For example, two populations with the same b A = K have the same individual risk of dying of COVID-19 outside of the healthcare system, but the system itself may have to deal with two very different spread of the disease in their population so that the same value of b A can be considered decent result in one setting but utterly unsatisfactory in others.
Overall, it can be interpreted as the damage or impact of under-reporting on the population's well-being, but it is not an appropriate indicator for the quality of the healthcare system.In order to make this metric comparable with the following one, scaling is applied to the original formula:

Multiplicative bias b M
The ratio between Y ij and D ij assesses the probability of a COVID-19-related death being officially reported.
In order to transform it into a bias metric, the complementary probability of not being reported is considered instead and called b (M ) ij .This is equivalent to the additive bias divided by the excess mortality rate. b This metric should also stay theoretically between 0 and 1, but for the same issues cited before, it can exceed these boundaries.Again, a larger bias indicates a bad situation.Regarding its meaning, b M measures the coverage of the pandemic by the healthcare system, thus it is an adequate indicator of its quality and a proxy for the efficacy of the crisis response.
The quality of the response is assumed to have been at least partially correlated to the evolving severity of the pandemic, e.g., a good emergency policy might nevertheless have performed poorly at the peak of the epidemic, while an ill-advised response policy may have exceeded expectations because of the negligibility of the crisis it had to face.Hence, the distribution of the multiplicative bias over time and space shall be assessed both in the absence and presence of a covariate approximating the intensity of the pandemic.
Processed data In this study, we consider the period starting from the first day of COVID-19 official data release (24th February) until the 11th May, after 11 weeks, when the national lockdown was lifted in Italy.We opted to consider this first wave window, as underreporting and overall data quality are at their lowest at the beginning of an epidemic for obvious reasons, and they tend to improve over time.

Model
Inspired by the work of Franco-Villoria et al. ( 2022), [15], the model aims to consider the spatial, temporal and spatio-temporal structure (as well as potential fixed effects) for b A ij and b M ij .The objectives of such models are: assessing the overall spatial distribution of these quality indicators on the Italian territory, as well as the temporal trend in the first wave of COVID-19; assessing the significance of the interaction structure in the model; estimating the importance of each of these components, including also the contribution of the potential fixed effects.
The chosen specification is a latent Gaussian model on the logit transformation of the response with random effects.The same modelling structure can be applied to both response variables.In matrix notation: where µ represents the the mean value all times and spaces.u i and v j respectively represent the spatial and temporal main effects, and w ij is the interaction random effect.This structure allows us to clearly distinguish between the three main sources of variation under investigation, i.e. a spatial association, a changing pattern over time, and an interaction between these two elements.
The parameters of this model are µ, β, σ u , σ v , σ w , σ , while the covariance matrices Σ u and Σ v respectively define appropriate spatial and temporal models and are treated as fixed and known.The covariance matrix of the interaction term is defined as the Kronecker product of the covariance matrix of the two main effects, following the work of Knorr-Held, [22].

Temporal component specification
The objective of the model is to retrieve the main temporal pattern behind the response and assess the impact of this factor.Hence, a RW 1 model could be a simple but flexible choice, appropriate for the goal of the project.Note that this is not assumed to be the actual temporal model behind the data but it is exclusively used because of its convenience.

Spatial component specification
One of the most popular models for lattice data is the ICAR model [3] [4], which consists of an improper Normal distribution for the spatial component, defined using a square symmetric matrix M with only nonnegative entries and a null diagonal, and a corresponding diagonal matrix D where d p,q = q m p,q .
The weights m p,q represent the neighbourhood structure between the spatial areas, and it is usually based on the adjacency matrix, according to which regions that share a border are assigned m p,q = 1 and 0 otherwise.This approach is appropriate when the adjacency is the most reliable factor providing information about the potential association between regions.However, it can be argued that other variables should be considered instead, whenever they are more reasonable estimators of the actual connection between different areas.Adjacency matrices often consider geographical factors inappropriately: e.g.islands are isolated from the mainland, while geographical barriers such as mountains are ignored.This is particularly relevant in the context of epidemiology, where the main factor causing spatial association is human mobility.This is considered the prevalent factor behind spatial association so traditional definitions of the weight matrix M , such as adjacency or distance-based matrices are unable to capture factors such as commuting routes, physical boundaries, air and highways traffic, metropolitan areas, etc.
Hence, it is proposed to replace the traditional adjacency matrix with smartphone location data, which actually estimates the average commuting of individuals between two provinces, no matter their actual geographical location.The spatial weight matrix used in this application was derived starting from data provided by Pepe et al. [29], who built a daily time series of an origin-destination flow matrix processing smartphones' location data across Italy.The daily matrices from 18th January to 21st February were averaged to create a mean matrix representing the average mobility across provinces before the beginning of the COVID-19 outbreak.The matrices were built so as to have rows summing to 1, including the within province mobility.In order to respect the conditions necessary for the CAR model, the diagonal entries were set to 0 and the remaining entries were normalised through a division by row sum.The origin-destination matrix was not symmetric as the flows were directional: hence, the inward and outward flow between two provinces were averaged in order to obtain an appropriate M. The single disadvantage of this approach consists of the loss of the sparsity of the matrix.In order to obtain a matrix sparse enough for computation efficiency, it was chosen to consider only the last quintile, i.e. the highest 20 % of the weights, assuming irrelevant the impact of smaller associations.The resulting weights can be visualized in Figure 1: it is clear how this weight matrix is an improvement with respect to the naive adjacency one, as it is able to recognize important industrial hubs and big cities, as well as air and sea traffic, giving an overall accurate representation of the spatial structure.
A snippet of the M for the provinces of Lombardy is represented in Figure 2: this plot shows how the smartphone location data precisely reconstruct the human mobility patterns not only at the national level but also locally.
Figure 2: Weights m ij between the provinces of Lombardy

Prior specification
In terms of prior specification, the same approach of Franco-Villoria et al. [15], as the innovative hierarchical variance decomposition method introduced by Fuglstad et al. [16] is employed in the context of spatiotemporal epidemiological data.The method is based on a reparametrization of the variance parameters in terms of a single total variance and a set of proportions, defined through a decomposition tree.In this setting, it becomes easier to translate prior assumptions on the parameters into prior distributions and hyperparameters.Moreover, it gives more intuitive results in the posterior analysis, as the parameters are already set as adimensional proportional contributions, rather than variances or standard deviations.As in Franco-Villoria et al. [15], here we choose to define the first split to separate the main from the interaction term, and secondly, the spatial and temporal effects are divided.This results in a new reparametrization of the three original variances σ 2 u , σ 2 v , σ 2 w into a total residual variance V , the proportion ψ of this V given by the interaction term, and the proportion φ of main effects variance imputable to the spatial effect.The prior specification is then chosen on this new set of parameters.Specifically, the INLA default prior on variance parameters on σ 2 , a Uniform on φ, and a Penalized Complexity (PC) prior [32] on ψ with base model φ 0 = 0, and a PC prior on V with base model V 0 = 0.

Fitting the models for b A and b M
The models were fitted in R using the INLA software [31].The same model was fitted for b A and b M , with a difference in the prior specification in the U upper bound parameter of the PC prior on ψ, to reflect the different variability levels found in the two responses.Specifically, U A = 0.1 and U M = 1.More details can be found in [32].First of all, the different contributions of the random components to the total variability can be assessed considering the marginal posterior distribution of proportions of total variance V , see Figure 3.For both b A and b M , it appears that the spatial component is the most relevant, followed by the interaction effect, while the temporal trend explains a much smaller share of the total variability.Specifically, the spatial component explains around 50 % of the structured variability of b A , while around 60 % for b M .This first result already shows how the multiplicative bias b M is more influenced by the spatial structure than b A , which means it is more correlated with the intrinsic characteristics of the local areas, e.g.health system quality.With respect to the temporal pattern, the two metrics also show differences.The average temporal trend for b A , shown in Figure 5, green curve, starts with an increasing part, up to the sixth week in the considered period, followed by a steady fall in the remaining weeks.This is a reasonable result, as it is expected that the indicator b A performed the worst at the peak of the "official" epidemic evolution, plus a delay due to the fact that deaths are considered instead of cases.Hence, this confirms the assumption that b A is related to the level of stress of the health system, rather than to the quality of its response to a certain amount of stress.The results for b M are again completely different as the posterior means, shown in Figure 5, red curve, show a steadily decreasing trend.It is important to note that the level is at its maximum at the beginning of the data reporting period and then, the quality level of the emergency response increases as the COVID-19 situation is taken more and more seriously and more effective policies and practices are put in place.Of course, this shows how the official data are unreliable, not only in magnitude but also in their trend, and should not be used raw to evaluate the evolution of an epidemic, particularly at the beginning of the reporting period, as the quality of the official data tends to significantly improve, i.e. a decrease in the b M indicator.While the temporal and spatial effects can offer general evaluations of the overall evolution, the interaction factor seems to play a vital role in the total variability, as shown in Figure 3. Hence, any conclusion on the single provinces' behaviour should be made on the whole linear predictor.In the following figure, we considered the posterior distribution of the fitted values reparametrized in the original scale of the b metrics, through a logistic transformation.The provinces were selected among the most populous areas in Italy, along with provinces particularly affected by COVID-19, and provinces displaying surprising behaviours in the results.The importance of the interaction term appears clearly in Figures 6 and 7 and, as the provinces display quite a variety of temporal evolution behaviours.In order to classify provinces on the basis of the b M metric, i.e. a data quality indicator, clustering was applied to the posterior means of the fitted values from Model b M .Specifically, partitional clustering using a dynamic time warping distance was chosen as an appropriate partitioning method for time series [13].The posterior mean of the fitted values from the model were chosen as the time series, in order to remove the residual noise present in the original raw data.Evaluating different performance metrics, the optimal cluster number was found to be 4. Figure 8 shows the 4 different classes and their centroid, ordered by best to worst overall performance.
The first cluster, coloured in green, includes all the provinces with a steady and relevant decrease over time; moreover, all of them except one also even start from a relatively low level already at the beginning of the reporting period.The red and last cluster identifies the worst-behaved class, as the index b M is already quite high at the beginning of the period is mostly stable, and does not display significant signs of improvement during the first wave.
The other two clusters represent intermediate behaviours and they are more ambiguous and heterogeneous.First, the orange cluster identifies provinces that performed better than the ones in the red cluster, but worse than the green provinces, because of a much lower rate of decrease.There are two main subtypes of trends in this cluster.Most of the provinces display an initial decrease, followed by a mild rise in the last few weeks, which could be a sign of a decline in alertness and a possible relaxation of some emergency measures initially put in place.Few provinces in the same cluster actually followed an inverse pattern, with an initial stable level of b M followed a quite steep improvement (e.g.Udine, Pordenone, Verona, Belluno, Vicenza, Padova, Venezia, Treviso), as if the severity of the epidemic was downplayed after the initial shock.The two groups have in common an insufficient degree of vigilance and promptness in dealing with an epidemic emergency, which became manifest as either a delay in putting in place effective response policies or a relaxation of the same way too early.
For the second group, the general trend, followed by the majority of the units in this cluster, is decreasing but more slowly than for the green provinces, hence indicating a worse performance in data reporting.However, looking more closely, we can identify a small subgroup of provinces in this cluster that displays an unusual concave U-shape trend: Sondrio, Verbano-Cusio-Ossola, Monza e della Brianza, Varese, Lecco, Como.The areas with this distinct behaviour are all located in North-West Italy and all fall into this class because they do not tend to reach extremes, either positive or negative, but are more or less stable around the 0.2 threshold.One hypothesis could be that these provinces might have experienced an actual delay in the spread of the epidemic with respect to the general evolution, which might have led the peak to be shifted later than the other province, for about 5 or 6 weeks.This delay was probably responsible for keeping b M low even in the worst moments for these provinces, as they had a temporal advantage consisting of all the experience and knowledge accumulated and shared by the other health system administrations over the first few weeks.The map presented in Figure 9 summarizes the cluster assignments over the Italian provinces, along with the location of the first COVID-19 hotspot in Italy in Codogno, Lodi.At first glance, it may seem counter-intuitive that areas closer to the first outbreak performed better than the ones far away.However, it may be that provinces in the North considered the emergency more seriously and sooner than the rest of Italy, as they perceived the risks earlier and may have put in place effective response protocols before the Southern areas.In this study, we considered two metrics to evaluate the quality of epidemic official data.This has been achieved using a measure of excess mortality registered during an epidemic period, with respect to the average level in the previous years.While the metric b A is a proxy for the underestimation of the epidemic severity, it has been shown how b M can be used as a quality indicator of the health system in terms of monitoring the epidemic situation and reporting accurately the data.Both metrics were smoothed through a spatio-temporal model, which highlighted the importance of the spatial components in the variability of both responses.Results for the b M are consistent with the assumptions about the data quality: in particular, the general temporal trend shows a steady decrease and it is in line with an improvement in the data quality in the first period of reporting.
Finally, the Italian provinces were grouped on the basis of the b M indicator fitted values, in order to classify them into clusters of different quality levels.This could be extremely useful for policymakers, since it offers an overview of the reporting performance of the different health systems and may suggest strategies employed locally, which may become best practices for future epidemics.
Future work will entail an implementation of the quality metrics to other epidemics and nations.We may be interested in comparing the overall data quality between nations facing the same epidemic, rather than within a single nation.With respect to the Italian case, the time series could be extended to the following waves if official data about COVID-19 were available at the provincial level.Covariates related to the health system could be included in the models as fixed effects, to check whether they reduce the spatial variability and may explain the b M metric.
In conclusion, the multiplicative bias is the main proposal of this work and it consists of a simple metric that could in theory be computed with a really short delay, i.e. as soon as registered deaths are aggregated by the governmental institutes for official statistics.Monitoring this indicator could help identify which areas display better results and performance in terms of the response to the epidemic, as well as areas where the severity of the situation is underestimated.This type of information could significantly accelerate the process of identification of effective and ineffective protocols and prevention measures, which consequently may save many lives down the line.

Figure 1 :
Figure 1: Mobility-based weight matrix in log scale

Figure 3 :
Figure 3: Posterior densities of the variance parameters as proportions of the total residual variance V in the b A and b M models

Figure 4 :
Figure 4: Posterior mean of the spatial random effects on b A and b M

Figure 5 :
Figure 5: Posterior mean of the temporal random effect on b A and b M

Figure 6 :
Figure 6: Posterior mean of the fitted values for b A in the original scale for selected provinces

Figure 7 :
Figure 7: Posterior mean of the fitted values for b M in the original scale for selected provinces

Figure 8 :
Figure 8: Posterior mean of the fitted values divided in 4 clusters with corresponding centroids in the provinces of Aosta, Rimini, Catanzaro, and Cosenza

Figure 9 :
Figure 9: Provinces by cluster and the city of Codogno (LO) marked by a black dot

Table 1 :
Disease outbreaks with at least 1 death reported by WHO.7/4/2022 -7/4/2023 (Source: Authors' calculations from WHO data) Therefore, as in the future new pandemics are going to spread due to new issues like global warming (see,