Bayesian nowcasting with Laplacian-P-splines ============================================ * Bryan Sumalinab * Oswaldo Gressani * Niel Hens * Christel Faes ## Abstract During an epidemic, the daily number of reported infected cases, deaths or hospitalizations is often lower than the actual number due to reporting delays. Nowcasting aims to estimate the cases that have not yet been reported and combine it with the already reported cases to obtain an estimate of the daily cases. In this paper, we present a fast and flexible Bayesian approach to do nowcasting by combining P-splines and Laplace approximations. The main benefit of Laplacian-P-splines (LPS) is the flexibility and faster computation time compared to Markov chain Monte Carlo (MCMC) algorithms that are often used for Bayesian inference. In addition, it is natural to quantify the prediction uncertainty with LPS in the Bayesian framework, and hence prediction intervals are easily obtained. Model performance is assessed through simulations and the method is applied to COVID-19 mortality and incidence cases in Belgium. Keywords * Nowcasting * Laplacian-P-splines * Epidemic * COVID-19 * Reporting delay ## 1 Introduction Nowcasting is a term used for estimating the occurred-but-not-yet-reported-events (Donker et al., 2011; Van de Kassteele et al., 2019). In epidemiology, real-time updates of new symptomatic/infected individuals are helpful to assess the present situation and provide recommendations for rapid planning and for implementing essential measures to contain an epidemic out-break. The exact number of new daily cases is frequently subject to reporting delays, resulting in underreporting of the real number of infected individuals for that day. Failing to account for the reporting delays will lead to possibly biased predictions that might have an effect on policy making (Gutierrez et al., 2020). The main goal of nowcasting is to estimate the actual number of new cases by combining the (predicted) not-yet-reported-cases with the already reported cases. Several early references that establish the statistical framework for this type of problem can be found in the paper of Lawless (1994) which demonstrates how nowcasting can be used not only in disease surveillance but also in other contexts such as warranty and insurance claims. More recently, Höhle and an der Heiden (2014) applied nowcasting to the outbreak of Shiga toxin-producing Escherichia coli in Germany and also to the SARS-CoV-2 outbreak (Glöckner et al., 2020; Günther et al., 2021). Their approach is formulated within a hierarchical Bayesian framework that consists of estimating the epidemic curve by using a quadratic spline based on a truncated power basis function and the time-varying reporting delay distribution is approximated by a discrete time survival model. Van de Kassteele et al. (2019) pointed out that a potential drawback of such a method is the long computation time required by the Markov chain Monte Carlo (MCMC) algorithm and therefore proposed an alternative fast and flexible modelling strategy based on bivariate P-splines (Penalized B-splines). P-splines (Eilers and Marx, 1996) provide a flexible smoothing tool used to describe trends in the data. It introduces a penalty parameter that controls the roughness of the fit and counterbalances the flexibility of a rich B-splines basis (Eilers and Marx, 1996; Eilers et al., 2015). Other attractive features of the P-splines smoother are the relatively simple structure of the penalty matrix that is effortlessly computed and a natural extension to the Bayesian framework (Lang and Brezger, 2004). Based on the approach of Van de Kassteele et al. (2019), the number of cases are structured in a two-dimensional table (with calendar time as the first dimension and delay time as the second dimension), yielding the data matrix used as an input in the model. The reporting intensity is assumed to be a smooth surface and is modelled using two-dimensional P-splines. In this paper, we build upon the work of Van de Kassteele et al. (2019) by proposing a new nowcasting methodology based on Laplacian-P-splines (LPS) in a fully Bayesian framework. A key advantage of working with the Bayesian approach is the ease to obtain the predictive distribution and quantify the uncertainty associated with the predictions. In addition, the posterior distribution of the penalty parameter can be explored, and hence its uncertainty can be accounted for. The Laplace approximation uses a second-order Taylor expansion to approximate the posterior distribution of the regression parameters by a Gaussian density. It is a sampling-free method with the major advantage of faster computational time as opposed to MCMC approaches that are commonly used in Bayesian inference. Therefore, given the flexibility of previously mentioned P-splines smoothers and the computational benefit of Laplace approximations, it can be a helpful tool in the daily monitoring of new cases during an epidemic period. Laplacian-P-splines already proved to be useful in survival models (Gressani and Lambert, 2018; Gressani et al., 2022a), generalized additive models (Gressani and Lambert, 2021) and also in epidemic models for estimating the effective reproduction number (Gressani et al., 2022b). We build on the work of Gressani and Lambert (2021) to extend the Laplacian-P-splines methodology to nowcasting, thereby providing a fast and flexible (fully) Bayesian alternative to Van de Kassteele et al. (2019). To evaluate the (predictive) performance of our method, a simulation study is implemented and several performance measures are reported such as the mean absolute percentage error (MAPE), prediction interval coverage, and prediction interval width. Finally, we apply our method to the COVID-19 mortality data in Belgium for the year 2021. The data and R codes used to implement the method are available on GitHub through the link [https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git](https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git). ## 2 Methodology ### 2.1 Bayesian model formulation In this section, we follow the work of Van de Kassteele et al. (2019) to formulate a fully Bayesian model based on P-splines. Let *y**t,d* denote the number of cases that occurred at time *t* = 1, 2, …, *T* (corresponding to the calendar day) and are reported with a delay of *d* = 0, 1, 2…, *D* days. The information on cases can be summarized in matrix form: ![Formula][1] with cases that have not yet been reported (at time *T*) highlighted in bold. The not-yet-reported cases correspond to (*t, d*) combinations satisfying *t* > *T* − *d*. The main objective is to predict the total number of cases, ![Graphic][2], for *t* = *T* − (*D* − 1), …, *T* for which the nowcasted and already reported cases can be combined. Let 𝒟 := ***y*** = (*y*1, *y*2 …, *y**n*)*′* denote the vector of the observed number of cases by stacking the columns of matrix Y for the reported cases, where each entry corresponds to each (*t, d*) combination of reported cases *y**t,d*. The model assumes that the number of cases either follows a Poisson or a negative binomial (NB) distribution, i.e., *y**t,d* ∼ Poisson(*µ**t,d*) or *y**t,d* ∼ NB(*µ**t,d*, *ϕ*) with mean *µ**t,d* > 0. For the negative binomial, *ϕ* > 0 is an overdispersion parameter and the variance is ![Graphic][3]. Following Van de Kassteele et al. (2019), the (log) mean number of cases is modeled with two dimensional B-splines: ![Formula][4] where *β* is the intercept; *b**j*(·) and *b**k*(·) are univariate B-splines basis functions specified in the time and delay dimensions, respectively; and *z**l*(*t, d*) represents additional covariates with regression coefficients *β**l*. In matrix notation: ![Formula][5] where the matrices *B* and *Z* correspond to the basis functions and covariates, respectively, and vectors ***θ*** and ***β*** are the associated parameters to be estimated (Details in Appendix A1). In the philosophy of P-splines (Eilers and Marx, 1996), we use a rich (cubic) B-splines basis and counterbalance the associated flexibility by imposing a discrete roughness penalty on contiguous B-spline coefficients. For the two dimensional P-splines, the penalty can be obtained based on row-wise (direction of calendar time) and column-wise (direction of reporting delay) differences for matrix Θ = (*θ**j,k*) with *j* = 1, …, *K**T* and *k* = 1, …, *K**D* (Appendix A1) (see Durbán et al. (2002) and Fahrmeir et al. (2013) (pp. 507-508)). Let ![Graphic][6] and ![Graphic][7] denote the *m*th order row-wise and column-wise difference matrix with dimensions (*K**T* − *m*) × *K**T* and (*K**D* − *m*) × *K**D*, respectively. In this paper, we use a second order (*m* = 2) difference penalty. For ease of notation, let ![Graphic][8] and ![Graphic][9]. The difference matrix for vector ***θ*** can be obtained by expanding the difference matrix into ![Graphic][10] and ![Graphic][11], where ⊗ denotes the Kronecker product. Using this notation, the row-wise and column-wise difference penalty can be written, respectively, as: ![Formula][12] Let *λ**t* > 0 and *λ**d* > 0 denote the row-wise and column-wise penalty parameter, respectively. The penalty for the two dimensional B-splines is then given by: ![Formula][13] Let us define the penalty matrices ![Graphic][14] and ![Graphic][15], where *δ* is a small number (say *δ* = 10−6), to ensure that the penalty matrices are full rank and thus invertible. Following Lang and Brezger (2004), the penalty can be translated in the Bayesian framework by specifying a Gaussian prior for each column and row vector of Θ. We write Θ in terms of its columns and rows: ![Formula][16] The priors are then given by ![Graphic][17] and ![Graphic][18] for *k* = 1, …, *K**D* and *j* = 1, …, *K**T*, respectively. It follows that the joint prior for ***θ*** can be written as: ![Formula][19] Gaussian priors are assumed for ***β*** and ***θ***, namely ![Graphic][20] with *V****β*** = *ζI**p*+1 (small *ζ*, e.g. *ζ* = 10−5) and (***θ***|***λ***) ∼ *𝒩**dim*(***θ***)(****, *𝒫*−1(***λ***)), where ***λ*** = (*λ**t*, *λ**d*)*′* is the penalty vector and ![Graphic][21] the global penalty matrix. Denote by *X* = (*B, Z*) the global design matrix, ***ξ*** = (***β****′*, ***θ****′*)*′* the latent parameter vector and ![Graphic][22] the precision matrix for ***ξ***. From here, we focus on the negative binomial model for the number of cases. The Poisson model is described in detail in Appendix A3. The full Bayesian (negative binomial) model is summarized as follows: ![Formula][23] where *𝒢*(*a, b*) denotes a Gamma distribution with mean *a/b* and variance *a/b*2. This robust prior specification on the penalty parameters follows from Jullion and Lambert (2007). They have shown that when the hyperparameters *a**δ*, *b**δ* are chosen to be equal and small enough (say 10−4), then the resulting fit is robust to the value of *ν* (e.g. *ν* = 3 in this paper). ### 2.2 Laplace approximation to the conditional posterior of *ξ* To obtain a Laplace approximation for the set of regression parameters ***ξ*** = (***β****′*, ***θ****′*)*′*, we use the posterior of ***ξ*** conditional on the vector of penalty parameters ***λ***. The gradient and Hessian of the (log) conditional posterior are analytically derived and used in a Newton-Raphson algorithm to obtain the Gaussian approximation to the conditional posterior distribution of ***ξ***. Note that for a negative binomial distributed *y**i* with mean *E*(*y**i*) = *µ**i*, the probability distribution can be written as an exponential dispersion family given by *p*(*y**i*; *γ, ϕ*) = exp{[*y**i**γ**i* − *b*(*γ**i*)]*/a*(*ϕ*)+*c*(*y**i*, *ϕ*)*}*, where ![Graphic][24], *a*(*ϕ*) = *ϕ* and Γ(·) is the gamma function. Thus, for fixed *ϕ*, a negative binomial regression model is a generalized linear model (GLM) where *µ**i* is linked on the linear predictor through the link function *g*(*µ**i*) such that ![Graphic][25] (see Agresti (2013) for a detailed account about GLMs). Here, we use the log-link function *g*(*µ**i*) = log(*µ**i*). The log-likelihood is given by log ![Graphic][26] and it can be shown that for a negative binomial model, the gradient and Hessian are: ![Formula][27] where *W* = diag(*w*1, …, *w**n*), *w**i* = [*V ar*(*y**i*)(*g**′*(*µ**i*))2]−1, *D* = diag(*g**′*(*µ*1), …, *g**′*(*µ**n*)), ![Graphic][28] and *M* = diag(*y*1 −*µ*1, …, *y**n* − *µ**n*) (see details in Appendix A2). Using Bayes’ rule, the posterior of ***ξ*** conditional on the penalty vector ***λ*** and overdispersion parameter *ϕ* is: ![Formula][29] The gradient and Hessian for the log-conditional posterior of ***ξ*** are: ![Formula][30] The above gradient and Hessian can then be used in a Newton-Raphson algorithm to obtain the mode of the conditional posterior of ***ξ***. After convergence, the Laplace approximation of the conditional posterior of ***ξ*** is a multivariate Gaussian density denoted by ![Graphic][31] where ![Graphic][32] is the mean/mode and ![Graphic][33] is the variance-covariance. ### 2.3 Hyperparameter optimization In this section, we derive the (approximate) posterior distribution of the hyperparameters belonging to the penalization part, i.e., ***λ*** and ***δ*** and overdispersion part *ϕ*. We first derive the joint posterior of ***λ, δ*** and *ϕ* and then proceed with an integration to isolate the marginal posterior distribution of ***λ*** and *ϕ*. Let ***η*** = (*λ**t*, *λ**d*, *δ**t*, *δ**d*, *ϕ*)*′* denote the vector of hyperparameters. Using Bayes’ theorem, the marginal posterior of ***η*** is: ![Formula][34] Following Rue et al. (2009), the above posterior can be approximated by replacing *p*(***ξ***|***η***, *𝒟*) with ![Graphic][35] obtained in Section 2.2 and by evaluating the latent vector at ![Graphic][36]. Note that ![Graphic][37], where ![Graphic][38] corresponds to the *i* th row of the design matrix *X*. Also, the determinant ![Graphic][39] is ![Graphic][40]. Hence, the approximated marginal posterior of ***η*** is: ![Formula][41] To obtain the joint marginal posterior of ***λ*** and *ϕ*, we need to integrate out the hyperparameters *δ**t* and *δ**d* from ![Graphic][42] as follows: ![Formula][43] Hence, ![Formula][44] The posterior mode (obtained via Newton-Raphson) is used as a point estimate for ***λ*** and *ϕ*. To ensure numerical stability, we log-transformed the penalty vector ***v*** = (*v**t*, *v**d*)*′* = (log(*λ**t*), log(*λ**d*))*′* and overdispersion parameter *v**ϕ* = log(*ϕ*). Using the method of transformations, we obtain the Jacobian of the transformation *J* = exp(*v**d*) exp(*v**t*) exp(*v**ϕ*). Hence, the joint posterior of ***v*** and *v**ϕ* is given by: ![Formula][45] Moreover, the joint log-posterior of ***v*** and *v**ϕ* is given by: ![Formula][46] Rather than maximizing the above joint posterior, we maximized the conditional posterior of *v**t* by keeping *v**ϕ* fixed at its method of moment estimate and setting *v**d* to zero. This approach offers faster computational speed while still producing nearly identical results as when maximizing the joint posterior. ### 2.4 Nowcasting with prediction interval To obtain the mean nowcast estimate with the prediction interval, note that ![Graphic][47] or equivalently ![Graphic][48], where ![Graphic][49]. Thus, ![Graphic][50]. The mean estimate for the not-yet-reported cases is calculated as ![Graphic][51] for all (*t, d*) combinations with *t* > *T* − *d*. Then, the estimate for the total number of cases for each *t* is obtained by summing the already reported cases and the mean estimate for the not-yet-reported cases, i.e. ![Graphic][52]. The prediction interval for the nowcast values is obtained by sampling from the posterior predictive distribution of the log-mean number of cases by following a five-step procedure: 1. For each (*t, d*) combinations with *t* > *T* − *d* (corresponding to the not-yet-reported cases), generate 1000 random samples ![Graphic][53] from a Gaussian distribution with mean ![Graphic][54] and variance ![Graphic][55]. 2. Exponentiate the sampled values from the previous step to obtain the average reporting intensities ![Graphic][56]. 3. Compute the average prediction for the not-yet-reported cases for each *t*. That is, compute ![Graphic][57]. 4. For each *t*, sample ![Graphic][58] containing 1000 values from ![Graphic][59]. 5. Finally, compute the quantiles of the sampled values ![Graphic][60] that correspond to the desired prediction interval. ### 2.5 Delay distribution To obtain the smooth estimate of the delay distribution, only the first term on the right hand side of equation (1) is used (excluding the day of the week effects), as explained in the paper of Van de Kassteele et al. (2019). Specifically, the procedure to compute the delay distribution is as follows: 1. Compute the contribution of the smoothing term in equation (1) to the reporting intensity for all (*t, d*) combinations: ***µ***smooth = exp(*B****θ***). 2. Arrange ***µ***smooth into a *T* × (*D* + 1) matrix with entries ![Graphic][61]. 3. For each *t* = 1, …, *T*, compute the reporting delay distribution given by: ![Graphic][62] ## 3 Simulations A simulation study is implemented in order to evaluate the predictive performance of the proposed method. The procedure to perform the simulations is as follows: 1. Consider a function *f* (*t*) that represents the mean epidemic curve of all cases such that *µ*(*t*) = exp(*f* (*t*)) for *t* = 1, …, *T*. 2. For each *t*, generate a random sample *y**t* from a negative binomial distribution with mean *µ*(*t*) and fixed overdispersion parameter (we choose the value of 10). 3. To account for possible delays *d* = 0, 1, 2, …, *D*, generate samples from a multinomial distribution with probabilities (*p*, *p*1, *p*2, …, *p**D* such that![Graphic][63]), i.e., ![Formula][64] This sample represents the reported number of cases for each (*t, d*) combination. 4. Repeat steps 1 to 3 for 500 times to generate 500 possible realizations. We consider two functions *f* (*t*) inspired from the paper of Noufaily et al. (2016) given by: ![Formula][65] for *t* = 1, …, 365. That is, we assume a one year (365 days) time window in the simulation. In terms of delay probabilities, we consider a maximum delay of *D* = 7 days with probabilities (*p*, *p*1, *p*2, …, *p*7) = (0.0, 0.1, 0.4, 0.2, 0.1, 0.1, 0.05, 0.05), as illustrated in Figure 1. For example, at time *T* − 6, only the cases that have a delay of 7 days (with probability *p*7 = 0.05) are not yet reported, or equivalently, 95% of the cases have already been reported. For time *T* − 5, only 90% of the cases have been reported, as cases with delays of 6 and 7 days are yet to be reported. On the nowcast day (time *T*), no case is reported, that is, the delay is 100%. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F1) Figure 1: Illustration of the delay probabilities considered in the simulation. Two scenarios are considered for the first function *f*1(*t*): (i) the first one is having a small number of cases with values of *θ*1 = 3 and *θ*2 = 1; (ii) the second scenario has relatively larger cases with *θ*1 = 3 and *θ*2 = 2 and an additional factor of 50 is added on the mean function such that *µ*(*t*) = 50 + exp(*f* (*t*)). These functions are denoted by *f*11(*t*) and *f*12(*t*), respectively. Similarly, two scenarios are also considered for the second function *f*2(*t*): (i) the first scenario has values *θ*1 = 0, *θ*2 = 0.4 and *θ*3 = 0.2; (ii) the second scenario has values *θ*1 = 1.5, *θ*2 = 0.4 and *θ*3 = 0.2. We denote these functions by *f*21(*t*) and *f*22(*t*), respectively. The first function has a symmetric curve with three peaks as shown in Figures 2a and 2b. On the other hand, the second function is not periodic as opposed to the first function (see Figures 2c and 2d). The plot for (one realization of) simulated cases based on these functions are shown in Figure 3 with (dashed) vertical lines corresponding to the different nowcast dates. ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F2) Figure 2: Mean epidemic curves considered in the simulations. ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F3) Figure 3: Simulated cases for different epidemic curves. The dashed vertical lines correspond to different nowcast dates. For each generated data set, we fit the proposed models (LPS-NB and LPS-Poisson) as well as the method proposed by Van de Kasteele et al. (VK) on the data, and compute the desired accuracy measures. The mean absolute percentage error (MAPE) and 95% prediction interval coverage are chosen to measure the predictive accuracy of our methodology. The MAPE is scale-independent and measures the error of prediction from the actual value in terms of percentage. The formula to compute the MAPE at a given calendar day *t* is given by: ![Formula][66] where *S* = 500 is the number of generated realizations, *y**st* is the true value (generated from step 2 of the simulation in Section 3 at iteration *s*) and ![Graphic][67] is the corresponding predicted value (mean nowcast estimate ![Graphic][68] in Section 2.4). The prediction interval coverage is obtained by determining the percentage of (true) unreported cases that fall within the interval. Moreover, for each simulation, the interval width is obtained, which is the difference between the upper and lower bound of the prediction interval. In the simulation, we fix the nowcast date at the end of the month (from March to November). For each nowcast date, the prediction measures are computed for the dates having unreported cases, that is, for *t* = *T* − (*D* − 1), …, *T*. As we have a maximum delay of 7 days, there will be seven dates (including the nowcast date) that involves prediction of unreported cases. All simulation results are presented in Appendix A5. Tables 1-4 provide a summary of the results for the prediction measures on the nowcast date (time *T*) for the different functions being considered. Results show that the MAPE values for LPS-NB and LPS-Poisson are relatively similar for the different function for each month. These MAPE values are a bit larger for functions *f*11(*t*) and *f*21(*t*) (that exhibit a smaller number of cases) ranging from 31% to 95%. However, for functions *f*12(*t*) and *f*22(*t*) (characterized by a larger number of cases), the MAPE for LPS-NB and LPS-Poisson are smaller, ranging between 32% and 47%. Van de Kassteele et al.’s method, on the other hand, has generally higher MAPE values (up to 1016%) compared to our LPS method, especially for functions *f*12(*t*) and *f*22(*t*). This means that LPS has better prediction accuracy, especially for scenarios with higher numbers of cases. View this table: [Table 1:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T1) Table 1: Performance measures on the nowcast date for function *f*11(*t*): LPS-NB - LPS model with a negative binomial distribution for the number of cases; LPS-Poisson - LPS model with a Poisson distribution for the number of cases; VK - Methodology of Van de Kassteele et al. (2019) View this table: [Table 2:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T2) Table 2: Performance measures on the nowcast date for function *f*12(*t*): LPS-NB - LPS model with a negative binomial distribution for the number of cases; LPS-Poisson - LPS model with a Poisson distribution for the number of cases; VK - Methodology of Van de Kassteele et al. (2019) View this table: [Table 3:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T3) Table 3: Performance measures on the nowcast date for function *f*21(*t*): LPS-NB - LPS model with a negative binomial distribution for the number of cases; LPS-Poisson - LPS model with a Poisson distribution for the number of cases; VK - Methodology of Van de Kassteele et al. (2019) View this table: [Table 4:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T4) Table 4: Performance measures on the nowcast date for function *f*22(*t*): LPS-NB - LPS model with a negative binomial distribution for the number of cases; LPS-Poisson - LPS model with a Poisson distribution for the number of cases; VK - Methodology of Van de Kassteele et al. (2019) In terms of prediction interval coverage, the LPS-NB has more stable coverage rates and is closer to the nominal 95% nominal prediction interval across the different functions and nowcast dates ranging between 88% and 96%. The method of Van de Kassteele et al. on the other hand, tends to have high coverage rates (up to 100%) resulting from wide prediction intervals, indicating higher uncertainty in its predictions. This is mainly the case when the number of cases is large, as for functions *f*12(*t*) and *f*22(*t*). The LPS-Poisson, however, generally has the narrowest prediction interval widths resulting in lower coverage rates. This is expected since we simulate the data from the negative binomial distribution with an overdispersion parameter which is not accounted for by the Poisson distribution. In general, these findings are also true for prediction at any other day (Appendix A5), such that the LPS-NB model performs better, particularly for larger cases, having lower MAPE and having coverage closer to the nominal 95% level. ## 4 Real Data application We apply our method to COVID-19 mortality data in Belgium for 2021 and incidence data for 2022. The raw data is available on the website of the Sciensano research institute (see sciensano.be/covid19data). The data contains the cumulative number of cases, reported up to the day of the file. The file is updated every day, and in this way, the number of cases and reporting delays are obtained. The data is structured in matrix format with the date of death/confirmed case as rows and number of days of reporting delay as columns. Figure 4 shows the total number of cases with (dashed) vertical lines corresponding to different nowcast dates used for illustration. As no cases are reported immediately on the day of death/confirmed case in the data, all cases have a delay of at least one day. The proportion of cases reported for each delay is shown in Table 5. Majority of the cases (66% for mortality and 61% for incidence) were reported with a delay of 2 days. Followed by a delay of 1 day (14%) for mortality and 3 days (37%) for incidence. The rest of the delays account for very small percentages of the cases. The reporting delay is truncated to a maximum number of seven days. View this table: [Table 5:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T5) Table 5: Proportion of cases reported for each delay in the mortality data. View this table: [Table 6:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T6) Table 6: Proportion of cases reported for each delay in the incidence data. ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F4) Figure 4: Covid-19 death (a) and incidence (b) cases in Belgium. Dashed vertical lines correspond to different nowcast dates. In this section, we present the model results at different nowcast dates. In the analysis, we choose *K**T* = 40 and *K**D* = 10 as the number of B-splines for the time and delay dimension, respectively. Moreover, a cubic B-splines basis with a second order penalty is used. In addition, we also consider the day of the week effect as an additional covariate in the model. Figure 5 presents the nowcasting results using the mortality data at different nowcast dates (namely, at every end of the month). The blue color represents the reported cases for the past 14 days on the nowcast date, the gray color are cases that have not yet been reported, and the orange line corresponds to the nowcast prediction with the 95% prediction interval. It can be seen that most of the nowcast predictions on the nowcast date are fairly close to the observed cases (gray). In addition, all of the observed cases fall within the prediction interval. For the incidence data, we used a maximum delay of 5 days. The nowcast plots are shown in Figure 6. The accuracy of nowcast estimates for incidence data is comparatively lower than that of mortality data in most cases. The nowcast results tend to exhibit varying degrees of underestimation. This means that some nowcast values are close to the observed values, while others are further away, and some are moderately close. ![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F5.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F5) Figure 5: Nowcast for mortality data with different nowcast dates. Blue - reported cases ; Gray - not-yet-reported cases; Orange - nowcast with prediction interval. ![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F6.medium.gif) [Figure 6:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F6) Figure 6: Nowcast for incidence data with different nowcast dates. Blue - reported cases ; Gray - not-yet-reported cases; Orange - nowcast with prediction interval. The plot for the estimated delay density is shown in Appendix A4. The steps to obtain the timevarying delay distribution is outlined in Section 2.5. The plots show that the delay distribution is fairly constant through time for both the incidence and mortality data. The density is highest for the delay of 2 days except in the beginning of the year for incidence data. This confirms the observed reporting intensity in our data that most cases are reported with a two-day delay. ## 5 Conclusion This paper presents a novel Bayesian approach of the nowcasting model proposed by Van de Kassteele et al. (2019). During a pandemic or outbreak, it is very important to get an idea on the number of actual cases in real-time. One advantage of using the Bayesian framework is that the prediction intervals, which is the main goal of nowcasting, can be naturally obtained. Our proposed methodology is based on a combination of Laplace approximations and P-splines, making it both fast and flexible. Based on the results presented in this paper, Bayesian nowcasting with Laplacian P-Splines seems to be a promising tool. We performed a simulation study to evaluate the (predictive) performance of our method. Simulation results show that our proposed method, when assuming a negative binomial distribution, performs better than Van de Kassteele et al. (2019), especially for larger cases. In particular, our model has lower MAPE values implying better prediction accuracy and the prediction interval coverage are fairly close to the nominal (95%) prediction interval. In addition, the method of Van de Kassteele et al. (2019) tends to have very wide prediction interval widths indicating higher uncertainty in its predictions. Moreover, we also applied our method to the COVID-19 mortality and incidence data in Belgium. The nowcast predictions for the mortality data seem to be fairly close to the actual observed values, and all of the 95% prediction intervals for the different nowcast dates being considered contained the observed values. In the case of incidence data with larger numbers of cases, the nowcast results tend to exhibit varying degrees of underestimation. However, it is important to take note that we have no data for the nowcast day, and there are very few reported cases (0.1%) with a one-day delay for the incidence data. This makes it much more difficult to produce accurate nowcast predictions. While nowcasting can provide valuable real-time information and predictions, it also has certain drawbacks. Nowcasting heavily relies on real-time data, which may not always be readily available, or of good quality. Biases in the the data (such as double-counting of cases) that are corrected at a later time, can impact the nowcasting prediction. In addition, gaps in data collection, processing, or dissemination can impact the timeliness and effectiveness of nowcasting predictions. Despite these drawbacks, nowcasting has been a very important tool for obtaining timely information and short-term predictions. While the proposed nowcasting models are complex, the R codes used to implement the methods are made freely available on a GitHub account, so that the proposed method is available for practical use. Finally, there are several potential extensions to the LPS methodology for nowcasting. One possible extension is to incorporate the LPS nowcasting model into the estimation of the reproduction number in the (LPS) method proposed by Gressani et al. (2022b). This would allow for a more comprehensive understanding of the underlying dynamics and transmission of the disease. Additionally, it would be beneficial to consider the correlation of cases for each time point (*t*). By doing this, the model can capture temporal dependencies and may provide better estimates of the variability. Furthermore, accounting for spatial correlation and other covariate effects would be valuable. Incorporating spatial information can account for local variations while considering additional covariates can provide a better understanding of the factors influencing disease spread. ## Data Availability The raw data is available on the website of the Sciensano research institute at [https://epistat.sciensano.be/covid/covid19\_historicaldata.html](https://epistat.sciensano.be/covid/covid19_historicaldata.html). In addition, the cleaned version of the data is available on GitHub at [https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git](https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git). [https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git](https://github.com/bryansumalinab/Laplacian-P-spline-nowcasting.git) [https://epistat.sciensano.be/covid/covid19\_historicaldata.html](https://epistat.sciensano.be/covid/covid19_historicaldata.html) ## 6 Appendix ### A1 In equation (1), we have a two-dimensional B-spline model: ![Formula][69] where * *B* = *B**D* *⊗ B**T* is the Kronecker product of B-splines matrices *B**T* and *B**D* with dimension *T* × *K**T* and *D* × *K**D*, given by: ![Formula][70] As a result, B has dimension *TD* × *K**T* *K**D*. * Z is the design matrix for additional covariates with dimension *TD* × (*p* + 1) and corresponding parameters ***β***. This is given by: ![Formula][71] * ***θ*** = vec(Θ) has dimension *K**T* *K**D* × 1 where vec(Θ) denotes the vectorization of a matrix Θ. It is obtained by vectorizing the matrix: ![Formula][72] of B-splines coefficients *θ**j,k* by stacking the columns of Θ, that is, ![Formula][73] ### A2. Gradient and Hessian of an negative binomial GLM #### A2.1 Negative Binomial as an exponential dispersion family Let *y ∼ NB*(*μ, ϕ*) then *y* has a probability distribution ![Formula][74] where *E*(*y*) = *μ* and ![Graphic][75]. The negative binomial probability distribution can be expressed as an exponential dispersion family: ![Formula][76] where *E*(*y*) = *b**′*(*γ*) = *μ* and *V ar*(*y*) = *a*(*ϕ*)*b**′′*(*γ*). We have: ![Formula][77] with ![Graphic][78] and *a*(*ϕ*) = *ϕ*. #### A2.2 Gradient and Hessian of an exponential dispersion family Suppose we have independent observations *y*1, …, *y**n* where each *y**i* has a probability distribution from an exponential dispersion family given by ![Formula][79] The GLM links the mean *E*(*y**i*) = *μ**i* to the linear predictor ![Graphic][80] through the link function *g*(*μ**i*), such that *ψ**i* = *g*(*μ**i*). The log-likelihood is given by ![Graphic][81] ##### Gradient We have: ![Formula][82] where ![Formula][83] Thus, ![Formula][84] Let ***ξ*** = (*ξ*1, …, *ξ**p*)*′*. The gradient is given by ![Formula][85] In matrix notation *∇****ξ****L* = *X**′**WD*(*y*−*μ*), where *W* = diag(*w*1, …, *w**n*), *w**i* = [*V ar*(*y**i*)(*g**′*(*μ**i*))2]−1 and *D* = diag(*g**′*(*μ*1), …, *g**′*(*μ**n*)). ##### Hessian We have ![Formula][86] Let *V* = diag(*v*1, *v*2, …, *v**n*) where ![Graphic][87] and *M* = diag(*y*1 − *μ*1, …, *y**n* − *μ**n*). Then, in matrix notation, the Hessian can be written as ![Formula][88] Using a log link function, i.e., *g*(*μ**i*) = log(*μ**i*), we recover ![Graphic][89] and ![Graphic][90]. Also, for the negative binomial model, we have ![Graphic][91]. Hence, ![Formula][92] Thus,![Graphic][93]. ### A3. Derivations for the LPS-Poisson model The full Bayesian LPS-Poisson model is given by: ![Formula][94] #### A3.1 Laplace approximation to the conditional posterior of *ξ* Note that for a Poisson distributed *y**i* with mean *μ**i*, we have *f* (*y**i*) ∝ exp (*y**i**γ**i* − *s*(*γ**i*)), where *γ**i* = log(*μ**i*) and *s*(*γ**i*) = exp(*γ**i*) = *μ**i*. Using Bayes’ rule, the posterior of ***ξ*** conditional on the penalty vector ***λ*** is: ![Formula][95] where **ℒ** (***ξ***; **𝒟**) denotes the likelihood function. Hence, the log-conditional posterior of ***ξ*** can be obtained as: ![Formula][96] where ≐ denotes equality up to an additive constant. It is known for Poisson generalized linear models that the gradient (score vector) and Hessian of the first term on the right hand side of (2) are given, respectively, by (see Agresti (2013), p. 136): ![Formula][97] where *W* = diag(*w*1, …, *w**n*) is a diagonal matrix with *w**i* = (*V ar*(*y**i*)[*g**′*(*μ**i*)]2)−1 = *μ**i* and *g*(·) is the log-link function such that *g*(*μ**i*) = log(*μ**i*). Hence, the gradient and Hessian for the log-conditional posterior of ***ξ*** are: ![Formula][98] #### A3.2 Hyperparameter optimization Let ***η*** = (*λ**t*, *λ**d*, *δ**t*, *δ**d*)*′* denote the vector of hyperparameters. The approximated marginal posterior of ***η*** is: ![Formula][99] The approximated marginal posterior of ***λ*** is given by: ![Formula][100] To ensure numerical stability, we log-transformed the penalty vector ***v*** = (*v**t*, *v**d*)*′* = (log(*λ**t*), log(*λ**d*))*′*. The log-posterior of ***v*** is then given by: ![Formula][101] ### A4. Estimated delay density plot ![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F7.medium.gif) [Figure 7:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F7) Figure 7: Estimated delay density for mortality data at different nowcast dates. ![Figure 8:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/06/28/2022.08.26.22279249/F8.medium.gif) [Figure 8:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/F8) Figure 8: Estimated delay density for incidence data at different nowcast dates. ### A5. Simulation results #### Results for function *f*11(*t*) View this table: [Table 7:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T7) Table 7: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-NB. View this table: [Table 8:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T8) Table 8: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using the method of Van de Kassteele et al. View this table: [Table 9:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T9) Table 9: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-Poisson. View this table: [Table 10:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T10) Table 10: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-NB. View this table: [Table 11:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T11) Table 11: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using the method of Van de Kassteele et al. View this table: [Table 12:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T12) Table 12: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-Poisson. View this table: [Table 13:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T13) Table 13: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-NB. View this table: [Table 14:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T14) Table 14: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using the method of Van de Kassteele et al. View this table: [Table 15:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T15) Table 15: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*11(*t*) using LPS-Poisson. #### Results for function *f*12(*t*) View this table: [Table 16:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T16) Table 16: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-NB. View this table: [Table 17:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T17) Table 17: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using the method of Van de Kassteele et al. View this table: [Table 18:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T18) Table 18: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-Poisson. View this table: [Table 19:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T19) Table 19: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-NB. View this table: [Table 20:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T20) Table 20: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using the method of Van de Kassteele et al. View this table: [Table 21:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T21) Table 21: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-Poisson. View this table: [Table 22:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T22) Table 22: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-NB. View this table: [Table 23:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T23) Table 23: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using the method of Van de Kassteele et al. View this table: [Table 24:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T24) Table 24: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*12(*t*) using LPS-Poisson. #### Results for function *f*21(*t*) View this table: [Table 25:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T25) Table 25: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-NB. View this table: [Table 26:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T26) Table 26: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using the method of Van de Kassteele et al. View this table: [Table 27:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T27) Table 27: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-Poisson. View this table: [Table 28:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T28) Table 28: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-NB. View this table: [Table 29:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T29) Table 29: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using the method of Van de Kassteele et al. View this table: [Table 30:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T30) Table 30: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-Poisson. View this table: [Table 31:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T31) Table 31: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-NB. View this table: [Table 32:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T32) Table 32: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using the method of Van de Kassteele et al. View this table: [Table 33:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T33) Table 33: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*21(*t*) using LPS-Poisson. #### Results for function *f*22(*t*) View this table: [Table 34:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T34) Table 34: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-NB. View this table: [Table 35:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T35) Table 35: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using the method of Van de Kassteele et al. View this table: [Table 36:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T36) Table 36: Mean absolute percentage error (MAPE) for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-Poisson. View this table: [Table 37:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T37) Table 37: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-NB. View this table: [Table 38:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T38) Table 38: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using the method of Van de Kassteele et al. View this table: [Table 39:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T39) Table 39: Prediction interval coverage for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-Poisson. View this table: [Table 40:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T40) Table 40: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-NB. View this table: [Table 41:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T41) Table 41: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using the method of Van de Kassteele et al. View this table: [Table 42:](http://medrxiv.org/content/early/2023/06/28/2022.08.26.22279249/T42) Table 42: Prediction interval width for the number of cases at day *T* − *t* (*t* = 0, …, 6), as well as the average over these days, for the function *f*22(*t*) using LPS-Poisson. ## Footnotes * The main difference between the earlier and updated edition of the manuscript lies in their assumptions regarding the distribution of the number of cases. In the previous version, a Poisson distribution was assumed, whereas the revised version assumes a negative binomial distribution for the number of cases. * Received August 26, 2022. * Revision received June 28, 2023. * Accepted June 28, 2023. * © 2023, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. Agresti, A. (2013). Categorical Data Analysis. John Wiley & Sons. 2. Donker, T., van Boven, M., van Ballegooijen, W. M., van’t Klooster, T. M., Wielders, C. C., and Wallinga, J. (2011). Nowcasting pandemic influenza A/H1N1 2009 hospitalizations in the Netherlands. European Journal of Epidemiology, 26(3):195–201. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s10654-011-9566-5&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21416274&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F06%2F28%2F2022.08.26.22279249.atom) 3. Durbán, M., Currie, I., and Eilers, P. H. C. (2002). Using P-splines to smooth two-dimensional Poisson data. In Proceedings of 17th International Workshop on Statistical Modelling, Chania, Crete, pages 207–214. 4. Eilers, P. H. C. and Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11(2):89–121. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1214/ss/1038425655&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1996WA79400001&link_type=ISI) 5. Eilers, P. H. C., Marx, B. D., and Durbán, M. (2015). Twenty years of P-splines. SORT: Statistics and Operations Research Transactions, 39(2):149–186. 6. Fahrmeir, L., Kneib, T., Lang, S., and Marx, B. (2013). Regression: Models, Methods and Applications. Springer Science & Business Media. 7. Glöckner, S., Krause, G., and Höhle, M. (2020). Now-casting the COVID-19 epidemic: The use case of Japan, March 2020. medRxiv. 8. Gressani, O., Faes, C., and Hens, N. (2022a). Laplacian-P-splines for Bayesian inference in the mixture cure model. Statistics in Medicine, 41(14):2602–2626. 9. Gressani, O. and Lambert, P. (2018). Fast Bayesian inference using Laplace approximations in a flexible promotion time cure model based on P-splines. Computational Statistics & Data Analysis, 124:151–167. 10. Gressani, O. and Lambert, P. (2021). Laplace approximations for fast Bayesian inference in generalized additive models based on P-splines. Computational Statistics & Data Analysis, 154:107088. 11. Gressani, O., Wallinga, J., Althaus, C. L., Hens, N., and Faes, C. (2022b). Epilps: A fast and flexible bayesian tool for estimation of the time-varying reproduction number. PLoS computational biology, 18(10):e1010618. 12. Günther, F., Bender, A., Katz, K., Küchenhoff, H., and Höhle, M. (2021). Nowcasting the COVID-19 pandemic in Bavaria. Biometrical Journal, 63(3):490–502. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F06%2F28%2F2022.08.26.22279249.atom) 13. Gutierrez, E., Rubli, A., and Tavares, T. (2020). Delays in death reports and their implications for tracking the evolution of COVID-19. Available at SSRN 3645304. 14. Höhle, M. and an der Heiden, M. (2014). Bayesian nowcasting during the STEC O104: H4 outbreak in Germany, 2011. Biometrics, 70(4):993–1002. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/biom.12194&link_type=DOI) 15. Jullion, A. and Lambert, P. (2007). Robust specification of the roughness penalty prior distribution in spatially adaptive Bayesian P-splines models. Computational Statistics & Data Analysis, 51(5):2542–2558. 16. Lang, S. and Brezger, A. (2004). Bayesian P-splines. Journal of Computational and Graphical Statistics, 13(1):183–212. 17. Lawless, J. (1994). Adjustments for reporting delays and the prediction of occurred but not reported events. Canadian Journal of Statistics, 22(1):15–31. 18. Noufaily, A., Farrington, P., Garthwaite, P., Enki, D. G., Andrews, N., and Charlett, A. (2016). Detection of infectious disease outbreaks from laboratory data with reporting delays. Journal of the American Statistical Association, 111(514):488–499. 19. Rue, H., Martino, S., and Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1467-9868.2008.00700.x&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000264374200002&link_type=ISI) 20. Van de Kassteele, J., Eilers, P. H. C., and Wallinga, J. (2019). Nowcasting the number of new symptomatic cases during infectious disease outbreaks using constrained P-spline smoothing. Epidemiology, 30(5):737. [1]: /embed/graphic-1.gif [2]: /embed/inline-graphic-1.gif [3]: /embed/inline-graphic-2.gif [4]: /embed/graphic-2.gif [5]: /embed/graphic-3.gif [6]: /embed/inline-graphic-3.gif [7]: /embed/inline-graphic-4.gif [8]: /embed/inline-graphic-5.gif [9]: /embed/inline-graphic-6.gif [10]: /embed/inline-graphic-7.gif [11]: /embed/inline-graphic-8.gif [12]: /embed/graphic-4.gif [13]: /embed/graphic-5.gif [14]: /embed/inline-graphic-9.gif [15]: /embed/inline-graphic-10.gif [16]: /embed/graphic-6.gif [17]: /embed/inline-graphic-11.gif [18]: /embed/inline-graphic-12.gif [19]: /embed/graphic-7.gif [20]: /embed/inline-graphic-13.gif [21]: /embed/inline-graphic-14.gif [22]: /embed/inline-graphic-15.gif [23]: /embed/graphic-8.gif [24]: /embed/inline-graphic-16.gif [25]: /embed/inline-graphic-17.gif [26]: /embed/inline-graphic-18.gif [27]: /embed/graphic-9.gif [28]: /embed/inline-graphic-19.gif [29]: /embed/graphic-10.gif [30]: /embed/graphic-11.gif [31]: /embed/inline-graphic-20.gif [32]: /embed/inline-graphic-21.gif [33]: /embed/inline-graphic-22.gif [34]: /embed/graphic-12.gif [35]: /embed/inline-graphic-23.gif [36]: /embed/inline-graphic-24.gif [37]: /embed/inline-graphic-25.gif [38]: /embed/inline-graphic-26.gif [39]: /embed/inline-graphic-27.gif [40]: /embed/inline-graphic-28.gif [41]: /embed/graphic-13.gif [42]: /embed/inline-graphic-29.gif [43]: /embed/graphic-14.gif [44]: /embed/graphic-15.gif [45]: /embed/graphic-16.gif [46]: /embed/graphic-17.gif [47]: /embed/inline-graphic-30.gif [48]: /embed/inline-graphic-31.gif [49]: /embed/inline-graphic-32.gif [50]: /embed/inline-graphic-33.gif [51]: /embed/inline-graphic-34.gif [52]: /embed/inline-graphic-35.gif [53]: /embed/inline-graphic-36.gif [54]: /embed/inline-graphic-37.gif [55]: /embed/inline-graphic-38.gif [56]: /embed/inline-graphic-39.gif [57]: /embed/inline-graphic-40.gif [58]: /embed/inline-graphic-41.gif [59]: /embed/inline-graphic-42.gif [60]: /embed/inline-graphic-43.gif [61]: /embed/inline-graphic-44.gif [62]: /embed/inline-graphic-45.gif [63]: /embed/inline-graphic-46.gif [64]: /embed/graphic-18.gif [65]: /embed/graphic-19.gif [66]: /embed/graphic-23.gif [67]: /embed/inline-graphic-47.gif [68]: /embed/inline-graphic-48.gif [69]: /embed/graphic-33.gif [70]: /embed/graphic-34.gif [71]: /embed/graphic-35.gif [72]: /embed/graphic-36.gif [73]: /embed/graphic-37.gif [74]: /embed/graphic-38.gif [75]: /embed/inline-graphic-49.gif [76]: /embed/graphic-39.gif [77]: /embed/graphic-40.gif [78]: /embed/inline-graphic-50.gif [79]: /embed/graphic-41.gif [80]: /embed/inline-graphic-51.gif [81]: /embed/inline-graphic-52.gif [82]: /embed/graphic-42.gif [83]: /embed/graphic-43.gif [84]: /embed/graphic-44.gif [85]: /embed/graphic-45.gif [86]: /embed/graphic-46.gif [87]: /embed/inline-graphic-53.gif [88]: /embed/graphic-47.gif [89]: /embed/inline-graphic-54.gif [90]: /embed/inline-graphic-55.gif [91]: /embed/inline-graphic-56.gif [92]: /embed/graphic-48.gif [93]: /embed/inline-graphic-57.gif [94]: /embed/graphic-49.gif [95]: /embed/graphic-50.gif [96]: /embed/graphic-51.gif [97]: /embed/graphic-52.gif [98]: /embed/graphic-53.gif [99]: /embed/graphic-54.gif [100]: /embed/graphic-55.gif [101]: /embed/graphic-56.gif