key: cord-0148415-44hle6mh
authors: Bonifazi, Gianluca; Lista, Luca; Menasce, Dario; Mezzetto, Mauro; Pedrini, Daniele; Spighi, Roberto; Zoccoli, Antonio
title: Study on the effects of the restrictive measures for containment of the COVID-19 pandemic on the reproduction number $R_t$ in Italian regions
date: 2021-06-01
journal: nan
DOI: nan
sha: 81f24be69c97a652ae6156b965b9956768434a64
doc_id: 148415
cord_uid: 44hle6mh

Since November 6$^{mathrm{th}}$, 2020, Italian regions have been classified according to four levels, corresponding to specific risk scenarios, for which specific restrictive measures have been foreseen. By analyzing the time evolution of the reproduction number $R_t$, we estimate how much different restrictive measures affect $R_t$, and we quantify the combined effect of the diffusion of virus variants and the beginning of the vaccination campaign upon the $R_t$ trend. We also compute the time delay between implementation of restrictive measures and the resulting effects. Three different models to describe the effects of restrictive measures are discussed and the results are cross-checked with two different algorithms for the computation of $R_t$.

The main purpose of this paper is to quantify the effects on the reproduction number R t of the Non-Pharmaceutical Interventions (NPI) put in place in the Italian regions since November 6 th , 2020.

Italian regions and the autonomous provinces of Trento and Bolzano have been classified into the four aforementioned colored levels: red, orange, yellow and white in decreasing levels of restrictions, corresponding to four risk scenarios, for which specific restrictive measures were foreseen, as reported in the Appendix. The scenarios are described in [1] and risk indicators are defined in the Ministerial Decree of April 30 th , 2020 [2] . The color assignments to the Italian regions along the period November 6 th , 2020 -April 26 th , 2021 are displayed in Figure 1 . We note that the change of NPI in this period was quite frequent. For instance, Lombardia, the most populated Italian region with about 10.0 million inhabitants, changed level 15 times in about six months.

Three possible models to evaluate the effects of the NPI on R t will be compared, and we will address two aspects that affect its computation:

• how can R t be properly synchronized with the daily number of positive swabs, from which it is computed;

• the value of the time delay between the implementation of a new level of NPI and its effects on R t .

We will compare the effects of the different NPI in the nine most populated Italian regions: Lombardia, Lazio, Campania, Sicilia, Veneto, Emilia-Romagna, Piemonte, Puglia, Toscana (in decreasing level of population). All together, these regions contain 80% of the total Italian population. We excluded regions with less than 2 million inhabitants, in order to ensure a sufficiently large number of cases and allow a precise R t computation.

In early January 2021 the vaccination campaign started in Italy, and in the same period novel virus variants spread across the country. These two important factors potentially affect R t in opposite ways and we will compute their combined effect. 2 R t and its synchronization with the new COVID-19 cases

We compute R t from the number of new COVID-19 daily cases reported by the Italian Dipartimento della Protezione Civile [3] . In reference [4] (CovidStat from here on) we developed a simple, yet precise, algorithm to compute R t . We publish daily real time estimates of R t , together with other indicators of the development of the Italian outbreak in [5] .

The official R t used by the Italian Government is computed from the date of onset of symptoms, not from the date of the positive swab as has been done in this paper. We verified in [6] that the two methods converge to essentially identical values of R t (see in particular Fig. 6 of reference [6] ), once the two estimations of R t are synchronized. Furthermore, since the R t used by the Italian Government is computed with the Cori et al. algorithm [7] (Cori from here on), we verified in [4] and [6] that our Covidstat algorithm computes R t in full agreement with it. In any case most of the results of this paper will be cross-checked with the Cori algorithm, by using the EpiEstim package [8] and assuming the same modeling of the generation time distribution as measured in [9] .

We prefer to keep our Covidstat algorithm as the default for the following computations because it's much faster to compute (as discussed in [4] ) and because we have full control of its computation details in all the necessary steps, in particular when we proceed to synchronize R t with the data.

The R t computation with the Covidstat algorithm requires a smoothing of the original data as a preliminary step and, subsequently, an exponential fit to the smoothed data, determining in this way the time-dependent growth rate λ(t) of the number of new daily cases. The smoothing is performed applying the Savitzky-Golay algorithm [10] in a time window of 21 days and with a second-order polynomial 1 . The fit to an exponential curve is performed within moving windows of 14 days.

R t is determined from λ = λ(t) using the equation 1:

where θ and κ are, respectively, the scale and shape parameters of the gamma distribution that models the distribution of the generation time, the time interval between infector-infected pair [4] . Estimates of the parameters θ and κ of the generation time distribution are taken from [9] .

To study the time-dependent effects of NPI upon the trends of R t , it is important to assess the correct synchronization of the latter with the daily number of new cases. In real-time applications, it's customary to assign both the smoothing and the resulting λ of the exponential fit to the last day of the time sequence. In this way, the value of R t can be assigned to the latest day. A more reasonable option for synchronization would be to assign the smoothing and the parameter of the exponential fit to the central day of the sequence.

The way the CovidStat algorithm is defined, provides an absolute mean to determine the optimal synchronization of R t : when the new cases are at a maximum or at a minimum, the growth rate λ becomes zero and, according to Eq.1, R t must therefore be equal to one. This feature has the additional advantage of being independent from the assumed distribution of the generation time.

In order to verify the synchronization of R t with the number of new daily cases, we determine the positions of the three prominent peaks in the period of time considered in this paper. To model them, we perform a fit to the sum of three Gompertz derivatives. Gompertz curves are widely and commonly used to describe the development of pandemics, as discussed for instance in [11] . Each derivative is defined as:

In this parametrization, a is the integral of the curve, T p it's the time when the peak is reached and k G is a growth-rate coefficient, modelling the shape, see also reference [11] .

To illustrate the procedure, Figure 2 .2 shows the distribution of the number of new daily cases in Lombardia. A smoothed distribution and the fit of the sum of Gompertz functions are superimposed on the data. We note the large dispersion of the number of new daily cases to be well beyond their expected Poisson statistical fluctuation. This dispersion makes the smoothing of the data necessary and sometimes produces false small secondary peaks that in turn are reflected in the evolution of R t . Table 1 displays the day in which the fitted curve reached the first peak in the nine Italian regions along with the dates when R t crossed the critical value of 1. The curves in Figure 2 .2 and the data in Table 1 show a good agreement between the dates in which a maximum of the new cases distribution is reached and the dates when R t = 1. However, we had to assign the R t values computed with Cori to the initial day of the 7-days time window under which it is computed, in contrast with the indications of [7] .

By adopting this synchronization (for the Covidstat algorithm the smoothing and the exponential fits are referred to the central day of the considered time window, while for Cori R t it is referred to the first day of the time window), the Covidstat R t can only be computed up to 8 days before the last day of available data, while Cori can be computed up to 7 days before the last day.

Data in Table 1 indicate a residual mismatch of a couple of days on average between the day in which a maximum is reached and R t crosses the value of 1. We decided not to adjust the synchronization for this amount of time, because it could be partially due to a bias induced by the fit to the Gompertz functions; in any case the procedure detailed in section 3.1 guarantees an overall correct synchronization.

Once R t is synchronized with the number of new daily cases, we can proceed to plot the R t trends for the nine Italian regions considered. Data are displayed in Figure 3 , the colors of the dots correspond to the three NPI levels: yellow, orange and red. None of the considered regions has ever been classified as white, so we can't compute the effects of Region Peak of new cases the white NPI level. We observe that all regions have the same descending trend at the beginning of the analyzed time range, because they were all subject, at the time, to the same common NPI. Subsequently, the trends began to differ from region to region, because of the uneven NPI levels and their enforcement date. We have to consider several factors that can influence the R t behaviour beyond what is induced by the three NPI levels. Noticeably: a) the contagion tracking and screening procedures differ from region to region; b) individual regions may have adopted specific additional restrictive provisions to those foreseen nation wide; in some occasions, limited parts of the regions, as single municipalities and provinces, where subject to stricter NPI;

c) some degree of freedom about the organization of schools of any degree along with Universities had been deferred to local level, more information about these points is reported in Appendix. Schools closures and partial openings, therefore, did not follow the NPI levels in the same homogeneous way. Protocols for tracking students and schools staff and the actions to be taken in case of positive cases during the school openings were left to decisions taken at regional level. d) mobility inside each region could significantly differ because of different work activities and different network of mobility connections; e) the vaccination campaign started in January 2021 with a potential impact on R t development;

f) virus variants appeared during the same period of time and became dominant towards its end with a potential impact on R t development.

For all these reasons, we believe that the effects induced by the NPI alone are not sufficient by themselves to fully account for the trends of R t . On the other hand, the purpose of this work is not to introduce a model capable of describing the development of the pandemic in all its aspects; we are rather interested in extracting the effects of the individual NPI. 

We introduce three different models to describe the effects of the NPI on the time evolution of R t . They compute a predicted value for R t ,R t , for all the time periods in which the NPI remained unchanged. At the beginning of each new time period we setR t = R t . The time periods do not begin on the calendar day in which the NPI changed, t, but at the later day t + t D , where t D is the time required for the effects of R t to become apparent. t D is an important parameter of this analysis and will be discussed in paragraph 3.1. The models we take into consideration are:

with i=yellow, orange, red. In this model, the effect of the NPIs consists in a multiplicative factor α on R t . This method has been adopted for instance in [13, 14] . Despite its popularity in literature, this model has severe limitations. The most important one is that the values of R t can only asympotically tend to infinite (when α is greater than one) or to zero (when α is smaller than 1).

To overcome the limitations of Model-a, we introduce a model where R t can asymptotically converge to any finite value R * t,i :R

The parameter N defines the number of steps needed to reach R * t . • Model-c

The data of Fig. 3 do not show discontinuities on the first derivative of R t , R t , in the regional trends. To guarantee the continuity of R t we additionally refine Model-b in the following way:

where the effects of the NPIs on the trend of R t are mediated by the time derivative R t (t + t D ); the weight β is an additional degree of freedom in the fits.

We determine estimates of the parameters α i , R * t , t D , N, β by minimizing a χ 2 defined globally over all considered regions, as described in the next section.

The delay time t D necessary to NPI changes to produce modifications in R t trends, is determined by a scan on the overall χ 2 computed for all the regions assuming Model-c in the period of time November 6 th , 2020 -January 15 th , 2021. The χ 2 is defined as

where t indicates the calendar days, T is the total number of days in the considered period of time, σ Rt is the error on R t . We have chosen Model-c because, as will be discussed in the following, it is the algorithm that best reproduces the data. Moreover, we limited the period of time to January 15 th , 2021 because, as also discussed in the following, it is the most stable period of time before the introduction of the vaccination campaign and the spread of virus variants.

The χ 2 scan reported in Figure 4 , obtained by varying the value of the time delay t D in the range from 0 to 22 days, indicates that an absolute minimum is reached for a delay of t D = 8 days. The χ 2 distribution is quite broad and irregular; on the other hand, the time delay is expected to have a quite broad distribution. Two contributions to the dispersion of t D are the finite width of the generation time distribution, and the collection and processing time of molecular swabs that introduce a dispersion of 7.7 days as we estimated in [6] .

Since we assigned the Covidstat R t to the center of the time interval of the fit (8 days before the last day) and the Cori R t to the start of the time interval (7 days before the last day), the value of t D = 8 days implies that changes on NPI can be detected not earlier than 16 days (Covidstat) and 15 days (Cori) after their enforcement. As a cross-check, we display the trend of R t in Sardegna, which is not one of our selected regions, but it is the only one that achieved a white color assignment on March 1 st , 2021, a decision which resulted in a pronounced rise of R t . The R t trend for Sardegna is shown in Figure 5 and is in good agreement with our evaluation of t D = 8 days. With an analogous χ 2 scan, the number of steps N for Model-b and Model-c is fixed to N = 9.

We estimate the parameters of each model for all the individual italian regions and for their combination (we will call it Italy from here on). As a cross-check, we perform all the calculations by using both the Covidstat and Cori algorithms. The results are displayed in Figure 6 . The values for the Italian sample are reported in Table 2 . the χ 2 is fully satisfactory, but, as already discussed in Section 2.3, these models cannot claim to completely describe the course of the pandemic. It's clear from the comparison of the χ 2 values that Model-a does not provide a fit quality to the data comparable to Model-b and Model-c. The latter provides the best fit to the data and in the following we will consider it as our reference.

In all the italian regions the red set of NPI turns out to provide better containment than the orange set and this latter is better than the yellow set, well beyond the errors of the fits. Absolute numbers in the Italian regions can differ (as discussed in Section 2.3, different regions implemented additional and localised containment measures beyond the NPI considered here), but the average values in Italy are in good agreement with the regional trends.

The results indicate that the yellow NPI level has been, in general, inadequate to guarantee the containment of the pandemic, since its fitted R * t value is bigger than 1 for all cases (with the only exception of Lazio with an estimated R * t = 1.0). The orange NPI level produces trends dangerously close to R * t = 1, while the red NPI level guarantees a good containment of the pandemics. We recall that, as reported in Appendix, the main additional restrictions set by the red level were: further restrictions to mobility, closure of all the stores with few exceptions and suspension of all sport competitions. Parameters computed with the CovidStat and Cori algorithms are in general in good agreement.

We display in Figure 7 the fitted valuesR t computed with the three models superimposed to the R t data in the case of Lombardia. A visual inspection of Figure 7 confirms the outcome of the χ 2 of the fits: Model-c is the model that better follows R t trends. It isn't fully satisfactory, but we have discussed in Section 2.2 that some features of the R t trends are due to fluctuations of the data on positive swabs and we discussed in Section 2.3 that the NPI we are modelling here are not the only factors that can influence the trends of R t . In Figure 8 we display the fits to R t ,R t computed with Model-c for all the regions. The colors of the R t dots in Figure 8 are different from those of 

During the period of time considered in this study, two new important facts happened that could potentially impact on the trends of R t . The vaccination campaign began in Italy early in January 2021. According to the data reported in the COVID-19 Opendata Vaccini repository [15] , the fraction of the Italian population vaccinated by January 15 th , 2021 was 4.3% (0%) with the first (second) dose respectively, while, as April 26 th , these fractions were increased to 20.1% (8.3%).

On the other hand, virus variants (variants of concerns or VOC), specifically lineage B.1.1.7 (English variant), began to spread out in Italy in the same period: according to [16] , 86.7% of the new cases in Italy were due to this variant at March 18 th , 2021, while this percentage was 17.8% at February 4 th , 2021 [17] .

Both VOC and vaccinations have a potentially important impact on the spread of the contagion, therefore we performed a fit with Model-c in two separate period of times: November 6 th ,2021 -January 15 th , 2021 and January 16 th , 2021 -April 26 th , 2021 to measure a possible overall effect of these two new factors in the progress of the pandemics. The results of the fits are reported in Table 3 .

The combined effect of virus variants and the vaccination campaign in the considered period of time, resulted in a worsening of NPIs effectiveness. The yellow NPI level, for instance, worsened its effectiveness in containing the spread of contagion by 15%. After this period of time, the effect of the variant B. the vaccination campaign continues. So, hopefully, if no further variants will influence the R t trends, the effect of the considered NPI will improve.

The goal of this paper is the computation of the effectiveness of the three NPI called yellow, orange and red introduced in the Italian regions since November 6 th , 2020.

We firstly synchronized the R t values with the dates in which the NPI changed in the Italian regions. The synchronization is a necessary step for any procedure that tries to compute changes in R t trends due to changes of NPI at given calendar dates. To this purpose we introduced two original methods, the first of which exploits the fact that our original Covidstat algorithm directly connects R t with the growth rate of an exponential fit.

We identified critical issues in the model often used in literature to quantify the effects of NPI measures. We introduced two methods to overcome them and the fit to the data shows that the two new models clearly outperform the original one. We selected the so called Model-c as a reference for all the results of the paper.

We demonstrated that our model is able to describe with just 4 parameters (3 asymptotic values R * t and a weight parameter β) the development of the pandemics in the 9 Italian regions considered in this paper. The development of the pandemics in the Italian regions is quite differentiated because they underwent different NPI restrictions at different times. The agreement between the model and the data is not complete because several local restrictive measures had been taken and because the poor quality of the positive swabs data introduces important noise in the R t trends. We discussed these features in the paper.

The model provides a coherent picture where in all the Italian regions the red set of NPI results to have smaller R * t values of the orange set and this latter has always R * t smaller than the yellow set, well beyond the statistical errors. We find, that averaged over the nine most populated Italian regions, only the red NPI level can produce a significant reduction of R t , below the threshold value of 1. The orange level keeps R * t around 1, while the yellow level targets R t values above one.

We measured a worsening of the containment effects in the period of time when the variant of concern lineage B.1.1.7 spread out in the country, even if, in the same period, the vaccination campaign started.

All the computations have been performed with two R t algorithms: the Covidstat algorithm that we developed [4] and the widely used standard Cori et al. algorithm [7] . We found a good agreement between the results obtained using the two algorithms.

Prevention and response to COVID-19: evolution of strategy and planning in the transition phase for the autumn-winter season, Ministero della Salute and Istituto Superiore Sanità

Dati COVID-19 Italia

A simplified estimate of the Effective Reproduction Number Rt using its relation with the doubling time and application to Italian COVID-19 data

A study on the possible merits of using symptomatic cases to trace the development of the COVID-19 pandemic

A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics

EpiEstim: Estimate Time Varying Reproduction Numbers from Epidemic Curves

The early phase of the COVID-19 outbreak in Lombardy, Italy

Smoothing and Differentiation of Data by Simplified Least Squares Procedures

The use of Gompertz models in growth analyses, and new Gompertzmodel approach: An addition to the Unified-Richards family

Imperial College COVID-19 Response Team, Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

Variant Of Concern) in Italia: lineage B.1.1.7, P.1, P.2,lineage B.1.351, lineage B.1.525

Prevalenza della variante VOC 202012/01, lineage B.1.1.7 Italia

We acknowledge the effort of ISS of making public the data of the symptomatic cases of COVID-19. The present work has been done in the context of the INFN CovidStat project that produces an analysis of the public Italian COVID-19 data. The results of the analysis are published and updated daily on the website https://covid19.infn.it/. The project has been supported in various ways by a number of people from different INFN Units. 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-Bicocca Unit), 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à Politecnica delle Marche) for his supportive contribution.

A brief description of the rules in place for the four NPI levels denominated as White, Yellow, Orange, Red is reported in table 4.Following the Ministerial Decrees of October 20 th , 2020 [18] and March 13 th , 2021 [19] , nursery and elementary schools had been always open, except the period March, 3 rd , 2021 -April 2 nd , 2021 when they closed. High schools and Universities were suggested to organize distance teaching, with some degree of freedom at local level. On April 6 th , 2021 it has been decided to bring the distance teaching at 50% and on April 26 th , 2021 this percentage become 70% to 100%. middle schools (three years, 12 to 14 years old students) had different implementations in the different regions, the most common was to follow for the first year the rules of elementary schools and for the following two years the rules of high-schools.

• Requirement to wear masks outdoors • Social distancing of 1 meter• Curfew: decision taken at the regional level, Sardegna chose a curfew from 11:30PM to 5:00AM.

All the restrictions of White, plus:• Curfew from 10:00PM to 5:00AM.• Shopping centers closed on holidays.• Museums and exhibitions closed.• Distance education for high schools; in-person teaching for preschools, elementary schools, and middle schools. Universities closed, except for laboratories• Public transportation must travel with max 50% occupancy• Coffee shops and restaurants open until 6:00PM. Take-away allowed until 10:00PM. Home deliveries without restrictions• Suspension of activities of amusement arcades, betting shops.• Pools, gyms, theaters, and movie theaters remain closed. Open sports centers.

All the restrictions of Yellow, plus:• Prohibited movement in and out of regions and municipalities (with some exceptions)• Coffee shops and restaurants closed 7 days a week.Take-away allowed until 10:00PM. Home deliveries without restrictions

All the restrictions of Orange, plus:• It is forbidden any movement, even within its own municipality, except for reasons of work, needs and health.• Closed all stores, except supermarkets, grocery stores, newsstands, tobacconists, pharmacies, laundries, hairdressers.• Suspended all sports competitions except those recognized as being of national interest. Closed sports centers.