key: cord-0739810-bb2xccl3
authors: Gressani, O.; Wallinga, J.; Althaus, C.; Hens, N.; Faes, C.
title: EpiLPS: a fast and flexible Bayesian tool for near real-time estimation of the time-varying reproduction number
date: 2021-12-05
journal: nan
DOI: 10.1101/2021.12.02.21267189
sha: e76b851ce6baee40ea7c15a6266310841d17b7dc
doc_id: 739810
cord_uid: bb2xccl3

In infectious disease epidemiology, the instantaneous reproduction number R(t) is a time-varying metric defined as the average number of secondary infections generated by individuals who are infectious at time t. It is therefore a crucial epidemiological parameter that assists public health decision makers in the management of an epidemic. We present a new Bayesian tool for robust estimation of the time-varying reproduction number. The proposed methodology smooths the epidemic curve and allows to obtain (approximate) point estimates and credible envelopes of R(t) by employing the renewal equation, using Bayesian P-splines coupled with Laplace approximations of the conditional posterior of the spline vector. Two alternative approaches for inference are presented: (1) an approach based on a maximum a posteriori argument for the model hyperparameters, delivering estimates of R(t) in only a few seconds; and (2) an approach based on a MCMC scheme with underlying Langevin dynamics for efficient sampling of the posterior target distribution. Case counts per unit of time are assumed to follow a Negative Binomial distribution to account for potential excess variability in the data that would not be captured by a classic Poisson model. Furthermore, after smoothing the epidemic curve, a ''plug-in'' estimate of the reproduction number can be obtained from the renewal equation yielding a closed form expression of R(t) as a function of the spline parameters. The approach is extremely fast and free of arbitrary smoothing assumptions. EpiLPS is applied on data of SARS-CoV-1 in Hong-Kong (2003), influenza A H1N1 (2009) in the USA and current SARS-CoV-2 pandemic (2020-2021) for Belgium, Portugal, Denmark and France.

health crisis. The reproduction number is also a good proxy for measuring the real-time growth phase of an epidemic and as such, constitutes a key signal about the severity of the outbreak and the required control effort. For this reason, having a robust, accurate and timely estimator of R(t) is a crucial matter that has attracted considerable interest in developing new statistical approaches during the last two decades as summarized in White et al. (2021) . The paper of Gostic et al. (2020) compares several methods for estimating R(t) and

gives clear insights about the main challenges and obstacles that have to be faced. They recommend the method of Cori et al. (2013) and its associated EpiEstim package (Cori, 2021) as an appropriate and accurate tool for near real-time estimation of the instantaneous reproduction number. Another recent approach is proposed in Parag (2021) , where a recursive Bayesian smoother based on Kalman filtering is used to derive a robust estimate of R(t) in periods of low incidence. The EpiNow2 package ) also provides interesting extensions and implementations of current best practices for precise estimation and forecast of the reproduction number using a Bayesian latent variable framework. Spline based approaches have shown to be an interesting tool for flexible modeling of the reproduction number. Azmon et al. (2014) use penalized radial splines for estimating R(t) under a Bayesian setting with misreported data and accelerated the computational implementation by replacing the MCMC scheme with Laplace approximations.

From a frequentist perspective, Pircalabelu (2021) uses truncated polynomials and radial basis splines to model the series of new infections and a derivative thereof as a candidate estimator for the reproduction number.

In this article, we propose a new Bayesian approach termed "EpiLPS" for estimating R(t) based on case incidence data and the serial interval distribution (the time elapsed between the onset of symptoms in an infector and the onset of symptoms in an infected case). Our estimator of R(t) is based on epidemic renewal-equations (Fraser, 2007; Wallinga and Lipsitch, 2007) and Laplacian-P-splines (LPS) smoothing of the mean number of new cases by day of reporting. Time series of new cases by day of reporting are assumed to follow a Negative Binomial distribution to account for potential excess variability in the 3 data that would not be captured by a classic Poisson model. Algorithms related to Laplace approximations and evaluations of B-spline bases are coded in C++ and embedded in the R language through the Rcpp package (Eddelbuettel et al., 2011) , making computational speed another key strength of EpiLPS as R(t) can be estimated in seconds. In addition, EpiLPS can also be used to obtain a smoothed estimate of the epidemic curve that can be of potential interest to further visualize an epidemic outbreak.

The proposed Bayesian methodology is based on a latent Gaussian model for the B-spline amplitudes and opens up two possible paths for inference. The first is called LPSMAP, a fully "sampling free" approach based on Laplace approximations to the conditional posterior of B-spline coefficients. The hyperparameter vector is fixed at its maximum a posteriori and credible envelopes of R(t) are computed via the "delta" method. The second path is called LPSMALA and is an MCMC-based approach with Langevin diffusions for efficient exploration of the posterior distribution of latent variables. The latter approach is computationally heavier than LPSMAP but has the merit of taking into account the uncertainty surrounding the hyperparameters. The underlying Metropolis-within-Gibbs structure keeps the practical implementation to a fairly simple level and the computational cost is reasonable even for long chain lengths.

Compared to existing methods, EpiLPS resembles EpiEstim from a methodological point of view in the sense that R(t) is estimated from incidence time series and a serial interval distribution, yet the two approaches fundamentally differ in many aspects. First, the methodology of Cori et al. (2013) assumes that incidence at time t is Poisson distributed, while EpiLPS assumes a Negative Binomial model. Second, as our approach uses penalized spline based approximations, the prior specifications are imposed on the roughness penalty parameter and not directly on R(t) as in EpiEstim. Third and most importantly, EpiLPS is free of any sliding window specification. An R package for EpiLPS has been developed and is available at https://github.com/oswaldogressani. The software also allows to compute the Cori et al. (2013) estimate of R(t) for the sake of comparison.

The manuscript is organized as follows. Section 2 aims at presenting the Laplacian-P-4 splines model for smoothing count data. We show how the Laplace approximation applies to the conditional posterior of the B-spline amplitudes and also derive the (approximate) posterior of the hyperparameter vector to be optimized. This yields the maximum a posteriori estimate of the spline vector via Laplacian-P-splines (LPSMAP). Section 3 uses LPSMAP and proposes a "plug-in" estimate of R(t) based on renewal equations. Approximate credible intervals for R(t) will also be shown in this section. Section 4 shows an alternative path for estimation of R(t) based on Markov chain Monte Carlo (MCMC). The latter approach uses Langevin dynamics for efficient sampling of the target posterior distribution.

This approach is termed LPSMALA for "Laplacian-P-splines with a Metropolis-adjusted Langevin algorithm". Section 5 is devoted to assess the performance of EpiLPS in various simulation scenarios and make comparisons with EpiEstim. In Section 6, we apply EpiLPS to real world epidemic outbreaks. Finally, Section 7 concludes with a discussion.

2 Methodology behind EpiLPS 2.1 Negative Binomial model for case incidence data Let {y t , t = 1, . . . , T } be a time series of counts during an epidemic of T days with y t ∈ N (set of non-negative integers) denoting the total number of new contaminations on day t.

We assume that the number of cases on day t follows a Negative Binomial distribution y t ∼ NegBin(µ(t), ρ), with µ(t), ρ ∈ R * + := {x ∈ R|x > 0} and probability mass function (see e.g. Anscombe, 1950; Piegorsch, 1990) :

where Γ(·) is the gamma function. Under the above parameterization, we have E(y t ) = µ(t)

and V(y t ) = µ(t) + µ(t) 2 /ρ, so that 1/ρ is the parameter responsible for overdispersion (variance larger than the mean) that is absent in a Poisson setting. In the limiting case lim ρ→+∞ V(y t ) = µ(t) = E(y t ) and we recover the Poisson model. We assume that µ(t)

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 evolves smoothly over the time course of the epidemic and model it with cubic B-splines (Eilers and Marx, 1996) :

where θ = (θ 1 , . . . , θ K ) ⊤ is the vector of B-spline amplitudes to be estimated and b(·) = (b 1 (·), . . . , b K (·)) ⊤ is a cubic B-spline basis defined on the domain T = [r l , T ], where r l is a lower bound on the time axis, typically the first day of the epidemic (i.e. r l = 1). The philosophy behind P-splines consists in specifying a "large" number K of basis functions together with a discrete roughness penalty λθ ⊤ P θ as a counterforce to the induced flexibility of the fit. The parameter λ > 0 acts as a tuning parameter calibrating the "degree" of smoothness and P = D ⊤ r D r + εI K is a penalty matrix built from rth order difference matrices D r of dimension (K − r) × K perturbed by an ε-multiple (here ε = 10 −6 ) of the K-dimensional identity matrix I K to ensure full rankedness. The reader is redirected to Eilers and Marx (2021) for a complete textbook treatment of P-splines. Following Lang and Brezger (2004) , we impose a Gaussian prior on the spline vector θ|λ ∼ N dim(θ) (0, Q −1 λ ), with precision matrix Q λ = λP . For full Bayesian inference, the following priors are imposed on the model hyperparameters. Following Jullion and Lambert (2007) , a robust Gamma prior is specified for the roughness penalty parameter λ|δ ∼ G (ϕ/2, (ϕδ)/2), where G(a, b)

is a Gamma distribution with mean a/b and variance a/b 2 , ϕ = 2 and δ is an additional dispersion parameter with hyperprior δ ∼ G(a δ = 10, b δ = 10). Finally, the following uninformative prior is imposed on the overdispersion parameter ρ ∼ G(a ρ = 0.0001, b ρ = 0.0001).

Let η := (λ, ρ) ⊤ denote the vector of hyperparameters. The full Bayesian model is thus:

6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint λ|δ ∼ G (ϕ/2, (ϕδ)/2) ,

The log-likelihood for the Negative Binomial model is given by:

with g(y t , ρ) = log Γ(y t + ρ) − log Γ(ρ) and= denoting equality up to an additive constant.

The gradient of the log-likelihood with respect to the spline coefficients is:

where:

The Hessian of the log-likelihood with respect to the B-spline amplitudes is:

7 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021.

Using Bayes' rule, the conditional posterior of θ for a given η is:

The gradient and Hessian of the log-likelihood (3) can be used to compute the gradient and

Hessian of the (log-)conditional posterior (4), namely:

The above two equations will be used in a Newton-Raphson algorithm to obtain the Laplace approximation to the conditional posterior of θ:

where θ * (η) and Σ * (η) is the mode and variance-covariance respectively after convergence of the Newton-Raphson algorithm. The latter two quantities are functions of the hyperparameter vector η. An intuitive choice for η is to fix it at its maximum a posteriori (MAP). This is the option retained here, although it is also possible to work with a grid-based approach (Rue et al., 2009; Gressani and Lambert, 2021 ).

The hyperparameter vector η = (λ, ρ) ⊤ will be calibrated by posterior optimization. Following Tierney and Kadane (1986) and Rue et al. (2009) , the hyperparameter vector can be approximated as follows:

8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ;

Approximation (6) can be written extensively as follows:

where the K/2 power of λ comes from the determinant |Q −1 λ | −1/2 = |λP | 1/2 ∝ λ K/2 . As δ ϕ 2 +a δ −1 exp −δ ϕλ 2 + b δ is the kernel of a Gamma distribution for the dispersion parameter δ, the following integral can be analytically solved:

Using the transformation of variables (ensuring numerical stability during optimization) w = log(ρ), v = log(λ), one can show that p(η|D) can be written as follows after using the multivariate transformation method:

The approximated log-posterior becomes:

9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint Equation (7) is numerically optimized and yields η * = argmax η log p( η|D). Plugging the latter vector into the Laplace approximation (5), one obtains the estimate θ = θ * ( η * ) of the spline vector. The latter can be seen as a maximum a posteriori (MAP) estimate of θ.

Thus, the approximated posterior of the spline vector is:

and can be used to construct credible intervals for functions that indirectly depend on θ, such as R(t) as shown in the following section.

In this section, we show how the Negative Binomial model for smoothing incidence counts can be used to estimate R(t) through the renewal equation. Let φ = {φ 1 , . . . , φ k } be a known k-dimensional vector representing the serial interval distribution. The renewal equation describes the link between the number of new cases at day t and past infected cases (up to time point t − 1) weighted by the serial interval distribution. Mathematically:

Rearranging (9) and taking the length k of the serial interval into account, we obtain an equation with the instantaneous reproduction number on the left-hand side:

10 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint A "plug-in" estimator of R(t) at any t ∈ T is obtained by replacing the number of new cases y t by the estimated mean number of cases µ(t) = exp θ ⊤ b(t) , yielding:

Using the indicator function I(·), i.e. I(A) = 1 if condition A is true and I(A) = 0 otherwise, the above "plug-in" estimate of R(t) can be written in a single line:

Replacing the number of new cases y t by the (theoretical) mean number of cases µ(t) = exp(θ ⊤ b(t)) in (10) and using the compact notation, one has:

Let us denote the log of the instantaneous reproduction number in (13) as:

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint

Note that h(θ|t) is seen here as a function of the spline vector θ for a given time point

method. Consider a first-order Taylor expansion of h(θ|t) around θ * ( η * ) (henceforth θ * for the sake of a light notation), the mean of the Laplace approximated posterior of the spline vector in (8):

where the kth entry of the gradient vector ∇h(θ|t) = (∂h(θ|t)/∂θ 1 , . . . , ∂h(θ|t)/∂θ K ) ⊤ is:

It follows that for k = 1, . . . , K we have:

12 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10. 1101 The Taylor expansion in (14) is a linear combination of the vector θ that is a posterior (approximately) Gaussian due to the Laplace approximation. As the family of Gaussian distributions is closed under linear combinations, it follows that h(θ|t) (and hence log R(t))

is a posteriori also (approximately) Gaussian with mean E(h(θ|t)) ≈ h(θ * |t) and variance- (8). This suggests to write:

Thus, from (15) a quantile-based (1 − α) × 100% approximate credible interval for R(t) is:

where z α/2 is the α/2-upper quantile of a standard normal variate. (Chib and Greenberg, 1996) . One of the most popular MCMC methods together with the Gibbs sampler (Geman and Geman, 1984) is the Metropolis-Hastings (MH) algorithm originally proposed by Metropolis et al. (1953) and later generalized by Hastings (1970) . In this section, we propose to implement a modified version of the Metropolis-adjusted Langevin algorithm (MALA) (Roberts and Tweedie, 1996) within the EpiLPS framework. The major advantage of MALA as compared to MH algorithms is that the proposal distribution is based upon a discretized approximation of the Langevin diffusion that uses the gradient of the target posterior distribution. These 13 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint "smarter" proposals make use of additional information about the target density so that algorithms based on Langevin dynamics can converge at sub-geometric rates and tend to be more efficient than naive random-walk Metropolis algorithms Rosenthal, 1998, 2001) . Let ζ = (θ ⊤ , ρ) ⊤ be the (K + 1)-dimensional vector gathering the B-spline coefficients θ and the overdispersion parameter ρ. Using Bayes' theorem, the joint posterior distribution for ζ, λ and δ is:

The analytical formulas of the chosen priors are:

14 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint

Injecting the above priors into (17) yields:

. (18) 4.1.2 Condtional posteriors of ζ, λ and δ

The following conditional posterior distributions can be directly obtained from (18):

(λ|ζ, δ, D) ∼ G 0.5(K + ϕ), 0.5(θ ⊤ P θ + δϕ) ,

(δ|ζ, λ, D) ∼ G (0.5ϕ + a δ , 0.5ϕλ + b δ ) .

4.2 Sampling from the joint posterior p(ζ, λ, δ|D)

As the full conditionals p(ζ|λ, δ, D), p(λ|ζ, δ, D) and p(δ|ζ, λ, D) are available, we follow a "Metropolis-within-Gibbs" strategy to sample the joint posterior p(ζ, λ, δ|D). In particular, the hyperparameters λ and δ will be sampled in a Gibbs step, while ζ will be sampled using a modified Langevin-Hastings algorithm. This approach is presented in Lambert and Eilers (2009) in the context of Bayesian density estimation (see also Lambert and Eilers (2005) for the use of MALA in a proportional hazards model). We adapt the algorithm of the latter reference to our EpiLPS methodology by incorporating the added value of our Laplacian-P-splines method. In particular, the variance-covariance matrix in the Langevin diffusion will be replaced by the variance-covariance matrix of LPSMAP. The correlation structure borrowed from LPSMAP improves convergence and chain mixing.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint

In what follows, we prefer to work under the log(·) parameterization for ρ, i.e. w = log(ρ) and denote by ζ = (θ ⊤ , w) ⊤ , the (K + 1)-dimensional vector of B-spline amplitudes and overdispersion w. Under this parameterization, the conditional posterior of ζ given λ and δ can be obtained from (19) by using the transformation method of random variables:

with the following log-likelihood under the reparameterization:

Let us denote by ζ (m−1)

∈ R (K+1) the state of the chain at iteration (m − 1). In the Langevin-Hastings algorithms, the proposal for the vector ζ at iteration m is a draw from the following multivariate Gaussian distribution:

where ϱ > 0 is a tuning parameter that has to be carefully chosen in order to reach a desired acceptance rate and Σ LH is the following block-diagonal variance-covariance matrix:

where Σ * is the K-dimensional covariance matrix obtained with LPSMAP. The gradient of log p( ζ|λ, δ, D) = ℓ( ζ; D) − 0.5λθ ⊤ P θ − b ρ exp(w) + a ρ w can be decomposed as follows:

16 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint and is analytically available (see Appendix A for more details). All the quantities related to the Langevin-Hastings proposal have been analytically derived, so that the draw in (24) can be obtained (for a given value of λ and δ). As in a classic MH algorithm, the next step consists in computing the acceptance probability:

where q(·, ·) denotes the (Gaussian) proposal distribution and p(·|λ, δ, D) the target (conditional) posterior distribution. Finally, we generate a uniform random variable u ∼ U(0, 1) and accept the proposed vector ζ

, ζ (prop) and reject it otherwise.

While iterating through the Metropolis-within-Gibbs algorithm, the tuning parameter ϱ is automatically adapted to reach the optimal acceptance rate of 0.57 (Haario et al., 2001; Atchadé and Rosenthal, 2005; Roberts and Rosenthal, 1998) . The pseudo-code below summarizes the LPSMALA algorithm.

LPSMALA algorithm to sample the posterior p(θ, ρ, λ, δ|D). , ϱ (m−1) Σ LH .

Compute acceptance probability π ζ (m−1)

, ζ

. 7: Draw u ∼ U(0, 1). 14: end for 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 The adaptive tuning part (on line 13) involves the step function H (z) = ϵI(z < ϵ) + zI(ϵ ≤ z ≤ A ) + A I(z > A ), with ϵ = 10 −4 and A = 10 4 , see Lambert and Eilers (2009) for details. Finally, the ratio q ζ (prop) , ζ (m−1)

, ζ (prop) entering the computation of the acceptance probability (line 6) is derived in Appendix B. 

Note also that S can be used to compute highest posterior density intervals of µ(t) at any point t. Using the renewal equation and the MCMC sample, one can apply the "plug-in" estimate of Section 3.1 and recover the following estimate of the instantaneous reproduction number at any given point t:

Also, using S , one can compute a highest posterior density interval of R(t) at any point t.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 

In this section, we implement a numerical study with four underlying epidemic scenarios to assess the accuracy with which EpiLPS is able to "track" the target reproduction number over time. EpiLPS results are compared with the estimate R() routine of the EpiEstim package using two sliding window options (the default weekly window of 7 days and a shorter window of 3 days). In each scenario, we simulate incidence data for an epidemic lasting 50 days according to a Poisson data generating process with mean counts governed by the renewal equation (see Eq.10). In total, S = 200 epidemics are simulated for each scenario and a cubic B-spline basis of size K = 20 with a second order penalty is used for inference in EpiLPS. In Scenario 1, a constant instantaneous reproduction number at R(t) = 2.5 is considered. In Scenario 2, an intervention strategy is replicated, so that R(t) = 2.5 and a sudden drop to R(t) = 0.7 occurs at day t = 20. The latter scenario allows to check whether EpiLPS is able to quickly react to such a brutal change in R(t). Scenario 3 implements a more wiggly structure for R(t) and finally, Scenario 4 considers the case of a vanishing epidemic with a monotonic decreasing reproduction number. Table 1 summarizes the target R(t) functions and the serial interval distribution for the considered scenarios. Table 2 reports for each scenario the average Bias, MSE, coverage of 90% and 95% credible envelopes and credible interval width for the R(t) estimator with EpiLPS (LPSMAP and LSPMALA) and EpiEstim. The average is computed over days t = 8, . . . , 50. For LPSMALA, a chain of length 5, 000 and a warm up period of 2, 000 is specified. Table 1 : Different scenarios for the simulation of epidemic data based on a Poisson data generating process for the incidence time series over T = 50 days.

Target R(t) Serial 19 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; For the first scenario, EpiEstim exhibits slightly lower bias than EpiLPS and the MSE is more or less of the same magnitude for both approaches. A general pattern that is apparent from Table 2 is that EpiLPS tends to overcover with LPSMAP (except for Scenario 2), while LPSMALA is closer to the nominal value in all scenarios. This phenomenon can be attributed to the following two main sources of approximations inherent to LPSMAP: (1) the Laplace approximation to the conditional posterior of the B-spline amplitudes and (2) the maximum a posteriori estimate of the hyperparameter vector that ignores the uncertainty related to the estimation of the roughness penalty parameter. In Scenario 2, EpiLPS has smaller bias than EpiEstim and the latter approach also suffers from mild undercoverage with a weekly sliding window. In Scenarios 3 and 4, where the target R(t) curve is less linear than in the previous scenarios, a more pronounced difference between EpiEstim and EpiLPS appears. EpiLPS shows better performance, especially in terms of coverage of credible intervals where EpiEstim exhibits severe undercoverage. Even when decreasing the sliding window to 3 days, the coverage probability is far from the nominal value.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101 https://doi.org/10. /2021 6 Application to observed case counts in infectious disease epidemics

In this section, the LPSMALA algorithm is applied on two historical outbreak datasets presented in Cori et al. (2013) . In particular, we consider the 2003 SARS outbreak in Hong Kong and the 2009 pandemic influenza in a school in Pennsylvania (USA). We use K = 40 B-splines with a second-order penalty and the serial interval distributions provided in the EpiEstim package Cori (2021) . The LPSMALA algorithm is implemented with a 22 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ;  chain of length 25,000. Acceptance rates for the generated chains are close to the optimal value of 57% and the posterior samples have converged according to the Geweke (1992) diagostic test (at the 1% level of significance). Figure 5 shows the smoothed epidemic curves and the estimated R(t) for the two outbreaks. Results for the SARS data show that the reproduction number reaches a first peak during the third week where R(t) = 9.66 (95% CI : 5.12 − 16.82) and a second more moderate peak around week 6 with R(t) = 2.79 (95% CI : 1.97 − 3.81). After day t = 43, the epidemic is under control and R(t) smoothly decays below 1. For the pandemic influenza in Pennsylvania, in the end of the second week R(t) is around 2.04 (95% CI : 1.24 − 3.07). During the middle of the third week, the situation is less severe and R(t) points below 1. As noted in Cori et al. (2013) , a few cases appeared in the last days of the epidemic generating an upward trend in R(t) estimates. 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10.1101/2021.12.02.21267189 doi: medRxiv preprint 6.2 Application on the SARS-CoV-2 pandemic

The EpiLPS methodology is illustrated on the SARS-CoV-2 pandemic using publicly available data from the Covid-19 Data Hub (Guidotti and Ardia, 2020) and its associated COVID19 package on CRAN (https://cran.r-project.org/package=COVID19). Countrylevel data on hospitalizations for Belgium, Denmark, Portugal and France from April, 5th 2020 to October 31st, 2021 is used with a uniform serial interval distribution over five days, i.e. φ = (0.2, 0.2, 0.2, 0.2, 0.2). In Figure 6 , the estimated reproduction number over the considered period is shown for the four countries. Results are obtained with the LPSMAP algorithm using K = 30 B-splines and a second order penalty. From a computational perspective, it takes less than 3 seconds to fit the EpiLPS model for the four countries with LPSMAP. The fitted reproduction numbers reflect the different waves of the COVID-19 pandemic and the recent rise in infections in the beginning of September 2021. 24 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10. 1101 EpiLPS (an acronym for Epidemiological modelling with Laplacian-P-Splines) is a fast and flexible tool for near real-time estimation of the instantaneous reproduction number R(t) during epidemic outbreaks. The tool is flexible in the sense that (penalized) spline based approximations provide smoothed estimates of R(t) with little computational effort and without the constraint of imposing any sliding window assumption that could potentially affect the timing and accuracy of the estimator. Moreover, the end user has the choice between a fully sampling-free approach (LPSMAP) or an efficient MCMC-based approach with Langevin diffusions (LPSMALA) for inference. The available EpiLPS package (https: //github.com/oswaldogressani) allows public health policy makers to analyze incoming data faster than existing methods relying on classic MCMC samplers, thus allowing them to be better informed when taking decisions on control measures for the ongoing SARS-CoV-2 epidemic. Simulation studies in this manuscript provide encouraging results and support EpiLPS as being a robust tool capable of a precise tracking of R(t) over time. The available EpiLPS software package is also straightforward to use with simple function calls that are specifically interesting for routine usage.

The EpiLPS project opens up several future research directions. A possible extension would be to formulate the EpiLPS model within a zero-inflated Poisson framework to cope with incindence time series characterized by an excess of zero counts. Another interesting extension would be to adapt the model to account for regional variation and imported cases.

In this appendix, we derive the analytical expressions of the gradient: ∇ ζ log p( ζ|λ, δ, D) = ∇ θ log p( ζ|λ, δ, D), ∂ log p( ζ|λ, δ, D) ∂w ⊤ , where the target function is log p( ζ|λ, δ, D) = ℓ( ζ; D) − 0.5λθ ⊤ P θ − b ρ exp(w) + a ρ w. Let us first concentrate on the partial derivatives with respect to the spline components:

∇ θ log p( ζ|λ, δ, D) = ∇ θ ℓ( ζ; D) − λP θ.

As already shown in Section 2.2, the gradient for the log-likelihood ∇ θ ℓ( ζ; D) is:

(y t + exp(w)) exp(θ ⊤ b(t)) exp(θ ⊤ b(t)) + exp(w) b k (t), k = 1, . . . , K.

The last term to be computed to recover the full gradient is:

where ∂ℓ( ζ; D) ∂w = T t=1 ∂ log Γ(y t + exp(w)) ∂w Term I − ∂ log Γ(exp(w)) ∂w Term II + ∂ exp(w)w ∂w Term III − ∂ ∂w (y t + exp(w)) log exp(θ ⊤ b(t)) + exp(w)

For Term I, using the chain rule, one recovers:

∂ log Γ(y t + exp(w)) ∂w = ∂ log Γ(y t + exp(w)) ∂Γ(y t + exp(w)) ∂Γ(y t + exp(w)) ∂(y t + exp(w)) ∂(y t + exp(w)) ∂w = Γ ′ (y t + exp(w)) Γ(y t + exp(w)) exp(w) = ψ(y t + exp(w)) exp(w),

27 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 5, 2021. ; https://doi.org/10. 1101 where ψ(·) is the digamma function. Using the same chain rule argument, one can easily show that for Term II:

∂ log Γ(exp(w)) ∂w = ψ(exp(w)) exp(w).

Term III is also trivial:

∂ exp(w)w ∂w = exp(w)(1 + w).

Term IV is as follows:

∂ ∂w (y t + exp(w)) log exp(θ ⊤ b(t)) + exp(w) = exp(w) log exp(θ ⊤ b(t)) + exp(w) +(y t + exp(w)) exp(w) exp(θ ⊤ b(t)) + exp(w) −1 = exp(w) log exp(θ ⊤ b(t)) + exp(w) +(y t + exp(w)) exp(θ ⊤ b(t)) + exp(w) −1

. (37) Gathering all the above intermediate results, we obtain the following derivative for (32):

∂ log p( ζ|λ, δ, D) ∂w = T t=1 exp(w) ψ(y t + exp(w)) − ψ(exp(w)) + (1 + w) − log exp(θ ⊤ b(t)) + exp(w) −(y t + exp(w)) exp(θ ⊤ b(t)) + exp(w)

28 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 5, 2021. ;

EpiNow2: Estimate Real-Time Case Counts and Time-Varying Epidemiological Parameters

Sampling theory of the Negative Binomial and logarithmic series distributions

On adaptive Markov chain Monte Carlo algorithms

On the estimation of the reproduction number based on misreported epidemic data

Markov chain Monte Carlo simulation methods in Econometrics

EpiEstim: estimate time varying reproduction numbers from epidemic curves

A new framework and software to estimate time-varying reproduction numbers during epidemics

Rcpp: Seamless R and C++ integration

Flexible smoothing with B-splines and penalties

Practical Smoothing: The Joys of P-splines

Estimating individual and household reproduction numbers in an emerging epidemic

Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images

Evaluating the accurating of sampling-based approaches to the calculation of posterior moments

Practical considerations for measuring the effective reproductive number

An approximate Bayesian approach for estimation of the reproduction number under misreported epidemic data

Laplace approximations for fast bayesian inference in generalized additive models based on P-splines

Covid-19 data hub

An adaptive Metropolis algorithm

Monte Carlo sampling methods using Markov chains and their applications

Robust specification of the roughness penalty prior distribution in spatially adaptive Bayesian P-splines models

Bayesian proportional hazards model with timevarying regression coefficients: A penalized Poisson regression approach

Bayesian density estimation from grouped continuous data

Bayesian P-splines

Equation of state calculations by fast computing machines

Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves

Maximum likelihood estimation for the Negative Binomial dispersion parameter

A spline-based time-varying reproduction number for modelling epidemiological outbreaks

Optimal scaling of discrete approximations to Langevin diffusions

Optimal scaling for various Metropolis-Hastings algorithms

Exponential convergence of Langevin distributions and their discrete approximations

Approximate Bayesian inference for latent Gaussian models by using Integrated nested Laplace approximations

Accurate approximations for posterior moments and marginal densities

How generation intervals shape the relationship between growth rates and reproductive numbers

Statistical estimation of the reproductive number from case notification data

To determine the analytical form of the ratio of proposal distributions in the Metropoliswithin-Gibbs algorithm, let us use the compact notation Υ (m−1) = ∇ ζ log p( ζ|λ, δ, D) ζ= ζ (m−1) ,, so that:The kernel of (39) is thus:Using a similar argument, one can show that the kernel of q ζ (m−1), ζ (prop) is:This can be used to compute the ratio: