Gabriel Frisch

Jean-Benoist Leger

Yves Grandvalet

Learning from missing data

HAL is a

Introduction

Biclustering or co-clustering simultaneously groups the rows and the columns of a data matrix. Co-clustering has found applications in many areas such as genomic analysis [START_REF] Pontes | Biclustering on expression data: A review[END_REF][START_REF] Kluger | Spectral biclustering of microarray data: Coclustering genes and conditions[END_REF]], text analysis [START_REF] Inderjit | Information-theoretic co-clustering[END_REF], Selosse et al., 2020b], collaborative filtering [START_REF] George | A scalable collaborative filtering framework based on co-clustering[END_REF]Merugu, 2005, Shan and[START_REF] Shan | Bayesian co-clustering[END_REF], or political analysis [Latouche et al., 2011, Wyse and[START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Co-clustering methods can be divided into categories such as, but not limited to, spectral methods [START_REF] Inderjit | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF][START_REF] Kluger | Spectral biclustering of microarray data: Coclustering genes and conditions[END_REF], mutual information methods [START_REF] Inderjit | Information-theoretic co-clustering[END_REF], modularity based methods [START_REF] Labiod | Co-clustering for binary and categorical data with maximum modularity[END_REF], non negative matrix tri-factorization [START_REF] Ding | Orthogonal nonnegative matrix t-factorizations for clustering[END_REF] or model-based methods. Among the model-based methods, the Latent Block Model [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF][START_REF] Nadif | Latent block model for contingency table[END_REF][START_REF] Lomet | Sélection de modèle pour la classification croisée de données continues[END_REF][START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] relies on mixtures, assuming that the observations are generated from finite mixture components in rows and columns.

Most standard methods of clustering or co-clustering presuppose complete information and cannot be applied with missing data, or may provide misleading conclusions when missingness is informative. A careful examination of the data generating process is necessary for the processing of missing values, which requires identifying the type of missingness [START_REF] Rubin | Inference and missing data[END_REF]: Missing Completely At Random (MCAR) refers to the mechanism in which the probability of being missing does not depend on the variable of interest or any other observed variable; whereas in Missing At Random (MAR) the probability of being missing depends on some observed data but is still independent from the non-observed data; and finally Missing Not At Random (MNAR) refers to the mechanism in which the probability of being missing depends on the actual value of the missing data. Under the MAR hypothesis, no information on the generation of data can be extracted from its absence, but under a MNAR assumption, this absence is informative, and ignoring this information in likelihood-based imputation methods may lead to strong biases in estimation [Little and Rubin, 1986]. Missing Not At Random is also known as non-ignorable missingness, in opposition to the ignorable missingness of MCAR and MAR settings, as the absence of data is assumed to convey some information.

In this paper, we aim at clustering the rows and columns of a data matrix whose entries are missing not at random. Equivalently, we consider the clustering of the vertices of a bipartite graph whose edges are missing not at random. For this purpose, we introduce a co-clustering model that combines a MNAR missingness model with the Latent Block Model (LBM).

In Section 2 we review the Latent Block Model introduced by Govaert and Nadif [2008]. In Section 3, we introduce our model, a LBM extended to a MNAR missingness process, and propose, in Section 4, a variational EM algorithm to infer its parameters. We also introduce, in Section 5, an Integrated Completed Likelihood (ICL) criterion to tackle model selection. We then conduct experiments on synthetic datasets in Section 6 to show that the overall approach is relevant to co-cluster MNAR data. Finally, an analysis of the voting records of the lower house of the French Parliament is presented in Section 7.

Related Works

Up to our knowledge, all existing co-clustering methods consider that missing data is either MCAR or MAR [Selosse et al., 2020a[START_REF] Jacques | Model-Based Co-clustering for Ordinal Data[END_REF][START_REF] Papalexakis | From kmeans to higher-way co-clustering: Multilinear decomposition with sparse latent factors[END_REF], except one proposed by [START_REF] Corneli | Co-clustering of ordinal data via latent continuous random variables and not missing at random entries[END_REF] used to co-cluster ordinal data. Their model is very parsimonious as it assumes that both data and missingness are only dependent on the row and column clusters. In this setting, they are able to consider MNAR data even if they suppose that missingness depends indirectly from the value of the data. The model we propose is less parsimonious, thus more flexible, as it supposes that missingness depends both on the value of the data and on the row and column indexes (not only on their respective cluster indexes). In addition to that, our missing data model can be easily re-used for any other statistical co-clustering model as it is weakly-coupled to the generative model of the full data matrix.

In the simple clustering framework, few mixture models handling MNAR data have been proposed. [START_REF] Marlin | Recommender systems, missing data and statistical model estimation[END_REF] combine a multinomial mixture clustering model, used as a complete data model, with a missingness model of type MNAR. They propose two versions of their missingness model. The first one, called CPT-v, models the data observation probability depending only on the underlying value of the data. The second one, called Logit-vd, allows the probability of a data entry to be missing to depend both on the value of the underlying data and the characteristics of the column, giving more flexibility to the model. Our missingness model respects the symmetry of the co-clustering problem by depending identically on the characteristics of the row and column. [START_REF] Kim | Bayesian binomial mixture model for collaborative prediction with non-random missing data[END_REF] propose Bayesian-BM/OR, a simple mixture model of binomials in a Bayesian formalism. The MNAR missingness is modeled by three factors, related to the row, the column and the data value, all three being modeled by Bernoulli variables combined together by a "or" logical operator. The choice of this missingness model is motivated by algorithmic considerations that are not relevant for co-clustering models. [START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF], in a graph perspective, deal with nonobserved dyads during the sampling of a network and consecutive issues in the inference of the stochastic block model. They propose three different MNAR sampling designs in which observing dyads depends either on their underlying value, or on the class or on the degree of the nodes. The Stochastic Block Model, though similar from the Latent Block Model we use, is not usable for co-clustering purposes.

Also related to missing data but not to clustering, MNAR is also considered in the Matrix Factorization framework. [START_REF] Steck | Training and testing of recommender systems on data missing not at random[END_REF] derives a weighted MF model and optimizes the parameters based on a metric that is robust to MNAR data. [START_REF] Miguel Hernández-Lobato | Probabilistic matrix factorization with non-random missing data[END_REF] use a double probabilistic MF model; one is for the complete data and one for the missing data, where users and items propensities are both modeled with low rank matrices. [START_REF] Schnabel | Recommendations as treatments: Debiasing learning and evaluation[END_REF] propose an empirical risk minimization framework to derive a propensity scored matrix factorization method that can account for selection bias.

The Latent Block Model

The Latent Block Model (LBM) is a co-clustering model that classifies jointly the rows and the columns of a data matrix [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF]. This probabilistic generative model assumes a double partition on the rows and the columns of a (n 1 × n 2 ) data matrix X that corresponds to a strong structure of the matrix in homogeneous blocks. This structure is unveiled by reordering the rows and columns according to their respective cluster index; for k 1 row clusters and k 2 column clusters, the reordering reveals k 1 ×k 2 homogeneous blocks in the data matrix. Note that we adopt here the original view where the data matrix is interpreted as a data table. The binary matrix X can also be interpreted as the biadjacency matrix of a bipartite graph, whose two sets of vertices corresponds to the rows and columns of the data matrix. In this interpretation, X ij = 1 if an edge is present between "row node" i and "column node" j, and X ij = 0 otherwise.

For the (n 1 × n 2 ) data matrix X, two partitions are defined by the latent variables Y and Z, with Y being the n 1 × k 1 indicator matrix of the latent row clusters (Y iq = 1 if row i belongs to group q and Y iq = 0 otherwise), and Z being the n 2 × k 2 indicator matrix of the latent column cluster. The group indicator of row i will be denoted Y i , and similarly, the group indicator of column j will be denoted Z j . The LBM makes several assumptions on the dependencies:

Independent rows and column clusters The latent variables Y and Z are a priori independent.

p(Y , Z) = p(Y )p(Z) .
Note that a priori independence does not imply a posteriori independence: given the data matrix X, the two partitions are (hopefully) not independent.

Independent and identically distributed row clusters The latent variables Y are independent and follow a multinomial distribution M(1; α), where α = (α 1 , ..., α k1 ) is the mixing proportions of rows:

p(Y ; α) = i p(Y i ; α) p(Y iq = 1; α) = α q , with α ∈ S (k1-1) = {α ∈ R k1 + ∑ q α q = 1}.
Independent and identically distributed column clusters Likewise, the latent variables Z are independent and follow a multinomial distribution M(1; β), where β = (β 1 , ..., β k2 ) is the mixing proportions of columns:

p(Z; β) = j p(Z j ; β) p(Z jl = 1; β) = β l ,
with β ∈ S (k2-1) .

Given row and column clusters, independent and identically distributed block entries Given the row and colum clusters (Y , Z), the entries X ij are independent and follow a Bernoulli distribution of parameter π = (π ql ; q = 1, ..., k 1 ; l = 1, ..., k 2 ): all elements of a block follow the same probability distribution.

p( X Y , Z; π ) = ij p( X ij Y i , Z j ; π ) p( X ij = 1 Y iq Z jl = 1; π ) = π ql .
To summarize, the parameters of the LBM are θ = (α, β, π) and the probability mass function of X can be written as: where φ(X ij ; π ql ) = π Xij ql (1-π ql ) 1-Xij is the mass function of a Bernoulli variable and where I (resp. J) denotes the set of all possible partitions of rows (resp. columns) into k 1 (resp. k 2 ) groups.

p(X; θ) = (Y ,Z)∈I×J ⎛ ⎝ i,q α q Yiq ⎞ ⎠ ⎛ ⎝ j,l β l Z jl ⎞ ⎠ ⎛ ⎝ i,j,q,l φ(X ij ; π ql ) YiqZ jl ⎞ ⎠ , Y i Z j X ij ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∀i, Y i iid ∼ M(1; α) ∀j, Z j iid ∼ M(1; β) ∀i, j, X ij Y iq = 1, Z jl = 1 ind ∼ B(π ql ) with α ∈ S k1-1 , β ∈ S k2-1 and π ql ∈ [0, 1]

Extension to Informative Missing Data

The standard Latent Block Model does not accommodate missing observations, that is, the data matrix X is fully observed. This section introduces our missingness model, which will be coupled to the LBM, thereby enabling to process missing data.

We start by introducing some notation: from now on, X (o) will denote the "partially observed" data matrix, with missing entries, whereas X denotes the "full" (unobserved) data matrix, without missing entries. The partially observed matrix X (o) is identical to the full matrix X except for the missing entries; X (o) takes its values in {0, 1, NA}, where NA denotes a missing value. It will be convenient to introduce a binary mask matrix M that indicates the non-missing entries of X (o)

: if M ij = 0, then X (o) ij = NA.

Models of Missingness

The three main types of missingness are Missing Completely At Random (MCAR), Missing At Random (MAR), and Missing Not At Random (MNAR).

We propose here a model for each missingness type. Instead of directly modeling the probability of being missing, we will model a real variable that defines the log-odds of this probability. This log-odds will be called here the "propensity" to be missing.

Missing Completely At Random (MCAR) Missingness does not depend on data, whether observed or not. A simple model of missingness is obtained by assuming that every entry of X (o) has the same propensity of being missing. This is modeled by a single propensity parameter µ. The graphical representation of this model is shown in Figure 2. ij is generated by the corresponding entries of (i) the full matrix X ij and (ii) the binary mask M ij . The binary mask M does not depend on X and is defined here from a single global effect parameter µ. ij is generated by the corresponding entries of (i) the full matrix X ij and (ii) the binary mask M ij . The binary mask M does not depend on X and is and is defined by a global effect parameter µ and two latent variables A and C that enable deviations from µ.

µ M ij X ij X (o) ij
µ M ij X ij X (o) ij A i C j
Missing At Random (MAR) Missingness depends on the observed data, but not on the unobserved data. The previous missingness model can be enlarged by allowing the propensity of missingness to depend on the row and column indexes. To do so, we can introduce a latent variable for every row, denoted A, and another one for every column, denoted C. For the sake of simplicity, all latent variables A i and C j are assumed independent. They allow deviations from the global propensity µ. The graphical representation of this model is shown in Figure 3.

Missing Not At Random (MNAR) Missingness here depends on unobserved data: the probability of observing the entries of the matrix depends on their values, whether observed or not. We equip the previous model with two additional latent variables to adapt the propensity of each entry of the data matrix to the unobserved data, that is, to X ij . These new row and column latent variables, B and D, adjust the propensity of missingness according to the actual value of X ij . The graphical representation of this model is shown in Figure 4. ij is generated by the corresponding entries of (i) the full matrix X ij and (ii) the binary mask M ij . The binary mask M depends on X and is defined by a global effect parameter µ, two latent variables A and C that enable deviations from µ, and two latent variables B and D, which drive the deviations from the MAR model.

A i B i C j D j µ X ij X (o) ij M ij
We model the latent variables A, B, C, and D with Gaussian distributions centered at zero with free variances σ 2 A , σ 2 B , σ 2 C , and σ 2 D , respectively:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∀i, A i iid ∼ N 0, σ 2 A , B i iid ∼ N 0, σ 2 B ∀j, C j iid ∼ N 0, σ 2 C , D j iid ∼ N 0, σ 2 D .
The global parameter µ and the latent variables define the propensity of missingness, that is, the log-odds of being missing as follows:

∀i, j P ij = µ + A i + B i + C j + D j if X ij = 1 µ + A i -B i + C j -D j if X ij = 0 . (1)
Then, given this propensity, every element M ij of the mask matrix is independent and follows a Bernoulli distribution:

∀i, j M ij A i , B i , C j , D j , X ij ind ∼ B(expit (P ij )) , (2)
with expit(x) = 1 (1 + exp(-x)).

Note that, if we omit the latent variables B i and D j , the missingness model follows the MAR assumption since P ij , and thus M ij , is then independent of X ij . If we also omit the latent variables A i and C j , the missingness model follows the MCAR assumption.

This model of missingness can be used for several applications. One of these, collaborative filtering, uses the history of user ratings to build a recommendation system. For this application, an MCAR modeling means that the probability of observing a rating for a particular item does not depend on the user nor the item; an MAR modeling means that missingness can depend on the user or the item; for example, some people give their opinion more often than others. The MAR simplifying assumption is often used in collaborative filtering. However, [START_REF] Marlin | Collaborative filtering and the missing at random assumption[END_REF] show that there is often a dependency between the rating frequency and the underlying preference level, lending support to the hypothesis that ratings are generated by a MNAR process, where missingness depends on the actual rating that would be given. Some people give their opinion more often when they are satisfied and other ones when they are dissatisfied. Most collaborative filtering methods do not have a principled method for extracting information from missing data, which can lead to strong biases in estimations that may in turn drastically affect predictions [START_REF] Miguel Hernández-Lobato | Probabilistic matrix factorization with non-random missing data[END_REF]. Our missingness model allows to account for the users' propensity to give their opinion, and for the items' propensity to be rated, that is, their notoriety. These propensities could also reflect exogenous factors such as price; for example, more expensive items could be evaluated more often.

LBM with MNAR data

We extend the standard LBM using the previous modeling to MNAR data. Given the full matrix X and the mask matrix M , all the elements of the observed matrix X (o) are independent and identically distributed:

X (o) ij X ij , M ij = X ij if M ij = 1 NA if M ij = 0 . (3)
Figure 5 summarizes the LBM extented to MNAR data. X o) taking its values in (0, 1, NA), the same model can be rewritten with a Categorial distribution using directly the latent variables of both models:

∀i, Y i iid ∼ M(1; α) ∀j, Z j iid ∼ M(1; β) ∀i, A i iid ∼ N 0, σ 2 A ∀i, B i iid ∼ N 0, σ 2 B (4) ∀j, C j iid ∼ N 0, σ 2 C ∀j, D j iid ∼ N 0, σ 2 D ∀i, j, X (o) ij Y iq = 1, Z jl = 1, A i , B i , C j , D j ind ∼ Cat ⎛ ⎜ ⎝ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 1 NA ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ p 0 p 1 1 -p 0 -p 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎞ ⎟ ⎠ with p 0 = (1 -π ql ) expit (µ + A i -B i + C j -D j )
(5)

p 1 = π ql expit (µ + A i + B i + C j + D j ) . ( 6
) Y i Z j A i B i C j D j LBM MNAR missingness model X ij M ij X (o) ij Latent Block Model Y i iid ∼ M(1; α), α ∈ S k1-1 Z j iid ∼ M(1; β), β ∈ S k2-1 ( X ij Y iq = 1, Z jl = 1 ) ind ∼ B(π ql ), π ql ∈ (0, 1) MNAR model A i iid ∼ N 0, σ 2 A , σ 2 A ∈ R * + B i iid ∼ N 0, σ 2 B , σ 2 B ∈ R * + C j iid ∼ N 0, σ 2 C , σ 2 C ∈ R * + D j iid ∼ N 0, σ 2 D , σ 2 D ∈ R * + ( M ij A i , B i , C j , D j , X ij = 1 ) ind ∼ B(expit (µ + A i + B i + C j + D j )) ( M ij A i , B i , C j , D j , X ij = 0 ) ind ∼ B(expit (µ + A i -B i + C j -D j )) Observations are generated according to: X (o) ij X ij , M ij = X ij if M ij = 1 NA if M ij = 0

Inference in the extented LBM

The dependency between the full data matrix X and the mask matrix M requires a joint inference of the LBM with the MNAR model. As the standard maximum likelihood approach cannot be applied directly, we adopt a strategy based on a variational EM.

During inference, we use the reformulation (Equation 4). We can split our random variables into two sets: the set of unobserved latent variables and the set of observed variables consisting of X (o) only. An observation of X (o) only is called the incomplete data, and an observation of X (o) together with the latent variables A, B, C, D, Y and Z is called the complete data. Given the incomplete data, our objective is to infer the model parameters θ via maximum likelihood θ = arg max θ p(X (o) ; θ). We resort to the Expectation Maximization (EM) algorithm to maximize p(X (o) ; θ) without explicitely calculating it. The EM algorithm iteratively applies the two following steps:

E-step Expectation step: from the current estimate θ (t) of θ, compute the cri- terion Q(θ θ (t)
) defined as the expectation of the complete log-likelihood, conditionally on the observations X (o) :

Q(θ θ (t) ) = E Y ,Z,A,B,C,D X (o) ,θ (t) log p X (o) , Y , Z, A, B, C, D; θ M-step Maximization step: find the parameters that maximize Q(θ θ (t) ). θ (t+1) = arg max θ Q(θ θ (t) )
The computation of the complete log-likelihood at the E-step requires the posterior distribution of the latent variables p(Y , Z, A, B, C, D X (o) ) which is intractable, because the search space of the latent variables is combinatorially too large. This problem is well known in the context of co-clustering; for the Latent Block Model, [START_REF] Celeux | The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem[END_REF], [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] propose a stochastic E-step with Monte Carlo sampling, but this strategy is not suited to large-scale problems. We follow the original strategy proposed by [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF], which relies on a variational formulation of the problem, since it is more efficient in high dimension.

Variational EM

The variational EM (VEM) [START_REF] Michael | An introduction to variational methods for graphical models[END_REF][START_REF] Jaakkola | Tutorial on variational approximation methods[END_REF] introduces q(⋅), a parametric inference distribution defined over the latent variables Y , Z, A, B, C, D and optimize the following lower bound on the log-likelihood of the incomplete data:

J (q, θ) = log p X (o) ; θ -KL q(⋅) ∥ p(⋅ X (o) ; θ) , (7)
where KL stands for the Kullback-Leibler divergence and q(⋅) denotes the variational distribution over the latent variables Y , Z, A, B, C and D. It can be shown that J (q, θ) is a concave function of the variational distribution q and that its maximum is reached for q(⋅) = p(⋅ X (o) ; θ). Thus, maximizing the criterion J is equivalent to minimizing the discrepancy between q(⋅) and p(⋅ X (o) ; θ), as measured by the Kullback divergence, and is also equivalent to maximizing the likelihood. The minimization of this Kullback divergence requires to explore the whole space of latent distributions; the difficulty of the problem is equivalent, in terms of complexity, to the initial problem. The criterion J (q, θ) can also be expressed as the sum of a negative "energy" and the entropy of q hence its name "negative variational free energy" in analogy with the thermodynamic free energy:

J (q, θ) = H(q) + E q log p X (o) , Y , Z, A, B, C, D; θ , (8)
where H(q) is the entropy of the variational distribution and E q is the expectation with respect to the variational distribution. The criteria J can become tractable if an exploration of a subspace, noted q γ , of the latent distributions is made. However, this solution comes with the cost that the maximum found, will be a lower bound of the initial criteria:

J (q, θ) ≥ J (q γ , θ) (9)
J (q γ , θ) is also known as the "Evidence Lower BOund" (ELBO) emphasizing the lower bound property on the evidence of the data.

A wise choice of the restriction on the variational distribution leads a feasible computation of the criterion. We choose to consider the following posterior shapes on the latent variables:

∀i Y i X (o) ∼ qγ M 1; τ (Y ) i ∀j Z j X (o) ∼ qγ M 1; τ (Y ) j ∀i A i X (o) ∼ qγ N ν (A) i , ρ (A) i ∀i B i X (o) ∼ qγ N ν (B) i , ρ (B) i ∀j C j X (o) ∼ qγ N ν (C) j , ρ (C) j ∀j D j X (o) ∼ qγ N ν (D) j , ρ (D) j .
We also impose the conditional independence of the latent variables to get a feasible computation of the entropy and of the negative "energy" (Equation 8) under q γ . This conditional independence is widely known as the "mean field Algorithm 1: Variational Expectation Maximization algorithm.

Input: observed data X (o) , k 1 and k 2 number of row groups and column groups ; -Initialize γ (0) and θ (0) ; while not convergence of criterion J do VE-step: find the variational parameters γ (t+1) that optimize

J γ, θ (t) γ (t+1) = arg max γ J γ, θ (t)
M-step: find the model parameters θ (t+1) that optimize J γ (t) , θ :

θ (t+1) = arg max θ J γ (t) , θ end
Result: θ and γ: model and variational parameters approximation" [START_REF] Parisi | Statistical field theory[END_REF]. We finally get the following fully factorized shape:

q γ = ∏ n1 i=1 M 1; τ (Y ) i × ∏ n2 j=1 M 1; τ (Z) j × ∏ n1 i=1 N ν (A) i , ρ (A) i × ∏ n1 i=1 N ν (B) i , ρ (B) i × ∏ n2 j=1 N ν (C) j , ρ (C) j × ∏ n2 j=1 N ν (D) j , ρ (D) j
, D) , ρ (D) ) denotes the parameters concatenation of the restricted variational distribution q γ .

where γ = (τ (Y ) , τ (Z) , ν (A) , ρ (A) , ν (B) , ρ (B) , ν (C) , ρ (C) , ν (
The new criteria J (γ, θ) that we want to optimize from now on is:

J (γ, θ) = H(q γ ) + E qγ log p X (o) , Y , Z, A, B, C, D; θ (10)
and the initial estimates of the model parameters θ are inferred as:

θ = arg max θ max γ J (γ, θ) . ( 11
)
This double maximization is realized with an iterative strategy and can be seen as an extension of the EM algorithm. The two steps are described in Algorithm 1.

Computation of the variational criterion

The restriction on the space of the variational distribution simplifies the computation of H(q γ ) as entropy is additive across independent variables:

H(q γ ) = - iq τ (Y ) iq log τ (Y ) iq - jl τ (Z) jl log τ (Z) jl + 1 2 i log 2πeρ (A) i + 1 2 i log 2πeρ (B) i + 1 2 j log 2πeρ (C) j + 1 2 j log 2πeρ (D) j
.

The independence of latent variables allows to rewrite the expectation of the complete log-likelihood as:

E qγ log p X (o) , Y , Z, A, B, C, D = E qγ [log p(Y )] + E qγ [log p(Z)] + E qγ [log p(A)] + E qγ [log p(B)
]

+ E qγ [log p(C)] + E qγ [log p(D)] + E qγ log p X (o) Y , Z, A, B, C, D . (12)
Despite the variational approximation, the expectation of the complete loglikelihood (12) can not be exactly computed as its last term involves an expectation under q γ of nonlinear functions:

E qγ log p X (o) Y , Z, A, B, C, D = ijql∶X (o) ij =0 τ (Y ) iq τ (Z) jl E qγ [log (p 0 )] + ijql∶X (o) ij =1 τ (Y ) iq τ (Z) jl E qγ [log (p 1 )] + ijql∶X (o) ij =NA τ (Y ) iq τ (Z) jl E qγ [log (1 -p 0 -p 1 )] , (13)
with p 0 and p 1 defined in Equations ( 5)-( 6).

These expectations can be approximated by the delta method [Wasserman, 2004, p. 79]. Using a first order Taylor expansion would lead to a criterion without maximum, so we use a second order Taylor expansion. The full expression of the criterion is given in Appendix A.

Maximization of the variational criterion

The VEM Algorithm 1 alternates between a maximization with respect to the variational parameters γ and a maximization w.r.t the model parameters θ. For our model, there is no explicit solution for the two maximizations of the criterion J (γ, θ), which are are carried out by the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm. We used automatic differentiation to compute the gradients needed for L-BFGS and for the Taylor series used in the variational criterion. We chose the Autograd library from HIPS and the submodule Autograd from PyTorch [Paszke et al., 2019]. These libraries rely on a reverse accumulation computational graph to compute exact gradients. Their high efficiency, even with large graphs, thanks to GPU acceleration, makes them particularly well adapted for the VEM algorithm.

Initialization

VEM does not ensure convergence towards a global optimum. The EM-like algorithms are known to be sensitive to the initialization, particularly when applied to models with discrete latent space, and may get stuck into unsatisfactory local maxima [Biernacki et al., 2003, Baudry andCeleux, 2015].

A simple solution consists in training for a few iterations from several random initializations, and pursue optimization with the solution with highest value of the variational criterion [see, e.g., small EM for mixtures [START_REF] Baudry | EM for mixtures[END_REF]. This exploration strategy spends a great deal of computing resources to bring out only a few good estimates. Another solution is to rely on simpler clustering methods, such as k-means or spectral clustering, to initialize the algorithm [START_REF] Shireman | Examining the effect of initialization strategies on the performance of Gaussian mixture modeling[END_REF].

For the Stochastic Block Model, a close relative of the Latent Block Model for graphs, [START_REF] Rohe | Spectral clustering and the highdimensional stochastic blockmodel[END_REF] prove the consistency of spectral clustering to identify the parameters of the Stochastic Block Model. Following this idea, we use a double spectral clustering (with absolute eigenvalues of the Laplacian as [START_REF] Rohe | Spectral clustering and the highdimensional stochastic blockmodel[END_REF]) on rows and columns on similarity matrices, to initialize our algorithm. Although this method is not designed for MNAR data, it can be expected to provide a satisfying initialization of the Latent Block Model if the missingness is not predominant. The parameters of our missingness model can not be initialized with this procedure; they are randomly initialized. The overall initialization procedure is described in Appendix B.

Model selection

Integrated Completed Likelihood criterion (ICL)

ICL, inspired by the Bayesian Information Criterion, was originally proposed to select a relevant number of classes for mixture models [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated classification likelihood[END_REF]]. It was extended to select an appropriate number of (row and column) clusters in the standard Latent Block Model [START_REF] Keribin | Model selection for the binary latent block model[END_REF]: for k 1 row classes and k 2 column classes, the criterion reads

ICL(k 1 , k 2 ) = log p(X, Y , Z) = log p( X, Y , Z θ; k 1 , k 2 )p(θ; k 1 , k 2 )dθ , (14)
with p(θ; k 1 , k 2 ) the prior distribution of parameters. By taking into account the latent variables Y , Z, ICL is a clustering-oriented criterion , whereas BIC or AIC are driven by the faithfulness to the distribution of X [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated classification likelihood[END_REF]].

For the LBM with MNAR missingness, ICL requires priors on the parameters of the missingness model. We chose independent and non informative InverseGamma(β, β) distribution (where β tends to zero) for the parameters

σ 2 A , σ 2 B , σ 2 C and σ 2 D .
As in [START_REF] Keribin | Model selection for the binary latent block model[END_REF], we use non-informative Dirichlet distribution priors on the parameters α and β of mixing proportions of classes. ICL reads

ICL(k 1 , k 2 ) = log p(X, Y , Z, A, B, C, D) (15) = log p( X, Y , Z, A, B, C, D θ; k 1 , k 2 )p(θ; k 1 , k 2 )dθ
Proposition 1. The ICL criterion for the LBM extended to the MNAR missingness process presented in Section 3.2 has the following asymptotic form for a data matrix of size n 1 × n 2 :

ICL(k 1 , k 2 ) = max θ log p X (o) , Y , Z, A, B, C, D; θ - k 1 k 2 2 log (n 1 n 2 ) - k 1 -1 2 log (n 1 ) - k 2 -1 2 log (n 2 ) + n 1 log (2π) -log (n 1 ) + n 2 log (2π) -log (n 2 ) + o(log n 1 ) + o(log n 2 ) .
See proof in Appendix C.

Since the maximum of the complete log-likelihood required to calculate the ICL is not available, in practice it is replaced by the lower bound provided by the variational approximation (see equation 10). An ICL criterion for the LBM with MAR missing data can be constructed in the same way, allowing for comparison with the MNAR model (see details in Appendix C).

Experiments on simulated data

Simulated data brings all the elements to assess clustering algorithms in controlled settings. Using controlled datasets provides the means to properly test the ability of an algorithm to recover the known underlying structure.

Assessing the difficulty of a co-clustering task

In co-clustering, several loss functions are suited for measuring the discrepancy between the underlying classes (Y , Z) and some predictions ( Ŷ , Ẑ). For our experiments, we will use the measure defined by [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF], that is, the ratio of misclassified entries in the data matrix:

l item (Y , Z, Ŷ , Ẑ) = 1-1 n 1 ∑ i δ Y i , Ŷi l row (Y , Ŷ ) + 1-1 n 2 ∑ j δ Z j , Ẑj l col (Z, Ẑ) -l row (Y , Ŷ ) l col (Z, Ẑ)
where δ is the Kronecker delta.

In standard clustering, the difficulty of a task is often assessed by its Bayes risk, that is, by the minimum of the expectation of the loss function, which is typically approximated by Monte Carlo on simulated data. Co-clustering poses specific difficulties. Adding more rows or more columns alter its difficulty because the dimensions of the spaces where the clustering is performed are expanded. The duality between the rows and the columns imply that the size of the matrix is a characteristic of a co-clustering problem. In other words, given a fixed generative distribution, as the matrix size increases, the difficulty of the task decreases, in contrast to simple clustering, where the difficulty, as measured by the Bayes risk, remains constant when more examples (that is, rows) are added.

A simple Monte Carlo approximation of the risk consists in averaging over many statistical units. In simple clustering, this means generating a great number of rows in a data matrix. In co-clustering, the statistical unit is the whole matrix, implying that a Monte Carlo approximation of the risk is obtained by generating a great number of data matrices; which then involves a great computational time. Furthermore, estimating the Bayes risk from a single data matrix is very inconstant; the risk may be very different between two data matrices of same size generated from the same distribution. Hence the usual notion of Bayes risk is not appropriate for co-clustering. [START_REF] Lomet | Design of artificial data tables for co-clustering analysis[END_REF] argue that conditioning the Bayes risk on the observed matrix is more appropriate. They give a protocol to simulate data matrices in which the difficulty of the clustering task is controlled by the following conditional Bayes risk:

r item ( Ŷ , Ẑ) = E l item (Y , Z, Ŷ , Ẑ) X (o) , (16)
where the expectation is taken over Y , Z only and Ŷ , Ẑ are the clusterings returned by the conditional Bayes classifier, that is, the maximum a posteriori : [START_REF] Lomet | Design of artificial data tables for co-clustering analysis[END_REF] released data sets, with different sizes and difficulties, simulated from the Latent Block Model. Using their protocol, we generated new data according the LBM with a MNAR missingness process. Data sets are generated according to the LBM with three row and column classes, with parameters

Ŷ , Ẑ = arg min (Y ,Z) r item (Y , Z) = arg max (Y ,Z) ij p Y i , Z j X (o) .
α = β = ⎛ ⎜ ⎝ 1 3 1 3 1 3 ⎞ ⎟ ⎠ and π = ⎛ ⎜ ⎝ 1 - 1 - 1 - 1 - 1 - ⎞ ⎟ ⎠ , ( 17
)
where defines the difficulty of the clustering task. The parameters of the MNAR process are

µ = 1, σ 2 A = 1, σ 2 B = 1, σ 2 C = 1, σ 2 D = 1 , (18)
which gives an average proportion of 35% of missing values. : Classification error with respect to the size of the data matrix (lower is better); ★ is the median of the conditional Bayes risk; ▲ is the median prediction error obtained by our algorithm.

Analyzing the classification of the proposed inference

We test here the ability of the proposed inference scheme to recover row and column classes. To conduct the experiments, we generate an initial data matrix of size n 1 = n 2 = 500 with a conditional Bayes risk of 5% set by choosing (17) by trial and error. The size of this matrix is then progressively reduced, removing rows and columns, to increase the difficulty of the classification task. The conditional Bayes risk is re-estimated on each sub matrix to provide a reference. Our algorithm is then run on these data matrices using 20 initializations for each run, as described in Section 4.4. We then predict the row and column classes (Y , Z) with their maximum a posteriori estimators on the variational distribution. This whole process is repeated 20 times, leading to the results presented in Figure 6.

As expected, the conditional Bayes risk decreases as the data matrices grow. The predictions returned by our algorithm follow the same pattern, with a diminishing gap to the conditional Bayes risk as the data matrices grow, which is consistent with our expectations. Appendix D provides additional experimental results that show consistent estimations of the model parameters.

Analyzing the benefit of a MNAR model versus a MAR model for MNAR data

The importance of using the right missingness model is tested by comparing the classifications returned by an LBM with and without an MNAR model. A data set is generated according to the LBM with MNAR values where the parameters α, β and π of the LBM are fixed as in (17), and is chosen in order to get a conditional Bayes risk of 12%, for data matrices of size n 1 = n 2 = 100; the MNAR model parameters µ, σ 2

A and σ 2 C all set to one which gives an average proportion of 35% of missing values. Several data matrices are generated using these parameters while varying the value of the σ 2 B and σ 2 D parameters that govern the MNAR effects; these variations do not affect the conditional Bayes risk nor the proportion of missing values. For each data matrix, we train the LBM with either the MAR or the MNAR model. This process is repeated 20 times, starting from the generation of a new fully observed data matrix.

The median of the classification errors l item are presented in Figure 7 as a function of the MNAR effect. They are essentially constant and close to the conditional Bayes risk for the LBM with the MNAR model, whereas the LBM with the MAR model is badly affected by MNAR data, eventually leading to a classification close to a totally random allocation 1 . Ignoring the nature of the missingness process leads here to strong biases in estimation that in turn drastically affect classification. Thankfully, the ICL criterion may be of great help to select the right missingness model as shown in Section 6.5.

Analyzing the ability of the model selection criterion to select the adequate number of classes

We reuse the parameters ( 17) and ( 18) to analyze the behavior of the ICL criterion. We consider different sizes of data matrices, between (30,30) and (150,150), with varying difficulty for each matrix size, with a conditional Bayes risk ( 16) of respectively 5%, 12% and 20%

The results in Figure 8 show that, as expected, the ICL criterion tends to select more often the right number of classes as the data matrices get larger and also when classes are more separated. We also observe that the ICL criterion tends to be conservative for small data matrices, by underestimating the number 1 With equal class proportions, the expected classification error of a random allocation is

k 1 -1 k 1 + k 2 -1 k 2 -k 1 -1 k 1 k 2 -1 k 2
, that is, 0.89 here where

k 1 = k 2 = 3. r item ( Ŷ , Ẑ) = 5% r item ( Ŷ , Ẑ) = 12% r item ( Ŷ , Ẑ) = 20% k 2 k 2 k 2 2 3 4 5 2 3 4 5 2 3 4 5 n 1 = n 2 = 30 k 1 2 9 2 13 3 14 2 3 1 7 2 2 2 1 4 5 1 1 n 1 = n 2 = 40 k 1 2 4 2 17 1 17 1 3 14 1 1 1 1 4 5 n 1 = n 2 = 50 k 1 2 1 11 2 15 1 2 3 2 17 7 1 1 4 5 n 1 = n 2 = 75 k 1 2 1 9 13 2 3 1 17 1 9 4 4 1 5 1 1 n 1 = n 2 = 100 k 1 2 2 11 3 19 18 1 7 1 4 5 1 n 1 = n 2 = 150 k 1 2 1 1 5 1 3 14 2 1 18 1 11 4 1 5 2 2
Figure 8: Count number of (k 1 , k 2 ) models selected by the ICL criterion among 20 data matrices for different difficulties, as measured by the conditional Bayes risk, and different matrix sizes. All matrices are generated with the same number of row and column classes:

k 1 = k 2 = 3.
of classes. It could come to the fact that the size of the matrix is not large enough to consider the asymptotic approximation as valid and/or it could come from the approximations used to compute the log-likelihood J (variational restriction and delta method).

Analysing the ability of the model selection criterion to select the adequate missingness model

We use the models fitted in Section 6.3 to analyze the ability of the ICL criterion to select the right missingness model (MNAR or MAR). The difference in ICL between the MAR and MNAR models is computed for each data matrix, assuming that the right numbers of classes (k 1 , k 2 ) are known.

The results, presented in Figure 9, show that ICL rightfully opts for the MNAR model almost everywhere, demonstrating the ability of this criterion to select the adequate missingness model. The MAR model is only chosen for some experiments with the lowest MNAR effect (σ2 B = σ 2 D = 0.01), where the prediction performances are almost identical (see Figure 7), with a median difference in ICL of -0.51 (the MAR model is chosen 13 times over the 20 repetitions).

Experiments on real data

We consider voting records 2 from the lower house of the French Parliament (Assemblée Nationale). This dataset gathers the results of the 1256 ballots of year 2018 of the 577 French members of parliament (MPs) for the procedural motions and amendments for the 15th legislature (June 2017). For each text, the vote of each MP is recorded as a 4-level categorical response: "yes", 'no", "abstained" or "absent". Using our model, we bring out some relevant groups of texts and MPs, as well as some structure in the behavior of nonvoters.

We gather the data in a matrix where each row represents an MP and each column represents a text. To use our model, we reduced the 4 response levels to 3 ("yes", 'no", "missing") assuming that merging the "abstained" and "absent" categories would not affect much the underlying missingness process ("abstained" votes represent about 4% of the expressed votes, "missing" responses represent 85% of all votes).

At the lower house of French Parliament, MPs may group together according to their political affinities. Groups with less that 15 members or MPs who choose to be independent are gathered under the "Non inscrits" (NI) label, giving a heterogeneous range of political hues inside it. The names of the groups and their frequency are detailed in Figure 10. Fewer classes lead to over-aggregated components hiding the subtleties of the network, but since they still correspond to well-identified groups and are more friendly to visual analysis, we provide them as additional material in Appendix E.

In Figure 11, classes of MPs are coherent to their political orientation: class 0 and 1 are mainly made up of left-wing MPs from the groups SOC, FI, GDR, LT, classes 2 and 3 are mainly made up of right-wing MPs from LR and the classes from 6 to 13 are mainly made up of centrist MPs from LaREM and MODEM who are known to be political allies. Classes of texts can be analyzed with the available metadata. A bipartite opposition system appears from classes A and C. Texts from class A are the original articles of law proposed by the government and are unsurprisingly voted positively by the MPs classes from 6 to 13 as they are from the same political mould as the French government. Texts from class C are mainly amendments proposed by minority and are voted positively by both the left wing (class 0 and 1) and the right wing (classes 2 and 3) and negatively by the MPs supporting the government (classes 6 to 13). The left and right wings are yet divided by usual issues such as immigration regulation amendments gathered in classes G and M or general economic matters gathered in classes H and I.

In our model, the latent variables A and B characterize the propensity of MPs to cast a vote. Figure 12 displays the scatter plot of ν

(A) i and ν (B) i
, the maximum a posteriori estimates of A i and B i for all MPs under the variational distribution. The abscissa represents the propensity to vote3 , with higher values of ν (A) corresponding to a higher propensity to vote, and the ordinate ν (B) represents the additional effect of casting a vote when supporting the text. The membership of MPs to their political group is indicated by the plotting symbol.

We see two obvious clusters separated by the vertical axis ν (B) : the bottom cluster is essentially formed by MPs from the LaREM and MODEM political groups, which support the government, whereas the top cluster is formed by the opposition political groups. The ν (B) estimates for the opposition cluster are positive, meaning that these MPs come to parliament to vote positively. This behavior is not surprising because the MPs of the opposition parties are outnumbered by the MPs supporting the government, so they must be diligent if they want their tabled motion or amendment passed. The dependency between the political groups and the MNAR effect encoded in the estimates ν (B) , which is confirmed by an ANOVA test (with a p-value smaller than numerical error), supports that the missingness patterns captured by our model are relevant for the problem at hand. A similar analysis is developed on texts in Appendix E.

Conclusion

In many estimation problems, the absence of data conveys some information on the underlying phenomenon that should be exploited for its modeling. We propose a co-clustering model that accounts for this absence of data; it aims at retrieving groups of rows and columns based on the complete data matrix instead of considering only the partitioning of the observed data matrix. This model consists of two building blocks: a co-clustering model (Latent Block Model) of the full data matrix, and a missingness model that manages the censoring that produces the observed data matrix. This missingness model preserves the symmetry of the co-clustering model by allowing two MNAR effects, one on the We used variational techniques and the Delta method to obtain a tractable approximation of the lower bound of the observed log-likelihood. We proposed a model selection criterion to select both the number of classes and the type of missingness (MAR versus MNAR).

-3 -2 -1 0 1 2 3 4 ν A -6 -4 -2 0
Our experiments on synthetic datasets show that ignoring an informative missingness can lead to catastrophic co-clustering estimates, supporting the value of using expressive missingness models on such type of data. We also illustrate the use of our model on a real-world case where the missingness model provides an interesting basis for analyzing and interpreting the motivations of nonvoters.

Our model should also be useful in other fields such as in ecology, where the probability of observing interaction between species derives from some factors that also explain the true interactions [START_REF] Vázquez | Uniting pattern and process in plant-animal mutualistic networks: a review[END_REF], or in collaborative filtering, where the probability of observing a rating depends on the actual rating that would be given by the user [START_REF] Marlin | Collaborative filtering and the missing at random assumption[END_REF]. In the latter application, the data sizes generally encountered in recommendation would require computational improvements in inference. Another useful future work is to extend our model to non-binary data.

with the following terms:

E qγ [log p(Y )] = iq E qγ Y iq log α q = iq τ (Y ) iq log α q E qγ [log p(Z)] = jl E qγ Z jl log β l = jl τ (Z) jl log β l E qγ [log p(A)] = - n 1 2 log 2π - n 1 2 log σ 2 A - 1 2σ 2 A i E qγ A 2 i = - n 1 2 log 2π - n 1 2 log σ 2 A - 1 2σ 2 A i ν (A) i 2 + ρ (A) i E qγ [log p(B)] = - n 1 2 log 2π - n 1 2 log σ 2 B - 1 2σ 2 B i ν (B) i 2 + ρ (B) i E qγ [log p(C)] = - n 2 2 log 2π - n 2 2 log σ 2 C - 1 2σ 2 C j ν (C) j 2 + ρ (C) j E qγ [log p(D)] = - n 2 2 log 2π - n 2 2 log σ 2 D - 1 2σ 2 D j ν (D) j 2 + ρ (D) j and E qγ log p X (o) Y , Z, A, B, C, D = ql,ij∶X (o) ij =1 τ (Y ) iq τ (Z) jl E qγ [log p 1 ] + ql,ij∶X (o) ij =0 τ (Y ) iq τ (Z) jl E qγ [log p 0 ] + ql,ij∶X (o) ij =NA τ (Y ) iq τ (Z) jl E qγ [log (1 -p 0 -p 1 )] , (20)
with p 0 and p 1 defined in Equations ( 5)-( 6).

Equation ( 20) involves the computation of the expectations of the following nonlinear functions:

f 1 (x, y) = log (π ql expit (µ + x + y)) f 0 (x, y) = log ((1 -π ql ) expit (µ + x -y)) f NA (x, y) = log (1 -π ql expit (µ + x + y) -(1 -π ql ) expit (µ + x -y)) .
The approximation of these expectations given by the second-order Delta method with independent random variables X and Y reads:

E[f (X, Y )] ≈ f (EX, EY ) + 1 2 var(X) ∂ 2 f (E[X], E[Y ]) ∂(X) 2 + 1 2 var(Y ) ∂ 2 f (E[X], E[Y ]) ∂(Y ) 2 ,
B Initialization of the VEM algorithm with spectral clustering.

Algorithm 2: Initialization of the VEM algorithm with spectral clustering. Input:

• observed data X (o)

• k 1 and k 2 number of row groups and column groups Function SpectralClustering(W adjacency matrix, k number of clusters):

• Define D the diagonal matrix, element of R n×n :

D ii = ∑ q W iq • Define L = D -1 2 W D -1 2
• Find the eigenvectors corresponding to the k eigenvalues of L that are largest in absolute value. Form the matrix U = [U 1 , ..., U k ] ∈ R n×k concatenating the eigenvectors into columns.

Return results of k -means with k clusters on U . begin

• Build Y the n 1 × k 1 indicator matrix of the row cluster memberships with results of SpectralClustering(XX T , k 1 )

• Build Z the n 2 × k 2 indicator matrix of the column cluster memberships with results of SpectralClustering(X T X, k 2 ).

• α, β and π estimated from Y and Z

• µ initialized such as expit(µ) is the global missingness rate

• σ 2 A , σ 2 B , σ 2 C and σ 2 D sampled from U ]0,1]
end Result: The ICL criterion of the LBM extended to the MNAR missingness process presented in Section 3.2 has the following asymptotic form for n 1 and n 2 :

• θ = α, β, π, µ, σ 2 A , σ 2 B , σ 2 C , σ 2
ICL(k 1 , k 2 ) =max θ log p X (o) , Y , Z, A, B, C, D; θ - k 1 k 2 2 log (n 1 n 2 ) - k 1 -1 2 log (n 1 ) - k 2 -1 2 log (n 2 ) (21) + n 1 log (2π) -log (n 1 ) + n 2 log (2π) -log (n 2 ) + o(log n 1 ) + o(log n 2 ) .
Proof. With independent latent variables and independent priors on the parameters, the ICL criterion reads

ICL = log p X (o) , Y , Z, A, B, C, D = log p X (o) , Y , Z, A, B, C, D θ p(θ)dθ = log p X (o) Y , Z, A, B, C, D, π p(π)p(µ)dπdµ (22)
+ log p( Y α )p(α)dα + log p( Z β )p(β)dβ + log p A σ 2 A p σ 2 A dσ 2 A + log p B σ 2 B p σ 2 B dσ 2 B + log p C σ 2 C p σ 2 C dσ 2 C + log p D σ 2 D p σ 2 D dσ 2 D .
As in the ICL developed by [START_REF] Keribin | Model selection for the binary latent block model[END_REF] for the standard LBM, we set non-informative Dirichlet distribution D(a, ..., a) priors on α and β:

log p(Y ) = log p( Y α )p(α; a)dα = i log q (α q ) Yiq 1 B(a) iq (α q ) a-1 dα = log B(a + i Y i ) -log B(a) = q log Γ(Y ∶q + a) + log Γ(k 1 a) -log Γ(n 1 + k 1 a) -k 1 log Γ(a) , where Y ∶q = ∑ i Y iq . The Stirling approximation log Γ(x) = x log x -x -1 2 log x +
o(log x) leads to the following asymptotic development of log p(Y ):

log p(Y ) = q log Γ(Y ∶q + a) -log Γ(n 1 + k 1 a) + o(log n 1 ) = q Y ∶q log Y ∶q -n 1 - 1 2 n 1 -n 1 log n 1 + k 1 a log n 1 -n 1 - 1 2 log n 1 + o(log n 1 ) .
With the non-informative Jeffrey prior a = 1 2 , this gives:

log p(Y ) = q Y ∶q log( 1 n 1 Y ∶q ) - k 1 -1 2 log n 1 + o(log n 1 ) = max α log p(Y ; α) - k 1 -1 2 log n 1 + o(log n 1 ) . (23)
Similarly, we get:

log p(Z) = l log Γ(Z ∶l + a) + log Γ(k 2 a) -log Γ(n 2 + k 2 a) -k 2 log Γ(a) = max β log p(Z; β) - k 2 -1 2 log n 2 + o(log n 2 ) , (24)
where Z ∶l = ∑ j Z jl . We set non-informative InverseGamma(β, β) distributions (as β tends to zero) as priors on σ 2 A , σ 2 B , σ 2 C and σ 2 D :

log p(A) = p A σ 2 A p σ 2 A ; β dσ 2 A = i log 2σ 2 A - n 1 2 exp - ∑ A 2 i 2σ 2 A β β Γ(β) exp - β σ 2 A σ 2 A -β-1 dσ 2 A = log β β Γ(β) 2 -n 1 2 σ 2 A - n 1 2 -β-1 exp - 2β + ∑ A 2 i 2 ⋅ 1 σ 2 A dσ 2 A = log β β Γ(β) 2 β 2β + A i - n 1 2 -β Γ( n 1 2 + β) .
To consider a non-informative InverseGamma(β, β) distribution, we realize a first order Taylor development as β tends to 0:

log p(A) ≈ log Γ( n 1 2 ) + log β - n 1 2 log( A 2 i ) .
Using the Stirling approximation of Γ(x) we get the following asymptotic development of log p(A):

log p(A) = n 1 2 log n 1 - n 1 2 - 1 2 log n 1 - n 1 2 log A 2 i + o(

C.2 ICL of the LBM with MAR data

We consider the following LBM extended with the MAR missingness process:

Latent Block Model

Y i iid ∼ M(1; α), α ∈ S k1-1 Z j iid ∼ M(1; β), β ∈ S k2-1 ( X ij Y i = q, Z j = l ) ind ∼ B(π ql ), π ql ∈ [0, 1] MAR data model A i iid ∼ N 0, σ 2 A , σ 2 A ∈ R * + C j iid ∼ N 0, σ 2 C , σ 2 C ∈ R * + ( M ij A i , C j )
ind ∼ B(expit (µ + A i + C j )) Observations are generated according to:

X (o) ij X ij , M ij = X ij if M ij = 1 NA if M ij = 0
The ICL of this model has the following asymptotic form for n 1 and n 2 :

ICL(k 1 , k 2 ) = max θ log p X (o) , Y , Z, A, C; θ - k 1 k 2 2 log (n 1 n 2 ) - k 1 -1 2 log (n 1 ) - k 2 -1 2 log (n 2 ) (28)
+ 1 2 (n 1 log (2π) -log (n 1 ) + n 2 log (2π) -log (n 2 ))

+ o(log n 1 ) + o(log n 2 ) .

D Supplemental figures for estimations

This section provides additional experimental results that show a consistent estimation of the model parameters. We reuse the data matrices generated by the LBM with missing data from Section 6.2. An initial data matrix of size n 1 = n 2 = 500 with a conditional Bayes risk of 5% was generated and progressively reduced, removing rows and columns, to increase the difficulty of the classification task.

Figure 13 displays the maximum absolute error made on the parameters π of the Bernoulli distributions that model the probability of X conditionally to the row and column classes. This error decreases as the size of the data matrices grows, which is consistent with our expectations. belonging to six text classes. This reflects the fact that our model does not link the MNAR effect to the LBM memberships.

Figure 1 :

 1 Figure 1: Summary of the standard Latent Block Model with binary data.

Figure 2 :

 2 Figure 2: Graphical representation of the MCAR model. The partially observed entry X (o)

Figure 3 :

 3 Figure 3: Graphical representation of the MAR model. The partially observed entry X (o)

Figure 4 :

 4 Figure 4: Graphical representation of the MNAR model. The partially observed entry X (o)

Figure 5 :

 5 Figure 5: Graphical view and summary of the Latent Block Model extended to MNAR missingness process. The observed data X (o) ij is generated by the necessary information carried by the class and propensity of row i and by the class and propensity of the column j.

  Figure6: Classification error with respect to the size of the data matrix (lower is better); ★ is the median of the conditional Bayes risk; ▲ is the median prediction error obtained by our algorithm.

Figure 9 :

 9 Figure9: Difference in ICL between the MAR and MNAR models with respect to an increase of the MNAR effect, ★ is the median. The MNAR model is selected when the difference in ICL is positive.

Figure 10 :

 10 Figure 10: Hemicycle of the political groups of the French National Assembly

Figure 11 :

 11 Figure 11: Left: matrix of votes reordered according to the row and column classes, for the MNAR LBM model selected by ICL, with 14 MP classes and 14 text classes. The red lines delineate class boundaries. The counts of MPs belonging to their political groups in each MP class is given on the left. Right: summary of the inferred opinions (expressed or not) for all classes of texts and MPs, as given by the estimated probability to support a text in each block of the reordered matrix.

Figure 12 :

 12 Figure 12: Left: maximum a posteriori estimates of the MPs propensities (ν (A) i , ν (B) i ), with their political group memberships. ν (A) i drives the MAR effect and ν (B) i drives the MNAR one. Right: ICL curve. Maximum is reached for k 1 =14 and k 2 =14

D

  the model parameters • Y and Z the row and column cluster memberships C Asymptotic form of the Integrated Completed Likelihood C.1 ICL of the MNAR model

Figure 13 :

 13 Figure13: Maximum error between the true (π) and the estimated ( π) probabilities associated to the blocks of the data matrix X as a function of its size.

Figure 14 displaysFigure 14 :

 1414 Figure 14 displays the mean squared error (MSE) between the generated and estimated values of the latent variables A, B, C, D responsible for the individual variability of missingness. The estimated values are given by the maximum a posteriori of their corresponding variational distribution. The MSE curves of the variables A and C are comparable as well as the curves of the variables B and D. This is expected as the data matrices are generated with symmetric characters in rows and columns.

Figure 15

 15 Figure 15 compares the estimated values of A, B, C and D to their true generated values for two different sizes of data matrices, all other parameters being equal. A linear trend is exhibited from these scatter plots showing a good aptitude of the proposed inference to recover extreme negative and positives values.

Figure 17 :

 17 Figure 17: maximum a posteriori estimates of the text propensities (ν (C) j , ν (D) j),

  + o(log n 2 ) .Using the standard BIC approximation, we havelog p X (o) Y , Z = log p X (o) Y , Z, π p(π)p(µ)dπdµ + o(log n 1 ) + o(log n 2 ) .The ICL criterion 21 is directly derived from equations 22, 23, 24, 25, 26 and 27.

	Similarly, we get:
	log p(B) = max σ 2 B	log p B; σ 2 B +	n 1 2	log(2π) -	1 2	log n 1 + o(log n 1 )
	log p(C) = max σ 2 C	log p C; σ 2 C +	n 2 2	log(2π) -	1 2	log n 2 + o(log n 2 )	(26)
	log p(D) = max σ 2 D	log p D; σ 2 D + log n 2 = max n 2 2 log(2π) -1 2 k 1 k 2 2	log(n 1 n 2 )	(27)
								log n 1 )
	= max σ 2 A	log p A; σ 2 A +	n 1 2	log(2π) -	1 2	log n 1 + o(log n 1 ) .	(25)

π log p X (o) Y , Z; π, µ +

Votes from the French National Assembly are available from http://data. assemblee-nationale.fr/travaux-parlementaires/votes.

More rigorously, the abscissa represents the global deviation from the average propensity to vote.

A Computing the criterion J (q γ , θ)

The criterion to be optimized is : J (q γ , θ) = H(q γ ) + E qγ [log p(X (o) , Y , Z, A, B, C, D; θ)] ,

(

where θ is the list of all model parameters: θ = (α, β, π, µ, σ 2 A , σ 2 B , σ 2 C , σ 2 D ). We restrict the form of the variational distribution q γ to get a fully factorized form:

, where γ denotes the list of parameters of the variational distribution:

). The entropy is additive across independant variables, so we get:

.

The independence of the latent variables allows to rewrite the expectation of the complete log-likelihood as:

which yields in our case:

The criterion is now fully computable.

E Supplemental figures for the French national assembly votes analysis

Figure 16 displays the reordered matrix of votes derived from a block clustering with a small number of classes. Such a simplification may be helpful for identifying global trends. With this model, the three MP classes are broadly identified as gathering the right-wing (first class) and left-wing (second class) opposition parties, the last class being formed of the political groups supporting the government. The opposition systems appear clearly: on the texts from classes A and E, the votes contrast the membership to the opposition parties versus the governmental alliance, whereas on texts from classes C and D, they separate the left-wing from the right-wing oppositions. Class B gathers various texts on topics of rather general agreement pertaining to social or health matters.

Going back to the model selected by ICL described in Section 7, we analyze the text propensities to be voted upon and to be positively perceived by nonvoters. These propensities are encoded in the values of the latent variables C and D. , the maximum a posteriori estimates of C j and D j under the variational distribution, for all voted texts. The abscissa ν (C) reflects the mobilization on the texts, with higher mobilization for higher values, and the ordinate ν (D) represents the additional effect of mobilizing specifically supporting voters. The fourteen-cluster membership of texts (there is no obvious relevant classification for texts) is indicated by the plotting symbol.

Some relationship between missingness and membership to text classes emerge from this plot. A first cluster of text appears in the positive quadrant, with texts mainly proposed by the government, categorized in text classes A and B. A second cluster, smaller, on the upper left, is mainly formed by texts categorized in class D, voted positively by few voters. All the texts are related to the same law project regarding housing and were voted over a short period (06/03/2018 and 06/08/2018). The largest cluster, on the lower part of the graph, gathers most of the remaining texts, that would have a tendency to be voted negatively by nonvoters. These texts were proposed by either the right-wing or left-wing opposition, and get little support from a vast majority of MPs. Note also that the small group of highly voted texts, on the right-hand side, is made of texts