The emergence of SARS-CoV-2, and its associated viral syndrome COVID-19, has raised important questions about the ways we analyse and identify dynamic temporal processes. In particular, by identifying similarities in principal timedependent indicators of epidemic dynamics, such as cumulative prevalence (the running total of confirmed cases over time), we can gain insight into similarities that are likely to emerge across various regions. Such similarities may be reflective of various hidden processes, be they related to the pathogen, to the response thereto or to various predisposing factors. By way of this, time series clustering has the potential to play a significant role in understanding the spatio-temporal factors governing the dynamic processes that drive an outbreak.
Vector quantisation and partitioning, or clustering, is the wider set of algorithms within unsupervised statistical learning that identify similar patterns among data in arbitrarily high-dimensional vector spaces, effectively taking a set P of N vectors in an n-dimensional vector space and assigning to each of these a label from the label set L, so that the assignment of each element of P to the groups defined by the labels comprising L minimise some objective function (typically referred to as the distance metric of the clustering). Cluster algorithms are widely used today and their practical applications are manifold, ranging from identifying clinical phenotypes in medicine and population health [1] [2] [3] [4] [5] through fraud detection [6] [7] [8] [9] [10] to image segmentation. [11] [12] [13] [14] [15] [16] Time series clustering presents a particular complication of this problem insofar as the subject of clustering is not a vector representing a single value, but rather a time series. These time series are typically not in synchrony, but rather exhibit a range of delays, lags and leads, and may depend on extrinsic and/or hidden variables. We may formulate the essential task of time series clustering as follows. Let T be a set of n time series x 1...tmax . Further, let k denote the cardinality of the label set L -in other words, the number of partitions we wish to split the data into, with k ≤ n. Then, the mapping f : T → l | l ∈ L is a clustering of the set of time series if it assigns to any element x i ∈ X | i ≤ n one (and only one) label l i ∈ L, so as to minimise an objective function (typically referred to in this context as a distance metric) J within each cluster C j = x p ∈ T | l p = j, j ∈ L defined by its label. This paper examines the use of two time series clustering algorithms -soft-DTW k-means clustering and k-shape clustering -to identify different patterns in COVID-19 prevalence in the continental United States, and comparing the results of the classifiers for inter-classifier consistency. By isolating the barycenters of the time-shifted clusters, we can identify consistent patterns in prevalence dynamics across multiple states. This in turn can be used to quantify the overall effect of pre-existing characteristics, population dynamics and nonpharmaceutical interventions (NPIs) between states.
Source data for the 48 states of the continental United States was obtained from the Starschema COVID-19 Data Set, [17] and filtered only for confirmed case counts. Data was loaded into Python 3.7 using pandas, [18] and values were scaled using tslearn.preprocessing's TimeSeriesScalerMeanVariance to µ = 0 and σ = 1. The results of this transformed raw data set are laid out, by state, in Figure 1.
Since first described by Sakoe and Chiba (1978) , [19] the dynamic time warping algorithm has been expressed in multiple formulations. The presentation below is based on Cuturi and Blondel's 2017 paper introducing Soft-DTW, with the marginal difference of using J(·, ·) instead of δ to represent the distance function. [20] Given two time series x : t x ∈ Z and y : t y ∈ Z, there exists a cost matrix ∆(x, y) for the distance function J, from which we can derive the cost matrix
For the two above-mentioned series, we may describe the set of matrices of all possible alignments as A tx,ty , which is a strict subset of {0, 1} tx×ty . Then, DTW can be defined as the function that for any pair of time series (x, y) identifies an alignment A ∈ A tx,ty so as to minimise the inner product of A with the cost matrix ∆(x, y) as
Thus, DTW can be conceived of as a search task, in which A tx,ty is the search space within which we search for an alignment A given x and y so as to minimise the inner product A, ∆(x, y) .
Soft-DTW universalises the notion underlying the DTW cost metric and the global alignment kernel metric
into a single metric. [21] Given the generalisation of the minimum metric with a smoothing factor γ ≥ 0 as
we may now define Soft-DTW as
Importantly, Soft-DTW -unlike the original DTW approach by Sakoe and Chiba [19] -is explicitly differentiable. In particular, as Saigo (2006) noted, [22] the gradient of Equation (3) can be calculated quite conveniently. LetÂ be the average alignment matrix following the Boltzmann distribution
and consequently
This can be easily calculated using backward recursion, as described in Algorithm 2 of Cuturi and Blondel (2017). [20] In addition, the notion of a clustering centroid can be generalised to the metric space comprising the time series to yield Frêchet means, also referred to in this context as barycenters. For a metric space (M, τ ), p ∈ M is a Frêchet mean of order q ≥ 1 of the time series
Thus, based on the soft-DTW metric, laid out above, we can extract from the time series of COVID-19 cumulative incidence a clustering that iteratively minimises the within-cluster sum of squares (k-means clustering). For the purposes of this paper, soft-DTW clustering was performed using tslearn 0.4.1 [23] using Python 3.7, with a γ parameter of 0.1.
k-shape clustering is a novel, robust clustering algorithm for time series that relies on iteratively refining clusters, with cross-correlation as the underlying distance metric. [24] Specifically, k-shape relies on a normalised version of crosscorrelation, referred to in this context as Shape Base Distance (SBD): time series are Z-normalised (i.e. µ = 0 and σ = 1), and the resulting cross-correlation sequence is divided by the geometric mean of the individual time series' autocorrelations. In this sense, k-shape can be understood as a k-means clustering that uses a cross-correlation based metric SBD(x, y). Let x s be the series x shifted, with zero-padding, by s, and the same be true for y s respectively, mutatis mutandis. For two time series of equal length x and y, we recursively define shift-wise cross-correlation for shifts in the range s ∈ [−m, m] as
Then, for the cross-correlation sequence, we obtain the cross-correlationρ w for any value of w ∈ 1, 2, . . . , 2m − 1 aŝ
Now, we can define the distance metric SBD(x, y) by
Because of the convolution theorem, which states that under certain conditions convolution in one domain of a time series (or more generally, any signal) is equivalent to elementwise multiplication in the other domain, [25] we can efficiently compute ψ(x, y) by taking the complex conjugate of the discrete Fourier transform of each series F(x) F(y), where is the complex conjugate operator. [24] Then, given the inverse discrete Fourier transform F −1 ,
and as Paparrizzos and Gravano showed, [24] Fast Fourier Transforms allow this to be calculated efficiently in O(|x|log(|x|)) time rather than O(|x| 2 ) time.
Similarly to the cluster analysis carried out in Subsection 2.2, k-shape clustering was performed using the tslearn package's clustering.KShape implementation, with an n init setting at 16 iterations for centroid seeds, using the result with the lowest inertia, random initialization and a convergence tolerance of 10 −6 .
After fitting the soft-DTW k-means and k-shape clustering models on the data set described in Section 2.1 with a label set cardinality (i.e. number of clusters) of 3, indicators for goodness of fit were obtained using sklearn.metrics. These show that the clustering is relatively sound. The silhouette scores (soft-DTW kmeans: 0.249, k-shape: 0.276) indicate while there is some chance of an overlap, the clustering is a relatively good fit. [26] This is confirmed by a strong Variance Ratio Criterion (Calinski-Harabasz score) of 18.809 and 20.155 for soft-DTW k-means and k-shape, respectively. [27] As Figure 2 clearly indicates, there are three distinctly characterisable patterns based on the barycenters: 1. Late peaking (k-means cluster 1, k-shape cluster 1): states in this cluster typically have a steady, consistent pattern affected only by weekly periodicities, and begin to surge around mid-June 2020.
2. Early peaking (k-means cluster 2, k-shape cluster 2): states in this cluster display a rapid-onset initial peak in April to May 2020, thereafter tapering off.
3. Bimodal (k-means cluster 3, k-shape cluster 3): within this cluster, states appear to exhibit a steady number of cases and the beginnings of a bimodal distribution over time, with a peak in April-May 2020 that subsides in June, then follows on to another rise in July and August. k-shape: 1, k-means: 3 2020-02 2020-03 2020-04 2020-05 2020-06 2020-07 2020-08 k-shape: 2, k-means: 3 2020-02 2020-03 2020-04 2020-05 2020-06 2020-07 2020-08 k-shape: 3, k-means: 3 Figure 3 : Combined time traces of k-means and k-shape classifications for the major consensus groups. Bimodal behaviour accounts for 21% of states, earlypeaking behaviour covers 17% and late-peaking, ascending behaviour accounts for over half (56%) of states. Three states do not fall within the major consensus groups.
algorithms (largely due to small differences in clustering, thus leading to different compositions for the barycenter calculation), they identify consistently the underlying pattern characteristic of the cluster. Notably, the barycenters calculated by the k-shape classification exhibit much stronger short-term (weekly) periodicity in all three clusters. At the same time, the second cluster's abnormal peak in June is much less reflected in the barycenter based on the k-shape clustering than it is on the k-means cluster, and the k-shape cluster presents a barycenter with a much flatter 'peak' in mid-April than the k-means barycenter. The distribution of time series (i.e. states) over the permutations of soft-DTW k-means and k-shape cluster assignments (see Figure 3) shows that the majority of states fall into matching soft-DTW k-means and k-shape categories, with only 3 states falling outside. Over half (56%) of states fall into the soft-DTW k-means cluster 3 and k-shape cluster 3, while 8 states (17%) fall into the soft-DTW k-means cluster 2 and k-shape cluster 2, and 9 states (21%) fall into the soft-DTW k-means cluster 1 and k-shape cluster 1. This distribution is displayed in the inter-classifier agreement matrix in Figure 4 .
In order to ascertain cross-cluster agreement, the Adjusted Rand Index (ARI) was used to quantify consensus between the k-shape and soft-DTW k-means classifiers. [28] This index, first proposed by Hubert and Arabie in 1985, is symmetric, thus it can be used to identify consensus between clusters with different metrics. At 0.864, the ARI indicates strong concurrence between the soft-DTW k-means and the k-shape classifiers. Cross-cluster agreement is illustrated in Figure 4 . As it is evident therefrom, over half of the states fall into the late-peaking (k-means cluster 3, k-shape cluster 3) category, with relatively few cases and no pronounced peaks until June 2020, after which the data evidences an oscillating but gradually increasing case count.
The strong cross-cluster agreement, covering 96% of all samples, indicates that despite their methodological differences, both the soft-DTW k-means clustering algorithm and the k-shape algorithm yield largely identical results when it comes to assigning states' time series to clusters. The strong concurrence and favourable ARI indicate that the cluster assignments are unlikely to be artefactual results of the underlying algorithms but rather reflect truly significantly distinct groupings of states by their case count time series.
k-shape and soft-DTW k-means classification strongly concur in identifying the three fundamental behavioural clusters of confirmed COVID-19 case count in the 48 states of the continental United States: a bimodal pattern, an early peaking k-shape: 1, k-means: 1 k-shape: 1, k-means: 2 k-shape: 2, k-means: 1 k-shape: 2, k-means: 2 k-shape: 3, k-means: 1 k-shape: 3, k-means: 3 pattern and a late, slower pattern that is largely stationary until approx. June 2020, then displays a rapid rise of cases.
The geographical distribution of these is worth noting. As Figure 5 shows, at the time of writing, most of the area of the continental United States follows the late peaking regime, and the calculated barycenters indicate these states are currently poised to experience further growth in case counts. Only a few states (green shades) have followed an early outbreak with a significant reduction in cases and no further resurgence, as may be considered evidence of successful mitigation/suppression efforts on their part. Finally, a number of states (blue shades) have experienced early outbreaks and are exhibiting a bimodal pattern, whereby an initial surge in April to late May 2020 has been followed not by successful suppression but a reduction followed by yet another rise in the number of reported cases of COVID-19.
As this paper has shown, time series clustering allows for finding commonalities between time series that are by necessity out of synchrony. In doing so, it can be helpful in illuminating geographical and regional patterns of disease dynamics. In particular, by using two different methods -a soft-DTW based, time-shifted k-means classifier and the correlation-based k-shape classifier -, a significant consensus between such classifications has been demonstrated where the number of confirmed COVID-19 cases in the continental United States is concerned. This lends credence to the hypothesis that epidemic dynamics of COVID-19 follow three distinct temporal patterns. These are in all likelihood conditioned by a combination of spatio-temporal factors (position along the epidemic's 'wavefront'), mitigation measures such as NPIs, their reltive effectiveness, as well as pre-existing factors of resilience and vulnerability.
Thus, by identifying the case count response, we can recognise different internally consistent clusters of case count progression over time. This may assist in understanding the governing patterns and dynamics of the SARS-CoV-2 pandemic, and assist in tailoring responses to the needs of individual areas and communities based on the temporal patterns of epidemic dynamics they exhibit.
The author declares no competing interests.
All simulations, code and data are available on Github and under the DOI 10.5281/zenodo.3970209. Shape files for the choropleth diagram in Figure 5 have been obtained from the United States Census Bureau, and are included in the data set noted above.