Francesco Amato
email: francesco.amato@univ-lyon2.fr

Julien Jacques
email: julien.jacques@univ-lyon2.fr

Clustering Longitudinal Ordinal Data via Finite Mixture of Matrix-Variate Latent Gaussians

Keywords: clustering, Données ordinal longitudinales, Tenseurs, Modèle de mélange, Loi gaussienne matricielle Keywords, Clustering, Ordinal longitudinal data, Three-way data, Mixture models, Matrix-variate Gaussians

Dans les sciences sociales ou la médecine, les études sont souvent basées sur des questionnaires demandant aux participants d'exprimer des réponses ordonnées à plusieurs reprises au cours d'une période d'étude. Nous présentons un modèle pour effectuer un clustering temporel sur ces données. Le modèle repose sur un mélange de distributions normales à variation matricielle, en tenant compte simultanément des structures de dépendance temporelle internes et intermédiaires. Un algorithme MC-EM pour l'estimation du modèle est utilisé. Des applications sur les données synthétiques et réelles sont présentées.

Mots-clés.

Context

In many areas of humanities and social sciences, the studies are based on questionnaires completed by participants several times over the study period. The researchers then analyse these questionnaires to determine typical behaviours within the studied population, being especially interested in their time evolution. Nonetheless, modelling temporal evolution is far from trivial. The most basic approach consists in performing analyses independently at each temporal phase, and then trying a posteriori to find links between these different analyses, by seeking from one phase to the other to find similar or different typical behaviours. An example is [START_REF] Selosse | Analysing a quality of life survey using a co-clustering model for ordinal data and some dynamic implications[END_REF], clustering of ordinal data for an application in psychology. The ideal way to cluster these data would be to account for the 1 temporal evolution, modelling all the responses to the questionnaires at the same time. We propose a model-based clustering technique aiming at facilitate such temporal analysis, by grouping together the units behaving similarly in time. Over the decades, research has produced a vast number of different approaches to clustering. From our prospective, probabilistic (or model-based) clustering offers the advantage of clearly stating the assumptions behind the clustering algorithm, and allows cluster analysis to benefit from the inferential framework of statistics to address some of the practical questions arising when performing clustering: determine the number of clusters, detecting and treating outliers, assessing uncertainty in the clustering [START_REF] Bouveyron | Model-based Clustering and Classification for Data Science: with applications in R[END_REF].

Related works

An approach to clustering longitudinal data consists in arranging the data in a threeway format and modelling them through a matrix-variate mixture model. This approach offers the advantage of accounting for the overall time-behavior, grouping together the units that have a similar pattern across and within time. While not being new [START_REF] Kaye | The mixture method of clustering applied to three-way data[END_REF], matrix-variate distributions have recently gained attention, and mixtures of matrix-normals (MMN) have been developed and applied both in a frequentist framework in Viroli, 2011a and within a Bayesian one by Viroli, 2011b, where it was used to cluster Italian provinces based on a longitudinal crime-related score. From a frequentist point of view, these models represent a natural extension of the multivariate normal mixtures to account for temporal (or even spatial) dependencies, and have the advantage of being also relatively easy to estimate by means of EM algorithm (a nice short description of the EM application to MNN is provided in §2.1 of Wang et al., 2020). More recently, in Gallaugher et al., 2018 and Melnykov et al., 2018, 2019 extensions for non-normal skewed cases have been proposed and applied. However, matrix-variate models suffer from over-parametrization that leads to estimation issues. To overcome this issue a more parsimonious model [START_REF] Sarkar | On parsimonious models for modeling matrix data[END_REF] and a new R package [START_REF] Zhu | MatTransMix: an R Package for Matrix Model-Based Clustering and Parsimonious Mixture Modeling[END_REF] has been proposed. Despite their efficacy, up to now these methods have only been applied to continous data. Our model expands the use to matrix-variate mixtures to ordinal data.

Model

Let denote by y i,j,t , i = 1, ..., N ; j = 1, .., J and t = 1, ..., T the observation of the j-th variable for the i-th unit at time t, that is: imagine to observe N units and measuring J different ordinal variables T times throughout the course of the study. Let us reorganize our data in a random-matrix form such that Y = {Y i } N i=1 is a sample of J × T -variate matrix observations (i.e. Y i = (y i,j,t ) ∈ R J×T ). Then, we can assume that each variable y i,j,t is the manifestation of an underlying latent variable z i,j,t which follows a Gaussian distribution, as done in the clustMD model [START_REF] Mcparland | Model based clustering for mixed data: clustMD[END_REF]. At this point, we can assume that each observed matrix Y i is indeed the manifestation of a latent random matrix Z i , which follows a matrix-normal distribution.

Y i =        y i,1,1 • • • y i,1,t • • • y i,1,T . . . . . . . . . • • • . . . y i,j,1 • • • y i,j,t • • • y i,j,T . . . • • • . . . . . . . . . y i,J,1 • • • y i,J,t • • • y i,J,T        ←-Z i =        z i,1,1 • • • z i,1,t • • • z i,1,T . . . . . . . . . • • • . . . z i,j,1 • • • z i,j,t • • • z i,j,T . . . • • • . . . . . . . . . z i,J,1 • • • z i,J,t • • • z i,J,T       
To map from Y i to Z i , let γ j denote a C j+1 -dimensional vector of thresholds that partition the real line for the j-th ordinal variable that has C j levels and let the threshold parameters be constrained such that -∞ = γ j,0 ≤ γ j,1 ≤ ... ≤ γ j,C j = ∞. If the latent z i,j,t is such that γ j,c1 < z i,j,t < γ j,c then the observed ordinal response, y i,j,t = c.

So, by assuming that each Z i follows a matrix-normal distribution, we can then cluster our data by means of finite Mixture of Matrix-Normals (MMN) (Viroli, 2011a). Indeed, let Z i ∼ MN (J×T ) (M, Φ, Σ), where M ∈ R J×T is the matrix of means, Φ ∈ R T ×T is a covariance matrix containing the variances and covariances between the T occasions or times and Σ ∈ R J×J is the covariance matrix containing the variance and covariances of the J variables. The matrix-normal probability density function (pdf) is given by

f (Z|M, Φ, Σ) = (2π) -T J 2 |Φ| -J 2 |Σ| -T 2 exp - 1 2 tr[Σ -1 (Z -M )Φ -1 (Z -M ) ⊺ ] . (1)
The matrix-normal distribution represents a natural extension of the multivariate normal distribution, since if Z ∼ MN (J×T ) (M, Φ, Σ), then vec(Z) ∼ MVN JT (vec(M ), Φ ⊗ Σ), where vec(.) is the vectorization operator and ⊗ denotes the Kronecker product. Then, the mean and the variance of the matrix-normal distribution are:

E(vec(Z)|M, Φ, Σ) = vec(M ) and V(vec(Z)|M, Φ, Σ) = Φ ⊗ Σ. ( 2
)
Being a special case of the multivariate normal distribution, the matrix-normal distribution shares the same various properties, like, for instance, closure under marginalization, conditioning and linear transformations [START_REF] Gupta | Matrix Variate Distributions[END_REF]. The separability condition of the covariance matrix has the twofold advantage of allowing the modeling of the temporal pattern of interest directly on the covariance matrix Φ and of representing a more parsimonious solution than that of the unrestricted Φ ⊗ Σ.

The pdf of the MMN model is given by

f (Z|π, Θ) = K k=1 π k ϕ (J×T ) (Z|M k , Φ k , Σ k ),
where ϕ (J×T ) represents the density function of a J ×T -dimensional matrix-variate normal, K is the number of mixture components, π = {π k } K k=1 is the vector of mixing proportions, subject to constraint K k=1 π k = 1 and

Θ = {Θ k } K k=1 is the set of component-specific parameters with Θ k = {M k , Φ k , Σ k }.
In addition to Z, we introduce a latent binary variable that indicate whether the unit i belongs to the k-th cluster, ℓ i = (ℓ i1 , ..., ℓ iK ), such that ℓ ik = 1 if the i-th unit belongs to the k-th cluster. Let the binary vector

Y R i = (Y R i1 , ..., Y R iR ) of length R
indicate the observed response patter for the i-th unit, such that if the r-th pattern is observed then Y R ir = 1 and any other entry in the vector equals zero. We can derive the joint density of

Z i , Y R i , ℓ i as: f (Y R i , Z i , ℓ i ) = f (Y R i |Z i , ℓ i )f (Z i |ℓ i )f (ℓ i ).
Assuming that:

ℓ i ∼ M(1, π), π = (π 1 , ..., π K ) Z i |ℓ ik = 1 ∼ MN (J×T ) (Z i |Θ k ), Θ k = {M k , Φ k , Σ k }, Y R i |Z i , ℓ ik = 1 ∼ M(1, ξ R i ), ξ R i = (1 Ω 1 (Z i ), . . . , 1 Ω R (Z i ))
we get:

f (ℓ i ) = K k=1 π ℓ ik k , f (Z i |ℓ i ) = K k=1 ϕ (J×T ) (Z i |Θ k ) ℓ ik , f (Y R i |Z i , ℓ i ) = R r=1 1 Ωr (Z i ) Y R ir ,
where M indicate the multinomial distribution and 1 Ωr (Z i ) is the indicator function that equals 1 when the elements in Z i have values that determine the r-th pattern: one can imagine Ω r as a J × T matrix whose elements are the indicator functions of the thresholds linked to the r-th pattern. Of course, when Z i is given, the the value of Y i is no more random.

Estimation

A key point is of course the choice of the thresholds γ = {γ j } J j=1 . In Corneli et al., 2020, thresholds are fixed in advance to avoid identifiability and computational complexity issues. While this can be a starting point in the estimation process, our model aims at treating them as parameters and estimating them, as in [START_REF] Mcparland | Clustering ordinal data via latent variable models[END_REF][START_REF] Mcparland | Model based clustering for mixed data: clustMD[END_REF] If the thresholds are assumed to be fixed in advance, the estimation process clearly simplifies. To estimate the model, since we do not observe neither Z nor ℓ, we resort to the EM algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. The complete log-likelihood can be then written as

log L(π, Θ|Y R , Z, ℓ) ∝ N i=1 R r=1 Y ir 1 Ωr (Z i )+ K k=1 ℓ ik log(π k )- T J 2 log(2π)- J 2 log(|Φ k |)- T 2 log(|Σ k |) - 1 2 tr[Σ -1 k (Z i -M k )Φ -1 k (Z i -M k ) ⊺ ]
The expectation step (E-step) of the EM algorithm consists of computing the expectation of the complete log-likelihood with respect to the latent data Z and the cluster labels ℓ.

In the first step we attempt to estimate the latent variables. For each response pattern r, we can approximate the value of Z i |ℓ i as the expected value of the truncated multivariate normal (using (2)), given the parameters Θ k of the assigned cluster. The latent variable ℓ can be computed by means of Bayes'theorem as: Mixture of matrix-variate normal distributions can be an efficient way to cluster longitudinal continuous data. Assuming that ordinal variables is a discretization of a latent continuous variable allows us to extend the use of these MMN to ordinal variables. Numerical study on synthetic data sets as well as real data application concerning diet choice during the pandemic (François-Lecompte et al., 2020) will be presented.

E(ℓ ik |Y R i = r, Θ, π) = π k Ωr f (Z|Θ k )dZ

K k=1 π k

 k Ωr f (Z|Θ k )dZ , which requires Monte-Carlo approximation on the multivariate reparametrization. By taking the first derivatives of the log-likelihood, the M-step goes by: ik(Z i -Mk ) ⊺ Σ-1 k (Z i -Mk ) J N i=1 ℓ ik , Σk = N i=1 ℓ ik (Z i -Mk ) Φ-1 k (Z i -Mk ) ⊺ T N i=1 ℓ ik 4 Conclusions