Matthieu Marbac

Christophe Biernacki

Vincent Vandewalle

Latent class model with conditional dependency per modes to cluster categorical data

Keywords: categorical data, clustering, integrated complete-data likelihood, MCMC algorithm, mixture models, model selection Mathematics Subject Classification (2000) 62H30, 62F15, 62-07, 62F07

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Clustering [START_REF] Jajuga | Classification, clustering and data analysis: recent advances and applications[END_REF] is a worthwhile method for analysing complex data sets. It aims to extract the main information from the data by grouping individuals into homogeneous classes. This paper focuses on the cluster analysis of categorical variables which occurs in many different fields (biology, sociology, marketing...). Clustering methods can be split into two approaches: the geometric ones based on distances between individuals and the probabilistic ones which model the data generation.

Geometric approaches are generally simpler and faster than the probabilistic ones, but they are sensitive to the choice of the distance between individuals. For the categorical data analysis, either they define a metric in the initial variable space like the k-means [START_REF] Huang | Automated variable weighting in k-means type clustering[END_REF], or they compute a metric on the factorial axes [START_REF] Chavent | A partitioning method for the clustering of categorical variables[END_REF]. Note that lots of geometric approaches can be interpreted as probabilistic ones [START_REF] Govaert | Data analysis[END_REF] revealing probabilistic hidden assumptions. Moreover, the probabilistic approaches allow to solve difficult questions, like the class number selection, in a rigorous mathematical framework.

The probabilistic approaches define a homogeneous class as the subset of the individuals arisen from the same distribution. In this context, the most classical approach models the data distribution by a finite mixture model of parametric distributions [START_REF] Mclachlan | Finite mixutre models[END_REF]. The resulting partition is meaningful since each class is described by the parameters of one mixture component.

The latent class model [START_REF] Goodman | Exploratory latent structure analysis using both identifiable and unidentifiable models[END_REF] which assumes the conditional independence between variables (further referred as CIM for conditional independent model) is a popular probabilistic approach to cluster categorical data. Its interpretation is easy since classes are explicitly described by the probability of each modality for each variable. Moreover, the sparsity involved by the conditional independence assumption is a great advantage since it circumvents the curse of dimensionality. In practice, this model obtains good results in lots of applications [START_REF] Hand | Idiot's Bayes Not So Stupid after All?[END_REF]. However, it leads to severe biases when its main assumption is violated (see Section 5), like an overestimation of the number of components [START_REF] Van Hattum | Market Segmentation Using Brand Strategy Research: Bayesian Inference with Respect to Mixtures of Log-Linear Models[END_REF]. Furthermore, the larger the number of variables, the higher the risk to observe conditionally correlated variables in a data set, and consequently the higher the risk to involve such biases by using CIM.

Different models relax the conditional independence assumption. Among them, the multilevel latent class model [START_REF] Vermunt | Multilevel latent class models[END_REF] assumes that conditional dependency between the observed variables can be explained by other unobserved variables. Another approach considers the intra-class dependencies by using a single latent continuous variable and a probit function [START_REF] Qu | Random Effects Models in Latent Class Analysis for Evaluating Accuracy of Diagnostic Tests[END_REF]. Recently proposed, the mixture of latent trait analyzers [START_REF] Gollini | Mixture of latent trait analyzers for modelbased clustering of categorical data[END_REF][START_REF] Bartholomew | Latent variable models and factor analysis: A unified approach[END_REF]) is a good challenger for CIM. It assumes that the distribution of the observed variables depends on many latent variables: one categorical variable (the class) and many continuous latent variables (modeling the intra-class dependencies between the observed categorical variables). However, its parameter inference is not so easy, so the authors advise to use a variational approach but they also propose to use the Gauss-Hermite quadrature to assess the log-likelihood. Although this model is very flexible, intra-class dependency is hardly interpretable since the intraclass correlations are interpreted throughout relationships with unobserved continuous variables.

The log-linear models' [START_REF] Agresti | Categorical data analysis[END_REF] purpose is to model the individual log-probability by selecting interactions between variables. Thus, the most general mixture model is the log-linear mixture model where all the kinds of interactions can be considered. It has been used for a long time [START_REF] Hagenaars | Latent structure models with direct effects between indicators local dependence models[END_REF] and it obtains good results in many applications [START_REF] Espeland | Using Latent Class Models to Characterize and Assess Relative Error in Discrete Measurements[END_REF][START_REF] Van Hattum | Market Segmentation Using Brand Strategy Research: Bayesian Inference with Respect to Mixtures of Log-Linear Models[END_REF]. However this model family is huge and the model selection stays a real challenge. In the literature, authors either fix the considered interactions in advance or they they perform a deterministic search like the forward method which is sub-optimal. Furthermore, the number of parameters increases with the conditional modality crossings, so there is an over-fitting risk and interpretation becomes harder. This paper presents the conditional modes model (referred by CMM) which groups the variables into conditionally independent blocks in order to consider the main conditional dependencies. Note that this idea was introduced in the Multimix software by [START_REF] Jorgensen | Mixture model clustering of data sets with categorical and continuous variables[END_REF] to cluster continuous and categorical data. Indeed, this software splits the variables into conditionally independent blocks but each block is allowed to contain one categorical variable at most (and a free number of continuous variables). Moreover, the categorical variable of a block is modelled by a full multinomial distribution. The CMM approach is developed only for categorical variables and considers blocks with more than one categorical variable. Moreover, each block of CMM follows a multinomial distribution per modes which assumes that few levels, named modes, are characteristic whereas the other ones follow a uniform distribution. Thus, the resulting multinomial distribution is parsimonious since its free parameters are limited to these few modes. Finally, the Multimix software carries out the model selection with a forward procedure while the model selection of CMM is achieved by a hybrid MCMC algorithm, inheriting thereby some optimality properties.

The resulting model is a good challenger for CIM since it preserves the sparsity and avoids many biases by considering the main conditional correlations. This model can also be interpreted as a parsimonious version of the log-linear mixture model. Indeed, the repartition of the variables into blocks defines the considered interactions while the multinomial distribution per modes defines specific intra-block interactions. In this way, it produces meaningful classes since the intra-class dependencies are brought out at two complementary levels: the block variable interaction level and the associated mode interaction level (through both locations and probabilities). Note that even if CMM can model some intra-class dependencies, it can also assume the conditional independence between variables. In this way, it generalizes lots of the parsimonious versions of CIM [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF].

For a fixed number of components, the model selection (repartition of the variables into blocks and mode numbers) is the most challenging problem since the number of competing models is huge. Therefore, the model selection is carried out by an MCMC algorithm whose the mode of the stationary distribution corresponds to the model having the highest posterior probability. This algorithm performs a random walk in the model space and requires the computation of the integrated complete-data likelihood. This quantity is not approximated by BIC-like methods since their results are biased (see our numerical experiments). Indeed, the integrated complete-data likelihood is accessible and non ambiguous through weakly informative conjugate prior. This approach provides an efficient model selection in a reasonable computational time since the parameters are estimated via an EM algorithm only for the single selected model. Thus, this approach is a possible answer to the combinatorial model selection problem which is known to be a real challenge for a log-linear mixture model. This paper is organized as follows. Section 2 presents the new mixture model. Section 3 presents the hybrid MCMC algorithm which performs the model selection. Section 4 presents the EM algorithm used to perform the maximum likelihood inference. In Section 5, numerical experiments show the relevance of the proposed criterion for model selection and the properties of the estimation algorithms. Section 6 presents three clusterings of real data sets performed by the R package CoModes. A conclusion is given in Section 7.

Conditional modes model

2.1 Mixture model of conditionally independent blocks Observations x = (x 1 , . . . , x n ) consist of n individuals x i = (x 1 i , . . . , x p i ) described by p categorical variables, called thereafter initial variables.

Running example

The Alzheimer data set [START_REF] Moran | Syndromes of behavioural and psychological symptoms in mild alzheimer's disease[END_REF] indicates the presence (yes) or absence (no) of six symptoms of Alzheimer's disease (AD) for n = 240 patients diagnosed with early onset AD conducted in the Mercer Institute in St. James's Hospital, Dublin. The p = 6 binary symptoms are, in this order: hallucination (Hal), activity (Act), aggression (Agg), agitation (Agi), diurnal (Diu) and affective (Aff).

The conditional modes model (CMM) assumes that individuals arise independently from a mixture with g components where variables are grouped into d conditionally independent blocks of variables. The repartition of the variables into the d blocks 2 is defined by the partition s = (s 1 , . . . , s d ) of {1, . . . , p}. The new variable resulting from the concatenation of the initial variables affiliated to block j (i.e. {x b i ; b ∈ s j }) is itself a categorical variable whose levels are defined by the Cartesian product of the variables affiliated to block j. This new (block dependent) categorical variable defined for block j is denoted by x j i = (x jh i ; h = 1, . . . , m j ), such as x jh i = 1 if individual i has level h and x jh i = 0 otherwise, where m j is the number of levels for block j. Since this mapping of the variables is bijective, defining a probability on x i = (x 1 i , . . . , x p i ) is equivalent to defining a probability on x i = (x 1 i , . . . , x d i ). The running example below provides an illustration of the relations between the initial variables x b i and the new (block dependent) variables x j i .

Running example If s = ({1}, {2}, {3, 5}, {4}, {6}), then the mixture components consider the dependency between the variable aggression and the variable diurnal (initial variables 3 et 5) while the other variables are conditionally independent. The repartition of the variables into blocks given by s provides five new categorical (block dependent) variables x j i (with j = 1, . . . , 5). Only one variable is affiliated to blocks 1, 2, 4 and 5, so the variables of these blocks are binary ones, the variable of block 3 takes 4 levels described in Table 1. The distribution of the (variables associated to the) block j under component k has u kj degrees of freedom and is parametrized by (δ kj , α kj ), k designating the index of the component among the g ones. We detail in Section 2.2 the exact expression of this distribution and in particular the specific definition of its two important parameters δ kj and α kj . The probability distribution function (pdf) of CMM is

p(x i |m, θ) = g k=1 π k d j=1 p(x j i |u kj , δ kj , α kj ), (1)
where m = (g, s, u) specifies the model with u = (u kj ; k = 1, . . . , g; j = 1, . . . , d), and where θ = (π, δ, α) denotes the parameters with δ = (δ kj ; k = 1, . . . , g; j = 1, . . . , d), α = (α kj ; k = 1, . . . , g; j = 1, . . . , d) and π = (π 1 , . . . , π g ) groups the mixture proportions with 0 < π k ≤ 1 and

g k=1 π k = 1.

Block of multinomial distribution per modes

The CMM model considers that each block follows a so-called multinomial distribution per modes. This distribution has only few free parameters corresponding to its modes, while the other parameters are equal. Thus, the free parameters are those of the levels having the greatest probabilities whereas uniformity holds for non-mode levels. The distribution of block j for component k has u kj degrees of freedom (so u kj modes) with 0 ≤ u kj < m j , and its mode locations are defined by the discrete parameter δ kj = {δ kjh ; h = 1, . . . , u kj } with δ kjh ∈ {1, . . . , m j } and δ kjh = δ kjh if h = h . Its probabilities are given by α kj = (α kjh ; h = 1, . . . , u kj +1) where α kjh denotes the probability of mode h for h = 1, . . . , u kj and where α kju kj +1 corresponds to the remaining probability mass. So, α kj is defined on a truncated simplex denoted by S(u kj , m j ) with

S(u kj , m j ) =    α kj : u kj +1 h=1 α kjh = 1 and for 1 ≤ h ≤ u kj , α kjh ≥ α kju kj +1 m j -u kj > 0    . (2)
Therefore, the pdf of block j for component k is

p(x j i |u kj , δ kj , α kj ) = u kj h=1 (α kjh ) x j(k,h) i α kju kj +1 m j -u kj 1-h∈δ kj x j(k,h) i , (3)
where the implicit notation x

j(k,h) i = x
jδ kjh i is used.

Running example

We consider the CMM model defined by m = (g, s, u) where g = 2, s = ({1}, {2}, {3, 5}, {4}, {6}) and where u is defined by Table 2. For instance, for component 1, uniformity holds in block 2 since u kj = 0

Block 1 Block 2 Block 3 Block 4 Block 5 Component 1 1 0 1 1 0 Component 2 1 1 1 0 1
Table 2: Example of the number of modes for the Alzheimer data set.

while the distribution of block 3 has one mode and can be summarized by Table 3. Table 3: Summary of the distribution of block 3 for component 1.

Model characteristics

The proposed model has two levels of interpretation. First, the intra-class dependencies of variables (equal over the mixture components) are brought out by the repartition of the variables into blocks given by s. Second, the intra-class and intra-block dependencies are summarized by the modes (possibly different over the mixture components) where the locations are given by δ kj and where the probabilities are given by α kj . The modes are interpreted as an overcontribution in comparison to the uniform distribution, since the probability of each mode is greater than the probability of the other locations (see constraints of ( 2)). A shorter summary for each distribution is also available by using the following two compact terms defined on [0, 1] and which reflect respectively the complexity and the strength of the intra-class and intra-block dependencies:

ūkj = u kj m j -1 and α kj• = u kj h=1 α kjh (4)
Thus, the smaller is ūkj and the larger is α kj• , the more massed in few characteristic levels is the distribution.

Running example For block 3 of component 1, ū13 = 1/3 and α 13• = 0.82, so this block distribution is strongly concentrated on one level.

The CIM model and its parsimonious versions implemented in Rmixmod [START_REF] Lebret | Rmixmod: The R Package of the Model-Based Unsupervised, Supervised and Semi-Supervised Classification Mixmod Library[END_REF] belong to the family of CMM. Indeed, if each block of CMM is composed by a single observed variable then the CMM model considers independence between all the variables within class (often called local independence). Table 4 presents the link between both model families (notations [ε jh k ] and [ε j k ] designate models used in [START_REF] Lebret | Rmixmod: The R Package of the Model-Based Unsupervised, Supervised and Semi-Supervised Classification Mixmod Library[END_REF]; see more details in this paper). It is also very interesting to notice that CMM

Rmixmod model CMM d s u kj ūkj α kj• Free [ε jh k ] p ({1}, . . . , {p}) m j -1 1 [1 -1 m j ; 1] Constrained [ε j k ] p ({1}, . . . , {p}) 1 1 m j -1 [ 1 m j ; 1]
Table 4: Link between models of Rmixmod and CMM. can be more parsimonious than the locally independent model although it takes into account the conditional dependencies. Finally, note that a sufficient condition for the generic identifiability of CMM (see Appendix A) results from the property of the generic identifiability of CIM [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF].

3 Model selection by a hybrid MCMC algorithm

Model selection by integrated likelihood

The number of mixture components for the competing models is usually bounded by a value g max . The goal of model selection is to find, among the set of competing models M, the model m having the largest posterior probability. Thus, m = arg max m∈M p(m|x).

(5)

We denote by M g the space of the models having g components, therefore m = arg max g∈{1,...,gmax}

p(m g |x) with m g = arg max m∈Mg p(m|x). (6)
By assuming that the model prior follows a uniform distribution,

p(m|x) ∝ p(x|m) with p(x|m) = p(x|m, θ)p(θ|m)dθ, (7)
where p(x|m, θ) = n i=1 p(x i |m, θ) is the likelihood function and p(θ|m) is the parameter distribution of the parameters. We assume independence between the prior distributions of the parameters, so

p(θ|m) = p(π|m) g k=1 d j=1 p(δ kj |m)p(α kj |m). (8)
The conjugate prior distribution of a multinomial distribution is a Dirichlet distribution. Thus, the Jeffreys non informative prior of the mixture proportions is the Dirichlet distribution D g ( 1 2 , . . . , 1 2 ). Moreover, δ kj |m follows a uniform distribution over the subset of {1, . . . , m j } of size u kj . Finally, α kj |m follows a Dirichlet distribution D u kj +1 (1, . . . , 1) restricted on the space S(u kj , m j ). Note that parameters of this last distribution are chosen equal to one in order to facilitate the computation of the integrated complete-data likelihood (see Section 3.4).

The number of competing models is huge, so an exhaustive approach is not doable. Therefore, the estimation of m g is performed by an MCMC algorithm for g = 1, . . . , g max .

Model selection via an ideal Gibbs sampler

For any g, the estimation of m g is equivalent to the estimation of the couple (s, u) maximizing p(s, u|x, g). This aim can be achieved by a Gibbs sampler having p(s, u|x, g) as marginal stationary distribution. Thus, we return the most generated couple (s, u).

We introduce the instrumental variable z = (z i ; i = 1, . . . , n) which indicates the class membership of the individuals, where z i = (z ik ; k = 1, . . . , g) with z ik = 1 if individual i belongs to component k and z ik = 0 otherwise. The Gibbs sampler alternates between the conditional sampling of the partition z and the conditional sampling of the couple (s, u). Thus, its stationary distribution is p(s, u, z|x, g). Its iteration [q] is written as

z [q+1] ∼ p(z|g, s [q] , u [q] , x) (9)
(s [q+1] , u [q+1] ) ∼ p(s, u|g, x, z [q+1] ). ( 10)

Both previous steps are performed with difficulty because independence between individuals does not hold in (9) and because the space of (s, u) is huge in (10). Thus, an hybrid algorithm derived from this Gibbs sampler is used. The mode of its marginal stationary distribution stays located at m g . The samplings from ( 9) and ( 10) are detailed in Section 3.3 and Section 3.4 respectively.

Gibbs algorithm for partition sampling

In this section, we show that the sampling from ( 9) can be achieved by a Gibbs sampler which alternates between both following samplings: z|m, x, θ and θ|m, x, z. Both previous samplings can be easily defined by introducing the complete-data likelihood function, which is written as

p(x, z|m, θ) = g k=1 π n k k d j=1 u kj h=1 (α kjh ) n kj(h) α kju kj +1 m j -u kj nkj 1 {α kj ∈S(u kj ;mj )} , (11)
where

n k = n i=1 z ik , n kj(h) = n i=1 z ik x j(k,h) i and nkj = n k - u kj h=1 n kj(h) .
The sampling from z|m, x, θ is performed by n independent samplings of the multinomial distributions M g (t i1 (θ|m), . . . , t ig (θ|m)) where

t ik (θ|m) = π k d j=1 p(x j i |u kj , δ kj , α kj ) g k =1 π k d j=1 p(x j i |u k j , δ k j , α k j ) . ( 12
)
The sampling from θ|m, x, z is performed as follows. The distribution of π|m, x, z is the Dirichlet D g 1 2 + n 1 , . . . , 1 2 + n g and the posterior distribution of (δ kj , α kj ) is written as

p(δ kj , α kj |m, x, z) = p(δ kj |m, x, z)p(α kj |m, x, z, δ kj ). ( 13
)
The distribution of δ kj |m, x, z is a multinomial one with too many values to be computable. Let the set δ kj = {δ kjh ; h = 1, . . . , u kj } containing the indices of the u kj largest values of

n kjh = n i=1 z ik x jh i such as ∀h ∈ δ kj , ∀h ∈ {1, . . . , m j } \ δ kj , n kjh ≥ n kjh . ( 14
)
We assume now that the difference between the mode probabilities and the non-mode probabilities are significant. So, the distribution of δ kj |m, x, z can be approximated by a Dirac in δ kj . This approximation is strengthened by the fast convergence speed of the discrete parameters [START_REF] Choirat | Estimation in discrete parameter models[END_REF]. Finally, α kj |m, x, z, δ kj follows the Dirichlet distribution D u kj +1 1+n kj(1) , . . . , 1+

n kj(u kj ) , 1 + n kj truncated on S(u kj ; m j ), where n kj(h) = n kjδ kjh and where

n kj = n k - u kj
h=1 n kj(h) .

MCMC algorithm for model sampling

The sampling from (10) requires the computation of the integrated completedata likelihood p(x, z|m) = p(x, z|m, θ)p(θ|m)dθ, since

p(s, u|g, x, z) ∝ p(x, z|m). ( 15
)
The integrated complete-data likelihood is written as

p(x, z|m) = Γ ( g 2 ) g k=1 Γ (n k + 1 2 ) Γ ( 1 2 ) g Γ (n + g 2 ) g k=1 d j=1 I kj (u kj ), ( 16
)
where I kj (u kj ) is the integral related to block j of component k defined by

I kj (u kj ) = δ kj S(u kj ,mj ) u kj h=1 (α kjh ) n kjδ kj (h) α kju kj +1 m j -u kj n k -h∈δ kj n kjh dα kj .
(17) The integral I kj (u kj ) has not a closed form. A BIC-like method could evaluate it by approximating the sum over the discrete parameters with its largest term and by approximating the integral with Laplace method. However, we propose an alternative which still neglects the sum over the discrete parameters of the modes locations but which now performs the exact computation on the continuous parameters. This approach is more precise and avoids the bias of the BIC criterion (see Section 5). The value of I kj (u kj ) is finally approximated by (proof is in Appendix B)

I kj (u kj ) ≈ 1 m j -u kj n kj u kj h=1 Bi 1 mj -h+1 ; n kj(h) + 1; n kjh + 1 m j -h , (18)
where n kjh = n k - It is not doable to compute the integrated complete-data likelihood for each (s, u) since the number of competing models is too huge. Thus, the sampling from ( 10) is performed by an MCMC algorithm, detailed in Appendix C, which has p(s, u|g, x, z) as stationary distribution and which also requires the computation of integral I kj defined by (17).

Maximum likelihood estimate

When model m has been assessed, its parameters θ m maximizing the likelihood function have to be obtained

θ m = arg max θ p(x| m, θ). ( 19
)
The direct optimization of the likelihood function involves to solve equations having no analytical solution. So, the parameter estimation is performed via an EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm[END_REF], which is often simple and efficient for missing data. This iterative algorithm alternates between two steps: the computation of the complete-data log-likelihood conditional expectation (E step) and its maximization (M step). Its iteration [r] is written as: E step: computation of the conditional probabilities t

[r] ik = t ik (θ [r] | m). M step: maximization of the complete-data log-likelihood π [r+1] k = t [r] •k n , δ [r+1] kj = arg max δ kj s.c. card(δ kj )=û kj h∈δ kj n i=1 t [r] ik xjh i , α [r+1] kjh =    1 t [r] •k n i=1 t [r] ik xjδ [r+1] kjh i if h = 1, . . . , ûkj 1 -ûkj h=1 α [r+1] kjh if h = ûkj + 1,
where xj i = (x jh i ; h = 1, . . . , m j ) is the variable resulting from the concatenation of the variables affiliated into block j by ŝ and where t

[r] •k = n i=1 t [r]
ik .

Simulations

Challenge of the mode number selection: Integrated likelihood vs BIC

Experiment design This experiment compares two approximations of I kj (u kj ): the proposed method defined by (18) and the BIC method [START_REF] Schwarz | Estimating the dimension of a model[END_REF] defined by

BIC(I kj (u kj )|z) = n i=1 g k=1 z ik ln p(x j i |u kj , δ kj , α kj ) - u kj 2 ln n k , (20)
where

α kjh = 1 n k n i=1 z ik x jδ kjh i if h = 1, .
. . , u kj , and α kju kj +1 = 1-ûkj h=1 α kjh . Note the approximations of I kj (u kj ) should be precise, since this quantity is used in the Gibbs sampler carrying out the model selection.

Samples composed of n i.i.d individuals are drawn from a multinomial distribution per modes M mj (

α kj• 3 , α kj• 3 , α kj• 3 , 1-α kj• mj -3 , . . . , 1-α kj• mj -3
) with m j levels and three modes having probability α kj• 3 . For different sizes of sample, 10 5 samples are generated with different values of (α kj• , m j ).

Results Figure 1 gives a comparison between (18) and the BIC-like approximation for the selection of the number of modes. The proposed criterion outperforms the BIC criterion in the four studied situations for large sample sizes. Indeed, its asymptotic behavior is better than the BIC criterion, since it rarely overestimates the mode number and since its variability is smaller than the BIC criterion. Since exact criteria are more efficient than approximated criteria (see [START_REF] Biernacki | Exact and Monte Carlo calculations of integrated likelihoods for the latent class model[END_REF]), this behavior was expected because the proposed criterion is closer to the exact value of the integrated likelihood than the BIC criterion. Indeed, the approximation of the integrated likelihood made by the proposed criterion is only made on the discrete parameters while the BIC criterion performs an approximation on both the discrete and the continuous parameters. We present now more specific comments. Fig. 1: Probability that the BIC criterion selects the true mode number (plain black line) and overestimates it (dotted black line) and probability that the proposed approach selects the true mode number (plain red line) and overestimates it (dotted red line).

In case (a), modes have a large probability mass and they are easily detected since there are few levels. Both criteria are efficient since they find the true number of modes with a probability close to one even for small samples. When the mode probabilities decrease (case (b)), it is more difficult to identify them. In such a case, the BIC criterion selects more often the true mode number with a moderate overestimation risk, for the small samples (size lower than 150), while the proposed approach can underestimate the number of modes. When the sample size is larger than 200, the proposed approach obtains better results since it finds a true number of modes almost always while the BIC criterion keeps an overestimation risk.

When the number of levels increases (case (c)), the problem is more complex and our approximation of the integrated complete-data likelihood given by ( 18) shows more clearly its interest. Indeed, the BIC criterion is strongly biased even for a large sample while the proposed approach almost always finds the true number of modes when the sample size is larger than 100. Finally, note that in the more complex situations like in case (d) (few probability mass for the modes and large number of levels), the approximation defined by ( 18) underestimates the mode number when the sample size is small. However, it converges to the true mode number when the sample size increases whereas the BIC criterion keeps significant bias even for a large data set.

Estimation accuracy with well specified model

Experiment design During this experiment, we highlight the relevance of the selected model m and of its parameter estimate θ m. For different values of n = (50, 100, 200, 400, 800), 25 data sets of six variables with three levels each are sampled from a bi-component CMM with the following parameters: s = ({1, 2}, {3, 4}, {5, 6}), u kj = 2, ∀(k, j) π = (0.5, 0.5), α kj = (0.4, 0.4, 0.2), δ 1j = {2, 4} and δ 2j = {6, 8}.

The estimation of the CMM model is performed with the R package CoModes (using the default option). Finally, the computing times are obtained on a 30 cores Intel(R) Xeon(R) CPU E5-4627 v2 @ 3.30GHz.

Results Table 5 presents the simulation results. The model selection works well since the true model is selected with a probability tending to one when n grows toward infinity. Moreover, the Kullback-Leibler divergence vanishes also when the sample size increases. Thus, the estimated distribution converges to the true one.

Estimation accuracy with misspecified model

Experiment design During this experiment, we underline the robustness of CMM. Thus, 25 samples of size 100 are generated by the following bi-component

; θ) = 0.5 3 h=1 p(x 2h-1 i , x 2h i ; θ) + 0.5 p(x 1 i ; θ)p(x 6 i ; θ) 2 h=1 p(x 2h i , x 2h+1 i ; θ), ( 21
) with p(x j i , x j+1 i ; θ) = p(x j i ; θ) λ1 {x j+1 i =x j i } +(1-λ)p(x j+1 i
; θ) and with x j i ∼ M 3 (1/3, 1/3, 1/3). Thus, when λ = 0, the sample is generated by a uniform distribution and classes are confused. The larger is the tuning parameter λ, the larger are the intra-class dependencies and the class separation.

Results Table 6 presents the Kullback-Leibler divergence of the model defined by ( 21) from CMM and CIM with their best number of classes (with g max = 4). λ 0.2 0.4 0.6 0.8 CMM KL-mean 0.08 0.24 0.44 0.64 KL-sd 0.02 0.01 0.07 0.11 g = 1 0.89 0.78 0.18 0.00 g = 2 0.11 0.22 0.72 0.68 g = 3 0.00 0.00 0.05 0.24 g = 4 0.00 0.00 0.05 0.08 CIM KL-mean 0.11 0.27 0.59 1.09 KL-sd 0.02 0.02 0.02 0.15 g = 1 1.00 1.00 0.89 0.44 g = 2 0.00 0.00 0.11 0.56 g = 3 0.00 0.00 0.00 0.00 g = 4 0.00 0.00 0.00 0.00

Table 6: Mean (KL-mean) and standard deviation (KL-sd) of the Kullback-Leibler divergence of the model defined by ( 21) from CMM and CIM and means of the selected number of classes (g).

The larger is λ, the larger is the Kullback-Leibler divergence for both models. However, the flexibility of CMM allows to keep an acceptable value of the Kullback-Leibler divergence while this divergence grows dramatically faster with CIM. Furthermore, when the classes are well separated (large value of λ), CMM finds more often the true class number than CIM.

Applications

This section illustrates the CMM model on three real data sets. Its results are obtained by the R package CoModes (using the default options) and are compared with the results of the CIM model obtained by the R package Rmixmod [START_REF] Lebret | Rmixmod: The R Package of the Model-Based Unsupervised, Supervised and Semi-Supervised Classification Mixmod Library[END_REF], using the default options too). Since, the first example is a real clustering problem, the partitioning accuracy of the competing models cannot be evaluated. For the other two data sets, a partition among the individuals is known and is used to compare the partitioning accuracy of the competing models. Finally, the computing times of CMM are obtained on a 30 cores Intel(R) Xeon(R) CPU E5-4627 v2 @ 3.30GHz. The computing times of CIM is not provided since it is very low compared to CMM (CIM is implemented in C++ for Rmixmod, CMM is only implemented in R).

Alzheimer clustering

Data We consider the Alzheimer data set introduced in the running example (see Section 2). We recall that this data set describes 240 patients with six binary variables (absence:no, presence:yes) indicating relevant symptoms of the Alzheimer disease. 7 summarizes the clustering results obtained by the CIM and the CMM models on the six features of the Alzheimer data set. The CMM model obtains a greater value of the BIC criterion than the CIM model independently of the number of components. These results were expected since the CIM model is included in the CMM model. Thus, the CMM fits the data better than the CIM model. Both models select two components but the CMM model also considers intra-component dependencies between the variables agression and diurnal. Table 8 shows that modelling the intra-component dependencies and introducing parsimonious constraints on the multinomial distributions can strongly impact the partition. 9. The model splits the population into two equally balanced classes (π 1 = 0.45 and π 2 = 0.55). The class 1 groups the individuals having few symptoms. Indeed, all the modes in this class indicate an absence of symptoms. Note that in this class, two blocks which are composed with a single variable (activity and affective) follow a uniform distribution. The class 2 contains the individuals presenting more symptoms. Indeed, two modes indicate a presence of symptoms activity and affective. However, the modes of blocks 1 and 3 are the same for both classes. However, their probabilities in class 2 are smaller than in class 1 (e.g. α 13• = 0.82 and α 13• = 0.44). Thus, the CMM model overall splits the population between two classes according to the stage of the disease. Table 10: Presentation of the five plumage and external morphological variables describing the puffins.

Model comparison Table

Model comparison Table 11 summarizes the clustering results obtained by the CIM and the CMM models on the five features of the Seabirds data set. Results

show that the CMM model allows to fit the data distribution better than the CIM model since it obtains greater values of the BIC criterion. Note that the parameter estimation of the CIM model is computed immediately but that the estimation of the CMM model is performed in a reasonable computing time (less than ten minutes for the whole analysis). Finally, the values of the adjusted Rand Index [START_REF] Hubert | Comparing partitions[END_REF] indicates that the partition resulting from the best CMM model (according to the BIC criterion) is closer to the partition defined by the subspecies than the partition resulting from the best CIM model. Table 12 presents the confusion matrices between the partitions resulting form the best CIM and CMM models and the subspecies. It shows that the Subalaris are more different than the two other subspecies. Indeed, both models assign all the Subalaris into the same class. Finally note that, the first principal correspondence axis allows to well-visualize the partition provided by the CMM model (see Figure 2). Best model interpretation The best CMM model and its parameters are displayed in Table 13. Thus, the variables collar, border and gender are considered as conditionally independent while the intra-component dependencies between the variables eyebrows and sub-caudal are modelled. Note that the block composed by these last two variables is strongly discriminative, since its modes Table 12: Confusion matrices between the subspecies and partitions with respect to the best CMM and CIM models. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 1.0 1.5

-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
Dim 1 (14.97%)

Dim 2 (10.63%) q q q dichrous lherminieri subalaris (a) Partition defined by the subspecies q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0.5 0.0 0.5 1.0 are different over the components and since they group a large probability mass for both components. Moreover, the variable gender follows a uniform distribution under each component (since there is an absence of mode), so it is not discriminative (the variable follows the same distribution over the mixture component). This results was expected since the gender should not provide any discriminative information regarding to the sub-specie of a seabird.

Acute inflammations clustering

Data The Acute data set [START_REF] Czerniak | Application of rough sets in the presumptive diagnosis of urinary system diseases. Artifical Intelligence and Security in Computing Systems[END_REF] Table 13: Parameters of the best CMM model for the seabirds data set. For each class, the first rows give the mode description (names of the levels h of the block j associated to their probabilty α kjh ) and the last row gives the indices of the block distribution. We recall that ūkj and α kj• stand respectively for the complexity and the strength as defined in ( 4).

(Lum), the urine pushing (Pus), the micturition pains (Mic) and the burning of urethra (Bur). Moreover, many patients have one of the following diseases of the urinary system: inflammation of urinary bladder (Uri) and Nephritis of renal pelvis origin (Ren). Therefore, there are three known partitions among the individuals: two two-group partitions defined by the variables Uri and Ren and one four-group partition defined by the Cartesian production between the variables Uri and Ren designed now by Uri/Ren. achieving the best value of the log-likelihood function produce systematically the same estimates.

Model comparison

Table 15: Parameters of the best CMM model for the Acute data set. For each class, the first rows give the mode description (names of the levels h of the block j associated to their probabilty α kjh ) and the last row gives the indices of the block distribution. We recall that ūkj and α kj• stand respectively for the complexity and the strength as defined in (4).

Table 15 summarizes the block distribution and lists the modes of each block. Based on these results, the class 1 is characterized by individuals with a high temperature, often lumbar pain and no nausea. On the other hand, the class 2 contains only individuals with a temperature equal or less than 38C. Note that the model provides strongly different classes. Indeed, block 1 provides different modes for both classes and the probability mass of the modes is high (α 11• = 0.90 and α 21• = 0.98). Moreover, the distributions of the two components for block 2 is mainly concentrated on the modes (α 12• = 0.98 and α 22• = 0.99) and only two modes appear in both components.

It is known [START_REF] Czerniak | Application of rough sets in the presumptive diagnosis of urinary system diseases. Artifical Intelligence and Security in Computing Systems[END_REF] that acute nephritis of renal pelvis origin begins with sudden fever, which reaches, and sometimes exceeds 40C. The fever is accompanied by shivers and one or both-side lumbar pains, which are sometimes very strong. Symptoms of acute inflammation of urinary bladder appear very often. Quite not infrequently there are nausea and vomiting and spread pains of whole abdomen. Therefore, the interpretation provided by the bi-component CMM model is relevant. As shown by Table 16, the class 1 mainly groups individuals with acute nephritis of renal pelvis origin.

CMM class1 class2

Ren -no 10 60 Ren -yes 50 0

Table 16: Confusion matrix between the acute nephritis of renal pelvis origin and the partition resuling from the best model.

Conclusion

In this article, we have presented a new mixture model (CMM) to cluster categorical data. Its strength is to relax the conditional independence assumption while staying parsimonious. A summary of the distribution of variables is given by both indices ūkj (complexity) and α kj• (strength) while each class can be summarized by the mode locations. As shown on the Seabirds application, CMM can improve the results of the classical latent class model even if the conditional independence assumption is true, thanks to its sparsity. The combinatorial problems of the block detection and of the modes number selection is solved by a hybrid MCMC algorithm which uses the computation of the integrated complete-data likelihood and which does not require estimates. Thus, this approach can be used to select the interactions of the log-linear mixture model per block. The parameters are only estimated for the single selected model. The R package CoModes allows to perform the model selection and the parameter estimation. Both data sets presented in this article are included in this package. To efficiently reduce the computing time, the functions of this package will be soon implemented in C++.

However, the model selection becomes difficult if the data set has a large number of variables, since the number of competing models becomes large. Some constraints on the block variables repartition could also be added (for instance the number of variables into blocks could be limited at three variables).

Another solution could be to estimate the model by a forward/backward strategy but it is known that these methods are sub-optimal.

Finally, we imposed the equality of the repartition of the variables into blocks for all the classes. This property allows us to prove the generic identifiability of CMM. This lack of flexibility is counterbalanced by flexible block distribution. However, one could try to relax the class-equality of s with the model non-identifiability risk.

A Generic identifiability of CMM

When d ≥ 3, CMM is generically identifiable (i.e. the parameter space where the model is not identifiable has a Lebesgue measure equal to zero). The demonstration is based on results of [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF] which use the Kruskal theorem [START_REF] Kruskal | Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics[END_REF][START_REF] Kruskal | More factors than subjects, tests and treatments: An indeterminacy theorem for canonical decomposition and individual differences scaling[END_REF]. The demonstration is cut into three steps: reminder of the Kruskal results for the three-way tables, demonstration of the model generic identifiability when d = 3, then its extension when model has more than three blocks.

Kruskal results

For a matrix M , the Kruskal rank of M , denoted by rank K M is the largest number I such that every set of I rows of M are linearly independent. Theorem 1 (Kruskal [START_REF] Kruskal | Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics[END_REF][START_REF] Kruskal | More factors than subjects, tests and treatments: An indeterminacy theorem for canonical decomposition and individual differences scaling[END_REF]). Let I j = rank K M j . If Corollary 1 The parameters of CMM with three blocs are generically identifiable, up to label swapping, provided: min(g, ξ 1 ) + min(g, ξ 2 ) + min(g, ξ 3 ) ≥ 2g + 2.

Generic identifiability of CMM with more than three blocks In the same way that [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF], we generalize the result with d blocks by observing that d blocks of categorical variables can be combined into three categorical variables. Thus, we can apply the Kruskal theorem.

Corollary 2 We consider a CMM with d blocks where d ≥ 3. If there exists a tri-partition of the set {1, . . . , d} into three disjoint non empty subsets S 1 , S 2 and S 3 , such that γ i = j∈S i ξ j with min(g, γ 1 ) + min(g, γ 2 ) + min(g, γ 3 ) ≥ 2g + 2, (

then the model parameters are generically identifiable up to label swapping.

B Approximation of I(u kj )

First, we define a new parametrization of the block distribution facilitating the integrated complete-data likelihood computation and the prior distribution related to this new block parametrization. Second, we underline the relationship between the embedded models. We conclude by the integrated complete-data likelihood computation, which is the target result.

B.1 New parametrization of the block distribution

Without loss of generality, we assume that the elements of δ kj are ordered by decreasing values of the probability mass associated to them and we introduce the new parametrization of α kj denoted by ε kj = (ε kjh ; h = 1, . . . , u kj ) where ε kj ∈ E kj = 1 m j ; 1 × . . . × 1 m j -u kj ; 1 and where ε kjh is defined by

ε kjh = α kj1 if h = 1 α kjh h-1 h =1
(1-ε kjh )

otherwise.

Each ε kjh follows a truncated beta distribution on the interval 1 m j -h+1 , 1 , so

p(ε kj |δ kj , m) = m j m j -u kj . ( 24
)
The conditional probability of x j = (x j i ; i = 1, . . . , n) is p(x j |z, u kj , δ kj , ε kj ) =

u kj h=1 (ε kjh ) n kj(h) (1 -ε kjh ) n kjh . (25)

B.2 Relation between embedded models

Let the model with u kj modes and the parameters (δ kj , ε kj ) and the model with u kj modes and the parameters (δ kj , ε kj ) such as u kj = u kj -1 and such as the u kj modes having the largest probabilities have the same locations (∀h ∈ δ kj , h ∈ δ kj ) and the same probability masses (ε kjh = ε kjh , h < u kj ). These embedded models follow this relation p(x j |z, u kj , δ kj , ε kj ) p(x j |z, u kj , δ kj , ε kj ) = (m j -u kj + 1)

B.3 Integrated complete-data likelihood

The integrated complete-data likelihood is finally approximated, by neglecting the sum over the discrete parameters of the modes locations and by performing the exact computation on the continuous parameters, by

I kj (u kj ) ≈ 1 m j -u kj n kj u kj h=1 Bi 1 m j -h+1 ; n kj(h) + 1; n kjh + 1 m j -h , (27)
Proof If, for the model with u kj -1 modes, the best modes locations are known and given by δ kj then the conditional probability of x j for a model with u kj modes is p(x j |z, u kj , δ kj , ε kj ) = 1 m j -u kj + 1 δ kj ∈∆ kj p(x j |z, u kj , δ kj , ε kj ),

where ∆ kj = {δ kj : δ kj ⊂ δ kj and card(δ kj ) = u kj }. Thus, by approximating this sum by its maximum element, we obtain that p(x j |z, u kj , δ kj , ε kj ) = 1 m j -u kj + 1 p(x j |z, u kj , δ kj , ε kj ), (29)

hh

  =1 n kj(h ) and Bi(x; a, b) = B(1; a, b) -B(x; a, b), B(x; a, b) denoting the incomplete beta function defined by B(x; a, b) = x 0 w a (1w) b dw.

  α kj• = 0.3, m j =27.

Fig. 2 :

 2 Fig.2: Seabirds on the first principal correspondence analysis map where an i.i.d. uniform noise on [0, 0.1] has be added on both axes for each individual to improve visualization.

I 1 +

 1 I 2 + I 3 ≥ 2g + 2, then the tensor [M 1 , M 2 , M 3 ] uniquely determines the M j , up to simultaneous permutation and rescaling rows. Generic identifiability of CMM with three blocks Let k 0 = argmin k u kj and the matrix M j where M j (k, h) = α kjh .(22)By denoting by ξ j = min k u kj + 1, generically, we have rank K M j = min(g, ξ j ).

Table 1 :

 1 Example of the new 4-level variable associated to block 3 and composed of two initial binary variables.

	observed variables	block 3
	x 3 i	x 5 i	x 3 i
	yes	yes	level 1
	yes	no	level 2
	no	yes	level 3
	no	no	level 4

Table 5 :

 5 Percentage (over the 25 samples) where ŝ corresponds to s and where û corresponds to u. Mean and standard deviation (sd) of the Kullback-Leibler (KL) divergence of the true model from the estimated distribution.

	mixture model of dimension six where the intra-class dependencies are different
	for both components
	p(x i

Table 7 :

 7 Results of CMM and CIM on the Alzheimer data set: BIC criterion (BIC) and computing time in minutes (Time). Values of the best model are in bold.

					number of components
			1	2	3	4	5	6
	CIM	BIC	-789.37 -785.13	-801.86 -818.02 -834.94 -855.88
	CMM	BIC	-779.22 -772.57	-785.38 -796.48 -810.81 -815.68
		Time	0.78	1.16	1.70	2.11	2.61	2.92

Table 8 :

 8 Confusion matrix between the partitions resulting from the CMM and the CIM models obtained on the Alzheimer data set.Best model interpretation The model selected by the BIC criterion is the model described in the running example (see Section 2): i.e. bi-component CMM model with s = ({1}, {2}, {3, 5}, {4}, {6}). Its parameters are summarized in Table

Table 11 :

 11 Results of CMM and CIM on the Seabirds data set: BIC criterion (BIC), adjusted Rand Index (ARI) and computing time in minutes (Time). Values of the best model are in bold.

					number of components
			1	2	3	4	5	6
	CIM	BIC	-711.59 -707.14	-725.21 -749.64 -782.26 -811.15
		ARI	0.00	0.13	0.19	0.30	0.27	0.26
		BIC	-701.98 -677.37	-693.74 -697.08 -705.46 -716.05
	CMM	ARI	0.00	0.21	0.10	0.07	0.04	0.04
		Time	0.56	0.95	1.30	1.64	1.93	2.22

  describes 120 patients by six features. The first variable is continuous and indicates the temperature of the patient (Temp), so only its integer part is used in order to set into a categorical variable with seven levels (integer from 36 to 42). The remaining five binary variables indicate the occurrence of nausea (Nau), the lumbar pain

	Class	Block 1		Block 2		Block 3	Block 4
		col.	α kjh	eye.	sub.	α kjh	bor.	α kjh	gen.	α kjh
	1	none	0.52	vis.	white	0.76	none	0.96
		dash	0.24	pron. white	0.10
		al.cont.	0.21
	ū1j α 1j•	0.75	0.97	0.13	0.86	0.50	0.96	0.00	0.00
	2	trace	0.46	vis.	black	0.43	none	0.74
		dash	0.31	trace	white	0.21	few	0.25
		al.cont.	0.18	trace	black	0.16
	ū2j α 2j•	0.75	0.95	0.20	0.80	1.00	0.99	0.00	0.00

Table 14 :

 14 Table14summarizes the clustering results obtained by the CIM and the CMM models on the six features of the Acute data set. The CMM model fits better the distribution of the observed variables since its BIC criterion is greater than this one of the CIM model. About the partitioning accuracy, both models provide relevant partitions and their interpretations are complementary. The best CIM model produces a partition strongly similar to the partition in four groups defined by the variable Uri/Ren (ARI=0.85). On the other hand, the best CMM model produces a partition more similar to the two-group partition defined by the absence/presence of Nephritis of renal pelvis origin (ARI=0.56). model interpretation According to the BIC criterion, the best model is the bi-component CMM model composed of two blocks where the first one only contains the variable Temp while the second one contains all the other variables. For both components, the block 1 follows a multinomial distribution with three modes while the block 2 follows a multinomial distribution with six modes. Thus, the model considers within-class dependencies between all variables, Temp excluded, and it requires only 19 parameters (while the bicomponent CIM model needs 23 parameters). However, since this model has only two blocks, we cannot prove its identifiability (see Appendix A). Nevertheless, this model seems to be identifiable because the runs of the EM algorithm Results of the Acute data set: BIC criterion (BIC), adjusted Rand Index (ARI) and computing time in minutes (Time). Values of the best models according to the BIC value are in bold.

	Best

Downloadable at https://r-forge.r-project.org/R/?group id=1809

Note that the repartition of the variables into blocks is identical between classes. This choice is motivated by reasons of identifiability and interpretation that we will detail later.

Class

Block For each class, the first row gives the mode description (names of the levels h of the block j associated to their probabilty α kjh ) and the last row gives the indices of the block distribution. We recall that ūkj and α kj• stand respectively for the complexity and the strength as defined in (4).

Seabirds clustering

Data The Seabirds data set is a biological data set describing 153 puffins (seabirds) by five plumage and external morphological characteristics presented in Table 10 ( Bretagnolle, 2007). These seabirds are divided into three subspecies dichrous (84 birds), lherminieri (34 birds) and Subalaris (35 birds).

As p(x j |z, u kj = 0) = (m j ) -n k , by applying recursively (26), we obtain that

C Details on the model sampling

The sampling of (s [s+1] , u [s+1] ) from ( 10) is performed in two steps. Firstly, a new repartition of the variables into blocks and the mode number of the modified blocks, respectively denoted by s [s+1] and u [s+1/2] , are sampled by one iteration of a Metropolis-Hastings algorithm. Secondly, the mode number of each block is sampled by one MCMC iteration. Thus, the sampling of m [s+1] is decomposed into the two following steps

Thus, this chain has p(s, u|g, x, z [s+1] ) as stationary distribution.

C.1 Metropolis-Hastings algorithm to sample from (31)

This sampling is performed by one iteration of the Metropolis-Hastings algorithm divided into two steps. Firstly, the proposal distribution q(.; m [s] ) generates a candidate m = (g, s , u ). Secondly m [s+1] is sampled according to the acceptance probability µ [s] defined by

Note that the computation of µ [s] involves to compute the integrated complete-data likelihood defined by ( 17). The sampling of m [s+1/2] is written as

m [s+1/2] = m with a probability µ [s] m [s] otherwise.

The proposal distribution q(.; m [s] ) samples m in two steps. The first step changes the block affectation of one variable. In practice, s is uniformly sampled in V (s [s] ) = {s : ∃!b as b ∈ s C.2 MCMC algorithm to sample u [s+1] This step allows to increase or decrease the mode number of each block by one at each iteration. So, u (35)