Decision support for the quickest detection of critical COVID-19 phases

During the course of an epidemic, one of the most challenging tasks for authorities is to decide what kind of restrictive measures to introduce and when these should be enforced. In order to take informed decisions in a fully rational manner, the onset of a critical regime, characterized by an exponential growth of the contagion, must be identified as quickly as possible. Providing rigorous quantitative tools to detect such an onset represents an important contribution from the scientific community to proactively support the political decision makers. In this paper, leveraging the quickest detection theory, we propose a mathematical model of the COVID-19 pandemic evolution and develop decision tools to rapidly detect the passage from a controlled regime to a critical one. A new sequential test—referred to as MAST (mean-agnostic sequential test)—is presented, and demonstrated on publicly available COVID-19 infection data from different countries. Then, the performance of MAST is investigated for the second pandemic wave, showing an effective trade-off between average decision delay \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta$$\end{document}Δ and risk \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document}R, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document}R is inversely proportional to the time required to declare the need to take unnecessary restrictive measures. To quantify risk, in this paper we adopt as its proxy the average occurrence rate of false alarms, in that a false alarm risks unnecessary social and economic disruption. Ideally, the decision mechanism should react as quick as possible for a given level of risk. We find that all the countries share the same behaviour in terms of quickest detection, specifically the risk scales exponentially with the delay, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \sim \exp {(-\omega \Delta )}$$\end{document}R∼exp(-ωΔ), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document}ω depends on the specific nation. For a reasonably small risk level, say, one possibility in ten thousand (i.e., unmotivated implementation of countermeasures every 27 years, on the average), the proposed algorithm detects the onset of the critical regime with delay between a few days to 3 weeks, much earlier than when the exponential growth becomes evident. Strictly from the quickest-detection perspective adopted in this paper, it turns out that countermeasures against the second epidemic wave have not always been taken in a timely manner. The developed tool can be used to support decisions at different geographic scales (regions, cities, local areas, etc.), levels of risk, instantiations of controlled/critical regime, and is general enough to be applied to different pandemic time-series. Additional analysis and applications of MAST are made available on a dedicated website.

with initial conditions r(0) = 0, s(0) = 1 − i(0); i(0) is a (fairly small) fraction of the total population that gives rise to the spread of the infection. Note that one of the three equations in (1) is redundant, in view of the conservation constraint r(t) + s(t) + i(t) = 1. Dividing the first of (1) by s(t) and solving for i(t) from the third equation, one gets d log s(t) = − β γ dr(t). The quantity β /γ goes also under the name of contact number and is a combined characteristic of the population and the disease. Also, using i(t) = 1 − r(t) − s(t) in the third equation, we eliminate i(t) and arrive at    s(t) = s(0)e − β γ r(t) , to be solved numerically. At the onset of the epidemic, i.e., for t → 0, the fractions of susceptible and recovered are approximately constant and equal to s(t) ≈ 1 and r(t) ≈ 0, respectively. As t → ∞, if a steady-state regime for r(t) emerges, we expect dr(t) dt = 0. Imposing this condition and denoting by r(∞) = lim t→∞ r(t) and s(∞) = lim t→∞ s(t), from the second of (2), we have r(∞) = 1 − s(0)e − β γ r(∞) ≈ 1 − e − β γ r(∞) , which can be numerically solved for r (∞). An analytical solution is also available, in terms of the so-called Lambert W function (also known as product logarithm): 4 From the first of (2), we see that the fractions of recovered and susceptible individuals reach steady-state constant values r(∞) (typically > 0) and s(∞) (< 1) verifying r( For instance, this happens at the beginning of the epidemic with s(t ⋆ ) = s(0) ≈ 1, and at the end of the epidemic, with s(t ⋆ ) = s(∞) < 1. Thus, at different stages of the epidemic, s(t) can be approximately assumed constant over short intervals of time. We make the assumption that variations of s(t) can indeed be neglected over short intervals of time. With s(t) = s(t ⋆ ), from the second equation in (1), we get initially grows exponentially at rate β − γ and eventually decreases exponentially to zero at rate β s(∞) − γ.

Proposed Model
Epidemic data are typically collected on a daily basis, and a discrete version of (4) with discretization step ∆t = 1 day can be obtained as follows. Let n be the day index. With the obvious notation i n = i(n∆t) and similarly for other quantities, the differential equation on the left-hand side of (4) is approximated by for some i 0 > 0, and therefore the ratio (also known as growth, for α > 0, or decay, for α < 0, rate) is constant and equal to (1 + α). It is useful to bear in mind that we typically have |α| ≪ 1. From (5) we see that the sequence i n grows or decays exponentially according to the sign of α, and its evolution can be locally assumed to be linear in n because (1 + α) n ≈ 1 + n α, in view of |α| ≪ 1. For later use, note that, from the third equation in (1), we also get ∆r n :=r n+1 − r n = γ i n .
It is evident that real-world data lead to "noisy" versions of the previous expressions. There exist two main sources of uncertainty. The first is implicit in the nature of the contagion, which obviously depends on a number of medical and social factors that are problematic to model in deterministic terms. Accordingly, the sequence {x n } n≥1 in (6) is better represented by a collection of random variables, leading to the following model for the fraction of infected individuals on day n: In (8), we assume that {x n } n≥1 is a sequence of nonnegative independent random variables with mean close to unity due to the fact that |α| ≪ 1. The independence assumption is crucial to ensure mathematical tractability; however, as we show in the section on "Performance Assessment -Synthetic Data," slight deviations from the condition of perfect independence do not significantly affect the performance of the proposed Mean-Agnostic Sequential Test (MAST). The second source of uncertainty is related to the way in which the epidemic data are collected, as we describe next. First, let us consider the values of "new positives" recorded on day n and its relationship to the "total cases" on days n and n − 1: p n := new positives n = total cases n -total cases n−1 .
Since the number of total cases on day n is equal to i n + r n , we have p n = i n + r n − (i n−1 + r n−1 ) = ∆i n−1 + ∆r n−1 = (α + γ)i n−1 , because of (5) and (7). Then, we get  Table S1. p-value of Kolmogorov-Smirnov goodness-of-fit test, and estimated value of σ . yielding for some initial value p 1 . This is important, as while there is a natural modeling in terms of the proportion of infected individuals i n , the time-sequence {i n } is difficult to obtain as the number that get removed from it (by recovery or death) is often not reported. On the other hand, the number of newly infected p n is easily available, and reported with vigor. Based on (11), the sequence {x n } is obtained from {p n } using the relationship x n = p n+1 /p n . Returning to the problem of the second source of uncertainty, the sequence of daily new positives that is available from the John Hopkins University 5 is expected to differ from the actual numbers of new positives {p n } n≥1 due to gross errors, missing values or delays in the reported data. For instance, it may happen that, on a given day, the fraction of new positives is not reported (some peripheral data collection unit did not communicate with the central unit in a timely manner) and such unreported fraction is then added to the value observed on a later day. Also, it typically happens that the number of reported cases over the weekend is systematically smaller than the values recorded on other days of the week, due to the smaller number of swabs tested during weekends. Similar periodic modulation effects are already reported in the literature. 3 For these reasons, we first eliminate from the data anomalous values (negative numbers) and then, in place of (11), we consider a more sensible observation model: where w n is a noise term taking into account the effect of gross errors. Due to their nature, we expect that gross errors are washed away by averaging the data over a time window of length L, where L is of the order of a week or longer. On the other hand, we have seen that the expected behavior of p n over a sufficiently small time window is approximately linear. For this reason, the downloaded sequence of new positives is first filtered by a moving average filter of length L days (in the paper we use L = 21, see below). Under these assumptions, we conclude from (12) that: where p n L = 1 L ∑ k p k denotes the L-th order uniformly-weighted moving average [MA(L), for short] of the sequence {p n }, where the sum involves L entries centered on the n-th sample. Admittedly, the approximation in (13) is rather sharp, but works for our purposes. Summarizing, the sequence of daily new positives (after removing negative entries) is first processed by an arithmetic moving average filter of order L, yielding {p n } n≥1 and then, from this, the sequence {x n = p n+1 /p n } is obtained as in (10).

Statistical Characterization
Recall that the subscript n to x n and other quantities denotes the index of the day. The statistical distribution of the sequence {x n } is derived from COVID-19 data made available by the Johns Hopkins University, 5 Figure S1. Example of the construction of the periodic counterparts of { µ 0,n } and { µ 1,n } from the Italian data, used to generate the synthetic grow rates under the controlled and the critical regime, respectively. From the growth rate sequence x n (green curve), the estimated mean sequence { µ n } (magenta) is obtained through a uniformly-weighted moving average filter of length 21 days. From { µ n }, the two sub-sequences { µ 0,n } (dark yellow) and { µ 1,n } (dark cyan) are extracted, so as to verify µ 0,n 1 and µ 1,n > 1. The two sub-sequences are then replicated, with the odd replicas flipped to preserve continuity, yielding the dark yellow and dark cyan solid curves, which are finally used as mean values to generate the synthetic data under the controlled and the critical regime, respectively. One realization of the resulting synthetic growth rates under H 0 and under H 1 is superimposed to the mean values, for illustration purpose.
let us refer to the data from Italy. Figure 2 Fig. 2(b) in green. Numerical investigations, not detailed for the sake of brevity, show that the growth rate x n can be modeled by a Gamma random variable with shape parameter κ and scale parameter θ , with κθ close to unity. Since κ is on the order of several hundreds, the Gamma distribution is practically indistinguishable from a Gaussian distribution with mean µ = κθ on the order of unity, and standard deviation σ = √ κθ ≪ 1. In light of these considerations, we adopt the Gaussian model because it leads to a detector with simple and intuitive structure.
Accordingly, let us reconsider the sequence {x n } shown in green in Figure 2(b). First, we disregard the initial portion of the sequence {x n } that corresponds to the first peak of {p n } (first wave of the epidemic). The omitted data corresponds to values of the day index n smaller than the index of the first passage (from the above) by one of the sequence {x n }. Then, we consider the smoothed version { µ n } of {x n }, shown in magenta in Figure 2(b). The smoothing is obtained by processing {x n } again by a MA(21) filter. The length L = 21 of the two MA filters used in this analysis have been selected by trial and error. Note that the filters are not causal: 6 to compute the output at time n, (L − 1)/2 successive samples of the input are needed. Note also that the filters are truncated at the endpoints where less than L samples are available.
By subtracting the smoothed version { µ n } from the sequence {x n }, a sequence {x n − µ n } with approximately zero-mean is obtained. Then, we verify that the random variables x n − µ n , n = 1, 2, . . . , have one and the same standard deviation σ ≪ 1, with good accuracy. In particular, there is no evidence of a substantial change in standard deviation between the regions with µ n ≤ 1 and those in which µ n > 1. The estimated values of σ for the 14 nations considered in the main document is reported in Table S1. To distinguish between these two regions, we introduce the notations µ 0,n :=µ n ≤ 1 and µ 1,n :=µ n > 1 for the actual (unknown) mean values. Note that, since µ 0,n and µ 1,n are close to unity and σ ≪ 1, we see that x k < 0 with negligible probability.
The quantity {x n − µ n } so obtained is modeled as a sequence of zero-mean independent Gaussian random variables with (nation-dependent but small) standard deviation σ . The independence assumption is necessary to ensure mathematical tractability and, in practice, slight deviations from the condition of perfect independence do not significantly affect our results. In this respect, it should be also noted that the smoothing operation shown in Eq. (13) enforces correlation over the sequence {p n }, and, in turn, over the sequence of ratios {x n = p n+1 /p n }. Finally, in Table S1, we verify that the p-value obtained by the Kolmogorov-Smirnov (KS) goodness-of-fit test 7 is larger than 0.01, meaning that the Gaussian assumption cannot be rejected at 1% significance level, for almost all the countries (except Spain and UK). We thus arrive at the following model of the observable growth rate sequence {x n }: x n ∼ N (µ 0,n , σ ) with µ 0,n ≤ 1, under the controlled regime, and x n ∼ N (µ 1,n , σ ) with µ 1,n > 1, under the critical one.

Real Data
The implementation of the MAST algorithm does not require knowledge of the mean sequences {µ 0,n } and {µ 1,n } and, in fact, it has been designed just to cope with this uncertainty. However, the performance of the test depends, of course, on the specific scenario under consideration. Performance of MAST run on publicly available COVID-19 infection data from different countries is reported in the main document and obtained as follows. For each country, we investigate the performance of MAST by using the estimates { µ 0,n } and { µ 1,n } of the mean sequences {µ 0,n } and {µ 1,n }. Since we need to simulate arbitrarily long data streams under both the controlled and critical regimes, we construct periodic counterparts of the sequences { µ 0,n } and { µ 1,n }, as illustrated in Fig. S1. Using these periodic counterparts, we obtain the performance of the test by standard Monte Carlo computer experiments, 2, 8, 9 involving 10 5 independent runs for each value of the threshold χ. This gives the relationship between R and χ, and between ∆ and χ, for relatively small values of χ. The former relationship is almost exactly exponential, and the latter almost exactly linear. Thus, the aforementioned relationships are extrapolated to arbitrarily large values of χ. This explains why the figures of the main document report values of R that would be difficult to obtain by standard computer experiments.

Synthetic Data
The performance of MAST is also evaluated in a synthetic scenario and compared with those of a naive Page's test. To emulate a realistic scenario, we assume that the mean sequence {µ n } has a periodic evolution not exceeding 1 in the controlled regime, and not falling below 1 in the critical regime. Specifically, we set  Figure S3. Performance of MAST and the naive Page's test in terms of risk versus mean delay. Red dashed lines and blue solid lines refer to the naive Page's test and MAST, respectively, when run on uncorrelated sequences; yellow dash-dotted lines and green dotted lines refer to the naive Paige's test and MAST, respectively, when applied on correlated sequences. The markers allow to distinguish between different values of σ : circles for σ = 0.065, squares for σ = 0.050, and triangles for σ = 0.035. and µ 1,n = 1 + ε 2 cos(2πnM −1 + φ 1 ) + 1 . Fig. S2(a) shows the mean sequences used for this analysis, with ε = 0.1, period M = 75 days, and phases φ 0 and φ 1 uniformly drawn from [0, 2π); in addition, one realization of the resulting synthetic growth rates under H 0 and under H 1 is provided, obtained by adding to each mean sequence a zero-mean Gaussian noise process with standard deviation σ = 0.050. Moreover, in order to numerically demonstrate that slight deviations from the condition of perfect independence of the sequence {x n − µ n } do not significantly affect the performance of the MAST (as claimed above), we also consider the case of a correlated zero-mean Gaussian noise process. The (normalized) correlation we enforce on the sequence {x n − µ n } is shown in blue in Fig. S2(b); this is obtained by averaging the correlations computed from the sequences {x n − µ n } of the 14 nations considered in the main document. The correlated Gaussian process is generated by filtering an uncorrelated Gaussian sequence as described in [10,Ch. 11.3.2].
MAST is compared with a naive Page's test that is aware of the bounds of the mean sequence under each hypothesis, namely, the lower bound 1 − ε and the upper bound 1 under H 0 , and the lower bound 1 and the upper bound 1 + ε under H 1 ; naive here implies that the test does not know the exact evolution of the mean sequence, but only its extreme values. The naive Page's test is based on Eq. (6), reported in the main document, imposing α = ε, and corresponds to the standard Page's test when M = 1, φ 0 = π, φ 1 = 0, and the data are uncorrelated. Fig. S3 presents the performance of MAST and the naive Page's test in terms of risk versus mean delay for different values of σ , obtained by standard Monte Carlo computer experiments 2, 8, 9 involving 10 5 independent runs for each value of the threshold χ. We observe that in the case of uncorrelated sequences, MAST outperforms the naive Page's test, even though the former does not have any prior information about the mean sequences, except for them being above or below 1, that is, µ 0,n 1 and µ 1,n > 1. By enforcing the correlation, the performance of both MAST and the naive Page's test improves; however, there are no remarkable differences, with the displacement in terms of mean delay being no more than 1 day for the MAST and 2 days for the naive Page's test at a risk level of 10 −4 .