Modern technologies have accelerated the pace at which data is collected, leading to new frontiers in quantitative research. For example, the proliferation of wearable devices, sensors that continuously capture detailed behavioral data at an individual level, has led to large numbers of observations and data types observed for each individual over time. When the dimension of the available data and the number of observations per person increases, over-parameterized models can lead to poor efficiency and overfitting. This problem can be addressed by combining information across similar individuals within a sample of heterogenous longitudinal time series data.
Longitudinal data can be characterize by repeated observations collected for multiple individuals over a period of time. Here we consider data collected at equally spaced time series for each individual. Individuals will inevitably have varying participation, differences in the sparsity of data, i.e. varying length of time series, and different patterns and frequency of missing data. In addition, time series for different individuals often differ in their trend as well as their lagged covariance structure. We consider these differences to be elements of heterogeneity that can be modeled. Centroid based clustering would separate differences in means, and may miss out on more nuanced differences that could be observed by evaluating the time series covariance structure.
Existing methods rarely accommodates these elements of heterogeneity. Classical functional clustering methods have difficulty handling some missing data mechanisms (Chiou and Li, 2007; James and Sugar, 2003) . Popular longitudinal clustering methods are based on calculating distances measures between individuals, which may not properly capture the true underlying signal (Montero et al., 2014; Ferreira and Zhao, 2016; Genolini et al., 2015) . We propose a method to cluster individuals that also addresses these elements of heterogeneity by evaluating the autocorrelations of time series data. Correlations are ancillary with respect to mean and variance and autocorrelations can be calculated pairwise when there is missing data. Furthermore, by clustering through a latent random effect, clusters will be identified during our model training and will therefore be selected so as to optimize model fit. This approach will be able to discern differences in underlying serial variation of time series between individuals while accounting for some variability intrinsic to an individual.
We use the Wishart mixture model (WMM), defined in Hidot and Saint-Jean (2010) , to cluster individuals by their K dimensional autocorrelation matrices while simultaneously estimating group specific autocorrelations matrices. This WMM incorporates the sparsity of the time series into the Wishart distribution as the degrees of freedom parameter, allowing the model to account for varying time series sample size across individuals. As a result, low degrees of freedom corresponds with reduced influence in the model. Together, the WMM of autocorrelations can accommodate these elements of heterogeneity: differences in mean, variance, and sparsity. We then use the respective estimated group autocorrelation matrix to estimate group specific autoregression (AR) coefficients using Yule-Walker equations (Yule, 1927; Walker, 1931) . The derivation of these AR estimates have the advantage of using pooled information across individuals, resulting in a robust and parsimonious population level autoregressive mixture model (ARMM).
In Section 2, we derive the WMM model under a variety of parametric assumptions and we detail an EM algorithm to estimate model parameters in each case. Next, we estimate the ARMM using the results from the WMM and outline an AIC model selection procedure as well as detail competing methods. In Section 3, we compare the WMM with competing approaches through simulation studies and an application with COVID-19 infection rate data. Ultimately we find WMM to be a powerful approach for clustering of longitudinal time series data, particularly in the presence of nuanced AR coefficients differences.
Our data consist of I individuals, with times series vector y i , with measurements equispaced over time for individual i ∈ {1, 2, 3, ..., I}. Individual i, has a time series vector of length n * i and we denote the time point index through the superscript in
. When there is no missing data n * i is equal to the number of observations, the sample size of time series y i .
Using the Wishart distribution, we evaluate C i , a square K dimension non-singular autocorrelation matrix with a Toeplitz, bi-symmetric structure such that
where r and c indicate the elements of matrix
, the autocorrelation between data points separated by lag k where y Under situations with missing data,ρ i (k) is calculated using pairwise autocorrelations to achieve the best possible estimation given the data available.
The degrees of freedom, n i paired with C i is given by number of rows without missing data
. When there is no missing data, n i = n * i − K + 1.
The corresponding group AR model is given by
g is are the AR coefficients shared by the group, g denotes the group index, and k is the lag index. The α i and the τ 2 i are the intercept and variance of the AR model specific to individual i, respectively.
Yule-Walker estimates are consistent method of moment estimates for causal AR processes, which rely on the autocovariance or autocorrelation matrices (Brockwell et al., 1991) . In a study population there may exist only a few types of autocorrelation matrices, which can be clustered together to form a population level model for each group, thereby combining signals across individuals. The WMM is a means to cluster individuals by their autocorrelation matrices and estimate group Wishart scale matrices. Using the group scale matrices, we can perform Yule-Walker estimation for group AR coefficients.
First we define the Wishart mixture model for G total number of groups or clusters,
where g ∈ {1, 2, . . . , G}, the missing group indicators Z ig follows a multinomial distribution with π 1 , π 2 , . . . , π G as the mixing probability and C i follows a Wishart distribution with Σ g as group scale matrix and n i degrees of freedom. n i controls the peakedness of the Wishart probability density function. A small n i corresponding to a flatter distribution; thus a lower contribution to the likelihood and less influence in the WMM. The complete data likelihood for Θ = {π 1 , . . . , π g , Σ 1 , . . . , Σ g } is given as
Using the EM algorithm, we estimate π g , Σ g , and impute z ig (Dempster et al., 1977) .
In the estimation step, the function Q Θ,Θ (t) , with current estimateΘ (t) , is given as
After conditioning and noting thatẑ
Here n i controls the degree of influence an individual has on the likelihood. A small n i , often associated with a noisy estimate of C i , results in a flat f W (C i |Σ g , n i ), leading to mixed soft clustering assignments through equation (3). This allows the data from noisy individuals to be dispersed throughout the each group's estimation rather than assumed by any individual group, contrary to hard clustering techniques.
In the maximization step, maximizing Q Θ,
yields our update for the mixing probabilityπ
The score function for Σ g is
which yields our update asΣ
The estimation ofΣ g , which maps to the group AR model, is done simultaneously along with clustering, allowing all individuals to be leveraged for the estimation of each group model under the soft clustering assignments. Repeating the EM algorithm, equations, (3),
(4), and (5) until convergence leads to the estimate of Θ. The final group indicator, to be used in the ARMM, is imputed as the index that maximizesẑ ig , also known as the maximum a posteriori rule for values z ig such that z ig * = 1 where g * = argmax g {ẑ ig }.
Computing Wishart densities can be numerically unstable when n i is large as the distribution function becomes very peaked, making convergence highly sensitive to initial parameter values. In addition, there may arise scenarios where n i is unknown despite the AR model fit being known, such as in meta-analyses. We address missingness as well as the numerical instability due to large n i with a modified version of the above proposed EM algorithm where we assume a constant degrees of freedom , n g , at a group level, where η = {n 1 , n 2 , . . . , n G }. While allowing n i to match the true number of data points for the ith individual seems intuitive, recall that the WMM is only an approximation as the observations forming the autocorrelation matrix are not independent. Thus allowing degrees of freedom to be flexible parameters may allow for a better model fit. We do this by solving the score function
where ψ is the digamma function. Updates for equation (4) remains the same, but the updates for (3), and (5) now becomê
We numerically update n g solvinĝ
Again, as our design matrix D i does not follow the classical definition of a Wishart distribution which relates to covariance matrices of independent multivariate normal observatinos, as we make no such assumptions on our data the degrees of freedom n i has a different meaning. Even with our modified algorithm equation (8) which allows for n g to be estimated at the cluster level instead of at the individual level, we may still want the degrees of freedom to somewhat reflect the amount of data different individuals have. As somewhat of a balance between the two approaches we propose a third extension of the EM algorithm which adjusts degrees of freedom at a cluster level by the adding a λ g adjustment term to n i , the degrees of freedom of the individual, where Λ = {λ 1 , λ 2 , . . . , λ G }.
The update of λ g is calculated by solving score function
Updates for equations (3), and (5) now becomê
We numerically update λ g usinĝ
Upper bound U prevents numerically instability incurred due to having a large degree of freedom. For this third algorithm, repeat, in order, (9), (4), (10), and (11) until convergence. The above three versions of the proposed EM algorithm allow for different assumptions on the Wishart distribution, namely with respect to how the degrees of the freedom parameter are handled, and as such we treat the three approaches as competing methods. Our subsequent analyses were conducted with U = 50 for both variations of EM.
The ARMM representation of the WMM through the use of the Yule-Walker equations can provide an alternative representation of the model that can provide some advantages.
The mean of Wishart distribution is n i Σ g , the mean of autocorrelation matrices, C i . The cluster level covarianceΣ g is an intuitive and proportional estimator of the autocorrelation matrix for each group. We can useΣ g in the Yule-Walker estimator for AR coefficients as it relies on a matrix that is proportional to the autocovariance matrix. The matrix can be blocked asΣ
where q g is scalar, u g is a K − 1 vector, and Q g is a K − 1 square matrix. The Yule-
Walker system for estimating AR coefficients is given as the method of moments estimator
The Yule-Walker estimation results in group AR models given in equation (1) where α i and τ i are the only parameters specific to individual i. Early on, we disregarded variance by working with autocorrelations, and now it reappears in the model as τ i , the parameter for individual specific variance. TheΣ g maps to the autoregressive parameters,φ (k) g and is estimated using data across individuals, creating a robust group model. Incorporating individual i's group indicator, z ig from the WMM, we arrive at the final ARMM
The z ig govern AR model membership for individual i; φ (k) g characterizes the autoregressive behavior; α i and τ 2 i allows for across individual heterogeneity; and n i is incorporated into the WMM likelihood in equation (2) and modulates the influence of individual i. In the first phase, four example time series are mapped to autocorrelation matrices, C i , and are processed by the WMM. In the second phase,Σ g is used to estimate the AR coefficients φ g via Yule-Walker estimation. In the third phase, group indicators z ig are combined with the group AR models to create the ARMM.
We can use a pre-specified number of clusters G and lag K, based off of prior knowledge, or it can selected based on a penalization criteria such as AIC (Akaike, 1974) . The complete data likelihood of ARMM model presented in equation (12) is
where Θ AR, = π 1 , . . . , π G ; α 1 , . . . , α I ; τ 2 1 , . . . , τ 2
is the conditional normal pdf of the AR model andπ g ,Φ g , and z ig are already calculated using the WMM. At this point, we only need to estimate α i and τ 2 i asα i = 1
There are 2I parameters for α i and τ 2 i but they do not depend on group number G. The only parameters that depend on G are Φ g and π g , which exist on a G dimension simplex.
Therefore, we propose using the AIC to identify G and K, more specifically by minimizing:
We consider a few other time series clustering methods to compare performance with. As one naive alternative approach to consider, each individual is fitted with a AR model and the AR coefficients are averaged across individuals to obtain the model average. We use a Gaussian mixture model (GMM) and treat individual i's Yule-Walker AR coefficients,
i , . . . ,φ (K−1) i T as observations of a multivariate normal distribution. We elect to use Yule-Walker coefficients due to the computational burden of finding maximum likelihood estimates for each individual. The GMM complete data likelihood is given as
. We calculateπ g ,Φ g , andẑ ig using the GMM, rather than the WMM, using R package mclust (Fraley et al., 2014) .
We also explore additional clustering methods that calculate pairwise distances between all individuals for hierarchical clustering, a variety of which are outlined in Montero et al. (2014) . Alternatively, these distances can be calculated using a network community detection approach outlined in Ferreira and Zhao (2016) . When working with concave distances, Valk and Pinheiro (2012) , outlined a quasi U-statistics approach. We compared hierarchical clustering based on three relevant distance measures. These methods fall under the class of hard clustering and do not estimate the ARMM model; while WMM and GMM are soft clustering methods and do estimate the ARMM model. Therefore, for the hard clustering methods after we obtain z ig we estimate π g , Σ g , using equations (4), (5) and estimate the ARMM coefficients using Yule-Walker equations so we may compare models using AIC.
First we consider the distance defined by Galeano and Peña (2001) between two time series derived using autocorrelations
13 Next, we define d P ACF (y i , y j ), as d ACF (y i , y j ), replacing the autocorrelation with partial autocorrelation as our second distance measure. Finally, we use Piccolo distance defined as
are Yule-Walker AR coefficients, for our third distance measure (Piccolo, 1990) .
These set of methods are compared through simulation to determine relative accuracy in identifying true cluster membership.
It is well known that EM converges to the local maximum; typically multiple runs from random starting guesses are executed to determine the MLE. In practice, it's difficult to recruit more individuals in a study. Computationally, the MLE of the ARMM likelihood given in equation (13) is a difficult optimization problem, while the WMM offers a simplified estimation procedure. We have found that, under certain conditions, the WMM is a consistent estimator as n i increases i.e. as more observations per individual are obtained:
Condition 1. The EM algorithm is initialized in a strictly concave region of the global maximum.
Condition 2. The y i is a causal AR(K − 1) process for all i.
Condition 3. The number of groups, G is correctly specified.
Theorem 3.1. Let Conditions 1-3 hold, and n i → ∞ for all i, then C i converges in probability to its true autocorrelation matrix,ẑ ig converges in distribution to an indicator function for true labels and AR coefficient estimates converge in distribution to the true
Proof. The proof is left to Appendix 5.
To evaluate the performance of competing methods, we conducted two simulation studies.
We have labeled the following algorithmic denotations: EM1: (3)-(5), EM2: (6) In Figure 2 , the first two columns of plots shows that the autocorrelation and coefficients exhibits distinct separation as time series sample size increases. The third column of plots shows WMM methods outperforms the competing methods at label recovery of the true z ig , with EM2 marginally out performs the other WMM methods at label recovery. Table   1 shows that as the time series sample size increases, overall model selection of the correct G improves. The WMM methods are favored by AIC selection, with EM2 out performing all other methods. EM2, which assumes group-level degrees of freedom (as opposed to individual-level) estimates n g and performs well in settings which the time series sample size is homogenous, such as our simulation study. However, in subsequent analysis, we find that EM1 and EM3 out perform EM2 when time series sample size is variable across individuals.
Outside of WMM methods, ACF performs best at label recovery and AIC model selection. ACF is closely relate to the WMM, as both utilize autocorrelations to cluster individuals. However, the WMM is better suited for clustering autoregressive time series, weight is increased relatively on the shorter lags as the jth lag only has 2(K − j) entries in the autocorrelation matrix, whereas each lag has equal weight in the AR coefficient-based models. As a result the GMM and AR coefficient-based methods are not able to discern the true difference between individuals as distance measures appear small, while differences across the respective correlation matrices are more easily detected.
The COVID-19 time series data for 50 U.S. states, Puerto Rico (PR), and Washington D.C (DC) was obtained from the The New York Times GitHub: https://github.com/ nytimes/covid-19-data (Times, 2020) . Time series data for each state consist of cumulative number of cases, recorded daily, from January 21, 2020 to May 22, 2020. In our analysis, we applied a log transformation on the time series data. We also omitted any time points with cumulative cases less than 100, as data collection often was inconsistent during the onset of an outbreak (COVID et al., 2020 ). An ARMM was fit in order to estimate varying rates of disease transmission over time in the U.S. We used a combination of G = 1, 2, 3, 4 clusters and AR(1)-AR(3) models, using WMM and competing methods with AIC, defined in equation (14), for model selection. Table 3 . The group ranks for transmission rates and state memberships are provided in Table   3 and Figure 3 . Estimation of AR models and clustering is agnostic of gross number of cases, but models lag conditional growth over time. Since we evaluate the estimator as all n i → ∞, for ease of notation, assume n i = n ∀i and n → ∞. In addition, assume that D i is zero centered. By the asymptotic behavior of the Yule-Walker estimates for causal AR processes, autocorrelation matrices as a sequence of n, given as C i (n) = (nVar(y i )) −1 D T i D i converges in probability to their correct group autocorrelation matrix, C i (n) P → V g , g = 1, 2, 3, . . . , G, and G ≤ I. V g and C i matrices are positive definite. EvaluatingΣ ĝ
Suppose labels z ig are correct for group g,Σ g becomes: Σ g = Vg n and nΣ g = V g .
Update for the label are given as z ig = π g f W (C i (n)|Σ g , n)
G g=1 π g f W (C i (n)|Σ g , n)
= π g f W (V j |Σ g , n)
G g=1 π g f W (V j |Σ g , n)
where C i (n) P → V j for j = 1, 2, 3, . . . , G. The pdf of the Wishart distribution is given as:
where M 1 is a constant w.r.t. n.
When j = g, the individual is in the correct group, then Σ −1 g V j = nV −1 g V g = nI, When j = g, the individual is in the incorrect group, then Σ −1 g V j = nV −1 g V j = nW and w i > 0 are the eigenvalues of W. The rest is given as A = n 2 Kn/2 K i=1 w i n/2 , and
. Evaluating the upper bound for C at k = K, C < 1 Γ( n−K+1 2 ) K . Thus, the upper bound for ABC is given as
where M 2 is a constant w.r.t. n. Finally, w i exp(w i −1) ≤ 1 and the equality only holds when w i = 1, but in our case w i = 1 for every i because W = I. Thus, r = K i=1 w i exp(w i −1) , and 0 < r < 1. The upper bound for ABC is given as O n M 2 r n/2 . Thus, lim n→∞ f W (V j |Σ g , n) = 0, when j = g.
In a single causal AR process, C i andΦ i converges in probability to the true values.
Under, strict concavity, iterative convergence of the EM algorithm, our label z ig given in (15) converge in distribution to the indicator function for the correct group, as n increases.
 is a consistent AR coefficient estimator through Yule-Walker equations, Φ g = Q −1 g u g and therefore also converges in distribution.