Framework for enhancing the estimation of model parameters for data with a high level of uncertainty

Reliable data is essential to obtain adequate simulations for forecasting the dynamics of epidemics. In this context, several political, economic, and social factors may cause inconsistencies in the reported data, which reflect the capacity for realistic simulations and predictions. In the case of COVID-19, for example, such uncertainties are mainly motivated by large-scale un- derreporting of cases due to reduced testing capacity in some locations. In order to mitigate the effects of noise in the data used to estimate parameters of models, we propose strategies capable of improving the ability to predict the spread of the diseases. Using a compartmen- tal model in a COVID-19 study case, we show that the regularization of data by means of Gaussian Process Regression can reduce the variability of successive fore- casts, improving predictive ability. We also present the advantages of adopting parameters of compartmental models that vary over time, in detriment to the usual approach with constant values.

Gaussian Process Regression· Regulariza- tion of data· Uncertainty quantification· Noisy data· Time-dependent parameters 1 Introduction Since the onset of the novel coronavirus (SARS-CoV-2) pandemic, at the end of 2019, a wealth of research has been carried out from across the globe, aiming to un- derstand the dynamics and transmission patterns of the disease. More than a year after the notification of the first case, the number of infected individuals keeps ris- National Laboratory for Scientific Computing, Petrópolis, Brazil E-mail: {glibotte, lanjos, rcca, smcm, rssr}@lncc.br ing significantly worldwide. In the meantime, the con- firmation of reinfections and identification of seasonal immunity [14] reinforces the need for actions to contain the disease, even in locations where the epidemic would be under control. Governmental decisions to mitigate the spread of the disease, such as the introduction of lockdown and social distancing measures, are usually based on com- putational simulations whose preeminent objective is to predict the way the disease spreads in the popula- tion, considering continuously reported data [12]. How- ever, there are several factors associated with natural, economic, and social aspects that make it difficult to adequately predict the spread of the disease and, con- sequently, the definition of a comprehensive policy for prevention and control of the disease [40,16,31]. The heterogeneity of the population in relation to attributes such as demographic diversity, age-dependent charac- teristics, and randomness related to the mobility and interaction of individuals makes it hardly possible to create a model capable of incorporating all these fea- tures together (see, for instance, Refs. [4,7,9,13]). Asymptomatic people also play a significant role in the ongoing pandemic. Oran and Topol [30] presented a comprehensive bibliographic review on the estimation ofasymptomaticcasesofCOVID-19indifferentpartsof the world and concluded that the proportion of asymp- tomatic individuals may vary from 40% to 45% in rela- tion to the total number of reported cases. The scenario of widespread underreporting coupled with a deficient screening and testing capacity leads to significant un- certainty in relation to the reported data of infected individuals. The impact of such uncertainties was ana- lyzed by Ioannidis [17], Li et al. [21], and Wu et al. [51]. Turning attention to the impact of the pandemic in Brazil, a country with a continental dimension, such . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. 2 Gustavo B. Libotte et al. problems tend to become even more evident [23]. So- cioeconomicinequalities[37,48]andculturalfactors[11, 33] have a direct impact on access to information and health services, which translates into high rates of infec- tion and, as a consequence, underreporting cases. Veiga e Silva et al. [49] also analyzed presumptive inconsis- tencies in the data collected and made available by the Ministry of Health in Brazil. They reported that there may be a difference of approximately 41% in the num- ber of deaths caused by COVID-10 related complica- tions. In an attempt to describe the spread of COVID- 19 on different population groups, several models have been proposed by means of integrating typical features related to the disease, such as the quarantine period, lockdown, social distancing, and hospitalization. Mas- sonis et al. [24] bring together several of these models, classifyingthemhierarchicallyinrelationtothenumber of coupled features. In general, even the most complex models, that is, those with supposedly more capacity toassociateknowledgeaboutthespreadingdynamicsof the disease, tend to experience some difficulties in iden- tifying the behavior of noisy data in the long-term, as shown by Alberti and Faranda [1] and Roda et al. [40]. The major motivation of this work is to provide al- ternatives to enhance the capacity of estimating pa- rameters in compartmental models and predicting the spread of COVID-19, taking into account data with a high level of uncertainty, such as the number of new casesreportedinthestateofRiodeJaneirosinceMarch 05, 2020. The data do not have well-defined behavior, in such a way that the variations in subsequent days are caused, in part, by the factors that deepen the dispar- ities in the notifications of cases of infection. In addi- tion,theBraziliangovernmentestimatestheCOVID-19 data considering the daily count in the municipalities— Brazil is made up of 5570 municipalities, of which 92 make up the state of Rio de Janeiro—which are au- tonomous in relation to population testing policy and do not follow a common strategy to prevent the disease. Asinsomelocationsdataarenotreportedonweekends, there are sudden drops in the number of new infections, which afterward cause unexpected increases when data are reported at the beginning of the following week. In this context, aiming at expanding the predictive capacity of compartmental models subject to scenar- ios of great uncertainty, the objectives of this work are to propose the use of strategies capable of reducing the variability of the estimation of model parameters and to discuss the advantages of considering time-de- pendent parameters. We propose that the noisy data set be regularized by means of Gaussian Process Re- gression (GPR). This approach allows a set of data to be smoothed, so as to decrease its noise level without significantly changing its behavior. To confirm this as- sumption, we compared a subgroup of reported data on dead and infected individuals in the state of Rio de Janeiro to simulations produced by the SEIRPD-Q model, whose parameters are estimated using regular- ized data. The results obtained are also compared to the simulations calculated in the usual way, without regularizing the data, in order to quantify the predic- tive capacity in relation to known data and the gain in relation to the parameter estimation approach with unchanged data. Arecentreviewoftheliteratureonthesubjectfound that, to date, few works that incorporate the concepts of Gaussian Process (GP) applied to the epidemiologi- cal modeling of COVID-19 have been published. These works are briefly presented below. Ketu and Mishra [19] proposed the Multi-Task GPR model, aiming to predict the COVID-19 outbreak worldwide. The authors com- pared the results obtained with other regression mod- els, in order to analyze the effectiveness of the proposed method. Zhou and Ji [52] proposed a model for trans- mission dynamics of COVID-19 considering underre- porting of cases (what they called undocumented in- fections) and estimated the time-varying disease rate of transmission using GPR and a Bayesian approach. Arias Velásquez and Mejía Lara [2] demonstrated the correlation between industrial air pollution and infec- tions by COVID-19 before and after the quarantine in Peru, by presenting a classification model using Re- duced-Space GPR. This methodology is used by the same authors in Ref. [2] to report a long-term forecast for COVID-19 in the USA. In turn, Ribeiro et al. [38] compared the predictive capacity of various machine learning regression and statistical models, considering short-term forecasting of COVID-19 cumulative cases in Brazil. As far as we are aware, this is the first time that GPR is employed for regularization of COVID-19 data, which are subsequently analyzed using a compart- mental model (although other methods for regulariza- tion have been applied [21]). In our analysis, we partition the data sets between training and test data, and carry out successive param- eter estimations varying the proportion between these types of data, in order to show the behavior of the ob- tained parameter set. We show that some of these pa- rameterscanbeapproximatedbyfunctionsand,consid- ering this possibility, we analyze the influence of adopt- ing variable parameters over time. We use both a de- terministic approach, in terms of least squares, to es- timate the gain related to the use of regularized data and time-varying parameters, and a Bayesian approach, in order to analyze parameter uncertainties, which are . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 3 propagated to the model to quantify uncertainties on the model outcomes. 2 Materials and methods 2.1 Compartmental model description Following the models proposed by Jia et al. [18] and Volpatto et al. [50], we develop an extension of the Susceptible-Exposed-Infected-Removed (SEIR) model, termed SEIRPD-Q model, which further considers Pos- itively Diagnosed (P) and Dead (D) groups of individ- uals. To account for social distancing measures, we im- plicitly consider a mechanism that isolates individuals from the virus transmission (−Q) rather than defining a specific compartment for individuals in quarantine. A schematic description of the SEIRPD-Q model is pre- sentedinFig.1.Weconsiderapopulationsusceptibleto a viral outbreak, whose rate of transmission per contact is given byβ. At the beginning of a COVID-19 infec- tion, infected individuals pass through a latency period in which they are not capable of transmitting the dis- ease to another individual, becoming infectious only af- ter this stage, even without showing any sign of disease. Individuals in such conditions are said to be exposed, with incubation period given by1/˜σ. The infectious in- dividuals are then divided into infected and positively diagnosed compartments, based on the premise that a large number of individuals who contract the virus are not diagnosed. This is due to the reduced testing ca- pacity in some locations so that the diagnosis is made, primarily, in individuals who have severe symptoms or are hospitalized. The proportion of infected individu- als, given by ρ, is related to the majority of people that only suffer mild symptoms and get recovered with- out significant complications. On the other hand, the complement of this group are those individuals who, in fact, have been positively diagnosed. In turn, indi- viduals who recover from the disease after being in the Susceptible Exposed Infected Removed Positively diagnosed Dead /uni0303(1 /uni2212) /uni0303 /uni2212 Fig. 1 Schematic description of the SEIRPD-Q model. infected compartment are moved to the removed com- partment at a rate ofγI. The same goes for positively diagnosed individuals, who are removed at a rate ofγP. In addition, it is reasonable to assume that most of the individuals who died from complications caused by the disease had severe symptoms and were tested or hos- pitalized. Therefore, we do not consider that individu- als in the infected compartment die from virus-related causes without being diagnosed, and the mortality rate of positively diagnosed individuals is given bydP. Quarantine measures are also incorporated in this model, affecting the susceptible, exposed, and infected compartments. Individuals in these compartments are kept in quarantine at a rate ofω, and are not assumed to be infectious considering restrictive quarantine mea- sures. The quarantine compartment is implicitly mod- eled and, therefore, the removed compartment includes individuals who have undergone quarantine, along with those who have been infected and have recovered from the disease. The formulation of the SEIRPD-Q model is given by the following system of ordinary differential equations: dS dt =−βSI−ωS dE dt =βSI− ˜σE−ωE dI dt = ˜σρE−γII−ωI dP dt = ˜σ (1−ρ)E−dPP−γPP dR dt =γII +γPP +ω (S +E +I) dD dt =dPP (1) The SEIRPD-Q model presents some fundamental differences in relation to those on which it was based: first, Jia et al. [18] consider that only susceptible in- dividuals are subject to quarantine measures, which is modeledusinganexplicitcompartmentandalsoconsid- ering an additional parameter that controls the social distancing relaxation; second, we disregard the asymp- tomatic compartment, according to Volpatto et al. [50], due to the lack of this information associated with lim- ited testing of the population; third, we consider only the mortality rate of positively diagnosed individuals, unlike Volpatto et al. [50]. 2.2 Data The data used in this work are the daily number of infected (positively diagnosed) and dead individuals in the state of Rio de Janeiro. The Brazilian Ministry of . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 4 Gustavo B. Libotte et al. Health reports the data daily, which are synthesized and made available as shown in Ref. [10]. The analyzed data refer to the period between March 10, 2020 and October5,2020,consistingof210records.Althoughthe number of recovered individuals is also available, it may be wise not to use this data to estimate the parameters of the model since, due to the unseemly policy of testing the population, there is much uncertainty about these data. In the following, the observable quantities are de- noted by the time seriesD(j) i =D(j) (ti), fori = 1, ..., p and j ∈ {P, D}. Overall, letD = {D1, ..., Dp} be the finite set of real-valued measurements collected from successive observations of the daily number of pos- itivelydiagnosed anddeadindividualsatdifferenttimes ti. 2.3 Deterministic approach for parameter estimation Here we want to solve the so-calledinverse problem, i.e., to determine the model parameters that generate out- puts as close as possible to the observable data. To this end, letθ be them-dimensional vector of model param- eters, andy(j) i =y(j) (ti, θ) be the model responses at different timesti,i = 1, ..., p andj∈{P, D}, as used previously. Additionally, letY (j) i denote the cumulative sum of all model responses given byy(j) i . The objective of the inverse problem is, therefore, to find the vectorˆθ (an estimate ofθ) that produces out- puts ˆy(j) i capable of fitting the available observations. The best fit between the responses of the modelˆy∈ Rp and the observed data can be estimated in terms of the residuals,thedifferencebetweenobservedandpredicted measurements, given byr (θ) = ˆy−D. The solution to an inverse problem is, in other words, the data fitting whose objective is to calculate an estimateˆθ that mini- mizes some error norm∥r∥ of the residuals [45,39]. The least squares fitting calculates the vector of optimal pa- rameters by taking the mean squared error, given by E (θ) = 1 p r (θ)T r (θ) = 1 p p∑ i = 1 ri (θ)2 , (2) and the estimate ˆθ is the vector that minimizes this quantity: ˆθ = arg min θ E (θ) . Equation (2) is usually called theobjective function (or cost function). WhenE→ 0, the estimate ˆθ generates an output vector ˆy that has a high level of agreement withtheobserveddata D,thatis,theresidualsaremin- imized. In general, real problems do not admitE = 0, since the noise that affects the model cannot be pre- dicted with such accuracy. 2.4 Bayesian approach for parameter estimation Bayesian inference provides another perspective for es- timating the value of a set of parameters that best char- acterizes the output of a model, given a set of data. Bayesian inference differs from the deterministic ap- proach because in addition to calibrate the parameter values, it measures their uncertainties, which is one of the focuses of this work. To conduct Bayesian infer- ence, we need some familiarity with a few basic con- cepts of probability. Here, we give a brief overview of such concepts. A more detailed description is provided by Refs. [47,22,3]. Bayes theorem provides a formulation to estimate the posterior probability of the model parameters given a set of observationsD, based on the likelihood of the event of interest occurring given the prior knowledge on the parameters. The theorem is stated as ppost (θ|D) = plike (D| θ)pprior (θ) pevid (D) , (3) whereplike (D| θ) is thelikelihood function,pprior (θ) is the prior informationorbeliefson θ,pevid (D)isthe ev- idence related to the observationsD, andppost (θ|D) is the posterior distribution associated withθ. Apriorknowledgecanbethoughtofastheprobabil- ity density function over the feasible values of the model parameters, the current knowledge on their values. In turn, the likelihood assumes the role of estimating the probability of characterizing the available data, given a set of parameters. In other words, the likelihood func- tion measures how good the data are being explained by the model. In this work we assume a Gaussian like- lihood, which has the form plike (D| θ) = ∏ j 1 σj √ 2π exp  − p∑ i = 1 ( D(j) i −y(j) i )2 2σ2 j  , in which σ2 j is a measure of the uncertainty (encom- passing data errors) related to each quantityj, forj∈ {P, D}. Here, bothσ2 P and σ2 D are considered hyper- parameters to be estimated together with the set of pa- rameters θ. The evidence, also referred to as marginal likelihood, is the integral of the likelihood over the prior and is considered as a normalization constant. Thus, we actually evaluate ppost (θ|D) ∝ plike (D| θ)pprior (θ) to produce the posterior distributionppost (θ|D) over . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 5 the parameters, that is, an updated belief about θ, givenD. WeconductedtheBayesianinferenceusingtheprob- abilistic programming library PyMC3 [41]. We adopt the Transitional Markov Chain Monte Carlo method proposed by Ching and Chen [8] for parameter esti- mation (which is not described here, for the sake of brevity). 2.5 Gaussian Process Regression Here we describe the general GP framework, with spe- cial attention to regression problems using noisy obser- vations. By definition, a GP is a collection of random variables, any finite number of which have a joint Gaus- sian distribution [35]. In other words, it is an extension of the multivariate Gaussian distributions to infinite di- mensionality. GPRs take place directly in the space of functions, defining priors over functions that, once we have seen some data, can be converted into posteriors over functions [27]. Thus, a GPR model is a Bayesian nonlinear regression model that takes into account the GP prior and whose posterior is the desired regression function that belongs to an infinite dimension random function space [43]. To introduce the GP, denote byt = (t1, ..., t p)T the time training points associated to a finite set of p observationsD = (D 1, ..., Dp)T. We assume that each observationD at locationt is a random variable associated to the GP stochastic process f (t)∼GP (m (t), k(t, t′)) , whichiscompletelydefinedbyits mean functionm (t) = E [f (t)], the expected value of all functions in the distri- bution evaluated for an arbitrary inputt, andk (t, t′), the covariance function of f (t), which describes the dependence between the function values for a pair of arbitrary input time pointst andt′, given byk (t, t′) = E [(f (t)−m (t)) (f (t′)−m (t′))]. We now consider the regression problem Di =f (ti) +ε (4) withεbeing an additive Gaussiannoise with zero mean and variance σ2. The GPR begins by assuming the vector-value function f∼N (0, K) as the prior distri- bution, where K is thep×p covariance matrix whose entries arek (t, t′). Considering this prior and noise in the time training set, as defined in Eq. (4), the joint distribution taking into account new input time points t∗ and their associated outputD∗ is given by [ D D∗ ] ∼N ( 0, [ K (t, t) +σ2I K (t, t∗) K (t∗, t) K (t∗, t∗) ]) , where I stands for thep×p identity matrix. Therefore, the predictive equations for GPR are derived from the conditional distribution property for the multivariate Gaussian distribution [35]. Considering the Schur com- plement (for more details on Schur complements, refer to Puntanen and Styan [34]), the posterior predictive distribution is the multivariate Gaussian distribution p (D∗| t∗, t,D) =N (µ∗, Σ∗) , (5) with mean µ∗ = K (t∗, t) ( K (t, t) +σ2I )−1 D (6) and covariance matrix Σ∗ = K (t∗, t∗)−K (t∗, t) ( K (t, t)+σ2I )−1 K (t, t∗) . (7) Therefore, the calculation of Eqs. (6) and (7) is suffi- cient to predictD∗. Note that this involves first calcu- lating the four covariance matrices. Furthermore, the covariance depends only on the time training set (t) and the new input points(t∗), and not on the observa- tion measures vector(D). By aplying the GPR to the observation setD(j), forj∈{P, D}, we then obtain p ( D(j)∗ ⏐⏐ t∗, t,D(j) ) , whose corresponding mean val- ues at the time training points are the regularized data used in this work. The covariance function is commonly called theker- nel of the GP. This function maps a pair of general input vectors t, t′ ∈ t into R. The idea behind the kernel is that ift and t′ are said to be similar, it is ex- pected that the function output (observations) at these points will also be similar. The main attribute of the kernel is to avoid the computation of an explicit non- linear mapping function that relates input and output data, obtaining the identification of the mapping in the space where the number of parameters to be optimized, the so-calledhyperparameters, is smaller [20]. Thus, the choice of an appropriate kernel is usually based on prior knowledge related to the behavior of the training data, as for example the occurrence of periodic oscillations, and assumptions such as smoothness [42]. Finding suit- able properties for the kernel function is one of the main tasks for defining an appropriate GP. The kernel can be any function that relates two in- put vectors, on the assumption that it can be formu- lated as an inner product, producing a positive semi- definite matrix [5]. Different functions can be combined to produce kernels with a variety of features, generally bybothaddingandmultiplyingkernels.Rasmussenand Williams [35] present a comprehensive analysis on the construction of kernel functions. Considering the ability . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 6 Gustavo B. Libotte et al. of GPR to approximate the behavior of the used data, we adopt the radial basis function (RBF) k (t, t′) = kRBF (t, t′), which is explicitly defined as kRBF (t, t′) = exp ( −|t−t′|2 2ℓ2 ) , where ℓ is the length scale of the kernel. Rasmussen and Williams [35] present an algorithm for GPR employing Cholesky factorization to solve the matrix inversion required by Eqs. (6) and (7). We use the Scikit-learn library [32] to implement the GPR. Hyperparameters are calculated using a computational routine internal to the library, which uses the L-BFGS- B algorithm [29] to obtain the optimal values. The op- timizer is restarted 50 times in order to increase the chances of convergence to the optimal set of hyperpa- rameters. In fact, the number of restarts is an arbitrary choice and, as the GPR is executed only once, the com- putational cost associated with this procedure is negli- gible. 2.6 Parameter and numerical experiment setups In models with more compartments and which, in gen- eral, have more parameters, finding a unique set of pa- rametervaluesthatbestfitssomedatamaybeunattain- able. Different combinations of parameters can produce similar results for data fitting. This fact is a charac- teristic of the non-identifiability of parameters [36]. To overcome this issue, some of the biologically-known pa- rameters are fixed, namely the incubation period˜σ = 1/5.8 day−1 [25], the proportion of symptomatic in- fected individualsρ = 0.6 [30], and both the recovery rate of infected and positively diagnosed individuals, defined asγI =γP = 1/16.7 day−1 [46]. Therefore, the parameters to be estimated are the rate of transmission β, rate of removal due to quarantine measuresω, and mortality ratedP. In order to define the initial conditions that make it possible to solve Eq. (1), first consider that the pop- ulation of the state of Rio de Janeiro is approximately equal to N = 17264943 individuals, according to the last demographic census conducted by the Brazilian In- stitute of Geography and Statistics [6]. Using the re- ported data, it is possible to define the initial condi- tions for the number of infected and positively diag- nosed individuals, which we assume to be the same at the beginning of the time series. The number of dead individuals on the first day considered in this analy- sis is also sufficient to define the initial condition for D. In addition, it is reasonable to assumeR (0) = 0, since it is not expected to have recovered individuals at the outbreak of COVID-19. Therefore, the only ini- tial condition for which there is no information from the reported data is related to the exposed individ- uals. Therefore, we assume E (0) as a parameter to be estimated and the initial condition for the num- ber of susceptible individuals is given byS (0) = N− (E (0) +I (0) +P (0) +R (0) +D (0)). We partition data sets into two subsets, which we call training data and test data. The key idea behind thisapproachistoanalyzethepredictivecapacityofthe model, by comparing test data with short-term simula- tions, which are calculated using parameters estimated with training data. In this way, it is possible to compare the gain of using regularized data in the parameter esti- mation procedure, analyzing the results considering an actual scenario. Furthermore, we only consider short- term predictions in our analyzes since, as the simula- tions are compared with original data. Therefore, we analyze scenarios in which we adopt 14 data points in the test set. Of note, other time windows could also be used. Arbitrarily,wechoosetheminimumsizeofthetrain- ing data set equal to 60. We set up 14 values in the test data set and calculate successive estimates of the parameters of the model, gradually increasing the pro- portion between training and test data. After each run, new data is added to the training set, so that the test set is composed of the next 14 values in the time se- ries. As data for 210 days are available, 136 sets of parameters are estimated, using the deterministic ap- proach described in Section 2.3. This procedure is per- formed both using training data as it stands, and after regularization. The simulations using the optimal pa- rameters are compared to the corresponding test set (without being regularized), by the computation of the normalized root-mean-square error (NRMSE) [15], con- sidering both the cumulative number of infected and dead individuals—in this step, the cumulative data is adopted, with the purpose of minimizing the influence of noise. The root-mean-square error is normalized by the difference between the highest and lowest values in the corresponding data set. All optimal parameters are calculated by combin- ing the Differential Evolution [44] and the Nelder-Mead Simplex [28] methods. For each problem, the solution is estimated by Nelder-Mead Simplex and refined by Differential Evolution, which searches for the optimal parameters in the vicinity of the previously obtained point. The solution to each problem—the best indi- vidual in the population—is taken as an initial esti- mate for the next problem. Nelder-Mead Simplex runs with coefficients of reflection, expansion, contraction, and shrinkage equal to 1, 2, 0.5, and 0.5, respectively . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 7 0 2 4 6Infected individuals ×103 Original Regularized 0.0 0.8 1.6 2.4 Cumulative number of cases ×105 19/03/2009/04/2030/04/2021/05/2011/06/2002/07/2023/07/2013/08/2003/09/2024/09/20 0 1 2 3Dead individuals ×102 19/03/2009/04/2030/04/2021/05/2011/06/2002/07/2023/07/2013/08/2003/09/2024/09/20 0.0 0.6 1.2 1.8 Cumulative number of deaths ×104 Fig. 2 Results of the regularization of data using GPR. The blue dots are the original data and the orange dots represent the resulting regularized data. Shaded areas indicate two standard deviations from the corresponding regularized data. We also show the cumulative data, for comparison purpose. (as commonly adopted in the literature). In turn, Dif- ferential Evolution runs with 20 individuals in the pop- ulation, amplification factor equal to 0.6, and crossover probability equal to 0.95. The search space is bounded by 0 ≤ β ≤ 10−6, 0 ≤ ω ≤ 1, 0 ≤ dP ≤ 1, and 0≤E (0)≤ 104 in appropriate units (a. u.). 3 Results and Discussion 3.1 Influence of regularization of data Initially, the data used in this analysis is regularized using GPR as described in Section 2.5. The values of the hyperparameters tunned for the RBF kernel are ℓ = 68.9andℓ = 50.2 fordailydataofinfectedanddead individuals, respectively. Figure 2 shows an illustrative comparison between the original data for daily infected and dead individuals and the corresponding regularized data. Shaded areas represent two standard deviations from the fitting data. We also show the respective cu- mulative data, in order to allow a visual inspection of the agreement of the data resulting from the regular- ization, but smoothing out the noise. Next, the influence of estimating the parameters of the SEIRPD-Q model, presented in Section 2.1, using theregularizeddataset,inrelationtotheapproachthat adopts the original data is analyzed. Our focus is to as- sess the gain related to the regularization of data within the scope of compartmental models and, therefore, we consider the simulations using the model adopted in this analysis. It is worth mentioning that the proposed analysis can be extended to any compartmental model with a structure similar to the model given by Eq. (1). The analysis is performed by evaluating the NRMSE between model predictions and test data.

concerning the optimal parameters calcu- lated by Differential Evolution and Nelder-Mead Sim- plex (Section 2.6) are shown in Fig. 3. For each set of parameters obtained in the 136 runs using the method- ology described in Section 2.6, considering each type of data (original and regularized), the model is simulated and we calculate the corresponding NRMSE. Then we compute the area under the curve, which is composed of the values of the NRMSEs. Analyzing Fig. 3, the ef- 0 25 50 75 100 125 Run 0.0 0.5 1.0 1.5 Normalized root-mean-square error ×102 Original Regularized Fig. 3 NRMSEs computed by comparing simulations of the SEIRPD-Q model, performed using parameters that best fit the training data, where θ = (β, ω, dP, E(0)), in relation to the test set (composed of 14 points). The area under the curve represents the total deviation relative to the test data in all runs, where the number of training data varies. . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 8 Gustavo B. Libotte et al. 0 2 4 6 Daily number of infected individuals ×103 Original data 0.0 0.6 1.2 1.8 ×103 Regularized data 0 1 2 3 Daily number of dead individuals ×102 0.0 0.4 0.8 1.2 ×102 102 104 106 Cumulative number of infected individuals 0.0 0.8 1.6 2.4 ×105 19/03/2009/04/2030/04/2021/05/2011/06/2002/07/2023/07/2013/08/2003/09/2024/09/20 101 104 Cumulative number of dead individuals 19/03/2009/04/2030/04/2021/05/2011/06/2002/07/2023/07/2013/08/2003/09/2024/09/20 0.0 0.5 1.0 1.5 ×104 Fig. 4 Simulations of the daily number of infected and dead individuals using optimal parameters obtained with the de- terministic approach, by varying the amount of data in the training set. The shaded areas represent the variation range of the simulations, whereas the points are the fitted data. Box-and-whisker diagrams are used to show the variability of the simulations on specific days. We also show the corresponding results for cumulative data, where the noise is less effective, for comparison purposes. The results follow the same color scheme as in Fig. 2, blue for original data, and orange for regularized data. fect of data regularization on parameter estimation is clear: for most of the estimated parameter sets, the cor- responding simulations have better agreement with the test data. Translating into numbers, regularized data resulted in more reliable predictions in 64.71% of runs. The area under the curve for original data is approxi- mately equal to 899.66 u. a. (unit of area), whereas for regularized data it is 705.70 u. a. This represents an av- erage improvement of approximately 21.56%. Note that the effect of the regularization is more prominent when the data set to be fitted is larger, since the influence of the noise tends to become more intensive. The implication of using regularized data is even more straightforward when we analyze the variability of simulations resulting from the optimal parameters corresponding to each point in Fig. 3. For this purpose, considertheresultsshowninFig.4,inwhichtheshaded area represents the range of the simulations related to the daily number of infected and dead individuals, for both original and regularized data, whose parameters are estimated by varying the amount of training data, as previously described. In turn, the points represent the corresponding data, which together with the box- and-whisker diagrams, aim to demonstrate the variabil- ity of the simulations on specific days. We also show . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 9 the corresponding cumulative values of the results ob- tained, in order to provide better conditions for com- parison. It is worth mentioning that the set of training data, even presenting variable size considering each set of parameters that is estimated, is completely shown in all cases. It is clear that simulations resulting from parame- ters estimated using data with no regularization have more variability than after regularization. Considering that on October 5, 2020 (the last day of the simulation), Rio de Janeiro had 273,335 confirmed cases, the bound- ary values of the shaded area on this day, for the cu- mulative number of infected individuals obtained with original data are 146,763.4 and 36,518,946.5, whereas the quantities obtained with regularized data vary be- tween 80,762.4 and 253,387.8. In the case of dead indi- viduals,suchvaluesvarybetween8,066.9and3,579,976.6 consideringoriginaldata,andbetween6,950.6and17,528.7 considering regularized data, with the actual cumula- tivenumberofdeadindividualsonthisdaybeing18,780. The dispersion of these results can be alternatively analyzed using a box-and-whisker diagram. The boxes areboundedbythelowerandupperquartiles(whichwe call Q1 and Q3, respectively) of the full model predic- tion values, and the horizontal line inside the box repre- sents the median. Whiskers, the vertical lines bounded by perpendicular dashes, extend from the minimum value of the data to the first quartile, and from the third quartile to the maximum value of the data. Whiskers expressthevariabilityoutsidethesequartiles.Theyran- ge from the lowest to the highest values of each set, both for original and regularized data, since no simula- tion can be considered an outlier. The statistical results associated with the box-and-whisker diagrams in Fig. 4 are summarized in Table 1. It is clear that regulariza- tion using GPR preserves the behavior of the data, al- though substantially reducing the variability of the sim- ulations. The reduction in variability provides evidence that regularizing the data with the GPR does not cause problems related to overfitting. Simulations influenced by overfitting would deviate significantly from the real data. In this case, parameter estimation using regular- ized data provides a means of reducing noise without negatively influencing forecasts [26]. Now consider a forecast of the epidemic in terms of the cumulative number of infected and dead indi- viduals, aiming to show the difference of the simula- tions in relation to the analyzed data. For this purpose, we adopt the Bayesian approach for parameter estima- tion, presented in Section 2.4. We arbitrarily choose the training data sets for the daily number of infected and dead individuals, both original and regularized, with 196 values, the maximum number of elements that the T able 1 Statistical results of all 136 simulations performed using the optimal parameters. The results refer to the cumu- lative number of infected and dead individuals, that is,Y(j) p for j∈{P, D}, obtained using both original and regularized data on the last day that the model was simulated, October 5, 2020. On this day, Rio de Janeiro accumulated 273,335 confirmed cases and 18,780 deaths. Original data Regularized data Y(P ) p Y(D) p Y(P ) p Y(D) p Min 146,763.4 8,066.9 80,762.4 6,950.6 Max 36,518,946.5 3,579,976.6 253,387.8 17,528.7 Median 259,406.5 17,313.5 217,621.5 16,196.7 Q1 172,065.9 14,796.5 163,064.3 15,249.4 Q3 3,970,094.4 342,994.5 243,274.8 17,059.5 test set can contain. In turn, the test set is made up of the next 14 values in the time series. By this analysis, we present a visual perspective of the benefit of using regularizeddataintheparameterestimationprocedure, providing a way of comparing the agreement between simulations and test data. The parameters to be estimated are the same as in the previous analysis, and are assumed to be uni- formly distributed, in such a way thatβ∼U ( 0, 10−6) , ω∼U (0, 1),dP∼U (0, 1), andE (0)∼U ( 0, 104) , in a. u. Table 2 shows themaximum a posterior (MAP) estimates of the parameters, as well as the correspond- ing 95% credible interval (CI) [3]. Likewise, all other parameters of the SEIRPD-Q model are taken as bio- logical parameters and have the same values previously reported (see Section 2.6). Figure 5 shows the poste- rior distribution of the estimated parameters of the SEIRPD-Q model, obtained when original (blue bins) and regularized (orange bins) data are employed for T able 2 MAP values and 95% CIs of the parameters esti- mated using Bayesian calibration (in a. u.). Data type Original Regularized β 1.5248× 10−8 1.2419× 10−8 (1.2795, 1.7991)× 10−8 (1.2206, 1.2639)× 10−8 ω 0.007121 0.004859 (0.005704, 0.008576) (0 .004738, 0.004996) dP 0.07126 0.0625 (0.06215, 0.08161) (0 .06030, 0.06475) E (0) 1190.8021 1236.5740 (655.0367, 1818.3565) (1168.6603, 1305.0816) . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 10 Gustavo B. Libotte et al. 0.0050 0.0075 0.0100 ω 0.060 0.075 0.090 dP 1.2 1.6 2.0 ×10−8 β 800 1600 2400 E(0) 0.45 0.60 0.75 0.90 1.05 ω ×10−2 1.00 1.25 1.50 1.75 2.00 β ×10−8 0.060 0.064 dP 0.0048 0.0050 ω 1200 1300 E(0) 1.20 1.23 1.26 ×10−8 β 6 7 8 9 dP ×10−2 0.5 1.0 1.5 2.0 2.5 E(0) ×103 Fig. 5 The left frame shows the posterior distribution of the parameters obtained in the Bayesian inference for original (in blue) and regularized (in orange) data, by fitting the daily number of infected and dead individuals; the right frame illustrates the variance of each parameter in a comparative way, where the same color scheme is adopted. calibration. In the right frame, violin plots express the variance of the inferred parameters using original and regularized data. Of note, the latter is much narrower than the former. This fact is reflected in the stochas- tic simulation of the SEIRPD-Q model shown in Fig. 6, along with the training and test data sets (note that the latter is never regularized). The solid curves represent the model responses when the free parameters are set to be the respective MAP values, shown in Table 2. In turn, the shading around the curves represents the 95% CI, the discrete points are the available data, and the ranges of the training and test sets are colored in gray for daily and cumulative data, respectively. Results re- lated to infected individuals are shown in red, whereas dead individuals are shown in green. Analyzing Fig. 6, it is clear that the Bayesian in- ference obtained suitable results, using both original and regularized data. The deviation between the simu- lations and the cumulative data in the test data range when regularized data are used for calibration is due to the regularization calculated using GPR, as can be seen in Fig. 2. Nevertheless, the regularization of the training data seems to make the predictive capacity of the model less unstable. Thus, the model could bet- ter capture the behavior of the data, eventually leading to more appropriate predictions. Of note, the narrower posterior distributions obtained by performing the cal- ibration with the regularized data resulted in less un- certainty about the values predicted by the model, as can be seen in Fig. 6. This is because the posterior dis- tributions obtained when data are regularized tend to have a smaller standard deviation, as can be seen in the right frame of Fig. 5 and in the CIs shown in Table 2. 3.2 Influence of time-varying parameters A further feature that may improve the predictive ca- pacityofcompartmentalmodelsistheadoptionoftime- varying parameters. In an effort to analyze this aspect, consider the 136 sets of optimal parameters, obtained using the deterministic approach that led to the results shown in Fig. 3. Figure 7 shows the values of the cali- brated parameters for each corresponding run, both us- ing original and regularized data. Initially, note that the variability of the set of parameters obtained by fitting the original data is much higher than those referring to . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 11 0 2 4 6 Daily number of infected individuals ×103 Original data Training data range 0 2 4 6 ×103 Regularized data 0 1 2 3 Daily number of dead individuals ×102 0 1 2 3 ×102 0 1 2 3 Cumulative number of infected individuals ×105 Test data range 0 1 2 3 ×105 29/03/2019/04/2010/05/2031/05/2021/06/2012/07/2002/08/2023/08/2013/09/2004/10/20 0.0 0.6 1.2 1.8 Cumulative number of dead individuals ×104 29/03/2019/04/2010/05/2031/05/2021/06/2012/07/2002/08/2023/08/2013/09/2004/10/20 0.0 0.6 1.2 1.8 ×104 Fig. 6 Simulations with parameters obtained by fitting the daily number of infected (in red) and dead (in green) individuals using the Bayesian approach. The shaded areas refer to the 95% CI (see Table 2). Hatched areas bound the range of training and test data. Training data is never regularized. regularized data. This is a consequence of the high level of noise in the data. Note, for instance, the gray shaded areas in Fig. 7. They refer to the same runs for all pa- rameters and correspond to the shaded area in the same color in Fig. 2. In this time window, the data show very large variations on subsequent days, especially those of the daily number of infected individuals. These varia- tions influence the values of the parameters obtained, giving rise to a great difference in the value of the set of parameters that best fits the data just by adding a new single value to the training set, which does not occur with regularized data. TheoptimalparametersinFig.7expresssomemean- ingful facts: first, the behavior of the initial conditions for exposed individuals reveals that defining this value just as a fixed proportion of the population size can undermine the capacity of the method for solving the system of differential equations; second, the parame- tersβ andω exhibit similar behaviors (especially when considering the results obtained with regularized data). Over time, the rate of contact between susceptible and infected individuals decreases (assuming the hypothesis that recovered individuals are immune for some time), so that quarantine measures end up being eased, which is reflected in the reduction of isolation measures. Def- initely, political and social interests also play a signifi- cant role in this behavior. The inspection of the behavior ofdP in Fig. 7 in- dicates that the mortality rate of positively diagnosed individuals basically only decreases after a certain run, . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 12 Gustavo B. Libotte et al. 1.5 3.0 4.5β ×10−8 Original Regularized 0.0 0.6 1.2 1.8ω ×10−2 0 25 50 75 100 125 Run 0.6 0.9 1.2dP ×10−1 0 25 50 75 100 125 Run 0.0 0.4 0.8 1.2E (0) ×103 Fig. 7 Optimal parameters obtained by fitting the daily data of infected and dead individuals, forθ = (β, ω, dP, E(0)). Each run is associated with a training set with a specific size, from 60 to 196 data in each set. The test data refer to the 14 subsequent data from the corresponding run. Parameters obtained by fitting original data are shown in blue, and with regularized data are shown in orange. around the gray shaded area. This behavior suggests that the parameter can be approximated by a function. In this case, we propose to representdP as a function of the form dP (t) = d0 exp (−d1t). (8) This approach assumes an inherent error related to the first runs, at expense of better representing the param- eter behavior in later times. Thus,dP becomes variable over time, and the vector of parameters to be estimated is formed byθ = (β, ω, d0, d1, E(0)). Considerthemethodologypreviouslyadopted,where the set of training data is gradually expanded with each parameter estimation, but now taking into ac- count the proposed approximation fordP (t) defined in Eq. (8). The search intervals for the new parameters are 0≤ d0≤ 1 and 0≤ d1≤ 1 and all other parameters follow the same specifications defined before. Figure 8 shows that this strategy is reflected in the reduction of the NRMSE considering the test data set with 14 values, in most simulated forecasts. The NRMSEs con- sidering regularized data and time-varyingdP (orange dots) represent a better approximation of the test data set in 73.53% of the analyzed runs, in relation to the re- sults considering original data (blue dots). The area un- der the curve obtained with regularized data anddP (t) is 733.06 u. a., whereas for original data the area is equal to 1427.28 u. a., which represents a reduction of approximately 48.64%. As in Fig. 7, we are interested in understanding the behavior of the parameters inherent to the function when dP varies over time. Figure 9 shows the parame- ters d0 and d1, as well as the other calibrated param- eters of the model, calculated in each parameter esti- mation procedure using both original and regularized data. The color scheme follows what has been adopted, blue for original data and orange for regularized data. In addition, Fig. 9 shows the range of Eq. (8), taking the optimal values ofd0 andd1, alongside the values of dP (shown in Fig. 7) for both original and regularized data. In this regard, we are interested in analyzing the behavior of dP (t) in terms of the values obtained for the case wheredP is constant. In the first runs,dP (t) is expressed by a nearly con- stant function, since the values ofd1 are very close to zero (see the behavior of d1 in Fig. 9). In this case, dP (t)≈d0. These curves practically coincide with the 0 25 50 75 100 125 Run 0.0 0.8 1.6 2.4 Normalized root-mean-square error ×102 Original Regularized Fig. 8 NRMSEs computed by comparing simulations of the SEIRPD-Q model, performed using parameters that best fit the training data, where θ = (β, ω, d0, d1, E(0)), in rela- tion to the test set (composed of 14 points). The area under the curve represents the total deviation relative to the test data in all runs, where the number of training data varies. . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 13 0 50 100 Run 0.0 0.8 1.6 2.4ω ×10−2 0 50 100 Run 1.5 3.0 4.5 6.0β ×10−8 0 50 100 Run 0.0 0.5 1.0 1.5E (0) ×103 0 25 50 75 100 125 Run 1.5 3.0 4.5 6.0d0, original data ×10−1 0 25 50 75 100 125 Run 0.0 0.8 1.6 2.4d1, original data ×10−2 0 25 50 75 100 125 Run, t 0 2 4 6 dP(t) original data ×10−1 0 25 50 75 100 125 Run, t 0.60 0.75 0.90 1.05 dP(t) regularized data ×10−1 0.75 0.90 1.05 d0, regularized data ×10−1 0.0 1.5 3.0 4.5 d1, regularized data ×10−3 Fig. 9 Optimal parameters obtained in a procedure similar to that of Fig. 7, forθ = (β, ω, d0, d1, E(0)). We also show the interval that includes the curves ofdP (t), for each optimal value ofd0 and d1 associated with the same run. In addition, we also show the optimal values ofdP presented in Fig. 7, for purposes of comparison with the curves ofdP (t). first points shown in the frame corresponding todP (t) obtained with regularized data in Fig. 9, which means that, in fact, the calibrations are quite similar. This be- havior occurs because, given the training set for such runs, the model identifies that the mortality rate is not decreasing and, therefore, the calibration procedure es- timates the most suitable function for such data (by means ofd0 and d1), which in this case is nearly con- stant, without loss of generality. As the training set gets larger,dP (t) starts to be- have as expected, exponentially decreasing. In this case, the first calibrations provide functions relatively distant from the corresponding points in Fig. 9, especially in the early times. In the last runs, the functions show good agreement with the compared points. However, it is important to note thatdP (t) is not expected to rep- resent the exact behavior of such data in the long run. In general, when the mortality rate varies over time, the compartmental model may be more capable of cap- turing the dynamics of the data, allowing for more ac- curate predictions. This hypothesis is supported by the last results obtained in Fig. 8, where the orange dots always represent the best approximation in relation to the test data. 4 Conclusions Our study provides a framework that aims to increase the predictive capacity of compartmental models. This is relevant from the point of view that the proposed strategies can be extended to other data sets and com- partmental models. Especially in the context of the epi- demiological modeling of COVID-19, approaches of this type can be useful, taking into account the wide range of existing compartmental models and the fact that the dataanalyzedherearesimilartoothersregardingnoise. This study has gone some way towards enhancing our understanding of the influence of noise on the es- timation of parameters of compartmental models. The work has revealed that the regularization of data by means of GPR can represent an alternative to miti- gate the effect of noise in a given parameter calibration. Since this procedure must not be repeatedly applied, in the context to which it is proposed, the computational cost of the GPR can be considered irrelevant. . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 14 Gustavo B. Libotte et al. Our research also suggests that it may be useful to use time-varying parameters, to the detriment of the usual approach that adopts constant parameters. It incorporates additional degrees of freedom into the model, in order to provide more flexibility to describe the behavior of the analyzed data. The choice of such a function depends on several factors, as for example the physical meaning of the parameter and the additional parameters inherent to the function. This analysis can be conducted for any model and, clearly, the gain is related to the appropriate choice of the function, espe- cially if many parameters to be estimated are included. Declarations Funding TheauthorswouldliketothanktheMinistryofScience, Technology, Innovation, and Communication (MCTIC) of Brazil. Gustavo Libotte and Lucas dos Anjos are supported by a postdoctoral fellowship from the In- stitutional Training Program (PCI) of the Brazilian National Council for Scientific and Technological De- velopment (CNPq), grant numbers 303185/2020-1 and 301327/2020-3, respectively. Conflict of interest The authors declare that the research was conducted in the absence of any commercial or financial relation- ships that could be constructed as a potential conflict of interest. The authors have no affiliation with any or- ganization with a direct or indirect financial interest in the subject matter discussed in the manuscript. Availability of data and material The data used in this work are publicly available ac- cording to Ref. [10]. Code availability The source code used to generate the results is pub- liclyavailableatgithub.com/gustavolibotte/enhancing- forecast-COVID-19.

1. Alberti, T., Faranda, D.: On the uncertainty of real-time predictions of epidemic growths: a COVID-19 case study for China and Italy. Communications in Nonlinear Sci- ence and Numerical Simulation90, 105372 (2020). DOI 10.1016/j.cnsns.2020.105372 2. Arias Velásquez, R.M., Mejía Lara, J.V.: Gaussian ap- proach for probability and correlation between the num- ber of COVID-19 cases and the air pollution in Lima. Urban Climate 33, 100664 (2020). DOI 10.1016/j.uclim. 2020.100664 3. Bailer-Jones, C.A.L.: Practical Bayesian Inference. Cam- bridge University Press, Cambridge (2017). DOI 10. 1017/9781108123891 4. Bhopal, S.S., Bhopal, R.: Sex differential in COVID-19 mortality varies markedly by age. The Lancet 396(10250), 532–533 (2020). DOI 10.1016/S0140-6736(20)31748-7 5. Bishop, C.M.: Pattern Recognition and Machine Learn- ing. Springer-Verlag, Berlin, Heidelberg (2006) 6. Brazilian Institute of Geography and Statistics: Demo- graphic Census (2020). URLhttps://www.ibge.gov.br/ cidades-e-estados/rj/rio-de-janeiro.html. Accessed November 20, 2020 7. Calvetti, D., Hoover, A.P., Rose, J., Somersalo, E.: Metapopulation network models for understanding, pre- dicting, and managing the coronavirus disease COVID- 19. Frontiers in Physics 8 (2020). DOI 10.3389/fphy. 2020.00261 8. Ching, J., Chen, Y.C.: Transitional Markov Chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging. Journal of Engineer- ing Mechanics 133(7), 816–832 (2007). DOI 10.1061/ (ASCE)0733-9399(2007)133:7(816) 9. Clark, A., Jit, M., Warren-Gash, C., Guthrie, B., Wang, H.H.X., Mercer, S.W., Sanderson, C., McKee, M., Troeger, C., Ong, K.L., Checchi, F., Perel, P., Joseph, S., Gibbs, H.P., Banerjee, A., Eggo, R.M., Nightingale, E.S., O’Reilly, K., Jombart, T., Edmunds, W.J., Rosello, A., Sun, F.Y., Atkins, K.E., Bosse, N.I., Clifford, S., Rus- sell, T.W., Deol, A.K., Liu, Y., Procter, S.R., Leclerc, Q.J., Medley, G., Knight, G., Munday, J.D., Kucharski, A.J., Pearson, C.A.B., Klepac, P., Prem, K., Houben, R.M.G.J., Endo, A., Flasche, S., Davies, N.G., Diamond, C., van Zandvoort, K., Funk, S., Auzenbergs, M., Rees, E.M., Tully, D.C., Emery, J.C., Quilty, B.J., Abbott, S., Villabona-Arenas, C.J., Hué, S., Hellewell, J., Gimma, A., Jarvis, C.I.: Global, regional, and national estimates of the population at increased risk of severe COVID-19 due to underlying health conditions in 2020: a modelling study. The Lancet Global Health 8(8), e1003–e1017 (2020). DOI 10.1016/S2214-109X(20)30264-3 10. Cota, W.: Monitoring the number of COVID-19 cases and deaths in Brazil at municipal and federative units level. SciELOPreprints:362 (2020). DOI 10.1590/ scielopreprints.362. URL https://doi.org/10.1590/ scielopreprints.362 11. Cupertino, G.A., Cupertino, M.d.C., Gomes, A.P., Braga, L.M., Siqueira-Batista, R.: COVID-19 and brazil- ian indigenous populations. The American Journal of Tropical Medicine and Hygiene103(2), 609–612 (2020). DOI 10.4269/ajtmh.20-0563 12. Currie, C.S., Fowler, J.W., Kotiadis, K., Monks, T., Onggo, B.S., Robertson, D.A., Tako, A.A.: How simu- lation modelling can help reduce the impact of COVID- 19. Journal of Simulation 14(2), 83–97 (2020). DOI 10.1080/17477778.2020.1751570 13. Davies,N.G.,Klepac,P.,Liu,Y.,Prem,K.,Jit,M.,Eggo, R.M.: Age-dependent effects in the transmission and con- . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint Framework for enhancing the estimation of model parameters for data with a high level of uncertainty 15 trol of COVID-19 epidemics. Nature Medicine 26(8), 1205–1211 (2020). DOI 10.1038/s41591-020-0962-9 14. Edridge, A.W.D., Kaczorowska, J., Hoste, A.C.R., Bakker, M., Klein, M., Loens, K., Jebbink, M.F., Matser, A., Kinsella, C.M., Rueda, P., Ieven, M., Goossens, H., Prins, M., Sastre, P., Deijs, M., van der Hoek, L.: Seasonal coronavirus protective immunity is short- lasting. Nature Medicine (2020). DOI 10.1038/ s41591-020-1083-1 15. Gujarati, D.N., Porter, D.C.: Basic Econometrics, 5 edn. McGraw-Hill Irwin, New York (2008) 16. Holmdahl, I., Buckee, C.: Wrong but useful—what COVID-19 epidemiologic models can and cannot tell us. New England Journal of Medicine 383(4), 303–305 (2020). DOI 10.1056/NEJMp2016822 17. Ioannidis, J.P.A.: Coronavirus disease 2019: The harms of exaggerated information and non-evidence-based mea- sures. European Journal of Clinical Investigation50(4) (2020). DOI 10.1111/eci.13222 18. Jia, J., Ding, J., Liu, S., Liao, G., Lin, J., Duan, B., Wang, G., Zhang, R.: Modeling the control of COVID-19: impact of policy interventions and meteorological factors. Electronic Journal of Differential Equations2020(23), 1– 24 (2020) 19. Ketu, S., Mishra, P.K.: Enhanced Gaussian process regression-based forecasting model for COVID-19 out- break and significance of IoT for its detection. Applied Intelligence (2020). DOI 10.1007/s10489-020-01889-9 20. Kocijan, J.: Modelling and Control of Dynamic Sys- tems Using Gaussian Process Models. Springer In- ternational Publishing, Cham (2016). DOI 10.1007/ 978-3-319-21021-6_2 21. Li, R., Pei, S., Chen, B., Song, Y., Zhang, T., Yang, W., Shaman, J.: Substantial undocumented infection fa- cilitates the rapid dissemination of novel coronavirus (SARS-CoV-2). Science368(6490),489–493(2020). DOI 10.1126/science.abb3221 22. Link, W.A., Barker, R.J.: Bayesian Inference. Academic Press,London(2010). DOI10.1016/B978-0-12-374854-6. 00004-1 23. Marson, F.A.L.: COVID-19 – 6 million cases worldwide andanoverviewofthediagnosisinBrazil:atragedytobe announced. Diagnostic Microbiology and Infectious Dis- ease 98(2), 115113 (2020). DOI 10.1016/j.diagmicrobio. 2020.115113 24. Massonis, G., Banga, J.R., Villaverde, A.F.: Structural identifiability and observability of compartmental models of the COVID-19 pandemic (2020) 25. McAloon, C., Collins, Á., Hunt, K., Barber, A., Byrne, A.W., Butler, F., Casey, M., Griffin, J., Lane, E., McEvoy, D., Wall, P., Green, M., O’Grady, L., More, S.J.: Incubation period of COVID-19: a rapid sys- tematic review and meta-analysis of observational re- search. BMJ Open10(8), e039652 (2020). DOI 10.1136/ bmjopen-2020-039652 26. Mohammed, R.O., Cawley, G.C.: Over-fitting in model selection with Gaussian Process Regression. In: P. Perner (ed.) Machine Learning and Data Mining in Pattern Recognition, pp. 192–205. Springer International Pub- lishing, Cham (2017) 27. Murphy, K.P.: Machine Learning: A Probabilistic Per- spective. The MIT Press, Cambridge (2012) 28. Nelder, J.A., Mead, R.: A simplex method for func- tion minimization. The Computer Journal7(4), 308–313 (1965). DOI 10.1093/comjnl/7.4.308 29. Nocedal, J., Wright, S.: Numerical Optimization. Springer New York, New York (2006). DOI 10.1007/ 978-0-387-40065-5 30. Oran, D.P., Topol, E.J.: Prevalence of asymptomatic SARS-CoV-2 Infection. Annals of Internal Medicine 173(5), 362–367 (2020). DOI 10.7326/M20-3012 31. Overton, C.E., Stage, H.B., Ahmad, S., Curran- Sebastian, J., Dark, P., Das, R., Fearon, E., Felton, T., Fyles, M., Gent, N., Hall, I., House, T., Lewkow- icz, H., Pang, X., Pellis, L., Sawko, R., Ustianowski, A., Vekaria, B., Webb, L.: Using statistics and mathemati- cal modelling to understand infectious disease outbreaks: COVID-19 as an example. Infectious Disease Modelling 5, 409–441 (2020). DOI 10.1016/j.idm.2020.06.008 32. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cour- napeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Ma- chine Learning Research12, 2825–2830 (2011) 33. Polidoro, M., de Assis Mendonça, F., Meneghel, S.N., Alves-Brito,A.,Gonçalves,M., Bairros,F., Canavese,D.: Territories under siege: risks of the decimation of indige- nous and quilombolas peoples in the context of COVID- 19 in south Brazil. Journal of Racial and Ethnic Health Disparities (2020). DOI 10.1007/s40615-020-00868-7 34. Puntanen, S., Styan, G.P.H.: Schur complements in statistics and probability. In: F. Zhang (ed.) The Schur Complement and Its Applications, pp. 163–226. Springer US, Boston (2005). DOI 10.1007/0-387-24273-2_7 35. Rasmussen, C.E., Williams, C.K.I.: Gaussian Processes for Machine Learning. Adaptive Computation and Ma- chine Learning. MIT Press, Cambridge (2006) 36. Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U., Timmer, J.: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics 25(15), 1923–1929 (2009). DOI 10.1093/ bioinformatics/btp358 37. Ribeiro, F., Leist, A.: Who is going to pay the price of COVID-19? Reflections about an unequal Brazil. Inter- national Journal for Equity in Health19(1), 91 (2020). DOI 10.1186/s12939-020-01207-2 38. Ribeiro, M.H.D.M., Silva, R.G., Mariani, V.C., Coelho, L.d.S.: Short-term forecasting COVID-19 cumulative confirmed cases: perspectives for Brazil. Chaos, Solitons & Fractals 135, 109853 (2020). DOI 10.1016/j.chaos. 2020.109853 39. de Ridder, D., Tax, D.M.J., Lei, B., Xu, G., Feng, M., Zou, Y., van der Heijden, F.: Parameter Estimation, chap. 4, pp. 77–113. John Wiley & Sons, Ltd (2017). DOI 10.1002/9781119152484.ch4 40. Roda, W.C., Varughese, M.B., Han, D., Li, M.Y.: Why is it difficult to accurately predict the COVID-19 epidemic? Infectious Disease Modelling 5, 271–281 (2020). DOI 10.1016/j.idm.2020.03.001 41. Salvatier, J., Wiecki, T.V., Fonnesbeck, C.: Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2, e55 (2016). DOI 10.7717/peerj-cs.55 42. Schulz, E., Speekenbrink, M., Krause, A.: A tutorial on Gaussian Process regression: Modelling, exploring, and exploiting functions. Journal of Mathematical Psychol- ogy 85, 1–16 (2018). DOI 10.1016/j.jmp.2018.03.001 43. Shi, J.Q., Choi, T.: Gaussian Process Regression Analysis for Functional Data, 1 edn. Chapman and Hall/CRC, New York (2011). DOI 10.1201/b11038 44. Storn, R., Price, K.: Differential Evolution—a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization 11(4), 341–359 (1997). DOI 10.1023/A:1008202821328 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint 16 Gustavo B. Libotte et al. 45. Sun, N.Z., Sun, A.: The Classical Inverse Problem, chap. 2, pp. 25–67. Springer New York, New York (2015). DOI 10.1007/978-1-4939-2323-6_2 46. Taghizadeh, L., Karimi, A., Heitzinger, C.: Uncertainty quantification in epidemiological models for the COVID- 19 pandemic. Computers in Biology and Medicine125, 104011 (2020). DOI 10.1016/j.compbiomed.2020.104011 47. Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial and Applied Mathematics, Philadelphia (2005). DOI 10.1137/1.9780898717921 48. Torres, T.S., Hoagland, B., Bezerra, D.R.B., Garner, A., Jalil, E.M., Coelho, L.E., Benedetti, M., Pimenta, C., Grinsztejn, B., Veloso, V.G.: Impact of COVID-19 pan- demic on sexual minority populations in Brazil: an anal- ysis of social/racial disparities in maintaining social dis- tancing and a description of sexual behavior. AIDS and Behavior (2020). DOI 10.1007/s10461-020-02984-1 49. Veiga e Silva, L., de Andrade Abi Harb, M.D.P., Teix- eira Barbosa dos Santos, A.M., de Mattos Teixeira, C.A., Macedo Gomes, V.H., Silva Cardoso, E.H., S da Silva, M., Vijaykumar, N.L., Venâncio Carvalho, S., Ponce de Leon Ferreira de Carvalho, A., Lisboa Frances, C.R.: COVID-19 mortality underreporting in Brazil: analysis of data from government internet portals. Journal of Medical Internet Research 22(8), e21413 (2020). DOI 10.2196/21413 50. Volpatto, D.T., Resende, A.C.M., Anjos, L., Silva, J.V.O., Dias, C.M., Almeida, R.C., Malta, S.M.C.: A generalized SEIRD model with implicit social distanc- ing mechanism: a Bayesian approach for the identifica- tion of the spread of COVID-19 with applications in Brazil and Rio de Janeiro state. medRxiv (2020). DOI 10.1101/2020.05.30.20117283 51. Wu, S.L., Mertens, A.N., Crider, Y.S., Nguyen, A., Pokpongkiat, N.N., Djajadi, S., Seth, A., Hsiang, M.S., Colford, J.M., Reingold, A., Arnold, B.F., Hubbard, A., Benjamin-Chung, J.: Substantial underestimation of SARS-CoV-2 infection in the United States. Nature Communications 11(1), 4507 (2020). DOI 10.1038/ s41467-020-18272-4 52. Zhou, T., Ji, Y.: Semiparametric Bayesian inference for the transmission dynamics of COVID-19 with a state- space model. Contemporary Clinical Trials 97, 106146 (2020). DOI 10.1016/j.cct.2020.106146 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 8, 2021. ; https://doi.org/10.1101/2020.12.17.20248389doi: medRxiv preprint