A statistical analysis of death rates in Italy for the years 2015-2020 and a comparison with the casualties reported for the COVID-19 pandemic
Gianluca Bonifazi, Luca Lista, Dario Menasce, Mauro Mezzetto, Alberto Oliva, Daniele Pedrini, Roberto Spighi, Antonio Zoccoli
AA statistical analysis of death rates in Italy for the years 2015–2020 anda comparison with the casualties reported for the COVID-19 pandemic.
Gianluca Bonifazi , Luca Lista , Dario Menasce , Mauro Mezzetto , Alberto Oliva , DanielePedrini , Roberto Spighi , and Antonio Zoccoli Universit`a Politecnica delle Marche INFN Sezione di Bologna Universit`a degli Studi di Napoli Federico II INFN Sezione di Napoli INFN Sezione di Milano Bicocca INFN Sezione di Padova Alma Mater Studiorum Universit`a di Bologna * Corresponding author , e-mail: [email protected]
Abstract
We analyze the data about casualties in Italy in the period 01/01/2015 to 30/09/2020 released by the ItalianNational Institute of Statistics (ISTAT). The data exhibit a clear sinusoidal behavior, whose fit allows for arobust subtraction of the baseline trend of casualties in Italy, with a surplus of mortality in correspondence tothe flu epidemics in winter and to the hottest periods in summer. While these peaks are symmetric in shape, thepeak in coincidence with the COVID-19 pandemics is asymmetric and more pronounced. We fit the former with aGaussian function and the latter with a Gompertz function, in order to quantify number of casualties, the durationand the position of all causes of excess deaths. The overall quality of the fit to the data turns out to be very good.We discuss the trend of casualties in Italy by different classes of ages and for the different genders. We finallycompare the data-subtracted casualties as reported by ISTAT with those reported by the Italian Department forCivil Protection (DPC) relative to the deaths directly attributed to COVID-19, and we discuss the differences.
Last year has seen an enhanced attention of the media concerning the death rates in various countries in the World,a fact due to the outbreak of the COVID-19 pandemics at the beginning of 2020 and the ensuing alarm it generatedworldwide. An analysis of the historical data archives[1] publicly made available by the Italian Istituto Nazionaledi Statistica (ISTAT)[2] shows a periodic variation of the death rate depending upon seasons, well represented by astable and regular sinusoid. Superimposed to the sinusoid trend there may be additional death excesses, most likelydue to seasonal diseases like influenza in winter or to very intense heat waves in the summer [3, 4].The approach adopted here for such an estimate is an interpolation of the data with a fit function exhibiting an ad hoc modeling of the main features of the curve. We present results of this method applied to the data providedby ISTAT in the period 2015–2020. A comparison is then carried out with another data sample[5], provided bythe Dipartimento della Protezione Civile (DPC) [6], which provides counts for deaths attributed to the COVID-19disease.
This study is based on publicly available data provided by ISTAT [1] as time series of recorded deaths by the NationalRegistry Office. The data, collected from all the 7903 districts located in the 20 Italian regions, covers the periodfrom January 1 st th a r X i v : . [ phy s i c s . s o c - ph ] F e b igure 1: The distribution of deaths collected by ISTAT [1] from 7903 districts in Italy between the years fromJanuary 1 st th . This feature remains partly confirmed also by disentangling the datausing age as a selection criteria. Fig. 2 shows the distribution of deaths for people in three different age classes: inblue those below 50 years, in green those in range 50 to 79; in red those above 80 and in magenta the sum of all thesethree classes. It is evident from these distributions that people with an age below 50 die, to a good approximation,with a constant average probability in any given day of the year while those above that age tend to have a varying,periodic probability of death with a maximum in winter and a minimum in summer. The older the age, the larger theexcess of death in particular periods of the year, appearing in the form of Gaussian-like excesses over the sinusoidalwave. Disentangling the data by gender, see Fig. 3, there seems to be a slight prevalence of female deaths withrespect to males, except for the COVID-19 peak, where the situation happens to be reversed. These are just rawvalues, though, not corrected to take into account the ratio between males and females in the Italian population.Later on we will quantify and appropriately weight these data. We perform a global fit over the data, where in a single step we estimate the sinusoidal baseline of the distribution,the seasonal excesses and the 2020 peak in correspondence of the COVID-19 pandemic. This method differs fromother methods often reported in literature, where the background is subtracted by computing the average numberof counts in the same period of the past 5 years, as in [8], or by fitting the sinusoid with the data during periods oftime where the excesses are not visible, like in spring or autumn [7, 9].We used a χ fit to interpolate the data with an appropriate function meant to model the data in order to determinethe value of the unknown parameters of the model along with their uncertainties. The actual minimization is carriedout by the MINUIT [10] package.The overall fit function has been defined as the sum of individual components in the following form: In the following we will discuss how we established that there is no significant slope of the average value of this wave. ( t ) = s ( t ) + k (cid:88) i =1 g i ( t ) + ˙ G ( t ) (1)where s ( t ), g i ( t ) end ˙ G ( t ) are defined and described below.The s ( t ) function is meant to model the wave-like variation of deaths with seasons, the g i ( t ) function describesthe excess peaks visible above the wave and ˙ G ( t ) represents the rightmost excess peak (spring 2020), which, unlikethe others, is asymmetrical. The index i runs from 1 up to k , the number of excess peaks featured by the datadistribution except the last one on the far right.The general wave-like behavior of the data is modeled by a sinusoidal function of the form: s ( t ) = c ( t ) + a sin (cid:18) πtT + ϕ (cid:19) (2)where t is the day number starting from t = 1 / / c ( t ) represents the slowly-varying offsetfrom zero deaths, a the amplitude of the oscillation, T is the period of variation (the time delay between consecutivemaxima) and finally ϕ the phase. We tried to model c ( t ) allowing for a linear dependence on t , as c ( t ) = c + c t ,but the fit determines a slope c compatible to zero within uncertainty. We therefore decided to maintain the c termindependent of time in the final fit.Each individual excess above this s ( t ) wave can be modeled by a Gaussian distribution of the canonical form: g i ( t ) = N i √ πσ i e − ( t − µ i ) / σ i (3)The choice of a Gaussian function here is only justified by being the simplest symmetrical function to describethese excesses representing, at the same time the distribution of a random variable. Modeling the excess peaks in thedescribed way has the advantage that the individual g i fit parameters correspond to a Gaussian with the backgroundcontribution already taken into account in the overall fit model. The N i parameter of each Gaussian, corresponds tothe number of excess deaths with respect to the wave-like background, whose values are also determined optimallyby the fit itself. An advantage of this approach is that in the case of adjacent, overlapping Gaussians (as can be seenin Fig. 4 in the case of the g end g peaks but also the g and the big peak on the far right of the distribution),each individual area is computed correctly by taking into account the nearby contributing ones.While the excess peaks look highly symmetrical around their maximum and can thus be reasonably well modeledwith Gaussians, as described before, the peak of the spring 2020, associated with the COVID-19 pandemic, is clearlyasymmetric. We have tried several possible parametrizations for that distribution, such as bifurcated Gaussianswith a common peak, generalized logistics, or else, to reflect the asymmetry, but in the end we resolved to adoptthe derivative of a Gompertz function[11] simply because it is customarily adopted by epidemiologists to describeepidemic evolution’s over time and we therefore considered it more suitable to our purpose.A Gompertz function is parametrized in the following way: G ( t ) = N G e − b e − ht (4)Equation 4 represents a cumulative distribution. Since our data represent instead daily counts, we used itsderivative, given by: ˙ G ( t ) = d G ( t )d t = N G bh e − b e − ht − ht (5)where the parameter N G is the value of the integral of this function. In Figure 4 and Tables 1 to 3 we report the results of a fit to the whole data sample, comprising both genders for allages in the six years from 2015 to 2020. The column labeled ‘ µ i ’, in Table 1, indicates the day when the maximumof an excess has been reached while those labeled ‘ µ i ± σ i ’ indicate, respectively, the day of onset and demise fromthe average background, a time interval in which occur 95% of the death cases (expressed with calendar dates). Thecolumn labeled ‘ Duration ’ is the time difference between onset and demise (namely 4 σ , expressed as number days).4igure 4: The whole data sample with a superimposed fit function obtained as specified in the text. The plot at thebottom shows the pulls (a quantity defined later on in the text) of the fit: their mean value, being compatible withzero, is a testimony of the appropriate choice of the particular model adopted. Numerical values for the integral ofeach individual Gaussian are provided in Table 1. g i N i µ i (days) σ i (days) µ i − σ i (date) µ i (date) µ i + 2 σ i (date) Duration (days)1 50 706 ± ± ± ±
361 201.1 ± ± ±
527 381.6 ± ± ±
534 743.5 ± ± ±
210 950.1 ± ± ±
704 1104.9 ± ± ±
616 1155.0 ± ± ±
256 1313.3 ± ± ±
685 1492.5 ± ± ±
504 1642.3 ± ± ±
613 1853.3 ± ± ±
217 2007.5 ± ± ±
362 2046.3 ± ± c a T (days) ϕ (rad)1678 ± ± ± ± pulls , p i , are defined as: p i = d i − F ( t i ) (cid:15) i (6)where d i is the number of death counts in a given day i and (cid:15) i the corresponding amount of statistical fluctuation.The data, being outcomes of counting values, are assumed to follow a Poisson distribution, hence (cid:15) i = √ d i .The χ /n DOF of the fit turns out to be 3.271.We report the distribution of the pulls in Fig. 5 fitted with a Gaussian function. The mean value of the fit is − . ± .
04, compatible with zero, while the standard deviation of the Gaussian fit turns out to be 1 . ± . i is given by the fit parameter N i defined in eq. 3, while the area of theGompertz derivative is the fit parameter N G in eq. 4. Yield From Peak To Duration (days)54 387 ±
557 20/2/2020 24/3/2020 11/5/2020 81Table 3: Results of the fit to the whole data set (no filters applied) for the Gompertz derivative function. Themeaning of the columns labeled
From , Peak and To is explained in the text.The width of the Gompertz is computed from the first day in which the integral of the function exceeds 2.5% ofthe total to the day in which the integral reaches 97.5% of the total. These two days are reported in Table 3 underthe columns labeled From and To .Figure 5: Distribution of the pulls (depicted as a time series in the bottom plot of Fig. 4) fitted with a Gaussianfunction. The peak of the Gaussian is at µ = − . ± .
04 (therefore compatible with zero) while the width id givenby σ = 1 . ± . T = 364 . ± . c = 1678 ± . st of everyyear.These results highlight an interesting feature of the COVID-19 deaths excess. As already noted, almost everywinter there is an excess of deaths with respect to the baseline, with the notable exception of the years 2015–2016(a period with a particularly balm winter[18], with a relatively small value of 4 455 excess of casualties). The peakin the spring of 2020, instead, shows characteristics markedly different from the winter excesses of previous year interms of amplitude, width and day of the year when the maximum is reached. In the following we will mention thepossible implications of these differences. 6 Age and gender mortality
We have also disentangled the data by age and gender and fit the distributions in these different categories to obtainaccurate numerical values.Figure 6: Casualties of people with age in the range 50 ≤ age ≤
59 with a superimposed fit based on function 6(blue points) while gray points are the category 0 < age ≤
49. Numerical results are listed in Table 4 and 5. Theplots at the bottom show the pulls of the two samples.We start with a cumulative plot for all people aged between 50 (included) and 60 years (excluded) who diedbetween 2015 and 2020, shown in Fig. 6. This plot shows that the average value of daily deaths for people in thisage range is about 70 casualties/day. In order to get a fit comparable with the one in Fig. 4, we are forced to adopta somewhat more stringent fit strategy. The wave parameter corresponding to the phase has been fixed to the value g i Yield ± ±
583 508 ±
714 910 ±
715 166 ±
346 724 ±
717 373 ±
608 80 ±
399 1105 ± ± ± ± ± s ( t ) function). This guarantees that maxima are reached in the winter and minima in the summer and no spurious time translation is introduced by thefit procedure when a local minima can eventually be found. Also the peak position and the width of the 13 Gaussians7ave been fixed to the values established for the whole data sample while the Gompertz parameters are all left free.The corresponding fit results are listed in Tables 4 and 5. Yield From Peak To Duration (days)1 373 ±
87 28/2/2020 24/3/2020 1/5/2020 63Table 5: Results of the fit to the whole data set (no filters applied) for the Gompertz derivative function. Thesevalues correspond to the fit depicted in Fig. 6.The picture shows two categories of age at the same time: those in the range 0 −
49 (in gray) do not show any signof seasonal variation around the mean value of ∼ /day (they were fit with a simple constant term). A sinusoidalvariation begins to be noticeable only in the range 50 −
59 (blue dots), along with the presence of the correspondingdeath excesses indicating a continuous increase in magnitude with age starting around 50. The results are affectedby larger uncertainties with respect to the full sample of Fig. 4, reflecting the smaller size of population in this range.The excess peaks and the sinusoid amplitude become more evident in a sample of even higher ages, namely60 ≤ age <
80. The average number of deaths in this category is much larger, due to an enhanced health fragilityfor people of progressively higher age, as seen in Fig. 7.Figure 7: Casualties of people with 60 ≤ age <
80 (with a superimposed fit based on eq.1). Numerical results arelisted in Table 6 and 7The fit is again pretty similar, in shape but not in amplitude of course, to the full sample shown in Fig. 4. The pulls feature a mean value compatible with zero also in this case. The fit strategy is the same as the one describedbefore for Fig. 6. Values obtained in this case are listed in Tables 6 and 7. The average death rate in this category is ∼ /day . Increasing the age threshold further up, by collecting deaths of people aged ≥
80, we get a sample withvery pronounced peaks, see Fig. 8. The average death rate in this last category reaches the high value of ∼ /day .Another information can be extracted from the data is the relative amount of deaths between genders. Fig 9show the distribution of males and females (summed over all ages) superimposed with the relative fits. In this case,since the two samples have a rather large statistical amount, both fits have been performed with all parameters freeto vary. These numbers need to be corrected by the relative number of males and females in the Italian population.The fraction of males in 2020 was 48 .
7% while females were 51 .
3% [12]: we compute a mortality factor (for eachgender) by normalizing the yields to 29 050 096 and 30 591 392 (the respective number of males and females of thetotal Italian population by January 1 st Mortality . While the absolute number of female deaths is higher than the males one every year, theopposite seems true for the 2020 peak. After re-weighting this small discrepancy between genders remains basically8 i Yield ± ± ±
84 8 432 ± ± ± ± ± ± ± ± ± ± ≤ age <
80 . These values correspond to the fit depicted in Fig. 7
Yield From Peak To Duration (days)15 951 ±
242 25/2/2020 22/3/2020 25/4/2020 60Table 7: Results of the fit for the category of 60 ≤ age <
80 for the Gompertz derivative function. These valuescorrespond to the fit depicted in Fig. 7Figure 8: Casualties of people with age ≥
80 (with a superimposed fit based on function 6). Numerical results arelisted in Table 8 and 9 9 i Yield Peak Duration (days)1 28 324 ±
721 10/1/2015 1762 8 844 ±
254 21/7/2015 453 630 ±
220 17/1/2016 424 23 732 ±
344 13/1/2017 605 4 655 ±
179 9/8/2017 206 11 978 ±
343 11/1/2018 637 4 367 ±
267 3/3/2018 418 3 607 ±
233 8/8/2018 339 18 217 ±
425 4/2/2019 10510 8 981 ±
379 4/7/2019 10711 6 648 ±
389 31/1/2020 9512 3 208 ±
238 4/7/2020 3713 6 396 ±
270 11/8/2020 53Table 8: Results of the fit to the data set of people aged over 80. These values correspond to the fit depicted inFig. 8
Yield From Peak To Duration (days)37 357 ±
365 22/2/2020 25/3/2020 5/5/2020 73Table 9: Results of the fit to the data set of people aged over 80 for the Gompertz derivative function. These valuescorrespond to the fit depicted in Fig. 8Figure 9: Number of daily casualties for males and females of all ages.10rue for all peaks except the 2020 one, where the mortality turns out to be larger for males than for females.The fraction of casualties for the two genders turn out to be about the same, at the level of one standard deviationin all the years till 2019 included.
Males Females g i Yield Mortality Peak Duration Yield Mortality Peak Duration ±
788 69.8 10/1/2015 178 32 405 ±
929 105.9 10/1/2015 2062 5 132 ±
258 17.7 21/7/2015 53 8 401 ±
266 27.5 21/7/2015 493 1 734 ±
394 6.0 17/1/2016 76 3 198 ±
441 10.5 17/1/2016 804 13 907 ±
357 47.9 13/1/2017 62 20 281 ±
404 66.3 13/1/2017 665 2 252 ±
156 7.8 9/8/2017 19 3 976 ±
174 13.0 9/8/2017 216 8 875 ±
468 30.6 11/1/2018 83 11 830 ±
415 38.7 11/1/2018 717 2 192 ±
268 7.5 3/3/2018 33 4 139 ±
305 13.5 3/3/2018 438 1 612 ±
194 5.5 8/8/2018 31 2 473 ±
203 8.1 8/8/2018 279 11 066 ±
462 38.1 4/2/2019 99 15 273 ±
509 49.9 4/2/2019 10910 3 728 ±
341 12.8 4/7/2019 99 5 224 ±
352 17.1 4/7/2019 8711 4 112 ±
443 14.2 31/1/2020 99 5 124 ±
439 16.7 31/1/2020 9312 964 ±
182 3.3 4/7/2020 29 1 764 ±
182 5.8 4/7/2020 2713 2 489 ±
236 8.6 11/8/2020 47 3 826 ±
253 12.5 11/8/2020 49Table 10: Results of the fit to the data set divided in a sample of males and another of females (of all ages) in Fig. 9.The meaning of the
Mortality column is described in the text.
Males Females
Yield Mortality Peak Duration Yield Mortality Peak Duration
27 240 ±
366 93.8 22/3/2020 63 26 079 ±
395 85.2 26/3/2020 72Table 11: Results of the fit to the data set (for the Gompertz peak only) divided in a sample of males and anotherof females (of all ages) in Fig. 9. The meaning of the
Mortality column is described in the text.
The data set provided by ISTAT [1] and used for the present analysis is not the only one publicly available: anotheris provided by the Dipartimento della Protezione Civile (DPC) [5], which provides a somewhat different kind ofinformation regarding the number of deaths in the context of the COVID-19 pandemic. In particular, the datarecord, which begins February 24 th directly attributed to the currentpandemic.A plot of the data from these two disparate sources is shown in Fig. 10. The magenta points (and the accompanyingfit result of a Gompertz function in red) correspond to the ISTAT data sample: these data are a subset of thosedisplayed in Fig. 4, specifically those between the dates of February 24 th and September 30 th th , due to the fact that a certain number of deaths werenot correctly reported in the preceding weeks and were recovered assigning that day as the actual death date. Inorder to compare the yield returned by the fit to the value provided by the ISTAT data we had to exclude thecontributions from the second pandemic peak: we decided to introduce a cutoff value while computing the sum ofentries of the DPC sample in correspondence to August 16 th , a day when the minimum number of casualties wasreached between the two pandemic waves, therefore including also the spike. The cutoff date is shown in Fig. 10 asa vertical green arrow.The sum in that period (the blue dots) results to be 35 468.Figure 10: Comparison between ISTAT and DPC data samplesOn the other hand, the yield obtained for the ISTAT sample is the one reported in Table 3, namely 54 387 ± ± The rich data sample provided by ISTAT allows for various additional visualizations. In Fig. 11 we display data forages in groups of 4 years to visualize the increase of death probability with age: it becomes more evident what wasalready shown in Fig. 2, namely the fact the young people tend to die with a rather flat probability along each year,while people of a progressively higher age tend to suffer more from illnesses in specific seasonal periods. Each bin inthis plot constitutes the number of deaths in six contiguous days. In Fig.12 we present a scatter plot of death ratesas a function of the day of the year (for the six years from 2015 to 2020) versus the age category. This graphicalrepresentation clearly illustrates the higher probability of death for the age category 70–95 with respect to the others.Each value in the ISTAT data sample comes with a geographical tagging marker, allowing for a categorization ofthe number of deaths in different parts of Italy.Fig. 13 shows the fits for each of the four zones in Italy, namely
North , Center , South and
Islands . The fits in The subdivision is arbitrary and we have defined
North as the sum of values for the following regions:
Piemonte, Valle d’Aosta,Liguria, Lombardia, Trentino-Alto Adige, Veneto, Friuli-Venezia Giulia, Emilia-Romagna . Center comprises T oscana, Umbria, Marche, g i North Central South Islands ±
518 13 293 ±
366 10 670 ±
333 5 730 ± ±
193 3 409 ±
139 2 995 ±
126 999 ±
913 1 333 ±
235 1 501 ±
166 1 250 ±
151 1 120 ± ±
254 9 184 ±
183 7 058 ±
165 3 135 ± ±
106 2 078 ±
85 1 716 ±
77 1 066 ±
586 9 393 ±
246 4 027 ±
170 3 676 ±
156 2 344 ± ±
203 2 233 ±
144 1 576 ±
129 1 451 ± ±
135 876 ±
93 511 ±
82 274 ±
619 10 692 ±
302 6 865 ±
213 6 620 ±
197 3 884 ± ±
255 2 631 ±
180 2 168 ±
161 1 374 ± ±
264 2 505 ±
186 2 033 ±
171 1 361 ± ±
119 898 ±
87 871 ±
79 300 ± ±
182 1 286 ±
131 1 188 ±
117 1 026 ± ±
33 6 063 ±
201 3 559 ±
205 2 158 ± Lazio, Abruzzo and Molise , South includes
Campania, Puglia, Basilicata and Calabria . Finally
Islands corresponds to
Sicilia andSardegna i North (%) Central (%) South (%) Islands (%) ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± The data provided by ISTAT allow for a detailed quantitative estimate of the number of deaths excesses with respectto an average baseline. This baseline is represented by a wave-like variation of the number of deaths which turnsout to be almost perfectly in phase with the yearly seasonal cycle. In this paper we present a study of these excessesevaluated by a statistical interpolation of the data based on a χ minimization method using a function which is thesum of a sinusoidal wave (to represent the repeated, seasonal variation of deaths) a number of Gaussian distributionsto represent the excesses above the sinusoid and, finally, a Gompertz derivative to model the asymmetric peak ofspring 2020. In this study we discussed the methodology adopted for the interpolations and analyze different samplesby disentangling genders, ages and locations. A comparison has also been carried out between the number of deathsprovided by ISTAT from the National Death Registry in the period corresponding to the first wave of the pandemicand the numbers provided by DPC in the same period for the deaths directly attributed to COVID-19. We found arather large discrepancy, amounting to ∼ ∼ covid19.infn.it/ . We are grateful to Mauro Albani, Marco Battaglini, Gianni Corsetti and SabinaPrati of ISTAT for useful insights and discussions. We wish also to thank Daniele Del Re and Paolo Meridiani foruseful discussions. The project has been supported in various ways by a number of people from different INFNUnits. In particular, we wish to thank, in alphabetic order: Stefano Antonelli (CNAF), Fabio Bredo (Padova Unit),Luca Carbone (Milano-Bicocca Unit), Francesca Cuicchio (Communication Office), Mauro Dinardo (Milano-BicoccaUnit), Paolo Dini (Milano-Bicocca Unit), Rosario Esposito (Naples Unit), Stefano Longo (CNAF), and Stefano Zani(CNAF). We also wish to thank Prof. Domenico Ursino (Universit`a Politecnica delle Marche) for his supportivecontribution. References [1] Public data on deaths in the Italian municipalities for the years 2015-2020, provided by ISTAT, [3] Hajat, S., & Gasparrini, A. (2016). The excess winter deaths measure: why its use is misleading for public healthunderstanding of cold-related health impacts. Epidemiology (Cambridge, Mass.), 27(4), 486.[4] Weinberger, K. R., Harris, D., Spangler, K. R., Zanobetti, A., & Wellenius, G. A. (2020). Estimating the numberof excess deaths attributable to heat in 297 United States counties. Environmental Epidemiology (Philadelphia,Pa.), 4(3).[5] Dati COVID-19 Italia,
Dipartimento della Protezione Civile , https://github.com/pcm-dpc/COVID-19[6] Dipartimento di Protezione Civile (DPC), [7] EuroMOMO, European MOrtality MOnitoring, [8] del Re, D., & Meridiani, P. (2020). Monitoring the Covid-19 epidemics in Italy from mortality data. medRxiv2020.05.07.20092775.[9] Rosano, A., Bella, A., Gesualdo, F., Acampora, A., Pezzotti, P., Marchetti, S., ... & Rizzo, C. (2019). Investigatingthe impact of influenza on excess mortality in all ages in Italy during recent seasons (2013/14–2016/17 seasons).International Journal of Infectious Diseases, 88, 127-134.[10] James, F., & Roos, M. (1975). MINUIT: a system for function minimization and analysis of the parameter errorsand corrections. Comput. Phys. Commun., 10(CERN-DD-75-20), 343-367.[11] Gompertz, B. (1825). XXIV. On the nature of the function expressive of the law of human mortality, and on anew mode of determining the value of life contingencies. In a letter to Francis Baily, Esq. FRS &c. Philosophicaltransactions of the Royal Society of London, (115), 513-583.[12] [13] Bonifazi, G., Lista, L., Menasce, D., Mezzetto, M., Pedrini, D., Spighi, R., & Zoccoli, A. (2020). A simplifiedestimate of the Effective Reproduction Number R t using its relation with the doubling time and application toItalian COVID-19 data. arXiv preprint arXiv:2012.05194., submitted to European Journal of Physics Plus.[14] Ohnishi, A., Namekawa, Y., & Fukui, T. (2020). Universality in COVID-19 spread in view of the Gompertzfunction. Progress of Theoretical and Experimental Physics, 2020(12), 123J01.[15] Johns Hopkins Coronavirus Resource Center, https://coronavirus.jhu.edu/ [16] CovidStat INFN, https://covid19.infn.it/ [17] Resident Italian population, ISTAT, http://dati.istat.it/Index.aspx?DataSetCode=DCIS POPRES1 [18]