Estimation and selection for the latent block model on categorical data

This paper deals with estimation and model selection in the Latent Block Model (LBM) for categorical data. First, after providing sufficient conditions ensuring the identifiability of this model, we generalise estimation procedures and model selection criteria derived for binary data. Secondly, we develop Bayesian inference through Gibbs sampling and with a well calibrated non informative prior distribution, in order to get the MAP estimator: this is proved to avoid the traps encountered by the LBM with the maximum likelihood methodology. Then model selection criteria are presented. In particular an exact expression of the integrated completed likelihood criterion requiring no asymptotic approximation is derived. Finally numerical experiments on both simulated and real data sets highlight the appeal of the proposed estimation and model selection procedures.

1. The data set is available from https://www.hds.utc.fr/coclustering/doku.php.

2. Congressional Voting Records data set is available from http://archive.ics.uci.edu/ml/datasets/Congressional+Voting+Records.

3. The code of Wyse and Friel is available at http://www.ucd.ie/statdept/jwyse.

The authors thank the reviewers and the Associate Editor for their very helpful comments which have greatly improved this paper. C. K. has been supported by LMH (Hadamard Mathematics Labex), backed by the foundation FMJH (Fondation Mathématique Jacques Hadamard).

Authors and Affiliations

1. Laboratoire de Mathématiques, Equipe Probabilités et Statistiques, UMR 8628, Bâtiment 425, Université Paris-Sud, 91405 , Orsay Cedex, France

   Christine Keribin & Vincent Brault

2. INRIA Saclay Île de France, Bâtiment 425, Université Paris-Sud,  91405 , Orsay Cedex, France

   Gilles Celeux

3. U.T.C, U.M.R. 7253, C.N.R.S. Heudiasyc, Centre de Recherches de Royallieu, B.P. 20529,  60205 , Compiègne Cedex, France

   Gérard Govaert

Corresponding author

Correspondence to Christine Keribin.

Appendices

Appendix A: Proof of Theorem 1

The identifiability of a parametric model requires that for any two different values \(\varvec{\theta }\ne \varvec{\theta }'\) in the parameter space, the corresponding probability distribution \({I\!P}_{\varvec{\theta }}\) and \({I\!P}_{\varvec{\theta }'}\) are different. We prove that, under assumptions of Theorem 1, there exists a unique parameter \(\varvec{\theta }=(\varvec{\pi },\varvec{\rho },\varvec{\alpha })\), up to a permutation of row and column labels, corresponding to \({I\!P}(\mathbf{y})\) the probability distribution function of matrix \(\mathbf{y}\) having at least \(2m-1\) rows and \(2g-1\) columns. The proof is adapted from Celisse et al. (2012) who set up a similar result for the Stochastic Block Model. Notice first that \(\varvec{\tau }= A \varvec{\rho }\) is the vector of the probability \(\tau _k\) to have \(1\) in a cell of a row of given row class \(k\):

As all the coordinates of \(\varvec{\tau }\) are distinct, \(R\) defined as the \(g\times g\) matrix of coefficients \(R_{ik}=(\tau _k)^i\), for \(0 \le i < g\) and \(1 \le k \le g\), is Vandermonde, and hence invertible. Consider now \(u_i\), the probability to have \(1\) on the \(i\) first cells of the first row of \(\mathbf{y}\):

With a given \({I\!P}(\mathbf{y}),\, u_1, \ldots ,u_{2g-1}\) are known, and we denote \(u_0=1\). Let now \(M\) be the \((g+1)\times g\) matrix defined by \(M_{ij}=u_{i+j-2}\), for all \(1\le i \le g+1\) and \(1\le j\le g\) and define \(M_i\) as the square matrix obtained by removing the row \(i\) from \(M\). The coefficients of \(M\) are

We can write, with \(A_{\pi }=\hbox {Diag}(\varvec{\pi })\)

Now, \(R\), unknown at this stage, can be retrieved by noticing that the coefficients of \(\varvec{\tau }\) are the roots of the following polynomial (Celisse et al. 2012) of degree \(g\)

where \(D_k=\hbox {det} M_k\) and \(D_g\ne 0\) as \(M_g\) is a product of invertible matrices. Hence, it is possible to determine \(\varvec{\tau }\), and \(R\) is now known. Consequently, \(\varvec{\pi }\) is defined in a unique manner by \(A_{\pi }=R^{-1} M_g R^{'-1}\).

In the same way, \(\varvec{\rho }\) is defined in a unique manner by considering the probabilities \(\sigma _\ell \) to have \(1\) in a column of column class \(\ell \) and the probabilities \(v_j\) to have \(1\) on the \(j\) first cells of the first column of \(\mathbf{y}\)

To determine \(\varvec{\alpha }\), we introduce for \(1\le i \le g\) and \(1\le j \le m\) the probabilities \(U_{ij}\) to have 1 in each \(i\) first cells of the first row and in each \(j\) first cells of the first column

These probabilities are known and we can write, with \(S_{j\ell }=(\sigma _\ell )^{j-1},\,j=1,\ldots m,\,\ell =1,\ldots ,m\)

and it leads to define \(A=A_{\pi }^{-1}R^{-1} U S^{'-1} A_{\rho }^{-1}\), and hence \(\varvec{\alpha }\), in a unique manner.

The proof is straightforwardly extended to the categorical case. The identification of \(\varvec{\pi }\) and \(\varvec{\rho }\) are obtained by considering the probabilities of \(y_{ij}=1\) where \(1\) is the first outcome of the multinomial distribution. Then, each \(A^h=(\alpha _{k\ell }^h)_{k=1,\ldots ,g;\ell =1,\ldots ,m}\) is successively identified by considering \(y_{ij}=\ell ,\,i=1,\ldots ,m,\,j=1,\ldots , g\).

Appendix B: VEM algorithm for categorical data

Govaert and Nadif (2008) proposed to use, for the binary case, a variational approximation of the EM algorithm by imposing that the joint distribution of the labels takes the form \(q^{(c)}_{zw}(\mathbf{z},\mathbf{w})=q^{(c)}_z(\mathbf{z})q^{(c)}_w(\mathbf{w})\). To get simpler formulas, we will denote

and

Using the variational approximation, the maximisation of the loglikelihood is replaced by the maximisation of the free energy

alternatively in \(q_{z},\,q_w\) and \(\varvec{\theta }\) (see Keribin 2010). The difference between the maximum loglikelihood and the maximum free energy is the Kullback divergence \(KL(q_{zw}||p(\mathbf{z},\mathbf{w}|\mathbf{y};\varvec{\theta }) )={I\!E}_{q_{zw}}\left[ \log \frac{q_{zw}(\mathbf{z},\mathbf{w})}{p(\mathbf{z},\mathbf{w}|\mathbf{y};\varvec{\theta })} \right] \). This can be extended to the categorical LBM, and leads to the following Variational EM (VEM) algorithm:

E step It consists of maximising the free energy in \(q_{z}\) and \(q_w\), and it leads to:

1. 1.

   Computing \(s_{ik}^{(c+1)}\) with fixed \(w^{(c)}_{jl}\) and \(\varvec{\theta }^{(c)}\)

   $$\begin{aligned} s_{ik}^{(c+1)}= \frac{\pi _k^{(c)} \psi _k(\mathbf{y}_{i\cdot };\varvec{\alpha }_{k\cdot }^{(c)})}{\sum _{k'} \pi _{k'}^{(c)} \psi _{k'}^{(c)}(\mathbf{y}_{i\cdot };\varvec{\alpha }_{k'\cdot }^{(c)}) }, k=1,\ldots , g \end{aligned}$$

   where \(\mathbf{y}_{i\cdot }\) denotes the row \(i\) of the matrix \(\mathbf{y}\), \(\varvec{\alpha }_{k\cdot }=(\alpha _{k1},\ldots ,\alpha _{km})\), and

   $$\begin{aligned} \psi _k(\mathbf{y}_{i\cdot };\varvec{\alpha }_{k\cdot }^{(c)})=\prod _{\ell ,h} ({\alpha _{k\ell }^{h}}^{(c)})^{\sum _j t_{j\ell }^{(c)} y_{ij}^{h}} \end{aligned}$$
2. 2.

   Computing \(t_{jl}^{(c+1)}\) with fixed \(s^{(c+1)}_{ik}\) and \(\varvec{\theta }^{(c)}\)

   $$\begin{aligned} t_{j\ell }^{(c+1)}= \frac{\rho _{\ell }^{(c)} \phi _{\ell }(\mathbf{y}_{\cdot j};\varvec{\alpha }_{\cdot \ell }^{(c)})}{\sum _{\ell '} \rho _{\ell '}^{(c)} \phi _{\ell '}(\mathbf{y}_{\cdot j};\varvec{\alpha }_{\cdot \ell '}^{(c)}) } , \ell =1,\ldots ,m \end{aligned}$$

   where \(\mathbf{y}_{\cdot j}\) denotes the column \(j\) of the matrix \(\mathbf{y},\, \varvec{\alpha }_{\cdot \ell }=(\alpha _{1\ell },\ldots ,\alpha _{g\ell })\) and

   $$\begin{aligned} \phi _{\ell }(\mathbf{y}_{\cdot j};\varvec{\alpha }_{\cdot \ell }^{(c)})=\prod _{k,h} ({\alpha _{k\ell }^{h}}^{(c)})^{\sum _i s_{ik}^{(c+1)} y_{ij}^{h}}. \end{aligned}$$

M step Updating \(\varvec{\theta }^{(c+1)}\). Denoting \(s_{.k}^{(c+1)}=\sum _i s_{ik}^{(c+1)}\), \(t_{.\ell }^{(c+1)}=\sum _j t_{j\ell }^{(c+1)}\), it leads to

Appendix C: SEM algorithm for categorical data

E and S step

1. 1.

   Computation of \( p(\mathbf{z}|\mathbf{y},\mathbf{w}^{(c)};\varvec{\theta }^{(c)})\), then simulation of \(\mathbf{z}^{(c+1)}\) according to \( p(\mathbf{z}|\mathbf{y},\mathbf{w}^{(c)};\varvec{\theta }^{(c)})\):

   $$\begin{aligned} p(z_i=k|\mathbf{y}_{i\cdot },\mathbf{w}^{(c)};\varvec{\theta }^{(c)})= \frac{\pi ^{(c)}_k \psi _k(\mathbf{y}_{i\cdot };\varvec{\alpha }^{(c)}_{k\cdot })}{\sum _{k'} \pi ^{(c)}_{k'} \psi ^{(c)}_{k'}(\mathbf{y}_{i\cdot };\varvec{\alpha }^{(c)}_{k'\cdot }) }, \end{aligned}$$

   for \(k=1,\ldots , g\) with

   $$\begin{aligned} \psi _k(\mathbf{y}_{i\cdot };\varvec{\alpha }_{k\cdot })=\prod _{\ell ,h} ({\alpha _{k\ell }^{h}}^{(c)})^{\sum _j w_{j\ell }^{(c)} y_{ij}^{h}} \end{aligned}$$
2. 2.

   computation of \(p(\mathbf{w}|\mathbf{y},\mathbf{z}^{(c+1)};\varvec{\theta }^{(c)})\), then simulation of \(\mathbf{w}^{(c+1)}\) according to \(p(\mathbf{w}|\mathbf{y},\mathbf{z}^{(c+1)};\varvec{\theta }^{(c)})\).

M step Denoting \(z_{.k}:=\sum _i z_{ik}\) and \(w_{.\ell }:=\sum _j w_{j\ell }\), it leads to

and

Note that the formulae of VEM and SEM-Gibbs are essentially the same, except that the probabilities \(s_{ik}\) and \(t_{j\ell }\) are replaced with binary indicator values \(z_{ik}\) and \(w_{j\ell }\).

Appendix D: Computing ICL

In this section, the exact ICL expression is derived for categorical data. Using the conditional independence of the \(\mathbf{z}\)s and the \(\mathbf{w}\)s conditionally to \(\varvec{\alpha }, \varvec{\pi }, \varvec{\rho }\), the integrated completed likelihood can be written

Thus

Now, according to the LBM definition,

Since the prior distribution of \(\varvec{\pi }\) is the Dirichlet distribution \(\mathcal{D}(a,\ldots ,a)\)

we have

We recognise a non-normalised Dirichlet distribution \(\mathcal{D}(z_{.1}+a,\ldots ,z_{.g}+a)\) with the normalising factor

The expression of \(p( \mathbf{z})\) directly follows from the Bayes theorem

In the same manner,

We now turn to the computation of \(p(\mathbf{y}|\mathbf{z}, \mathbf{w} )\). Since \(\varvec{\alpha }\) and \((\mathbf{z}, \mathbf{w})\) are independent, we have

and, using the conditional independence of \(y_{ij}\) knowing \(\mathbf{z},\,\mathbf{w}\) and \(\varvec{\alpha }\),

Therefore

because \(p(\alpha _{k\ell })\) is the density of a Dirichlet distribution \(\mathcal{D}(b,\ldots ,b)\). Each term \(k\ell \), is the density of a non-normalised Dirichlet distribution \({\mathcal {D}}(b+N_{k\ell }^1,\ldots ,b+N_{k\ell }^r)\) with the normalising factor

Thus, by the Bayes formula,

And, the ICL criterion, presented in Section 4.1, is straightforwardly derived from equations (3), (8), (9) and (10).

Reprints and permissions