Modeling COVID-19 Incidence by the Renewal Equation after Removal of Administrative Bias and Noise

Abstract

Simple Summary

In the past two years, the COVID-19 incidence curves and reproduction number $R_t$ have been the main metrics used by policy makers and journalists to monitor the spread of this global pandemic. However, these metrics are not always reliable in the short term, because of a combination of delay in detection, administrative delays and random noise. In this article, we present a complete model of COVID-19 incidence, faithfully reconstructing the incidence curve and reproduction number from the renewal equation of the disease and precisely estimating the biases associated with periodic weekly bias, festive day bias and residual noise.

Abstract

The sanitary crisis of the past two years has focused the public’s attention on quantitative indicators of the spread of the COVID-19 pandemic. The daily reproduction number $R_t$, defined by the average number of new infections caused by a single infected individual at time t, is one of the best metrics for estimating the epidemic trend. In this paper, we provide a complete observation model for sampled epidemiological incidence signals obtained through periodic administrative measurements. The model is governed by the classic renewal equation using an empirical reproduction kernel, and subject to two perturbations: a time-varying gain with a weekly period and a white observation noise. We estimate this noise model and its parameters by extending a variational inversion of the model recovering its main driving variable $R_t$. Using $R_t$, a restored incidence curve, corrected of the weekly and festive day bias, can be deduced through the renewal equation. We verify experimentally on many countries that, once the weekly and festive days bias have been corrected, the difference between the incidence curve and its expected value is well approximated by an exponential distributed white noise multiplied by a power of the magnitude of the restored incidence curve.

1. Introduction

The renewal equation, first formulated for birth-death processes in a 1907 note of Alfred Lotka [1], establishes a model for epidemic propagation based on the individual infectiousness. The infectiousness of individuals at time t is characterized by the reproduction number $R_t$, defined as the average number of cases generated by an infected person at time t, and by the generation time [2,3] defined as the probability distribution of the time between infection of a primary case and infections in secondary cases. This probability distribution depends on the incubation time (a permanent biological factor) and on the detection time (which we assume stationary). For these reasons, the distribution of the generation time is supposed to be independent of t. In practice, the generation time is replaced by the observable serial interval ${\Phi }_s$ which represents the time distribution of the delay of the onset of symptoms between primary and secondary cases. In Figure 1, we show the serial interval obtained in [4] using 689 observed pairs of primary and secondary cases.

The case renewal Equation [5,6] is a classic equation linking $R_t$, $\Phi$ and the incidence $i_t$ of new daily cases,

(1)$i_t=\sum _s{i}_{t-s}{R}_{t-s}{\Phi }_s\mathrm{for}t=0,\dots ,t_c,$

where $t_c$ is the current time. This equation does not account for several strong perturbations of $i_t$. Government statistics of the observed incidence curve are indeed affected by changes in testing and polling policies and by weekend reporting delays. These recording delays and subsequent rash corrections result in impulse noise, and in a strong weekly periodic bias observable on the observed incidence curve $i_t^0$. In [3] this bias is corrected by a seven days sliding average and in [7] it is corrected by multiplying $i_t^0$ by a 7-day periodic factor $q_t$. These bias correcting coefficients $q_t$ are learned by a variational method that we describe below. Our first purpose in this note is to resolve the festive day problem. We denote by $\mathbf{F}$ the set of festive days t, at which the $i_t^0$ curve is strongly affected by the reduction in the number of registered cases. This reduction is compensated by an increase in the number of registered cases the following days. No model has been proposed so far to address this problem, which creates strong impulse noise in any estimation of $i_t$ and $R_t$. We tackle this problem by a variational method computing $R_t$, where both $i_t$ and $R_t$ are considered unknown on festive days and in the next few days. To that purpose, we shall denote by ${\mathbf{F}}_{+}$ the union of festive days and the ones following them affected by the festive day (typically 2 or 3 days after the festive day).

Our second purpose is to provide a noise model for the difference ${\widehat{i}}_t-i_t^r$ between the signal ${\widehat{i}}_t$ corrected of the week-end and festive effects, and its restored version $i_t^r$ using the renewal equation, defined by

(2)$i_t^r=\sum _s{\widehat{i}}_{t-s}{R}_{t-s}{\Phi }_s.$

We provide strong experimental evidence that the relation between ${\widehat{i}}_t$ and $i_t^r$, can be empirically modeled by

(3)${\widehat{i}}_t=i_t^r+{\epsilon }_t{\left(i_t^r\right)}^a,$

where $a>0$ and ${\epsilon }_t$ is a white noise.

This leads us to propose a signal processing version of the renewal equation model taking into account noise and bias and justifying a posteriori the variational method. The proposed observation model linking the observed signal $i_t^0$ to the ground truth incidence $i_t$ is

(4)$q_t{i}_t^0=i_t+{\epsilon }_t{\left(i_t\right)}^a\mathrm{for}t\in [0,t_c]\backslash {\mathbf{F}}_{+},$

where $q_t$ is a quasi-periodic gain with period 7, ${\epsilon }_t$ is a white noise. The exponent a can be estimated for each country and varies between 0.6 and 0.9. The exceptional set ${\mathbf{F}}_{+}$ is introduced because festive days provoke perturbations of the observation model (4). Specifically, the 7 days period of $q_t$ is broken for these groups of days.

We shall verify experimentally on 38 countries (and detail the results on USA, France and Germany) that the normalized error ${\epsilon }_t$ is indeed a white noise with a distribution that is well described by an exponential distribution. This a posteriori noise model contradicts the classic a priori stochastic formulation of the renewal equation where the first member $i_t$ of Equation (1) is assumed to be a Poisson variable, and the second member of this equation is interpreted as the expectation of this Poisson variable. Using this Poisson model leads maximum likelihood estimation strategies to compute $R_t$ [3,8,9,10]. As we shall see, the Poisson model is not verified. Indeed, as we mentioned, the empirically observed standard deviation of the noise follows a power law with exponent a significantly larger than 0.5, which is incompatible with the Poisson model.

The proposed observation model (4) of the pandemic’s incidence curve provides a simple framework enabling:

• a computation of the reproduction number $R_t$;

• a correction of the weekend and festive days bias on $i_t$;

• a verification that the difference between the observed incidence curve after bias correction and its expected value using the renewal equation is a white noise, the parameters of which can be estimated.

Paper organization:

In Section 2, we describe an anterior variational method [7] and point out its main three limitations: its weekly bias correction is strongly periodic, which does not work on long periods; the festive days cause strong perturbations in the inversion, finally no residual noise model is proposed. We therefore modify its variational formulation. In Section 3, we present the results of the statistical analysis of the residual noise on many countries. These examples lead to specify the noise model and to validate a posteriori the proposed inversion model. In Section 4, we discuss the a priori noise models proposed in the literature. Finally, in Section 5, we present the conclusions of this work.

Timely estimates of restored versions of $i_t$ and $R_t$ are extremely useful to tame a pandemic. The proposed restoration and inversion algorithm can be run through an online demo [11] for every day in every country and U.S. state. The demo plots the objects of this paper, namely the incidence curve $i_t^0$, its bias corrected version ${\widehat{i}}_t$, its fully restored version $i_t^r$, finally the main pandemic index, the time-dependent reproduction number $R_t$. Figure 2 illustrates the application of the variational method of Section 2 to USA on 1 February 2022, as displayed by the online demo. Figure 3 compares the results of this inversion method, applied with and without festive day bias correction, obtained for France on 6 January 2022.

We can summarize the main contributions of this paper in the following way:

• Based on the case renewal equation, we propose a new variational model which estimate:

  • A time varying reproduction number $R_t$

  • A restored incidence curve with the weekly and festive day biases corrected.

  • The weekly seasonality profile of the incidence curve.

• We verify experimentally, on many countries, that, once the weekly and festive days biases have been corrected, the difference between the incidence curve and its expected value using the renewal equation is well approximated by an exponential distributed white noise multiplied by a power of the magnitude of the restored incidence curve.

2. The Proposed Variational Model

The EpiInvert method proposed in [7] is a deconvolution + denoising procedure to solve the functional Equation (1) using the Tikhonov–Arsenin [12,13] variational approach. EpiInvert estimates both $R_t$ and a restored $i_t$ corrected for the weekend bias. To remove the weekend effect, it computes a 7-day periodic multiplicative factor $q=(q_0,q_1,q_2,q_3,q_4,q_5,q_6)$. From the observed incidence curve and the serial interval, $R_t$ and q are jointly estimated by minimizing

(5)$\begin{array}{cc}\hfill & E(R,q)=\sum _{t=0}^{t_c}{\left(\frac{q_{t\%7}{i}_t^0-{\sum }_s{R}_{t-s}{i}_{t-s}^0{q}_{(t-s)\%7}{\Phi }_s}{media{n}_{(t-\tau ,t]}\left(i^0\right)}\right)}^2+w\sum _{t=1}^{t_c}{(R_t-R_{t-1})}^2\hfill \end{array}$

where $t\%7$ denotes the remainder of the Euclidean division of t by 7 and $media{n}_{(t-\tau ,t]}\left(i^0\right)$ is the median of $i_t^0$ in the interval $(t-\tau ,t]$ used to normalize the energy with respect to the size of $i_t$ (the value of $\tau$ is fixed to 21 (3 weeks) in the experiments). The total number of cases is preserved by adding to (5) the constraint on $q_t$:

(6)$\sum _{t=t_c-T+1}^{t_c}{i}_t^0=\sum _{t=t_c-T+1}^{t_c}{q}_{t\%7}{i}_t^0,$

where T is a period of analysis empirically fixed to $T=56$ days. The minimization of the above energy yields estimates of $R_t$, q and a restored incidence curve.

One limitation of using a 7-day periodic formulation to model the weekend effect is that it does not take into account the variation over time of the seasonal profile. To deal with this issue, we consider $q_t$ for $t=0,\dots ,t_c$ allowing different correction factors $q_t$ for every day but keeping the values $q_t-q_{t-7}$ small which forces $q_t$ to be quasi-periodic. A regularity assumption for the seasonality is commonly used in the study of time series as it is the case of the standard Holt–Winters’ seasonal method [14].

In addition to the weekend bias, festive days can introduce a strong bias in the incidence values. On a festive day $t\in \mathbf{F}$, a sharp decrease in the number of registered incident cases is generally observed. This is compensated by increased incidence numbers in the next few days. Assuming that each festive day, $t\in \mathbf{F}$, mainly affects the value of the incidence curve in the festive day and in the next $M_t$ days (where $M_t$ is an algorithm parameter (by default we fix $M_t=min\{2,t_c-t\}$)), we consider the values of $i_t^0,i_{t+1}^0,\dots ,i_{t+M_t}^0$ as unknown. We denote by ${\mathbf{F}}_{+}$ the union of the festive days $t\in \mathbf{F}$ and the $M_t$ days following them. We set $i_t^f=i_t^0$ for $t\notin {\mathbf{F}}_{+}$ and consider the values ${\left(i_t^f\right)}_{t\in {\mathbf{F}}_{+}}$ as unknowns. Then the new proposed inversion functional is

(7)$\begin{array}{cc}\hfill & E(R,q,{\left(i_t^f\right)}_{t\in {\mathbf{F}}_{+}})=\sum _{t=0}^{t_c}{\left(\frac{q_t{i}_t^f-{\sum }_s{R}_{t-s}{i}_{t-s}^f{q}_{t-s}{\Phi }_s}{media{n}_{(t-\tau ,t]}\left(i^0\right)}\right)}^2+w_R\sum _{t=1}^{t_c}{(R_t-R_{t-1})}^2+\hfill \\ \hfill & \sum _{t\in \mathbf{F}}{\lambda }_t{\left(\frac{{\sum }_{k=0}^{M_t}{i}_{t+k}^f-{\sum }_{k=0}^{M_t}{i}_{t+k}^0}{media{n}_{(t-\tau ,t]}\left(i^0\right)}\right)}^2+w_q\sum _{t=7}^{t_c}{(q_t-q_{t-7})}^2,\hfill \end{array}$

The values $i_t^f$ for $t\in {\mathbf{F}}_{+}$ are set free in the minimization. Yet the third term in the functional ensures that the overall number of cases in the affected days remains unchanged. For each $t\in \mathbf{F}$, ${\lambda }_t\ge 0$ represents the weight we assign to this constraint for each festive day. We fix, experimentally, ${\lambda }_t=2^{t_c-t-2}$ if $t_c>t$ and ${\lambda }_t=0$ if $t_c=t$. In other terms, the value of ${\lambda }_t$ is adjusted according to the number of days that have passed since the festive day. To keep a smooth seasonality we add to the energy a regularization term where we penalize high values of $q_t-q_{t-7}$. The parameters $w_R$ and $w_q$ are regularization weights with default values $w_R=w_q=5$. Their values are proven in [7] to be nearly optimal for COVID-19 incidence curves.

By minimizing this energy we obtain the reproduction number $R_t$, the seasonality $q_t$ and $i_t^f$, which corresponds to the original incidence $i_t^0$ but with the optimized values in the festive days. The bias corrected incidence ${\widehat{i}}_t$ defined in model (3) is given by ${\widehat{i}}_t=q_t{i}_t^f$.

The estimated incidence curve must preserve the number of cases. In the original EpiInvert formulation this constraint is enforced by (6) on its analysis interval $(t_c-T,T]$. In the new formulation, the interval time of analysis is the whole time interval $[0,t_c]$. Extra conditions are required to keep $i_t^0$ close to ${\widehat{i}}_t$ and $i_t^r$. Therefore, to preserve the number of cases we add to the energy (7) the constraints on $q_t$:

(8)$\sum _{t=0}^{t_c}{i}_t^f=\sum _{t=0}^{t_c}{q}_t{i}_t^f;\sum _{t_c-14(k+1)}^{t_c-14k}{i}_t^f=\sum _{t_c-14(k+1)}^{t_c-14k}{q}_t{i}_t^f\mathrm{for}k=0,1,2,\dots .$

The first constraint corresponds to a global preservation of the number of cases in the whole period and the second one corresponds to a local preservation of the number of cases every 2 weeks. In particular, the second constraint ensures a good agreement between the epidemiological indicator given by the accumulated number of cases in the last 14 days of the original incidence curve and the estimated ones using the proposed method. This indicator is currently widely used to evaluate the current epidemic transmission.

The minimization of the energy (7) is obtained by alternating steps computing in turn $R_t$, $q_t$, and then $i_t^f$ (for $t\in {\mathbf{F}}_{+}$) until convergence. The above constraints are added to the minimization by the Lagrange multiplier technique.

3. Results

We used the incidence data published in [15] for France, [16] for Germany, [17] for Spain and [18] for the rest of countries. We checked the observation model and its inversion on the 626 daily incidence data from 24 March 2020 to 9 December 2021 for 38 countries and will detail the results for France, Germany, and the USA. In general, for the festive days we fixed $M_t=2$, so the method estimated the incidence value of the festive day and of the next 2 days. However, not all festive days disturb the incidence in the same way. Parameter $M_t$ allows us to adapt the number of days affected. To illustrate this option we set $M_t=5$ for Thanksgiving in the USA in 2021 because this festive day causes in 2021 a longer perturbation in the number of registered cases. Figure 4, Figure A1 and Figure A2 show the minimization results for the energy (7). They display for each country (i) the original incidence curve $i_t^0$, (ii) the incidence curve after bias correction ${\widehat{i}}_t$, (iii) the restored incidence curve $i_t^r$ using the renewal Equation (1), (iv) the weekly bias correction factors $q_t$, (v) the reproduction number estimation $R_t$ and (vi) the normalized error defined by

(9)${\epsilon }_t=\frac{{\widehat{i}}_t-i_t^r}{{\left(i_t^r\right)}^a}.$

The power a was obtained through log-log linear regression. Indeed, if $|{\widehat{i}}_t-i_t^r|$ is proportional to ${\left(i_t^r\right)}^a$, then $log\left(\right|{\widehat{i}}_t-i_t^r\left|\right)\approx a·log\left(i_t^r\right)+b$, and a and b can be estimated by a linear regression between $log\left(\right|{\widehat{i}}_t-i_t^r\left|\right)$ and $log\left(i_t^r\right)$. Its results are illustrated for 38 countries in Figure 5 and Table A2. The Pearson correlation p-values in this table confirm the linear relation. The estimated exponent a varies between 0.7 and 0.9, and the constant coefficient b varies between −0.11 and −2.6. For the world we have $a=0.76$ and $b=-1.16$.

We performed a control test on a Brownian motion simulated by starting from 10,000 and sampling $i_{t+1}-i_t\simeq \mathcal{N}(0,100)$. The obtained exponent a is negative ($a=-1.01$) and we have $b=13.4$. Both values are far away from the group of coefficients of real incidence curves. The p-value for the control is anyway non significant (0.0844), compared to the extremely small p-values for the real incidence curves. Figure A3 shows the results of the variational inversion method on the Brownian control. For this control, both $R_t$ and the weekly seasonality correction coefficients stay very close to 1 as should be expected, with means 1.001 and 1.00002, and standard deviations 1.7% and 0.3% respectively.

Next, we looked for a stochastic model of the normalized error ${\epsilon }_t$ defined by (9). Figure 4, Figure A1 and Figure A2 visually support a stationarity assumption for ${\epsilon }_t$ in France, Germany and USA.

In Figure 6 we show the autocorrelation function for these three countries. For most non-zero shifts, its value stays inside the 95% confidence interval for the stationarity assumption. (This interval is indicated by horizontal blue lines in the plot.) Similar results were obtained on 33 more countries, as illustrated in Figure A5. These results support a white noise assumption for ${\epsilon }_t$.

We finally estimated the parameters of the distribution of ${\epsilon }_t$ assuming an exponential power distribution with density

(10)$\frac{\beta }{2\alpha \Gamma (1/\beta )}{e}^{-{\left(\frac{|x-\mu |}{\alpha }\right)}^{\beta }},$

where $\mu$ is the location, $\alpha$ the scale and $\beta$ the shape. These parameters to approximate ${\epsilon }_t$ by an exponential power distribution were estimated by the R-package normalp [19].

In Figure 7, we plot for these three countries the histogram of the distribution of ${\epsilon }_t$ and its approximation by a normal ($\beta =2$) and by the obtained optimal exponential distribution. We display the same result for 33 more countries in Figure A4.

Table A1 provides the results for all countries. Columns 5 to 8 in the table provide the parameters of the optimal exponential law: location, scale, shape. In all cases the exponent remains close to 1. Figure 8 displays a quantile-quantile plot comparing ${\epsilon }_t$ with the estimated exponential distribution for three countries: France, Germany, USA. The linear fit is excellent, and this goodness of fit is confirmed for 33 more countries in Figure A6.

4. Discussion of Previous Models

4.1. The Fraser Renewal Equation

In our proposed incidence model, we used the general integral Equation (1), which is a functional equation in $R_t$. Integral equations have been previously used to estimate $R_t$: in [20], the authors estimate $R_t$ as the direct deconvolution of a simplified integral equation where $i_t$ is expressed in terms of $R_t$ and $i_t$ in the past, without using the serial interval. A simpler functional equation than (1) was proposed in Fraser [21] (Equation (9)),

(11)$i_t=R_t\sum _s{i}_{t-s}{\Phi }_s.$

This equation is derived from the general case renewal Equation (1) by assuming that $R_t$ is constant in the serial interval support. It computes the “instantaneous reproduction number” and represents the number of secondary cases arising from an individual showing symptoms at a particular time, assuming that conditions remain identical after that time, in contrast with the case renewal Equation (1). This last equation applied to the incidence curve is coherent if ${\Phi }_s$ denotes the serial interval between two cases, which can have negative dates, because an infectious may be detected after the infection cases she caused. Using (11) requires that ${\Phi }_s$ only has positive dates. This explains why [22] proposed to estimate the generation time, namely the (always positive) time between two infections, before using it in (11). The advantage of Equation (11) is that $R_t$ is estimated at time t from the past incidence values $i_{t-s}$ by a simple division, provided that ${\Phi }_s=0$ for $s<0$:

(12)$R_t=\frac{i_t}{{\sum }_s{i}_{t-s}{\Phi }_s}.$

4.2. Deterministic Implementations Using Fraser’s Renewal Equation and Other Models

Many papers estimating $R_t$ use the deterministic causal renewal Equation (11). This is the case of [23,24,25]. This last paper also involves the Wallinga–Teunis formulation [2], also based on the renewal equation but only allowing a backward estimate of $R_t$ (see the discussion in [7]). Some papers such as [26] propose a simplified version of (11). See also [27], who use this equation but estimate the probability distribution ${\Phi }_s$ by a maximum entropy method. A few papers use another deterministic model, the Wallinga–Teunis formulation, to compute $R_t$ [28], or a SIR model, such as in [29], where the time variable parameter $\beta \left(t\right)$ of the three ODE’s of a SIR model is estimated from incidence data in a seven days sliding window.

4.3. Stochastic Observation Models for $i_t$ and $R_t$

The renewal Equation (11) is often endowed with an a priori stochastic Poisson model as

(13)$i_t=\mathcal{P}\left(R_t\sum _s{i}_{t-s}{\Phi }_s\right).$

In this stochastic formulation, the first member $i_t$ of Equation (11) is assumed to be a Poisson variable, and the second member of this equation is interpreted as the expectation of this Poisson variable. This leads to a maximum likelihood estimation strategy to compute $R_t$ (see [3,8,9,10,30]). This form of the renewal equation is proposed and used in [3] and in the EpiEstim software. It is highly recommended in a recent review [31] signed by representatives from ten different epidemiological labs from several continents. Many papers dedicated to the computation of $R_t$ use this model, for example [32,33,34], who also assume that $R_t$ is a Poisson variable, and [35] who also assume that $R_t$ also is a random variable following a Gamma distribution. In [36], the authors use the stochastic form of the renewal Equation (13) where they call ${\Phi }_s$ causal serial interval. Then $R_t$ is estimated jointly on all regions of a country by a variational model containing a spatial total variation regularization to ensure that $R_t$ is piecewise constant, and the $L^1$ norm of its time Laplacian to ensure time regularity. The functional also penalizes outliers, typically Sundays and holidays by assuming a sparse structure of such events. See also [37] for an exposition of the application of this method.

In [38], the method Epifilter is introduced as an extension of EpiEstim and of the Wallinga–Teunis formulation. Epifilter has been applied in practical studies such as [39]. The core of Epifilter is again the causal renewal equation in Poisson form (13). Yet, the author proposes a doubly stochastic model, as $R_t$ is assumed to follow a recursive discrete Brownian motion of the sort

(14)$R_t=R_{t-1}+\eta \sqrt{R_{t-1}}{ϵ}_{t-1},$

where $ϵ\simeq \mathcal{N}(0,1)$ and $\eta$ is a user parameter, that we can interpret as a regularity control on $R_t$. Then ${\left(R_s\right)}_{s\le t}$ is computed from the incidence data ${\left(i_t\right)}_{s\le t}$ by recursive filtering. The method is complemented by Bayesian (backward) recursive smoothing that brings a better estimate on low incidence periods.

Similarly, in [40], a parametric model with a stochastic multiplicative term is proposed for $R_t$ where the stochastic term is a Gamma law with prescribed standard deviation. The parameters are estimated in several prefectures in interaction to provide the best fit to incidence data linked to $R_t$ through the causal renewal Equation (11).

A few papers assume a negative binomial a priori for the incidence [41]. Nevertheless, the equations given in the paper indicate the adoption of the renewal Equation (11) and put the stochastic process on $R_t$ by assuming $R_t\simeq R_{t-1}GP$ where $GP$ is a squared exponential kernel. The very same model is used in [42], and is based on the authors’ software EpiNow2. Similarly in [43], incidence $i_t$ and reproduction number $R_t$ are linked through the classic SIR model; a parametric piecewise linear model for $R_t$ is estimated by fitting the parameters to real incidence data. Here, the daily incidence data are modeled as a negative binomial, with mean given by the deterministic solution of the SIR equations and unknown dispersion.

In [44], a direct stochastic model is proposed for $R_t$, assuming that its log derivative is Brownian, namely

$d\left(log\right(\beta \left(t\right))=\nu dB(t)$

where $\nu$ is the volatility of the Brownian process to be estimated. Then we have

$R_t=C\beta \left(t\right)s\left(t\right),$

where C is a constant depending on steady transmission characteristics and $s\left(t\right)$ is the proportion of the population that is susceptible. The case incidence is then estimated through an SEIR model. We refer to [45] for a still more complex stochastic model for $R_t$, depending on three stochastic parameters.

5. Conclusions

In [7], we have proven extensively by simulations and experiments on live worldwide COVID-19 incidence data that using the simplified causal renewal Equation (11) incurs in a five days delay in the estimation of $R_t$, compared to the Nishiura renewal Equation (1). This is why we used here this second model.

All of the stochastic models mentioned in Section 4.3 are formulated a priori. To the best of our knowledge, no there has been no a posteriori verification of their noise models on $i_t$ or $R_t$. In contrast, we have proposed to learn the noise model from data and to verify a posteriori that the noise model is correct. Our experiments show that the weekly and festive administrative perturbations are more important than the noise. Hence, they must be corrected first to enable a proper noise analysis.

These experiments seem to confirm the validity of the observation model (4). As we saw, this model can be inverted by minimizing the energy (7). This minimization yields three signals: a restored incidence on the festive days, the administrative bias correcting coefficients $q_t$ that are quasi-periodic with period 7, and the time varying reproduction number $R_t$, arguably the pandemic’s most useful control parameter. Last but not least, the renewal equation deduces a restored incidence $i_t^r$ by (2) from the bias compensated incidence ${\widehat{i}}_t$. The modeling loop was closed by verifying that the normalized error defined by (9) is a white noise. We also found that this noise follows an exponential distribution. This analysis discards the Poisson model for the pandemic’s case count $i_t$. A pure case count should be a Poisson noise, but we saw that the main perturbation was an administrative bias which, once compensated, leaves behind a noise with standard deviation proportional to a power larger than 0.5 of the case count $i_t$. Under the Poisson model this standard deviation would have been equal to the square root of $i_t$.

In summary, based on the renewal equation inversion, this work contributes to a better understanding of the dynamic of the registered administrative observation of the incidence curve, its weekly seasonality, the influence of the festive days and the expected noise model in the observation of the incidence curve.

Author Contributions

Conceptualization, L.A. and J.-M.M.; methodology, L.A., J.-D.M. and J.-M.M.; software, L.A.; validation, L.A., J.-D.M. and J.-M.M.; formal analysis, L.A., J.-D.M. and J.-M.M.; investigation, L.A., J.-D.M. and J.-M.M.; resources, L.A., J.-D.M. and J.-M.M.; data curation, L.A., J.-D.M. and J.-M.M.; writing—original draft preparation, L.A. and J.-M.M.; writing—review and editing, L.A., J.-D.M. and J.-M.M.; visualization, J.-D.M.; supervision, L.A. and J.-M.M.; project administration, L.A. and J.-M.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets of the incidence curve were analyzed in this study. These data can be found in [15] for France, [16] for Germany, [17] for Spain and [18] for the rest of countries.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. The serial interval [Forumla omitted. See PDF.] obtained by [4]. The bars represent the observed number of cases in function of the number of days between the onset of symptoms in primary and secondary cases. The dotted line is its approximation by a scaled and shifted log-normal distribution.

Figure 2. Illustration of the online inversion method [11]. On the left in red, the obtained reproduction number [Forumla omitted. See PDF.] and in black its estimate obtained by the classic EpiEstim method. On the right in green, the original incidence curve [Forumla omitted. See PDF.] of new cases, in blue the incidence curve [Forumla omitted. See PDF.] corrected of the weekend and festive biases, and in red the final reconstructed incidence curve [Forumla omitted. See PDF.] obtained from [Forumla omitted. See PDF.] by the application of the renewal equation. Estimate obtained for USA on 1 February 2022.

Figure 3. Incidence curve (in green) of France up to 6 January 2022. In blue, the incidence corrected of the weekly bias, in cyan the incidence corrected of the weekly and festive day. The Christmas holidays introduce a distortion in the weekly bias corrected incidence that is corrected by the festive day bias correction.

Figure 4. From top to bottom: (i) the original incidence curve [Forumla omitted. See PDF.] of France, (ii) the incidence curve after bias correction [Forumla omitted. See PDF.], (iii) the restored incidence curve using the renewal equation [Forumla omitted. See PDF.], (iv) the weekly bias correction factors [Forumla omitted. See PDF.], (v) the reproduction number estimation [Forumla omitted. See PDF.] and (vi) the normalized error [Forumla omitted. See PDF.], where a is the optimal exponent obtained by regression (see Table A2).

Figure 5. Worldwide log-log correlations between restored incidence [Forumla omitted. See PDF.] and the residual [Forumla omitted. See PDF.] (defined as restored incidence - bias-corrected incidence). The plot presents the log(error) as a function of the log(incidence). The regression parameters were computed through robust linear regression by the R package MASS. (A): Correlation in France, Germany, and USA, with festive day correction. (B): Spread of the values for 38 countries, without festive corrections. (C): Robust linear regression curves for all countries. The linear regression coefficients a and b can be found in Table A2. The worldwide coefficients are [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 6. For France, Germany and USA, autocorrelation of the normalized error [Forumla omitted. See PDF.], using the festive day correction, obtained with the R-software functionalities (acf() function). The orange dotted line provides the 95% confidence interval for non-correlation. Similar plots for the same countries and 33 more countries, without using the festive day correction, are displayed in Figure A5.

Figure 7. For France, Germany and USA, histogram of the normalized error [Forumla omitted. See PDF.], using the festive day correction, its normal approximation (blue line) and its optimal approximation using an exponential distribution (red line) (we use the R-package normalp to approximate [Forumla omitted. See PDF.] by an exponential distribution). See Figure A4 for the results for the same countries and 33 more countries without using the festive day correction.

Figure 8. Quantile-quantile plot with France, Germany and the USA comparing [Forumla omitted. See PDF.], using the festive day correction, with the optimal exponential distribution using the R-package normalp.

Appendix A

Figure A1. From top to bottom: (i) the original incidence curve of Germany [Forumla omitted. See PDF.], (ii) the incidence curve after bias correction [Forumla omitted. See PDF.], (iii) the restored incidence curve using the renewal equation [Forumla omitted. See PDF.], (iv) the weekly bias correction factors [Forumla omitted. See PDF.], (v) the reproduction number estimation [Forumla omitted. See PDF.] and (vi) the normalized error [Forumla omitted. See PDF.], where a is the optimal exponent obtained by regression (see Table A2).

Figure A2. From top to bottom: (i) the original incidence curve [Forumla omitted. See PDF.] of USA, (ii) the incidence curve after bias correction [Forumla omitted. See PDF.], (iii) the restored incidence curve using the renewal equation [Forumla omitted. See PDF.], (iv) the weekly bias correction factors [Forumla omitted. See PDF.], (v) the reproduction number estimation [Forumla omitted. See PDF.] and (vi) the normalized error [Forumla omitted. See PDF.], where a is the optimal exponent obtained by regression (see Table A2).

Figure A3. Control test, from top to bottom: (i) the test incidence curve [Forumla omitted. See PDF.] which is a Brownian motion, (ii) the test curve after bias correction [Forumla omitted. See PDF.], (iii) the restored incidence curve using the renewal equation [Forumla omitted. See PDF.], (iv) the weekly bias correction factors [Forumla omitted. See PDF.], (v) the reproduction number estimation [Forumla omitted. See PDF.] and (vi) the relative error [Forumla omitted. See PDF.]. Both [Forumla omitted. See PDF.] and the weekly seasonality correction coefficients stay very close to 1, with means 1.001 and 1.00002, and standard deviations 1.7% and 0.3% respectively.

Table A1

Table with the mean and standard deviation of ${\epsilon }_t$ and the parameters of the best fit to the exponential distributions for 36 countries. The data of starred countries in the first three rows have undergone the festive bias correction.

Country	Mean	Std	Location	Scale	Shape ($\mathit{\beta }$)
			Exponential	Exponential	Exponential
FRA *	−0.0283	0.8290	−0.0286	0.5394	1.0000
DEU *	−0.0178	0.4785	−0.0135	0.3433	1.0144
USA *	−0.0044	0.2169	−0.0059	0.1537	1.0000
FRA	0.0109	1.0024	−0.0316	0.6026	1.0000
DEU	0.0091	0.5143	0.0050	0.3458	1.0000
USA	0.0032	0.4779	−0.0097	0.3003	1.0000
ARG	0.0025	0.4430	−0.0286	0.3153	1.0000
AUT	0.0419	1.1030	−0.0041	0.9035	1.2701
BEL	0.0413	1.2175	−0.0366	0.8304	1.0000
BRA	−0.0018	0.4825	−0.0368	0.3312	1.0000
CAN	0.0068	1.2720	−0.0252	0.8290	1.0000
CHL	0.0019	0.2960	−0.0082	0.2138	1.0252
COL	−0.0026	0.2006	−0.0107	0.1490	1.0751
CZE	0.0116	0.5671	−0.0415	0.3755	1.0000
DNK	0.0278	1.2446	−0.0298	0.8126	1.0000
GRC	0.0218	1.2764	−0.0410	0.8847	1.0000
HUN	0.0069	0.6600	−0.0267	0.4410	1.0000
IND	0.0419	0.9891	−0.0084	0.6786	1.0000
IDN	−0.0015	0.3374	−0.0140	0.2607	1.1466
IRL	0.0030	1.1778	−0.0748	0.8252	1.0000
ITA	0.0368	1.1441	0.0141	0.7130	1.0000
JPN	0.0243	0.6647	−0.0254	0.4515	1.0000
MEX	−0.0318	1.7329	−0.0955	1.1091	1.0000
NPL	0.0035	0.8994	0.0005	0.5652	1.0000
NLD	0.0437	0.7185	−0.0404	0.4910	1.0000
PHL	−0.0196	2.0401	−0.0930	1.4011	1.0000
POL	−0.0017	0.1911	−0.0043	0.1268	1.0000
ROU	0.0063	0.9465	−0.0011	0.5798	1.0000
RUS	0.0107	0.3383	0.0066	0.2270	1.0000
SRB	0.0675	1.0140	0.0758	0.7932	1.1728
SVK	0.0024	1.3671	−0.0778	0.8194	1.0000
ZAF	0.0139	0.9110	−0.0320	0.7059	1.1497
ESP	0.0637	1.6068	−0.0047	1.0840	1.0000
CHE	0.0528	1.2228	0.0017	0.8667	1.0000
THA	0.0299	1.3738	−0.0312	0.9374	1.0000
TUN	0.0123	1.3033	−0.0845	0.9224	1.0000
UKR	0.0034	0.4117	−0.0215	0.2586	1.0000
ARE	0.0108	0.4192	−0.0127	0.3265	1.1588
GBR	0.0085	0.3304	−0.0171	0.2163	1.0000

Table A2

Coefficients a and b for 38 countries of the log-log linear regression $ax+b$ between restored incidence $i_t^r$ and the residual $|{\widehat{i}}_t-i_t^r|$ as displayed in Figure 5. The Pearson correlation p-values given by the stats R package confirm a linear relation. The exponent a varies between 0.7 and 0.9. Stars * indicate countries with festive correction. The pvalues are slightly better with festive correction than without. The last row shows the results on the control curve, simulated as a Brownian process. Its large p-value discards a linear log-log relation, and the estimated values of a and b also stand far away from the estimated values for real incidence curves.

Country	a	b	p-Value	Country	a	b	p-Value
FRA *	0.8074272	−1.164141	2.01 × 10${}^{-75}$	FRA	0.8136197	−1.1710322	2.76 × 10${}^{-71}$
DEU *	0.8233846	−1.496739	5.99 × 10${}^{-95}$	DEU	0.8235076	−1.5057318	3.01 × 10${}^{-92}$
USA *	0.9076139	−2.264255	3.16 × 10${}^{-42}$	USA	0.8638492	−1.7287377	6.37 × 10${}^{-37}$
ARG	0.8340299	−1.5574878	1.71 × 10${}^{-101}$	AUT	0.6628437	−0.5661912	3.45 × 10${}^{-86}$
BGD	0.9104934	−2.5672893	6.14 × 10${}^{-56}$	BEL	0.7184413	−0.6589731	3.65 × 10${}^{-61}$
BRA	0.8906214	−1.536314	1.03 × 10${}^{-58}$	CAN	0.7240632	−0.6726824	2.96 × 10${}^{-44}$
CHL	0.8349688	−1.9543089	2.64 × 10${}^{-40}$	COL	0.9175985	−2.2638884	3.03 × 10${}^{-112}$
CZE	0.8520268	−1.4708978	2.88 × 10${}^{-133}$	DNK	0.6900743	−0.6284769	2.78 × 10${}^{-68}$
GRC	0.6555842	−0.5683038	2.58 × 10${}^{-102}$	HUN	0.7838904	−1.3618843	4.47 × 10${}^{-142}$
IND	0.7042499	−0.8457334	8.20 × 10${}^{-68}$	IDN	0.8406915	−1.7674138	5.30 × 10${}^{-97}$
IRL	0.7043354	−0.5484242	1.35 × 10${}^{-89}$	ITA	0.6964125	−0.8659193	2.53 × 10${}^{-71}$
JPN	0.7222903	−1.2445353	5.65 × 10${}^{-85}$	MEX	0.725394	−0.4661005	1.76 × 10${}^{-32}$
NPL	0.7548857	−1.0559482	1.42 × 10${}^{-55}$	NLD	0.7494921	−1.1280471	3.35 × 10${}^{-96}$
PHL	0.6715338	−0.1103984	1.90 × 10${}^{-47}$	POL	0.9306078	−2.6041615	3.02 × 10${}^{-133}$
ROU	0.6920366	−1.0282145	4.11 × 10${}^{-77}$	RUS	0.7212814	−2.0048746	4.05 × 10${}^{-26}$
SRB	0.628712	−0.65103	6.76 × 10${}^{-92}$	SVK	0.7381511	−0.7853881	8.53 × 10${}^{-164}$
ZAF	0.7275793	−0.7811203	9.48 × 10${}^{-69}$	ESP	0.6806819	−0.3916179	2.03 × 10${}^{-42}$
CHE	0.6138378	−0.5491828	1.38 × 10${}^{-75}$	THA	0.7110685	−0.4672682	1.63 × 10${}^{-222}$
TUN	0.7539949	−0.503523	3.03 × 10${}^{-163}$	TUR	0.8998264	−2.658924	1.32 × 10${}^{-68}$
UKR	0.8172308	−1.8996555	1.75 × 10${}^{-70}$	ARE	0.7511088	−1.5460453	3.80 × 10${}^{-52}$
GBR	0.8705096	−1.9395546	1.70 × 10${}^{-96}$	World	0.7631129	−1.1389749	0.00
Brownian	−1.0155743	13.3969412	0.0844

Figure A4. Fit of exponential distributions for 36 countries. In red, the best fitting exponential distribution with shape larger or equal to 1. In black, the best fitting normal law.

Figure A5. Autocorrelation of the normalized error [Forumla omitted. See PDF.] using the R-software functionalities (acf() function) for 36 countries. The dotted lines give the 95% confidence interval for non-correlation.

Figure A6. Quantile-quantile plot with 36 countries comparing [Forumla omitted. See PDF.] (without using the festive day correction) with the optimal exponential distribution using the R-package normalp. In the horizontal axis we show the theoretical quantiles and in the vertical axis, the sample quantiles. Note that the exponential distribution shape parameter [Forumla omitted. See PDF.], indicated on the graphs can have values >1.

References

1. Lotka, A.J. Relation between birth rates and death rates. Science; 1907; 26, pp. 21-22. [DOI: https://dx.doi.org/10.1126/science.26.653.21.b] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17754777]

2. Wallinga, J.; Teunis, P. Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures. Am. J. Epidemiol.; 2004; 160, pp. 509-516. [DOI: https://dx.doi.org/10.1093/aje/kwh255] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15353409]

3. Cori, A.; Ferguson, N.M.; Fraser, C.; Cauchemez, S. A new framework and software to estimate time-varying reproduction numbers during epidemics. Am. J. Epidemiol.; 2013; 178, pp. 1505-1512. [DOI: https://dx.doi.org/10.1093/aje/kwt133] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24043437]

4. Ma, S.; Zhang, J.; Zeng, M.; Yun, Q.; Guo, W.; Zheng, Y.; Zhao, S.; Wang, M.H.; Yang, Z. Epidemiological parameters of coronavirus disease 2019: A pooled analysis of publicly reported individual data of 1155 cases from seven countries. medRxiv; 2020; [DOI: https://dx.doi.org/10.1101/2020.03.21.20040329]

5. Nishiura, H. Time variations in the transmissibility of pandemic influenza in Prussia, Germany, from 1918–19. Theor. Biol. Med. Model.; 2007; 4, 20. [DOI: https://dx.doi.org/10.1186/1742-4682-4-20]

6. Nishiura, H.; Chowell, G. The Effective Reproduction Number as a Prelude to Statistical Estimation of Time-Dependent Epidemic Trends. Mathematical and Statistical Estimation Approaches in Epidemiology; Chowell, G.; Hyman, J.M.; Bettencourt, L.M.A.; Castillo-Chavez, C. Springer: Dordrecht, The Netherlands, 2009; pp. 103-121.

7. Alvarez, L.; Colom, M.; Morel, J.D.; Morel, J.M. Computing the daily reproduction number of COVID-19 by inverting the renewal equation using a variational technique. Proc. Natl. Acad. Sci. USA; 2021; 118, e2105112118. [DOI: https://dx.doi.org/10.1073/pnas.2105112118]

8. Thompson, R.; Stockwin, J.; van Gaalen, R.D.; Polonsky, J.; Kamvar, Z.; Demarsh, P.; Dahlqwist, E.; Li, S.; Miguel, E.; Jombart, T. et al. Improved inference of time-varying reproduction numbers during infectious disease outbreaks. Epidemics; 2019; 29, 100356. [DOI: https://dx.doi.org/10.1016/j.epidem.2019.100356]

9. Liu, Q.H.; Ajelli, M.; Aleta, A.; Merler, S.; Moreno, Y.; Vespignani, A. Measurability of the epidemic reproduction number in data-driven contact networks. Proc. Natl. Acad. Sci. USA; 2018; 115, pp. 12680-12685. [DOI: https://dx.doi.org/10.1073/pnas.1811115115]

10. Obadia, T.; Haneef, R.; Boëlle, P.Y. The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks. BMC Med. Informatics Decis. Mak.; 2012; 12, 147. [DOI: https://dx.doi.org/10.1186/1472-6947-12-147]

11. Alvarez, L.; Colom, M.; Morel, J.D.; Morel, J.M. EpiInvert Online Interface, IPOL: Image Processing On Line. Available online: http://www.ctim.es/epiinvert (accessed on 9 December 2021).

12. Tikhonov, A.N.; Arsenin, V.Y. Solutions of Ill-Posed Problems; Wiley: New York, NY, USA, 1977.

13. Benning, M.; Burger, M. Modern regularization methods for inverse problems. Acta Numer.; 2018; 27, pp. 1-111. [DOI: https://dx.doi.org/10.1017/S0962492918000016]

14. Hyndman, R.; Athanasopoulos, G. Forecasting: Principles and Practice; 2nd ed. OTexts: Melbourne, Australia, 2018; ISBN 1886529043 Available online: OTexts.com/fpp2 (accessed on 9 December 2021).

15. Government of France. Informations COVID-19, Carte et Données. Available online: https://www.gouvernement.fr/info-coronavirus/carte-et-donnees (accessed on 9 December 2021).

16. Robert Koch-Institut. COVID-19-Dashboard. Available online: https://experience.arcgis.com/experience/478220a4c454480e823b17327b2bf1d4 (accessed on 9 December 2021).

17. Spanish Goverment. Situación actual COVID-19. Available online: https://www.sanidad.gob.es/en/profesionales/saludPublica/ccayes/alertasActual/nCov/situacionActual.htm (accessed on 9 December 2021).

18. Ritchie, H.; Ritchie, H.; Mathieu, E.; Rodés-Guirao, L.; Appel, C.; Giattino, C.; Ortiz-Ospina, E.; Hasell, J.; Macdonald, B.; Beltekian, D. et al. Coronavirus Pandemic (COVID-19), OurWorldInData.org. Available online: https://ourworldindata.org/coronavirus-source-data (accessed on 9 December 2021).

19. Mineo, A.M. On the estimation of the structure parameter of a normal distribution of order p. Statistica; 2003; 63, pp. 109-122. [DOI: https://dx.doi.org/10.6092/issn.1973-2201/342]

20. Demongeot, J.; Oshinubi, K.; Seligmann, H.; Thuderoz, F. Estimation of Daily Reproduction rates in COVID-19 Outbreak. medRxiv; 2021; [DOI: https://dx.doi.org/10.1101/2020.12.30.20249010]

21. Fraser, C. Estimating Individual and Household Reproduction Numbers in an Emerging Epidemic. PLoS ONE; 2007; 2, e758. [DOI: https://dx.doi.org/10.1371/journal.pone.0000758]

22. Knight, J.; Mishra, S. Estimating effective reproduction number using generation time versus serial interval, with application to COVID-19 in the Greater Toronto Area, Canada. Infect. Dis. Model.; 2020; 5, pp. 889-896. [DOI: https://dx.doi.org/10.1016/j.idm.2020.10.009]

23. Bonifazi, G.; Lista, L.; Menasce, D.; Mezzetto, M.; Pedrini, D.; Spighi, R.; Zoccoli, A. A simplified estimate of the effective reproduction number R_t R t using its relation with the doubling time and application to Italian COVID-19 data. Eur. Phys. J. Plus; 2021; 136, pp. 1-14. [DOI: https://dx.doi.org/10.1140/epjp/s13360-021-01339-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33868891]

24. Flaxman, S.; Mishra, S.; Gandy, A.; Unwin, H.J.T.; Coupland, H.; Mellan, T.A.; Zhu, H.; Berah, T.; Eaton, J.W.; Guzman, P.N.P. et al. Estimating the Number of Infections and the Impact of Nonpharmaceutical Interventions on COVID-19 in 11 European Countries. Imperial College COVID-19 Response Team. Available online: https://www.imperial.ac.uk/media/imperial-college/medicine/sph/ide/gida-fellowships/Imperial-College-COVID19-Europe-estimates-and-NPI-impact-30-03-2020.pdf (accessed on 30 September 2020).

25. Koyama, S.; Horie, T.; Shinomoto, S. Estimating the time-varying reproduction number of COVID-19 with a state-space method. PLoS Comput. Biol.; 2021; 17, e1008679. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1008679] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33513137]

26. Drewes, H.; Flaeschner, G.; Moeller, P. Improving the reproduction number calculation by treating for daily variations of SARS-CoV-2 cases. medRxiv; 2021; [DOI: https://dx.doi.org/10.1101/2021.08.15.21262071]

27. Tao, Y. Maximum entropy method for estimating the reproduction number: An investigation for COVID-19 in China and the United States. Phys. Rev. E; 2020; 102, 032136. [DOI: https://dx.doi.org/10.1103/PhysRevE.102.032136]

28. Wang, K.; Zhao, S.; Li, H.; Song, Y.; Wang, L.; Wang, M.H.; Peng, Z.; Li, H.; He, D. Real-time estimation of the reproduction number of the novel coronavirus disease (COVID-19) in China in 2020 based on incidence data. Ann. Transl. Med.; 2020; 8, 689. [DOI: https://dx.doi.org/10.21037/atm-20-1944]

29. Shapiro, M.B.; Karim, F.; Muscioni, G.; Augustine, A.S. Adaptive Susceptible-Infectious-Removed Model for Continuous Estimation of the COVID-19 Infection Rate and Reproduction Number in the United States: Modeling Study. J. Med. Int. Res.; 2021; 23, e24389. [DOI: https://dx.doi.org/10.2196/24389]

30. Boulmezaoud, T.Z.; Alvarez, L.; Colom, M.; Morel, J.M. A Daily Measure of the SARS-CoV-2 Effective Reproduction Number for all Countries. Image Process. Line; 2020; 10, pp. 191-210. [DOI: https://dx.doi.org/10.5201/ipol.2020.304]

31. Gostic, K.M.; McGough, L.; Baskerville, E.B.; Abbott, S.; Joshi, K.; Tedijanto, C.; Kahn, R.; Niehus, R.; Hay, J.A.; De Salazar, P.M. et al. Practical considerations for measuring the effective reproductive number, Rt. PLoS Comput. Biol.; 2020; 16, e1008409. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1008409] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33301457]

32. You, C.; Deng, Y.; Hu, W.; Sun, J.; Lin, Q.; Zhou, F.; Pang, C.H.; Zhang, Y.; Chen, Z.; Zhou, X.H. Estimation of the time-varying reproduction number of COVID-19 outbreak in China. Int. J. Hyg. Environ. Health; 2020; 228, 113555. [DOI: https://dx.doi.org/10.1016/j.ijheh.2020.113555] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32460229]

33. Chintalapudi, N.; Battineni, G.; Sagaro, G.G.; Amenta, F. COVID-19 outbreak reproduction number estimations and forecasting in Marche, Italy. Int. J. Infect. Dis.; 2020; 96, pp. 327-333. [DOI: https://dx.doi.org/10.1016/j.ijid.2020.05.029]

34. Hong, H.G.; Li, Y. Estimation of time-varying reproduction numbers underlying epidemiological processes: A new statistical tool for the COVID-19 pandemic. PLoS ONE; 2020; 15, e0236464. [DOI: https://dx.doi.org/10.1371/journal.pone.0236464]

35. Salas, J. Improving the estimation of the COVID-19 effective reproduction number using nowcasting. Stat. Methods Med. Res.; 2021; 30, pp. 2075-2084. [DOI: https://dx.doi.org/10.1177/09622802211008939]

36. Pascal, B.; Abry, P.; Pustelnik, N.; Roux, S.G.; Gribonval, R.; Flandrin, P. Nonsmooth convex optimization to estimate the COVID-19 reproduction number space-time evolution with robustness against low quality data. arXiv; 2021; arXiv: 2109.09595

37. Abry, P.; Pustelnik, N.; Roux, S.; Jensen, P.; Flandrin, P.; Gribonval, R.; Lucas, C.G.; Guichard, É.; Borgnat, P.; Garnier, N. Spatial and temporal regularization to estimate COVID-19 reproduction number R (t): Promoting piecewise smoothness via convex optimization. PLoS ONE; 2020; 15, e0237901. [DOI: https://dx.doi.org/10.1371/journal.pone.0237901]

38. Parag, K.V. Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves. PLoS Comput. Biol.; 2021; 17, e1009347. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1009347]

39. Mee, P.; Alexander, N.; Mayaud, P.; Gonzalez, F.d.J.C.; Abbott, S.; de Souza Santos, A.A.; Acosta, A.L.; Parag, K.V.; Pereira, R.H.; Prete, C.A., Jr. et al. Tracking the emergence of disparities in the subnational spread of COVID-19 in Brazil using an online application for real-time data visualisation: A longitudinal analysis. Lancet Reg.-Health-Am.; 2021; 5, [DOI: https://dx.doi.org/10.1016/j.lana.2021.100119]

40. Jung, S.m.; Endo, A.; Akhmetzhanov, A.R.; Nishiura, H. Predicting the effective reproduction number of COVID-19: Inference using human mobility, temperature, and risk awareness. Int. J. Infect. Dis.; 2021; 113, pp. 47-54. [DOI: https://dx.doi.org/10.1016/j.ijid.2021.10.007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34628020]

41. Abbott, S.; Hellewell, J.; Thompson, R.N.; Sherratt, K.; Gibbs, H.P.; Bosse, N.I.; Munday, J.D.; Meakin, S.; Doughty, E.L.; Chun, J.Y. et al. Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts. Wellcome Open Res.; 2020; 5, 112. [DOI: https://dx.doi.org/10.12688/wellcomeopenres.16006.1]

42. Sherratt, K.; Abbott, S.; Meakin, S.R.; Hellewell, J.; Munday, J.D.; Bosse, N. CMMID COVID-19 Working Group Jit, M.; Funk, S. Exploring surveillance data biases when estimating the reproduction number: With insights into subpopulation transmission of COVID-19 in England. Philos. Trans. R. Soc. B; 2021; 376, 20200283. [DOI: https://dx.doi.org/10.1098/rstb.2020.0283] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34053260]

43. Karnakov, P.; Arampatzis, G.; Kičić, I.; Wermelinger, F.; Wälchli, D.; Papadimitriou, C.; Koumoutsakos, P. Data-driven inference of the reproduction number for COVID-19 before and after interventions for 51 European countries. Swiss Med. Wkly.; 2020; 150, w20313. [DOI: https://dx.doi.org/10.4414/smw.2020.20313] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32677705]

44. Cazelles, B.; Champagne, C.; Nguyen-Van-Yen, B.; Comiskey, C.; Vergu, E.; Roche, B. A mechanistic and data-driven reconstruction of the time-varying reproduction number: Application to the COVID-19 epidemic. PLoS Comput. Biol.; 2021; 17, e1009211. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1009211]

45. Mellan, T.A.; Hoeltgebaum, H.H.; Mishra, S.; Whittaker, C.; Schnekenberg, R.P.; Gandy, A.; Unwin, H.J.T.; Vollmer, M.A.; Coupland, H.; Hawryluk, I. et al. Subnational analysis of the COVID-19 epidemic in Brazil. MedRxiv; 2020; [DOI: https://dx.doi.org/10.1101/2020.05.09.20096701]