Analysis of COVID-19 evolution based on testing closeness of sequential data

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.).


Introduction
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 (Chan et al. 2014;Canonne 2020;Daskalakis et al. 2018). Wolfer and Kontorovich, for example, developed an identity tester that determines whether sequential data represented by two Markov chains are identical (Wolfer and Kontorovich 2020). 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 (Canonne and Wimmer 2020). The main criticism of these algorithms is that generally the domain of the distribution is taken as [n] = {1, . . . , n}, which is not always realistic or representative of the true data. To overcome this limitation, a suitable quantization is needed, as suggested by Canonne and Wimmer (2020). In this work, since we did not have a prior information on the data distribution, we adopted a uniform discretization, for a number of bins which was determined empirically via numeric search.
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 (Wolfer and Kontorovich 2020). We tested it using it to analyze the evolution of COVID-19 during a certain time period (week, month, etc.). Although Markov switching models and Markov agent models have been widely used for general compartmental models in epidemiology such as for SIR (susceptible-infectiousrecovered) and SEIR (susceptible-exposed-infectious-recovered) models to represent the state transition (Bestehorn et al. yyy;Boukanjime et al. 2020;Gribaudo et al. 2021;Larsen et al. 2020;Raherinirina et al. 2021), there has been little use of Markov chains to model COVID-19 data. Ma et al. (2021) recently proposed using a Markov process combined with LSTM (long short-term memory) model to categorize the reported COVID-19 cases. To our knowledge, there have been no reports of applying Markov chains to COVID-19 data using testing techniques.
In the following section, we briefly describe related work on distribution testing and Markov chain testing. Our analysis methods are described in Sect. 3, and their usage for analyzing spatio-temporal data like that for COVID-19 is described in Sect. 4. We discuss the testing sensitivity in Sect. 5 and conclude with a summary of the key points in Sect. 6.

Distribution testing
Distribution testing is a field of computer science concerned with statistical (composite) hypothesis testing questions, with a focus on finite sample guarantees and efficient algorithms. While the findings of many distribution tests have been reported, the main focus has been on three 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) (Batu et al. 2001;Goldreich and Ron 2011;Paninski 2008). The identity testing problem is to determine whether D = D * (a fixed distribution over Ω) or D is ε-far from D * (Valiant and Valiant 2011;Acharya et al. 2015;Valiant and Valiant 2017). The closeness testing problem is to determine whether D and D (another distribution on Ω) are equal or ε-far from each other (Batu et al. 2013;Diakonikolas and Kane 2016;Valiant 2011). Here, we focus on tolerant closeness testing as it is useful for analyzing the COVID-19 situation. The resulting tolerant closeness tester is as follows.
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 (Daskalakis et al. 2018). The χ 2 -type statistics defined by Chan et al. (2014) are used here.

Markov chain testing
Learning and testing discrete distributions has been a hot research area, especially for sample complexity problems in identity testing and closeness testing (Canonne 2020). 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., Wolfer andKontorovich 2019, 2020). 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 i j 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 testing 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 i B ) and q i. = (q i1 , . . . , q i B ) 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 (Wolfer and Kontorovich 2020) constructs a tester T that can determine whether a given trajectory was generated from an unknown ergodic Markov chain M having B states. The following distance between Markov chains M 1 and M 2 is used.
where . T V stands for the total variation norm (see Wolfer and Kontorovich 2020). 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. This issue has also been studied by Dikkala and Gravin (2018), who, inspired by the early work of Kazakos (1978), proposed a difference measure that captures the scaling behavior of the total variation distance between growing trajectories of the Markov chains. They then presented efficient identity testers and gave its information lower bounds. Recently, (Cherapanamjeri and Bartlett 2019) succeeded to remove a dependency in the hitting time of the sample complexity for symmetric chains. Fried and Wolfer (2021)

Analysis methods using distribution testing and Markov chain testing
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.

Closeness analysis 2. Periodical evolution analysis 3. Key factor analysis
In the following sections, we first formulate the problem and then describe these analysis methods.

Observation model formulation
Let us consider a population P and suppose that P = P , where {P } =1,...,L is a partition of the population and P 's are disjoints. This segmentation 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 (t) of the number of infected people in segment P 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 π t of infected people in segment P at time t is assigned a state s i if s i < π t ≤ s i+1 . The observed proportion isπ t = n t /N , where n t is the number of infected people in population P at time t, and N is the size of the population segment P . For each t and , the applicationπ t :−→ M[s] is to take a random variable in M[s], which is the set of discrete probability measures on [s].

Closeness analysis
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. (2014) and Theorem 5.9 of Canonne (2020), tight upper O and lower Ω bounds for sample complexity with the total variation distance in Eq. (2) are given by The algorithm we designed for closeness analysis satisfies the following two conditions under the assumption of oracle access (Canonne 2020;Chan et al. 2014). On input ε ∈ (0, 1) (a constant), C ∈ R + (an absolute constant) and B ∈ N (the number of states), it takes C · max( B 2/3 ε 4/3 , B 1/2 ε 2 ) 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 Sect. 2.2. For each state b, the discrete conditional probability distributions , . . . , , . . . , )) are compared. In accordance with Theorem 1 of Chan et al. (2014) and Theorem 5.9 of Canonne (2020), m 0 is sampled from a Poisson distribution with mean m (line 21), and m 0 samples are sampled from the distributions (lines 23 and 24). For the acceptance probability, the χ 2 -type statistic z(n) defined by Chan et al. is calculated for each sample n (line 28) and compared with a threshold (Canonne 2020) (line 30). 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 29) 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 ( Daskalakis et al. 2018 for instance). These inequalities show that the χ 2 -divergence 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.
Note that the distance also depends on the mixing properties of the Markov chains and the stationary distribution, particularly when the number of states is small (Wolfer and Kontorovich 2020). For such a case, the mixing time should be estimated, for example, according to Algorithm 1 in Wolfer (2020) and confirmed to be smaller than m (line 21) in Algorithm 1.

Periodical evolution analysis
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 using Markov chain modeling.
Input: ε ∈ (0, 1), C ∈ R + , N ∈ N, B ∈ N, μ ∈ N Data: x = (x 1 , x 2 , . . . , x I ), y = (y 1 , y 2 , . . . , y J ) ∈ R d Output: acceptance probability P A , reject probability P R , χ 2 -type statistic Z , total variation distance D for each state 1 /* Quantize x and y into B d bins (or B d states of Markov chains) Sample variables m 0 from Poisson distribution with mean m 23 Sample a set S x of m 0 samples from Markov chain with transition probability Sample a set S y of m 0 samples from Markov chain with transition probability Algorithm 2: Periodical evolution analysis using testing closeness.

Key factor analysis
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) (T.J. Hastie 1990) 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 (Hastie 1992).

COVID-19 sequential data
We used reported data for the number of newly infected people n t 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 of each city. Segmentation {P } =1,...,L (described in Sect. 3.1) was linked to each city in Tokyo (which is a prefecture, not a city). The observed proportionπ t (= n t /N ) was quantized into B-states, and B was set to 20. Figure 1 shows 53 cities ×53 cities matrices of acceptance probabilities (the mean of P A (b) over all states in Algorithm 1) and distances of χ 2 -type statistics (the mean of Z (b) over all states in Algorithm 1) between all pairs of 53 cities in Tokyo for each month from April 2020 to April 2021, calculated using Algorithm 1. C and μ were 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 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 Fig. 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. Figure 3 shows the matrices of acceptance probabilities, distances of χ 2 -type statistics, reject probabilities (mean of P R (b) over all states in Algorithm 1), and total variation distances (mean of D(b) over all states in Algorithm 1) between all pairs of 13 months for Shinjuku and Tachikawa calculated using Algorithm 2. C and μ were chosen empirically and set to 100 and 3, respectively. Tachikawa-shi is located in the middle west of Tokyo, in a suburban area. For Shinjuku-ku (in the heart of Tokyo), as in Fig.  1, almost all the pairs are different while the May-October 2020 pair are similar. For Tachikawa-shi, the pairs from April to November 2020 and for February and March 2021 are similar. The number of infected people for these months was relatively and stably small. This figure illustrates the characteristics of monthly COVID-19 evolution for both cities.  For all cities in Tokyo, 57 weeks × 57 weeks matrices of acceptance probabilities, distance of χ 2 -type statistic, reject probability, and total variation distance Figure 4 shows the matrices of acceptance probabilities, distances of χ 2 -type statistic, reject probabilities, and total variation distances between all pairs of 57 weeks from 1 April 2020 to 5 May 2021 for all of Tokyo calculated using Algorithm 2 and all the numbers accumulated for all the cities in Tokyo. C and μ were chosen empirically and set to 100 and 3, respectively. The acceptance probabilities show that the weeks from April to June, 2020 and for August and September, 2020, tended to be similar among the cities. The distances show that the weeks in January, April, and May 2021 were very different. This indicates that the number of infected people for the weeks in January 2021 dynamically changed, probably because of an increase in contacts between people due to year-end and beginning-of-year parties and meetings. In April and May 2021, variants of the 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.

Key factor analysis for COVID-19 evolution
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 Fig. 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. (3). The s(term) indicates a smoothed transform in which term is computed using a smoothing spline, as mentioned in Sect. 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) 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. (4). 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 Sect. 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. (5). All the terms were significant: 0.001 significance level for week, vehicle, s(deathTokyo), s(deathWorld), and s(patientHospital), For the total variation distances, the best model for a delay of two weeks was selected; it is shown in Eq. (6). 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.

Discussion
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 A). 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 21 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

Conclusions
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.