Block clustering of Binary Data with Gaussian Co-variables

TL;DR

This paper introduces a new co-clustering method for binary data that incorporates Gaussian co-variables, improving the grouping of rows and columns in large-scale data analysis.

Contribution

It presents a novel latent block model that integrates co-variables into the co-clustering process, enhancing clustering accuracy.

Findings

Effective on simulated datasets

Successful application to genetic data

Improved clustering performance

Abstract

The simultaneous grouping of rows and columns is an important technique that is increasingly used in large-scale data analysis. In this paper, we present a novel co-clustering method using co-variables in its construction. It is based on a latent block model taking into account the problem of grouping variables and clustering individuals by integrating information given by sets of co-variables. Numerical experiments on simulated data sets and an application on real genetic data highlight the interest of this approach.

Figures

Equations

$$f(\mathbf{x};\mathbf{\theta })=\sum _{\mathbf{u}\in \mathcal{U}}p(\mathbf{u};\mathbf{\theta })f(\mathbf{x}|\mathbf{u};\mathbf{\theta })$$

$$f(\mathbf{x};\mathbf{\theta })=\sum _{(\mathbf{z},\mathbf{w})\in \mathcal{Z}\times \mathcal{W}}p(\mathbf{z};\mathbf{\theta })p(\mathbf{w};\mathbf{\theta })f(\mathbf{x}|\mathbf{z},\mathbf{w};\mathbf{\theta })$$

$$f(\mathbf{x},\mathbf{y};\mathbf{\theta })=\sum _{(\mathbf{z},\mathbf{w})\in \mathcal{Z}\times \mathcal{W}}p(\mathbf{z};\mathbf{\theta })p(\mathbf{w};\mathbf{\theta })f(\mathbf{x}|\mathbf{y},\mathbf{z},\mathbf{w};\mathbf{\theta })f(\mathbf{y}|\mathbf{z};\mathbf{\theta }).$$

$$f(\mathbf{x},\mathbf{y}|\mathbf{z},\mathbf{w};\mathbf{\theta })=\prod _{i,j}f(x_{ij},{\mathbf{y}}_i|{\mathbf{z}}_i,{\mathbf{w}}_j;\mathbf{\theta }).$$

$$f(x_{ij}|{\mathbf{y}}_i,{\mathbf{\beta }}_{z_i{w}_j})=\mathrm{logis}({\beta }_{0,z_i{w}_j}+{\mathbf{\beta }}_{z_i{w}_j}^T{\mathbf{y}}_i{)}^{x_{ij}}{\left(1-\mathrm{logis}({\beta }_{0,z_i{w}_j}+{\mathbf{\beta }}_{z_i{w}_j}^T{\mathbf{y}}_i)\right)}^{1-x_{ij}}$$

$$f(\mathbf{y}|\mathbf{z};\mathbf{\theta })=\prod _i\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i})$$

$$f(\mathbf{x},\mathbf{y}|\mathbf{\theta })=\sum _{(\mathbf{z},\mathbf{w})\in \mathcal{Z}\times \mathcal{W}}\prod _{i,j}{\pi }_{z_i}{\rho }_{w_j}\mathrm{logis}({\mathbf{y}}_i^T{\mathbf{\beta }}_{z_i{w}_j}{)}^{x_{ij}}{\left(1-\mathrm{logis}({\mathbf{y}}_i^T{\mathbf{\beta }}_{z_i{w}_j})\right)}^{1-x_{ij}}\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i}).$$

$$l(\mathbf{\theta })=\log f(\mathbf{x},\mathbf{y};\mathbf{\theta })$$

$$L_C(\mathbf{z},\mathbf{w},\mathbf{\theta })=\sum _k{z}_{.k}\log {\pi }_k+\sum _{\ell }{w}_{.\ell }\log {\rho }_{\ell }+\sum _{i,j,k,\ell }{z}_{ik}{w}_{j\ell }\log f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}_{k\ell }).$$

$$Q(\mathbf{\theta },{\mathbf{\theta }}^{(c)})=\mathbb{E}\left[L_C(\mathbf{z},\mathbf{w},\theta )|\mathbf{x},\mathbf{y},{\mathbf{\theta }}^{(c)}\right]=\sum _{i,k}{t}_{ik}^{(c)}\log {\pi }_k+\sum _{j,\ell }{r}_{j\ell }^{(c)}\log {\rho }_{\ell }+\sum _{i,j,k,\ell }{e}_{ikj\ell }^{(c)}\log f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}_{k\ell })$$

$$t_{ik}^{(c)}=P(z_{ik}=1|\mathbf{x},\mathbf{y},{\mathbf{\theta }}^{(c)}),r_{jl}^{(c)}=P(w_{j\ell }=1|\mathbf{x},\mathbf{y},{\mathbf{\theta }}^{(c)}),e_{ikj\ell }^{(c)}=P(z_{ik}{w}_{j\ell }=1|\mathbf{x},\mathbf{y},{\mathbf{\theta }}^{(c)})$$

$$l(\mathbf{\theta })={\tilde{F}}_C(\mathbf{t},\mathbf{r};\mathbf{\theta })+KL(q(\mathbf{z},\mathbf{w})\parallel p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},\mathbf{\theta }))$$

$${\tilde{F}}_C(\mathbf{t},\mathbf{r};\mathbf{\theta })=\sum _k{t}_{.k}\log {\pi }_k+\sum _{\ell }{r}_{.\ell }\log {\rho }_l+\sum _{i,j,k,\ell }{t}_{ik}{r}_{j\ell }\log f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}_{k\ell })+H(\mathbf{t})+H(\mathbf{r})$$

$$H(\mathbf{t})=\sum _{ik}{t}_{ik}\log t_{ik},H(\mathbf{r})=\sum _{jl}{r}_{jl}\log r_{jl}.$$

$${\mathrm{argmax}}_{\mathbf{t},\mathbf{r},\mathbf{\theta }}{\tilde{F}}_C(\mathbf{t},\mathbf{r},\mathbf{\theta }).$$

$${\pi }_k=\frac{t_{.k}}{n}\text{ and }{\rho }_{\ell }=\frac{r_{.\ell }}{m}.$$

$${\mathbf{\beta }}_{kl}={\mathrm{argmax}}_{\mathbf{\beta }}\sum _{ij}{t}_{ik}{r}_{jl}\log f(x_{ij}|{\mathbf{y}}_i;\mathbf{\beta }),k=1,\dots ,g,l=1,\dots ,d.$$

$${\mathbf{\mu }}_k=\frac{1}{t_{.k}}\sum _i{t}_{ik}{\mathbf{y}}_i\text{mbox}and{\mathbf{\Sigma }}_k=\frac{1}{t_{.k}}\sum _i{t}_{ik}({\mathbf{y}}_i-{\mathbf{\mu }}_k)({\mathbf{y}}_i-{\mathbf{\mu }}_k{)}^T.$$

$$t_{ik}^{(c+1)}=\frac{{\pi }_k^{(c)}\prod _{jl}{\left(f(x_{ij}|{\mathbf{y}}_i;{\mathbf{\beta }}_{kl}^{(c)})\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_k^{(c)},{\mathbf{\Sigma }}_k^{(c)})\right)}^{r_{jl}^{(c)}}}{{\sum }_k{\pi }_k^{(c)}\prod _{jl}{\left(f(x_{ij}|{\mathbf{y}}_i;{\mathbf{\beta }}_{kl}^{(c)})\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_k^{(c)},{\mathbf{\Sigma }}_k^{(c)})\right)}^{r_{jl}^{(c)}}}.$$

$$r_{jl}^{(c+1)}=\frac{{\rho }_l^{(c)}\prod _{ik}f(x_{ij}|{\mathbf{y}}_i;{\mathbf{\beta }}_{kl}^{(c+1/2)}{)}^{t_{ik}^{(c+1)}}}{{\sum }_l{\rho }_l^{(c)}\prod _{ik}f(x_{ij}|{\mathbf{y}}_i;{\mathbf{\beta }}_{kl}^{(c+1/2)}{)}^{t_{ik}^{(c+1)}}}.$$

$$\mathrm{BIC}=-2\underset{\mathbf{\theta }}{\max }l(\mathbf{\theta })+D\log (N)$$

$$\mathrm{BIC}(g,d)=-2\underset{\mathbf{\theta }}{\max }\log f(\mathbf{x},\mathbf{y};\mathbf{\theta })+(g-1)\log n+\lambda \log n+(d-1)\log m+gd(p+1)\log (mn)$$

$$f(\mathbf{y}|\mathbf{x},\mathbf{z},\mathbf{w},\mathbf{\theta })$$

$$\propto \prod _{i=1}^n{\pi }_{z_i}\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i})\prod _{l=1}^d{\rho }_l^{n_l}\frac{e^{n_{il}{\mathbf{y}}_i^T{\mathbf{\beta }}_{z_{i}l}}}{{\left(1+e^{{\mathbf{y}}_i^T{\mathbf{\beta }}_{z_{i}l}}\right)}^{m_l}}$$

$$f(\mathbf{y}|\mathbf{x},\mathbf{z},\mathbf{w},\mathbf{\theta })\propto \prod _{j=1}^m{\rho }_{w_j}\prod _{k=1}^g{\pi }_k^{m_j}\prod _{i:z_i=k}\mathrm{logis}({\mathbf{y}}_i^T{\mathbf{\beta }}_{k{w}_j}{)}^{x_{ij}}{\left(1-\mathrm{logis}({\mathbf{y}}_i^T{\mathbf{\beta }}_{k{w}_j})\right)}^{1-x_{ij}}\varphi ({\mathbf{y}}_i;{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k).$$

$$I(j)$$

$$=\log {\rho }_{w_j}+\sum _{i=1}^n\left(x_{ij}{\mathbf{y}}_i^T{\mathbf{\beta }}_{z_i{w}_j}-\log (1+\exp ({\mathbf{y}}_i^T.{\mathbf{\beta }}_{z_i{w}_j}))\right)$$

$${\hat{w}}_j=\arg \underset{l=1}{\overset{d}{\max }}{r}_{jl}.$$

$$\frac{\partial {\tilde{F}}_C(\mathbf{t},\mathbf{r};\mathbf{\theta })}{\partial t_{ik}}=\log {\pi }_k+\sum _{j,\ell }{r}_{j\ell }\log f_{k\ell }(x_{ij},{\mathbf{y}}_i;\mathbf{\theta })-\log t_{ik}-1.$$

$$t_{ik}^{(c+1)}=\frac{{\pi }_k^{(c)}{\prod }_{j,l}{\left[f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}^{(c)})\right]}^{r_{jl}^{(c)}}}{{\sum }_k{\prod }_{j,l}{\left[f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}^{(c)})\right]}^{r_{jl}^{(c)}}}.$$

$$\log (t_{ik}^{(c+1)})\propto \log ({\pi }_k^{(c)})+\sum _{j,l}{r}_{jl}^{(c)}\log f(x_{ij},{\mathbf{y}}_i;{\mathbf{\theta }}^{(c)}).$$

Full text

1 Introduction

Classification is a method of data analysis that aims to group together a set of observations into homogeneous classes. It plays an increasingly important role in many scientific and technical fields. Its aim is the automatic resolution of problems by decision-making based on the observations and to define the rules for classifying objects depending on qualitative or quantitative variables.

Clustering is the most popular technique for data analysis in many disciplines. In recent years, co-clustering has been increasingly used. Unlike classical clustering, which groups similar objects from a single collection of objects, co-clustering or bi-clustering [1] aims at simultaneously grouping objects from two disjoint sets, thus revealing interactions between elements of two sets.

It is most often used with bipartite spectral graphing partitioning methods in the field of extracting text data [2] by simultaneously grouping documents and content (words) and analyzing huge corpora unlabeled documents [3] to simultaneously understand aggregates of subsets of web users sessions and information from the page views. Co-clustering algorithms have also been developed for computer vision applications. It is used for grouping images simultaneously with their low-level visual characteristics and for content-based search [4].

In this paper we extend co-clustering methods allowing simultaneous detection of associations between variables and individuals by taking into account co-variables. These co-variables can be additional measures of interest. Consideration of a co-variable is expected to provide better separation of groups of variables and especially groups of individuals. Classification quality is determined by general validation measures specific to the co-clustering method. This approach can be useful when co-clustering a set $\mathbf{X}$ of variables and individuals in coherence with an independent $\mathbf{Y}$ variable measured on these same individuals. For example, in the co-clustering of several SNP (Single-Nucleotide Polymorphism) variables on different patients with respect to a measured phenotype (see application in section 3).

The paper is organized as follows. In the first part, we explain the principle of block mixture models through section 2. The latent block model for binary variable takes into account co-variables and the model parameters estimation is proposed in Section 2.2. The parameter estimation method is described in section 2.3. The choice of the optimal number of blocks and the measure of influence of each variable on the co-variable $\mathbf{Y}$ is presented in the second part (section 2.5 and 2.6). The method is illustrated on simulated and real genetic data in the last part (section 3).

2 Block mixture models

2.1 Classical latent block model

Let $\mathbf{x}$ be a data set doubly indexed by a set $I$ with $n$ elements (individuals) and a set $J$ with $m$ elements (variables). We represent a partition of $I$ into $g$ clusters by $\mathbf{z}=(z_{11},\ldots,z_{ng})$ with $z_{ik}=1$ if $i$ belongs to cluster $k$ and $z_{ik}=0$ otherwise, $z_{i}=k$ if $z_{ik}=1$ and we denote by $z_{.k}=\sum_{i}z_{ik}$ the cardinality of row cluster $k$. Similarly, we represent a partition of $J$ into $d$ clusters by $\mathbf{w}=(w_{11},\ldots,w_{md})$ with $w_{j\ell}=1$ if $j$ belongs to cluster $\ell$ and $w_{j\ell}=0$ otherwise, $w_{j}=\ell$ if $w_{j\ell}=1$ and we denote $w_{.\ell}=\sum_{j}w_{j\ell}$ the cardinality of column cluster $\ell$.

The block mixture model formulation is defined in [5] and [6] (among others) by the following probability density function

[TABLE]

where ${\cal U}$ denotes the set of all possible labels of $I\times J$ and $\boldsymbol{\theta}$ contains all the unknown parameters of this model. By restricting this model to a set of labels of $I\times J$ defined by a product of labels of $I$ and $J$, and further assuming that the labels of $I$ and $J$ are independent of each other, one obtain the decomposition

[TABLE]

where ${\cal Z}$ and ${\cal W}$ denote the sets of all possible labellings $\mathbf{z}$ of $I$ and $\mathbf{w}$ of $J$. Equation (1) define a Latent Block Model (LBM).

2.2 LBM for binary variables with co-variables: General formulation

From now, we assume that $\mathbf{x}$ is a binary data set. Let $\mathbf{y}$ represents a data-set (co-variables) of $\mathbb{R}^{p}$ indexed by $I$. In order to take into account this set of co-variables the classical block model formulation is extended to propose a block mixture model defined by the following probability density function

[TABLE]

By extending the latent class principle of local independence to our block model, each data pair $(x_{ij},\mathbf{y}_{i})$ will be independent once $z_{i}$ and $w_{j}$ are fixed. Hence we have

[TABLE]

We choose to model the dependency between $x_{ij}$ and $\mathbf{y}_{i}$ using the canonical link for binary response data

[TABLE]

with $(\beta_{0},\boldsymbol{\beta}_{k,l})\in\mathbb{R}^{p+1}$ and $\operatorname{logis}(x)=e^{x}/(1+e^{x})$. Each data point $\mathbf{y}_{i}$ will be independent once $z_{i}$ are fixed. In the examples presented in section 3, we choose

[TABLE]

with $\phi$ denoting the multivariate Gaussian density in $\mathbb{R}^{p}$.

In order to simplify the notation, we add a constant coordinate $1$ to vectors $\mathbf{y}_{i}$ and write $\boldsymbol{\beta}_{k,l}$ in the latter rather than $(\beta_{0,k,l},\boldsymbol{\beta}_{k,l})$.

The parameters are thus $\boldsymbol{\theta}=(\boldsymbol{\pi},\boldsymbol{\rho},\boldsymbol{\beta},\boldsymbol{\mu},\boldsymbol{\Sigma})$, where $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{g})$, $\boldsymbol{\rho}=(\rho_{1},\ldots,\rho_{d})$ are the vectors of probabilities $\pi_{k}$ and $\rho_{\ell}$ that a row and a column belong to the $k$th row component and to the $\ell$th column component respectively, $\boldsymbol{\beta}=(\boldsymbol{\beta}_{kl})$ are the coefficients of the logistic function, $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are the means and variances of the Gaussian density. In summary, we obtain the latent block mixture model with pdf

[TABLE]

Using the above expression, the randomized data generation process can be described by the four steps row labellings (R), column labellings (C), co-variable data generation (Y) and data generation (X) as follows:

1. (R)

Generate the labellings $\mathbf{z}=(z_{1},\ldots,z_{n})$ according to the distribution $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{g})$. 2. (C)

Generate the labellings $\mathbf{w}=(w_{1},\ldots,w_{m})$ according to the distribution $\boldsymbol{\rho}=(\rho_{1},\ldots,\rho_{d})$. 3. (Y)

Generate for $i=1,...,n$ vector $\mathbf{y}_{i}$ according to the Gaussian distribution $\mathcal{N}_{p}(\boldsymbol{\mu}_{z_{i}},\boldsymbol{\Sigma}_{z_{i}})$. 4. (X)

Generate for $i=1,...,n$ and $j=1,...,m$ a value $x_{ij}$ according to the Bernoulli distribution $f(x_{ij}|\mathbf{y}_{i};\boldsymbol{\beta}_{z_{i}w_{j}})$ given in (3).

2.3 Model Parameters Estimation

The complete data is represented as a vector $(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{w})$ where unobservable vectors $\mathbf{z}$ and $\mathbf{w}$ are the labels. The log-likelihood to maximize is

[TABLE]

and the double missing data structure, namely $\mathbf{z}$ and $\mathbf{w}$, makes statistical inference more difficult than usual. More precisely, if we try to use an EM algorithm as in standard mixture model [7] the complete data log-likelihood is found to be

[TABLE]

The EM algorithm maximizes the log-likelihood $l(\boldsymbol{\theta})$ iteratively by maximizing the conditional expectation $Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(c)})$ of the complete data log-likelihood given a previous current estimate $\boldsymbol{\theta}^{(c)}$ and $(\mathbf{x},\mathbf{y})$:

[TABLE]

where

[TABLE]

Unfortunately, difficulties arise due to the dependence structure in the model, in particular to determine $e_{ikj\ell}^{(c)}$. The assumed independence of $\mathbf{z}$ and $\mathbf{w}$ in (1) is not conserved by the posterior probability.

To solve this problem an approximate solution is proposed in [5] using the [8] and [9] interpretation of the VEM algorithm. Consider a family of probability distribution $q(z_{ik},w_{j\ell})$ verifying $q(z_{ik},w_{j\ell})>0$ and the relation $q(z_{ik},w_{j\ell})=q(z_{ik})q(w_{j\ell})$, for all $i,j,k,l$. Set $t_{ik}=q(z_{ik})$ and $r_{jl}=q(w_{j\ell})$, $\mathbf{t}=(t_{ik})_{ik}$ for $i=1,\ldots,n$, $k=1,\ldots,g$ and $\mathbf{r}=(r_{jl})_{jl}$ for $j=1,\ldots,m$ and $l=1,\ldots,d$. One shows easily that

[TABLE]

with $KL(q\parallel p)$ denoting the Kullback-Liebler divergence of distribution $p$ and $q$,

[TABLE]

and $H(\mathbf{t})$, $H(\mathbf{r})$ denoting the entropy of $\mathbf{t}$ and $\mathbf{r}$, i.e.

[TABLE]

$\tilde{F}_{C}$ is called the free energy or the fuzzy criterion. As the Kullback-Liebler divergence is always positive, the fuzzy criterion is a lower bound of the log-likelihood and is used as a replacement for it. Doing that, the maximization of the likelihood $l(\boldsymbol{\theta})$ is replaced by the following problem

[TABLE]

This maximization can be achieved using the BEM algorithm detailed hereafter.

2.4 Block expectation maximization (BEM) Algorithm

The fuzzy clustering criterion given in (8) can be maximized using a variational EM algorithm (VEM). We here outline the various expressions evaluated during E and M steps.

E-Step:

we compute either the values of $\mathbf{t}$ (respectively $\mathbf{r}$) with $\mathbf{r}$ (respectively $\mathbf{t}$) and $\boldsymbol{\theta}$ fixed (formulas (12), (13) hereafter). Details are given in appendix A.

M-Step:

we calculate row proportions $\boldsymbol{\pi}$ and column proportions $\boldsymbol{\rho}$. The maximization of $\tilde{F}_{C}$ w.r.t. $\boldsymbol{\pi}$, and w.r.t $\boldsymbol{\rho}$, is obtained by maximizing $\sum_{k}t_{.k}\log\pi_{k}$, and $\sum_{\ell}r_{.\ell}\log\rho_{\ell}$ respectively, which leads to

[TABLE]

Also, for $\mathbf{t}$, $\mathbf{r}$ fixed, the estimate of model parameters $\boldsymbol{\beta}$ will be obtained by maximizing

[TABLE]

Detail are given in appendix B. Finally parameters of the Gaussian density are given by the usual formulas

[TABLE]

Putting everything together, we obtain the BEM algorithm.

BEM algorithm:

Using the E and M steps defined above, BEM algorithm can be enumerated as follows:

Initialization

Set $\mathbf{t}^{(0)},\mathbf{r}^{(0)}$ and $\boldsymbol{\theta}^{(0)}=(\boldsymbol{\pi}^{(0)},\boldsymbol{\rho}^{(0)},\boldsymbol{\beta}^{(0)},\boldsymbol{\mu}^{(0)}$, $\boldsymbol{\Sigma}^{(0)}).$

(a) Row-EStep

Compute $\mathbf{t}^{(c+1)}$ using formula

[TABLE]

(b) Row-MStep

Compute $\boldsymbol{\pi}^{(c+1)}$, $\boldsymbol{\mu}^{(c+1)}$, $\boldsymbol{\Sigma}^{(c+1)}$ using equations (9) and (11) and estimate $\boldsymbol{\beta}^{(c+1/2)}$ by solving maximization problem (10).

(c) Col-EStep

Compute $\mathbf{r}^{(c+1)}$ using formula

[TABLE]

Observe that $r_{jl}$ does not depend of the density of $\mathbf{y}$.

(d) Col-MStep

Compute $\boldsymbol{\rho}^{(c+1)}$ using equations (9) and estimate $\boldsymbol{\beta}^{(c+1)}$ by solving maximization problem (10).

Iterate

Iterate (a)-(b)-(c)-(d) until convergence.

2.5 Selecting the number of blocks

BIC is an information criterion defined as an asymptotic approximation of the logarithm of the integrated likelihood ([10]). The standard case leads to write BIC as a penalised maximum likelihood:

[TABLE]

where $N$ is the number of statistical units and $D$ the number of free parameters and $l(\boldsymbol{\theta})$ defined in (5). Unfortunately, this approximation cannot be used for LBM, due to the dependency structure of the observations $(\mathbf{x},\mathbf{y})$. However, a heuristic have been stated to define BIC in [11] and [12]. BIC-like approximations ICL lead to the following approximation as $n$ and $m$ tend to infinity

[TABLE]

with $\lambda$ the number of parameters of the $\mathbf{y}$ distribution. For LBM, the intractable likelihood $f(\mathbf{x},\mathbf{y};\boldsymbol{\theta})$ is replaced by the maximized free energy $\tilde{F}_{C}$ in (8) obtained by the BEM algorithm.

2.6 Measuring Influence of a Variable

Let $j$ be fixed (a column of the matrix $\mathbf{x}$). We would like to measure the effect of the variable $\mathbf{x}^{j}=(x_{ij})_{i=1}^{n}$ on $\mathbf{y}$. It is possible to obtain a measure of this effect by looking to the posterior probability of $\mathbf{y}$.

Lemma 1

Let $(\mathbf{x},\mathbf{z},\mathbf{w})$ fixed. For $l=1,\ldots,d$ let $m_{l}$ denotes the number of columns with label $l$, i.e $m_{l}=\#\{w_{jl}=1,\;j=1,\ldots m\}$ and for a row $i$ fixed let $m_{il}$ denotes the number of elements such that ${w_{jl}=1}$ and $x_{ij}=1$, i.e. $m_{il}=\#\{w_{jl}x_{ij}=1,\;j=1,\ldots n\}$. The posterior probability of the co-variable $\mathbf{y}$ is

[TABLE]

Alternatively, for $k=1,\ldots,g$, let $n_{k}$ denotes the number of rows with label $k$, i.e. $n_{k}=\#\{z_{ik}=1,i=1,\ldots,m\}$. The posterior probability of the co-variable $\mathbf{y}$ is

[TABLE]

The proof of this lemma is straightforward and therefore omitted.

Assuming $\mathbf{z}$ and $\mathbf{w}$ known, we measure the influence of a variable using its contribution to the posterior probability. Fixing $j$, taking the logarithm and eliminating terms independent of $\mathbf{x}^{j}$, we obtain the influence measure criteria

[TABLE]

which is interpreted as the $\log$-contribution to the posterior distribution (16) of the variable $\mathbf{x}^{j}$. Replacing the unknown labels $w_{j}$ and $z_{i}$ by their MAP estimators $\hat{w}_{j}$ and $\hat{z}_{i}$, we are able to sort the variables from the most to the less influential.

3 Examples

3.1 Simulated data

3.1.1 Computational time

We compute 80 times the elapsed time of the model for various configurations of the parameter on a HP Zbook G3. The (averaged) computing time as a function of $m$ when $g=2$ for different values of $m$ (the number of columns) and when $d$ (the number of cluster in columns) take values 2 and 6 is plotted in figure 1 below

We can observe that as $n$ grows the elapsed time grows linearly, but that the slope increases as $d$ (the number of class in columns) is increased.

3.1.2 Error rate

Next we simulate 80 times the number of columns well classified when $g=2$ and for various configurations of $m$ and $d$. The cluster of a column is estimated using the maximum a posterior (MAP) estimator

[TABLE]

From these partial results, we see that the number of bad classified columns labels increases as $d$ increases while it remains relatively constant with $m$. An other salient feature is that when the number of individuals ($n$) is greater, this error rate is lower. The number of well classified rows is stable near 0.9 for all tested configurations of the parameters and is not displayed.

3.2 Real Data Analysis

Here, we study data from an epidemiological and genetic survey of malaria disease in Senegal. Data were collected between 1990 to 2008. We worked on a dataset including $n=885$ individuals with measured malaria risk score (phenotype) and genotype available on several candidate genes for susceptibility/resistance to the disease. A total of $m=45$ Single Nucleotide Polymorphisms (SNPs) was considered across these genes and was used as genetic variables. The malaria risk score was a quantitative measure normally distributed and was considered as a co-variable for this co-clustering method. The SNPs are coded in dominant effect on the disease risk. Using the BIC criteria (see graph 3), we choose to focus on the model $d=2$ groups of individuals and $g=11$ groups of SNPs.

3.2.1 Analysis for phenotype data

The choice of a mixture model or not depends on the application context. In the case of genetic data, we are often interested in the comparison of the susceptible and the resistant to a given phenotype. In this application, we look for genes to explain the difference between susceptible and resistant which justifies the use of a mixture model on the target variable. After block-clustering, we find that the individuals are divided in two groups: the susceptibility category composed of a group of individuals with a value of phenotype essentially greater than zero and the resistant category composed of a group of individuals with a value of phenotype essentially less than zero (see figure 4).

Observe how the marginal distribution of the phenotype, which is uni-modal, becomes multi-modal when conditioned by ($\mathbf{x},\mathbf{z}$).

3.2.2 Analysis for genotypes data

We looked at the SNPs to determine which ones would potentially be involved in malaria susceptibility / resistance.

The proposed methodology allowed the selection of the most significant SNPs according to the influence measure proposed in section 2.6. The most frequent SNPs are grouped into the following classes: class 1 and 9. It is noted that the SNPs of these classes have been shown in the literature to have a high significance effect on malaria. Most G6P and hemoglobin SNPs are grouped into these 2 classes. Reviews from exiting literature gives us: Glucose-6-phosphate dehydrogenase (G6PD) deficiency is prevalent in sub-Saharan African populations and has been associated with protection against severe malaria [13, 14, 15, 16]. Studies above haplotype analysis reveal that the G6PD locus is an under-balanced selection, suggesting a malaria protection mechanism based on modest frequency alleles and avoiding parasite attachment [14]. Hemoglobins S and C (HbS and HbC respectively) are known to be two structurally variant forms of normal adult hemoglobin (HbA) resulting from distinct mutations in the $\beta$-globin gene. The protective effect of HbS against Plasmodium falciparum malaria has been shown by several authors [17, 18, 19]. In the case of HbC, the protection is highest in homozygous individuals with HbCC. The proposed model confirmed the strong link between sickle cell polymorphism (HBS), blood group ABO (HBC) and falciparum malaria in the West African population.

3.2.3 Association between phenotype and genotypes

The most common approach used in genetic data is the GWAS method (Genome Wide Association Studies). This method makes a linear regression of the quantitative phenotype on each genotype variable. By applying co-clustering with the phenotype as co-variable, we could obtain a dichotomy of the phenotype. This dichotomy allows us to divide individuals into two categories: susceptible and resistant. In this part, we compare the results of GWAS studies between the quantitative phenotype, the binary phenotype ($\mathbf{1}_{y_{i}\leq 0}$) and the (co-)clustered phenotype. Figure 6 shows that there are more signals at the 5% threshold for the clustered phenotype compared to the two other phenotypes. In summary the proposed methodology allows to detect more significant SNPs compared to the quantitative and binary phenotype.

4 Conclusion

In this article, our main contribution has been to develop a co-clustering model taking into account a (mixture of) Gaussian co-variable. Applications have been made on simulated and real data sets. Our preliminary results are confirmed in previous studies in Africa. The method offers good classification performance on complex data sets (large number of variables and classes). This method can be useful in a wide variety of classification problems with Gaussian predictors and will allow us to discover new patterns of genes allowing to understand and evaluate the mechanism existing between genetics and malaria in an African population particularly in a Senegalese rural area. Further analysis could be done with more SNPs in another paper in preparation. Estimation is performed using a R package (with computational part in C++) that will be soon be available on the CRAN website https://cran.r-project.org/. Meanwhile the package is available on demand to the authors.

Appendix A Computing the (rows and columns) E-Step

For the E-Step $t_{ik}$ value maximize the fuzzy criterion given in equation (8). Derivative with respect to $t_{ik}$ gives

[TABLE]

Equating this equation to zero, taking exponential and recalling that $\sum_{k}t_{ik}=1$, we obtain that $t_{ik}$ is updated as

[TABLE]

For numerical reason, we prefer to compute the logarithm of this expression which is

[TABLE]

Recall that (see equation 3)

[TABLE]

giving

[TABLE]

Similar computation gives for $r_{jl}$

[TABLE]

Observe that the Gaussian distribution does not depend of $j$ nor $l$. This term become constant when summing over $i$ and $k$ and disappears when $r_{jl}$ values are normalized.

Appendix B Computing the M-Step

For the M-Step, we use a Newton-Raphson algorithm in order to solve the equation (10). For each pair $(k,l)$ the function to maximize can be written

[TABLE]

The first derivative with respect to the d-th coordinate $\beta_{d}$ is

[TABLE]

giving the following expression for the gradient

[TABLE]

with $Y=\left[\mathbf{y}_{i}\right]_{i=1}^{N}$, $X=\left[\sum_{j}r_{jl}x_{ij}\right]_{i=1}^{N}$, $\boldsymbol{\mu}=\left[r_{.l}\frac{\exp(\mathbf{y}_{i}^{T}.\boldsymbol{\beta})}{1+\exp(\mathbf{y}_{i}^{T}\boldsymbol{\beta})}\right]_{i=1}^{N}$, $D=\mathop{\mathrm{diag}}(t_{ik})_{i=1}^{N}$ The second derivative with respect to $\beta_{d}$ and $\beta_{d^{{}^{\prime}}}$ is

[TABLE]

giving the following expression for the hessian

[TABLE]