Cluster-based dual evolution for multivariate time series: Analyzing COVID-19.

This paper proposes a cluster-based method to analyze the evolution of multivariate time series and applies this to the COVID-19 pandemic. On each day, we partition countries into clusters according to both their cases and death counts. The total number of clusters and individual countries' cluster memberships are algorithmically determined. We study the change in both quantities over time, demonstrating a close similarity in the evolution of cases and deaths. The changing number of clusters of the case counts precedes that of the death counts by 32 days. On the other hand, there is an optimal offset of 16 days with respect to the greatest consistency between cluster groupings, determined by a new method of comparing affinity matrices. With this offset in mind, we identify anomalous countries in the progression from COVID-19 cases to deaths. This analysis can aid in highlighting the most and least significant public policies in minimizing a country's COVID-19 mortality rate.

COVID-19 has resulted in a global pandemic with severe human, social, and economic costs. In order to manage the economic ramifications of prioritizing citizen safety, policymakers have sought a multi-level approach involving social distancing, business closures, and movement restrictions. For this purpose, a careful identification of the most and least successful countries at responding to the spread of COVID-19 is of great relevance. This paper meets such a demand by developing a new method to analyze multivariate time series, in which the variables are the cumulative cases and death counts of each country on each day. We have three goals: first, we analyze the cases and death counts on a country by country basis; second, we analyze the two multivariate time series in conjunction to elucidate their similarity further; and third, we determine anomalous countries relative to cases and deaths.

INTRODUCTION

Understanding the trajectories of COVID-19 cases and death counts assists governments in anticipating and responding to the impact of the pandemic. As the disease spreads, the timely identification of anomalous countries, both successful and unsuccessful, provides opportunities to determine effective response strategies. This analysis can be difficult as death counts naturally lag behind case counts.

This paper builds on the extensive literature of multivariate time series analysis, developing a new mathematical method and a more extensive analysis of COVID-19 dynamics than previously performed. Existing methods of time series analysis include parametric models, such as exponential or power-law models, and nonparametric methods, such as distance analysis, distance correlation, and network models. Both parametric and nonparametric methods have been used to model COVID-19.

Cluster analysis is another common statistical method with successful applications to COVID-19 and more broadly, epidemiology. Designed to group data points according to similarity, cluster analysis has been used to study non-communicable diseases, infectious diseases, and epidemic outbreaks such as Ebola, SARS, and COVID-19. Clustering algorithms are highly varied—common examples are K-means and spectral clustering, which partition elements into discrete sets, and hierarchical clustering, which does not specify a precise number of clusters. In this paper, we will use hierarchical clustering, K-means, and its optimal one-dimensional variant Ckmeans.1d.dp. K-means and Ckmeans.1d.dp require an initial choice of the number of clusters . We draw upon several methods to address the subtle question of how to select this . The goal of this paper is to use a dynamic and smoothed implementation of cluster analysis to study the worldwide spread of COVID-19, track the relationships between different countries’ cases and death counts, and make inferences regarding the most successful strategies in managing the progression from cases to deaths.

This paper is structured as follows: in each of the following three sections, we introduce portions of our methodology and present our results. Section II investigates the multivariate time series of cases and deaths individually. Section III analyzes the two time series in conjunction, determining suitable offsets for the number of clusters and the cluster memberships. Section IV determines anomalous countries with respect to cases and deaths. Section V summarizes the results and the new findings regarding COVID-19.

INDIVIDUAL ANALYSIS OF COVID-19 CASES AND DEATHS

Time-varying cluster analysis methodology

The most general setup of our methodology is as follows: let be a multivariate time series over an interval of length , for and , with each belonging to a common normed space . Slightly different procedures apply if is one-dimensional, namely, , or higher-dimensional.

In this paper, the two multivariate time series we present are the cumulative daily counts of cases and deaths on a country by country basis. We order the countries by alphabetical order and denote these counts by , respectively. We choose cumulative counts to best analyze the evolution of the disease over time. Our data spans 12/31/2019 to 04/30/2020, a period of days across countries.

Given the exponential nature of the data, we choose a logarithmic difference as our metric. First, we do the following data preprocessing: any entry in the data that is empty or 0—before any cases are detected—we replace with a 1, so that the log of that number is defined. Then, we define a distance on case and death counts by . Effectively, this pulls back the standard metric on under the homeomorphism and makes the positive real numbers a one-dimensional normed space.

The goal is to partition the counts into a certain number of clusters at each time . We wish to carefully choose the number of clusters in such a way that provides us meaningful inference on how the data change. A wildly varying number of clusters would obscure inference on individual countries’ cluster memberships changing with time. Thus, we combine several methods of choosing this number to reduce the bias in our estimator and perform additional exponential smoothing to yield a suitably changing number with time. In our experiments, we use six methods outlined in Appendix A. These have been chosen after experimentation and consultation with the literature, but our method is flexible and could use any combination of methods. Given cluster numbers offered by these methods, we compute the average . This is not necessarily an integer; we do not compute clusters directly with this value.

In our implementation, this average value exhibits itself as approximately locally stationary. Thus, we apply exponential smoothing to to produce a smoothed integer value . We use this value at each to obtain a clustering at that time. As the daily case and death data are one-dimensional, the most appropriate clustering method is the optimal implementation of K-means specific to one-dimensional data, Ckmeans.1d.dp. We implement this algorithm to group daily counts into clusters and sort the clusters according to the ordering on .

Similar experiments can also be performed for higher-dimensional data. Analyzing three-day rolling counts of cases and deaths requires the use of standard K-means clustering. These yield similar results to the daily analysis and can be seen in Appendix B.

Matrix analysis of multivariate time series

We record the results of this analysis in several sequences of matrices. Having performed the data preprocessing described above, first let be the matrix of (logarithmic) distances between counts at time , that is, . Next, let and be two different affinity matrices defined as follows: We term and standard and Gaussian affinity matrices, respectively. These definitions are motivated by standard constructions, but we appropriately normalize for subsequent analysis. We vary in experiments so that the matrix entries mimic Gaussian spreads over standard deviations, respectively. Then, let be an adjacency matrix defined as follows:

Finally, we define a distance on the set of dates . Let the Frobenius norm of an matrix be defined as . Given , let . Performing hierarchical clustering on these distances produces a dendrogram on the set of dates that we term the cluster evolution dendrogram. This groups moments in time according to similarity in the evolving cluster structures. In Appendix C, we include an algorithmic presentation of the steps taken in Secs. II A and II B. In Appendix D, we include a list of mathematical objects and their respective definitions used in this paper.

Results for time series of cases

In this section, we implement Ckmeans.1d.dp on daily counts of cases. Experiments using standard K-means on three-day rolling counts of cases produce similar results included in Appendix B. Our analysis supports several aspects of the empirically observed natural history regarding the spread of COVID-19 cases. The smoothed number of clusters , depicted in Fig. 1(a), ranges between . Until the end of January, there were only two clusters, with China being the only country severely impacted by the virus. However, as the virus has spread around the world, reported counts have changed day by day, with the number of clusters increasing rapidly toward a peak in early March. As depicted in Fig. 2(a), Italy was the first country to join the most severely impacted cluster, with the United States (US), Spain, France, Germany, Iran, and the United Kingdom (UK) all joining by late March. Subsequently, cluster numbers slowly declined until the end of our analysis window and appear to have stabilized. Indeed, the ranking of worst affected countries has largely stabilized in April, producing more consistent clustering results.

FIG. 1.

Smoothed number of clusters as a function of time, defined in Sec. II A. In (a), the blue and orange curves track the number of clusters for cases and deaths, respectively, from 12/31/2019 to 04/30/2020. In (b), the curves are shown after translation by the optimal series evolution offset, defined in Sec. III, computed to be . There is a strong similarity between the two curves up to this offset: both peak at 17 clusters before declining, suggesting reduced spread in the data.

FIG. 2.

Heat maps track the changing cluster membership of the 15 most severely impacted countries with respect to their counts of (a) cases and (b) deaths, respectively. Cluster membership, determined by Ckmeans.1d.dp, depicts COVID-19 severity relative to the rest of the world. Clusters are ordered with 1 being the worst impacted at any time. Darker and lighter colors correspond to smaller and greater numbered cluster labels and represent worse and less affected clusters, respectively.

In Fig. 3(a), we depict the cluster evolution dendrogram for the daily cases, defined in Sec. II B, to study the evolution of the cluster structure. This uses hierarchical clustering to determine similarity between adjacency matrices at different times, which encode the cluster structure on each day. We exclude the first 50 days, in which the cluster structure and associated adjacency matrices are all identical, with only China in its own cluster. The dendrogram identifies two distinct clusters, the larger of which contains two meaningful sub-clusters. All three (sub-)clusters identified are contiguous intervals of dates, 02/19–03/01, 03/02–03/14, and 03/15–04/30. This reveals a marked transition in cluster behavior on 03/02 for the case counts, with a smaller transition on 03/15.

FIG. 3.

Cluster evolution dendrograms, defined in Section II B for (a) cases and (b) deaths. These apply hierarchical clustering to the distance between adjacency matrices at varying times , thereby grouping different dates according to the cluster structures at these times. The -axis excludes the first 50 days for cases and 66 days for deaths, as the cluster structure of counts is trivial before these periods, respectively. Each cluster is an unbroken interval of dates. There is a clear break in the cluster structure between 03/01 and 03/02 for cases, and 03/18 and 03/19 for deaths, with a 17-day difference.

Results for time series of deaths

In this section, we implement Ckmeans.1d.dp on daily counts of deaths. The smoothed number of clusters , depicted in Fig. 1(a), ranges between . The trajectory for number of death clusters follows a similar pattern to that of cases, with a lag of approximately one month. As with the case counts, our analysis highlights the key takeaways in severely impacted countries. Although we have highlighted a one-month offset in the general evolution of COVID-19 cases and deaths, there are dissimilarities regarding the membership of the worst affected cluster. In mid-March, China moved out of the worst cluster into the second death cluster, demonstrating its relative success in responding to the pandemic. On the other hand, the US, Spain, Italy, France, and the UK have recently moved into the worst cluster, as depicted in Fig. 2(b). Examining cluster constituencies of cases and deaths over time confirms that China has managed potential COVID-19 deaths relatively effectively, while Italy, Spain, the UK, and the US have been ineffective.

In Fig. 3(b), we depict the cluster evolution dendrogram, defined in Sec. II B, for the daily deaths. We exclude the first 66 days, in which the cluster structure and associated adjacency matrices are all identical. Figures 3(a) and 3(b) show near-identical hierarchical clustering results for cases and deaths, respectively. Again, two distinct clusters are identified, with two meaningful sub-clusters within the larger cluster. All three (sub-)clusters are again contiguous intervals of dates, 03/06–03/18, 03/19–03/30, and 03/31–04/30. This reveals there is a marked transition in cluster behavior on 03/19 for the death counts, with a smaller transition on 03/31. These are 17 and 16 days later than the corresponding breaks for the case counts.

SERIES OFFSET ANALYSIS

In this section, we describe further analysis on two related multivariate time series and valued in a common normed space . With the application to COVID-19 in mind, we develop a new method that can determine if there is an appropriate time offset between the two time series. We perform several analyses for this purpose; in Sec. IV, we can subsequently study anomalous individual countries. We adopt our notation from Sec. II, using subscripts or to refer to mathematical objects pertaining to the cases or deaths counts.

First, we have already observed a clear offset in the evolution of for the time series of cases and deaths and wish to determine it precisely. We define the series evolution offset with respect to the changing number of clusters as follows: let and be the smoothed number of clusters for each time series. Given an offset , let be the translated function defined by Let the series evolution offset be the integer that minimizes the distance between functions, For our application, this offset is , confirming the one-month offset observation in Fig. 1(a).

Next, we determine the offset that minimizes the discrepancy between affinity matrices and of the two time series. Given an offset , let the normalized total offset difference between affinity matrices be defined as follows: We normalize by the number of terms in this sum, which varies with , for an appropriate comparison. When we can rewrite this as follows: Let the cluster consistency offset be the integer that minimizes the normalized total offset difference. We can also do the same for the offset with respect to the Gaussian affinity or adjacency matrices and , respectively. All these matrices are normalized, so a comparison of their values is appropriate. We choose the normalization parameter of the Gaussian affinity matrix in Eq. (2) for this purpose. We standardize notation such that always refers to the series evolution offset, while refers to the cluster consistency offset.

Results are displayed in Table I, with the optimal affinity matrix offset determined in Fig. 4. To illustrate the flexibility of the method, we choose different start dates for our offset analysis. The first 30 days carry some triviality in the cluster structure, with few cases observed outside China, so it may be desirable to exclude them from the analysis. Fortunately, the optimal offset differs only slightly with different start dates.

TABLE I.

Cluster consistency offset for various adjacency and affinity matrices at different starting dates. These are determined by minimizing the normalized total offset difference in Eq. (3), as well as its analog for Gaussian and adjacency matrices. The parameter m is defined in Eq. (2).

Optimal cases vs deaths offset
Start date	Gaussian m = 1	Gaussian m = 2	Gaussian m = 3	Adj	Aff
12/31/2019	16	16	16	20	16
01/13/2020	12	13	14	20	15
01/21/2020	12	13	14	19	15
01/31/2020	12	13	14	19	15

FIG. 4.

The normalized total offset difference as a function of the offset , defined in Eq. (3). The convex nature of this plot indicates that is a globally optimal value.

The optimal cluster consistency offset is overwhelmingly around 16. This confirms known medical findings indicating time from diagnosis to death has generally been around 17 days. Moreover, this is consistent with the results of Fig. 3, where two breaks in the cluster behavior occurred 17 and 16 days later in the death counts relative to the case counts. This is quite different from the series evolution offset of 32 days. While the cluster consistency offset seeks to align the similarity of case and death counts among individual countries, the series evolution offset seeks to quantify the overall spread of the data as a function of time.

ANOMALY ANALYSIS

Having identified a suitable offset between two multivariate time series, one can then investigate the existence of any anomalies. In this case, we use as the cluster consistency offset relative to affinity matrices, as depicted in Table I and then perform a closer analysis of the affinity matrices to identify anomalous countries. Let be the inconsistency matrix defined entry-wise by , where the absolute value of each entry is taken. Smaller entries indicate greater consistency between cases and deaths, while greater entries indicate anomalous (inconsistent) countries. Let the anomaly score of any individual country be defined as . Larger values indicate more anomalous countries and the sequence of anomaly scores can reveal the emergence and disappearance of anomalies over time. Let the lag-adjusted death rate for each country be defined as follows: These ratios may be orders of magnitude higher than standard reported death rates and are no longer bound between 0 and 1. This measure provides insight into the rate of spread and a country’s success in minimizing the number of deaths, conditional on a given number of cases days prior.

In Table II, we depict the results of ordering the ten most anomalous countries, by anomaly score, from 01/28/2020 to 04/27/2020. In Fig. 5, we display the affinity matrices for cases and deaths and the inconsistency matrix for 04/27/2020, with an offset of from Table I. We only include countries that had at least 5000 cases as of 04/30/2020. Anomalies may signify either disproportionately high or low number of deaths relative to the number of cases.

TABLE II.

The ten most anomalous countries in progression from cases to deaths as defined by their anomaly score from Sec. IV and a lag of τ = 16. AE: United Arab Emirates, AT: Austria, AU: Australia, BD: Bangladesh, BE: Belgium, BY: Belarus, CA: Canada, CN: China, DE: Germany, DO: Dominican Republic, ES: Spain, FR: France, ID: Indonesia, IE: Ireland, IL: Israel, IN: India, IR: Iran, IT: Italy, JP: Japan, KR: South Korea, ME: Mexico, MY: Malaysia, NL: Netherlands, NO: Norway, QA: Qatar, SG: Singapore, SW: Sweden, TR: Turkey, UA: Ukraine, UK: United Kingdom, US: United States, ZA: South Africa.

Ten most anomalous countries: inconsistency matrix analysis
Date	A1	A2	A3	A4	A5	A6	A7	A8	A9	A10
01/28/2020	US	UK	IT	IL	IE	IR	ID	IN	DE	FR
02/07/2020	US	DO	IT	IL	IE	IR	ID	IN	DE	FR
02/17/2020	SG	JP	KR	AU	MY	US	DE	FR	AE	CA
02/27/2020	IR	SG	MY	IT	AU	US	DE	UK	AE	CA
03/08/2020	IT	IR	SG	MY	DE	AE	CA	JP	ES	US
03/18/2020	ES	SG	IT	IR	AE	UK	NL	FR	US	KR
03/28/2020	QA	ES	TR	UK	SG	KR	AE	BY	US	IT
04/07/2020	QA	SG	KR	UK	CN	UA	NO	ZA	AU	TR
04/17/2020	BD	QA	SG	UK	AU	KR	BE	ZA	AT	FR
04/27/2020	QA	SG	BD	ME	AU	UK	SW	BE	DE	IL

FIG. 5.

(a) depicts the affinity matrix for case counts at 04/27/2020, (b) depicts the deaths affinity matrix for 04/11/2020, and (c) depicts the inconsistency matrix with an offset of from Table I. Only countries with greater than 5000 cases at 04/30 are included and ordered alphabetically along the axes. The more prominent the respective row and column in the inconsistency matrix, the more anomalous the country. The three most prominent anomalies in (c) are Qatar, Singapore, and Bangladesh.

This analysis supports several aspects of the empirically observed spread of COVID-19, identifying the most and least successful countries in the progression of cases to deaths. Early in the global spread of COVID-19, Iran and Italy were internationally known as countries that were struggling to contain the number of deaths. Table II identifies both as anomalous on 02/27/2020 and 03/08/2020, reflecting their sharp rise in deaths even before other severely impacted countries. On the other hand, Singapore is identified as anomalous during this period due to its relatively small number of deaths. As at 03/07/2020, Singapore had 130 COVID-19 cases and 0 deaths.

A similar trend continued until late March, during which Spain and Italy are identified as the most consistently anomalous countries due to their high death rates. The lag-adjusted death rates for Spain and Italy are and , respectively. Indeed, the number of deaths in Spain on 03/28/2020 was more than two times greater than the number of cases 16 days earlier. This confirms the severity of the COVID-19 pandemic: Spain and Italy suffered a large number of deaths within a short window. As of late March, Singapore was still identified as anomalous due to the relatively small number of deaths. Toward the end of our analysis window, Qatar and Australia are also identified as anomalous with low death rates, while the UK is identified as anomalous due to a high death rate. The lag-adjusted death rates for Qatar and Australia as of 04/27/2020 are and , respectively. The lag-adjusted death rate for the UK is .

CONCLUSION

In this paper, we introduce a new method of analyzing a multivariate time series via cluster analysis. Unlike typical applications of time series analysis to epidemiology, it is nonparametric; and unlike existing applications of cluster analysis to time series, we produce a dynamically smoothed number of clusters that changes over time. The analysis of case and death counts over time produces two multivariate time series, which we partition into clusters on each day. While previous studies examine fewer countries over shorter time windows, we study 208 countries over 4 months. Individual countries’ cluster membership tracks their severity of counts relative to the rest of the world, while the number of clusters reflects the overall spread of the data.

The high degree of similarity between the two time series facilitates the identification of anomalous countries in the progression of cases to deaths. We introduce another method herein, using inconsistency matrices and lag-adjusted death rates to highlight the sequential emergence and disappearance of such anomalies over time. These may be used to evaluate a country’s effectiveness at handling the pandemic, taking into account an appropriate time offset in mortality due to the disease. Our inconsistency matrices provide a multivariate method with greater generality than the included application. For this reason, they do not identify high or low mortality rates, which are only applicable in a one-dimensional context. The lag-adjusted death rate meets this purpose in our application and any other one-dimensional setting. Last, this methodology is flexible: different metrics between data, clustering methods, and means of learning offset could all be used to study related multivariate time series and identify changing similarity and anomalies.

Our analysis also provides new insights into the spread of COVID-19 across countries and over time. We show a strong similarity between the evolution of case and death counts, identifying a suitable time offset of 16 days for cluster membership between the two time series. This confirms known medical findings, indicating time from diagnosis to death as approximately 17 days. The cluster evolution dendrograms provide further support of a distinct lag between cases and deaths. These dendrograms are highly similar, also up to an offset of 16 days, and demonstrate sharp transition points at 03/02/2020 and 03/19/2020 for cases and deaths, respectively, again with a 17-day difference. These transitions reflect the natural history of the spread of COVID-19 cases and deaths, respectively. On 03/02/2020, numerous countries began to report their first instances of COVID-19 cases, predominantly imported from Iran and Italy. On 03/19/2020, Italy’s death toll surpassed that of China. Less pronounced transitions exist on 03/15/2020 and 03/31/2020 for cases and deaths, respectively. Again, a 16-day offset is observed.

On the other hand, the time offset between the evolution of the number of clusters is 32 days. One explanation for the series evolution offset being longer is that there is an additional delay between cluster membership changes with respect to cases and deaths that can be attributed to stresses on a country’s healthcare resources. First, the number of cases may increase significantly, placing a country into a different cluster relative to cases and overwhelming its healthcare resources, thereby leading to a greater number of death counts. That is, the progression from elevation in case clusters to death clusters is not necessarily due to a natural progression from infection to death, but involves mediating factors like stresses on hospital capacity. Perhaps the initial wave of patients can be treated with ventilators, but these may quickly run out, causing more deaths from later instances of cases. Regardless, it is an interesting observation that the offset of 32 days in the number of clusters does not minimize the offset in affinity or adjacency matrix norm differences.

This analysis may assist in identifying the characteristics of the most and least successful government strategies for managing COVID-19. In particular, Singapore, Qatar, Australia, and South Korea are four countries whose policies have been most successful in minimizing COVID-19 mortality. Each of these countries provided a substantial amount of easily accessible testing in the early stages of COVID-19 development. Singapore and Australia also closed their borders to travel before a critical mass in total case counts was established and were early to implement strict lockdown procedures.

By contrast, Italy, Spain and the UK are three countries whose policies managed the progression from COVID-19 cases to deaths least effectively. Many argue that lockdown procedures in Italy and Spain, although severe once in place, were implemented too late. Similarly, the UK initially elected not to shut down large gatherings or introduce social distancing measures in an attempt to build herd immunity among the community. Ultimately, however, the UK did implement strict lockdown policies as mortality rates rose.

These findings suggest that the timeliness of various lockdown procedures is perhaps more important than their severity. Countries with easy access to early testing also appear to manage the progression from cases to deaths more effectively. Conversely, countries that struggled to minimize their COVID-19 mortality rate also exhibit some general similarities. First, these countries were slow to implement measures that would restrict people’s movements. Second, many of these countries carried an early high case burden, suggesting that mediating factors such as undue stress from finite healthcare resources may contribute to the mortality rate.

Overall, this paper introduces a new method for analyzing multivariate time series individually and in conjunction, thereby providing new insights into the caseload and mortality rate affecting different countries. As the pandemic evolves, it is the objective of emerging research to facilitate timely and appropriate means of producing effective government strategies for minimizing the extensive human, social, and cultural costs of COVID-19.

Algorithm

Cluster-based evolution analysis

TABLE III.

Mathematical objects and definitions.

Mathematical objects glossary
Object	Description
D(t)	Distance matrix between log counts
Aff(t)	Standard affinity matrix
G(t)	Gaussian affinity matrix
kav(t)	Unsmoothed number of clusters obtained as average of six methods
k^^(t)	Smoothed number of clusters
Adj(t)	Adjacency matrix coding cluster outputs for k^^(t) clusters
d(s, t)	Frobenius distance between adjacency matrix of various dates
δ	Series evolution offset with respect to number of clusters
τ	Cluster consistency offset with respect to cluster membership
Inc(t)	Lag-adjusted inconsistency matrix
aj(t)	Anomaly score of country j
LDRj(t)	Lag-adjusted death rate of country j
∥f_δ−g∥_L^1	L1 norm between functions