Data-driven inference of the reproduction number for COVID-19 before and after interventions for 51 European countries

Summary The reproduction number is broadly considered as a key indicator for the spreading of the COVID-19 pandemic. Its estimated value is a measure of the necessity and, eventually, effectiveness of interventions imposed in various countries. Here we present an online tool for the data-driven inference and quantification of uncertainties for the reproduction number, as well as the time points of interventions for 51 European countries. The study relied on the Bayesian calibration of the SIR model with data from reported daily infections from these countries. The model fitted the data, for most countries, without individual tuning of parameters. We also compared the results of SIR and SEIR models, which give different estimates of the reproduction number, and provided an analytical relationship between the respective numbers. We deployed a Bayesian inference framework with efficient sampling algorithms, to present a publicly available graphical user interface (https://cse-lab.ethz.ch/coronavirus) that allows the user to assess and compare predictions for pairs of European countries. The results quantified the rate of the disease’s spread before and after interventions, and provided a metric for the effectiveness of non-pharmaceutical interventions in different countries. They also indicated how geographic proximity and the times of interventions affected the progression of the epidemic.


Introduction
The forecasting of the evolution of the ongoing coronavirus disease 2019  and the effects of nonpharmaceutical interventions are critical components for decision makers across the world.A broad range of data analysis tools and forecasting models have been deployed since the beginning of 2020 to assess the spread of the disease, as well as the expected number of infections and numbers of deaths [1][2][3].A metric that is often deployed to quantify the progress of the disease is the reproduction number, which is the expected number of secondary cases caused by a single infected individual over the infectious period.The beginning of the epidemics is characterised by the basic reproduction number R 0 , whereas the effective reproduction number R t describes its progression in time.R t exhibits significant complexity [4], but it is broadly considered that values of R t above 1.0 indicate a rapid expansion of the infections.To reduce this number below 1.0, governments have applied various kinds of measures, such as social distancing, travel restrictions, closing of public places, schools and nonessential production.The estimation of R t hinges on the forecasting model and the data used to infer their parameters [5].Moreover, an often overlooked fact is that the method by which inference is made plays an important role in these estimates.A well-established technique for such inference is Bayesian inference [6].However, even though Bayesian inference is a very potent method to estimate uncertainties in model structure and parameters, model parameters may be unidentifiable for the chosen type and amount of data, as well as the choice of priors [7,8].In particular for epidemic models, we believe that inferences for which the sampling algorithm has not converged may indicate a mismatch between the model and the data and need to be reported.
In this work we deployed Bayesian inference to quantify the evolution of R t , as well as the time-points of interventions for 51 European countries.The study relied on the Bayesian calibration of the simple, well-established SIR model [9] extended to account for interventions, with data from reported daily infections.We present an online interface that allows for entry of customised data and comparisons between countries.The parameters of the model included the reproduction number, the day of the first intervention and a reduction factor for the reproduction number.By inferring these parameters from data, we identified when interventions became effective and determined the reproduction number before and after the measures had become effective.This indicates how well the imposed restrictions were able to slow down the spread of the disease.
Very recently a related work [10] inferred the intervention points and the evolution of R t in Germany.In the present work we applied our model to 51 European countries.Our model was able to fit the data for the majority of the countries without individual tuning of the parameters or the inference method for each country.Our results are in excellent agreement with the work previously mentioned for Germany [10].Moreover, our results indicated that capturing accurately the results for one country does not generalise well with the same accuracy to all countries.We visualised the inferred quantities on a map and ranked the countries by the disease's spread rate and the effectiveness of the interventions.We found that the results of our model were consistent with those from related studies in Switzerland [11].The results suggest that at the current stage of the pandemic, countries can be categorised according to the inferred R t .Most countries imposed their restrictions within the first weeks of the epidemic starting in March 2020, and since then about 80% of them have managed to bring their country specific R t below 1.0, which indicates a decaying epidemic.Other countries, such as Armenia, Sweden and Moldova, remain at a stage with larger R t .

Data: daily confirmed cases
We calibrated the well established SIR model using data from daily confirmed cases reported in the open source repository Humanistic GIS Lab, University of Washington [12], whose main data source for European countries is the World Health Organization [13].To avoid the stochastic regime of the system, for each country we considered data only after the total number of confirmed cases exceeded five and the number of cases per million of population exceeded two.Population figures for countries are retrieved from the open repository [14].The analysis in this paper considers data up to 18 May 2020.Official times of interventions that we use for comparison in the section "Inference of interventions for 51 European countries" are taken from reference [15], which aggregates the official beginning of non-pharmaceutical interventions by various governments.

Epidemic modelling using the SIR
The SIR model describes the evolution of susceptible (S), infected (I), and removed (R) population where R t is the reproduction number, γ is the removal rate (including recovery and mortality), I 0 is the initial number of infected individuals and N is the size of the population.To model an intervention or a series of measures taken by the government, we consider R t to be a piece-wise linear function of time split into three phases , (2) The three phases are: (i) uncontrolled disease outbreak before the first intervention takes place (t ≤ t int − ½ δ int , (ii) an adoption phase during which one or multiple measures take place (t int − ½ δ int ≤ t ≤ t int + ½ δ int ), and (iii) the period after all interventions have become effective (t ≥ t int + ½ δ int ).
In the first and last regime we assumed that the reproduction number R t remained constant, and during the adoption phase we assumed a linear transition from R 0 to k int R 0 .The intervention time t int , the reduction factor k int ∈ (0,1) and the duration of transition δ int were inferred from the data.

Bayesian inference and sampling algorithms
We quantified the uncertainty in the extended SIR model, as described above, using data of daily confirmed cases for each country.We denoted this data with Î i corresponding to day t i .From the SIR model, the daily incidence was computed as (3)   where ϑ = (R 0 , γ, k int , t int , δ int ) are the model parameters of the SIR model including an intervention as described above.The initial number of infected individuals was set from the available data so that I 0 = Î 0 .We consider the following generative model for the daily incidence data where NB is the negative binomial parametrised by the mean f(ϑ,t) and the dispersion g(ϑ,t).Notice that using this parametrisation, the variance of an observation Î i is given by f(ϑ,t i ) + f(ϑ,t i ) 2 /g(ϑ,t i ).Here, the dispersion function is constant, i.e., g(ϑ,t) = r, where the parameter vector has been extended to include the dispersion parameter r, ϑ = (R 0 , γ, k int , t int , δ int , r).
We denoted the collection of observations by the vector Î = (Î 1 ,...,Î N ) and the corresponding observation times by t = (t 1 ,...,t N ).For the estimation of the probability of ϑ conditioned on the data, we apply Bayes' theorem, p(Î; t) (5) where p(Î|ϑ;t) is the likelihood function, p(ϑ) is the prior probability distribution and p(Î;t) is the model evidence.
Assuming the observations are independent, the likelihood function is given by (6)   We assumed uninformative priors for the model parameter The lower and upper bounds of the priors are summarised in table 1.We assumed that the removal rate γ was the same in all countries and set it to γ −1 = 5.2 days following [16].The inference of ϑ was performed using Korali [17,18], a high-performance framework for uncertainty quantification and optimisation of computational models.The posterior distribution was estimated with Bayesian Annealed Sequential Importance Sampling (BASIS), a reduced variant of the Transitional Markov Chain Monte Carlo (TM-CMC) algorithm [19], run with 5000 samples and default parameters taken from [20].

Graphical user interface for model predictions
We provide an online interface for real time evaluation of our predictive model.The interface allows for sideby-side comparison of two countries with default input data obtained from [12].Modification of individual data points or substitution of different data sets can be achieved by changing the data in the input form.The accuracy of the Bayesian inference can be improved by increasing the number of samples used by the solver.User requests are forwarded from the web server to a remote workstation where the runs are executed.The back-end code of our model evaluation was implemented in Python and performance-critical parts were implemented in C++.The mean and confidence interval of the computed prediction for total and new daily infected cases were displayed together with input data.The inferred intervention date, t int , was highlighted by the vertical line with corresponding reproduction number before and after the intervention.Figure 1 shows a screenshot of our online interface comparing the predictions for Switzerland and Sweden.Interventions introduced by the Swiss government starting on 17 March 2020 have been more successful in stopping the epidemic than the less strict approach of Sweden.The samples drawn from the posterior distribution are further shown in the columns of the selected countries.

SIR model verification
We verified the presented intervention model using the results of a recent related study [10] that reported the results of a similar SIR model with interventions.The investigators performed Bayesian inference of the model parameters applied to Germany.That model had more parameters to described multiple interventions and considered a periodic modulation of the daily reported cases to account for trends such as under-reporting at weekends.We compared our results to this study in figure 2, where the reported quantities are the effective growth rate λ * (t) = γ(R t (t)−1), the number of new daily cases and the total number of cases.The effective growth rate relates to one limiting case of the SIR model where the total number of infected individuals is small compared to the population such that S/N → 1 and the model reduces to an exponential growth I(t)

Original article
∝ exp(λ * t), which is valid at early stages of the epidemic.Figure 3 shows the corresponding parameters drawn from the posterior distribution.
We want to highlight that our model was simpler and had fewer parameters than the one suggested in [10].Also we observed that the fitted model agreed well with the reported data and provides an estimate for the intervention time in Germany.

Inference of interventions for 51 European countries
We applied our framework to 51 countries of the European continent and infer the model parameters for each country separately.The inferred reproduction number before and after intervention is shown in figures 4 and 5.The map and the scatter plot reveal the differences in the values of R t and effects of interventions among various countries.Certain countries showed strong similarity in the inferred parameters.For example, the neighbouring Switzerland and Austria both have R t dropping from 2.3 to 0.7 after interventions, even though the latter imposed more constraining measures.The same was evident for other pairs such as Italy and France, and also Poland and Ukraine.However, no single explanation for these similarities and differences can be suggested.Physical proximity of the countries did not always lead to similar progression of the epidemic.One example is Poland, where the epidemic started slower but the imposed restrictions proved less efficient than in adjacent Lithuania.The disease's development depended on the types of measures taken by governments, as well as the demographics and the particular features of the social interactions in the population.For instance, the UK and Sweden initially decided against imposing strict measures and therefore underwent a longer epidemic with more infections than their neighbors, Ireland and Norway, which corresponded to larger inferred values of R t .
Another inferred parameter was the intervention time.The countries are ordered by the inferred intervention time, which helps to identify countries with similar inferred parameters.For certain countries, such as Switzerland, Germany and Greece, the inferred intervention time matched the first official announcements with a delay of at most one week.In other cases, the delays reached up to one month.This may indicate that the measures were imposed gradually and did not show an immediate effect.

A comparison of the SIR and SEIR models
A limitation of the SIR model is that its calibration is performed only on confirmed cases, which may be only a fraction of those infected in a certain population.More complex models also account for hospitalised cases and reported deaths.However, such data were not be readily available for all 51 countries examined in this study.Nevertheless, we provide a comparison of the dynamics of the SIR model with the more complex SEIR model that also accounts for an incubation period, during which an infected individual is not yet infectious [23].We used the data of infected individuals in Switzerland and the results indicated that the SIR and the SEIR model have very similar fits for the daily infections and exhibit similar dynamics.At the same time, the SEIR model resulted in a value of R t larger than the one estimated by SIR.This difference can also be quantified by a simple relationship (equation 7) that hinges on the fact that while models have the same definition of R t they estimate different growth rates of I(t) and therefore of N -S(t).Compared with SIR (equation 1), the SEIR model adds a parameter α  Fitting these linearised models to the same data implies equal growth rates, i.e., λ SIR = λ SEIR .In turn this leads to a simple analytical relationship between the respective R t R t SEIR = R t SIR ( 1 + ( R t SIR − 1 ) γ α ) (7)   An example in figure 8 shows that both models can achieve the same growth rates before and after interventions, but for the SEIR model that requires a larger R 0 (2.25 for SIR and 3.82 for SEIR) and a smaller reduction factor k int (0.27 and 0.124, respectively).Likewise, the inference with both models on data for Switzerland in figure 9 gives similar predictions and intervention times, but different estimates of R 0 : 2.25 for SIR and 4.0 for SEIR.For the incubation period in the SEIR, we used α −1 = 2.9 days [24].Results in figure 9 were obtained using the nested sampling algorithm [25] implemented in Korali.

Outlook
The present Bayesian inference framework integrated data from reported daily infection cases with a simple, wellestablished epidemiology model, to provide forecasts for the evolution of the disease with quantified uncertainties.Moreover, it provided data-driven estimates for the time-    Georgia, Ireland, Poland, Russia, Serbia).Understanding the success and failures of the model is a subject of our ongoing investigations.Such differences may lead to a better integration of models and data that may need to account for factors such as demographics, social interactions and geographical proximity between countries.
The present study relied on the SIR model using as the only type of data the reported infections.Models with more parameters fitted on heterogeneous datasets, such as hospitalisations and reported deaths, may provide a more detailed description of the epidemic at the cost of added complexity.At the same time, such data are not readily available for all countries examined in this study.Moreover, increasing the complexity of models does not necessarily lead to more useful predictions [26].
The present framework is modular and can be easily modified to incorporate advanced epidemiology models [11,27], different types of data [28] and advanced sampling methodologies [29].Our ongoing efforts include the handling of heterogeneous data to inform agent-based models for the evolution of COVID-19 as well as studies for optimal testing in order to increase the quantity and veracity of the data regarding the reported cases.
We believe that the present study indicated the need for Bayesian inference, with properly integrated data, models and sampling algorithms, to forecast of the evolution of COVID-19.Forecasts based on models with quantified uncertainties help to better assess the risk involved by the spreading of the disease and guide future interventions.

Disclosure statement
No financial support and no other potential conflict of interest relevant to this article was reported.

Figure 1 :
Figure 1: Screenshot of the online interface comparing predictions for Switzerland and Sweden on 2 July 2020 with 2000 samples https://cse-lab.ethz.ch/coronavirus.

Figure 2 :
Figure 2:The effective growth rate (top), the number of new infections per day (middle) and the total number of infections in Germany (bottom) fitted by present model (blue lines and shades) compared to median predictions from[10] (green crosses).The solid lines show the mean prediction and the shades are the 90% confidence intervals.

Figure 3 :
Figure 3: Inferred parameters (the basic reproduction number R 0 , intervention time t int , reduction factor after intervention k int , duration of the transition δ int and the dispersion parameter r of the negative binomial distribution [4] for Germany: samples drawn from the posterior distribution (upper triangle), marginal distributions of the individual parameter (diagonal) and likelihood heat map (lower triangle).
shows the relation between the reproduction number and the intervention time measured from the beginning of the epidemic in each country.One trend is apparent from the values of R 0 before intervention: countries with larger R 0 tend to introduce the restrictions earlier.The inferred times of the start of interventions t int − ½ δ int are shown in figure7and listed in table 2 together with the times when the restrictions were officially announced.

Figure 4 :
Figure 4: Inferred reproduction number R t before (right) and after interventions (left) shown in color on a map of Europe [21].

Figure 5 :
Figure 5: Inferred reproduction number R t .Mean values (circles) and 90% confidence intervals (gray shades).The countries are ordered by the value of R t after interventions.

Figure 6 :
Figure 6: Inferred mean reproduction number before (R 0 , right) and after intervention (R t , left) vs the time of the start of interventions t int − ½ δ int measured from the beginning of the epidemic in each country.The dots have the same colours as in figure 5 and the country codes are listed in table 2. The red dashed line shows the linear regression of the intervention time with respect to R 0 and R t .

Figure 7 :
Figure 7: Inferred mean time of the start of interventions t int − ½ δ int (circles) with 90% confidence intervals (gray shades) compared with the offcial beginning (vertical bars) of non-pharmaceutical interventions [15] in each country.Empty circles indicate missing data for the offcial time.The countries are ordered by the inferred intervention time such that countries with similar parameters are located closer.

Figure 9 :
Figure 9: Inference with SIR (left) and SEIR (right) on data for Switzerland: predictions (top) and sampling (bottom).Inferred parameters are R 0 , t int , δ int and k int .The incubation period in the SEIR model is 2.9 days.

Table 1 :
Parameters and the prior distributions for the model discussed in the section "Epidemic modelling using the SIR".

Table 2 :
Inferred time of the start of interventions t int − ½ δ int (mean and standard deviation) compared to the official beginning of non-pharmaceutical interventions [15] in each country.
point of interventions across different countries that are in agreement with the first announcement of interventions for a number of European countries.Although the model was able to infer accurately timepoints of interventions for some countries (for example, Austria, Azerbaijan, Denmark, Greece, France, Luxembourg, Norway, Switzerland), it did not capture the correct date for others (for example,