Aude Sportisse

Matthieu Marbac

Fabien Laporte

Gilles Celeux

Claire Boyer

Julie Josse

Christophe Biernacki

Model-based Clustering with Missing Not At Random Data

Keywords: Model-based Clustering, Informative Missing Values, EM and Stochastic EM Algorithms, Medical Data

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Clustering remains a crucial tool for the comprehensive analysis of large datasets, providing a concise summary through the grouping of observations. Notably, the model-based paradigm [START_REF] Mclachlan | Mixture Models: Inference and Applications to Clustering[END_REF][START_REF] Bouveyron | Model-based Clustering and Classification for Data Science: with Applications in R[END_REF] facilitates clustering by yielding interpretable models that enhance our understanding of the relationships between the formed clusters and the features at play. This parametric framework exhibits flexibility in handling highdimensionality problems [START_REF] Bouveyron | High-dimensional data clustering[END_REF][START_REF] Bouveyron | Model-based clustering of high-dimensional data: A review[END_REF], mixed datasets [START_REF] Marbac | Model-based clustering of gaussian copulas for mixed data[END_REF], and even time series and dependent data [START_REF] Ramoni | Bayesian clustering by dynamics[END_REF][START_REF] Xiong | Time series clustering with arma mixtures[END_REF]. However, the challenge accompanying this multifaceted model-based clustering lies in the requisite modeling efforts to design mixture models tailored to the data structure.

In large-scale data analysis, the problem of missing data is ubiquitous, data collection being never perfect (e.g. machines which fail, non-responses in a study). Classical approaches for dealing with missing data consist of working on a complete dataset [START_REF] Little | Statistical Analysis with Missing Data[END_REF], either by using only complete individuals, or by imputing missing values. However, both methods can cause huge problems in the analysis, either by reducing too drastically the dataset to a possibly biased subsample, or by distorting the distribution of the completed samples, respectively. Furthermore, it's crucial to note that both of these strategies are essentially preprocessing steps and are not specifically tailored for the final clustering task. Alternatively, one can explore likelihood-based approaches, employing methods like Expectation Maximization (EM) algorithms [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. In this paper, we adopt such an approach to enable model-based clustering to effectively handle informative missing data in an efficient manner.

We assume that the missing data are missing not at random (MNAR) [START_REF] Rubin | Inference and missing data[END_REF][START_REF] Ibrahim | Missing responses in generalised linear mixed models when the missing data mechanism is nonignorable[END_REF][START_REF] Mohan | Estimation with incomplete data: The linear case[END_REF]. More specifically, we consider that the cause of the missingness can be explained by the membership to a class, which is not observed, and we call this specific MNAR model MNARz. As the MNAR mechanism is neither ignorable for the density estimation (parameters estimation), nor for the clustering (partition estimation), dealing with such data does require the specific modeling effort for the distribution of the missingdata pattern, indicating where are the missing values in the data. An example of MNAR data includes clinical data collected in emergency situations, where doctors may choose to treat patients suffering from severe trauma before taking measurements. This intrinsically leads to more missing data in this class.

Related works on clustering despite missing values

In order to handle missing values in a model-based clustering framework, Hunt and Jorgensen [START_REF] Hunt | Mixture model clustering for mixed data with missing information[END_REF] implement the standard EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] based on the observed likelihood and Serafini et al. [START_REF] Serafini | Handling missing data in model-based clustering[END_REF] propose to perform multiple imputations (with Monte Carlo methods) in the E-step. However, both works only consider MCAR data, when the cause of the missingness is completely independent from the data values.

Different clustering methods have been developed to deal with MNAR mechanisms. In a partition-based framework, Chi et al. [START_REF] Chi | k-pod: A method for k-means clustering of missing data[END_REF] propose an extension of k-means clustering for missing data, called k-Pod, without requiring the missing-data pattern to be modelled. However, like k-means clustering, the k-Pod algorithm relies on strong assumptions as equal proportions between clusters. Du Roy De Chaumaray and Marbac [START_REF] Roy De Chaumaray | Clustering data with nonignorable missingness using semi-parametric mixture models[END_REF] perform clustering via a semiparametric mixture model using the patternmixture approach to formulate the joint distribution, which makes the method not suitable for estimating the density parameters or imputing missing values. For longitudinal data, Beunckens et al. [START_REF] Beunckens | A latent-class mixture model for incomplete longitudinal gaussian data[END_REF] and Kuha et al. [START_REF] Kuha | Latent variable modelling with nonignorable item nonresponse: multigroup response propensity models for crossnational analysis[END_REF] jointly model the measurements and the dropout process by using an extension of the shared-parameter model.

Contributions

We first present the considered MNARz mechanism, in which the missingness depends on class membership, within the context of unsupervised classification based on mixture models for different types of data (continuous, count, categorical, and mixed). We then demonstrate that, under MNARz, statistical inference can be conducted on the augmented matrix formed from the concatenation of the data matrix and the missingdata pattern (binary mask indicating where are the missing values) by considering a missing at random (MAR) mechanism instead. In other words, the cause of the missingness can be explained by the observed variables, making it ignorable. This offers a significant advantage, as there is no need to model the missing-data mechanism in such cases. This also provides theoretical insights into an approach commonly used in practice without theoretical foundations, wherein working on the augmented data matrix under a MAR assumption is usually suggested to facilitate efficient learning, despite a more complex underlying missing mechanism [START_REF] Josse | On the consistency of supervised learning with missing values[END_REF]. We then propose an EM algorithm, which has been implemented and is available for reproducibility 1 . More general MNAR models are also discussed, for which the cause of the missingness can depend on the values of the missing variable themselves, as well as identifiability and estimation strategies. However, they are empirically shown to be well apprehended by the more simple MNARz model. Finally, we assess the numerical performance of our method on synthetic data and the real medical registry TraumaBase ® .

Missing-data in model-based clustering

Set the dataset Y = (y 1 | . . . |y n ) T consisting of n individuals, where each observation y i = (y i1 , . . . , y id ) T belongs to a space Y, depending on the type of data, defined by d features. The pattern of missing data is denoted by

C = (c 1 | . . . |c n ) T ∈ {0, 1} n×d , with c i = (c i1 , . . . , c id ) T ∈ {0, 1} d : c ij = 1 indicates
that the value y ij is missing and c ij = 0 otherwise. The values of the observed (resp. missing) variables for individual i are denoted by y obs i (resp. y mis i ). The objective of clustering is to estimate an unknown partition Z = (z 1 | . . . |z n ) T ∈ {0, 1} n×K that groups the full dataset Y into K classes, with z i = (z i1 , . . . , z iK ) T ∈ {0, 1} K and where z ik = 1 if y i belongs to cluster k, z ik = 0 otherwise. Consequently, in a clustering context, the missing data are not only the values y mis i but also the partition labels z i .

Mixture models

Mixture models allow for clustering by modeling the distribution of the observed data (y obs i , c i ). Assuming an underlying mixture model with K components, the probability distribution function (pdf) of the couple (y i , c i ) reads as

f (y i , c i ; θ) = K k=1 π k f k (y i ; λ k )f k (c i | y i ; ψ k ), (1)
where θ = (γ, ψ) gathers all the model parameters, γ = (π, λ) groups the parameters related to the marginal distribution of y i , π = (π 1 , . . . , π K ) is the vector of proportions with K k=1 π k = 1 and π k > 0 for all k ∈ {1, . . . , K}. Given λ = (λ 1 , . . . , λ K ), f k (• ; λ k ) is the pdf of the k-th component parameterized by λ k , ψ = (ψ 1 , . . . , ψ K ) groups the parameters of the missingness mechanisms and f k (c i | y i ; ψ k ) is the pdf related to the missingness mechanism under component k

(i.e., f k (c i | y i ; ψ k ) = f (c i | y i , z ik = 1; ψ k ))
. In many cases, the parameter ψ is interpreted as a nuisance parameter. However, when the mechanism is not ignorable, i.e. can not be ignored when performing inferences for λ, we need to consider the whole parameter θ to achieve clustering, since the pdf of the observed data is

f (y obs i , c i ; θ) = f (y i , c i ; θ)dy miss i . (2)
Different types of pdf f k (• ; λ k ) can be considered, depending on the types of features at hand. Thus, if y i is a vector of continuous variables, the pdf of a d-variate Gaussian distribution [START_REF] Mclachlan | Mixture Models: Inference and Applications to Clustering[END_REF][START_REF] Banfield | Model-based gaussian and non-gaussian clustering[END_REF] can be considered for f k (y i ; λ k ) and thus λ k groups the mean vector and the covariance matrix. Moreover, if some components of y i are discrete or categorical, the latent class model (see [START_REF] Geweke | Alternative computational approaches to inference in the multinomial probit model[END_REF][START_REF] Mcparland | Model based clustering for mixed data: clustmd[END_REF]) defining f k (y i ; λ k ) = d j=1 f kj (y ij ; λ kj ) can be used, with λ k = (λ k1 , . . . , λ kd ). In such a case, f kj could be the pdf of a Poisson (resp. multinomial) distribution with parameter λ kj if y ij is an integer (resp. categorical) variable. The choice of the modeling for the missingness mechanism (i.e., the distribution f k (c i | y i ; ψ k )) is discussed in the following.

To formulate the joint distribution of the data and the missing-data pattern, we consider in this paper the selection model [START_REF] Heckman | Sample selection bias as a specification error[END_REF], which factorizes it into the product of the marginal data density and the missing-data mechanism [START_REF] Mclachlan | Mixture Models: Inference and Applications to Clustering[END_REF]. This approach has the great advantage of allowing imputation of the missing values and density estimation throughout the parameter estimation of the mixture model. Another approach, called the pattern-mixture model [START_REF] Little | Pattern-mixture models for multivariate incomplete data[END_REF], can be used, involving the product of the marginal density of C and the conditional density of Y given C; it has been considered by [START_REF] Chaumaray | Clustering data with non-ignorable missingness using semi-parametric mixture models assuming independence within components[END_REF] for a clustering purpose.

The MNARz model

To handle MNAR data in selection models, the distribution of the missing-data pattern given the data and the partition should be specified. We consider that the elements of c i are conditionally independent given (y i , z i ). By the categorical nature of the mask c i , this independence assumption is a quite natural hypothesis in the context of clustering [START_REF] Chi | k-pod: A method for k-means clustering of missing data[END_REF][START_REF] Roy De Chaumaray | Clustering data with nonignorable missingness using semi-parametric mixture models[END_REF].

In the MNARz model, we consider that the only effect of missingness is on the class membership k, being the same for all variables. More specifically, the conditional distribution of c ij given (y i , z i ) is assumed to be a generalized linear model with link function ρ, so that finally

f k (c i | y i ; ψ k ) = d j=1 ρ(α k ) cij (1 -ρ(α k )) 1-cij , (3)
where ψ k = α k , in this case. The MNARz model is the simplest of the MNAR models we can propose (see Section 3.3 for more general ones). Roughly speaking, this model assumes that the proportion of missing values can vary among the clusters. Although MNARz does not directly involve y i in its ground definition (3), it does not mean that the pattern c i does not depend on y i since z i depends itself on y i ; see Figure 1 for an illustration.

Proposal

Reinterpretation of the MNARz model as a MAR strategy

Interestingly, the MNARz model can be turned into a MAR-like one by working on the augmented matrix formed from the concatenation of the data matrix and the missingdata pattern. This is the purpose of the next theorem, proven in Appendix A.

; θ) = K k=1 π k f k (y obs i ; λ k ) d j=1 ρ(α k ) cij (1 -ρ(α k )) 1-cij . ( 4
)
Then for fixed parameter θ = (π, λ, ψ), the mixture model for ỹobs i is the same as the distribution for y obs i with the mixture model (1) under the MNARz mechanism (3).

Theorem 1 implies that the maximum likelihood estimate of θ is the same considering ỹobs i under the MAR assumption and y obs i under the MNARz assumption (3). This implies that if the mechanism is MNARz, an (EM) algorithm designed for MAR data can be used on the augmented data set instead, capitalizing on efficient implementations dedicated to such a well-studied setting (see Section 4). In fact, Theorem 1 is the first theoretical result in unsupervised learning in line with the intuition developed in [START_REF] Josse | On the consistency of supervised learning with missing values[END_REF] for supervised learning and in [START_REF] Sportisse | Imputation and low-rank estimation with missing not at random data[END_REF] for estimation in low-rank models, that working with MAR strategies on the data set augmented by the missing pattern can actually tackle certain types of MNAR settings.

Furthermore, the identifiability of the model parameters when considering MNARz data follows directly from this reinterpretation as a MAR model. Indeed, identifiability in the complete case implies identifiability when MCAR or MAR values occur.

Associated EM algorithm

Assuming identifiability, we estimate parameters via likelihood maximization using the EM algorithm specifically designed for Gaussian, Poisson, multinomial and mixed data with MNAR data. Details of the algorithm are given in Appendix B. Assuming that the number K of clusters is known (its choice in practice is discussed in Section 4) and that the samples (y i , z i , c i ) i=1,...,n are i.i.d., the complete-data log-likelihood can be written as

comp (θ; Y, Z, C) = n i=1 K k=1 z ik log (π k f k (y i ; λ k )f k (c i | y i ; ψ k )) . (5)
The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] is an iterative algorithm that permits to maximize the likelihood function under missingness. Initialized at the point θ [0] , its iteration [r] consists, at the E-step, in computing the expectation of the complete-data log-likelihood

Q(θ; θ [r-1] ) = E θ [r-1] comp (θ; Y, Z, C) | Y obs , C
, then, at the M-step, updating the parameters by maximizing this function θ [r] = arg max θ Q(θ; θ [r-1] ). Note that

Q(θ; θ [r-1] ) = n i=1 K k=1 t ik (θ [r-1] ) log(π k ) + τ y (λ k ; y obs i , c i , θ [r-1] ) + τ c (ψ k ; y obs i , c i , θ [r-1] ) , where t ik (θ [r-1] ) = 1 f (y obs i , c i ; θ [r-1] ) π [r-1] k f k (y i ; λ [r-1] k )f k (c i | y i ; ψ [r-1] k )dy miss i , τ y (λ k ; y obs i , c i , θ [r-1] ) = E θ [r-1] log f k (y i ; λ k ) | y obs i , c i , z ik = 1 , τ c (ψ k ; y obs i , c i , θ [r-1] ) = E θ [r-1] log f k (c i | y i ; ψ k ) | y obs i , c i , z ik = 1 .
Thus, the iteration [r] of the EM algorithm is defined by

• E-step: Computation of t ik (θ [r-1]
), τ y (λ k ; y obs i , c i , θ [r-1] ) and τ c (ψ k ; y obs i , c i , θ [r-1] ). • M-step: Updating the parameters 1] ).

λ [r] k = arg max λ k n i=1 t ik (θ [r-1] )τ y (λ k ; y obs i , c i , θ [r-1] ) π [r] k = 1 n n i=1 t ik (θ [r-1] ) ψ [r] k = arg max ψ k n i=1 t ik (θ [r-1] )τ c (ψ k ; y obs i , c i , θ [r-
The E-step requires to be able to integrate the distribution of y mis i given (y obs i , z ik = 1, c i ) and the M-step requires to maximize the resulting function. This is straightforward under the MNARz model, because the effect of the missingness does not depend on y i (see Appendix B.1 for computation details in the case of Gaussian or categorical data).

Beyond the MNARz model

In Section 2.2, we proposed the MNARz model in [START_REF] Bouveyron | High-dimensional data clustering[END_REF]. A more general model, called MNARy k z j , can be considered, when the effect of missingness is on both the class membership and the variable itself:

f k (c i | y i ; ψ k ) = d j=1 (ρ(α kj + β kj y ij )) cij (1 -ρ(α kj + β kj y ij )) 1-cij , (6)
where ψ k = (α k1 , β k1 , . . . , α kK , β kK ). The parameter α kj represents a mean effect of missingness on the k-th class membership for the variable j (note that within a same class k, α kj is not necessarily equal to α kj for j = j ). The parameter β kj represents the direct effect of missingness on the variable j which depends on the class k as well.

Simpler models can also be derived from (6) by imposing equal parameters either across the class membership, or across the variables likely to be missing. First, we introduce three models, with a lower complexity than [START_REF] Ramoni | Bayesian clustering by dynamics[END_REF], that still allow the probability of being missing to depend on both the variable itself and the class membership:

• MNARyz j model: when β 1j = . . . = β Kj , ∀j, the effect of missingness on a variable is the same regardless of the class (while keeping different mean effects α kj on the class membership). • MNARy k z model: when α k1 = . . . = α kd , ∀k, the missingness has a same mean effect on class membership shared by all variables (while allowing different selfmasked and class-wise parameters β kj ). • MNARyz model: when β 1j = . . . = β Kj , ∀j and α k1 = . . . = α kd , ∀k, the effects on a particular variable and on the class membership can be respectively the same for all the classes and for all the variables.

Secondly, the probability to be missing can also depend only on the variable itself:

• MNARy model2 : when α 11 = . . . = α 1d = α 21 = . . . = α Kd and β 1j = . . . = β Kj ∀j, the only effect of missingness is thus on the variable j, being the same regardless of the class membership. • MNARy k model: when α 11 = . . . = α 1d = α 21 = . . . = α Kd , the effect of missingness on the variable j depends on the class k.

Thirdly, the probability to be missing can also depend only on the class membership, so that the missingness is class-wise only. This is the case of the MNARz model given in 3, but we can also consider a slightly more general case:

• MNARz j model: when β kj = 0, ∀(k, j), the effect of missingness on the class membership k is not the same for all the variables. Finally, the simplest model is the missing completely at random (MCAR) one, characterized by no dependence on variables, neither on class membership, i.e., each variable has the same probability of missing, MCAR:

β kj = 0, ∀(k, j) and α 1j = . . . = α Kj , ∀j. (7)
For each of these MNAR models, we have studied the identifiability and have proposed a specific algorithm (EM or Stochastic EM); all the details are given in the accompanying note available here. Note that for MNAR data, beyond the clustering task, the main challenge to overcome consists in proving the identifiability of the parameters of the data and the missing-data pattern distributions [START_REF] Molenberghs | Every missingness not at random model has a missingness at random counterpart with equal fit[END_REF]. The identifiability study showed that the most general models lead to non-identifiable parameters for categorical data (but the identifiability holds only for the MCAR, MNARz and MNARz j mechanisms).

Despite the possibility of defining a large number of MNAR models, we have chosen to focus on the MNARz mechanism, which is a good compromise, clearly outperforming methods that do not consider MNAR data, while limiting the computational cost of the estimation in regard of more general MNAR mechanisms. Moreover, the MNARz model is robust to model mispecification (see Section 4, Figure 3).

Numerical experiments on synthetic data

To assess the quality of the clustering, it is possible to use an information criterion such as the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF] or the Integrated Completedata Likelihood (ICL) [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. The BIC criterion is expected to select a relevant mixture model from a density estimation perspective, while the ICL is expected to select a relevant mixture model for a clustering purpose [START_REF] Baudry | Estimation and model selection for model-based clustering with the conditional classification likelihood[END_REF]. Thus, we consider the latter in the following. As the ICL involves an integral which is generally not explicit, we can use an approximate version [START_REF] Baudry | Estimation and model selection for model-based clustering with the conditional classification likelihood[END_REF] that we adapt with missing data. For a model M with ν M parameters, the ICL reads as where θM is a maximum likelihood estimator, (θ; Y obs , C) is the observed loglikelihood, and with

ICL(M) = ( θM ; Y obs , C) - ν M 2 log n + n i=1 K k=1 z MAP ik ( θM ) log(P(z ik = 1|y obs i , c i ; θM )),
z MAP ik (θ) = argmax k∈{1,...,K} P(z ik = 1|y obs i , c i ; θ). (8)
In addition, the Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] can be computed between the true partition Z and the estimated one.

Leveraging from MNAR data in clustering

MNAR data are often considered as a real obstacle for statistical processing. Yet, the following numerical experiment illustrates that the MNAR mechanism may help performing the clustering task. Let us consider a bivariate isotropic Gaussian mixture model with two components and equal mixing proportions, under the MNARz mechanism (3) with a probit link function. The difference between the centers of both mixture components is taken as ∆ µ = µ 21 -µ 11 = µ 22 -µ 12 ∈ {0.5, 1, . . . , 3}. This cluster overlap controls the mixture separation, which can vary from a low separation (∆ µ = 0.5) to a high separation (∆ µ = 3). We also make the discrepancy between inter-cluster missing proportions ∆ perc = |perc 2 -perc 1 |, vary in {0, 0.1, 0.2, 0.3}3 corresponds to emphasize the MNAR evidence: indeed, ∆ perc = 0 corresponds to a MCAR model, whereas a high value of ∆ perc corresponds to a high difference of missing pattern proportions between clusters. For all possible values of (∆ µ , ∆ perc ), 15% missing values are introduced. Figure 2 gives the theoretical ARI (i.e., we compute the ARI with the theoretical parameters) as a function of the cluster overlap ∆ µ and the MNAR evidence ∆ perc . Although the good classification rate is mainly influenced by center separation ∆ µ , it also increases with the MNAR evidence ∆ perc . When classification is difficult because the mixture is not well separated (∆ µ = 0.5), the fact that the data is MNAR helps clustering: the theoretical ARI for ∆ perc = 30 (MNAR data) is significantly higher than the one for ∆ perc = 0 (MCAR data). This toy example illustrates how clustering can leverage from MNAR values, rather generally considered a true hindrance for any statistical analysis.

Generic experiments

We consider a Gaussian mixture with three components having unequal proportions (π 1 = 0.5, π 2 = π 3 = 0.25) and independent variables:

∀j ∈ {1, . . . , d}, y ij = δ 3 k=1 ϕ kj z ik + ij , (9)
with ij ∼ N (0, 1) the noise term, ϕ k ∈ {0, 1} d and δ > 0. Thus, each entry y ij follows a Gaussian distribution with variance 1 and mean δ 3 k=1 ϕ kj z ik . The values of ϕ kj are arbitrary chosen and highlight the interations between the variable j and the class membership k. This formulation allows to control, in any scenario, the theoretical rate of misclassification through the value of δ (and hence the theoretical ARI). We introduce missing values with a MNAR model ( 6), using a probit link function and control the rates of missingness through the value of ψ k . For each experiment, the values of δ, ψ and φ are given in Appendix D. When not specified, the simulations have been performed for a theoretical rate of misclassification of 10% and a theoretical missing rate in the whole dataset of 30%.

Comparison of MNARz with other MNAR settings

We first vary the number of variables (d = 3, 6, 9) and consider n = 100 observations. The missing values are sequentially introduced with a MNAR setting. We compare the method considering the true mechanism (the one used to generate the missing values) with the EM algorithm for MCAR and MNARz values and the two-step heuristic based on Mice. This latter consists of first imputing the missing values using multiple imputations by chained equations [START_REF] Buuren | mice: Multivariate imputation by chained equations in r[END_REF] to get M completed datasets. Then, classical model-based clustering is performed on each completed dataset, for which the ARI is computed, Figure 3 shows the boxplot of the ARI for each scenario. First, the methods that consider a MNAR mechanism (MNAR ) always outperform those that consider the MCAR mechanism and the two-step procedure based on Mice. Finally, the MNARz model remains a good compromise, clearly outperforming methods that do not consider MNAR data.

Moreover, in Appendix C, Figure C1 and C2 show the computation times for these numerical experiments; while the MNAR models considering that the probability of being missing depends on the variable itself are computationally very costly, MNARz model clearly limits the computational cost of the estimation.

Focus on the MNARz mechanism

Considering the setting (9) and under a MNARz mechanism, we then evaluate the impact of misspecification of the link function (Figure 4 In Figure 4(a), our algorithm always considering a probit function gives the best ARI (outperforming strategies assuming only MCAR data) regardless of the link function (Laplace distribution, logit, probit) used to introduce missing values under a MNARz model. This highlights the robustness of the MNARz setting to the link function. In Figure 4(b), we consider a three-component Gaussian mixture with non-diagonal covariance matrices. For each component, the diagonal terms of the covariance matrix are Σ ii = 1 and the other terms Σ ij = , i = j, with ∈ {0, 0.1, 0.25, 0.5}, while the algorithms assume = 0. If the EM algorithm designed for MNARz data suffers from a huge deviation ( = 0.5) regarding the data distribution, it remains competitive for smaller ones ( = 0.1, 0.25). Finally, Figure 4(c) shows the boxplots of the ARI for 10%, 30% and 50% of missing values in the entire dataset. As the percentage of missing data increases, the gap between algorithms considering MCAR and MNARz data is widening, proving the relevancy of our algorithm even with high missing-data rates (50%). In Appendix C, we also provide the experiments for a theoretical rate of misclassification of 15%. Same conclusions hold.

When the number K of clusters is not known a priori, it can be automatically chosen using the ICL criterion: the idea is to run algorithms with several values for K (K = 1, 2, 3, 4 here), and to choose the model with the highest resulting ICL. To our knowledge, no method proposes an automatic choice of the number of clusters in unsupervised classification for the two-step heuristics, which is also a major drawback. We also illustrate in Appendix C (Figure C2) the findings of Theorem 1, by comparing our algorithm considering MCAR or MNARz data with the algorithm of the RMixtComp package [START_REF] Biernacki | Model-based clustering with mixed/missing data using the new software mixtcomp[END_REF] considering MCAR data and using the augmented data matrix (Y |C). As expected, both approaches give similar results.

TraumaBase ® dataset

In this section, we illustrate our approach on a public health application with the TraumaBase ® Group (https://www.traumabase.eu/en US) on the management of traumatized patients. This dataset contains 41 mixed variables (continuous, quantitative) on 8, 248 polytraumatized patients who suffer from a major trauma (injuries from cycle or car accident). Data have been collected from 15 different hospitals. In this dataset, 11% of the data are missing and only 1.4% of the individuals are fully observed. More information on the variables can be found in Appendix E. The purpose of this real data analysis is twofold: (i) we want to know if considering the missingness process has an impact on the estimated partition, (ii) we compare our method with the classical imputation methods in Appendix E. After discussion with doctors, some variables can be considered to have MNAR values, such as Shock.index.ph, which denotes the ratio between heart rate and systolic arterial pressure. In fact, if this rate has a value that indicates that the patient's condition is critical, doctors cannot measure heart rate or systolic arterial pressure in emergency situations. Therefore, we expect that considering a MNAR mechanism can improve the classification.

We compare our algorithm designed for the MNARz data (3) and the MCAR data [START_REF] Xiong | Time series clustering with arma mixtures[END_REF]. Figure 5a presents the ICL values in the Traumabase dataset for different numbers of classes. If both algorithms select K = 3 number of classes, the ICL of the algorithm which considers MNARz data is nonetheless always significantly higher than that of the algorithm for MCAR data. Their corresponding ARI between classifications obtained assuming either MNARz or MCAR mechanisms is about 0.90. Thus, both partitions are close but not equal, which may reflect the influence of the mechanism. To deepen this issue, we focus on the variable Shock.index.ph. Table 5b and 5c compare the performances of the algorithm considering MNARz data with the one considering MCAR data in terms of modelling of the marginal distribution of Shock.index.ph and partition estimation. As the values can be compared only up to label swapping, we notice that the minimum values (on the diagonals) are significantly higher than zero, which indicates that there is an influence of the MNARz mechanism on the modeling of the data and on the classification rules.

In Appendix E, we assess also the results by using the function catdes of the R package FactoMineR [START_REF] Lê | FactoMineR: An R package for multivariate analysis[END_REF] which allows to see how the cluster of the classification is described by the variables. The three groups described with our algorithm assuming MNARz data seem to be described by the same characteristics than those given by the doctors: the first group is formed by patients with a higher mortality rate and more severe injuries than the average population and the third group by a lower mortality rate and less serious injuries, whereas the second group may correspond to other cases. Note that this classification was done without using variables related to patient death and that it is quite striking to retrieve the same characteristics. This reinforces the idea that the classification obtained makes sense and may provide other information than the one of the doctors, taking into account more variables.

Concluding remarks

This paper addresses model-based unsupervised learning when MNAR values occur. We propose to cluster individuals via an estimation of the mixture model parameters in play. A by-product of such an approach is that the missing values can be also imputed, once the distribution is estimated. To this end, we have proposed an approach which embeds MNAR data directly within model-based clustering algorithms, in particular the EM algorithm. We have discussed several possible MNAR specifications. However, the numerical experiment leads us to recommend using algorithms considering a simple missing-data mechanism, the MNARz mechanism, which models the probability of being missing only depending on the class membership. By its very simplicity, the latter is indeed able to straightforwardly deal with any kind of data. In addition to being interpretable (which is especially important for real applications), this MNARz mechanism can be apprehended as a MAR one on the augmented matrix [Y |C], including the missing-data pattern C (Theorem 1). This echoes a widely-used approach in practice, not theoretically studied so far.

The seminal motivation of this work was clustering patients of the Traumabase dataset, in particular to assist doctors in their medical care. After a first conclusive application, there are still key challenges to make this work entirely applicable to real datasets. First, if our methodology can be applied to mixed data (categorical/quantitative), a straightforward extension of the proposed approach should be doable to handle variables that are not necessarily of the same type (MCAR, MAR and MNAR variables are indeed often coupled). Without any prior help from experts, this actually remains an open question to automatically evaluate the missing type of variables. Note however that one can arbitrate between the presented MNAR mechanisms using the ICL criterion, at the price of running multiple times the algorithm for the different MNAR scenarios.

Appendix A Proof of Theorem 1

Proof of Theorem 1. We denote by (c 1 , . . . , cn ) the patterns of missing data associated to the observed data ỹobs . It is thus the concatenation ci = (c i , 0 d ) of c i with the zero vector 0 d = (0, . . . , 0) of length d. Since all c i values are observed in ỹobs i , it is the reason why the last d values in ci are fixed to zero. Then, the MAR assumption indicates that P(c i | ỹi , z i ; ζ) = P(c i | ỹobs i ; ζ), with ζ the related parameter. Consequently, using the MAR assumption and the i.i.d. assumption of all uplets (ỹ i , z i , ci ), the whole likelihood can be decomposed into two likelihoods, one has

L(θ, ζ; Ỹ obs , C) = n i=1 f (ỹ i , ci ; θ, ζ)dỹ mis i = n i=1 f (ỹ i ; π, λ, ψ)f (c i |ỹ i ; ζ)dỹ mis i = n i=1 f (c i |ỹ obs i ; ζ) × Y mis i f (ỹ i ; π, λ, ψ)dỹ mis i = n i=1 L(ζ; ci | ỹobs i ) × n i=1 L(π, λ, ψ; ỹobs i ). L(π, λ, ψ, ζ; ỹobs i , ci ) = L(ζ; ci | ỹobs i ) × L(π, λ, ψ; ỹobs i ). (A1)
Providing that (π, λ, ψ) and ζ are functionally independent (ignorability of the MAR mechanism), the maximum likelihood estimate of θ = (π, λ, ψ) is obtained by maximizing only L(π, λ, ψ; ỹobs i ), and does not depend on L(ζ; ci | ỹobs i ). Finally, by using (4), the observed likelihood L(π, λ, ψ; ỹobs i ) is

L(π, λ, ψ; ỹobs i ) = K k=1 π k f k (y obs i ; λ k ) d j=1 ρ(α k ) cij ρ(α k ) (1-cij ) (A2) = K k=1 π k f k (y obs i ; λ k ) d j=1 f (c ij | z ik = 1; ψ). ( A3
)
As P(c ij | z ik = 1; ψ) corresponds to the MNARz definition (3), the observed likelihood L(π, λ, ψ; ỹobs i ) is equal to the full observed likelihood L(π, λ, ψ; y obs i , c i ) associated to the MNARz model,

L(π, λ, ψ; y obs i , c i ) = K k=1 π k f k (y obs i ; λ k ) d j=1 f (c ij | z ik = 1; ψ).

It leads to the decomposition

Q(θ; θ [r-1] ) = n i=1 K k=1 t ik (θ [r-1] ) log(π k ) + τ y (λ k ; y obs i , c i , θ [r-1] ) + τ c (ψ k ; y obs i , c i , θ [r-1] ) ,
where the terms involved in this decomposition are now detailed.

(a) the expectation of the data mixture part over the missing values given the available information (i.e., the observed data and the indicator pattern), the class membership and the current value of the parameters:

τ y (λ k ; y obs i , c i , θ [r-1] ) = E θ [r-1] log f k (y i ; λ k ) | y obs i , z ik = 1, c i , (b)
the expectation of the missing mechanism part over the missing values given the available information, the class membership and the current value of the parameters:

τ c (ψ k ; y obs i , c i , θ [r-1] ) = E θ [r-1] log f k (c i | y i ; ψ k ) | y obs i , z ik = 1, c i .
(c) the conditional probability for an observation i to belong to the class k given the available information and the current value of the parameters:

t ik (θ [r-1] ) = f (z ik = 1 | y obs i , c i ; θ [r-1] ).
Terms (a) and (b) require to integrate over the distribution f (y mis i | y obs i , z ik = 1, c i ; θ [r-1] ). For Term (a), one has

f (y mis i | y obs i , z ik = 1, c i ; θ [r-1] ) = f (y mis i , y obs i , z ik = 1, c i ; θ [r-1] ) f (y obs i , z ik = 1, c i ; θ [r-1] ) (B4) = f (c i | y mis i , y obs i , z ik = 1; ψ [r-1] )f (y mis i , y obs i , z ik = 1; λ [r-1] ) Y mis i f (c i | y mis i , y obs i , z ik = 1; ψ [r-1] )f (y mis i , y obs i , z ik = 1; λ [r-1] )dy mis i . (B5)
Term (c) corresponds to the conditional probability for an observation i to arise from the kth mixture component with the current values of the model parameter. More particularly, one has

t ik (θ [r-1] ) = f (z ik = 1, y obs i , c i ; θ [r-1] ) f (y obs i , c i ; θ [r-1] ) = f (z ik = 1, y obs i , c i ; θ [r-1] ) K h=1 f (z ih = 1, y obs i , c i ; θ [r-1] ) = π [r-1] k f (y obs i | z ik = 1; λ [r-1] k )f (c i | y obs i , z ik = 1; θ [r-1] ) K h=1 π [r-1] h )f (y obs i | z ih = 1; λ [r-1] h )f (c i | y obs i , z ih = 1; θ [r-1] ) = π [r-1] k f k (y obs i ; λ [r-1] k )f (c i | y obs i , z ik = 1; θ [r-1] ) K h=1 π [r-1] h f h (y obs i ; λ [r-1] h )f (c i | y obs i , z ih = 1; θ [r-1] ) (B6)

B.1 Gaussian mixture for continuous data

The pdf f k (y i ; λ) = φ(y i ; µ k , Σ k ) is assumed to be a Gaussian distribution with mean vector µ k and covariance matrix Σ k . First, let us detail the terms of the E-step. Term (a) is written as follows:

E log(φ(y i ; µ k , Σ k )) | y obs i , z ik = 1, c i ; θ [r-1] = - 1 2 [n log(2π) + log((| Σ k |))] - 1 2 E (y i -µ k ) T (Σ k ) -1 (y i -µ k ) | y obs i , z ik = 1, c i ; θ [r-1] .
This last term could be expressed using the commutativity and linearity of the trace function:

E (y i -µ k ) T (Σ k ) -1 (y i -µ k ) | y obs i , z ik = 1, c i ; θ [r-1] = tr(E (y i -µ k )(y i -µ k ) T | y obs i , z ik = 1, c i ; θ [r-1] (Σ k ) -1 ).
Finally note that only 1] has to be calculated.

E (y i -µ k )(y i -µ k ) T | y obs i , z ik = 1, c i ; θ [r-
For the MNARz model, the effect of the missingness is only due to the class membership.

• For Term (a), note that

f (y mis i | y obs i , z ik = 1, c i ; θ [r-1] ) = f (y mis i | y obs i , z ik = 1; λ [r-1] ),
which makes the computation easy. Indeed, using (B5),

f (y mis i | y obs i , z ik = 1, c i ; θ [r-1] ) = d j=1 ρ(α [r-1] k ) cij (1 -ρ(α [r-1] k )) 1-cij f (y mis i , y obs i , z ik = 1; λ [r-1] ) Y mis i d j=1 ρ(α [r-1] k ) cij (1 -ρ(α [r-1] k )) 1-cij f (y mis i , y obs i , z ik = 1; λ [r-1] )dy mis i = f (y mis i , y obs i , z ik = 1; λ [r-1] ) Y mis i f (y mis i , y obs i , z ik = 1; λ [r-1] )dy mis i = f (y mis i | y obs i , z ik = 1; λ [r-1] ), since d j=1 ρ(α [r-1] k ) cij (1 -ρ(α [r-1] k
)) 1-cij does not depend on y mis i and is simplified with the numerator. The law of (y mis i | y obs i , z ik = 1) is Gaussian (up to a reorganization of the variables associated to individual i). Noting that

y i | z ik = 1; λ [r-1] = y obs i y mis i | z ik = 1; λ [r-1] ∼ N (µ obs ik ) [r-1] (µ mis ik ) [r-1] ,
(Σ obs,obs ik

) [r-1] (Σ obs,mis ik ) [r-1] (Σ mis,obs ik ) [r-1] (Σ mis,mis ik ) [r-1] ,
one obtains

y mis i | y obs i , z ik = 1; λ [r-1] ∼ N (μ mis ik ) [r-1] , ( Σmis ik ) [r-1] . (B7)
with (μ mis ik ) [r-1] and ( Σmis ik ) [r-1] the standard expression of the mean vector and covariance matrix of a conditional Gaussian distribution (see for instance [START_REF] Anderson | An Introduction to Multivariate Statistical sAnalysis[END_REF]) detailed as follows

(μ mis ik ) [r-1] = (µ mis ik ) [r-1] + (Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1] -1 y obs i -(µ obs ik ) [r-1] , (B8) ( Σmis ik ) [r-1] = (Σ mis,mis ik ) [r-1] -(Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1] -1 (Σ obs,mis ik ) [r-1] . (B9)
Note also that we have

(y i -µ k )(y i -µ k ) T = (y obs i -µ obs ik ) T (y obs i -µ obs ik ) (y obs i -µ obs ik ) T (y mis i -µ mis ik ) (y mis i -µ mis ik ) T (y obs i -µ obs ik ) (y mis i -µ mis ik ) T (y mis i -µ mis ik )
.

Therefore, the expected value of each block for the current parameter value is

E (y obs i -µ obs ik ) T (y obs i -µ obs ik ) | y obs i , z ik = 1; λ [r-1] = (y obs i -µ obs ik ) T (y obs i -µ obs ik ) E (y obs i -µ obs ik ) T (y mis i -µ mis ik ) | y obs i , z ik = 1; λ [r-1] = (y obs i -µ obs ik ) T ((μ mis ik ) [r-1] -µ mis ik ) E (y mis i -µ mis ik ) T (y mis i -µ mis ik ) | y obs i , z ik = 1; λ [r-1] = ((μ mis ik ) [r-1] -µ mis ik ) T ((μ mis ik ) [r-1] -µ mis ik ) + ( Σmis ik • For Term (b), f (c i | y i , z ik = 1; ψ) is independent of y i , which implies log(f (c i | z ik = 1; ψ)) = d j=1 c ij log ρ(α k ) + (1 -c ij ) log(1 -ρ(α k )) (B10)
• For Term (c), one first remark that

P(c i | y obs i , z ik = 1; θ [r-1] ) = d j=1 P(c ij = 1 | y obs i , z ik = 1; θ [r-1] ) cij P(c ij = 0 | y obs i , z ik = 1; θ [r-1] ) 1-cij .
In particular, for MNARz, by independence of y i , one has

P(c ij = 1 | y obs i , z ik = 1; θ [r-1] ) = P(c ij = 1 | z ik = 1; θ [r-1] ) = ρ(α k )
Using (B6), one obtains

t [r-1] ik (θ [r-1] ) = π [r-1] k φ(y obs i ; (µ obs ik ) [r-1] , (Σ obs,obs ik ) [r-1] ) d j=1 ρ(α [r-1] k ) cij (1 -ρ(α [r-1] k )) 1-cij K h=1 π [r-1] h φ(y obs i ; (µ obs ih ) [r-1] , (Σ obs,obs ih ) [r-1] ) d j=1 ρ(α [r-1] h ) cij (1 -ρ(α [r-1] h )) 1-cij (B11)
If ρ is the logistic distribution, the expression can be written more simply

t ik (θ [r-1] ) ∝ π [r-1] k φ(y obs i ; λ [r-1] k ) d j=1 1 + exp(-δ ij α [r-1] k ) -1
where

δ ij = 1 if c ij = 1 -1 otherwise.
Finally, the E-step and the M-step can be sketched as follows in the Gaussian mixture case.

E-step

The E-step for Term (a) consists of computing for k = 1, . . . , K and i = 1, . . . , n

(μ mis ik ) [r-1] = (µ mis ik ) [r-1] + (Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1] -1 y obs i -(µ obs ik ) [r-1] ( Σmis ik ) [r-1] = (Σ mis,mis ik ) [r-1] -(Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1] -1 (Σ obs,mis ik ) [r-1] (ỹ i,k ) [r-1] = (y obs i , (μ mis ik ) [r-1] ) Σ[r-1] ik = 0 obs,obs i 0 obs,mis i 0 mis,obs i ( Σmis ik ) [r-1] .
Note that whenever the mixture covariance matrices are supposed diagonal then ( Σmis ik ) [r-1] is also a diagonal matrix. Term (c) also requires the computation of t ik (θ [r-1] ) given in (B11) for k = 1, . . . , K and i = 1, . . . , n.

M-step

The maximization of Q(θ; θ [r-1] ) over (π, λ) leads to, for k = 1, . . . , K,

π [r] k = 1 n n i=1 t ik (θ [r-1] ) µ [r] k = n i=1 t ik (θ [r-1] )(ỹ k,i ) [r-1] n i=1 t ik (θ [r-1] ) Σ [r] k = n i=1 t ik (θ [r-1] ) (ỹ i,k ) [r-1] -µ r k )((ỹ i,k ) [r-1] -µ r k ) T + Σ[r-1] ik n i=1 t ik (θ [r-1] ) .
Then, the maximization of Q(θ; θ [r-1] ) over ψ can be performed using a Newton Raphson algorithm. For k = 1, . . . , K, it remains to fit a generalized linear model with the binomial link function for the matrix (J MNARz k ) [r] and by giving t ik (θ [r-1] ) as prior weights to fit the process.

(J MNARz k ) [r] =    c .1 1 . . . . . . c .d 1    . (B12)
The EM algorithm for the MNARz model is described in Algorithm 1 for Gaussian mixture.

B.2 Latent class model for categorical data

For categorical data, we have φ(

y i ; λ k ) = d j=1 φ(y ij ; λ kj ) = d j=1 j =1 (λ kj ) y ij . Term (a) is E log(φ(y i ; p k )) | y obs i , z ik = 1, c i ; λ [r-1] = j,cij =0 j =1 y ij + j,cij =1 j =1 log(λ y ij kj ) (B13) Term (b)
is the same as in the Gaussian case given in (B10). Finally, the EM algorithm can be summarized as follows E step: For k = 1, . . . , K and i = 1, . . . , n, compute

t ik (θ [r-1] ) = π [r-1] k j,cij =0 j =1 (λ kj ) y ij ρ(α k ) K h=1 π [r-1] h j,cij =0 j =1 (λ hj ) y ij ρ(α h ) (ỹ ij,k ) [r-1] = c ij (θ kj ) [r-1] + (1 -c ij )y ij , ∀j = 1, . . . , d, ∀ = 1, . . . , j .
for i = 1 to n, k = 1 to K do 6: (μ mis ik ) [r-1] = (µ mis ik ) [r-1] +(Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1]
-1

y obs i -(µ obs ik ) [r-1] . 7: ( Σmis ik ) [r-1] = (Σ mis,mis ik ) [r-1] -(Σ mis,obs ik ) [r-1] (Σ obs,obs ik ) [r-1] -1 (Σ obs,mis ik ) [r-1] .
8:

(ỹ i,k ) [r-1] = (y obs i , (μ mis ik ) [r-1]
).

9: 1] , where 0 obs,obs

Σ[r-1] ik = 0 obs,obs i 0 obs,mis i 0 obs,mis i ( Σmis ik ) [r-
t ik (θ [r-1] ) ∝ π [r-1] k φ(y obs i ; (µ obs ik ) [r-1] , (Σ obs,obs ik ) [r-1] ) d j=1 ρ(α [r-1] k ) cij (1 - ρ(α [r-1] k )) 1-cij 11:
end for 12:

M-step:

13:

for k = 1 to K do 14: π [r] k = 1 n n i=1 t ik (θ [r-1] ), µ [r] k = n i=1 t ik (θ [r-1] )(ỹ k,i ) [r-1] n i=1 t ik (θ [r-1] ) 15: Σ [r] k = n i=1 t ik (θ [r-1] ) (ỹ i,k ) [r-1] -µ [r] k )((ỹ i,k ) [r-1] -µ [r] k ) T + Σ[r-1] ik n i=1 t ik (θ [r-1] ) 16: Let ψ [r]
k be the coefficients of a GLM with a binomial link function, by giving prior weights t ik (θ [r-1] ). In particular, the optimization problem is

(J MNARz k ) [r] ψ [r] k = log 1 -E[c|(J MNARz k ) [r] ]) E[c|(J MNARz k ) [r] ] ,
for a matrix (J MNARz k ) [r] given in (B12) and c = (c .1 , . . . , c .d ) the concatenated missing data patterns for the variables 1, . . . , d.

17:

end for 18: end for M step: The maximization of Q(θ; θ [r-1] ) over θ leads to, for k = 1, . . . , K, The M-step for ψ consists of performing a GLM with a binomial link for the following matrix: Fig. E6: Plot of the catdes output for the EM algorithm applied to the Traumabase dataset considering MNARz data. We consider K = 3 clusters. A cell colored in blue (resp. red) means that the average value of the variable in the given cluster is significantly lower (resp. higher) than in the overall data. In addition, the more transparent the cell is, the less significant the difference between the behavior of the variable in the given cluster and in the overall data is.

π r k = 1 n n i=1 t ik (θ [r-1] ) (θ kj ) r = n i=1 t ik (θ [r-1] )(ỹ ij,k ) [r-1] n i=1 t ik (θ [r-1] ) , ∀j = 1, . . . , d, ∀ = 1, . . . , j . d ϕ 3 ϕ 11 = ϕ 22 = ϕ 33 = 1 6 ϕ 11 = ϕ 22 = ϕ 33 = ϕ 14 = ϕ 36 = 1 9 ϕ 11 = ϕ 22 = ϕ 33 = ϕ 14 = ϕ 36 = ϕ 17 = ϕ 27 = ϕ 39 = 1
(H MNARz ) [r] =    c .1 z .1 . . . z .K . . . . . . . . . . . . c .d z .1 . . . z .K    =               c 11 z [r] 11 . . . z [r] 1K . . . . . . . . . . . . c n1 z [r] n1 . . . z [r] nK . . . . . . . . . . . . c 1d z [r] 11 . . . z [r] 1K . . . . . . . . . . . . c nd z [r] n1 . . . z [r] nK .               (

E.2 Imputation performances in the Traumabase dataset

We perform now simulations on the real dataset in order to be able to measure the quality of the imputation of our method compared to the multiple imputation [START_REF] Buuren | mice: Multivariate imputation by chained equations in r[END_REF] (Mice). We introduce some additional missing values in three quantitative variables (TCD.PI.max, Shock.index.ph, FiO2 ) by using the MNARz mechanism (3). The variables contain initially 51%, 31%, 7% and finally 63%, 50% and 32% missing values. The algorithm for continuous data specifically designed for MNARz data for K = 3 classes is compared with mean imputation and multiple imputation in terms of mean squared error (MSE). Denoting by Ŷ ∈ R n×d the imputed dataset and C ∈ R n×d the indicator pattern of missing data newly introduced, the mean squared error is given by

E[( Ŷ -Y ) C] 2 F E[Y C] 2 F ,
where is the Hadamard product and E[] 2 F = E[ . 2 F ] denotes the expectation of the Frobenius norm squared. In particular, to impute missing values using our clustering algorithm, we use the conditional expectation of the missing values given the observed ones, given that the data are assumed to be Gaussian and that all the parameters of the distribution are given by our algorithm. Imputation is carried out by taking the mean over 10 4 draws. In Figure E7, our clustering algorithm, designed for the MNAR setting, gives a significantly smaller error than other methods.

E.3 Description of the variables in the Traumabase dataset

A description of the variables which are used in Section 5 is given. Figure E8 gives the percentage of missing values per variable. The indications given in parentheses ph (pre-hospital) and h (hospital) mean that the measures have been taken before the arrival at the hospital and at the hospital. surgical intervention to reduce intracranial hypertension (h). • Death (categorical, binary variable, no missing values): death of the patient (h).

• AIS.head, AIS.face (ordinal, discrete, integers between 0 and 6 and 4 1.7% NA & 1.7% NA): Abbreviated Injury Score, describing and quantifying facial and head injuries (AIS = 0: no injury; the higher the more critical) (h). • AIS.external (continuous, discrete, integers between 0 and 5, 1.7% NA): Abbreviated Injury Score for ex-ternal injuries, here it is assumed to be a proxy of information avail-able/visible during pre-hospital phase (ph/h). • ISS (continuous, discrete, integers between 0 and 75, 1.6% NA): Injury Severity Score, sum of squares of top three AIS scores (h). • Activation.HS.procedure (categorical, binary variable, 3.7% NA): activation of hemorragic shock procedure in case of HS suspicio (h).

Fig. 1 :

 1 Fig. 1: Illustration of dependency between c i and y i in a MNARz model by drawing P(c i | y i ; π, λ, ψ) for a three component univariate Gaussian model with mixing proportions π 1 = π 2 = 0.3 and π 3 = 0.4, with centers µ 1 = µ 3 = -5 and µ 2 = 0, and with variances σ 2 k = k (k ∈ {1, 2, 3}). The MNARz parameters are fixed to α 1 = 2, α 2 = 0 and α 3 = 1.

Fig. 2 :

 2 Fig. 2: Relative effect of both the separation strength ∆ µ of the mixture component and the MNAR evidence ∆ perc on theoretical ARI (e.g. if ∆ perc = 10% (green line), the second class has 10% more missing values).

Fig. 3 :

 3 Fig.3: Boxplot of the ARI obtained for 50 samples composed of d = 3, 6, 9 variables (rows) and n = 100 observations. Missing values are introduced with MNARy, MNARy k z j MNARyz or MNARz j settings (columns). The boxplot in green is the one for the algorithm considering the true MNAR setting; the boxplot in blue (resp. in gray) is the one for the EM algorithm considering the MNARz setting (resp. the MCAR setting); the boxplot in red in the two-step heuristic (Mice). The red dashed line indicates the theoretical ARI.

Fig. 4 :

 4 Boxplot of the ARI obtained for 50 samples of dimension d = 6 variables. The missing values are introduced using a MNARz setting. The red dashed line indicates the theoretical ARI.

  ICL values for different K in the mixed case.

Fig. 5 :

 5 Fig. 5: Comparison of the MCAR and MNARz mechanism on Traumabase dataset.

20 Algorithm 1 2 :

 2012 EM algorithm for Gaussian mixture and MNARz model 1: Input: Y ∈ R n×d (matrix containing missing values), K ≥ 1, r max . Initialize π 0 k , µ 0 k , Σ 0 k and ψ 0 k , for k ∈ {1, . . . , K}. 3: for r = 0 to r max do

  of observed (reps. missing) variables for individual i.

  Choice of the values of δ and α for all the experiments of Section 4.2 for the MNARz j mechanism.

  1 -1.2 0.4 -0.1 -1.3 0.3 -0.1 -1

  Fig. C1: Boxplot of the computational times (in seconds) obtained for 50 samples composed of d = 6 variables.

  Fig. E7: Mean squared error of the imputation task for the Traumabase dataset.

Table 1 :

 1 Proportion of good selections of K (K = 3) using the ICL criterion for the EM algorithm, over 50 repetitions (d = 6).

		MCAR	MNARz
	Sample size n	100	500	100	500
	10 % NA	94% 100%	94% 100%
	30 % NA	8%	96%	56% 100%
	50 % NA	0%	0%	20% 98%

Table 1

 1 gathers the percentages of times (over 50 repetitions) the correct number of classes (K = 3) is chosen by the ICL criterion for different missing-data rates (10, 30, 50%) and different sample sizes (n = 100, 500). In any case, the EM algorithm for MNARz data always outperforms the algorithm for MCAR data in terms of accurate model selection. The EM algorithm for MNARz data manages also to select the best model despite a high percentage of missing data (50%) provided that the sample size is large enough (n = 500).

  Shock.index.ph considering MNARz data and the one considering MCAR data.

		MCAR	MNARz Class 1 Class 2 Class 3
		Class 1	0.03	0.47	0.61
		Class 2	0.45	0.05	0.25
		Class 3	0.63	0.23	0.03
	(b)	Total variation distance between
	the marginal distribution of the variable
		MCAR	MNARz Class 1 Class 2 Class 3
		Class 1	2.43	26.5	37.6
		Class 2	26.2	3.40	20.1
		Class 3	39.3	19.2	2.05
	(c)	Euclidean distance between condi-
	tional probabilities of the cluster member-
	ships given observed values of the vari-
	able Shock.index.ph, considering MNARz or
	MCAR data (more details in Appendix E).

Table D1 :

 D1 Choice of the values of ϕ and α for all the experiments of Section 4.2. Other values ϕ kj are null.

  If the data are mixed (continuous and categorical), the formulas can be extended straightforwardly if the continuous and the categorical variables are assumed to be independent knowing the latent clusters.

	K	% NA	link	rate of misclassification	l	δ	α
	3	30%	probit	90%	0	2.6	-1 -0.3 0
	3	30%	logit	90%	0	2.76	-1.5 -0.8 0.1
	3	30%	Laplace	90%	0	2.85	-1.1 0.3 0
	3	30%	probit	85%	0	2.27	-1 -0.3 0
	3	30%	logit	85%	0	2.44	-1.5 -0.8 0.1
	3	30%	Laplace	85%	0	2.46	-1.1 0.3 0
	3	30%	probit	90%	0.1	2.3	-1.16 0.3 -0.42
	3	30%	probit	90%	0.25 2.17	-1.16 0.3 -0.4
	3	30%	probit	90%	0.5	1.85	-1.16 0.3 -0.4
	3	30%	probit	85%	0.1	1.97	-1.16 0.3 -0.42
	3	30%	probit	85%	0.25 1.86	-1.16 0.3 -0.4
	3	30%	probit	85%	0.5	1.57	-1.16 0.3 -0.4
	3	10%	probit	90%	0	2.18	-1.65 -1.2 -0.9
	3	50%	probit	90%	0	3.3	-0.55 0.25 1.7
	3	10%	probit	85%	0	1.95	-1.65 -1.2 -0.9
	3	50%	probit	85%	0	2.62	-0.55 0.25 1.7	B14)
	B.3 Combining Gaussian mixture and latent class model for
		mixed data

Table D2 :

 D2 Choice of the values of δ and α for all the experiments of Section 4.2 and Appendix C for the MNARz mechanism. K denotes the number of class, the column denoted as % NA gives the rate of missingness, the column called link gives the link function of the missing-data mechanism used in the introduction of the missing values, l is the coefficient of correlation (antidiagonal terms), δ is given in[START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] and α in[START_REF] Ramoni | Bayesian clustering by dynamics[END_REF].

Table D4 :

 D4 Choice of the values of δ, α and β for all the experiments of Section 4.2 for the MNARy mechanism.

	d	δ	α	β
	3	4.72	-1.2 -0.8 -0.5	-3 0.3 1
	6	2.12	-1.35 -0.29 0	-3 0.3 -3 -3 -2 1
	9	1.71	-1.34 -0.34 0	-3 0.3 -3 -2.8 -2 1 0.2 0.1 0.4

Table D5 :

 D5 Choice of the values of δ, α and β for all the experiments of Section 4.2 for the MNARyz mechanism.

  Table5cgives the Euclidean distance between the conditional probabilities of the cluster memberships given the observed values of the variable Shock.index.ph obtained using the algorithm considering MNARz data and those obtained using the algorithm considering MCAR data. For clarity, the latter quantity is reported here, Boxplot of the ARI obtained for 50 samples composed of d = 6 variables.The missing values are introduced using a MNARz setting. The misclassification rate is of 15%. Impact of the misspecification of the data distribution. Boxplot of the ARI obtained for 50 samples composed of d = 6 variables. The missing values are introduced using a MNARz setting. The misclassification rate is of 15%. Impact of the percentage of missing values.with s the index of the variable Shock.index.ph, θ MCAR (resp. θ MNAR ) the estimator returned by the algorithm considering MCAR data (resp. MNAR data).

	l=0	l=0.1	l=0.25	l=0.5
	1.00
	MCAR MNARz Fig. C4: 10% NA 0.00 0.25 0.50 0.75 ARI MCAR MNARz MCAR MNARz 30% NA MCAR MNARz MCAR MNARz MCAR 50% NA MCAR MNARz MNARz 0.25 0.50 0.75 ARI 2 IG S .II D ea th D ea th _3 0d V as op re ss or .th er ap y T B I_ D ea th _3 0d T B I_ D ea th IS S Fig. C5: 1 F iO 2 A IS .h ea d T C D .P I.m ax A ct iv at io n. H S .p ro ce du re Tr an ex am ic .a ci d C ris ta llo id .v ol um e IIC P T B I C ol lo id .v ol um e Tr au m a. ce nt er S ho ck .in de x. ph D ec om pr es si ve .c ra ni ec to m y N eu ro su rg er y. da y0 A IS .fa ce E V D D el ta .h em oC ue m aj or E xt ra cr an ia l A IS .e xt er na l H R .p h H R .p h. m ax H em oC ue .in it S B P .p h D B P .p h S pO 2. ph .m in D B P .p h. m in S B P .p h. m in G C S .in it G C S .m ot or .in it G C S G C S .m ot or P up il. an om al y= M yd ria se (b ila té ra le ) C ar di ac .a rr es t.p h= O ui Im pr ov .a no m al y. os m o= N on O sm ot he ra py .p h= M an ni to l P up il. an om al y. ph =M yd ria se B ila té ra le P up il. an om al y. ph =P as pr éc is é O sm ot he ra py .p h= R ie n O sm ot he ra py =M an ni to l P up il. an om al y. ph =A ni so co rie (u ni la té ra le ) O sm ot he ra py =S S H P up il. an om al y= A ni so co rie (u ni la té ra le ) O sm ot he ra py .p h= N A P up il. an om al y. ph =N A Im pr ov .a no m al y. os m o= O ui O sm ot he ra py .p h= S S H A nt ip la te le t.t he ra py =N A A nt ic oa gu la nt .th er ap y= N A C ar di ac .a rr es t.p h= N A P up il. an om al y= N A A nt ip la te le t.t he ra py =O ui A nt ic oa gu la nt .th er ap y= N on A nt ip la te le t.t he ra py =N on O sm ot he ra py =R ie n C ar di ac .a rr es t.p h= N on Im pr ov .a no m al y. os m o= N on te st é P up il. an om al y= N on O sm ot he ra py .p h= P as de m yd ria se P up il. an om al y. ph =N on A nt ic oa gu la nt .th er ap y= O ui	Mechanism MCAR MNARz Mechanism MCAR MNARz 3
	n
	(P(z ik = 1|y obs is ; θ MCAR ) -P(z i k = 1|y obs is ; θ MNAR )) 2 , ∀k, k ∈ {1, 2, 3}
	i=1

  • Trauma.center (categorical, integers between 1 and 16, no missing values): name of the trauma center (ph & h). • Anticoagulant.therapy (categorical, binary variable, 4.3% NA): oral anticoagulant therapy before the accident (ph). • Antiplatelet.therapy (categorical, binary variable, 4.4% NA): anti-platelet therapy before the accident (ph). • GCS.init, GCS (ordinal, integers between 3 and 15, 2% NA & 42% NA): Initial Glasgow Coma Scale (GCS) on arrival on scene of enhanced care team and on arrival at the hospital (GCS = 3: deep coma; GCS = 15: conscious and alert) (ph & h). • GCS.motor.init, GCS.motor (ordinal, integers between 1 and 6, 7.6% NA & 43%): Initial Glasgow Coma Scale motor score (GCS.motor = 1: no response; GCS.motor = 6: obeys command/purposeful movement) (ph % h). • Pupil.anomaly.ph, Pupil.anomaly (categorical, 3 categories: Non, Anisocoire (unilaterale), Mydriase Bilaterale, 2% NA & 1.7%): pupil dilation indicating brain herniation (ph & h). • Osmotherapy.ph, Osmotherapy (categorical, 4 categories: Pas de mydriase, SSH, Mannitol, Rien, 1.7% NA and no missing values): administration of osmotherapy to alleviate compression of the brain (either Mannitol or hypertonic saline solution) (ph & h) • Improv.anomaly.osmo (categorical, 3 categories: Non testé, Non, Oui, no missing values): change of pupil anomaly after ad-ministration of osmotherapy (ph). • Cardiac.arrest.ph (categorical, binary variable, 2.3% NA): cardiac arrest during pre-hospital phase (ph). • SBP.ph, DBP.ph, HR.ph (continuous, 29.3% NA & 29.6% NA & 29.5% NA): systolic and diastolic arterial pressure and heart rate during pre-hospital phase (ph). • SBP.ph.min, DBP.ph.min (continuous, 12.8% NA & 13% NA): minimal systolic and diastolic arterial pressure during pre-hospital phase (ph). • HR.ph.max (continuous, 13.7 % NA): maximal heart rate during pre-hospital phase (ph). • Cristalloid.volume (continuous, positive values, 30% NA): total amount of prehospital adminis-tered cristalloid fluid resuscitation (volume expansion) (ph). • Colloid.volume (continuous, positive values, 31.3% NA): total amount of prehospital administered colloid fluid resuscitation (volume expansion) (ph). • HemoCue.init (continuous, 34.9% NA): prehospital capillary hemoglobin concentration (the lower, the more the patient is probably bleeding and in shock); hemoglobin is an oxygen carrier molecule in the blood (ph). • Delta.hemoCue (continuous, 37.2% NA): difference of hemoglobin level between arrival at the hospital and arrival on the scene (h). • Vasopressor.therapy (continuous, no missing values): treatment with catecholamines in case of physical or emotional stress increasing heart rate, blood pres-sure, breathing rate, muscle strength and mental alertness (ph). • SpO2.min (continuous, 11.7% NA): peripheral oxygen saturation, measured by pulse oxymetry, to estimate oxygen content in the blood (95 to 100%: considered normal; inferior to 90% critical and associated with considerable trauma, danger and mortality) (ph). • TCD.PI.max (continuous, 51.2% NA): pulsatility index (PI) measured by echodoppler sonographic examen of blood velocity in cerebral arteries (PI ¿ 1.2: indicates altered blood flow maybe due to traumatic brain injury) (h). • FiO2 (categorical, in {1, 2, 3, 4, 5}, 6.8% NA): inspired concentration of oxygen on ventilatory support (the higher the more critical; Ventilation = 0: no ventilatory support) (h). • Neurosurgery.day0 (categorical, binary variable, no missing values): neurosurgical intervention performed on day of admission (h). • IGS.II (continuous, positive values, 2% NA): Simplified Acute Physiology Score (h). • Tranexomic.acid (categorical, binary variable, no missing values): administration of the tranexomic acid (h). • TBI (categorical, binary variable, no missing values): indicates if the patient suffers from a traumatic brain injury (h). • IICP (categorical, binary variable, 70.9% NA): at least one episode of increased intracranial pressure; mainly in traumatic brain injury; usually associated with worse prognosis (h). • EVD (categorical, binary variable, no missing values): external ventricular drainage (EVD); mean to drain cerebrospinal fluid to reduce intracranial pressure (h). • Decompressive.craniectomie (categorical, binary variable, no missing values):

  • TBI Death (categorical, binary variable, no missing values): death of the patients suffering from a traumatic brain injury (h). • TBI Death 30d (categorical, binary variable, no missing values): death of the patients suffering from a traumatic brain injury in the 30 days (h). • TBI 30d (categorical, binary variable, no missing values): traumatic brain injury in the 30 days (h). • Death 30d (categorical, binary variable, no missing values): death in the 30 days (h). • Shock.index.ph (continuous, positive values, 30.5% NA): ratio of heart rate and systolic arterial pressure during pre-hospital phase (ph). • majorExtracranial (categorical, binary variable, no missing values): major extracranial lesion (h).

The code is available on https://github.com/AudeSportisse/Clustering-MNAR/tree/main.

This is actually a particular case of MNAR mechanims, widely used in practice[START_REF] Mohan | On handling self-masking and other hard missing data problems[END_REF].

The value ∆perc means that if the percentage of missing values in the first cluster is perc 1 , the percentage of missing values in the second cluster is perc 2 = (perc 1 + ∆perc).

Note that in FigureC2, the differences for n = 100 can be explained by the difference in initialization of the algorithms, which can play an important role for small sample sizes.

Appendix B Detail on EM algorithm

The EM algorithm consists on two steps iteratively proceeded: the E-step and M-step.

For the E-step, one has

).

Appendix D Complements on generic experiments

This section gives the values of δ (see ( 9)) ψ (see [START_REF] Ramoni | Bayesian clustering by dynamics[END_REF]) and ϕ (see ( 9)) used during the different experiments. As explained in Section 4.2, their choice allows to control the rates of misclassification and missingness, as well as the interation between the variables and the class membership. To estimate these values, we have generated a large sample (n = 10 5 ) and compute the misclassification rate and the missingness rate for several values of δ and ψ and pick the ones which correspond to the setting of the experiment.