key: cord-0628620-elfueecx
authors: Braca, Paolo; Gaglione, Domenico; Marano, Stefano; Millefiori, Leonardo M.; Willett, Peter; Pattipati, Krishna
title: Quickest Detection of Critical COVID-19 Phases: When Should Restrictive Measures Be Taken?
date: 2020-11-23
journal: nan
DOI: nan
sha: b460b23da7db50d8382634f9f83955aea86fd1e2
doc_id: 628620
cord_uid: elfueecx

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 $Delta$ and risk $R$, where $R$ is inversely proportional to the time required to declare the need to take unnecessary restrictive measures. We find that all the countries share the same behaviour in terms of quickest detection, specifically the risk scales exponentially with the delay, $R sim exp{(-omega Delta)}$, where $omega$ 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 three weeks, much earlier than when the exponential growth becomes evident.

With more than 57 million cases worldwide and over one and 1.36 million deaths as of November 20, 2020, the outbreak of coronavirus disease (COVID- 19) is undeniably one of the worst global crises since World War II. In March 2020, the exponential increase of individuals needing hospitalization in intensive care units, combined with the lack of effective cures and vaccines, pushed many governments to take extraordinary measures aimed at "flattening the curve" of infections. 1-3 The adopted measures included the limitation of mobility and social activities, closure of schools, universities, shops, factories, and so forth, up to the extreme act of national lockdowns. Evidence that such measures achieved a reduction of the rate of new infections are gradually appearing in the scientific literature. [4] [5] [6] These measures contributed to keep the spread of COVID-19 under control for some time, but we are now, in November 2020, experiencing the onset of a new exponential growth of confirmed cases, the "second wave," with severe risks for personal health and healthcare systems under severe stress. Governments and authorities are facing again the difficult task of deciding if and when new containment measures may be needed. In absence of limitations to mobility and social activities, the pandemic spreads exponentially in time, so that any delay in applying restrictions may lead to severe consequences. On the other hand, accounting for the social and economic impact of the possible countermeasures, already observed in the first half of 2020, 7-10 restrictions should be taken only if and when it is strictly necessary. Managing the trade-off between these contrasting requirements is extremely challenging.

To address this challenge, we leverage sequential detection theory, 11-13 and specifically quickest detection schemes [14] [15] [16] to propose a rigorous methodology aimed at identifying as quickly as possible the onset of an exponential growth of the pandemic Each curve is normalized by the initial value and shifted to the same initial time. The initial point of the first wave for a specific nation is when the first positives are reported, while the initial point of the second wave is when the daily positives start to increase again -growth rate larger than one. The first wave, upper region (shaded in red), was more aggressive than the second wave, mostly confined to the lower region (shaded in blue). These two regions are delimited by different values of α in Eq. (1). Curves of each nation are dashed if related to the first wave otherwise solid if related to the second wave.

evolution. The proposed procedure -referred to as MAST (mean-agnostic sequential test) -is designed to minimize the average time to detect a change in regime: [14] [15] [16] from a situation in which the pandemic is under control (H 0 regime) to the onset of an exponential growth (H 1 regime).

Quickest detection theory has a long history, 17 and has been successfully applied in several fields, including quality control, climate modeling, remote sensing, financial analysis, image analysis, security, signal and speech processing, and biomedical applications. 15, 16 In the context of public health surveillance, where the timely detection of various types of adverse health events is crucial, quickest detection techniques have found several applications, 18 e.g., to reveal the onset and the peak of the epidemic period 19 or that the peak is over. 20 The approach pursued in this article is substantially different from most of the epidemiological models, based for instance on stochastic evolution of epidemic compartments 4, 6, 21-23 and metapopulation networks, 24, 25 where the goal is to predict the mid/long-term behavior of the outbreak. For instance, in stochastic compartmental models, given an initial condition, the epidemic can have two outcomes: the number of infected individuals can increase, in which case we have a major outbreak, or decrease. The probability of a major outbreak can be computed, 26 but it is of limited use in taking timely on-line decisions.

The importance of taking a decision as quickly as possible in an epidemic scenario can be understood by looking at the curve of daily cases of COVID-19 infection, reported in Fig. 1 for several nations. 27 In the same figure we also report the deterministic curves of daily cases, with an exponential growth described by the following equation p n+1 = p n (1 + α) = p 1 (1 + α) n , n = 1, 2, . . . ,

where n is the time index (day), 1 + α is the growth rate, and α > 0, which corresponds to the H 1 regime of a major outbreak. This exponential behavior is equivalent to a recently-proposed disease-transmission model, 28 and the reproduction number defined therein is equivalent to the growth rate. All the curves in Fig. 1 are normalized to the initial value p 1 and shifted to the same initial time. Figure 1 shows that first wave was noticeably more aggressive than the second one in terms of growth rate; specifically, in the first wave α was varied by country but ranged between 0.06 and 0.40, while in the second wave it is between 0.01 and 0.06. In the latter scenario, any delay in revealing the onset of an exponential growth phase and the consequent implementation of containment measures produces a costly exponential increase of the number of new cases. In other words, it is essential to reduce as much as possible the decision delay to level off the curve of infection as early as possible. At the same time, it is important to enforce restrictions only if essential, in order to avoid unnecessary social unrest and economic cost. And, inevitably, any detection procedure under the controlled regime H 0 can produce false alarms, i.e., it can wrongly declare the upcoming onset of an exponential growth phase. The risk of taking a wrong decision can be quantified by the inverse of the mean time between false alarms. Balancing between detection delay and risk represents a fundamental system trade-off. In this paper, we consider a relatively simple, but effective, mathematical model of the pandemic and develop a decision tool to quickly detect the passage from a controlled regime to a critical one. Its effectiveness in terms of delay/risk trade-off is demonstrated on publicly available COVID-19 data from several countries. It is our hope that the proposed MAST procedure can be useful in making timely and rational decisions to control the pandemic evolution.

The main result of this research is the design of the MAST quickest-detection procedure, which is specifically tailored to epidemic scenarios, and its application on COVID-19 data. In this context, we show that the mean time ∆ required to reveal the onset of an exponential phase is in the order of a few days (or few weeks), with a risk R that scales exponentially with the delay:

where the symbol "∼" refers to asymptotic behaviour (small R, large ∆), and ω is a parameter whose value varies from country to country. This exponential relationship holds for most optimal quickest detection procedures under specific mathematical conditions and optimality criteria. 14 In this work, we show empirically that this optimality condition is verified by the MAST procedure when applied to the COVID-19 data, despite the fact that MAST is designed to work in the presence of uncertainty and non-stationarity of the data statistics.

To illustrate these concepts, let us refer to Fig. 2 , where data from Italy are considered. Panel (a) shows the curve of daily positives from Johns Hopkins University COVID-19 data, 27 along with the smoothed version thereof, shown in green. Using the latter, we obtain the sequence of growth rates x 1 , x 2 , . . . , shown in panel (b) (green curve) as the ratio of two consecutive values of daily positives, with its smoothed version (in magenta); see the data model in Eq. (3).

In Fig. 2 (a), we note that the peak of the first wave is smaller than the second one, which is likely due to the smaller number of swabs made during the first wave. In other words, in the two waves we are observing different fractions of the total infected populations. This difference does not affect our analysis because the sequence x 1 , x 2 , . . . , of growth rates are ratios of successive daily positives.

As discussed in the Methods section and illustrated in the notional scheme of Fig. 5 , the proposed MAST procedure takes as input the sequence x 1 , x 2 , . . . and recursively computes the sequence of decision statistics T 1 , T 2 , . . . , see Eq. (7), where T n only depends on the observed growth rates x 1 , . . . , x n , up to day n. It is worth stressing that the growth rate process x 1 , . . . , x n , is assumed to be Gaussian distributed with unknown and time-varying mean value sequence. The MAST decision statistic is "mean-agnostic," because it does not require knowledge of the time evolution of such a mean sequence, and it is therefore robust to its deterministic fluctuations. . Application of the MAST procedure on the COVID-19 pandemic data from Italy. On the left-side vertical axis, we select a decision threshold to correspond to a desired risk, e.g., R = 10 −4 . Then, the blue curve indicates the stopping day (about July 18, in the example) corresponding to the selected value of risk. Finally, the red curve referred to the right-side vertical axis shows the mean delay ∆ corresponding to the selected risk R (about 3 days). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range.

The onset of the critical regime is declared by the MAST at the smallest index n such that T n > χ, where χ is a threshold level. Panel (c) of Fig. 2 shows the corresponding decision performance. The left axis reports the risk R as a function of the threshold value χ, and the right axis reports the mean delay of decision ∆, again as a function of χ. We see that the function R(χ) is exponentially decreasing, while ∆(χ) is linearly increasing, consistent with Eq. (2) and with known relationships for Page's test and other quickest-detection procedures. 14, 15 Combining the evidence shown in Fig. 2 , we obtain Fig. 3 that can be interpreted as follows. The blue curve refers to the left-side vertical axis, while the red curve refers to the right-side vertical axis. Starting from a given level of risk selected on the left-side vertical axis, we find the stopping day of MAST, which is the day on which the onset of the critical regime is declared (first threshold crossing). Then, referring to the red curve, one obtains the corresponding delay ∆ on the right-side vertical axis. This is the mean delay incurred by the MAST procedure in declaring the onset of the critical regime. The interpretation of ∆ is that the passage to the exponential growth of the pandemic takes place, on the average, ∆ days before the alert produced by MAST. In Fig. 3 , values of the risk smaller than 10 −9 are collapsed to R = 10 −9 , because for R ≤ 10 −9 we can safely assume that risk is essentially negligible. It is worth noting that, even at negligible risk level, the stopping day is approximately July 27, and the corresponding mean delay is less than 8 days.

The analysis in Figs. 2-3 can be repeated for other regions, yielding similar insights. We report the cases of: United States of America in Fig. 7 ; United Kingdom in Fig. 8 ; France in Fig. 9 ; and Germany in Fig. 10 . Note that in the United States, there exists also a third wave of the pandemic, which we analyze separately in Fig. 7(d) , by restarting the test after declaring the onset of the second wave. (Actually, such a third wave is likely to be a delayed second wave in different geographic regions of the U.S., as a state-by-state analysis seems to imply.) Additional analysis on data from more geographical regions is available online 29 and updated regularly.

A comparison of the MAST performance for 14 different nations is addressed in Fig. 4 . Different decision performances reflect different values of the parameter ω governing the relationship between risk and delay, shown in Eq. (2) . We see that, accepting a risk of R = 10 −4 , the mean detection delay ∆ is about 3 days for Netherlands and below 20 days for Spain. Thus, the MAST procedure detects an outbreak very expeditiously at very low risk level.

The outer black curves labeled by ω = 0.32 and ω = 11.52 represent the ideal performance of the Page's test (see Eq. (6)), in which the mean values of the growth rates x 1 , x 2 , . . . , are assumed constant, known in advance, and equal to (1 + α) under H 1 and (1 − α) under H 0 . For such a scenario, the parameter ω appearing in Eq. (2) is the Kullback-Leibler distance 2 (α/σ ) 2 between the distributions under the two regimes. 14 Figure 4 . Operational curve -risk versus mean delay for decision -for 14 Countries. For large ∆, the operational curve is described by Eq. (2), namely R ∼ exp(−ω∆). The outermost black curves correspond to the ideal performance of the Page's test, assuming known and constant growth rate, with α = 0.01 and 0.06, respectively (extreme values of α for the second wave, see Fig. 1 ). Each nation is characterized by a specific value of ω (reported between brackets on the legend), and all values of ω lie in the range between ω = 0.32 and ω = 11.52. and α = 0.06 corresponding to the extreme growth rates of the second wave reported in Fig. 1 , and set σ = 0.025, which is the arithmetic mean of the estimated standard deviations of the 14 countries. For the two values of α, this yields 2 (α/σ ) 2 = 0.32 and 11.52, respectively, which are the values reported in Fig. 4 . By setting ω = 2 (α/σ ) 2 in Eq. (2) we observe that the higher is the growth rate, the better is the detection capability for a given risk. In other words, consistent with intuition, the more aggressive is the outbreak, the more quickly it can be detected. On the other hand, the larger is σ , the higher is the delay for a given risk. Also this effect is intuitive, because the standard deviation σ measures the entity of random fluctuations in the data x 1 , x 2 , . . . , and reliable decisions require more time because the data is uncertain. These trade-offs are also observed for the MAST test run over COVID-19 data for different nations, and are captured by the parameter ω in Eq. (2).

Let us focus again on Italy, and set R = 10 −4 . Accepting a risk in the order of R = 10 −4 means that there is one chance in ten thousand that the countermeasures are taken too early or, in other words, unjustified actions against the pandemic are adopted every 27 years, on the average. By taking drastic countermeasures on July 18 -the stopping day prescribed by our MAST procedure for R = 10 −4 -we would have left only ∆ ≈ 3 days of uncontrolled exponential growth of the pandemic, before addressing it. In this respect, it should be noted that the adoption of severe countermeasures in Italy has been decided only at the beginning of November 2020. The delayed decision situation is analogous for other nations.

It is clear that managing an unprecedented pandemic event is a huge and a multifaceted problem, which can only be addressed by taking into account many different perspectives. It is also clear that the decision when to take pandemic countermeasures depends on a large number of societal factors. The contribution of this article is limited to the analysis of the pandemic strictly from a quickest-detection viewpoint and, from this perspective, we obtain useful insights and quantitative analyses. One evidence, as just pointed out, is that critical regimes of many nations began dramatically earlier than when countermeasures were taken.

Aside from the above retrospective analysis, a major contribution of the MAST quickest detection tool developed in this paper consists of providing proactive decision support for detecting future waves of the COVID-19 outbreak, and the onset of future pandemics. In these cases, precise nation-dependent performance predictions, such as those given in Fig. 4 , cannot be available in advance because the curves have been obtained by exploiting estimates of the mean value of the growth rate (see e.g., the curve in magenta in Fig 2) , derived by forensic inspection of the data. However, good approximate performance bounds can be obtained by assuming constant growth rates (1 ± α) in the critical and controlled regime, respectively, for "reasonable" values of α, for instance, the faster and the slower rates reported in Fig. 2 and used to compute the outermost black curves in Fig. 4 Figure 5 . Flowchart of the proposed MAST procedure. Input data are the daily infected. These noisy data are filtered to mitigate imperfections in data collection, randomness, and delays. Filtered daily positive p n are used to compute the growth rate x n , which is used to compute the MAST statistic T n . The statistic is then compared with the threshold χ; if it is larger than the threshold, the outbreak is declared. Otherwise, the procedure continues to collect and process the data.

governs the performance of a clairvoyant Page's test that knows exactly the constant growth rate. Figure 5 shows the flowchart of the MAST procedure. The filtering operation in Fig. 5 is important to mitigate gross errors and lack or delayed reporting of the input data. To compute p n , the filtering operation used in this paper requires one to observe data samples beyond the current day n, which causes delays for on-line implementations. In these cases, alternative causal filtering strategies 31 are more appropriate.

The quickest-detection tool developed can also be applied to different time-series, other than the sequence of growth rate of daily new positives {x n }, for instance to the daily new hospitalized individuals addressed, for the Italian case, in Fig. 6 . The sequence of daily new hospitalized individuals does not require the smoothing operation shown in the flowchart of Fig. 5 , because the collected data are less affected by gross errors. These data are also normally distributed with good accuracy, as confirmed by goodness-of-fit analysis (not reported). In Fig. 6 we see that the decision taken by using the hospitalized sequence is delayed as compared to that obtained by the sequence of daily new positives. This behaviour is intuitive because the hospitalized individuals are a subset of the positive ones, and a possible hospitalization follows the onset of clinical symptoms.

The possibility of processing different sequences of data opens the way to the design of more sophisticated decision rules based on joint processing of multiple time-series. In addition, the MAST procedure can easily accommodate different definitions of critical regime. For instance, if a pandemic growth at rate (1 + α * ) is considered acceptable, for some 0 < α * ≪ 1, the critical regime can be characterized by rates exceeding (1 + α * ), rather than 1. The necessary modifications to the MAST statistic are straightforward. We also envision that considering both the passage from a controlled to a critical regime and vice versa can be addressed by minor modifications to the MAST procedure. All these extensions are left for future work. An extended analysis of COVID-19 infection data from more countries than those covered in this paper is available on the web. 29 Finally, we hope that in the near future the publicly available data can be: (i) more reliable so as to mitigate bias effects due to, e.g., false positives, contrasting multiple test outcomes for the same individual, markedly different contagion incidence in close geographical areas, etc.; and (ii) released with finer granularity so as to allow for analyses stratified by population age, comorbidity, etc. These aspects are also relevant to the design of effective vaccination policies.

The observation model used in this paper can be formally obtained by replacing the constant growth rate (1 + α) appearing in Eq. (1) by the sequence of random variables x 1 , x 2 , . . . , yielding

where p n is the number of new cases on day n, while x n = p n+1 /p n , a time-series that makes explicit the growth rate that we seek to coopt. The model in Eq. (3) Fig. 3 ) and on the growth rate sequence of the daily new hospitalized individuals (dotted lines). On the left-side vertical axis we select a desired risk, e.g., R = 10 −4 . Then, the blue curves indicate the stopping day (about July 18 if the growth rate sequence of the daily new positive individuals is used, and August 10 if the growth rate sequence of the daily new hospitalized individuals is used) corresponding to the selected value of risk. Finally, the red curves referred to the right-side vertical axis show the mean delay ∆ corresponding to the selected risk R (about 3 days if the growth rate sequence of the daily new positive individuals is used, and below 5 days if the growth rate sequence of the daily new hospitalized individuals is used). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range.

sequence of daily new positives, for a certain region of interest; the analysis presented here relies on the data provided by the Johns Hopkins University. 27 Referring for instance to Fig. 2 , such a sequence is shown in gray in the left panel. To address gross errors, missing values and delays in reporting the data, the sequence is smoothed by a moving average filter with uniform weights (the moving average filter of length L is always assumed to have equal 1/L weights). By trial and error, we select the length L of the filter to be 21 days, and the smoothed sequence {p n } so obtained is shown in green in Fig. 2(a) . The growth rate process used as observable is then computed as {x n = p n+1 /p n }, and is represented by the green curve in Fig. 2(b) . The time-varying statistical mean {µ n } of the sequence can be estimated by low-pass filtering of the sequence {x n }, and for this we use again a moving average filter of length L = 21 days. By subtracting from each x n the estimated mean value µ n , that is, the curve in magenta in Fig. 2(b) , one obtains the sequence x n − µ n . Statistical analysis conducted by Kolmogorov-Smirnov goodness-of-fit test 32 reveals that x n − µ n , for each n = 1, 2, . . . , can be modeled by a zero-mean Gaussian random variable with (country-dependent) standard deviation σ ≪ 1. Since µ n is close to unity, we see that x k < 0 with negligible probability, hence the Gaussian approximation should not be problematic.

By observing the sequence {x n } as time n elapses, we want to detect the passage from the controlled regime H 0 , to the critical regime H 1 , when there is the exponential growth. In the controlled regime, the mean value of the random variable x n is below one, i.e., µ n = µ 0,n ≤ 1, meaning that the number of new positives remains approximately stable or decreases. Conversely, in the critical regime, the mean value is greater than unity, µ n = µ 1,n > 1, which implies an explosion of daily new positives in the long run. Formally:

critical regime H 1 ⇒ x n ∼ N (µ 1,n , σ ), µ 1,n > 1.

We assume that the x k 's are mutually independent under either regime. 1 Different nations are characterized by different values of σ , hence for each country σ is assumed known, because in practice it can be estimated on-line from the data. Conversely, the quantities {µ 0,n } and {µ 1,n } are modeled as deterministic but unknown sequences.

To detect the change of regime, we rely on the Generalized Likelihood Ratio Test (GLRT) approach, 11-13 a milestone of decision theory in scenarios where the statistical distributions of the data contain unknown parameters -in our case, the sequences of mean values {µ 0,n } and {µ 1,n }, and the time at which the passage of regime occurs. If we assume the sequences of mean values in the two regimes to be constant and known, say µ 0,n = 1 − α and µ 1,n = 1 + α (cf. equation (1)), the GLRT solution to the quickest-detection problem would be the celebrated Page's test: compare to an appropriate threshold the CUSUM statistic 14, 15, 17  

In the presence of unknown and time-varying sequences complying with the constraints on µ 0,n and µ 1,n , shown in (4)-(5), the GLRT solution to the quickest-detection problem amounts to compare the MAST statistic T n to a threshold χ:

and declaring the change of regime at the first occurrence of the threshold crossing. It is worth noting that the MAST statistic is formally obtained by replacing the unknown value of α appearing in the CUSUM statistic, with the estimate α n = |x n − 1| (constant factors can be incorporated in the threshold). The reader educated in detection theory will recognize the analogy with the energy detector arising in testing the presence of an unknown time-varying deterministic signal buried in Gaussian noise. 34 The threshold χ employed in the MAST procedure is selected to trade-off decision delay ∆ and risk R, two quantities that are defined as follows. The mean delay ∆ is the mean value of the difference between the time at which the MAST statistic crosses the threshold χ and the time of passage from the controlled to the critical regime. The risk R is defined as reciprocal of the mean time between successive false alarms (in fact, "false alarm probability" is perhaps a more familiar jargon to readers with background in detection theory), where the false alarm is defined as a threshold crossing during the controlled regime (i.e., the regime in which there is no need for intervention). In this paper, the functions R(χ) and ∆(χ) are obtained by standard Monte Carlo simulations 12, 26, 34 for relatively small values of χ, see e.g., Fig. 2(c) . We found that the first function is essentially exponential, and the second is essentially linear, consistent with known expressions for the Page's test. 14, 15 This allows us to extrapolate the behaviour of the two functions for values of χ that would be problematic to obtain from real data or by computer experiments, such as those used in Fig. 3 and similar figures for other nations. pandemic data from US. Both the second wave (left side) and the third wave (right side) are analysed. For each wave, we select a desired risk on the left-side vertical axis, e.g., R = 10 −4 . Then, the blue curve indicates the stopping day (about June 6 for the first wave and September 10 for the third wave) corresponding to the selected value of risk. Finally, the red curve referred to the right-side vertical axis shows the mean delay ∆ corresponding to the selected risk R (approximately 4 days for both the second and the third waves). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range. On the left-side vertical axis, we select a desired risk, e.g., R = 10 −4 . Then, the blue curve indicates the stopping day (about July 11, in the example) corresponding to the selected value of risk. Finally, the red curve referred to the right-side vertical axis shows the mean delay ∆ corresponding to the selected risk R (below 6 days). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range. On the left-side vertical axis, we select a desired risk, e.g., R = 10 −4 . Then, the blue curve indicates the stopping day (about July 7, in the example) corresponding to the selected value of risk. Finally, the red curve referred to the right-side vertical axis shows the mean delay ∆ corresponding to the selected risk R (below 20 days). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range. On the left-side vertical axis, we select a desired risk, e.g., R = 10 −4 . Then, the blue curve indicates the stopping day (about July 19, in the example) corresponding to the selected value of risk. Finally, the red curve referred to the right-side vertical axis shows the mean delay ∆ corresponding to the selected risk R (below 13 days). For clarity, note that the right-side scale for the delay is split into two linear ranges, for a better rendering of the small-∆ range.

Let us start, for the sake of simplicity, from the classical Susceptible-Infectious-Recovered (SIR) model of epidemic spread; a similar discussion holds for more complicated compartmental models. There are P individuals (for concreteness, this may be the population of a state or a country), grouped into three classes, or compartments: susceptible S(t), infected I(t), and removed (or recovered) R(t). Let s(t) = S(t)/P, i(t) = I(t)/P and r(t) = R(t)/P, denote the corresponding fractions of individuals. Each susceptible individual infects, on the average, β randomly-chosen other individuals per unit of time. Each infected individual will then recover (or, unfortunately, pass away) at a constant average rate γ. The popular SIR model 6, 21, 26 is formalized mathematically as follows

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 (8) is redundant, in view of the conservation constraint r(t)+ s(t)+ i(t) = 1. Dividing the first of (8) 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

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 (9), 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): 35

From the first of (9), we see that the fractions of recovered and susceptible individuals reach steady-state constant values r(∞) (typically > 0) and s(∞) (< 1) verifying r(∞) = 1 − s(∞), while i(∞) = lim t→∞ i(t) = 0. Suppose s(t) ≈ s(t ⋆ ), for some t ⋆ . 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 (8), we get

with α:=β s(t ⋆ ) − γ. Usually, β − γ > 0, while β s(∞) − γ < 0, so that i(t) initially grows exponentially at rate β − γ and eventually decreases exponentially to zero at rate β s(∞) − γ.

Epidemic data are typically collected on a daily basis, and a discrete version of (11) 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 (11) is approximated by i n+1 −i n ∆t = i n+1 − i n = α i n . This yields ∆i n :=i n+1 − i n = α i n ⇒ i n = i 0 (1 + α) n , n ≥ 1,

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 (12) 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 (8), 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 (13) is better represented by a collection of random variables, leading to the following model for the fraction of infected individuals on day n:

In (15), 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 (12) and (14) . Then, we get

yielding p n+1 = p n x n = p 1

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 (18) , 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 27 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. 6 For these reasons, we first eliminate from the data anomalous values (negative numbers) and then, in place of (18), 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 (19) 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 (20) 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 (17) .

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, 27 as described next. For concreteness, 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: 31 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 1 . 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. (20) Figure 11 . 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.

sequence {p n }, and, in turn, over the sequence of ratios {x n = p n+1 /p n }.

Finally, in Table 1 , we verify that the p-value obtained by the Kolmogorov-Smirnov (KS) goodness-of-fit test 32 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.

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. 11 . Using these periodic counterparts, we obtain the performance of the test by standard Monte Carlo computer experiments, 12, 26, 34 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.

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 µ 0,n = 1 + ε 2 cos(2πnM −1 + φ 0 ) − 1 , Figure 13 . 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. 12(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. 12 (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 [36, 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. 13 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 12, 26, 34 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 .

Coronavirus disease (COVID-19) pandemic

How will country-based mitigation measures influence the course of the COVID-19 epidemic?

Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts

Effective containment explains subexponential growth in recent confirmed COVID-19 cases in China

Inferring change points in the spread of COVID-19 reveals the effectiveness of interventions

Adaptive Bayesian learning and forecasting of epidemic evolution -Data analysis of the COVID-19 outbreak

The socio-economic implications of the coronavirus pandemic (COVID-19): A review

COVID-19 pandemic, oil prices, stock market, geopolitical risk and policy uncertainty nexus in the US economy: Fresh evidence from the wavelet-based approach

Global supply-chain effects of COVID-19 control measures

Poor, H. V. An Introduction to Signal Detection and Estimation

Testing Statistical Hypotheses

Detection of abrupt changes: theory and application

Quickest Detection

Selective review of offline change point detection methods. Signal Process

Continuous inspection schemes

Optimal sequential surveillance for finance, public health, and other areas

Robust outbreak surveillance of epidemics in sweden

Statistical surveillance of epidemics: Peak detection of influenza in sweden

A contribution to the mathematical theory of epidemics

Monitoring and prediction of an epidemic outbreak using syndromic observations

Evaluation and prediction of the COVID-19 variations at different input population and quarantine strategies, a case study in Guangdong province

Substantial undocumented infection facilitates the rapid dissemination of novel Coronavirus (SARS-CoV-2)

The effect of travel restrictions on the spread of the 2019 novel Coronavirus (COVID-19) outbreak

A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis

COVID-19 data repository

Reporting, epidemic growth, and reproduction numbers for the 2019 novel coronavirus (2019-nCoV) epidemic. Annals Intern

MAST: COVID-19 pandemic onset test -multi-country analysis and visualization

Elements of Information Theory

Signals & Systems

On the Kolmogorov-Smirnov limit theorems for empirical distributions

Quickest detection of covid-19 pandemic onset

Detection Theory

Application of the Lambert W function to the SIR epidemic model. The Coll

PW was supported by AFOSR under contract FA9500-18-1-0463.

SM developed the approach and supervised the work. PB, LMM, and DG did most of the analysis. DG prepared the figures. LMM prepared the website. PW and KP supervised the work. All authors contributed in writing the manuscript.

A visualization of MAST performance on updated COVID-19 infection data from an extended set of countries is available at https://covid-mast.github.io.

The authors declare no competing interests.