A stochastic Bayesian bootstrapping model for COVID-19 data.

We provide a stochastic modeling framework for the incidence of COVID-19 in Castilla-Leon (Spain) for the period March 1, 2020 to February 12, 2021, which encompasses four waves. Each wave is appropriately described by a generalized logistic growth curve. Accordingly, the four waves are modeled through a sum of four generalized logistic growth curves. Pointwise values of the twenty input parameters are fitted by a least-squares optimization procedure. Taking into account the significant variability in the daily reported cases, the input parameters and the errors are regarded as random variables on an abstract probability space. Their probability distributions are inferred from a Bayesian bootstrap procedure. This framework is shown to offer a more accurate estimation of the COVID-19 reported cases than the deterministic formulation.

Introduction

COVID-19 is an infectious disease caused by coronavirus SARS-CoV-2. It was detected for the first time in Wuhan, China, in December 2019, and quickly spread around the globe becoming an ongoing pandemic. The virus is rapidly transmitted between persons through small droplets and aerosols. The most common symptoms of the disease are fever, dry cough and fatigue. Data to date suggest that 80% of infections are mild or asymptomatic, 15% are severe, and 5% are critical. Lethality strongly depends on age and comorbidities. As of July 2021, more than 190 million people have been infected and more than 4 million people have died. To contain the spread of the virus and alleviate the pressure on the health systems, governments put several restrictions such as city or country lockdowns, quarantines, social distancing and hygiene measures, curfews, mandatory masks, etc. (World Health Organization 2021; Wu et al. 2020; Kraemer et al. 2020; Wu et al. 2020).

The use of mathematical models is an effective tool to describe and predict the evolution of epidemics and to propose targeted measures (Chitnis et al. 2010; Remuzzi and Remuzzi 2020). Due to the fast transmissibility of SARS-CoV-2 and the containment measures frequently implemented by governments, the modeling of its spread is a difficult problem. For example, as already noticed by other researchers (Acedo et al. 2016; Muñoz-Fernández et al. 2021; Moein et al. 2021), the usual autonomous SIR (susceptible-infected-recovered) model cannot capture the quick variations of COVID–19 reported infections.

When aggregated time-series data are present, the logistic differential equation model may be useful to fit the measurements on COVID-19 infections (Wang et al. 2020). The curve is characterized by an increasing growth in the beginning period, combined with a decreasing growth at a later stage. Generalizations of the logistic growth curve, to allow for more general sigmoid shapes, have been employed to model COVID-19 (Pelinovsky et al. 2020; Wu et al. 2020; Aviv-Sharon and Aharoni 2020; Lee et al. 2020). Those approaches seem to be adequate when fitting a single wave of the COVID-19 epidemic. For multiple phases of growth, as occurs with COVID-19 infection waves, a single S-shaped trajectory is not applicable. In this line, ideas from (Meyer 1994) on the bi-logistic model have been used to fit cumulative cases of COVID-19 along two different outbreaks or waves (Zhang et al. 2020; Salpasaranis and Stylianakis 2020; da Silva et al. 2020). These types of models, usually called phenomenological or statistical, are often useful to reproduce and forecast the course of an epidemic (Chowell et al. 2016; Pell et al. 2018), when the insight is limited, treatments and interventions are rapidly changing, data vary abruptly, and mechanistic models (compartmental models with laws of transmission) present difficulties.

The incidence of an epidemic and its modeling have intrinsic uncertainties that are irreducible (random uncertainty). Thus, to better mimic reality, models should incorporate stochastic components (Smith 2013). For example, reference (Lee et al. 2020) also included stochastic effects through a Bayesian formalization. Bayesian inference has been suggested for other models to accommodate COVID-19 data, such as the Gompertz curve (Berihuete et al. 2021). However, the use of the Bayesian bootstrap (Rubin 1981) to infer input uncertainties does not seem to have been investigated so far.

In this paper, we investigate the use of four combined generalized logistic differential equations to model the four waves of the COVID-19 epidemic in Castilla-Leon (a Spanish autonomous region with 2.5 million inhabitants) at once. The temporal period runs from March 1, 2020 to February 12, 2021. The calibration of model parameters is conducted by a least-squares minimization procedure. Due to data measurement errors, some form of stochasticity is introduced into the model. The input parameters and the errors are considered as random variables, whose probability distributions are inferred from a Bayesian bootstrap technique.

The plan of the paper is the following. Section 2 provides the methodology, including a brief description of the data, and the deterministic and stochastic approaches. Section 3 presents the results (numerical calculations, fittings and plots). A discussion comes in Sect. 4. The paper ends with some final conclusions in Sect. 5.

Methods

Data

We have data on the number of new daily COVID-19 reported infections in the Spanish autonomous region of Castilla-Leon. This region is the largest community in Spain by area, it is located in the northwest of Spain, and it has a population of around 2.5 million. The data correspond to the temporal period of almost one year, from March 1, 2020 to February 12, 2021, which encompasses four waves of the epidemic. The cases have been retrieved from the open data portal of Castilla-Leon1. This dataset only captures a small fraction of the true burden, due to asymptomatic cases, lack of resources and omission of suspected but not confirmed cases. In this paper, we treat Castilla-Leon as a whole where people interact homogeneously.

In Fig. 1, the number of daily new reported infections, for 349 consecutive days, is depicted in the left panel, while the accumulated number of daily reported infections is shown in the right panel. We note that there are four clear waves of the epidemic. The first wave corresponds to the first entrance of the virus in Spain, which ended up in summer 2020 due to the severe lockdown imposed by authorities. The second wave started after summer 2020 due to relaxation of measures and overlapped with a larger third wave along autumn. Finally, the fourth wave began after the end-of-year vacations and ended in February 2021 due to some restrictions and the vaccination program. A significant variability in the daily data, with abruptly increasing and decreasing magnitudes, is observed between nearby days, highlighting some sort of uncertainty entailing some stochastic nature in the data. This may be due to highly variable factors, such as the quantity of tests available and performed, symptomatology, etc. The implementations and computations are performed with Mathematica®, version 12.0, and are included as supplementary material, where the data are available (variable vtotal).

Fig. 1

Left panel: number of new daily reported infections. Right panel: accumulated daily number of reported infections

Multiple generalized logistic growth curves

In this subsection, the reported infections of the COVID-19 disease presented in Fig. 1 are modeled. Each wave of the epidemic—the accumulated version from the right panel of Fig. 1—is described by means of a generalized logistic growth curve. The four waves are then modeled by juxtaposing four generalized logistic growth curves. The input parameters of the complete model are calibrated by certain optimization procedures. This framework provides a smooth curve that fits the data from both panels of Fig. 1. Note that in this subsection, stochastic effects are not taken into account yet.

A generalized logistic differential equation

The Malthusian model, proposed by T.R. Malthus in 1798 in an essay (Malthus 1999), describes an exponential growth in a population through the ordinary differential equationwhere a is a positive parameter defined as the intrinsic growth rate. Given an initial condition , the solution of the above equation is given byA modern formulation of the Malthusian growth model can be read at any introductory text (Murray 2002). In the field of population ecology, it is considered as the first law of population dynamics (Turchin 2001).

In order to capture the decrease in the growth rate with time, P.F. Verhulst proposed in 1838 (Verhulst 1838; Kingsland 1982) the logistic model, given bywhere is the carrying capacity (the limit of y(t) as t tends to infinity). This is a Bernoulli-type ordinary differential equation. Given an initial condition , the solution is known in closed formThis function was employed to forecast an Ebola epidemic (Chowell et al. 2014). The saturation effect implicitly captures public health interventions, without complex mechanistic assumptions about the transmission process.

To allow for more flexible S-shaped curves to model growth phenomena over time, the following modification of the logistic differential equation has been suggested in the literatureOriginally, (1) was used for the analysis of tumor growth (Marusic et al. 1994; Spratt et al. 1996; Birch 1999; Sachs et al. 2001), though applications in epidemiology are also found. Examples of diseases include SARS (Hsieh et al. 2004; Hsieh 2009), dengue fever (Hsieh 2005), influenza H1N1 (Hsieh 2010), Zika (Chowell et al. 2016), Ebola (Pell et al. 2018), and COVID-19 (Pelinovsky et al. 2020; Wu et al. 2020; Aviv-Sharon and Aharoni 2020; Lee et al. 2020). Here is a power that controls the asymmetry of the curve and how fast the limiting number K is approached. It endows the model with higher flexibility. When , the classical logistic differential equation is obtained, and when b tends to 0, the Gompertz equation is recovered. This generalization belongs to the class of Bernoulli differential equations too. If is the initial condition, the solution takes the formThis is called a generalized logistic growth curve (or sometimes Richards’ curve). It may be appropriate to model the aggregated cases of a single wave of the COVID-19 epidemic. For new cases (not accumulated), consecutive differences are considered.

Combination of growth curves

For multiple phases of growth, a single sigmoid curve is not appropriate to describe such data. Thus, we propose a combination of generalized logistic growth curves of the form (2). This is an extension of the work initiated by P.S. Meyer for the bi-logistic model (Meyer 1994), with subsequent applications in sociology (Fenner et al. 2013), agriculture (Shehu 2015) or epidemiology (Zhang et al. 2020; Salpasaranis and Stylianakis 2020; da Silva et al. 2020; Lavrova et al. 2017), for instance.

Mathematically, a combination of generalized logistic growth curves takes the following formThis sum of trajectories, supplied with four terms , allows modeling the accumulated cases of the four concatenated COVID-19 waves. In practice, one should try fitting several concatenated models and compare their goodness-of-fit. In our case, we tried three terms, but the fit was not good. For new cases (not accumulated), consecutive differences are considered. It is important to note that the four generalized logistic growth curves are not independent (the four waves are not treated independently). Finally, notice that captures the beginning of the i-th wave. Each measures the highest infection level of the i-th wave.

A deterministic fit

Let be the number of new reported cases at time . Let be the number of accumulated reported cases at time , scaled by the total population in Castilla-Leon ( inhabitants). The following simple relations holdA combination of generalized logistic growth curves as in (3) is used to model . The parameters of (3) are estimated at once by a deterministic least-squares procedure (Stanescu et al. 2009),In most of the cases, the model curve will not go through all the data, so this minimum value will not be zero. The minimum gives a measure of how good the fit is.

Parametric randomization of the model

We now proceed with modeling the daily random variability of the data through a probabilistic setting. Following a Bayesian formalism, the parameters and the errors will be regarded as random variables on an abstract probabilistic space. Their probability distributions will be then inferred by means of Bayesian bootstrapping. In this context, the output of the model is a stochastic process, which will render a stochastic fit of the data, rather than an averaged estimation.

Randomization

Part of the variability of the data is not captured by the deterministic model. We assume the errors are random variables naturally defined in a probabilistic framework. We also consider the lack of knowledge on the input parameters, which prescribe the constitutive laws of the system, represented within a probabilistic framework (Smith 2013, chapter 1), (Le Maître and Knio 2010, chapter 1). Thus, we consider in this new setup the parameters and the errors of the model as random variables.

The field of uncertainty quantification studies the impact of random uncertainties on models (Smith 2013). This quantification is necessary to evaluate the discrepancies between the model predictions and the current system behavior. Inverse uncertainty quantification deals with inference of the probability distributions of the parameters from the data. These probability distributions are not, in general, independent. Forward uncertainty quantification extracts the main statistical content of the model output, once the probability distributions of the parameters are fixed (Le Maître and Knio 2010; Xiu 2010). The various stages of the probabilistic modeling process are schematically illustrated in Fig. 2. Inverse uncertainty quantification is not an easy task. Here we rely on the Bayesian bootstrap technique. Forward uncertainty quantification will be conducted via Monte Carlo simulation.

Fig. 2

Schematic illustration of the various stages of the probabilistic modeling process

Bayesian bootstrap

Given a mathematical model, the model error coming from a deterministic least-squares optimization technique may be regarded as a random variable X. The error varies with time, and thus the errors at different instants of time may be seen as copies of X. In this context, the errors of the model are identically distributed and independent random variables , where m is the length of the data. Actually, these random variables are unknown; only different realizations are available from the data and the deterministic fit. The bootstrap methodology assumes that the observed residuals are all possible distinct values of X, based on the principle that all observed variables are discrete.

The Bayesian bootstrap, developed by D.B. Rubin (Rubin 1981), infers the distribution of X by resampling with Dirichlet weights of coefficients . This procedure corresponds to the following hierarchical Bayesian statistical modelwithBayes’ theorem yields the posterior distributionbecause, for ,That is, the posterior proportions are uniformly distributed on the simplex. Here is the Categorical distribution on and stands for the Dirichlet distribution, which is a conjugate prior. Parameters and come from the frequencies in before and after observing them, respectively. Note that is an improper prior.

The distribution can be sampled as follows. From independent realizations , assume that these values are ordered as . Set and . Define , . Then are m independent realizations of .

For each resampling of with Dirichlet weights of coefficients , the input parameters of the model are determined by a least-squares fitting. Formally, it is assumed thatThis gives rise to samples of the input parameters, so their (posterior) probability distributions may be inferred asThis allows solving the problem of inverse uncertainty quantification.

Monte Carlo simulation

Monte Carlo simulation is a popular method for forward uncertainty quantification. It is simple to implement and robust. It uses a collection of independent random realizations of the model response, usually obtained from deterministic numerical techniques. The statistics of the model response are derived from the statistics of that sample. For example, the mean of a function h of the model output is estimated by , by the law of large numbers. Monte Carlo simulation essentially amounts to conducting M deterministic resolutions, where M is generally large. The robustness of the method is due to its independence of the random dimensionality, the variable t, or regularity issues. Thus, in contrast to spectral methods, its use is advantageous when there is a large number of input random parameters or the variable t may be large. Further, if the inverse parameter estimation method generates realizations of the parameters (such as Bayesian methods), then Monte Carlo simulation seems the logic option for forward uncertainty quantification. The convergence of the Monte Carlo estimate behaves as due to the central limit theorem. It is assessed based on the estimated statistics. The reader is referred to (Smith 2013; Le Maître and Knio 2010; Xiu 2010) for similar discussions on Monte Carlo sampling.

A stochastic fit

The aim here is to randomize the deterministic model. The input parameters and the model errors are assumed to be random variables. Note that to apply the Bayesian bootstrap methodology, one needs errors that are identically distributed and independent. However, the fit for the accumulated infections does not yield independent residuals. Indeed, two consecutive residuals are correlated, because of the increasing character of the curve. Also, in addition, the fit for the daily new infections does not yield identically distributed errors. Certainly, a resampling of residuals may give rise to negative data points, which does not make sense. To fix these issues and achieve, as far as possible, an independent and identically distributed sample, the (natural) logarithms of the daily new infections will be considered. While the parameters of model (3) for cumulative infections are calibrated by least-squares fitting, the resampling is performed for the residuals obtained for the logarithms of the daily new infections. For each resampling of these residuals, we come back to cumulative infections and refit model (3). In this way, realizations of the model parameters are obtained for each refit.

The following steps summarize the procedure for estimating the probability distributions of the parameters: Once and are available, the following steps outline the procedure to calculate the realizations of the model output:

Determine the residuals for the logarithms of the daily new infections: . The parameters of y were previously determined by the deterministic least-squares fitting (4).

Start a FOR loop to generate M bootstrap samples.

Resample the residuals with Dirichlet weights of coefficients . Keep the resampling as a vector with components .

Define . These are the new generated data for the logarithms of the daily new infections.

Let . These are the new generated data for the accumulated infections.

Fit the deterministic model (3) to by least-squares fitting. Keep the parameters as a vector .

End the FOR loop.

An ensemble of realizations for the residuals (errors) of the logarithms of the daily new infections. An ensemble of realizations for the input parameters.

Start a FOR loop over .

With parameters , evaluate . This gives a vector , a discrete sample path of the model for the logarithms of the daily new infections.

Incorporate the error: . This is a a discrete sample path of the model for the logarithms of the daily new infections, taking into account the random errors.

Let . Here is a discrete sample path for the accumulated infections.

End the FOR loop.

The M discrete sample paths of the model output for the accumulated infections. From them, statistics such as the mean, the variance, quantiles, etc. may be determined by means of Monte Carlo simulation.

Results

Deterministic fit

We have used the built-in function FindFit in Mathematica® to obtain the optimal parameters in (4), which are (up to four significant digits) as follows:The fitted model is depicted in Fig. 3. In the left panel, the aggregated reported infections in percentages, , are shown, while the right panel shows the estimated new daily infections, . The combination of generalized logistic growth curves allows for a good fit of the COVID-19 incidence, at least from an averaged (smoothed) point of view. In the following subsection, stochasticity will be taken into account to deal with the daily random variability of the measurements. In Fig. 4, the predictability of the model is quantified. We use two train sets, up to and up to , and predict the epidemic size a few days later. The forecast is reasonably good for two weeks. It seems to be better at than at , possibly due to the fact that corresponds to the peak of the second wave, while corresponds to the very early growth phase of the fourth wave. Predictions seem to be better when a larger dataset of the wave is available.

Fig. 3

Left panel: deterministic fit for the daily accumulated reported infections in percentages. Right panel: deterministic fit for the new daily reported infections

Fig. 4

Left panel: deterministic prediction for the daily accumulated reported infections in percentages. Right panel: deterministic prediction for the new daily reported infections. The vertical dashed line indicates the end of the calibration period

Stochastic fit

Table 1 reports the estimated marginal posterior statistics (mean, median, standard deviation, and quantiles 0.025 and 0.975) of the corresponding parameters. Figure 5 plots histograms for some marginal posterior distributions. Realizations are employed, for bootstrap samples (larger bootstrap samples do not render significant differences at the scale of the figures). In Fig. 6, the fit of the randomized model is shown, both for the accumulated reported infections (left panel) and the new reported infections per day (right panel). We draw mean values and regions of probability 0.95 for . Optical inspection shows that the estimated mean values for the model response are very similar to the deterministic fit from Fig. 3. Thus, the deterministic fit is extended by incorporating probabilistic features. In the stochastic approach, the probability regions must contain the variability of the data, but in a correct way, in the sense that the realizations generated should resemble the pattern of the data. A stochastic method whose probability regions are unnecessarily wide is not good, despite containing all recorded measurements. To better appreciate the similarity between the real data and the stochastic model, compared to the deterministic counterpart, some realizations (, and ) are plotted in Fig. 7. In Fig. 8, the predictability of the model is illustrated. We use two sets of incidence data, up to and up to as in the deterministic subsection.

Table 1

Estimated marginal posterior statistics (mean, median, standard deviation, and quantiles 0.025 and 0.975) of the input random parameters

	Mean	Median	SD	Quantile 0.025	Quantile 0.975
b_1	0.350	0.283	0.224	0.103	0.892
a_1	0.580	0.513	0.191	0.361	0.962
K_1	0.018	0.018	0.004	0.010	0.024
y_0,2	0.011	0.010	0.003	0.005	0.018
b_2	0.564	0.674	0.316	0.074	0.930
a_2	0.248	0.137	0.203	0.104	0.769
K_2	0.011	0.010	0.005	0.006	0.019
t_0,2	115.1	114.8	2.839	111.7	118.4
y_0,3	0.019	0.019	0.002	0.016	0.021
b_3	4.183	4.490	0.857	2.146	4.985
a_3	0.025	0.022	0.043	0.020	0.027
K_3	0.055	0.055	0.004	0.049	0.060
t_0,3	212.6	212.4	1.412	211.5	213.7
y_0,4	0.001	0.001	0.0003	0.0004	0.002
b_4	0.712	0.814	0.269	0.131	0.984
a_4	0.220	0.168	0.138	0.122	0.690
K_4	0.032	0.032	0.003	0.026	0.037
t_0,4	304.4	304.2	1.138	302.5	306.8

Fig. 5

Histograms for some marginal posterior distributions

Fig. 6

Left panel: stochastic fit for the accumulated daily reported infections in percentages. Right panel: stochastic fit for the new daily reported infections

Fig. 7

Left panel: some realizations of the stochastic model for the accumulated daily reported infections in percentages. Right panel: some realizations of the stochastic model for the new daily reported infections

Fig. 8

Left panel: stochastic prediction for the daily accumulated reported infections in percentages. Right panel: stochastic prediction for the new daily reported infections. The vertical dashed line indicates the end of the calibration period

Discussion

The COVID-19 spread in Castilla-Leon is modeled along the four waves in a single equation, both from the deterministic and the stochastic points of view. Generalized logistic differential equations are concatenated four times to deal with the four waves at once. This procedure extends the use of the generalized logistic differential equation to more than one wave, compared to (Pelinovsky et al. 2020; Wu et al. 2020; Aviv-Sharon and Aharoni 2020; Lee et al. 2020), and to more than two waves, compared to (Zhang et al. 2020; Salpasaranis and Stylianakis 2020; da Silva et al. 2020) from the work (Meyer 1994) on the bi-logistic model. The deterministic model gives a good fit for COVID-19 reported cases, especially for accumulated cases. Daily new reported infections present significant variability, perhaps due to highly variable factors, such as amount of tests available and performed, symptomatology, etc. Thus, stochasticity is incorporated, to extend the deterministic model and to obtain realizations that resemble the irregular dynamics of the reported new cases closer (Smith 2013). The Bayesian bootstrap (Rubin 1981), together with a trick to deal with residuals, is employed to set probability distributions for the input parameters and the errors. The potential of the Bayesian bootstrap for mathematical models with uncertainties does not seem to have been investigated. Our approach resembles the previous use of the frequentist bootstrap (Efron 1979) for mathematical models with deterministic parameters and random errors (Dogan 2007). Apart from fitting the available data, predictions have been performed by using two train sets, for a quite large forecast period of 15 days. Phenomenological models like the one proposed in this paper may be useful to generate reasonable forecasts in near time of the incidence of starting and advanced epidemic outbreaks, without entering in accounting for possible underlying mechanisms of the studied phenomenon (such as temporal or spatial dependencies), which would clearly be the base for further analysis. This idea was also discussed in (Chowell et al. 2016; Pell et al. 2018) for a single generalized logistic model.

We justify the use of the Bayesian bootstrap by commenting the unfeasibility of other techniques:Some limitations of the present work, which define potential avenues for future research, are the following: (a) Ideally, the probabilistic interval for the model output should be narrower. It would be of interest to investigate alternative deterministic models (the averaged, smooth curve) or Bayesian approaches. (b) We needed to consider the logarithms of the daily new infections to have an adequate resampling of residuals for which increasing sample paths (accumulated infections) or negative values (daily new infections) were not a problem. It would be interesting to directly define a model error for the accumulated infections that takes into account the correlation between successive days. This correlation is due to the increasing character of the curve. (c) Spatial, behavioral, or environmental effects have been neglected. These conditions may provide a more faithful fit of the COVID-19 data, albeit at the expense of higher complexity. Nonetheless, when uncertainty is present on the phenomenon itself and the data, sometimes it may be preferable to simply consider the model as uncertain, rather than augmenting its complexity.

Maximum entropy principle. This principle has been used in the literature to infer consistent probability distributions for input parameters (Dorini and Sampaio 2012). The probability density function of the parameter is taken by maximizing the Shannon entropy functional, often restricted to a certain support, to a mean value equal to the deterministic estimate, and rarely to a variance if available. In the case studied in the present paper, reliable supports of the parameters are unknown a priori.

General Bayesian model. Not restricted to the Bayesian bootstrap, Markov Chain Monte Carlo algorithms may solve the problem for any set of prior distributions (Smith 2013; Lesaffre and Lawson 2012). However, the large amount of input random parameters in our case study prevented us from using this option.

Itô-type stochastic differential equations. One could naturally ask about the incorporation of a white noise (formal derivative of Brownian motion) into the generalized logistic differential equation (1) (Mao 2007; Allen 2007). However, accumulated infections give rise to increasing sample paths, which would contradict the everywhere non-differentiability of Itô processes.

Conclusion

We have shown in the paper that the sum of four generalized logistic growth curves allows for a proper fit of the accumulated reported infections along the four waves of the COVID-19 epidemic in Castilla-Leon (Spain). Daily new reported infections are described by consecutive differences. The input parameters are pointwise calibrated by least-squares fitting. However, this calibration lacks of probabilistic interpretations.

Taking into account the significant variability in the daily reported data, with noisy features, stochasticity is incorporated into the model by treating the input parameters and the model errors as random variables. This conception of uncertainty matches with the Bayesian formalism of Statistics. The Bayesian bootstrap is an adequate approach for inverse uncertainty quantification and infers the probability distributions of the parameters. The model response is stochastic and includes realizations that permit a more reliable fit of the daily new COVID-19 reported cases, compared to the smooth deterministic counterpart.