Clustering data with non-ignorable missingness using semi-parametric mixture models assuming independence within components

Abstract

We propose a semi-parametric clustering model assuming conditional independence given the component. One advantage is that this model can handle non-ignorable missingness. The model defines each component as a product of univariate probability distributions but makes no assumption on the form of each univariate density. Note that the mixture model is used for clustering but not for estimating the density of the full variables (observed and unobserved). Estimation is performed by maximizing an extension of the smoothed likelihood allowing missingness. This optimization is achieved by a Majorization-Minorization algorithm. We illustrate the relevance of our approach by numerical experiments conducted on simulated data. Under mild assumptions, we show the identifiability of the model defining the distribution of the observed data and the monotonicity of the algorithm. We also propose an extension of this new method to the case of mixed-type data that we illustrate on a real data set. The proposed method is implemented in the R package MNARclust available on CRAN.

Appendices

Appendix A: Weak and strong ignorability for clustering

In a likelihood-based estimation, the missingness mechanism is said to be ignorable for likelihood inference if the missing data are missing at random and if the distinctness property is satisfied by the parameters (see Definition 6.4 in Little and Rubin 2002). These conditions ensure that it is appropriate to ignore the missingness mechanism, especially for parameter estimation. Indeed, when ignorability holds true, the parameters of the distribution of $\varvec{X}_i$ can be consistently estimated from the observed sample $(\varvec{x}^{\text {obs}}_1,\ldots ,\varvec{x}^{\text {obs}}_n)$ without modelling the missingness process. This framework has been extended to the case of ignorability for a subset of the parameters (Little et al. 2017). In such a case, despite the fact that the missingness mechanism is MNAR, the condition of ignorability for a subset of the parameters ensures that a subset of the parameters can be consistently estimated by ignoring the missingness mechanism. Thus, this section introduces the notion of ignorability for clustering that defines the conditions for the clustering to be consistently achieved without modeling the missingness process.

In clustering, there are two quantities of interest: the partition and the posterior probabilities of classification (see (4) and (5)). Thus, we introduce the notions of weakly and strongly ignorable mechanisms for clustering that allow the missingness mechanism to be neglected for estimating the partition and the posterior probabilities of classification respectively. In such a case, the marginal pdf of the observed variable $f(\varvec{x}^{\text {obs}}_i)=\sum _{k=1}^{K} \pi _k f_k(\varvec{x}^{\text {obs}}_i)$ where $f_k(\varvec{x}^{\text {obs}}_i)=\int \lambda _k(\varvec{x}_i) d\varvec{x}^{\text {miss}}_i$ can be used to compute the posterior probabilities of classification and the classification rule.

Definition 1

The missingness mechanism is said to be strongly ignorable for clustering if it can be neglected for computing posterior probabilities of classification as follows

$$\begin{aligned} \forall (\varvec{x}^{\text {obs}}_i,\varvec{r}_i),\; \frac{\pi _k g_{k}^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i) }{\sum _{\ell =1}^K \pi _\ell g_{\ell }^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i) }=\frac{\pi _k f_{k}(\varvec{x}^{\text {obs}}_i) }{\sum _{\ell =1}^K \pi _\ell f_{\ell }(\varvec{x}^{\text {obs}}_i) }. \end{aligned}$$

The missingness mechanism is said to be weakly ignorable for clustering if it can be neglected for computing the classification rule, which means that

$$\begin{aligned} \forall (\varvec{x}^{\text {obs}}_i,\varvec{r}_i),\; \zeta (\varvec{x}^{\text {obs}}_i) = \eta (\varvec{x}^{\text {obs}}_i), \end{aligned}$$

where $\eta (\varvec{x}^{\text {obs}}_i)={{\,\mathrm{arg\,max}\,}}_{k=1,\ldots ,K} \pi _k f_{k}(\varvec{x}^{\text {obs}}_i) /\left( \sum _{\ell =1}^K \pi _\ell f_{\ell }(\varvec{x}^{\text {obs}}_i) \right)$.

The strong ignorability for clustering requires that the missingness mechanism can be ignored for evaluating, for each individual, the posterior probabilities of each class given the observed values, while weak ignorability is only concerned with the value of k which maximizes those probabilities, i.e., the class to which the individual will be affected via the maximum a posteriori rule.

As illustrated by the following example, ignorability for likelihood inference implies strong ignorability that implies weak ignorability.

Example 2

(Ignorability for clustering) Let $\varvec{X}_i\in {\mathbb {R}}$ be a univariate random variable following a bi-component Gaussian mixture model with equal proportions where the means of the components are $\mu _1$ and $\mu _2$ and where the variances of the components are equal to one. The conditional probability of observing $\varvec{X}_i$ is given by ${\mathbb {P}}(\varvec{R}_i=1\mid \varvec{X}_i,Z_{ik}=1)=(2\pi )^{-1/2}\exp (-(\varvec{x}_i-\mu _k)^2/2)$. The probability of observing the realization of the random variable depends on the value of this realization itself. Thus, the missingness mechanism is not MAR and thus the mechanism is not ignorable for likelihood inference. Noting that ${\mathbb {P}}(\varvec{R}_i=1\mid Z_{ik}=1)=(2\sqrt{\pi })^{-1}$, the following probabilities of classification computed on the observed data by considering the missingness mechanism are

$$\begin{aligned}{} & {} {\mathbb {P}}(Z_{ik}=1\mid \varvec{X}_i=\varvec{x}_i,\varvec{R}_i=\varvec{r}_i)\\{} & {} \quad =\left\{ \begin{array}{ll} \frac{\exp (-(\varvec{x}_i-\mu _k)^2)}{\exp (-(\varvec{x}_i-\mu _1)^2)+\exp (-(\varvec{x}_i-\mu _2)^2)} &{}\quad \text {if } \varvec{r}_i=1\\ 1/2 &{} \quad \text {if } \varvec{r}_i=0\\ \end{array} \right. . \end{aligned}$$

Moreover, the observed values $\varvec{X}^{\text {obs}}_i$ follow bi-component Gaussian mixture with equal proportions where the means of the components are $\mu _1$ and $\mu _2$ and where the variances of the components are equal to $2^{-1/2}$. Depending on the restrictions made on the model, the mechanism is either weakly or strongly ignorable:

• if the clustering is achieved by considering a bi-component Gaussian mixture with equal proportions and variances within component equal to 1, then the classification rule can be consistently estimated but not the probabilities of classification. Therefore, the mechanism is weakly ignorable.

• if the clustering is achieved by considering a bi-component Gaussian mixture with equal proportions and equal variances within component, then the probabilities of classification computed on the observed data are equal if the mechanism is ignored or not. This implies that the mechanism is strongly ignorable for clustering.

The previous example shows that weak ignorability can be obtained when the missingness mechanism leads to a misspecification of the model used to fit the densities $f_k(\varvec{x}^{\text {obs}}_i)$. In this particular example, the specific choice of the missingness mechanism leads the distributions of $\varvec{X}_i$ and $\varvec{X}^{\text {obs}}_i$ within components to belong to the same parametric family (i.e., Gaussian distribution). In the general case, these distributions do not belong to the same parametric family.

Despite the fact that, for some specific distributions, the ignorability for clustering could hold under MNAR scenarios, in the general case, the missingness mechanism has to be taken into account for clustering under the MNAR scenarios. Moreover, in the general case, the conditional distributions of $\varvec{X}_i\mid Z_{ik}$ and $\varvec{X}_i\mid Z_{ik}, \varvec{R}_i=1$ do not belong to the same parametric family. As an example, if $\varvec{X}_i\in {\mathbb {R}}$ follows a Gaussian mixture and if the conditional probability of $\varvec{R}_i \mid (\varvec{X}_i,\varvec{Z}_i)$ is given by a logistic link, then the conditional distribution of the observed values within a cluster (e.g., $\varvec{X}_i\mid Z_{ik},\varvec{R}_i=1$) is no longer Gaussian. Thus, if parametric assumptions are relevant for modelling the distribution within clusters, it is more hazardous to make parametric assumptions for the distribution of the observed values within the clusters. Thus, the next section introduces a semi-parametric mixture model that takes into account the missingness process to also handle mechanisms that are not ignorable for clustering.

Appendix B: Proofs

Proof of Lemma 1

The model defined by (8)–(9) is identifiable, if

$$\begin{aligned} \forall (\varvec{r}_i^\top ,\varvec{x}_i^\top )^\top \in \{0,1\}^d\times {\mathbb {R}}^d,\; g^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })=g^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\tilde{\varvec{\theta }}) \; \Rightarrow \; \varvec{\theta }={\tilde{\varvec{\theta }}},\quad \end{aligned}$$

(B1)

where $\varvec{\theta }$ groups the finite dimensional parameters, $\pi _k$ and $\tau _{kj}$, and the infinite dimensional parameters $p_{kj}$ for $k=1,\ldots ,K$ and $j=1,\ldots ,d$. Thus, considering the case where all the variables are observed (i.e, $r_{ij}=1$, for $j=1,\ldots ,d$), the left hand side of (B1) implies

$$\begin{aligned} \forall \varvec{x}_i\in {\mathbb {R}}^d,\; \sum _{k=1}^K \rho _k\prod _{j=1}^dp_{kj}(x_{ij})=\sum _{k=1}^K {\tilde{\rho }}_k\prod _{j=1}^d{\tilde{p}}_{kj}(x_{ij}), \end{aligned}$$

(B2)

where $\rho _k=\pi _k\prod _{j=1}^d \tau _{kj}$ and ${\tilde{\rho }}_k={\tilde{\pi }}_k\prod _{j=1}^d {\tilde{\tau }}_{kj}$. Theorem 8 in Allman et al. (2009) states that a mixture whose components are defined as a product of univariate densities, is identifiable if all the univariate densities are linearly independent and if $d\ge 3$, up to label swapping. Thus, under the conditions of Lemma 1, Theorem 8 in Allman et al. (2009) implies that $\forall k=1,\ldots ,K,$ and $\forall j=1,\ldots ,d,$

$$\begin{aligned} \pi _k\prod _{j=1}^d\tau _{kj}={\tilde{\pi }}_k\prod _{j=1}^d{\tilde{\tau }}_{kj} \text { and } p_{kj}={\tilde{p}}_{kj}. \end{aligned}$$

(B3)

This results and the left-hand side of (B1) imply

$$\begin{aligned}{} & {} \forall (\varvec{r}_i^\top ,\varvec{x}_i^\top )^\top \in \{0,1\}^d\times {\mathbb {R}}^d,\; \sum _{k=1}^K \pi _k\prod _{j=1}^d\tau _{kj}\prod _{j=1}^d p_{kj}(x_{ij})\nonumber \\ {}{} & {} \quad = \sum _{k=1}^K {\tilde{\pi }}_k\prod _{j=1}^d{\tilde{\tau }}_{kj}\prod _{j=1}^d p_{kj}(x_{ij}). \end{aligned}$$

(B4)

Considering the marginal distribution of $(r_{ij}^\top ,x_{ij}^\top )^\top$ with $r_{ij}=1$, for any $j=1,\ldots ,d$, we have from (B4)

$$\begin{aligned} \forall j=1,\ldots ,d,\; \forall x_{ij}\in {\mathbb {R}},\; \sum _{k=1}^K (\pi _k\tau _{kj} - {\tilde{\pi }}_k{\tilde{\tau }}_{kj}) p_{kj}(x_{ij}) =0. \end{aligned}$$

(B5)

The densities $p_{kj}$ are linearly independent, so $\forall (\alpha _1,\ldots ,\alpha _K)^\top \in {\mathbb {R}}^K{\setminus }\{{\varvec{0}}\}$, $\sum _{k=1}^K \alpha _k p_{kj}$ is not the zero function. Thus, (B5) implies that

$$\begin{aligned} \forall j=1,\ldots ,d,\; \forall k=1,\ldots ,K,\; \pi _k\tau _{kj} = {\tilde{\pi }}_k{\tilde{\tau }}_{kj}. \end{aligned}$$

(B6)

From (B3) and (B6), recalling that $\pi _k>0$ and $\tau _{kj}>0$, for $k=1,\ldots ,K$ and $j=1,\ldots ,d$, we obtain that

$$\begin{aligned} \forall k=1,\ldots ,K,\; {\varvec{M}}{\varvec{u}}_k={\varvec{0}}, \end{aligned}$$

where

$$\begin{aligned} {\varvec{M}}=\begin{bmatrix} 1 &{} \quad 1 &{} \quad \ldots &{}\quad 1 \\ 1 &{} &{}&{} \\ \vdots &{} \quad &{}\quad {\varvec{I}}_d &{} \\ 1 &{} &{} &{} \end{bmatrix}\quad \text {and}\quad {\varvec{u}}_k= \begin{bmatrix} \log (\pi _k/{\tilde{\pi }}_k)\\ \log (\tau _{k1}/{\tilde{\tau }}_{k1})\\ \vdots \\ \log (\tau _{kd}/{\tilde{\tau }}_{kd})\\ \end{bmatrix}, \end{aligned}$$

where ${\varvec{I}}_d$ is the identity matrix of size d. As ${\varvec{M}}$ has full rank for $d\ge 2$, we deduce that ${\varvec{u}}_k={\varvec{0}}$ and thus $\pi _k={\tilde{\pi }}_k$ and $\tau _{kj}={\tilde{\tau }}_{kj}$ for $k=1,\ldots ,K$ and $j=1,\ldots ,d$. $\square$

Proof of Lemma 2

This proof is similar to the proof of Theorem 1 of (Levine et al. 2011) and is only given for ease of reading. We have

$$\begin{aligned} \ell _n(\varvec{\theta }) - \ell _n(\varvec{\theta }^{[r]}) \ge b^{[r]}(\varvec{\theta }) - b^{[r]}(\varvec{\theta }^{[r]}), \end{aligned}$$

where $b^{[r]}(\varvec{\theta })=\sum _{i=1}^n \sum _{k=1}^K t_{ik}(\varvec{\theta }^{[r]})\ln \left( \pi _k{\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })\right)$. Indeed, using the concavity of the logarithm,

$$\begin{aligned} \ell _n(\varvec{\theta }) - \ell _n(\varvec{\theta }^{[r]})&= \sum _{i=1}^n \ln \left( \sum _{k=1}^K t_{ik}(\varvec{\theta }^{[r]}) \frac{\pi _k{\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })}{\pi _k^{[r]}{\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta }^{[r]})}\right) \\&\ge \sum _{i=1}^n \sum _{k=1}^K t_{ik}(\varvec{\theta }^{[r]}) \ln \frac{\pi _k{\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })}{\pi _k^{[r]}{\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta }^{[r]})}\\&= b^{[r]}(\varvec{\theta }) - b^{[r]}(\varvec{\theta }^{[r]}). \end{aligned}$$

To prove the monotonicity of the algorithm, it suffices to show that $\varvec{\theta }^{[r+1]}$ is such that $b^{[r]}(\varvec{\theta }) - b^{[r]}(\varvec{\theta }^{[r]})\ge 0$. Note that the following decomposition holds

$$\begin{aligned} b^{[r]}(\varvec{\theta }) = b_1^{[r]}(\varvec{\theta }) +b_2^{[r]}(\varvec{\theta }) + \sum _{k=1}^K\sum _{j=1}^d b^{[r]}_{3kj}(\varvec{\theta }) \end{aligned}$$

where

$$\begin{aligned} b_1^{[r]}(\varvec{\theta }) = \sum _{i=1}^n \sum _{k=1}^K t_{ik}(\varvec{\theta }^{[r]})\ln \pi _k, \end{aligned}$$

$$\begin{aligned} b_2^{[r]}(\varvec{\theta }) = \sum _{i=1}^n \sum _{k=1}^K \sum _{j=1}^d t_{ik}(\varvec{\theta }^{[r]}) \left( r_{ij}\ln \tau _{kj} + (1-r_{ij})\ln (1-\tau _{kj})\right) , \end{aligned}$$

and

$$\begin{aligned} b_{3kj}^{[r]}(\varvec{\theta })=\int _{\Omega _j} \sum _{i=1}^n t_{ik}(\varvec{\theta }^{[r]}) r_{ij} \frac{1}{h_j} {\mathcal {K}}\left( \frac{x_{ij} - u}{h_j}\right) \ln p_{kj}(u)du. \end{aligned}$$

Maximizing $b^{[r]}(\varvec{\theta })$ on the proportions $\pi _1,\ldots ,\pi _k$ is equivalent to maximizing $b_1^{[r]}(\varvec{\theta })$ on the proportions. Similarly, maximizing $b^{[r]}(\varvec{\theta })$ on the probabilities $\tau _{kj}$ is equivalent to maximizing $b_2^{[r]}(\varvec{\theta })$ on the $\tau _{kj}$’s. Thus, one can check that the estimators $\pi _k^{[r+1]}$’s and $\tau _{kj}^{[r+1]}$’s maximize $b^{[r]}(\varvec{\theta })$ on the $\pi _k$’s and on the $\tau _{kj}$’s. Finally, note that we have

$$\begin{aligned} b_{3kj}^{[r]}(\varvec{\theta })= & {} - c_{kj}^{[r]} \int _{\Omega _j} p^{[r+1]}_{kj}(u) \ln \frac{p^{[r+1]}_{kj}(u)}{p_{kj}(u)} du \\ {}{} & {} +c_{kj}^{[r]} \int _{\Omega _j} p^{[r+1]}_{kj}(u) \ln p^{[r+1]}_{kj}(u) du, \end{aligned}$$

where $c_{kj}^{[r]}=\sum _{i=1}^n t_{ik}(\varvec{\theta }^{[r]})r_{ij}$. The second term of the right-hand side of the equation does not depend on $p_{kj}$. The first term of the right-hand side of the equation is based on Kullback–Leibler divergence from $p_{kj}$ to $p^{[r+1]}_{kj}$. Thus, noting that $c_{kj}^{[r]}>0$ (because $\sum _{i=1}^n r_{ij}\ge 1$), $p^{[r+1]}_{kj}$ is unique, up to changes on a set of Lebesgue measure zero, density function maximizing $b_{3kj}^{[r]}(\varvec{\theta })$. This proof is concluded by noting that $\varvec{\theta }^{[r+1]}={{\,\mathrm{arg\,max}\,}}_{\varvec{\theta }}b^{[r]}(\varvec{\theta })$ yielding that $b^{[r]}(\varvec{\theta }^{[r+1]}) \ge b^{[r]}(\varvec{\theta }^{[r]})$ and thus $\ell _n(\varvec{\theta }^{[r+1]}) \ge \ell _n(\varvec{\theta }^{[r]})$. $\square$

Appendix C: Simulation

This section gives the values of $\delta$ and $\gamma$ used during the different experiments. These values have been estimated by generating a large sample ($n=10^4$ observations) where the theoretical rate of misclassification is computed between the true partition and the partition given by the rule of the maximum a posteriori computed with the true parameters. Note that, because the missingness mechanism impacts the distribution within components and thus the overlaps between components, changing the value of $\gamma$ allows the rate of missingness to be changed but implies changing the value of $\delta$ to retain the rate of misclassification. Table 6 presents the parameters used to generate the data with $K=3$ components, $d=6$ variables and a theoretical rate of misclassification of $10\%$ (related to Figs. 1 and  3).

Table 6 Parameters used to generate the data with $K=3$ components, $d=6$ variables and a theoretical misclassification of $10\%$

Table 7 presents the parameters used to generate the data with $K=3$ components, a theoretical missing rate per variable of $30\%$ and a theoretical rate of misclassification of $10\%$ (related to Fig. 4).

Table 7 Parameters used to generate the data with $K=3$ components, a theoretical missing rate per variable of $30\%$ and a theoretical rate of misclassification of $10\%$

Table 8 presents the parameters used to generate the data with $K=3$ components, $d=6$ variables and a theoretical missing rate per variable of $30\%$ (related to Fig. 5).

Table 8 Parameters used to generate the data with $K=3$ components, $d=6$ variables and a theoretical missing rate per variable of $30\%$

Appendix D: Extension of the approach to mixed-type data

This sections considers that $\varvec{X}_i=(\varvec{X}_{i1}^\top ,\ldots ,\varvec{X}_{i d_c}^\top ,X_{i d_c +1},\ldots ,X_{id})^\top$ is a d-variate vector of mixed-type such that the first $d_c$ elements are categorical and the last $d-d_c$ elements are continuous. Each categorical variable $\varvec{X}_{ij}=(X_{ij1},\ldots , X_{ij m_j } )^\top$, with $1 \le j \le d_c$, has $m_j$ levels and $X_{ijh}=1$ if subject i takes level h for variable j and $X_{ijh}=0$ otherwise. The definition of $\varvec{R}_i=(R_{i1},\ldots ,R_{id})^\top$ is unchanged yielding that $R_{ij}=1$ if variable j is observed for subject i and $R_{ij}=0$ otherwise. Similarly to Sect. 2, we consider that the couples $(X_{ij},R_{ij})^\top$ are conditionally independent given $\varvec{Z}_i$. Thus, the conditional distribution of $R_{ij}$ given $Z_{ij}=1$ is a Bernoulli distribution with parameter $\tau _{kj}$. The conditional distributions of a continuous variable (i.e., $X_{ij}$ with $d_c + 1 \le j \le d$) given $Z_{ij}=1$ and $R_{ij}=1$ and given $Z_{ij}=1$ and $R_{ij}=0$ are defined by the densities $p_{kj}$ and $q_{kj}$ respectively. Finally, the conditional distributions of a categorical variable (i.e., $X_{ij}$ with $d_c \le j \le d_c$) given $Z_{ij}=1$ and $R_{ij}=1$ and given $Z_{ij}=1$ and $R_{ij}=0$ are defined by two multinomial distributions. We denote $\beta _{kj}=(\beta _{kj1},\ldots ,\beta _{kjm_j})^\top$, the vector defining the multinomial distribution of variable j (with $1\le j \le d_c$) given $Z_{ik}=1$ and $R_{ij}=1$. Thus, $0<\beta _{kjh}$ is the probability that subject i takes level h for variable j under component k when this variable is observed and $\sum _{h=1}^{m_j} \beta _{kjh}=1$. Similarly to (8) and (9), the distribution of the observed variables $(\varvec{x}_i^{\text {obs}\top },\varvec{r}_i^\top )^\top$

$$\begin{aligned} g^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })=\sum _{k=1}^K \pi _k g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta }), \end{aligned}$$

where the pdf of component k is a specific version of (3) defined by

$$\begin{aligned} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta })= & {} \left( \prod _{j=1}^d\tau _{kj}^{r_{ij}}(1-\tau _{kj})^{1-r_{ij}} \right) \left( \prod _{j=1}^{d_c} \prod _{h=1}^{m_j} \left( \beta _{kjh}\right) ^{r_{ij} x_{ijh}}\right) \\{} & {} \quad \times \left( \prod _{j=d_c+1}^d p_{kj}^{r_{ij}}(x_{ij})\right) . \end{aligned}$$

Thus, a sufficient condition for identifiability, is that the proportions are not zero (i.e., $0<\pi _k$ for any k) and that there are at least three continuous variables (i.e., $d-d_c\ge 3$) that have linearly dependent densities $p_{kj}$ and non-zero probability of observing these variables under each component (i.e., $0<\tau _{kj}$).

Fig. 9

Maximum of the smoothed log-likelihood with respect to the number of clusters

The estimation is performed by maximizing the smoothed log-likelihood where the smoothing is performed only on the continuous variables. Therefore, the smoothed log-likelihood function is defined by

$$\begin{aligned} \ell _n(\varvec{\theta })= & {} \sum _{i=1}^n \ln \left[ \sum _{k=1}^K \pi _k \left( \prod _{j=1}^d\tau _{kj}^{r_{ij}} (1-\tau _{kj})^{1-r_{ij}} \right) \left( \prod _{j=1}^{d_c} \prod _{h=1}^{m_j} \left( \beta _{kjh}\right) ^{r_{ij} x_{ijh}} \right) \right. \times \\{} & {} \quad \left. \left( \prod _{j=d_c+1}^d \exp \left( r_{ij} {\mathcal {S}} \ln p_{kj}(x_{ij}) \right) \right) \right] . \end{aligned}$$

In this context, the maximization on $\varvec{\theta }$ of the smoothed log-likelihood function is performed via an MM algorithm. This iterative algorithm starts at the initial value of the parameters $\varvec{\theta }^{[0]}$. At iteration [r], it performs the following two steps

• Computing the smoothed probabilities of subpopulation memberships

  $$\begin{aligned} t_{ik}(\varvec{\theta }^{[r]}) = \frac{\pi _k^{[r]} {\mathcal {N}} g_k^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta }^{[r]})}{\sum _{\ell =1}^K\pi _\ell ^{[r]} {\mathcal {N}} g_\ell ^{\text {obs}}(\varvec{x}^{\text {obs}}_i,\varvec{r}_i;\varvec{\theta }^{[r]})}. \end{aligned}$$

• Updating the estimators

  • Updating of the proportions

    $$\begin{aligned} \pi _k^{[r+1]}= \frac{1}{n} \sum _{i=1}^n t_{ik}(\varvec{\theta }^{[r]}). \end{aligned}$$

  • Updating of the parameters of the missingness mechanism

    $$\begin{aligned} \tau _{kj}^{[r+1]} = \frac{\sum _{i=1}^n r_{ij} t_{ik}(\varvec{\theta }^{[r]})}{\sum _{i=1}^n t_{ik}(\varvec{\theta }^{[r]})}. \end{aligned}$$

  • Updating of the parameters of the categorical variables, for $1\le j \le d_c$

    $$\begin{aligned} \beta _{kjh}^{[r+1]} = \frac{\sum _{i=1}^n r_{ij} x_{ijh} t_{ik}(\varvec{\theta }^{[r]})}{\sum _{i=1}^n r_{ij} t_{ik}(\varvec{\theta }^{[r]})}. \end{aligned}$$

  • Updating of the conditional distribution of the continuous variables, for $d_c+1\le j \le d$

    $$\begin{aligned} p_{kj}^{[r+1]}(u) = \frac{\sum _{i=1}^n r_{ij} t_{ik}(\varvec{\theta }^{[r]}) \frac{1}{h_j}{\mathcal {K}}\left( \frac{x_{ij} - u}{h_j}\right) }{\sum _{i=1}^n r_{ij} t_{ik}(\varvec{\theta }^{[r]})}. \end{aligned}$$

Fig. 10

Rate of missingness and empirical counterpart of ${\mathbb {E}}[\max _k \ln {\mathbb {P}}(Z_{ik}\mid R_{ij})]$ for each variable

Appendix E: Echocardiogram Data Set

Figure 9 helps in selecting a suitable number of components by presenting the evolution of the maximum smoothed log-likelihood with respect to the number of clusters

Fig. 11

Rate of missingness and empirical counterpart of ${\mathbb {E}}[\max _k \ln {\mathbb {P}}(Z_{ik}\mid X_{ij})]$ for each variable

Table 9 p-values obtained for testing the significance of the correlation coefficient computed on the couple of observed variables for observations assigned to cluster 1

Table 10 p-values obtained for testing the significance of the correlation coefficient computed on the couple of observed variables for observations assigned to cluster 3

Figure 10 illustrates the relation between the rate of missingness and how the missingness of a variable is informative for the partition. Moreover, Fig. 11 illustrates the relation between the rate of missingness and how the observed values of a variable is informative for the partition.

Tables 9 and 10 presents the p-values obtained by testing the nullity of the correlation coefficient of the conditional distribution of couple of variables conditionally on component 1 and 3 respectively. The high values of the p-values suggest that the assumption of conditional independence given the component membership is appropriate. Note that results related to component 2 are not presented due to a lack of subjects affected to this class.