Advertisement

Penalized splines model to estimate time-varying reproduction number for Covid-19 in India: A Bayesian semi-parametric approach

Open AccessPublished:November 02, 2022DOI:https://doi.org/10.1016/j.cegh.2022.101176

      Abstract

      Statistical modelling is pivotal in assessing intensity of a stochastic processes. Novel Corona virus disease demanded proactive measures to understand the severity of disease spread and to plan its control accordingly. We propose estimation of reproduction number as a crucial factor to monitor the random dynamics of Covid-19 in India. In the present paper, semi-parametric regression based on penalized splines embedded under Bayesian formulation is utilised to estimate reproduction number while incorporating effects of underreporting and delay in reporting for the actual number of daily occurrences. Monte Carlo Markov Chain approximations are utilised to perform simulation study and thereby to assess the impact of the reporting probability and misspecification of delay pattern on potential for further substance of the pandemic. For a cycle of reporting on weekly basis, the proposed penalized spline Bayesian framework fits closest to the empirical data drawn for a two-day delay in reporting with approximately half of the actual cases being reported. The present paper is a contribution towards estimation of the true daily reproduction number of Covid-19 incidences in its next generation cycle.

      Keywords

      1. Introduction

      It is difficult to capture the exact evolution and dynamics of any pandemic through deterministic modelling. Appearance of novel pathogen leading to sustained pandemic, therefore calls for accounting of temporal piecewise changes in short time-spans. Existing count of infectees and transmission time of infection are the two-key factors in assessment of growth rate of the infection in a specified population. Transmission and propagation of any novel communicable pathogen is quantified through statistical models which are validated through short-history data curated from the current time-chain of the pandemic. Average infection transmissions generated from a currently active single infected unit in a given specified time t is called reproduction number (Rt). In simple words, Rt represents number of secondary infections that one infected individual can spread further on time t. Rt1 indicates that the epidemic is under control and approaching a disease free state while Rt>1 indicates that the control measures lack efficacy leading to the possibility of endemic transforming into epidemic. During progression of any epidemic the expected rise or decline in the reproduction number explains the time-based renewals. These estimates are therefore vital in future medical-assistance and planning in terms of manpower, medical machinery, medicines and health care infrastructure.
      Spline is a flexible mathematical function that represents changing data through disaggregated polynomials. Splines have been an important key for addressing various mathematical problems in approximation theory and in numerical analysis.
      • Agarwal G.G.
      Splines in statistics.
      Incorporating splines in modelling of Rt as an alternative to other popular mathematical models like SIR
      • Ganyani T.
      • Faes C.
      • Chowell G.
      • Hens N.
      Assessing inference of the basic reproduction number in an SIR model incorporating a growth-scaling parameter.
      and SIERD
      • Korolev I.
      Identification and estimation of the SEIRD epidemic model for COVID-19.
      portray challenges
      • Dolton P.
      The statistical challenges of modelling COVID-19.
      in modelling of Covid-19 data and create essential need to appropriately trace the sample-path of its evolution. Rapid changes in trends on daily incidence of Covid- 19 is modelled using semiparametric regression in conjunction with versatile penalized splines to estimate Rt.
      • Mullah M.A.S.
      • Yan P.
      A semi-parametric mixed model for short-term projection of daily COVID-19 incidence in Canada.
      Renewal equation-based estimates of Rt have been explored under closed population assumption which mirrors epidemic persistence under lockdown.
      • Cauchemez S.
      • Boelle P.Y.
      • Thomas G.
      • Valleron A.J.
      Estimation in real time the efficacy of measures to control emerging communicable diseases.
      This approach is motivated by the fact that Rt is viewed as an autoregressive entity determined completely by its present value and the transmission period only. This idiom is also studied through likelihood approach.
      • White L.F.
      • Wallinga J.
      • Finelli L.
      • et al.
      Estimation of the reproductive number and the serial interval in early phase of the 2009 influenza A/H1N1 pandemic in the USA.
      Rt is also investigated under Bayesian lens for efficacy of control measures of communicable diseases in real time.
      • Wallinga J.
      • Lipsitch M.
      How generation intervals shape the relationship between growth rates and reproductive numbers.
      Underreporting of Rt is undertaken for H1N1 pandemic using Bayesian data augmentation and renewal process.
      • Hens N.
      • Van Ranst M.
      • Aerts M.
      • Robesyn E.
      • Van Damme P.
      • Beutels P.
      Estimating the effective reproduction number for pandemic influenza from notification data made publicly available in real time: a multi-country analysis for influenza A/H1N1v 2009.
      Estimation of Rt based on various delay scenarios and misreported H1N1data is studied in context of Mexico and USA.
      • Azmon A.
      • Faes C.
      • Hens N.
      On the estimation of the reproduction number based on misreported epidemic data.
      During any epidemic outbreak, regional governmental bodies are the sole source of reporting of the related incidence data on infection and mortality. Potential for misreporting exists in such data sets due to uncertainties in the real-time reporting of the test outcome as well as due to the limitations of the capacity of the testing facilities usually when pandemic is close to its peak. These limitations have been classified as delay in reporting and underreporting respectively. Length of diagnosis method, limited testing centres, hospital holidays including its restricted-capacity functioning during weekends are multiple sources which cause delays in reporting, while misreporting may arise due to false negative test and ignorance among the population units during early stages of the pandemic. Distribution of the time-elapsed since onset of infection till its occurrence in the infectee is assumed to be known and is termed as serial (or generation) interval distribution. Use of Bayesian inference for semi-parametric regression model brings flexibility in the estimation process through sequential updation of the new infectee counts. In the present work, we adopt Bayesian semiparametric regression model with splines to estimate Rt along with delay and reporting probability. The obtained Bayes estimates are expected to determine the degree of success of the epidemic control strategies.
      Our paper is organised as follows. Section 2 describes the statistical model for step-wise inclusion of the underreporting active cases and delay in the reporting process. In Section 3, we analyse corona virus disease (COVID-19) incidence data for India, estimate the daily reproduction number and the delay parameter for the duration of the outbreak from 15th March 2020 up to 13th April 2020(https://www.covid19india.org/). We conduct a simulation study to assess and validate the proposed method of Bayes estimation of Rt, after adjusting for the impact of misspecification of delay patterns and different reporting probabilities. Section 4 is attributed to discussion. Section 5 concludes the study and explores avenues for further research.

      2. Methods

      The manuscript comprises of majorly four broad methodological concepts. First, observed cases and actual cases are assumed to follow Poisson distribution. Second, the epidemic renewal equation is incorporated for reflecting the progression of disease spread. Third, semi-parametric spline regression which is used to model Rt incorporating effect of underreporting and delay structures fitted under Bayesian paradigm and different scenarios created are compared through Deviance Information Criterion (DIC). Fourth, the simulated Rt under various scenarios is compared with estimated Rt with Mean average square error (MASE) and its components.

      2.1 Notations and assumptions

      Let At={A1,A2,,AT} denote the actual cases of new disease counts during T days of an epidemic which consists of a reported part St={S1,S2,,ST} and an unreported part Ct={C1,C2,,CT} such that At=St+Ct for t = 1, 2, 3 … T. During any epidemic, the total count of reported cases (St) are always less than the actual counts (At). This underreporting at time t and delay in reporting are responsible for underestimation of Rt.
      Reporting probability (τ) may be fixed or time dependent and is referred to as thinning parameter. Suppose maximum length of generation interval is r then the distribution of time interval between the infection times of an infected case and its infector is represented by the ordered set g={g1,g2,,gr}and termed as the generation interval probability. Hence, count of individuals infected on day ti is given as giAti for i = 1, 2, …, r.
      Assuming that all the infected individuals on day ti have the same capacity of further transmission (Rt) and the transmission capability Rt changes with each generation, the average infections at time t is denoted by λt=Rt(i=1t1giAti) for t = 2,3, ….,r while average infected in the first generation remains deterministic at λ1. Additionally, dit represents delay probabilities which account for the percentage of cases on day i reported on day t. More specifically, dit captures the delay structure or delayed cases and is the proportion of total cases which are reported on day t but actually belong to day i or in other words if total number of reported or observed cases are 100 on day t and dit is 0.4 then 40 cases actually belong to day i and are reported on day t. Observed count of infectees, including delaying and underreporting is denoted by Ht={H1,H2,,HT}. Due to delay structures, Ht would exceed At for some t. For constant value of τ and dit over time T, shape of epidemic curve remains same accompanied by its positional shift only. For dynamic τ and dit over a time period, epidemic curve experiences change in its shape. Each of the two types of cases, described observed and actual are assumed to follow Poisson distribution with different means. Observed cases which comprise of underreporting and delay structures ultimately are of great analytical importance for this study and are considered for estimation and inferential purposes.

      2.2 Building distributional structures

      Since Poisson process is a renewal counting process, therefore, total count of infected individuals on day t is assumed to follow Poisson law with parameter λt. Probability of observing At(t>r) conditional on the past prevalence Atr=(At1,At2,.,Atr) is given as,
      P(At/Atr,g,Rt)=exp(λt)λtAt(At)!,
      (1)


      More specifically Atr represents the process prior to the current renewal point. Probability of reported cases is assumed to follow Binomial law, while accounting for data augmentation of actual prevalence due to underreporting is expressed as,
      P(St/At,τ)=(AtSt)τSt(1τ)AtSt
      (2)


      Thus, the observed effective mean prevalence reduces to λtτt without disturbing the existing correlation dynamics in the chain of the infectee and the infected as follows,
      P(St/Atr,Rt,τ,g)=exp(λtτ)(λtτ)St(St)!,
      (3)


      Thus, (3) includes the impact of underreporting in the data. Further, incorporating impact caused by dit, in (3), we have
      P(Ht/λtτ,dit)=exp(ηt)(ηt)Ht(Ht)!
      (4)


      where, ηt=i=1t(λiτ)dit, represents the average number of infectees or the mean number of observed cases (Ht). Also, it is the ηt which includes the three components as reproduction number, underreporting parameter and delay parameter, which are to be estimated.
      Symbolically, the effective mean due to underreporting is weighed by the accumulated cases in the time interval [i,t] being reported at time t. Physically, often the cases from temporal underreporting and from the delay in reporting are not distinguishable. Hence, to resolve such identifiability
      • Azmon A.
      • Faes C.
      • Hens N.
      On the estimation of the reproduction number based on misreported epidemic data.
      ,
      • Thompson R.
      • Baker R.J.
      Composite link functions in generalized linear models.
      usage of composite link function L is made to map transition from {λt} to {λtτ}. Next, we describe the different delay structures as follows.
      • (i)
        One-day delay: A fraction of the new cases on day t is reported on day t + 1, denoted by dtt+1
      • (ii)
        Two-day delay: Fraction of new cases on day t is respectively reported partially on days t +1 and t +2 respectively.
      • (iii)
        Weekend Delay: Due to weekly off on Saturday and Sunday at the reporting health centres, no reporting is recorded on these days. Wednesday, Thursday and Friday will have the same delay pattern as in case (i) above. Monday and Tuesday will have reporting of additional fraction which is carried over from the preceding Saturday and Sunday.
      We use a penalized spline model to allow Rt to change over time t. Furthermore, it is well known that penalized splines can be embedded in the linear mixed model framework where the selection of the smoothing parameter is provided by estimating the random-effects variance component.
      • Crainiceanu C.M.
      • Ruppert D.
      • Wand M.P.
      Bayesian analysis for penalized spline regression using WinBUGS.
      ,
      • Ruppert D.
      • Wand M.P.
      • Carrol R.J.
      Semiparametric Regression.
      We model stochastic Rt in the epidemic renewal equation (5) by using a penalized spline as,
      P(Ht/ηt)Poisson(ηt)


      ηt=i=1t(λiτ)dit=τi=1tditRi(s=1kgsAis)
      (5)


      log(Rt)=α0+α1t+j=1Euj|tεj|3
      (6)


      where θ=(α0,α1,u1,..,uE) is the vector of regression coefficients and εj:E is an ordered set of fixed knots such that ε1<ε2<..<εE and j = 1,2,3 …,E. Spline component in (6) represents the non-parametric part, which is amenable to mathematical computations. The statistical description of the considered model and its estimation through computational software is well documented
      • Crainiceanu C.M.
      • Ruppert D.
      • Wand M.P.
      Bayesian analysis for penalized spline regression using WinBUGS.
      in literature and application of the adapted model is being utilised well for H1N1 pandemic in USA.
      • Azmon A.
      • Faes C.
      • Hens N.
      On the estimation of the reproduction number based on misreported epidemic data.
      To ensure the desired flexibility, the number of knots E should be sufficiently large enough. The choice for number of knots is fixed to 15 and is based on past implementations
      • Azmon A.
      • Faes C.
      • Hens N.
      On the estimation of the reproduction number based on misreported epidemic data.
      ,
      • Crainiceanu C.M.
      • Ruppert D.
      • Wand M.P.
      Bayesian analysis for penalized spline regression using WinBUGS.
      of similar model. The vector of random coefficients U={u1,u2,,uE} is assumed to be independent and normally distributed with E(U)=0 and cov(U)=σu2ΛE1 where ΛE is a matrix in which the (m,n)th entry is |εmεn|3. The class of splines used in the present research is well formulated and decoded
      • Crainiceanu C.M.
      • Ruppert D.
      • Wand M.P.
      Bayesian analysis for penalized spline regression using WinBUGS.
      in past studies. We are modelling log(Rt) instead of Rt because Rt cannot be negative and can take any value between 0 to infinity. The right hand side of equation (6) can be negative after estimation of model parameters. Hence, log(Rt) instead of Rt.
      We assign diffused informative priors as α0,α1Normal(0,104); σu2Gamma(105,105) and postulating flat prior Beata(1,1) for τ and dij each. The choice of priors has been considered after understanding the application of the considered model under Bayesian setup in past studies
      • Azmon A.
      • Faes C.
      • Hens N.
      On the estimation of the reproduction number based on misreported epidemic data.
      ,
      • Crainiceanu C.M.
      • Ruppert D.
      • Wand M.P.
      Bayesian analysis for penalized spline regression using WinBUGS.

      3. Results

      Considering, E = 3 with fixed generation distribution g={0.2,0.3,0.5} we estimate time varying Rt, with fixed τt for t = 1,…,30 days. Smoothed estimate of Rt thus obtained are plotted in Fig. 1. We use R version 4.1.0 and OpenBUGS to execute the proposed theoretic formulations for estimation of dij and Rt and used Deviance Information Criterion
      • Spiegelhalter D.J.
      • Best N.G.
      • Carlin B.P.
      • Van Der Linde A.
      Bayesian measures of model complexity and fit.
      (DIC) to assess model adequacy and complexity (Table 1, Table 2, Table 3) under fixed τ=0.1,0.5,0.9. Posterior mean and posterior standard deviation are obtained based on Markov Chain Monte Carlo (MCMC) approximations. We run MCMC simulations for 105 iterations and attain convergence with respect to MCMC error (<0.05). Each simulation is repeated 1000 times to obtain the corresponding posterior estimate. First, we generate epidemic data on the basis of a “reference”, which we will refer to as the true underlying delay pattern, under τ = 0.5. Next, we estimate Rt,t=1,2,3,,30, on the basis of four other delay patterns. Five combinations or scenarios of underreporting and delay structures are considered.
      • (i)
        Underreporting, one-day-delay (DUO).
      • (ii)
        Underreporting, two day-delay (DUT).
      • (iii)
        Underreporting weekend-delay (DUW).
      • (iv)
        Underreporting, no-delay (DUN).
      • (v)
        No-underreporting with no-delay (DNN).
      Fig. 1
      Fig. 1Smoothed mean estimates of the reproduction numbers for the one-day, two-day and weekend delay patterns with reporting probability of τ = 0.5.
      Table 1Posterior Summary for delay parameters under the case DUO.
      Delay PatternOne-Day Delay
      Reporting Probabilityτ = 0.1τ = 0.5τ = 0.9
      Delay ParametersMeanSDMeanSDMeanSD
      dMTu0.2193240.1478980.1660860.1219790.1513880.097722
      dTuW0.2988630.1455660.2597860.1227430.2660190.106074
      dWTh0.2148090.143090.2043690.1185020.2139840.103941
      dThF0.1247430.1397310.1304010.1155150.1314810.098246
      dFSa0.1611930.1406940.1648480.1154840.152240.095382
      dSaSu0.2681780.1417420.2570180.1155370.234050.093469
      dSuM0.1599340.146530.1221530.1185080.0968410.092972
      τ0.8631450.1201710.8684460.1097870.8914940.088073
      DIC = 402.9DIC= 396.5DIC=530.8
      Table 2Posterior Summary for the delay parameters under the case DUT.
      Delay PatternTwo-Day Delay
      Reporting Probabilityτ = 0.1τ = 0.5τ = 0.9
      ParametersMeanSDMeanSDMeanSD
      dMTu0.8222160.0779090.8516640.059090.3640460.197813
      dMW0.1038380.0786860.0750730.0590530.4578520.190282
      dTuW0.8506030.1014630.8463670.0985350.3693430.20879
      dTuTh0.6229340.1021480.5644620.1008540.2724690.153726
      dWTh0.2617710.183070.312790.1960050.4358110.237962
      dWF0.5920060.1772130.5447390.1799810.3177140.219379
      dThF0.5315260.1883380.5550020.196510.5817970.210321
      dThSa0.1576240.1257820.1612490.1233210.0659320.062097
      dFSa0.8679980.1071570.8578390.124850.1090260.095322
      dFSu0.1424620.1033170.1560390.1081980.7995820.149343
      dSaSu0.4026180.1140680.3541120.1127020.1092460.088297
      dSaM0.6897880.0516350.7205710.050940.2009290.112594
      dSuM0.0300940.0284550.0284530.0267090.5754380.128044
      dSuTu0.619710.0628310.5652320.0595120.1686540.107355
      τ0.8578950.0937550.8847830.0666890.877760.063583
      DIC = 388.6DIC = 359.0DIC = 347.5
      Table 3Posterior Summary for the delay parameters under the case DUW.
      Delay PatternWeek-Day Delay
      Reporting Probabilityτ = 0.1τ = 0.5τ = 0.9
      ParametersMeanSDMeanSDMeanSD
      dMTu0.8125550.0603030.7934820.0632610.8252310.057408
      dTuW0.9481640.0392350.9415620.0409130.9725140.023348
      dWTh0.5546710.0486580.5923350.0468730.6236210.037133
      dThF0.2389580.0425920.2904530.0469450.2829720.035502
      dFSa0.0397450.0331620.0580160.0426540.0246230.021552
      dSaSu0.0707550.0671460.0916190.0803960.063220.059453
      dSuM0.6463070.104530.6379990.1195410.6027090.092993
      τ0.84350.0936070.8612990.0731370.8763650.070275
      DIC = 6.914e+13DIC = 6.914e+13DIC=6.914e+13
      We perform the next stage of simulation to assess the impact of assuming a wrong delay pattern on the estimation of the Rt,t=1,2,3,,T. We calculate the Mean Average Squared Error (MASE) with N = 1000 replications for T = 30 days to assess model efficacy as under,
      MASE=1Nn=1N[1Tt=1T(RˆtRt)2]


      MASE and its bias-variance components for Rt for each true delay pattern vis-à-vis other considered delay-underreporting combinations (Table 4 and Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6). We conduct sensitivity analysis to examine the impact of varying reporting probability τ= 0.15, 0.2, 0.4 and 0.6. on the Rt estimates (Table 5). Estimates of τ are plotted under one-day, two-day and weekend delay misspecifications in Fig. 7, Fig. 8, Fig. 9.
      Table 4MASE and its bias–variance decomposition for the estimates of the reproduction numbers based on the sensitivity analysis of different delay patterns.
      True

      Delay

      Pattern
      Fitted Pattern
      DUODUTDUWDUNDNN
      DUOMASE3.96E-050.0021950.0005973489.81E-059.78E-05
      Var1.3882690.042434.0848021.6392462.187404
      (Bias)27.08E-060.3110710.038451180.0018790910.00971007
      DUTMASE0.0035396.27E-050.0158430.007740840.008196599
      Var3.20956211.204480.1594351.0399511.325875
      (Bias)20.349430.010764762.0021580.62810180.5307091
      DUWMASE0.0067872.31E-050.008088790.0040915650.00454935
      Var4.8847140.0447226.0907932.9661993.347738
      (Bias)21.10E-010.4318511.34E-010.0018080020.006190591
      DUNMASE0.0007710.0017360.0007471936.08E-058.00E-05
      Var4.355360.0379454.3255292.214971.887889
      (Bias)20.036460.2490050.032461080.005959660.003908461
      DNNMASE0.0002460.0015860.0035606945.88E-040.000183114
      Var2.4288340.0305268.7550433.4949222.384256
      (Bias)20.021970.2309760.2029110.025730870.00244879
      Fig. 2
      Fig. 2Estimation of Rt taking DUN as reference scenario.
      Fig. 3
      Fig. 3Estimation of Rt taking DNN as reference scenario.
      Fig. 4
      Fig. 4Estimation of Rt taking DUO as reference scenario.
      Fig. 5
      Fig. 5Estimation of Rt taking DUT as reference scenario.
      Fig. 6
      Fig. 6Estimation of Rt taking DUW as reference scenario.
      Table 5Sensitivity analysis for non-constant reporting probability (τ).
      Delay PatternOne-Day
      Reporting Probabilityτ = 0.15τ = 0.2τ = 0.4τ = 0.6
      MASE1.38E-050.0009798.57E-055.01E-05
      Var1.5592311.9961152.1968161.856806
      (Bias)20.0007460.0648640.0022980.003574
      Delay PatternTwo-Day
      Reporting Probabilityτ = 0.15τ = 0.2τ = 0.4τ = 0.6
      MASE0.0002185.89E-053.54E-053.24E-05
      Var8.7369629.8103910.380212.301
      (Bias)20.0153070.0014467.41E-060.000477
      Delay PatternWeekend
      Reporting Probabilityτ = 0.15τ = 0.2τ = 0.4τ = 0.6
      MASE0.0092350.0080480.0073510.007217
      Var7.7889366.5511355.8166345.847377
      (Bias)21.0671531.0942361.101530.971853
      Fig. 7
      Fig. 7Estimated Rt under different reporting probabilities.
      Fig. 8
      Fig. 8Estimated Rt under different reporting probabilities.
      Fig. 9
      Fig. 9Estimated Rt under different reporting probabilities.
      Table 1 exhibited overall decrease in estimated one-day-delay probability for increasing reporting probability. DIC shows best, when reporting probability is half the actual occurrence (i.e., τ = 0.5)
      Table 2 displayed mixed trend of estimated delay in probabilities due to increase in reporting probability under two day set-up. Again, DIC is observed to be lowest when number of reported cases are half of the actual cases.
      Estimates of delay in probabilities are same for all the considered reporting probabilities. Overall, among all delay structures considered, the lowest DIC is seen for τ = 0.5 under two-day delay setup for the considered data set.
      Table 4 shows the corresponding MASE and its bias–variance decomposition for estimates of the reproduction numbers based on different delay and reporting patterns. Accounting for delay and underreporting is seen to increase the accuracy for estimates of Rt. Lowest MASE is obtained for the correct delay pattern, which validates the model. Largest MASEs are obtained for misspecification in DUT. In general, Models with underreporting show less bias than those which incorrectly ignore underreporting. MASE's for the two considered scenarios of no-delay structures are closest to each other.
      Table 5 yields the smaller MASE's for high reporting probabilities (0.4 and 0.6) for two-day and weekend delay patterns as compared to one-day delay patterns. Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6 depict the estimated daily reproduction numbers are shown for different reference and fitted delay pattern combinations. Fig. 7, Fig. 8, Fig. 9 show the estimated reproduction numbers under different delay patterns for different values of reporting probability. MASE is seen to be highest for DUW and lowest for DUO for all the reporting probability situations. MASE's are seen to decrease with increase in reporting probability for DUT and DUW.

      4. Discussion

      Estimation of reproduction number, usually referred to as “Ro” in most of the studies, comprises of sufficient number of classical statistical methods utilised in recent past in India. Exponential growth rate model
      • Marimuthu S.
      • Joy M.
      • Malavika B.
      • Nadaraj A.
      • Asirvatham E.S.
      • Jeyaseelan L.
      Modelling of reproduction number for COVID-19 in India and high incidence states.
      in India for Covid-19 has gained decent popularity in estimation reproduction number by using “Ro pacakage”
      • Marimuthu S.
      • Joy M.
      • Malavika B.
      • Nadaraj A.
      • Asirvatham E.S.
      • Jeyaseelan L.
      Modelling of reproduction number for COVID-19 in India and high incidence states.
      ,
      • Mitra A.
      • Pakhare A.P.
      • Roy A.
      • Joshi A.
      Impact of COVID-19 epidemic curtailment strategies in selected Indian states: an analysis by reproduction number and doubling time with incidence modelling.
      and “incidence package”
      • Mitra A.
      • Pakhare A.P.
      • Roy A.
      • Joshi A.
      Impact of COVID-19 epidemic curtailment strategies in selected Indian states: an analysis by reproduction number and doubling time with incidence modelling.
      in R software. Estimates for reproduction number ranged from 1.6 to 2.7 in studies for Covid-19 in India and its states.
      • Marimuthu S.
      • Joy M.
      • Malavika B.
      • Nadaraj A.
      • Asirvatham E.S.
      • Jeyaseelan L.
      Modelling of reproduction number for COVID-19 in India and high incidence states.
      • Mitra A.
      • Pakhare A.P.
      • Roy A.
      • Joshi A.
      Impact of COVID-19 epidemic curtailment strategies in selected Indian states: an analysis by reproduction number and doubling time with incidence modelling.
      • Patrikar S.R.
      • Kotwal A.
      • Bhatti V.K.
      • et al.
      Incubation period and reproduction number for novel coronavirus (COVID-19) infections in India.
      • Rai B.
      • Shukla A.
      • Dwivedi L.K.
      COVID-19 in India: predictions, reproduction number and public health preparedness.
      A series of studies
      • Liu Ying
      • Gayle Albert A.
      • Wilder-Smith Annelies
      • Rocklöv Joacim
      The reproductive number of COVID-19 is higher compared to SARS coronavirus.
      was carried out for COvid-19 in China for estimation of reproduction number using all possible statistical techniques based on maximum likelihood estimation, SEIRD, and MCMC. The estimated reproduction number for Covid-19 ranged from 2 to 7 in most of the studies.
      • Liu Ying
      • Gayle Albert A.
      • Wilder-Smith Annelies
      • Rocklöv Joacim
      The reproductive number of COVID-19 is higher compared to SARS coronavirus.
      Renewal equations and their role in stochastic process
      • Karlin S.
      A First Course in Stochastic Processes.
      is well known in understanding the time-based random phenomenon observed in real life applications. In context of epidemic modelling, renewal theory based
      • Mishra S.
      • Berah T.
      • Mellan T.A.
      • et al.
      On the Derivation of the Renewal Equation from an Age-dependent Branching Process: An Epidemic Modelling Perspective.
      applications are limited in number. Embedding renewal equation into spline regression for capturing the dynamics of pandemic caused by the novel corona virus under different scenarios is the soul of this research. More specifically, the present research was to estimate reproduction number incorporating the effect of underreporting and delay in reporting through innovative approach of penalized spline regression through Bayesian toolkit which isn't being explored yet in context of Covid-19 in India.
      We note that there is a large variability in the estimation of the reproduction number for the first few observations. Often Information is scarce at the beginning of any epidemic break-out, hence Rt curve shows higher probability for the initial recorded cases. It can also imply that the initial infectees are responsible for larger number of secondary cases in absence of either awareness or control measures or both.
      DIC being a measure of model adequacy and complexity gave us the answers for identifying the best case scenario for underreporting and delay parameters. Evidence from the analysis (Table 2) shows that among all the considered delay structures for different reporting probabilities two-day delay pattern with a reporting probability of 50% was the most suitable scenario which follows the considered dataset closely. Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6 show that mis-specifying the delay pattern for DUO, DUT, DUW, DUN and DNN has a moderate impact on the estimated trend for the reproduction number. Note that the structure of the DUT pattern is an extension of the DUO pattern. There is a substantial impact when misspecifying the DUW pattern. MASE with the help of simulation gave us the answers on how model estimates vary on mis-specifying a particular scenario and reflected the stability and suitability of the adopted modelling procedure. Simulation study (Table 4, Table 5) justified the modelling strategy adopted for the stated problem under different cases of misspecification of underreporting and delay structures.

      5. Conclusion

      Reporting of Rt instead of daily cases, alerts both the government and public to start and regulate preventive measures accordingly. In the present paper, explanation of the logarithmic transform of Rt through penalized splines under Bayesian setup have generated stronger evidence in favour of the proposed model as is evident from following discussion. The present paper follows Bayesian paradigm under diffused informative priors by treating the reporting fractions, reposting delays and their interactive influences distinctly. Time dependent Rt are seen to be influenced by reporting fractions and systematic delays in reporting. The present study conclude that observed data accompanied by 50% underreporting for a two day lag in recording the incidence of covid-19 cases are closest to the empirical data. A precise estimate of Rt, therefore, is crucial in formulation of corrective or preventive policies towards handling the epidemic and towards efficacy of intervention indicators such as lockdown, quarantine and other related control measures. Future scope of study could include a stochastic serial distribution, correlated reproduction numbers and a different delay pattern than the weekly logistics studied in the present paper. Also, propagation of pandemic is considered for closed population in the present work, which could be deviated from, to understand the advantage and necessity of imposition of movement across spatial contours. Such a study could be useful in understanding cost-benefit bargains in economic activities. Rt can be computed geographically i.e., for each state, district or any other administrative boundary to understand the similarities and differences in dynamics of pandemic like Covid-19 for geospatial analysis and local level covariates responsible for pandemic situation.

      Sources of funding

      The first author gratefully acknowledges IoE grant from University of Delhi.

      Declaration of competing interest

      All the authors declared to have no conflict of interest.

      Acknowledgements

      Authors express their gratitude to all the anonymous reviwers and the editor for their constructive comments and suggestions which resulted in immense improvement in the original manuscript.

      References

        • Agarwal G.G.
        Splines in statistics.
        Bull Allahabad Math Soc. 1989; 4: 1-5
        • Ganyani T.
        • Faes C.
        • Chowell G.
        • Hens N.
        Assessing inference of the basic reproduction number in an SIR model incorporating a growth-scaling parameter.
        Stat Med. 2018; : 1-17
        • Korolev I.
        Identification and estimation of the SEIRD epidemic model for COVID-19.
        J Econom. 2021; 220: 63-85
        • Dolton P.
        The statistical challenges of modelling COVID-19.
        Natl Inst Econ Rev. 2021; 257: 46-82https://doi.org/10.1017/nie.2021.22
        • Mullah M.A.S.
        • Yan P.
        A semi-parametric mixed model for short-term projection of daily COVID-19 incidence in Canada.
        Epidemics. 2022; 38100537
        • Cauchemez S.
        • Boelle P.Y.
        • Thomas G.
        • Valleron A.J.
        Estimation in real time the efficacy of measures to control emerging communicable diseases.
        Am J Epidemiol. 2006; 164: 591-597
        • White L.F.
        • Wallinga J.
        • Finelli L.
        • et al.
        Estimation of the reproductive number and the serial interval in early phase of the 2009 influenza A/H1N1 pandemic in the USA.
        Influ Other Respir Virus. 2009; 3: 267-276
        • Wallinga J.
        • Lipsitch M.
        How generation intervals shape the relationship between growth rates and reproductive numbers.
        Royal Soc. 2007; 274: 599-604
        • Hens N.
        • Van Ranst M.
        • Aerts M.
        • Robesyn E.
        • Van Damme P.
        • Beutels P.
        Estimating the effective reproduction number for pandemic influenza from notification data made publicly available in real time: a multi-country analysis for influenza A/H1N1v 2009.
        Vaccine. 2011; 29: 896-904
        • Azmon A.
        • Faes C.
        • Hens N.
        On the estimation of the reproduction number based on misreported epidemic data.
        Stat Med. 2014; 33: 1176-1192
        • Thompson R.
        • Baker R.J.
        Composite link functions in generalized linear models.
        Appl Stat. 1981; 30: 125-131
        • Crainiceanu C.M.
        • Ruppert D.
        • Wand M.P.
        Bayesian analysis for penalized spline regression using WinBUGS.
        J Stat Software. 2005; 14: 1-24
        • Ruppert D.
        • Wand M.P.
        • Carrol R.J.
        Semiparametric Regression.
        first ed. Cambridge University Press, Cambridge2003
        • Spiegelhalter D.J.
        • Best N.G.
        • Carlin B.P.
        • Van Der Linde A.
        Bayesian measures of model complexity and fit.
        J R Stat Soc B Stat Meth. 2002; 64: 583-639
        • Marimuthu S.
        • Joy M.
        • Malavika B.
        • Nadaraj A.
        • Asirvatham E.S.
        • Jeyaseelan L.
        Modelling of reproduction number for COVID-19 in India and high incidence states.
        Clin Epiderm Global Health. 2021; 9: 57-61
        • Mitra A.
        • Pakhare A.P.
        • Roy A.
        • Joshi A.
        Impact of COVID-19 epidemic curtailment strategies in selected Indian states: an analysis by reproduction number and doubling time with incidence modelling.
        PLoS One. 2020; 15e0239026https://doi.org/10.1371/journal.pone.0239026
        • Patrikar S.R.
        • Kotwal A.
        • Bhatti V.K.
        • et al.
        Incubation period and reproduction number for novel coronavirus (COVID-19) infections in India.
        medRxiv. 2020;
        • Rai B.
        • Shukla A.
        • Dwivedi L.K.
        COVID-19 in India: predictions, reproduction number and public health preparedness.
        medRxiv. 2020;
        • Liu Ying
        • Gayle Albert A.
        • Wilder-Smith Annelies
        • Rocklöv Joacim
        The reproductive number of COVID-19 is higher compared to SARS coronavirus.
        J Trav Med. March 2020; 27 (taaa021,)https://doi.org/10.1093/jtm/taaa021
        • Karlin S.
        A First Course in Stochastic Processes.
        Academic press, 2014
        • Mishra S.
        • Berah T.
        • Mellan T.A.
        • et al.
        On the Derivation of the Renewal Equation from an Age-dependent Branching Process: An Epidemic Modelling Perspective.
        2020 (arXiv preprint arXiv:2006.16487)