Measuring COVID-19 spreading speed through the mean time between infections indicator

Introduction

The COVID-19 outbreak, declared pandemic in 2020, attracted the attention of scientists from different domains (biologists, physicists, engineers and mathematicians, among others). The urgent need to control, predict and monitor the disease progress made it essential to count with mathematical models and algorithms (see for example Cori et al. (2013)), both to manage data on infections and deaths and to perform calculations and predictions. Not only were the well known SIR (Susceptible – Infectious – Removed) compartmental model and its derivations vastly applied (Ali and Khan 2020; Cao et al. 2020; Gleeson et al. 2022; Huang et al. 2021; Kermack and McKendrick 1927; Liu, Zhang, and Wang 2020; Rojas 2020; Simon 2020) but also many new models (Al-Ani 2021) and AI algorithms (Fokas, Dikaios, and Kastis 2021; Gomes da Silva et al. 2020; Silva et al. 2020; Xiong et al. 2020) were proposed. The BPM model, named after Barraza, Pena, and Moreno (2020) and Moreno, Pena, and Barraza (2021), is a stochastic Markovian contagion model described by a nonhomogeneous birth process (NHBP). The authors were able to obtain the functional form of the cumulative infection cases and deaths curves, which have a sub-exponential shape, a behaviour previously pointed out for epidemics (Chowell et al. 2016; Ganyani, Faes, and Hens 2020; Triambak et al. 2021; Viboud, Simonsen, and Chowell 2016). The model has two parameters: one that represents the power of the outbreak and another that models the immunization rate. The expected time between events can be obtained as a standard calculation in counting stochastic processes (Pena, Moreno, and Barraza 2022). Thus, a useful indicator to measure the outbreak speed is obtained: the mean time between infections (MTBI). This indicator expresses the spreading speed, since less time between successive infections implies a more rapid disease spread. Therefore, the MTBI indicator is expected to decrease as the peak gets closer during the initial stage of a disease wave. As we will show, the model parameters taken individually do not suffice to gain insight on the progress of the epidemic, but the MTBI actually does, which makes it a quite robust epidemiological indicator. Thus, the MTBI can be used to evaluate the impact of the actions taken by health institutions.

To calibrate the model parameters, we fit the data of the total reported cumulative infection cases to the mean value function of the process. Since this function does not have concavity changes, the data curve needs to be split into sections (which we will call stages) with a fixed concavity, and then the model can be separately applied to each of them. To achieve this, we look for local maximums and minimums of the daily data curve, which cannot be directly observed due to sharp variations in the reported data. Those jumps are commonly associated to high frequency additive noise in signal processing. Thus, we apply a filtering routine to smooth the daily data curve and then a standard maximum and minimum detection algorithm. Once the stages are properly separated, we fit the model for some specific stages in Argentina, Germany and the United States by the same method described in (Barraza, Pena, and Moreno 2020) to obtain the MTBI indicator, which depends on the model parameters.

An interesting result of the filtering algorithm is the comparison between the filtered curves of daily cases and daily deaths, which shows that both have rather similar shapes, with a short time delay in between. This behaviour was observed in the three considered countries.

In the BPM model the parameter estimation at a given time is performed using the whole history of the corresponding stage. A natural question is whether it is correct to use the full history, or just some part of it. Different approaches to this question were recently presented in the context of the SIR model (Capobianco et al. 2021; Cordelli et al. 2020). All these cases contain some arbitrariness regarding the choice of the data used to calibrate the models. In this work, we perform a comparative study between the possibilities of using different data inputs, for example using the whole stage history or short time windows of a certain number of previous days. We show statistically that the MTBI does not change significantly, whereas the model parameters do. Evolution over time of the MTBI and those parameters is also shown. The most important result obtained is that the MTBI indicator is invariant with respect to the size of the data used for calibration, which lets us conclude that the MTBI is a robust indicator with respect to the data used to estimate it.

Materials and methods

Filtering algorithm

Since the M(t) curve does not have inflection points, the model needs to be applied separately in each stage having either positive or negative concavities, which correspond to the stages (early/mitigation) of each outbreak wave. This implies that a criterion for the wave and stage separation needs to be defined. We choose a rather classic approach:

1. Local minimums in the daily data report represent time instants where a wave finishes and a new wave starts.

2. Local maximums in the daily data report represent time instants where, within a single wave, the initial stage ends and the mitigation stage begins.

The main difficulty regarding these definitions is that the measured data (daily cases and daily deaths) is extremely noisy. This implies that finding local extrema by direct observation is not straightforward. To overcome this problem we performed a filtering procedure to suppress the noise and smooth the curve, followed by a maximum and minimum detection algorithm. In this work we considered Argentina, Germany and the United States as our study cases, but the same results shown here can be achieved for any other dataset by properly setting the filter parameters.

Some interesting observations regarding the noise can be made by looking at the data spectrum. Figure 1 depicts the absolute value of the Discrete Fourier Transform (DFT) of the daily reported cases in the three mentioned countries. Significant peaks can be clearly distinguished at frequency 1/7 and its second harmonic (and its third in the United States case), which is explained by the “once per week” periodic variations in the testing amount on weekends, which are known to be lower. An analogous behaviour can be seen in the deaths dataset. These observations evidence the need of a preprocessing, namely a low-pass filtering, in order to suppress measurement noise. Consequently, we use moving-average (MA) filters, also known as sliding windows, which are the simplest low-pass Finite Impulse Response (FIR) filters; see (Smith 2003) for technical details about MA filters. We proposed to apply n MA filters of length L in series; L and n are the algorithm hyperparameters.

Figure 1:

Spectrum of the daily cases data, as obtained by applying the fast fourier transform (FFT) algorithm. Frequency axis has 1/day units, and is shown up to 1/2, the highest observable frequency (due to the data being sampled at one datum per day). (a) Argentina, (b) Germany and (c) United States.

The reported cases from Argentina, Germany and the United States and their respective filtered curves are shown in Figure 2. Since we want to remove the 1/7 frequency component at least, L must be chosen to be ≥ 7 . In order to avoid harming the signal too much, we chose to fix L=7 and increase n as much as necessary, achieving a nearly Gaussian behaviour (Smith 2003) without suppressing useful signal frequencies.

Figure 2:

Daily cases, reported (purple) and filtered (red). (a) Argentina, (b) Germany and (c) United States.

Having a smooth curve, detecting maximums and minimums can be done by straightforward methods. It must be remarked that the algorithm is not completely automatic: proper hyperparameter values must be carefully chosen by a human user, and since filtering is not perfect, in some cases human interpretation is necessary to determine which minimums and maximums are “false” and which are not.

In Figure 3 we show the filtered curves corresponding to the daily cases and daily deaths in a single plot in order to compare them. It is interesting to note that the cases and deaths curves have similar shapes, with the deaths curve delayed by a short period (approximately two weeks) as expected due to the length of the infectious period. This also shows that the propagation speed of deaths is directly related to the propagation speed of infections. However, due to the existence of measurement noise, they do not always exhibit the same amount of waves.

Figure 3:

Daily cases (blue, scale shown on the left side axis) and deaths (red, scale shown on the right side axis) filtered curves. Filled circles indicate local maximums; × marks indicate local minimums. (a) Argentina, (b) Germany and (c) United States.

Results

Once the data is properly segmented into stages, the disease’s propagation speed can be assessed by analyzing the evolution over time of the MTBI indicator within each of them. To achieve this, we performed the following procedure:

Step 1: Fix a day T=T ₀ to start the procedure.

Step 2: Estimate the model parameters using data from the beginning of the stage and up to T. This yields two estimates γ ̂ ( T ) and ρ ̂ ( T ) of the value of these parameters at day T.

Step 3: Use these estimates on Eq. (4) to obtain M T B I ̂ ( T ) , an estimate of the MTBI on day T.

Step 4: Take T=T ₀+1 and repeat Steps 2–3.

Step 5: Repeat Step 4 until T is the last day of the stage.

This yields three curves γ ̂ ( t ) , ρ ̂ ( t ) and M T B I ̂ ( t ) , which represent estimations of the parameters γ and ρ, and the MTBI, as functions of time, ranging from t=T ₀ to the end of the stage. We typically use T ₀=20 or T ₀=30 if possible to ensure good enough estimates, but it may be necessary to take smaller values in short stages. Figure 4a shows a flow diagram describing this procedure. The fittings were performed by the nonlinear Least-Squares method, using the Levenberg-Marquardt algorithm implemented in the Python language; the source code provided by the authors of Barraza, Pena, and Moreno (2020), where they implemented the BPM model, is available at Pena (2021). The experiments were done using the real (unfiltered) data, since the cumulative curves do not differ much after filtering. Datasets were obtained from Our World In Data (2021).

Figure 4:

Flow diagrams describing estimation and the statistical test. (a) Procedure to obtain the estimated curves γ ̂ ( t ) , ρ ̂ ( t ) and M T B I ̂ ( t ) . (b) Procedure to statistically test MTBI invariance with respect to the window size.

A natural question arises from this study: is the whole history of a stage significant to estimate the parameters, at any given time instant, or should a shorter time period be considered instead? This matter has already been addressed in the context of the SIR model (Capobianco et al. 2021; Cordelli et al. 2020). Then, we intend to answer how much the parameters ρ ̂ ( t ) , γ ̂ ( t ) , and the indicator M T B I ̂ ( t ) change if, at each t, we fit the model using a fixed number W of previous days instead of the whole stage history. We call this a ”time window”, and the whole procedure ”windowing” (the usual terminology in signal processing). Then, for each day t, we considered different window sizes W and calculated the γ ̂ ( t ) , ρ ̂ ( t ) and M T B I ̂ ( t ) curves as previously described (Figure 4a). Our main assertion is that the choice of W does not produce a significant impact on the MTBI estimation, which makes it a quite robust epidemiologic indicator compared to the γ and ρ parameters, which are highly dependant on W. In this section we present statistical evidence supporting this assertion.

For the analysis, we considered the following study cases:

1. Argentina: initial stage of the first wave, from March 03, 2020 to October 17, 2020.

2. Germany: initial stage of the fourth wave, from July 02, 2021 to November 08, 2021.

3. United States: initial stage of the third wave, from September 07, 2020 to December 30, 2020.

As it can be seen in Figure 5, the obtained ρ ̂ ( t ) and γ ̂ ( t ) / ρ ̂ ( t ) curves are indeed very different according to the value of W used in each case. The resulting ρ ̂ ( t ) curves change in both shape and scale (i.e., the numerical values have different orders of magnitude, which is why we plot each curve on its own y axis), as it can be seen in Figure 5a–c. The γ ̂ ( t ) / ρ ̂ ( t ) curve changes in shape, but the numerical values do not differ much (the pairs of curves in Figure 5d–f are plotted in the same scale); this implies that even if γ ̂ ( t ) and ρ ̂ ( t ) change individually between the two methods, they do it proportionally. In any case, the estimation of either ρ ̂ ( t ) or γ ̂ ( t ) / ρ ̂ ( t ) is severely impacted by the choice of W, i.e., to the size of the input data.

Figure 5:

Evolution over time of ρ ̂ ( t ) (upper row) and γ ̂ ( t ) / ρ ̂ ( t ) (lower row) parameters for each of the three countries in its chosen stage. Green curves were estimated using the whole stage history. Red curves were obtained using a time window of size W=30. In the case of the ρ ̂ ( t ) estimations (upper row), numerical values attain different orders of magnitude on each case, which is why the green curve scale is shown in the left side axis whereas the red curve scale is shown in the right side axis. The curves being so different show explicitly how highly the values of γ and ρ depend on the choice of W. (a) Argentina ρ, (b) Germany ρ, (c) United States ρ, (d) Argentina γ/ρ, (e) Germany γ/ρ and (f) United States γ/ρ.

Now we turn to the analysis of the MTBI indicator. In contrast to the behaviour of the ρ and γ/ρ parameters, the M T B I ̂ ( t ) did not change much when estimating with different sizes of window W. This means that the MTBI captures the epidemic behaviour more efficiently than any of the individual parameters, and that it can be evaluated without risk even if the available data is short. The rest of this section is dedicated to provide statistical support to this assertion. To do this, we propose the following procedure (see Figure 4b).

Step 1: Choose n different window sizes W _i .

Step 2: Compute the M T B I i ̂ ( t ) curve for each window W _i . This yields n different estimators of MTBI(t).

Step 3: For each index (day) t _j in the curve, form a sample set with every M T B I i ̂ ( t j ) . This yields a total of L sets S _j of size n, where L is the length of the curves.

Step 4: Test normality on each S _j .

Step 5: Compute the average S j ̄ for each j as a point estimator of MTBI(t _j ).

Step 6: Use each S _j to produce confidence intervals around S j ̄ .

This algorithm can be interpreted as follows: produce n different estimators for the MTBI(t) curve and test for normality pointwise. Then, take the average of all the estimator as a ”best guess” of the real MTBI(t) and measure its precision via confidence intervals. The normality of the samples, if achieved, means that observations are disposed around its mean in a Gaussian manner and its differences are mostly due to random estimation noise. Figure 6a shows the results of applying a Kolmogorov-Smirnov test with 0.05 significance level: the green points indicate the indexes (days) where no evidence against normality can be found. The test p-values for the three countries are plotted in Figure 6b. In all three cases there are some days when the test fails (more precisely: 40 out of 185 (21.6%) for Argentina, three out of 86 (3.5%) for Germany and 11 out of 71 (15.5%) for USA), but it is clear than most of the samples can be safely considered normal.

Figure 6:

Results of the 95% Kolmogorov-Smirnov normality test applied pointwise to the sample of M T B I ̂ ( t ) estimated curves with W ranging from 15 to 45. The x-axis shows the index of the day in the curve (recall that the three stages have different lengths). The left figure shows in which indexes normality can be rejected (red dots) and in which the corresponding sample is normal (green dots). The right figure shows the pointwise p-values of the test as well as the 0.05 rejection threshold. (a) Test results and (b) test P-values.

The last step is to show that the estimated M T B I ̂ i ( t ) curves are ”sufficiently close to each other” to be considered identical, except for random measurement noise. In order to do this, we consider the average curve M T B I ̄ ( t ) , which is itself a point estimation of MTBI(t). At each day t=t _j we have a Gaussian sample S _j with sample mean S j ̄ = M T B I ̄ ( t j ) . Then, a straightforward confidence interval can be built around S j ̄ using Student’s t distribution. The already proven normality turns this into a robust procedure, but since we used quite large samples (n=31) it is not a strict requirement, so this analysis is valid even in the points where the Kolmogorov-Smirnov test allowed to reject normality. The results are shown in Figure 7 for the three analysed countries in its chosen stage. We plotted M T B I ̄ ( t ) (solid blue line) and two confidence intervals of 95 and 99% (dashed green and red lines respectively). The light blue shadow represents the standard deviation of the sample. Simple inspection of these plots shows a remarkable fact: the red dashed line is located around half way between the average curve and the outer edge of the blue shadow (this is especially clear in the Germany curve and the first half of the USA curve). Since that outer edge is M T B I ̄ ( t j ) + σ ( t j ) , the location of the red dashed line can be interpreted as follows: we have 99% confidence that the ”unknown” MTBI(t) does not differ from M T B I ̄ ( t ) in more than a half standard deviation. The intervals being so small is notorious evidence that all the M T B I i ̂ ( t ) estimates are extremely close to each other. It is also important to remark that all the estimates come from a well fit model, achieving R² coefficients of over 0.97 in all cases. All the facts previously exposed allows us to reach the conclusion that the MTBI values does not change significantly with respect to the size of the input data used to fit the model. The analyzed samples suggest that time windows of size W=15 can be safely used and that the ”far” history of the stage has no severe impact on the estimations.

Figure 7:

Statistical analysis of M T B I ̂ ( t ) . The solid blue curves are the average curves, M T B I ̄ ( t ) , of 31 estimations made using different window sizes. Red and green dashed lines are two pairs of confidence intervals around this average. The light blue shadow indicates the pointwise standard deviation within the 31 averaged curves. In the case of Argentina, the image is shown in a log scale to cover the high range of values. (a) Argentina, (b) Germany and (c) United States.

Discussion/conclusions

In this work we have analyzed the evolution over time of the MTBI indicator given by a recently proposed stochastic nonhomogeneous Markov model when applied to the COVID-19 pandemic, which is useful to measure the propagation speed of the epidemic. Three countries were considered as study cases: Argentina, Germany and the United States. For each of these countries, a particular stage was chosen and the model was fitted using different sets of input data. On the one hand, the full history of the stage until the day t was considered to calibrate the model parameters at that day, and on the other hand, the calibration at t was made using only the data from several time windows of different lengths. We show the evolution over time of the averaged MTBI indicator, as well as the evolution of the model parameters. The most important conclusion was that the MTBI indicator (which depends on those parameters) does not change significantly using different sizes of input data, which shows this indicator is robust with respect to the dataset size. Of all the three study cases, only the one from Germany has a vaccination process, rising from 35 to 50% over that period. The consequently sudden slow down in the speed propagation can be seen in the picture. Due to the population size and social mobility (between 20 and 30% under normal) in the considered stage of the United States, the propagation speed is much higher than in the other two study cases, which is reflected in the extremely low numerical values of the MTBI (Figure 7c). In the case of Argentina, the strict lockdown imposed during the first wave of the epidemic is reflected on high M T B I ̂ ( t ) values at the beginning of the first wave (i.e., a slow propagation). As future work, we propose to study the predictive power of the MTBI indicator.