Ph.D Matthieu Marbac
email: matthieu.marbac-lourdelle@inria.fr

Christophe Biernacki
email: christophe.biernacki@math.univ-lille1.fr

Vincent Vandewalle
email: vincent.vandewalle@univ-lille2.fr

Model-based clustering for conditionally correlated categorical data

Keywords: Categorical data, Clustering, Correlation, Expectation-Maximization algorithm, Gibbs sampler, Mixture model, Model selection

An extension of the latent class model is presented for clustering categorical data by relaxing the classical "class conditional independence assumption" of variables. This model consists in grouping the variables into inter-independent and intra-dependent blocks, in order to consider the main intra-class correlations. The dependency between variables grouped inside the same block of a class is taken into account by mixing two extreme distributions, which are respectively the independence and the maximum dependency. When the variables are dependent given the class, this approach is expected to reduce the biases of the latent class model. Indeed, it produces a meaningful dependency model with only a few additional parameters. The parameters are estimated, by maximum likelihood, by means of an EM algorithm. Moreover, a Gibbs sampler is used for model selection in order to overcome the computational intractability of the combinatorial problems involved by the block structure search. Two applications on medical and biological data sets show the relevance of this new model. The results strengthen the view that this model is meaningful and that it reduces the biases induced by the conditional independence assumption of the latent class model.

Introduction

Nowadays practitioners often face very large data sets, which are difficult to analyze directly. In this context, clustering [START_REF] Jajuga | Classification, clustering and data analysis: recent advances and applications[END_REF] is an important tool providing a partition among individuals. Other approaches cluster both individuals and variables simultaneously [START_REF] Govaert | Clustering with block mixture models[END_REF]. Furthermore, with the increasing number of variables at hand, the risk of observing correlated descriptors, even within the same class, is often high. In view of these difficulties, the practitioner can choose between two approaches. The first approach consists in performing a selection among the observed variables [START_REF] Maugis | Variable selection in model-based clustering: A general variable role modeling[END_REF] in order to extract uncorrelated data, thereby losing some potentially crucial information. The second approach consists in applying a method for modeling the conditional dependencies on the whole set of variables.

Clustering methods can be divided into two kinds of approaches: geometric ones based on the distances between individuals and probabilistic ones modelling the data generation process.

Geometric methods are generally faster than probabilistic ones, but they are often quite sensitive to the choice of distance between individuals. Furthermore, as the probabilistic tools are not available for these methods, difficult questions, e.g. the selection of the number of clusters cannot be addressed rigorously. For categorical data, geometric methods either define a metric in the initial variables space e.g. k-means [START_REF] Huang | Automated variable weighting in k-means type clustering[END_REF], or compute their metric on the axes of the multiple correspondence analysis [START_REF] Chavent | A partitioning method for the clustering of categorical variables[END_REF][START_REF] Guinot | Use of multiple correspondence analysis and cluster analysis to study dietary behaviour: food consumption questionnaire in the SU. VI. MAX. cohort[END_REF].

Many geometric approaches can also be interpreted in a probabilistic way. Thus, for the continuous data, the classical k-means algorithm can be identified as an homoscedastic Gaussian mixture model [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF][START_REF] Celeux | Gaussian parsimonious clustering models[END_REF] with equal proportions. For the categorical variables, [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF] show that the CEM algorithm [START_REF] Mclachlan | The EM algorithm[END_REF], applied to a classical latent class model, maximizes a classical information criteria close to a χ 2 metric. Other links between the two approaches are described in Govaert (2010), Chapter 9. Let us now introduce our proposal for this problem.

In the categorical case, the latent class model also known as naive Bayes belongs to the folklore [START_REF] Goodman | Exploratory latent structure analysis using both identifiable and unidentifiable models[END_REF][START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF]. In this article, we refer to this model as the conditional independence model (further denoted by cim). Classes are explicitly described by the probability of each modality for each variable under the conditional independence assumption.

The sparsity of the model, implied by this assumption, is a great advantage since it restricts the curse of dimensionality. cim has obtained quite good results in practice in different areas [START_REF] Hand | Idiot's Bayes Not So Stupid after All[END_REF]), e.g. behavioral science [START_REF] Reboussin | A latent class analysis of underage problem drinking: Evidence from a community sample of 16-20 year olds[END_REF] and medicine [START_REF] Strauss | Using latent class analysis to identify patterns of hepatitis C service provision in drug-free treatment programs in the US[END_REF].

However, cim may suffer from severe biases when the data are intra-class correlated. For instance, an application presented by [START_REF] Van Hattum | Market Segmentation Using Brand Strategy Research: Bayesian Inference with Respect to Mixtures of Log-Linear Models[END_REF] shows that cim over-estimates the number of clusters when the conditional independence assumption is violated. For a long time, people have tried to relax the conditional independence assumption by modeling conditional interactions between variables using an additive model [START_REF] Harper | Local dependence latent structure models[END_REF]. The main drawback of this approach is that its number of parameters becomes huge, making their estimation intractable. Some other methods take into account the intra-class dependencies as mixtures of Bayesian networks [START_REF] Cheng | Comparing Bayesian network classifiers[END_REF]. Conditionally on each class, a directed acyclic graph is built with a set of nodes representing each variable. However, if no constraint is added, the network estimation is also quite complex. By constraining the network to be a tree, the model selection and the parameter estimation can be easily performed. Moreover the correlation model enjoys great flexibility. The extension of the dependency tree of [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF] was done by [START_REF] Friedman | Bayesian network classifiers[END_REF] for the supervised classification and by [START_REF] Meila | Learning with mixtures of trees[END_REF] for the clustering.

However the main problem of these models is that they too often require too many parameters.

When covariates are available, the conditional dependencies between the categorical ones can be modeled by a logistic function [START_REF] Formann | Linear logistic latent class analysis for polytomous data[END_REF][START_REF] Reboussin | Locally dependent latent class models with covariates: an application to under-age drinking in the USA[END_REF]. By assuming that these covariates are unobserved, the multilevel latent class model (Vermunt, 2003(Vermunt, , 2007) ) naturally incorporates the intra-class dependencies. This model has connections with the approach of [START_REF] Qu | Random Effects Models in Latent Class Analysis for Evaluating Accuracy of Diagnostic Tests[END_REF] where the intra-class dependencies are modeled by a latent continuous variable with a probit function. The hybrid model [START_REF] Muthén | Latent variable hybrids: Overview of old and new models[END_REF]) is more general approach in which, for each class, a factor analysis model is fitted to either all categorical variables or to those categorical variables having dependencies. Recently, [START_REF] Gollini | Mixture of latent trait analyzers for model-based clustering of categorical data[END_REF] proposed the mixture model of latent traits analyzers which assumes that the distribution of the categorical variables depends on both a categorical latent variable (the class) and many continuous latent traits variables. The parameter estimation is also a difficult point which is solved via a variational approach. All these models consider the intra-class dependencies, but their main drawback is that these dependencies must be interpreted among relations with a latent variable. As a result, pertinent interpretation can be difficult.

The log-linear models [START_REF] Agresti | Categorical data analysis[END_REF][START_REF] Bock | Loglinear models and entropy clustering methods for qualitative data. Classification as a tool of research[END_REF] were originally proposed to model the logprobability of each individual by selecting interactions between variables. The most general mixture model is the log-linear mixture model as it is able to incorporate many forms of interactions.

It has been used at least since [START_REF] Hagenaars | Latent structure models with direct effects between indicators local dependence models[END_REF]. [START_REF] Espeland | Using Latent Class Models to Characterize and Assess Relative Error in Discrete Measurements[END_REF] used it to cluster radiographic cross-diagnostics and [START_REF] Van Hattum | Market Segmentation Using Brand Strategy Research: Bayesian Inference with Respect to Mixtures of Log-Linear Models[END_REF] in a market segmentation problem. However this model family is huge and the model selection is a real challenge. In the literature, authors often require the modeled interactions ahead of time. Another option is to perform a deterministic search e.g. the forward method which is sub-optimal. Furthermore, the number of parameters to estimate increases with the conditional modalities interactions, thus implying potential over-fitting and more difficult interpretation. The latent class model (cim) can be seen as a particular log-linear mixture model, in which interactions are discarded. Our aim is to present a version of the log-linear mixture model which takes into account the interactions of order one or more while keeping the number of unknown parameters to a reasonable amount.

We propose to extend the classical latent class model (cim) for categorical data, by a new latent class model which relaxes the conditional independence assumption. We refer to this new model as the conditionally correlated model (denoted by cmm). This model is a parsimonious version of the log-linear mixture model, and thus benefits from its interpretative power. Furthermore, we propose a Bayesian approach to automatically perform model selection.

The ccm model groups the variables into conditionally independent blocks given the class.

The main intra-class dependencies are thus shown by the repartition of the variables into blocks. This approach, allowing modeling of the main conditional interactions, was first proposed by [START_REF] Jorgensen | Mixture model clustering of data sets with categorical and continuous variables[END_REF] in order to cluster continuous and categorical data. For cmm, each block follows a particular dependency distribution which corresponds to our main contribution. This distribution consists in a bi-component mixture of an independence and a maximal dependency distribution according to the Cramer's V criterion. This specific distribution of the blocks allows summarizing the conditional dependencies of the variables with only one continuous parameter: the maximum dependency distribution proportion. Thus, the model underlines the main conditional dependencies and their strength.

The new model is a two-degree parsimonious version of a log-linear mixture model. A first degree of parsimony is introduced by grouping the variables which are conditionally dependent into the same block. This repartition of the variables per blocks defines the interactions considered by the model for each class. Moreover, the strength of the correlation is reflected by the proportion of maximum dependency distribution. A second degree of parsimony is induced by the specific distribution of the blocks. As for all log-linear mixture models, the selection of the pertinent interactions is a combinatorial problem. We propose to perform this model selection via a Gibbs sampler in order to overcome the enumeration of all the models. Thus, this general approach could also select the interactions of any log-linear mixture model. This paper is organized as follows. Section 2 reviews the principles of the latent class model. Section 3 presents the new mixture model taking into account the intra-class correlations. Section 4 is devoted to parameter estimation in the case where the number of classes and the blocks of variables are supposed to be known. Section 5 presents a Gibbs algorithm for avoiding the combinatorial difficulties inherent to block selection. Section 6 presents results on simulated data.

Section 7 presents a comparison between two main model-based clustering approaches and our proposal on a classical medical data set, and presents another application on a larger real data set.

A tutorial of the R package Clustericat1 which performs the model selection and the estimation of the parameters of cmm is given with the first application (see Appendix A). Section 8 presents our conclusions.

Classical models

Latent class model: intra-class independence of variables

Observations to be classified are described with d discrete variables x = (x 1 , . . . , x d ) defined on the probabilistic space X . Each variable j has m j response levels with m j ≥ 2 and is written x j = (x j1 , . . . , x jm j ) where x jh = 1 if variable j takes modality h and x jh = 0 otherwise. In the standard latent class model (cim), the variables are assumed to be conditionally independent knowing the latent cluster. Furthermore, data are supposed to be drawn independently from a mixture of g multivariate multinomial distributions with probability distribution function (pdf)

p(x ; θ) = g k=1 π k p(x ; α k ) with p(x ; α k ) = d j=1 m j h=1 (α jh k ) x jh , (1)
with θ = (π 1 , . . . , π g , α 1 , . . . , α g ), π k being the proportion of the component k in the mixture where π k > 0 and g k=1 π k = 1, and α k = (α jh k ; j = 1, . . . , d; h = 1, . . . , m j ) where α jh k denotes the probability that variable j has level h if the object is in cluster k and satisfies the two following constraints: α jh k > 0 and m j h=1 α jh k = 1. The classical latent class model is much more parsimonious than the saturated log-linear model, which requires ( j m j ) -1 parameters, since it only requires ν cim parameters with

ν cim = (g -1) + g g k=1 (m j -1).
(2)

Its maximum likelihood estimator is easily computed via an EM algorithm [START_REF] Mclachlan | The EM algorithm[END_REF]. For the cluster analysis, the mixture identifiability up to a permutation of the class is generally necessary [START_REF] Mclachlan | Finite mixutre models[END_REF]. However, there are mixtures, such as the products of Bernoulli distributions, which are not identifiable but which produce good results in several applications. In order to relax this too stringent concept of identifiability, the notion of generic identifiability was introduced by [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF]: a model is generically identifiable if it is identifiable except for a subset of the parameter space with Lebesgue measure zero.

Latent class model extension: intra-class independence of blocks

Despite its simplicity, the latent class model s to good results in many situations [START_REF] Hand | Idiot's Bayes Not So Stupid after All[END_REF]. However, in the case of intra-correlated variables, it can entail severe biases in the partition estimation and it may also overestimate the number of clusters. In order to reduce these biases, a classical extension of the latent class model was introduced by [START_REF] Jorgensen | Mixture model clustering of data sets with categorical and continuous variables[END_REF] for conditionally correlated mixed data. This model is implemented in the Multimix software [START_REF] Hunt | Theory & Methods: Mixture model clustering using the MUL-TIMIX program[END_REF]. For each component k, each block b follows a specific parametric distribution denoted as p(x {kb} ; θ kb ) where θ kb groups the parameters of this distribution. The model pdf is written as

p(x ; σ, θ) = g k=1 π k p(x ; σ k , θ k ) with p(x ; σ k , θ k ) = B k b=1 p(x {kb} ; θ kb ), ( 3
)
where θ is redefined as θ = (π 1 , . . . , π g , θ 1 , . . . , θ g ) with θ k = (θ k1 , . . . , θ kb k ). Figure 1 independence, would be represented by white cells off the diagonal and black cells on the diagonal.

X 5 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 5 x {11} x {12} (a) X 3 X 4 X 2 X 5 X 1 X 1 X 5 X 2 X 4 X 3 x {21} x {22}
x {23}

(b) This approach is very general, since any distribution can be chosen for each block as soon as it is different from the distribution of independence. The mixture model by conditional independent blocks is a parsimonious version of the log-linear mixture model. Indeed, the distribution of variables in blocks determines which interactions are modeled. The interactions between variables of different blocks will be zero and those between variables of the same block can be modeled by the specific distribution of the block. The limiting case of this model where b k = d for each class is equivalent to the latent class model with the conditional independence assumption.

The generic identifiability of the mixture model with conditionally independent blocks follows, under specific constraints, from Theorem 4 of [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF] by assuming that the distribution of each block is itself identifiable. This proof is given in Appendix A of [START_REF] Marbac | Model-based clustering for conditionally correlated categorical data[END_REF].

Intra-block parsimonious distribution

The goal is now to define a parsimonious distribution for each block that takes into account the dependency between variables. Furthermore, the parameters of the distribution inside the block must be meaningful for the practitioner. In

Maximum dependency distribution

The maximum dependency distribution is defined as the "opposite" distribution of independence according to the Cramer's V criterion computed on all the variable couples. Indeed, this latter minimizes this criterion while the maximum dependency distribution maximizes it. Under this distribution, the modality knowledge of one variable provides the maximum information on all the subsequent variables. Note that it is a non-reciprocal functional link between variables. Indeed, if x {kb} arises from this distribution, the knowledge of the variable having the largest number of modalities determines exactly the others but the reverse does not necessarily apply. So, this distribution defines successive surjections from the space of x {kb}j to the space of x {kb}j+1 with j = 1, . . . , d {kb} -1 (recall that the variables are ordered by decreasing number of modalities in each block). In fact, it is a reciprocal functional link only when m h=1 τ h kb = 1, and the product between the conditional distributions defined as specific multinomial distributions. So, conditionally on x {kb}1h = 1, x {kb}j is deterministic for each j = 2, . . . , d {kb} . Indeed, in such a case, x {kb}j follows a specific multinomial distribution whose parameters are 0 and 1. More precisely, this distribution is parametrized by

δ hj kb = (δ hjh ′ kb ; h ′ = 1, . . . , m {kb} j
) with the following constraints defining the successive surjections:

δ hjh ′ kb ∈ {0, 1}, m {kb} j h ′ =1 δ hjh ′ kb = 1 (multinomial distribution) and m {kb} 1 h=1 δ hjh ′ kb ≥ 1 (surjections).
By denoting δ kb = (δ hj kb ; h = 1, . . . , m {kb} 1

; j = 2, . . . , d {kb} ), the pdf of maximum dependency distribution is defined as:

ṕ(x {kb} ; τ kb , δ kb ) = p(x {kb}1 ; τ kb ) d {kb} j=2 p(x {kb}j |x {kb}1 ; {δ hj kb } h=1,...,m {kb} 1 ) = m {kb} 1 h=1 τ h kb d {kb} j=2 m {kb} j h ′ =1 (δ hjh ′ kb ) x {kb}jh ′ x {kb}1h . (4)
Figure 2 shows two examples of the maximum dependency distributions. The probabilities of the joint distribution are represented by the area of dark boxes. Notice that δ kb defines the position where the probabilities are non-zero (location of a dark boxes) and τ kb defines the probabilities of this non-zero cells (area of the dark boxes). For the example illustrated in Figure 2.a, the probability that the first variable takes the modality one is 0.1 (τ 1 11 = 0.1). Moreover, conditionally on x {11}11 = 1, variable x {11}2 is deterministic since the probability that x {11}21 = 1 given x {11}11 = 1 is one. x {11}12

x {11}13 x {11}14

x {11}23

x {11}22

x {11}21 . This distribution has very limited interest because it is so unrealistic that it can almost never be used alone. We will see in the next section how to use it in a more efficient way.

A new block distribution: mixture of two extreme distributions

We proposed to model the distribution of each block by a bi-components mixture between an independence distribution and a maximum dependency distribution. For block b of component k, the block distribution is modeled by:

p(x {kb} ; θ kb ) = (1 -ρ kb )p(x {kb} ; α kb ) + ρ kb ṕ(x {kb} ; τ kb , δ kb ), (5)
where θ kb = (ρ kb , α kb , τ kb , δ kb ) and where ρ kb is the proportion of the maximum dependency distribution in this mixture with 0 ≤ ρ kb ≤ 1. The proposed model requires few additional parameters compared with the conditional independence model. In addition, it is easily interpretable as explained in the next paragraph. Note that the limiting case where ρ kb = 0 considers that the block follows an independence distribution. In this particular case, the parameters of the maximum dependency distribution are no longer defined.

Under this distribution, the proportion of the maximum dependency distribution reflects the deviation from independence under the assumption that the other allowed distribution is the maximum dependency distribution. The parameter ρ kb gives an indicator of the inter-variables correlation of the block. It is not here a pairwise dependency among variables but a dependency between all variables of the block. Furthermore, it stays bounded when the number of variables is larger than two while the Cramer's V is non-upper-bounded in this case. The intra-variables dependencies between the variables are defined by δ kb . The strength of these dependencies is explained by τ kb since it gives the weight of the over-represented modalities crossing compared with the independence distribution.

Above, we interpreted the distribution by conditionally independent blocks as a parsimonious version of the log-linear mixture model because it determines the interactions to be modeled for each class. By choosing the proposed distribution for blocks, a second level of parsimony is added.

Indeed, among the interactions allowed by this distribution with independent blocks, only those corresponding to the maximum dependency distribution will be modeled. Other interactions are considered as null.

Properties:

• The ccm, stays parsimonious compared with cim since, for each block with at least two variables, the number of the additional parameters depends only on the number of modalities of the first variable of the block and not on the number of variables in the block. By using ν cim defined in Equation ( 2), the number of parameters of ccm is denoted ν ccm by:

ν ccm = ν cim + {(k,b)|d {kb} >1} m {kb} 1 . (6)
• The proposed distribution is identifiable under the condition that the block is composed by = 2 then a new constraint is added. In order to have the most meaningful parameters, the chosen value of ρ kb is the largest value maximizing the log-likelihood. This additional constraint does not falsify the definition of ρ kb as an indicator of the dependency strength between the variables of the same block. Furthermore, this constraint is natural since blocks with the biggest dependencies are wanted. Note that ρ kb seems to be correlated with the Cramer's V. An example is given in Section 3 of [START_REF] Marbac | Model-based clustering for conditionally correlated categorical data[END_REF].

• Note that the marginal probabilities given class can be straightforwardly deduced from the parameters of the model.

Estimation of the parameters

For a fixed model (g, σ), the parameters must be estimated. Since the proposed distribution ccm has two latent variables (the class membership and the intra-block distribution membership), two algorithms derived from the EM algorithm perform the estimation of the associated continuous parameters. The combinatorial problems arising from the consideration of the discrete parameters are avoided by using a Metropolis-Hastings algorithm.

Global GEM algorithm

The whole data set consisting of n independent and identically distributed individuals is denoted by x = (x 1 , . . . x n ) where x i ∈ X . The objective is to obtain the maximum log-likelihood estimator θ defined as (g is now implicit in each expression)

θ = argmax θ L(θ; x, σ) with L(θ; x, σ) = n i=1 ln p(x i ; σ, θ) . (7)
The search for maximum likelihood estimates for mixture models entails solving equations having no analytical solutions. For the mixture models, the assignments of the individuals into the classes can be considered as missing data. So, the tool generally used is the Expectation-Maximization algorithm (denoted EM algorithm) and its extensions [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm[END_REF]. We denote the unknown indicator vectors of the g clusters by z = (z ik ; i = 1, . . . , n; k = 1, . . . , g) where z ik = 1 if x i arises from cluster k, z ik = 0 otherwise.

Thus, the mixture model distribution corresponds to the marginal distribution of the random variable X obtained from the couple distribution of the random variables (X, Z). In order to maximize the log-likelihood, the EM algorithm uses the complete-data log-likelihood which is defined as

L c (θ; x, z, σ) = n i=1 g k=1 z ik ln π k p(x i ; σ k , θ k ) . (8)
The EM algorithm is an iterative algorithm which alternates between two steps: the computation of the complete-data log-likelihood conditional expectation (E step) and its maximization (M step). Many algorithms are derived from the EM algorithm and among them the Generalized EM algorithm (GEM) is of interest to us. It works on the same principle as the EM algorithm, but the maximization step is replaced by a GM step where the proposed parameters increase the expectation of the complete-data log-likelihood according to its previous value without necessarily maximizing it.

We prefer to use the GEM algorithm, since the maximization step in the EM algorithm requires estimating the continuous parameters for too many possible values of the discrete parameters in order to warrant the maximization of the complete-data log-likelihood expectation. Indeed, exhaustive enumeration for estimating the discrete parameters is generally impossible when a block contains variables with many modalities and/or many variables, as detailed below. If S(a, b)

is the number of possible surjections from a set of cardinal a into a set of cardinal b, then δ kb is defined in the discrete space of dimension d {kb} -1 j=1 S(m {kb} j

, m {kb} j+1 ). For example, a block with three variables and m {kb} = (5, 4, 3) implies 51, 840 possibilities for δ kb . Thus, a stochastic approach is proposed in Section 4.2 to overcome this problem. Then, the estimation of the continuous parameters conditionally on the discrete parameters is performed via the classical EM algorithm presented in Section 4.3 since their estimation cannot be obtained in closed form. At iteration (r), the steps of the global GEM can be written as:

• E global step: z (r) ik = π (r) k p(x i ;σ k ,θ (r) k ) g k ′ =1 π (r) k ′ p(x i ;σ k ′ ,θ (r) k ′ ) , • GM global step: π (r+1) k = n (r) k n where n (r) k = n i=1 z (r) ik and ∀(k, b) θ (r+1) kb
is updated under the constraint that the conditional expectation of complete-data log-likelihood increases (see Sections 4.2 and 4.3).

Initialization of the algorithm: Since this algorithm is performed in an stochastic algorithm used for the model selection (see Section 5) and since this latter has an influence on the GEM initialization, this point will be detailed in Section 5.2.

Stopping criterion:

The GEM algorithm is stopped after r max iterations and we fix θ = θ (rmax) .

Details of the GM global step of the GEM

The maximization of the expected complete-data log-likelihood is done by optimizing its terms for each (k, b). Thus, the determination of θ Then the continuous parameters (ρ

(r,s+ 1 2 ) kb , α (r,s+ 1 2 ) kb , τ (r,s+ 1 2 ) kb
) are computed, conditionally on the value of δ (r,s+ 1 2 ) kb , in order to maximize the expected complete-data log-likelihood of block b for component k:

n i=1 z (r) ik ln p(x {kb} i ; ρ kb , α kb , τ kb , δ (r,s+ 1 2 ) kb
).

(10)

The candidate parameters are now denoted by θ

(r,s+ 1 2 ) kb = (ρ (r,s+ 1 2 ) kb , α (r,s+ 1 2 ) kb , τ (r,s+ 1 2 ) kb , δ (r,s+ 1 2 ) kb
). The whole block parameters θ (r,s+1) kb of the next step are then defined as θ (r,s+ 1 2 ) kb with the acceptance probability µ (r,s+1) and θ (r,s) kb otherwise, where: ; θ (r,s) kb ). Thus, the proposed initialization ensures an increased likelihood at each iteration of the GEM algorithm.

µ (r,s+1) = min n i=1 p(x {kb} i ; θ (r,s+ 1 2 ) kb ) z (r) ik |∆(δ (r,s+ 1 2 ) kb )| n i=1 p(x {kb} i ; θ (r,s) kb ) z (r) ik |∆(δ (r,s) kb )| , 1 , (11)
x {kb}11 x {kb}12 x {kb}13

x {kb}22

x {kb}21 (a)

x {kb}11 x {kb}12 x {kb}13

x {kb}22

x {kb}21 (b.1) Remark: When the space of possible δ kb is small (for example when the block groups a small number of binary variables), an exhaustive approach obtains the same results as the proposed algorithm with less computation time. Thus, the retained approach (exhaustive or stochastic) depends on the number of variables and modalities.

x {kb}11 x {kb}12 x {kb}13 x {kb}22 x {kb}21 (b.2) x {kb}11 x {kb}12 x {kb}13 x {kb}22 x {kb}21 (b.3) x {kb}11 x {kb}12 x {kb}13

Details of the M θ step of the GM global step

As there is a second level of mixing, another EM algorithm can be performed for the contin-

L full c (θ; x, y, z, σ) = n i=1 g k=1 z ik ln π k + B k b=1 (1 -y {kb} i ) ln(1 -ρ kb ) + (1 -y {kb} i ) ln p(x {kb} i ; α kb ) + y {kb} i ln ρ kb + y {kb} i ln ṕ(x {kb} i ; τ kb , δ kb ) . ( 12
)
At iteration (t), the local EM algorithm estimates the continuous parameters of block b, with fixed values of z (r) and δ (r,s+ 1 2 ) kb , by the following two steps:

• E local step: y {kb}(r,s+ 1 2 ,t) i = ρ (r,s+ 1 2 ,t) kb ṕ(x {kb} i ;τ (r,s+ 1 2 ,t) kb ,δ (r,s+ 1 2 ) kb ) p(x {kb} i ;ρ (r,s+ 1 2 ,t) kb ,α (r,s+ 1 2 ,t) kb ,τ (r,s+ 1 2 ,t) kb ,δ (r,s+ 1 2 ) kb ) , • M local step: ρ (r,s+ 1 2 ,t+1) kb = n (r,s+ 1 2 ,t) kb n (r) k , τ (r,s+ 1 2 ,t+1) kb = n i=1 z (r) ik y {kb}(r,s+ 1 2 ,t) i x {kb}1h i n (r,s+ 1 2 ,t) kb , α (r,s+ 1 2 ,t+1) kb = n i=1 z (r) ik (1-y {kb}(r,s+ 1 2 ,t) i )x {kb}jh i n (r) k -n (r,s+ 1 2 ,t) kb
, where n (r,s+ 1 2 ,t) kb

= n i=1 z (r)
ik y {kb}(r,s+ 1 2 ,t) i .

Conjecture: During our numerous experiments, we empirically noticed that the log-likelihood function of the mixture between the independence and the maximum dependency distributions had a unique optimum. We conjecture that this function has indeed a unique maximum.

Initialization of the algorithm: The previous conjecture allows to perform only one initialization of the EM algorithm fixed to: (ρ

(r,s+ 1 2 ,0) kb , α (r,s+ 1 2 ,0) kb , τ (r,s+ 1 2 ,0) kb ) = (ρ (r,s) kb , α (r,s) kb , τ (r,s) kb ).
Stopping criterion: This algorithm is stopped after a number of iterations denoted by t max and returns the value of block parameters θ

(r,s+ 1 2 ) kb defined as θ (r,s+ 1 2 ) kb = θ (r,s+ 1 2 ,tmax) kb .
Remark: In the specific case where δ kb are known for each (k, b), the estimation of all the continuous parameters could be performed by a unique EM algorithm where, at iteration (r), the E step would compute both z (r) and y (r) while the M step would estimate all the parameters maximizing the expectation of the full complete-data log-likelihood.

Model selection

Gibbs algorithm for exploring the space of models

Since the number of components g determines the dimension of σ, the model construction is done in two steps. Firstly, the selection of the number of components and, secondly, the determination of the variable repartition per blocks for each component. In a Bayesian context, the best model (ĝ, σ) is defined as [START_REF] Robert | Le choix bayesien: principes et pratique[END_REF]:

(ĝ, σ) = argmax g,σ p(g, σ|x). ( 13
)
Thus, by considering that p(g) = 1 gmax if g ≤ g max and 0 otherwise, where g max is the maximum number of classes allowed by the user, and by assuming that p(σ|g) follows a uniform distribution, the best model is also defined as:

(ĝ, σ) = argmax g argmax σ p(x|g, σ) . ( 14
)
To obtain (ĝ, σ), a Gibbs algorithm is used for estimating argmax σ p(x|g, σ), for each value of g ∈ {1, . . . , g max }. Indeed, this method limits the combinatorial problem involved by the detection of the block structure of variables. A reversible jump method could be used [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components (with discussion)[END_REF], but this approach is rarely performed with mixed parameters (continuous and discrete). Indeed, in such a case, it is difficult to define a mapping between the parameters space of two models. So, we propose to use an easier Gibbs sampler-type having p(σ|x, g) as stationary distribution. It alternates between two steps: the generation of a stochastic neighborhood Σ [q] conditionally on the current model σ [q] by a proposal distribution and the generation of a new pattern σ [q+1] included in Σ [q] with a probability proportional to its posterior probability. At iteration [q], it is written as:

• Neighborhood step: generate a stochastic neighborhood Σ [q] by a proposal distribution given below conditionally on the current model σ [q] ,

• Pattern step: σ [q+1] ∼ p(σ|x, g, Σ [q] ) with p(σ|x, g,

Σ [q] ) =    p(x|g,σ) σ ′ ∈Σ [q] p(x|g,σ ′ ) if σ ∈ Σ [q]
0 otherwise.

A possible deterministic neighborhood of σ [q] could be defined as the set of models where, at most one variable is affected, for one component, in another block (possibly creating a new block):

σ : ∃!(k, b, j) j ∈ σ [q]
kb and j / ∈ σ kb ∪ σ [q] . However, as this deterministic neighborhood can be very large, our proposal distribution allows reducing it to a stochastic neighborhood Σ [q] by reducing the number of (k, b) where σ kb could be different to σ

[q]

kb . Thus, one component k [q] is randomly sampled in {1, . . . , g} then one block b

[q] f rom is randomly sampled in {1, . . . , B [q] k [q] }. Another block b [q] is randomly sampled in {1, . . . , B [q] k [q] } \ b [q]
f rom and the set b

[q] to = {b [q] , B [q] k [q] + 1}
is built. The stochastic neighborhood Σ [q] is then defined as:

Σ [q] = σ : ∃!(k, b, j) j ∈ σ [q] kb , j / ∈ σ kb and j ∈ σ kb ′ with k = k [q] , b = b [q] from , b ′ ∈ b [q]
to ∪ σ [q] . (15)

We denote the elements of Σ [q] as σ [q+ε(e)] where ε(e) = e |Σ [q] |+1 and e = 1, . . . , |Σ [q] |. Figure 4 shows an illustration of this definition.

(c)

Figure 4: Example of the support of Σ [q] in a case of four variables. If the variable of row j and the variable of column j ′ are in the same block then cell (j, j ′ ) is painted in black. This cell is painted in white otherwise. (a) Graphical representation of σ

[q] k = ({1, 2}, {3, 4}); (b) Elements of Σ [q] if b [q] f rom = 1; (c) Elements of Σ [q] if b [q] f rom = 2.
At the generation pattern step, the previous algorithm needs the value of p(x|g, σ) ∀σ ∈ Σ [q] . By using the BIC approximation [START_REF] Schwarz | Estimating the dimension of a model[END_REF][START_REF] Lebarbier | Une introduction au critre BIC : fondements thoriques et interprtation[END_REF], this probability is approximated by:

ln p(x|g, σ) ≃ L( θ; x, g, σ) - ν ccm 2 log(n), (16)
θ being the maximum likelihood estimator obtained by the GEM algorithm previously described in Section 4. Thus, at the iteration [q], for each e = 1, . . . , |Σ [q] |, estimator θ[q+ε(e)] associated to element σ [q+ε(e)] is computed by the GEM algorithm.

Initialization: Whatever the initial value selected for σ [0] , the algorithm converges to the same value of σ. However, this convergence can be very slow when the initialization is poor. Since blocks consist of the most correlated variables, a Hierarchical Ascendant Classification (HAC) is used on the matrix of Cramer's V distances on the variable couples. The partition produced by the HAC minimizing the block number without blocks consisting of more than four variables is chosen for each σ

[0]

k . For the initialization, the variables number of a block is limited to four, because very few blocks having more than four variables were exhibited in the course of our experiments.

Obviously, the Gibbs algorithm can then violate this initial constraint if necessary.

Stopping criterion:

The algorithm is stopped when q max successive iterations have not discovered a better model.

Consequences of the Gibbs algorithm on the GEM algorithm

Initialization of the GEM algorithm: At iteration [q] of the Gibbs algorithm, the GEM algorithm estimates θ[q+ε(e)] associated to model σ [q+ε(e)] for e = 1, . . . , |Σ [q] |. Since these models are close to σ [q] , their maximum likelihood estimators should be closed to θ [q] . The GEM algorithm initialization is also done by the value of θ[q] for the not modified blocks. Thus, θ

[q+ε(e)](0) kb = θ[q]
kb if the blocks are not modified (σ

[q+ε(e)] kb = σ [q]
kb ). For the other blocks, the continuous parameters are randomly sampled. For those blocks, in order to avoid the combinatorial problems, we use a sequential method to initialize δ [q+ε(e)](0) kb : the surjections from x {kb}1 to x {kb}j are sampled, according to x and to the continuous parameters previously sampled (ρ

[q+ε(e)](0) kb , α [q+ε(e)](0) kb , τ [q+ε(e)](0) kb ),
for each j = 2, . . . , d {kb} as follows:

δ .j[q+ε(e)](0) kb ∝ n i=1 p(x {kb}1 i , x {kb}j i ; ρ [q+ε(e)](0) kb , α 1[q+ε(e)](0) kb , α j[q+ε(e)](0) kb , τ [q+ε(e)](0) kb , δ .j kb ) z [q] ik , (17) where δ
.j[q+ε(e)] kb

= (δ

hj[q+ε(e)] kb

; h = 1, . . . , m {kb} 1
) and where z

[q] ik = E Z ik |x i , θ [q] .
Remark about r max : As said in Section 4.1, the algorithm is stopped after a fixed number of iterations r max . If the algorithm is stopped before its convergence, the proposed initialization limits the problems. Indeed, if the model has a high a posteriori probability, it will stay in neighborhood Σ [q] during some successive iterations, so its log-likelihood will increase.

Simulations

Table 1 presents the adjustment parameters values used for all the simulations.

Algorithms Gibbs GEM Metropolis-Hastings EM Criteria q max = 20 × d r max = 10 s max = 1 t max = 5
Table 1: Values of the different stopping criteria. for controlling the overlapping of classes: when it is close to one their overlapping (Bayes error) is close to one. These parameters set the error rate to 0.10 for each studied situation: (0.60, 0.20, 0.20), (0.20, 0.20, 0.60),

σ kb = (d/b, 1 + d/b) ρ kb = 0.6(1 -u) τ kb =
δ h2h ′ 1b = 1 iff h = h ′ δ 122 1b = δ 223 1b = δ 321 1b = 1 α j 1b =
α 1 2b = α 1 1b (1 -u) + (0.

Application

Dentistry clustering

The Handelman's dentistry data [START_REF] Handelman | Clinical radiographic evaluation of sealed carious and sound tooth surfaces[END_REF] display the evaluation of 3869 dental

x-rays (sound or carious) performed by five dentist and possibly showing incipient caries. This data set has been clustered by several models in the past. It is suggested that there are two main classes: the sound teeth and the carious ones.

According to the BIC criterion, data are split into three classes by cim. Furthermore, dependencies are observed between the variables into classes since the Cramer's V computed per class is not close to zero. Thus, [START_REF] Espeland | Using Latent Class Models to Characterize and Assess Relative Error in Discrete Measurements[END_REF] apply a log-linear mixture model to fit the data. The authors set the model, by adding some assumptions to better fit the data. More precisely, they consider a mixture with four components. The first two components take into account the interactions between dentists 3 and 4. The last two components are specific since they allow only one modality interaction, when all the diagnoses are respectively carious and sound. Note that these assumptions are required by the above authors due to their realistic nature. Indeed, this model fits the data better than cim. On the other hand, its interpretation needs the analysis of four classes.

As the last two classes seem artificial, [START_REF] Qu | Random Effects Models in Latent Class Analysis for Evaluating Accuracy of Diagnostic Tests[END_REF] prefer to use the random effects models in a latent class analysis with two classes. They assume that the conditional dependencies can be modeled by a single continuous latent variable varying among the individuals. According to the authors, one class represents the sound teeth and the other represents the carious ones, while the random effect represents all the patient-specific unrecorded characteristics of the x-ray images.

Their model does not need the two additional artificial classes. Thus their interpretation is easier.

We now display the results of the proposed model ccm estimated with the R package Clustericat (the code is presented in Appendix A). The BIC criterion selects two classes with a value of -7473.

It claims that cmm better fits the data than the model of [START_REF] Qu | Random Effects Models in Latent Class Analysis for Evaluating Accuracy of Diagnostic Tests[END_REF] since their BIC criterion value is -7487. The BIC criterion values for cim and cmm are displayed in Table 3. We indicate the computing time (in seconds), obtained with an Intel Core i5-3320M processor, to estimate cmm where 20 MCMC chains were started with a stopping rule q max = 100 while cim needs less than 0.1 sec with 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]. Table 3: BIC criterion values for the cim and the cmm according to different numbers of classes for the dentistry data set. For each model, the best results according to the BIC criterion are in bold. The computing time in seconds is indicated for cmm where 20 MCMC chains were started with a stopping rule q max = 500.

We note that cmm obtains better values for the BIC criterion than cim when g = 1, 2. When the number of classes is larger (g ≥ 3) the best model of cmm assumes the conditional independence between variables.

The BIC criterion selects two classes for cmm and this is coherent with a clustering of the teeth between the sound and the carious ones. Furthermore, the two main characteristics of the model fixed by [START_REF] Espeland | Using Latent Class Models to Characterize and Assess Relative Error in Discrete Measurements[END_REF] are automatically detected by the model: importance of the two modality crossings where all the dentists have the same diagnosis and a dependency between the diagnosis of dentists 3 and 4. Thus, the estimated model is coherent with the imposed model of [START_REF] Espeland | Using Latent Class Models to Characterize and Assess Relative Error in Discrete Measurements[END_REF] while no information was given a priori.

The fitted model can be interpreted as:

• the majority class (π 1 = 0.86) mainly gathers the sound teeth. There is a strong dependency between the five dentists (σ 1 = ({1, 2, 3, 4, 5}) and ρ 11 = 0.35). The dependency structure of the maximum dependency distribution indicates an over-contribution of both modality interactions where the five dentists have the same diagnosis, especially when they claim that the teeth are sound (τ all sound 11 = 0.93 and τ all carious 11 = 0.07). This class could be interpreted as grouping the teeth who these diagnosis is obvious.

• the minority class (π 2 = 0.14) groups principally the carious teeth. There is a dependency between dentists 3 and 4 while the diagnosis of the other ones are independent given the class (σ 2 = ({3, 4}, {1, 2, 5}), ρ 21 = 0.31 and ρ 22 = 0). This class could be interpreted as grouping the teeth whose diagnosis is more complex.

Figure 6 helps the interpretation of the clusters for the ccm with two components (best model according to the BIC criterion). On ordinates, the estimated classes are represented with respect to their proportions in decreasing order. Their corresponding area depends on their proportion.

The cumulated proportions are indicated on the left side. On abscissae, three indications are given.

The first one is the inter-variables correlations (ρ kb ) for all the blocks of the class ordered by their strength of correlation (in decreasing order). The second one is the intra-variables correlations (τ kb ) for each block drawn according to their strength dependencies (in decreasing order). The third is the variables repartition per block. A black cell indicates that the variable is assigned to the block and a white cell indicates that, conditionally on this class, the variable is independent of the variables of this block. For example, this figure shows that the first class has a proportion of 0.86 and that all the variables are assigned into the same block.

Calves clustering

In this section, the results obtained by the ccm are compared to those obtained for the cim by the RMixmod software [START_REF] Lebret | Rmixmod: The R Package of the Model-Based Unsupervised, Supervised and Semi-Supervised Classification Mixmod Library[END_REF]. The "Genes Diffusion" company has collected information from French breeders in order to cluster calves. The 4270 studied calves are described by nine variables related to behavior (aptitude for sucking Apt, behavior of the mother just before The computing time for the ccm estimation is given in minutes.

For the cim, the BIC criterion selects a high number of classes, since it selected eight classes.

The interpretation of the clusters is also difficult and we can assume that the estimator's quality is very poor. Figure 7 helps the interpretation for the ccm with five components (best model according to the BIC criterion). Its interpretation is the same as the interpretation of Figure 6.

For example, this figure shows that the first class has a proportion of 0.29 and that it is composed of four blocks. The most correlated block of the first class has ρ kb ≃ 0.80 and the strength of the biggest modalities link is also close to 0.85. This block consists of the variables TDC and TRM.

Here is now a possible interpretation of Class 1 (note that the others classes are also meaningful; • General: This class has a proportion equal to 0.29 and consists of three blocks of dependency and one block of independence.

A p t T O C T R C T D C I s o D i s E m p T R M T D M A p t T O C T R C T D C I s o D i s E m p T R M T D M A p t T O C T R C T D C I s o D i s E m p T R M T D M A p t T O C T R C T D C I s o D i s E m p T R M T D M Class 2 A p t T O C T R C T D C I s o D i s E m p T R M T D M A p t T O C T R C T D C I s o D i s E m p T R M T D M A p t T O C T R
• Block 1: There is a strong correlation (ρ 11 ) between the variables diarrhea treatment of the calf and mother preventive treatment against respiratory disease, especially between the modality no treatment against the calf diarrhea and the absence of preventive treatment against respiratory disease of its mother (τ 11 and δ 11 ).

• Block 2: There is a strong correlation (ρ 12 ) between the variables treatment against respiratory illness of the calf and mother preventive treatment against diarrhea, especially between the modality preventive treatment against respiratory illness of the calf and the presence of diarrhea preventive treatment of its mother (τ 12 and δ 12 ).

• Block 3: There exists another strong link between the behavior of mother, the emptying of the umbilical and its disinfection (τ 13 and δ 13 ).

• Block 4: This block is characterized by absence of preventive treatment against omphalitis and having 50% of the calves infected by this illness (α 14 ).

Conclusion

By using the block extension of the latent class model, a new mixture model is proposed for clustering categorical data by taking into account the intra-class correlation. The block distribution is defined as a mixture between an independent distribution and a maximum dependency distribution. This specific distribution, which remains parsimonious, is compared to the full latent class model and allows different levels of interpretation. The blocks of variables detect the conditional dependencies between variables while their strengths are reflected by the proportions of maximum dependency distribution. The parameters of the block distribution reflect the links and the strength between modalities.

The parameter estimation and the model selection are simultaneously performed via a Gibbs sample-type algorithm. It allows reducing the combinatorial problems of the block structure detection and the links between modalities search for the estimation of the maximum dependency distribution. The results are good when the number of modalities is small for each variable. For more than six modalities, the detection of other links encounters some persistent difficulties. So the algorithm can be slow in this case. The proposed approach to estimate the block structure is not adapted for data sets with many variables. A deterministic but sub-optimal solution could be used to perform a forward algorithm.

The R package Clustericat allows clustering categorical data sets by using cmm. This package is available on Rforge at the following url https://r-forge.r-project.org/R/?group id=1803.

The proposed model can be easily extended to the case of ordinal data. For this purpose, additional constraints on the dependency structure of each distribution of maximum dependency need to be added. A Dentistry clustering with the R package Clustericat

The R package Clustericat is available on Rforge website at the following url: https://r-forge.rproject.org/R/?group id=1803. This section presents the code used to cluster the dentistry data # a plot summarizing the best model like Figure 6 > plot(res)

  It considers that conditionally on the class k, variables are grouped into b k independent blocks and that each block follows a specific distribution. The repartition in blocks of the variables determines a partition σ k = (σ k1 , . . . , σ kb k ) of {1, . . . , d} in b k disjoint non-empty subsets where σ kb represents the subset b of variables in the partition σ k . This partition defines x {kb} = x σ kb = (x {kb}j ; j = 1, . . . , d {kb} ) which is the subset of x associated to σ kb . The integer d {kb} = card(σ kb ) is the number of variables in block b of component k and x {kb}j = (x {kb}jh ; h = 1, . . . , m {kb} j ) corresponds to variable j of block b for component k with x {kb}jh = 1 if the individual takes modality h for variable x {kb}j and x {kb}jh = 0 otherwise and where m {kb} j represents the number of modalities of x {kb}j . Note that different repartitions of the variables into blocks are allowed for each component and they are grouped into σ = (σ 1 , . . . , σ g ).

  presents an example of the distribution with conditional independent blocks for a mixture with two components described by five variables. Blank cells indicate that the intra-class correlation is neglected and black cells indicate that this correlation is taken into account. For instance, Figure 1.a illustrates the distribution of the class 1 where two blocks (b 1 = 2) are considered. More precisely, the first block is composed by the first two variables (σ 11 = {1, 2}) and the second block is composed by the last three variables (σ 12 = {3, 4, 5}). Note that the classical latent class model with conditional

Figure 1 :

 1 Figure 1: Example of the conditional independent blocks mixture model with g = 2 and d = 5, (a) k = 1, B 1 = 2 and σ 1 = ({1, 2}, {3, 4, 5}), (b) k = 2, B 2 = 3 and σ 2 = ({1, 5}, {2, 4}, {3}).

  this context, we propose to model the distribution of each block by a mixture of the extreme distributions according to the Cramers V criterion computed on all the variable couples. The model results in a bi-component mixture between an independence distribution and a maximum dependency distribution which can be easily interpreted by the user. The maximum dependency distribution is introduced first. The resulting conditional correlated model (CCM) is also defined as a block model extension of the latent class model where the distribution inside the block is modeled by this bi-component mixture. Remark: Without loss of generality, the variables are considered as ordered by decreasing number of modalities in each block: (k, b) mkb j mkb j The goal is now to define a parsimonious distribution for each block that takes into account the dependency between variables. Furthermore, the parameters of the distribution inside the block must be meaningful for the practitioner. In this context, we propose to model the distribution of each block by a mixture of the extreme distributions according to the Cramers V criterion computed on all the variable couples. The model results in a bi-component mixture between an independence distribution and a maximum dependency distribution which can be easily interpreted by the user. The maximum dependency distribution is introduced first. The resulting conditional correlated model (ccm) is also defined as a block model extension of the latent class model where the distribution inside the block is modeled by this bi-component mixture. Remark: Without loss of generality, the variables are considered as ordered by decreasing number of modalities in each block: ∀(k, b) m {kb} j ≥ m {kb} j+1 where j = 1, . . . , d {kb} -1.

  variable determines the other ones, this distribution is defined by a product between the multinomial distribution of the first variable parametrized by τ kb = (τ h kb ;

Figure 2 :

 2 Figure 2: Two examples of the maximum dependency distributions for the first component of the mixture illustrated by Figure 1(a). (a) The first block is displayed with m {11}1 = 4, m {11}2 = 3, δ h1h 11 = 1 for h = 1, 2, 3, δ 413 11 = 1 and τ 11 = (0.1, 0.3, 0.2, 0.4); (b) The second block is displayed with m {12}1 = m {12}2 = m {12}3 = 2, δ hjh ′ 12 = 1 iff (h = h ′ ) and τ 12 = (0.5, 0.5).

  at least three variables (d {kb} > 2) or that the modality number of the last variable of the block is more than two (m {kb} 2 > 2). This result is demonstrated in Appendix B of Marbac et al. (2013). The parameter ρ kb is a new indicator allowing measuring the correlation between variables, not limited to correlation between variable couples. In the case where the identifiability conditions cannot be met, we distinguish two cases. If d {kb} = 1, then block b contains only one variable, and the proposed model is reduced to model a multinomial distribution, ρ kb = 0 and the maximum dependency distribution is not defined. If d {kb} = 2 and m {kb} 2

  of the parameters of the other blocks. A Metropolis-Hastings algorithm[START_REF] Robert | Monte Carlo statistical methods[END_REF] is also performed, for each (k, b), to avoid the combinatorial problems induced by the detection of the discrete parameters δ kb . It performs a random walk over the discrete parameters space and computes the maximum likelihood estimators of the continuous parameters (ρ kb , α kb , τ kb ) associated with them. This stochastic algorithm allows finding the estimator maximizing the expected complete-data log-likelihood of block b for component k: each iteration (s) of this Metropolis-Hastings algorithm, a discrete parameter denoted by δ (r,s+ 1 2 ) kb is sampled with a uniform distribution in a neighborhood of δ (r,s) kb denoted as ∆(δ (r,s) kb ).

••

  kb )| denoting the cardinal of ∆(δ (r,s) kb ). Thus, at iteration (s), the algorithm performs the three following steps: • Stochastic step on δ kb : generate δ Maximization step on the continuous parameters (M θ step): compute the continu-Stochastic step on θ kb : sample θ kb ) is defined as the set of parameters where at most two surjections are different from that of δ (r,s) kb . Figure 3 illustrates this definition.Initialization of the algorithm: The initialization of the algorithm is done by θ This algorithm is stopped after a number of iterations s max . The parameter θ

Figure 3 :

 3 Figure 3: Example of ∆(δ kb ) with d {kb} = 2 and m {kb} = (3, 2). For row h ′ and column h, a black cell indicates that δ h2h ′ kb = 1 and a white cell that δ h2h ′ kb = 0: (a) δ kb ; (b.1), (b.2), (b.3), (b.4) are

  uous parameters (ρ kb , α kb , τ kb ) estimation by introducing other unknown vectors corresponding to the indicator of the blocks distributions conditionally on z. These vectors are denoted by y = (y {kb} ; k = 1, . . . , g; b = 1, . . . , b k ) with y {kb} = (y maximum dependency distribution for block b of cluster k and y {kb} i = 0 if x {kb} i arises from the independence distribution for block b of cluster k. The whole mixture model distribution corresponds to the marginal distribution of the random variable X obtained from the triplet distribution of the random variables (X, Y, Z). Since the blocks are independent conditionally on Z, the full complete-data log-likelihood (both in Y and Z) is defined as:

Figure 5 :

 5 Figure 5: Box-plots of the number of iterations required by the Metropolis-Hastings algorithm for finding the best links between modalities, according to the number of modalities when datasets are simulated with a proportion of maximum dependency distribution equal to 0.5. (a) Three variables with the proposed initialization; (b) Three modalities per variable with the proposed initialization; (c) Three variables with a random initialization; (d) Three modalities per variable with a random initialization.

Figure 6 :

 6 Figure 6: Summary of the best ccm according to BIC for the dentists data set.

Figure 7 :

 7 Figure 7: Summary of the best ccm according to BIC for the calves data set.

#

  Loading of the data set > data("dentist") # to define the parameters of the algorithm performing the estimation # here 25 MCMC are performed with a stopping criterion equals to # 200 successive iterations having not found a better model > st <strategycat(dentist, nb init=25, stop criterion=200) # estimation of the model for a class number equal to 2. # for the data set with five binary variables (modal) > res <clustercat(dentist, 2, modal=rep(2,5), st) # presentation of the best model > summary(res) # presentation of the parameters of the conditional dependencies for the best model > summary dependencies(res)

Table 2 :

 2 Table2shows the mean and the standard deviation of the Kullback-Leibler divergence between the parameters used for the dataset generation and the estimated parameters according to the number of variables. When n increases, the Kullback-Leibler divergence converges to zero. It confirms the good behavior of the proposed algorithm. Mean (standard deviation) of the Kullback-Leibler divergence.

	d \ n	100	200	400	800
	4	0.77 (1.34 ) 0.26 (0.26 ) 0.15 (0.05 ) 0.12 (0.05 )
	6	1.22 (1.77 ) 0.27 (0.14 ) 0.09 (0.07 ) 0.05 (0.05 )
	8	1.72 (2.50 ) 0.41 (0.20 ) 0.09 (0.05 ) 0.05 (0.03 )
	10	1.73 (4.06 ) 0.52 (0.14 ) 0.10 (0.03 ) 0.04 (0.03 )

Table 4

 4

	g		1	2	3	4	5	6	7	8
	cim	BIC	-28589 -26859 -26526 -26333 -26238 -26235 -26226 -26185
		ν cim	17	35	53	71	89	107	125	143
	ccm	BIC	-26653 -26289 -26173 -26038 -26025 -26059 -26045	-26058
		ν ccm	24	48	80	89	112	131	148	163
	time (min)	0.97	3.32	6.16	6.56	10.03	11.76	12.31	14.92

displays the BIC criterion values and the number of parameters for the cim and ccm models. Furthermore, the computing time in minutes (obtained with an Intel Core i5-3320M processor) to estimate ccm by starting 20 MCMC chains with a stopping criterion of q max = 180 while cim needs 3 seconds with 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]

.

Table 4 :

 4 Results for the cim and the cmm according to different class numbers. For both models, first row corresponds to the BIC criterion values and the second row indicates the continuous parameter number. For each model, the best results according to the BIC criterion are in bold.

  Vermunt, J.K. 2003. Multilevel latent class models. Sociological methodology, 33(1), 213-239. Vermunt, J.K. 2007. Multilevel mixture item response theory models: an application in education testing. Proceedings of the 56th session of the International Statistical Institute. Lisbon,

	Portugal, 22-28.

The R package Clustericat is available on Rforge website at the following url: https://r-forge.rproject.org/R/?group id=1803

X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 (b) X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4 X 4 X 3 X 2 X 1 X 1 X 2 X 3 X 4

In fact, the algorithm is stopped as soon as it finds a discrete estimator involving a likelihood higher than or equal to the likelihood obtained with the true discrete parameters used for the simulation.

Acknowledgments: The author thank the Associate Editor and the three anonymous referee for their useful comments and references. The authors are grateful to the Genes Diffusion corporation for the provision of the data set and especially its members: Amélie Vallée, Julie Hamon and Claude Grenier. We are grateful to Parmeet Bhatia, Modal Team engineer, and Stéphane Chrétien for their precious assistance. This work was financed by DGA and Inria.

As these algorithms are interlocked, the iterations number of the most internal algorithms are small. Since the number of possible models increases with d, we propose to set: q max = 20 × d.

When the best model is selected by the Gibbs algorithm, this latter will stay in this model during many iterations so the Metropolis-Hastings and the EM algorithm are performed many times.

Thus, it is not necessary to have a large number of iterations as stopping criterion.

6.1 Study of the algorithm for the δ kb estimation

In this section, we illustrate the performance of the Metropolis-Hastings algorithm used for the δ kb estimation (see Section 4.2) and the relevance of its initialization (see Equation ( 17)). Since this algorithm is interlocked in the Gibbs and in the GEM algorithm, we need it to converge quickly.

The following simulations show that the algorithm stays relevant up to six modalities per variable and up to six variables per block. These conditions hold in most situations.

Samples of size 200 described by variables having the same number of modalities are generated by a mixture between an independence distribution and a maximum dependency distribution.

The parameter estimation is also performed by the Metropolis-Hastings algorithm, described in Section 4.2, since only one class is generated. The discrete parameters initializations are performed according to Equation ( 17) with z i1 = 1 for i = 1, . . . , 200.

Figure 5 shows the box-plots of the iterations number needed by the Metropolis-Hastings algorithm for finding the true links between modalities maximizing the likelihood 2 . According to these simulations, one observes that the results of this algorithm are good as a result of its initialization which allows significantly reducing the number of iterations needed in order to find the maximum likelihood estimators.

Study of the algorithm for model selection

In order to illustrate the efficiency of the algorithm for the model selection (and also the included estimation process), we want to study the evolution of the Kullback-Leibler divergence according to the number of variables and to the size of the data set. Thus, 100 samples are generated for many situations according to the ccm with two components. Note that the parameter u is introduced