key: cord-0534870-pwdysbpg
authors: Matsui, Tomoko; Azzaoui, Nourddine; Murakami, Daisuke
title: Analysis of COVID-19 evolution based on testing closeness of sequential data
date: 2021-06-30
journal: nan
DOI: nan
sha: 7bf043e86a9039d086032ea5bf6759a8c0222aa4
doc_id: 534870
cord_uid: pwdysbpg

A practical algorithm has been developed for closeness analysis of sequential data that combines closeness testing with algorithms based on the Markov chain tester. It was applied to reported sequential data for COVID-19 to analyze the evolution of COVID-19 during a certain time period (week, month, etc.).

The COVID-19 coronavirus has spread worldwide, and as of May 31, 2021, the number of confirmed cases was 170M, and the number of deaths was 3.54M. A fourth wave of infections due to the emergence of variants with strong infectivity began hitting a number of countries in Spring 2021. Coping with a worldwide pandemic like the COVID-19 one requires understanding the infection situation. This requires development of techniques for analyzing the various types of sequential data that are available. These data include the number of confirmed infections, the number of deaths, and the number of polymerase chain reaction tests and rapid antigen tests by location and time.

As the availability of various types of data has increased in recent years, faster and more sample-efficient algorithms have been developed for statistical testing. In particular, for data collected by sensors, closeness testing of distributions to infer information from the underlying probability distributions is rapidly evolving [5, 3, 6] . Wolfer and Kontorovich, for example, developed an identity tester that determines whether sequential data represented by two Markov chains are identical [15] . Although the theory is quite rich in this area, there have been few reports of proposed algorithms being tested on actual applications or of simulation studies. Moreover, the algorithms are suitable only for discrete distributions, so a quantization technique is needed to transform continuous distributions into discrete ones. Canonne and Wimmer discussed the difficulties inherent in binning and segmentation and their limitations [4] .

We have developed a practical algorithm for closeness analysis of sequential data by combining distribution testing and algorithms based on Wolfer and Kontorovich's identity tester [15] . We tested it by using it to analyze the evolution of COVID-19 during a certain time period (week, month, etc.).

In the following section, we briefly describe related work on distribution testing and Markov chain testing. Our analysis methods are described in section 3, and their usage for analyzing spatio-temporal data like that for COVID-19 is described in section 4. We discuss the testing sensitivity in section 5 and conclude with a summary of the key points in section 6.

2 Related work 2 

Distribution testing typically involves three types of problems: the uniform testing problem, the identity testing problem, and the closeness testing problem. Let D be a distribution over a (countable) domain Ω. The uniform testing problem is to determine whether D = U Ω (the uniform distribution on Ω) or the distance between D and U Ω is far from ε ∈ (0, 1) (ε-far) [1, 7, 9] . The identity testing problem is to determine whether D = D * (a fixed distribution over Ω) or D is ε-far from D * [12, 11] . The closeness testing problem is to determine whether D and D ′ (another distribution on Ω) are equal or ε-far from each other [2, 13] . Here, we focus on closeness testing as it is useful for analyzing the the COVID-19 situation. The resulting problem is as follows.

Given sample access to distributions D and D ′ over Ω, and bounds η 1 ≥ 0, η 2 > 0, δ ∈ (0, 1), distinguish with probability at least 1 − δ between d 1 (D, D ′ ) ≤ η 1 and d 2 (D, D ′ ) ≥ η 2 whenever D, D ′ satisfy one of these two inequalities.

Here, d 1 and d 2 are the distances between two distributions. Depending on the purpose of the analysis, the total variation distance, l 2 , the χ 2 distance, or the Hellinger distance are generally used as d 1 and d 2 in distribution testing. The total variation distance is standard, and the properties of the other two distances have been theoretically and comparatively studied [6] . The χ 2 -type statistics defined by Chan et al. [5] are used here.

Learning and testing discrete distributions has been a hot research area, especially for sample complexity problems in identity testing and closeness testing [3] . Most of the work in this area has relied on independent and identically distributed (iid) sample testing, which is based on an unrealistic assumption. Emergent work has started to address the three testing problems described above, especially for data generated from a finite Markov chain (e.g., [14, 15] ). Since COVID-19 data observations are obviously not iid in time and space, we assume here that the observed proportions π (where the distribution D is estimated by π) are generated by a Markov chain over a discrete state space [s] = {s 1 , . . . , s B }; this means that it verifies the Markovian property

where p ij denotes the transition probability from state s i to state s j . Given an observed trajectory π = (π 0 , . . . , π T ) from some unknown Markov chain up to time T , we are interested in learning the transition probabilities from only this trajectory. Two strategies can be adopted for Markov chain testing: (i) naive use of distribution testing techniques (closeness testing, identity testing, and so on) for conditional transition probability comparison and (ii) less obvious comparison of the stationary distributions of the two Markov chains. With the first strategy, the discrete conditional probability distributions p i. = (p i1 , . . . , p iB ) and q i. = (q i1 , . . . , q iB ) as defined in (1) are compared for each fixed state s i . With the second strategy, this technique needs existence conditions through mixing time concept.

Wolfer and Kontorovich's identity tester [15] constructs a tester T that can determine whether a given trajectory was generated from an unknown ergodic Markov chain M having B states. They showed that the tester can determine with a probability of at least 1 − δ whether the sample trajectory was generated from M or ε-far from M .

Focusing on COVID-19, we investigated whether the pandemic evolved in the same way in different regions and for different segments of the population. We tested three analysis methods based on distribution testing and Markov chain testing that can be applied to the spatio-temporal data of COVID-19 and potentially any novel coronavirus.

In the following sections, we first formulate the problem and then describe these analysis methods.

Let us consider a population P and suppose that P = P l , where {P l } l=1,...,L is a certain segment of the population. This segment can be linked to geographic regions, socio-demographics categories, age, and other relevant auxiliary variables. We are interested in monitoring the dynamic distribution of a coronavirus like COVID-19. We are especially interested in the evolution of the distribution D l (t) of the number of infected people in population P l at time t.

Our testing framework is applicable to only discrete distributions, so we need to quantize the state space into B bins. Let us denote the discretized states as [s] = {s 1 , . . . , s B } (in the univariate case), and discretization of the interval [0, p max ], where p max is the maximum allowed proportion (in the experiments, the segmentation is uniform and p max is less than 1). To investigate the severity of COVID-19, the proportion π l t of infected people in segment P l at time t is assigned a state s i if s i < π l t ≤ s i+1 . The observed proportion isπ l t = n l t /N l , where n l t is the number of infected people in population P l at time t, and N l is the size of the population segment P l . For each t and l, the application π l t :−→ M[s] is to take a random variable in M[s], which is the set of discrete probability measures on [s].

We designed an algorithm for closeness analysis by combining distribution testing (closeness testing) and Markov chain testing in order to analyze the closeness of two sequential data. In distribution testing, there is generally assumed to be oracle access to the distributions. For closeness testing, according to Theorem 1 of Chan et al. [5] and Theorem 3.2.9 of Canonne [3] , tight upper O and lower Ω bounds for sample complexity are given by

and Ω(max(

). The algorithm we designed for closeness analysis satisfies the following two conditions under the assumption of oracle access [3, 5] . On input ε ∈ (0, 1) (a constant), C ∈ R + (an absolute constant) and B ∈ N (the number of states), it

samples from the distributions and,

• if the distributions are equal, it outputs ACCEPT with probability at least 2/3; • if the total variation distance between the distributions is greater than ε, it outputs REJECT with probability at least 2/3.

As shown in Algorithm 1, five parameters are input: ε, C, B, N ∈ N (the number of testing iterations) and µ ∈ N (the minimum number of samples for testing). The sequential data (x and y with d-dimension) are first quantized into B bins (or B states). Algorithm 1 follows the naive use strategy described in section 2.2. For each state b, the discrete conditional probability distributions (

T y (k) )) are compared. In accordance with Theorem 1 of Chan et al. [5] and Theorem 3.2.9 of Canonne [3] , m 0 is sampled from a Poisson distribution with mean m (line 19), and m 0 samples are sampled from the distributions (lines 21 and 22). For the acceptance probability, the χ 2 -type statistic z(n) defined by Chan et al. is calculated for each sample n (line 26) and compared with a threshold [3] (line 28). The statistic can be viewed as a modification of the empirical triangle distance applied to c x and c y . For the reject probability, the total variation distance d(n) is calculated for each sample n (line 27) and compared with a threshold ε.

After application of Algorithm 1, the acceptance P A and reject P R probabilities, the distance of the χ 2 -type statistic Z, and the total variation distance D for closeness testing between x and y can be calculated as the mean, median, or minimum value over all states. The minimum value is the most conservative; the mean value was used in the experiments. The χ 2 -type statistic is an estimate of χ 2 -divergence. The relation between the divergence and the total variation distance is as follows; for distributions p and q, the following inequalities hold.

Additional details and discussion can be found elsewhere ( [6] for instance). These inequalities show that the χ 2divergence d χ 2 is more conservative than the Hellinger distance d H and the total variation distance d TV . This motivated our use of the χ 2 -type statistic.

For a sequential data such as COVID-19 data, it is often demanded to analyze the evolution situation. Here, we investigate a method of periodical evolution analysis with closeness analysis. As shown in Algorithm 2, input sequence x is first segmented into L segments. Then, for each pair of segments, closeness of the pair is tested using Algorithm 1. We can analyze the periodical properties on the resulting L × L matrices for the acceptance probabilities and the distances.

Algorithm 1: Closeness analysis of sequential data. Algorithm 2: Periodical evolution analysis using testing closeness. Input: ε ∈ (0, 1), C ∈ R + , N ∈ N, B ∈ N, L ∈ N Data: x = (x 1 , x 2 , . . . , x I ) ∈ R d Output: acceptance probability P A , reject probability P R , χ 2 -type statistic Z, total variation distance D for each pair of segments 1 /* Segment x into L segments */ 2 x = {x 1 , x 2 , . . . , x L } 3 /* Test closeness for each pair of segments */ 4 for i ← 1 to L do 5 for j ← 1 to L do 6 if i = j then 7 /* Test closeness of segments x i and x j using Algorithm 1 */ 8 end 9 end 10 end

When planning measurements such as those for COVID-19, it is important to analyze the key factors, i.e., the factors that correlate with changes in, for example, the number of infections. We investigated a method for analyzing the key factors that uses a generalized additive model (GAM) [10] in which the response variable depends linearly on the unknown smooth functions of some predictor variables and the focus is on making inferences about the smooth functions. The benefit of GAM is that it takes advantage of the smoothed transforms of the predictor variables using basis functions such as smoothing splines. The distances obtained by the closeness analysis are used as the response variables. The data for the key factor candidates, e.g., vehicle and public transport increase rates, are used as predictor variables. The best model is then selected in a step-wise fashion using either Akaike Information Criterion or model residual deviance [8] .

We used reported data for the number of newly infected people n t l for each of the 53 cities on the main island of Japan as reported daily by the Tokyo metropolitan government from April 1, 2020, to May 6, 2021, along with the population N l of each city. Segmentation {P l } l=1,...,L (described in section 3.1) was linked to each city in Tokyo (which is a prefecture, not a city). The observed proportionπ l t (= n l t /N l ) was quantized into B-states, and B was set to 20. For the acceptance probabilities, the matrices between the waves tend to be darker; that is, many cities are considered to have had similar characteristics of the changes in the number of infected people for each of the months. In fact, for such cities, the number of infected people was relatively and stably small during those months.

For the distances, the overall matrix color is the darkest for January 2021, when the third wave peaked and the number of infected people was the largest. Many cities experienced an explosion of infections and different characteristics of the changes in the number of infected people for the month. Figure 2 shows the k-means clustering for the distance matrices in Figure 1 . To facilitate recognition of the differences in the level of increases in infection, the number of color codes was set to three: red indicates relatively high level, yellow indicates moderate level, and blue indicates low level. For April 2020, two cities in the heart of Tokyo, Shinjuku-ku and Minato-ku, had the highest level. This is attributed to Shinjuku-ku and Minato-ku having a popular entertainment district. Until October 2020, most cities had the lowest level. Starting with the third wave, roughly from December 2020 to February 2021, the levels of the nearby cities increased to moderate and then to high. These figures illustrate how the characteristics of the changes in the number of infected people were transformed. COVID-19 virus with higher infectivity began to gradually spread, so the characteristics of the changes in the number of infected people differed from those in previous weeks.

For the key factor analysis, we used the distances of the χ 2 -type statistic Z and the total variation distances D between all pairs of 52 weeks from 6 May 2020 to 4 May 2021 for all of Tokyo, which are included in figure 4 in which 57 weeks were used. Table 1 lists the key factor candidates used in the experiments such as vehicle and public transport increase rates and average temperature in Tokyo, which are considered to affect the rate of new infections. We set a delay of zero (no delay), one week, or two weeks between the distances.

For the distances of the χ 2 -type statistic, the R-squared (adjusted) values are listed in Table 2 . R-squared is a statistical measure of the success in explaining the response by the model, and R-squared (adjusted) is a version adjusted for the number of predictors in the model for parsimony. The table shows that the fitting was fairly accurate. The best model for a delay of two weeks was selected; it is shown in eq. (2). The s(term) indicates a smoothed transform in which term is computed using a smoothing spline, as mentioned in section 3.4. All the terms were significant: 0.001 significance level for vehicle, s(temperature), and s(deathTokyo), 0.01 for s(week), s(patientHospital), and s(roomHospital), and 0.05 for pedestrian and s(deathWorld).

Z ∼s(week) + vehicle + pedestrian + s(tepmerature) + s(deathTokyo) + s(deathWorld) + s(patientHospital) + s(roomHospital)

For the total variation distances, the fitting accuracy on the R-squared (adjusted) values was fairly good, as shown in Table 2 . The best model for a delay of two weeks was selected; it is shown in eq. (3). All the terms were significant except for s(patientHospital): 0.001 significance level for s(week), vehicle, s(temperature), s(deathTokyo), and s(infectedWorld) and 0.01 for pedestrian and s(roomHospital).

D ∼s(week) + vehicle + pedestrian + s(temperature) + s(deathTokyo)

Moreover, we divided the 52 weeks from 6 May 2020 to 4 May 2021 into two periods: (i) the 30 weeks from May to November 2020 and (ii) the 22 weeks from December 2020 to May 2021. For the first period, the R-squared (adjusted) values for both the χ 2 -type statistic and total variation distance in Table 2 were low, making it is difficult to find correlation between the distances and the key factors. For the second period, the R-squared (adjusted) values for both distances were high. As mentioned in section 4.2, the third wave roughly started in December 2020 in Tokyo, and stronger correlations between the distances and the key factors are evident for the second period.

For the distances of the χ 2 -type statistic, the best model for a delay of two weeks was selected; it is shown in eq. (4). All the terms were significant: 0.001 significance level for week, vehicle, s(deathTokyo), s(deathWorld), and s(patientHospital), 0.01 for s(transport) and infectedWorld, and 0.05 for s(temperature).

Z ∼week + vehicle + s(transport) + s(temperature) + s(deathTokyo) + infectedWorld + s(deathWorld) + s(patientHospital) For the total variation distances, the best model for a delay of two weeks was selected; it is shown in eq. (5). All the terms were significant except for s(patientHospital)): 0.001 significance level for week, vehicle, s(deathTokyo), and s(patientHospital) and 0.05 for s(transport).

These results indicate that the increase rates for vehicles and public transport can be used in the COVID-19 measurements, especially for the second period. The temperature, numbers of deaths, and number of patients in hospitals in Tokyo should be considered key factors that can be correlated with a change in COVID-19 infection rates.

We first discuss the properties of Algorithm 1 as a Markov chain tester and the sensitivity of its parameters. We do this using simulated data: (i) sequence Q x randomly generated from a transition probability matrix with 5 states (Markov chain), (ii) sequence Q y generated using sorting sequence X, and (iii) sequence Q z consisting of (100 − α)% sequences (the same as for Q x ) and an α% sequence (different from Q x ). All sequences had a length of 100 with state components s 1 = 1, . . . , s 5 = 5 (see appendix 7). Note that although sequences Q x and Q y included the same portion of each state, Q y had no Markovian property. Figure 5 shows the acceptance probabilities, the distances of the χ 2 -type statistic, and the threshold values of closeness analysis between two sequences (Q x and Q y ) with and without the Markovian property and with various values of ε and C in Algorithm 1. When ε was smaller than 0.3, the algorithm could accurately distinguish Q x and Q y for all values of C. However, when ε was 0.4 or 0.5 and C was 1 or less, the test results were incorrect although the inaccuracy was less than 4%. These results show that strict testing can be conducted with small values of ε and large values of C although with these setting, m (line 19 in Algorithm 1) becomes large and the computation cost is higher. However, the required level of strictness in closeness analysis should differ between applications, meaning that the values can be set accordingly, especially that of ε. Moreover, both C and ε should be set in accordance with the available computation power. Figure 6 shows the acceptance probabilities, the distances of the χ 2 -type statistic, and the threshold values of closeness analysis between two identical sequences (Q x and Q x ) with various values of ε and C in Algorithm 1. For ε from 0.1 to 0.9 and C from 1 to 100, the algorithm correctly determined that the two sequences were the same. Table 3 lists the acceptance probabilities, the distances of the χ 2 -type statistic, the reject probabilities, and the total variation distances of closeness analysis between (100 -α)% similar sequences (Q x and Q z ) with ε = 0.1 and C = 100 in Algorithm 1. α was varied from 0 to 5%. The algorithm was able to distinguish the similar sequences when α = 2% or more. In contrast, the classical hypothesis tests for two distributions (Wilcoxon rank-sum test and Kolmogorov-Smirnov test) could not reject the null hypothesis for all values of α. The proposed algorithm thus has strong testing power for sequential data. Table 3 : Acceptance probability, distance of χ 2 -type statistic, reject probability, and total variation distance of closeness analysis between (100 -α)% similar sequences (Q x and Q z ) with ε = 0.1 and C = 100 in Algorithm 1. 

We have designed a practical algorithm for testing the closeness of sequential data by combining distribution testing and Markov chain testing. We used it to analyze the closeness, the periodical evolution, and the key factors for the number of people infected with COVID-19 for each city in Tokyo. The results showed that whether or not the epidemic evolves in the same way in different cities or in different months or weeks with numerical indicators of the acceptance and reject probabilities and the significance levels. Examination of the properties of the algorithm as a Markov chain tester and the sensitivity of the parameters showed that strict testing can be conducted with small values of ε and large values of C under the constraint of the available computation power. Comparison with the classical Wilcoxon rank-sum test and Kolmogorov-Smirnov test demonstrated that the algorithm has a strong testing power for sequential data.

Testing random variables for independence and identity

Testing closeness of discrete distributions

A Survey on Distribution Testing: Your Data is Big. But is it Blue?

Testing data binnings

Optimal algorithms for testing closeness of discrete distributions

Which Distribution Distances are Sublinearly

On testing expansion in bounded-degree graphs

Generalized additive models

A coincidence-based test for uniformity given very sparsely sampled discrete data

Generalized Additive Models

An automatic inequality prover and instance optimal identity testing

The power of linear estimators

Testing symmetric properties of distributions

Estimating the mixing time of ergodic markov chains

Minimax testing of identity to a reference ergodic markov chain

The simulated data, Q x , Q y and Q z are as follows. 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5) Q z = (1 4 1 2 2 5 1 2 2 5 5 5 1 2 5 5 3 3 4 5 4 2 4 4 5 3 4 4 5 5 5 5 4 3 2 2 5 1 4 3 2 4 5 3 5 5 1 5 2 3 5 3 2 4 1 2 4 4 The transition probability matrix used to generate Q x is as follows.