# Clustering : une vision unifiée pour une utilisation éclairée - Partie 3 : Formalisation par modèles de mélange

C Biernacki

► To cite this version:

C Biernacki. Clustering : une vision unifiée pour une utilisation éclairée - Partie 3 : Formalisation par modèles de mélange. Doctoral. France. 2023. hal-04370761

HAL Id: hal-04370761

https://inria.hal.science/hal-04370761v1

Submitted on 3 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Clustering : une vision unifiée pour une utilisation éclairée

—

Partie 3
Formalisation par modèles de mélange

C. Biernacki

Ateliers statistiques de la SFdS, 15 et 16 juin 2023, IHP, Paris

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

Université
de Lille
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

# Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
Clustering: an ill-posed problem

■ We have seen how to perform and evaluate clustering...

■ ... but we do not know what is a cluster

■ Thus we have built something without defining it!

It is a serious problem

Impossible to provide guarantees on by-products of clustering (ex.: some user decisions) since no guarantees on clustering itself is really available

$$x \rightarrow \hat{z} = f(x) \rightarrow \underbrace{\text{decision}}_{\text{need guarantees}} = g(\hat{z})$$
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

Expected guarantees on clustering

■ $\hat{z}, \hat{\mu}, \hat{K}, \text{etc.}$ are estimates of theoretical quantities $z, \mu, K, \text{etc.}$

■ It can be thus expected classical guarantees provided for any estimate in statistics

■ Typically: consistency, bias, variance

■ Examples:

$$p(\hat{K} = K) = 1 \quad \text{as } n \to \infty$$

$$p(\hat{z} = z) \text{ for finite } n$$

■ In the previous lessons we were very far from such a requirement...

Key idea

Formalize the rigorous definition of a cluster
The model-based clustering paradigm

a cluster ⇐⇒ a distribution

■ It recasts all previous/next questions into model design/estimation/selection

■ It takes benefits from all theoretical statistics environment

■ How to choose the best metric $M_{(k)}$?

■ How to choose the number K of clusters?

■ Clusters of different sizes are they well estimated?

■ How to choose the data unit?

■ How to select features?

■ How to deal with mixed data?

■ How to deal with missing data?

■ How to deal with outliers?

■ ...
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

What about empirical clustering?

Somewhere it works pretty well even if it has the previous mentioned limits

1. Interesting to understand why

2. Interesting to overcome their limits then

In fact many empirical methods are hidden model-based clustering ones!
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

Reformulate K-means: the hidden Gaussian assumption

$$
\begin{align*}
W_I(\mathbf{z}) &= \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \|\mathbf{x}_i - \boldsymbol{\mu}_k\|_1^2 \\
&= -2 \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \ln \left[ \underbrace{\frac{1}{K} \frac{1}{(2\pi)^{d/2} \|\mathbf{I}\|^{1/2}} \exp \left( -\frac{1}{2} (\mathbf{x}_i - \boldsymbol{\mu}_k)' (\mathbf{x}_i - \boldsymbol{\mu}_k) \right)}_{N_d(\boldsymbol{\mu}_k, \mathbf{I})} \right] + \text{cst}
\end{align*}
$$

Model

d-variate Gaussian with variance matrix I and same cluster sample size (see later)
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

Reformulate K-means: the hidden estimate choice

$$W_1(z) = -2\ell_c(\mu; x, z) + \text{cst}$$

<figure><img src="image_41.png" /></figure>

$$\hat{z}^W = \arg \min_z W_1(z) \\ (\hat{z}^{\ell_c}, \hat{\mu}^{\ell_c}) = \arg \max_{(z, \mu)} \ell_c(\mu; x, z) \} \Rightarrow \hat{z}^W \equiv \hat{z}^{\ell_c}$$

Estimate

Maximum of the so-called complete-likelihood (see later for its statistical properties)
<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

## Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

Parametric mixture model

■ Parametric assumption:

thus

$$p_k(\mathbf{x}_1) = p(\mathbf{x}_1; \alpha_k)$$

$$p(\mathbf{x}_1) = p(\mathbf{x}_1; \theta) = \sum_{k=1}^{K} \pi_k p(\mathbf{x}_1; \alpha_k)$$

■ Mixture parameter:

$$\theta = (\pi, \alpha) \text{ with } \alpha = (\alpha_1, \dots, \alpha_K)$$

■ Model: it includes both the family $p(\cdot; \alpha_k)$ and the number of groups $K$

$$\mathbf{m} = \{p(\mathbf{x}_1; \theta) : \theta \in \Theta\}$$

The number of free continuous parameters is given by

$$\nu = \dim(\Theta)$$
The clustering process in mixtures

1. Estimation of $\theta$ by $\hat{\theta}$

2. Estimation of the conditional probability that $x_i \in G_k$

$$t_{ik}(\hat{\theta}) = p(Z_{ik} = 1 | \mathbf{X}_i = x_i; \hat{\theta}) = \frac{\hat{\pi}_k p(\mathbf{x}_i; \hat{\alpha}_k)}{p(\mathbf{x}_i; \hat{\theta})}$$

3. Estimation of $z_i$ by maximum a posteriori (MAP)

$$\hat{z}_{ik} = \mathbb{I}_{\{k=\arg\max_{h=1,...,K} t_{ih}(\hat{\theta})\}}$$

<figure><img src="image_54.png" /></figure>
<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

Geometric interpretation of $\Sigma_k$

$$ \Sigma_k = \underbrace{\lambda_k}_{\text{volume}} \cdot \underbrace{\mathbf{D}_k}_{\text{orientation}} \cdot \underbrace{\mathbf{A}_k}_{\text{shape}} \cdot \mathbf{D}_k' $$

<figure><img src="image_65.png" /></figure>
<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

Estimation of $\theta$ by complete-likelihood

Maximize the complete-likelihood over $(\theta, z)$

$$ \ell_c(\theta; \mathbf{x}, \mathbf{z}) = \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \ln \{ \pi_k p(\mathbf{x}_j; \alpha_k) \} $$

* Equivalent to traditional methods

<table><thead><tr><th>Metric</th><th>M = I</th><th>M free</th><th>M<sub>k</sub> free</th></tr></thead><tbody><tr><td>Gaussian model</td><td>[πλI]</td><td>[πλC]</td><td>[πλ<sub>k</sub>C<sub>k</sub>]</td></tr></tbody></table>

* Bias of $\hat{\theta}$: heavy if poor separated clusters

* Associated optimization algorithm: CEM (see later)

* CEM with $[\pi\lambda I]$ is strictly equivalent to $K$-means

* CEM is simple et fast (convergence with few iterations)
<figure><img src="image_72.png" /></figure>

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /></figure>

<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

<figure><img src="image_77.png" /></figure>

Estimation of $\theta$ by observe-likelihood

Maximize the observe-likelihood on $\theta$

$$ \ell(\theta; \mathbf{x}) = \sum_{i=1}^{n} \ln p(\mathbf{x}_i; \theta) $$

* Convergence of $\hat{\theta}$, asymptotic efficiency, asymptotically unbiased

* General algorithm for missing data: EM

* EM is simple but slower than CEM

* Interpretation: it is a kind of fuzzy clustering
<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>

<figure><img src="image_84.png" /></figure>

<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /></figure>

## Principle of EM and CEM

* Initialization: $\theta^0$

* Iteration $n^q$:

* Step E: estimate probabilities $t^q = \{t_{ik}(\theta^q)\}$

* Step C: classify by setting $t^q = \text{MAP}(\{t_{ik}(\theta^q)\})$

* Step M: maximize $\theta^{q+1} = \arg\max_{\theta} \ell_c(\theta; x, t^q)$

* Stopping rule: iteration number or criterion stability

### Properties

* $\oplus$: simplicity, monotony, low memory requirement

* $\ominus$: local maxima (depends on $\theta^0$), linear convergence (EM)
<figure><img src="image_88.png" /></figure>

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /></figure>

<figure><img src="image_91.png" /></figure>

<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

Gaussian M-step

$$n_k^{(q)} = \sum_{i=n'+1}^{n} t_{ik}(\theta^{(q)})$$

$$\pi_k^{(q+1)} = \frac{n_k^{(q)}}{n}$$

$$\mu_k^{(q+1)} = \frac{1}{n_k^{(q)}} \left( \sum_{i=1}^{n} t_{ik}(\theta^{(q)}) x_i \right)$$

$$\Sigma_k^{(q+1)} = \frac{1}{n_k^{(q)}} \left( \sum_{i=1}^{n} t_{ik}(\theta^{(q)})(x_i - \mu_k^{(q+1)})(x_i - \mu_k^{(q+1)})' \right)$$
<figure><img src="image_94.png" /></figure>

<figure><img src="image_95.png" /></figure>

<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /></figure>

<figure><img src="image_98.png" /></figure>

<figure><img src="image_99.png" /></figure>

<figure><img src="image_100.png" /></figure>

<figure><img src="image_101.png" /></figure>

<figure><img src="image_102.png" /></figure>

<figure><img src="image_103.png" /></figure>

Example of an EM in the univariate case

<figure><img src="image_104.png" /></figure>

Note: low at the beginning but increase of the log-likelihood
<figure><img src="image_105.png" /></figure>

<figure><img src="image_106.png" /></figure>

<figure><img src="image_107.png" /></figure>

<figure><img src="image_108.png" /></figure>

<figure><img src="image_109.png" /></figure>

<figure><img src="image_110.png" /></figure>

Local maxima

<figure><img src="image_111.png" /></figure>

<figure><img src="image_112.png" /></figure>
<figure><img src="image_113.png" /></figure>

<figure><img src="image_114.png" /></figure>

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>

<figure><img src="image_118.png" /></figure>

Comparison EM-CEM in practice

<table><thead><tr><th rowspan="2">n</th><th rowspan="2">overlap<fnref n="1" /></th><th colspan="2">K̂L(θ, ̂θ)</th><th colspan="2">err(z, ŷ)</th></tr><tr><th>EM</th><th>CEM</th><th>EM</th><th>CEM</th></tr></thead><tbody><tr><td rowspan="3">20</td><td>low</td><td>0.2770</td><td>0.2771</td><td>0.3383</td><td>0.3217</td></tr><tr><td>middle</td><td>0.4916</td><td>0.3699</td><td>0.2050</td><td>0.1700</td></tr><tr><td>high</td><td>0.4108</td><td>0.3132</td><td>0.0983</td><td>0.0667</td></tr><tr><td rowspan="3">200</td><td>low</td><td>0.0209</td><td>0.0822</td><td>0.3342</td><td>0.3188</td></tr><tr><td>middle</td><td>0.0187</td><td>0.0425</td><td>0.1638</td><td>0.1587</td></tr><tr><td>high</td><td>0.0172</td><td>0.0209</td><td>0.0530</td><td>0.0500</td></tr><tr><td rowspan="3">2000</td><td>low</td><td>0.0014</td><td>0.0454</td><td>0.3112</td><td>0.3113</td></tr><tr><td>middle</td><td>0.0017</td><td>0.0246</td><td>0.1620</td><td>0.1619</td></tr><tr><td>high</td><td>0.0017</td><td>0.0059</td><td>0.0509</td><td>0.0510</td></tr></tbody></table>

<footnote n="1">¹high: 30%, middle: 15%, low:5%</footnote>
<figure><img src="image_119.png" /></figure>

<figure><img src="image_120.png" /></figure>

<figure><img src="image_121.png" /></figure>

<figure><img src="image_122.png" /></figure>

<figure><img src="image_123.png" /></figure>

<figure><img src="image_124.png" /></figure>

Categorical variables: latent class model

■ Categorical variables: $d$ variables with $m_j$ modalities each, $x_i^j \in \{0, 1\}^{m_j}$ and
$x_i^{jh} = 1 \Leftrightarrow$ variable $j$ of $x_i$ takes level $h$

■ Conditional independence:

$$p(\mathbf{x}_i; \boldsymbol{\alpha}_k) = \prod_{j=1}^{d} \prod_{h=1}^{m_j} (\alpha_k^{jh})^{x_i^jh}$$

and

$$\alpha_k^{jh} = p(\mathbf{x}_i^{jh} = 1 | z_{ik} = 1)$$

$$\text{with } \boldsymbol{\alpha}_k = (\alpha_k^{jh}, j = 1, \dots, d; h = 1, \dots, m_j)$$
<figure><img src="image_125.png" /></figure>

<figure><img src="image_126.png" /></figure>

<figure><img src="image_127.png" /></figure>

<figure><img src="image_128.png" /></figure>

<figure><img src="image_129.png" /></figure>

<figure><img src="image_130.png" /></figure>

Integer: Poisson mixture model

■ integer variables: $d$ variables $\mathbf{x}_i^j \in \mathbb{N}$

■ Intra conditional independence:

$$p(\mathbf{x}_i^{int}; \alpha_k^{int}) = \prod_{j=1}^{d} \frac{(\alpha_k^j)^{\mathbf{x}_i^j}}{\alpha_k^j!} e^{-\alpha_k^j}$$
<figure><img src="image_131.png" /></figure>

<figure><img src="image_132.png" /></figure>

<figure><img src="image_133.png" /></figure>

<figure><img src="image_134.png" /></figure>

<figure><img src="image_135.png" /></figure>

<figure><img src="image_136.png" /></figure>

SPAM E-mail Database<fnref n="3" />

■ $n = 4601$ e-mails composed by 1813 "spams" and 2788 "good e-mails"

■ $d = 48 + 6 = 54$ continuous descriptors<fnref n="2" />

■ 48 percentages that a given word appears in an e-mail ("make", "you"...)

■ 6 percentages that a given char appears in an e-mail (";", "$"...)

■ Transformation of continuous descriptors into binary descriptors

$$x_i^j = \begin{cases} 1 & \text{if word/char } j \text{ appears in e-mail } i \\ 0 & \text{otherwise} \end{cases}$$

<footnote n="2">²There are 3 other continuous descriptors we do not use</footnote>

<footnote n="3">³https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/</footnote>
<figure><img src="image_137.png" /></figure>

<figure><img src="image_138.png" /></figure>

<figure><img src="image_139.png" /></figure>

<figure><img src="image_140.png" /></figure>

<figure><img src="image_141.png" /></figure>

<figure><img src="image_142.png" /></figure>

An EM run with a binary data set

Initial binary data

<figure><img src="image_143.png" /></figure>
<figure><img src="image_144.png" /></figure>

<figure><img src="image_145.png" /></figure>

<figure><img src="image_146.png" /></figure>

<figure><img src="image_147.png" /></figure>

<figure><img src="image_148.png" /></figure>

<figure><img src="image_149.png" /></figure>

An EM run with a binary data set

Iteration 1

<figure><img src="image_150.png" /></figure>
<figure><img src="image_151.png" /></figure>

<figure><img src="image_152.png" /></figure>

<figure><img src="image_153.png" /></figure>

<figure><img src="image_154.png" /></figure>

<figure><img src="image_155.png" /></figure>

<figure><img src="image_156.png" /></figure>

An EM run with a binary data set

Iteration 2

<figure><img src="image_157.png" /></figure>
<figure><img src="image_158.png" /></figure>

<figure><img src="image_159.png" /></figure>

<figure><img src="image_160.png" /></figure>

<figure><img src="image_161.png" /></figure>

<figure><img src="image_162.png" /></figure>

<figure><img src="image_163.png" /></figure>

An EM run with a binary data set

Iteration 3

<figure><img src="image_164.png" /></figure>
<figure><img src="image_165.png" /></figure>

<figure><img src="image_166.png" /></figure>

<figure><img src="image_167.png" /></figure>

<figure><img src="image_168.png" /></figure>

<figure><img src="image_169.png" /></figure>

<figure><img src="image_170.png" /></figure>

An EM run with a binary data set

Iteration 4

<figure><img src="image_171.png" /></figure>
<figure><img src="image_172.png" /></figure>

<figure><img src="image_173.png" /></figure>

<figure><img src="image_174.png" /></figure>

<figure><img src="image_175.png" /></figure>

<figure><img src="image_176.png" /></figure>

<figure><img src="image_177.png" /></figure>

An EM run with a binary data set

Iteration 5

<figure><img src="image_178.png" /></figure>
<figure><img src="image_179.png" /></figure>

<figure><img src="image_180.png" /></figure>

<figure><img src="image_181.png" /></figure>

<figure><img src="image_182.png" /></figure>

<figure><img src="image_183.png" /></figure>

<figure><img src="image_184.png" /></figure>

An EM run with a binary data set

Iteration 6

<figure><img src="image_185.png" /></figure>
<figure><img src="image_186.png" /></figure>

<figure><img src="image_187.png" /></figure>

<figure><img src="image_188.png" /></figure>

<figure><img src="image_189.png" /></figure>

<figure><img src="image_190.png" /></figure>

<figure><img src="image_191.png" /></figure>

An EM run with a binary data set

Iteration 7

<figure><img src="image_192.png" /></figure>
Need to formalize

<figure><img src="image_193.png" /></figure>

Formalizing estimation

<figure><img src="image_194.png" /></figure>

Formalizing selection

<figure><img src="image_195.png" /></figure>

More advanced formalizing

<figure><img src="image_196.png" /></figure>

Experiments

<figure><img src="image_197.png" /></figure>

To go further

<figure><img src="image_198.png" /></figure>

An EM run with a binary data set

Iteration 8

<figure><img src="image_199.png" /></figure>
<figure><img src="image_200.png" /></figure>

<figure><img src="image_201.png" /></figure>

<figure><img src="image_202.png" /></figure>

<figure><img src="image_203.png" /></figure>

<figure><img src="image_204.png" /></figure>

<figure><img src="image_205.png" /></figure>

An EM run with a binary data set

<figure><img src="image_206.png" /><figcaption>Iteration 9</figcaption></figure>
<figure><img src="image_207.png" /></figure>

<figure><img src="image_208.png" /></figure>

<figure><img src="image_209.png" /></figure>

<figure><img src="image_210.png" /></figure>

<figure><img src="image_211.png" /></figure>

<figure><img src="image_212.png" /></figure>

An EM run with a binary data set

Iteration 10

<figure><img src="image_213.png" /></figure>
<figure><img src="image_214.png" /></figure>

Formalizing estimation

<figure><img src="image_215.png" /></figure>

Formalizing selection

<figure><img src="image_216.png" /></figure>

More advanced formalizing

<figure><img src="image_217.png" /></figure>

Experiments

<figure><img src="image_218.png" /></figure>

To go further

<figure><img src="image_219.png" /></figure>

An EM run with a binary data set

<figure><img src="image_220.png" /></figure>
<figure><img src="image_221.png" /></figure>

<figure><img src="image_222.png" /></figure>

<figure><img src="image_223.png" /></figure>

<figure><img src="image_224.png" /></figure>

<figure><img src="image_225.png" /></figure>

<figure><img src="image_226.png" /></figure>

# Mixed data: classical approaches

Usually, unify data type by transformation :

* Quantify continuous variables: **loose some information**

* MCA dof categorical variable: **loose the meaning**

* ...

## Proposal

Model-based directly on **raw data**
<figure><img src="image_227.png" /></figure>

<figure><img src="image_228.png" /></figure>

<figure><img src="image_229.png" /></figure>

<figure><img src="image_230.png" /></figure>

<figure><img src="image_231.png" /></figure>

<figure><img src="image_232.png" /></figure>

## Mixed data: conditional independence everywhere

The aim is to combine continuous, categorical and integer data

$$x_1 = (x_1^{\text{cont}}, x_1^{\text{cat}}, x_1^{\text{int}})$$

The proposed solution is to mixed all types by inter-type conditional independence

$$p(x_1; \alpha_k) = p(x_1^{\text{cont}}; \alpha_k^{\text{cont}}) \times p(x_1^{\text{cat}}; \alpha_k^{\text{cat}}) \times p(x_1^{\text{int}}; \alpha_k^{\text{int}})$$

In addition, for symmetry between types, intra-type conditional independence

Only need to define the univariate pdf for each variable type!

* Continuous: Gaussian

* Categorical: multinomial

* Integer: Poisson
<figure><img src="image_233.png" /></figure>

<figure><img src="image_234.png" /></figure>

<figure><img src="image_235.png" /></figure>

<figure><img src="image_236.png" /></figure>

<figure><img src="image_237.png" /></figure>

<figure><img src="image_238.png" /></figure>

Missing data: a seminal paper

Biometrika (1976), 63, 3, pp. 581-92

Printed in Great Britain

581

Inference and missing data

BY DONALD B. RUBIN

*Educational Testing Service, Princeton, New Jersey*

SUMMARY

When making sampling distribution inferences about the parameter of the data, $\theta$, it is appropriate to ignore the process that causes missing data if the missing data are ‘missing at random’ and the observed data are ‘observed at random’, but these inferences are generally conditional on the observed pattern of missing data. When making direct-likelihood or Bayesian inferences about $\theta$, it is appropriate to ignore the process that causes missing data if the missing data are missing at random and the parameter of the missing data process is ‘distinct’ from $\theta$. These conditions are the weakest general conditions under which ignoring the process that causes missing data always leads to correct inferences.

*Some key words:* Bayesian inference; Incomplete data; Likelihood inference; Missing at random; Missing data; Missing values; Observed at random; Sampling distribution inference.
<figure><img src="image_239.png" /></figure>

<figure><img src="image_240.png" /></figure>

<figure><img src="image_241.png" /></figure>

<figure><img src="image_242.png" /></figure>

<figure><img src="image_243.png" /></figure>

<figure><img src="image_244.png" /></figure>

Missing data: current solutions

<table><thead><tr><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>Cluster</th></tr></thead><tbody><tr><td>1.23</td><td>?</td><td>3.42</td><td>?</td></tr><tr><td>?</td><td>?</td><td>4.10</td><td>?</td></tr><tr><td>4.53</td><td>1.50</td><td>5.35</td><td>?</td></tr><tr><td>?</td><td>5.67</td><td>?</td><td>?</td></tr></tbody></table>

Discarded solutions

■ Suppress units and/or variables with missing data ⇒ loss of information

■ Imputation of the missing data by the mean or more evolved methods ⇒ uncertainty of the prediction not taken into account

Retained solution

Use an integrated approach which allows to take into account all the available information to perform clustering
<figure><img src="image_245.png" /></figure>

<figure><img src="image_246.png" /></figure>

<figure><img src="image_247.png" /></figure>

<figure><img src="image_248.png" /></figure>

<figure><img src="image_249.png" /></figure>

<figure><img src="image_250.png" /></figure>

<figure><img src="image_251.png" /></figure>

<figure><img src="image_252.png" /></figure>

<figure><img src="image_253.png" /></figure>

<figure><img src="image_254.png" /></figure>

Missing data: notations and MNAR assumption

■ $O_i \subseteq \{1, \dots, d\}$ the set of the observed variables from sample $i$

■ $x_i^O$ the observed data from sample $i$

■ $M_i$ the set of the missing variables for sample $i$

■ $\mu_{ik}^O$ the sub-vector of $\mu_k$ associated to index $O_i$ (the same for $M_i$)

■ $\Sigma_{ik}^{OM}$ the sub-matrix of $\Sigma_k$ associated to row $O_i$ and columns $M_i$ (the same for any other combination)

Assumption on the missingness mechanism

Missing At Randon (MAR): the probability that a variable is missing does not depend on its own value given the observed variables.
<figure><img src="image_255.png" /></figure>

<figure><img src="image_256.png" /></figure>

<figure><img src="image_257.png" /></figure>

<figure><img src="image_258.png" /></figure>

<figure><img src="image_259.png" /></figure>

<figure><img src="image_260.png" /></figure>

Missing data: maximum likelihood estimator

Observed log-likelihood...

$$ \ell(\theta; \mathbf{x}^O) = \sum_{i=1}^{n} \log \left( \sum_{k=1}^{K} \pi_k p(\mathbf{x}_i^O; \alpha_k) \right) = \ln \left[ \sum_{k=1}^{K} \pi_k \underbrace{\int_{\mathbf{x}_i^M} p(\mathbf{x}_i^O, \mathbf{x}_i^M; \alpha_k) d\mathbf{x}_i^M}_{\text{MAR assumption}} \right] $$
<figure><img src="image_261.png" /></figure>

<figure><img src="image_262.png" /></figure>

<figure><img src="image_263.png" /></figure>

<figure><img src="image_264.png" /></figure>

<figure><img src="image_265.png" /></figure>

<figure><img src="image_266.png" /></figure>

EM and Gaussian case: E step

$\theta$ and $\theta^+$ the parameters for two successive steps (idem for missing data)

$$z_{ik}^{+} = P(Z_{ik} = 1 | \mathbf{x}_{i}^{O}; \theta) = \frac{\pi_{k} \phi(\mathbf{x}_{i}^{O}; \mathbf{\Sigma}_{k})}{\sum_{\ell=1}^{K} \pi_{\ell} \phi(\mathbf{x}_{i}^{O}; \mathbf{\Sigma}_{\ell})}$$

$$\mathbf{x}_{ik}^{M+} = E[\mathbf{x}_{i}^{M} | \mathbf{x}_{i}^{O}, Z_{ik} = 1; \theta] = \mu_{ik}^{M} + \mathbf{\Sigma}_{ik}^{MO} (\mathbf{\Sigma}_{ik}^{OO})^{-1} (\mathbf{x}_{i}^{O} - \mu_{ik}^{O}).$$

Interpretation

- $z_{ik}^{+}$: class posterior probability membership given the available information $\mathbf{x}_{i}^{O}$.

- $\mathbf{x}_{ik}^{M+}$: conditional imputation of the missing data given the cluster.
<figure><img src="image_267.png" /></figure>

<figure><img src="image_268.png" /></figure>

<figure><img src="image_269.png" /></figure>

<figure><img src="image_270.png" /></figure>

<figure><img src="image_271.png" /></figure>

<figure><img src="image_272.png" /></figure>

EM and Gaussian case: M step

$$ \pi_k^+ = \frac{1}{n_k^+} \sum_{i=1}^{n} z_{ik}^+, \quad \mu_k^+ = \frac{1}{n_k^+} \sum_{i=1}^{n} z_{ik}^+ \mathbf{x}_{ik}^+ $$

$$ \Sigma_k^+ = \frac{1}{n_k^+} \sum_{i=1}^{n} z_{ik}^+ [(\mathbf{x}_{ik}^+ - \mu_k^+)(\mathbf{x}_{ik}^+ - \mu_k^+)' + \Sigma_{ik}^+] $$

where $n_k^+ = \sum_{i=1}^n z_{ik}^+$, $\mathbf{x}_{ik}^+ = (\mathbf{x}_i^O)$, $\Sigma_{ik}^+ = (\begin{pmatrix} 0_i^O & 0_i^{OM} \\ 0_i^{OM} & \Sigma_{ik}^{M+} \end{pmatrix})$ with 0 the $d \times d$ null matrix, and $\Sigma_{ik}^{M+} = \Sigma_{ik}^{MO} (\Sigma_{ik}^{O})^{-1} \Sigma_{ik}^{OM}$.

Interpretation of $\Sigma_{ik}^{M+}$

Variance correction due to the under-estimation of variability caused by the imputation of missing data.
<figure><img src="image_273.png" /></figure>

<figure><img src="image_274.png" /></figure>

<figure><img src="image_275.png" /></figure>

<figure><img src="image_276.png" /></figure>

<figure><img src="image_277.png" /></figure>

<figure><img src="image_278.png" /></figure>

# Missing data: SEM algorithm

A SEM algorithm to estimate $\theta$ by maximizing the observed-data log-likelihood

* Initialisation: $\theta^{(0)}$

* Iteration nb $q$:

* E-step: compute conditional probabilities $p(x^M, z|D; \theta^{(q)})$

* S-step: draw $(x^{M(q)}, z^{(q)})$ from $p(x^M, z|x^0; \theta^{(q)})$

* M-step: maximize $\theta^{(q+1)} = \arg\max_{\theta} \ln p(x^O, x^{M(q)}, z^{(q)}; \theta)$

* Stopping rule: iteration number

## Properties: simpler than EM and interesting properties!

* Avoid possibly difficult E-step in an EM

* Classical M steps

* Avoids local maxima

* The mean of the sequence $(\theta^{(q)})$ approximates $\hat{\theta}$

* The variance of the sequence $(\theta^{(q)})$ gives confidence intervals
<figure><img src="image_279.png" /></figure>

<figure><img src="image_280.png" /></figure>

<figure><img src="image_281.png" /></figure>

<figure><img src="image_282.png" /></figure>

<figure><img src="image_283.png" /></figure>

<figure><img src="image_284.png" /></figure>

Missing data: SE algorithm

A SE algorithm estimates then ($x^M, z^M$)

■ Iteration nb $q$:

■ E-step: compute conditional probabilities $p(x^M, z^M|D; \hat{\theta})$

■ S-step: draw ($x^{M(q)}, z^{M(q)}$) from $p(x^M, z^M|D; \hat{\theta})$

■ Stopping rule: iteration number

Properties

■ simplicity because of conditional independence

■ the mean/mode of the sequence ($x^{M(q)}, z^{M(q)}$) estimates ($x^M, z^M$)

■ confidence intervals are also derived
<figure><img src="image_285.png" /></figure>

<figure><img src="image_286.png" /></figure>

<figure><img src="image_287.png" /></figure>

<figure><img src="image_288.png" /></figure>

<figure><img src="image_289.png" /></figure>

<figure><img src="image_290.png" /></figure>

Missing data: illustration with the cancer data set (1/2)

■ Strategy "mice<fnref n="4" /> + mixture": mixture on the dataset completed by mice

> data.imp=mice(data)

> data.comp.mice=complete(data.imp)

■ Strategy "full mixture": mixture on the observed (no completed) dataset

<footnote n="4">⁴http://cran.r-project.org/web/packages/mice/mice.pdf</footnote>
<figure><img src="image_291.png" /></figure>

<figure><img src="image_292.png" /></figure>

<figure><img src="image_293.png" /></figure>

<figure><img src="image_294.png" /></figure>

<figure><img src="image_295.png" /></figure>

<figure><img src="image_296.png" /></figure>

MMissing data: illustration with the cancer data set (2/2)

<table><thead><tr><th>Strategy</th><th>mice + mixture</th><th>full mixture</th></tr></thead><tbody><tr><td>% misclassified</td><td>12.8</td><td>8.1</td></tr></tbody></table>

Avoid to complete missing data (imputation depends on the purpose)
<figure><img src="image_297.png" /></figure>

<figure><img src="image_298.png" /></figure>

<figure><img src="image_299.png" /></figure>

<figure><img src="image_300.png" /></figure>

<figure><img src="image_301.png" /></figure>

<figure><img src="image_302.png" /></figure>

## Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
<figure><img src="image_303.png" /></figure>

<figure><img src="image_304.png" /></figure>

<figure><img src="image_305.png" /></figure>

<figure><img src="image_306.png" /></figure>

<figure><img src="image_307.png" /></figure>

<figure><img src="image_308.png" /></figure>

Keep in mind

George E.P. Box (1987)

"Essentially, all models are wrong, but some are useful"

- $\mathcal{M} = \{\mathbf{m}\}$ will denote the set of competing models

- The true distribution $\mathbf{p}$ is **not** necessarily in $\mathcal{M}$

- **Density estimation**: AIC, BIC

- **Clustering**: ICL, CL, NEC
<figure><img src="image_309.png" /></figure>

<figure><img src="image_310.png" /></figure>

<figure><img src="image_311.png" /></figure>

<figure><img src="image_312.png" /></figure>

<figure><img src="image_313.png" /></figure>

<figure><img src="image_314.png" /></figure>

Importance of model selection: example

Model = number of clusters + parametric structure of clusters

Too simple model: bias

<figure><img src="image_315.png" /></figure>

true modèle: $[\pi \lambda_k l]$

too simple model: $[\pi \lambda l]$

Too complex model: variance

<figure><img src="image_316.png" /></figure>
<figure><img src="image_317.png" /></figure>

<figure><img src="image_318.png" /></figure>

<figure><img src="image_319.png" /></figure>

<figure><img src="image_320.png" /></figure>

<figure><img src="image_321.png" /></figure>

<figure><img src="image_322.png" /></figure>

<figure><img src="image_323.png" /></figure>

<figure><img src="image_324.png" /></figure>

<figure><img src="image_325.png" /></figure>

Importance of model selection: bias/variance trade-off

Partition error rate: $err(z_1, z_2) \ge 0$ a distance-like between two partitions $z_1, z_2$

Gap between true and model partition:

$$ \theta_m^* = \arg \min_{\theta \in \Theta_m} err(z, z(\theta)) $$

MLE:

$$ \hat{\theta}_m = \arg \max_{\theta \in \Theta} \ell(\theta; D) $$

Fundamental decomposition of $err(z, z(\hat{\theta}_m))$:

$$ \begin{aligned} & \mathrm{err}(z, z(\hat{\theta}_m)) \\ &= \{\mathrm{err}(z, z(\theta_m^*)) - \mathrm{err}(z, z)\} + \{\mathrm{err}(z, z(\hat{\theta}_m)) - \mathrm{err}(z, z(\theta_m^*))\} \\ &= \{\mathrm{bias}_m\} + \{\mathrm{variance}_m\} \end{aligned} $$
<figure><img src="image_326.png" /></figure>

<figure><img src="image_327.png" /></figure>

<figure><img src="image_328.png" /></figure>

<figure><img src="image_329.png" /></figure>

<figure><img src="image_330.png" /></figure>

<figure><img src="image_331.png" /></figure>

<figure><img src="image_332.png" /></figure>

<figure><img src="image_333.png" /></figure>

<figure><img src="image_334.png" /></figure>

<figure><img src="image_335.png" /></figure>

<figure><img src="image_336.png" /></figure>

<figure><img src="image_337.png" /></figure>

<figure><img src="image_338.png" /></figure>

<figure><img src="image_339.png" /></figure>

Importance of model selection: illustration of the variance effect

30 samples from a bivariate mixture with two components

$$ \pi_1 = \pi_2 = 0.5, \quad \mu_1 = (0, 0)', \quad \mu_2 = (2, 2)', \quad \Sigma_1 = \Sigma_2 = I $$

$\mathcal{M} = \{\text{spherical, general}\}$

<table><thead><tr><th>n</th><th>m</th><th>err(z, ŵ<sub>m</sub>)</th></tr></thead><tbody><tr><td rowspan="2">40</td><td>spherical</td><td>0.0967</td></tr><tr><td>general</td><td>0.1100</td></tr><tr><td rowspan="2">200</td><td>spherical</td><td>0.0840</td></tr><tr><td>general</td><td>0.0872</td></tr></tbody></table>
<figure><img src="image_340.png" /></figure>

<figure><img src="image_341.png" /></figure>

<figure><img src="image_342.png" /></figure>

<figure><img src="image_343.png" /></figure>

<figure><img src="image_344.png" /></figure>

<figure><img src="image_345.png" /></figure>

<figure><img src="image_346.png" /></figure>

<figure><img src="image_347.png" /></figure>

<figure><img src="image_348.png" /></figure>

<figure><img src="image_349.png" /></figure>

<figure><img src="image_350.png" /></figure>

<figure><img src="image_351.png" /></figure>

## Some heuristics entropy-based criteria: examples

* A fundamental decomposition of $\ell(\theta; \mathbf{x})$: for any "fuzzy partition" $\mathbf{c} = \{c_{ik}\}$

$$ 
\begin{aligned}
\ell(\theta; \mathbf{x}) &= \sum_{i=1}^{n} \sum_{k=1}^{K} c_{ik} \ln \{\pi_k p(\mathbf{x}_i; \alpha_k)\} - \sum_{i=1}^{n} \sum_{k=1}^{K} c_{ik} \ln t_{ik}(\theta) \\
&= \ell_c(\theta; \mathbf{x}, \mathbf{c}) + \xi(\theta; \mathbf{c}) \\
&= \text{complete-data log-likelihood} + \text{entropy}
\end{aligned}
$$

* NEC criterion (Normalized Entropy Criterion): retain $\mathbf{m}$ minimizing

$$ \text{NEC}_{\mathbf{K}} = \begin{cases} \frac{\xi(\hat{\theta}_{\mathbf{K}}; t(\hat{\theta}_{\mathbf{K}}))}{\ell(\hat{\theta}_{\mathbf{K}}; \mathbf{x}) - \ell(\hat{\theta}_{1}; \mathbf{x})} & \text{if } \mathbf{K} > 1 \\ 1 & \text{if } \mathbf{K} = 1 \end{cases} $$

* CL criterion (Completed Likelihood): retain $\mathbf{m}$ maximizing

$$ \text{CL} = \ell_c(\hat{\theta}; \mathbf{x}, \hat{\mathbf{z}}) = \underbrace{\ell(\hat{\theta}; \mathbf{x})}_{\text{model adequacy}} - \underbrace{\xi(\hat{\theta}; \hat{\mathbf{z}})}_{\text{partition evidence}} $$

* Behaviour: not completely satisfactory but something happens...
<figure><img src="image_352.png" /></figure>

<figure><img src="image_353.png" /></figure>

<figure><img src="image_354.png" /></figure>

<figure><img src="image_355.png" /></figure>

<figure><img src="image_356.png" /></figure>

<figure><img src="image_357.png" /></figure>

Some heuristics entropy-based criteria: NEC illustration

<figure><img src="image_358.png" /></figure>

<figure><img src="image_359.png" /></figure>
<figure><img src="image_360.png" /></figure>

<figure><img src="image_361.png" /></figure>

<figure><img src="image_362.png" /></figure>

<figure><img src="image_363.png" /></figure>

<figure><img src="image_364.png" /></figure>

<figure><img src="image_365.png" /></figure>

## Some heuristics entropy-based criteria: CL illustration

Interpretation as a penalized within-cluster inertia criterion:

$$CL_{([p \lambda l], K)} = -\frac{nd}{2} \ln(W_K) - n \ln(K) + \text{cst}$$

<figure><img src="image_366.png" /></figure>

<figure><img src="image_367.png" /></figure>

<figure><img src="image_368.png" /></figure>
<figure><img src="image_369.png" /></figure>

<figure><img src="image_370.png" /></figure>

<figure><img src="image_371.png" /></figure>

<figure><img src="image_372.png" /></figure>

<figure><img src="image_373.png" /></figure>

<figure><img src="image_374.png" /></figure>

## Theoretical model selection criteria

The most widespread principle

$$ \underbrace{\text{Criterion}}_{\text{to be maximized}} = \underbrace{\text{maximum log-likelihood}}_{\text{model-data adequacy}} - \boxed{\text{penalty}} $$
“cost” of the model

<table><thead><tr><th>criterion</th><th>penalty</th><th>interpretation</th><th>user purpose</th></tr></thead></table>

general criteria in statistics

<table><tr><td>AIC</td><td>ν</td><td>model complexity</td><td>prediction</td></tr><tr><td>BIC</td><td>0.5ν ln(n)</td><td>model complexity</td><td>identification</td></tr></table>

specific criterion for the clustering aim

<table><tr><td>ICL</td><td>0.5ν ln(n) - Σ<sub>j,k</sub> ẑ<sub>jk</sub> ln t<sub>jk</sub>(θ̂)</td><td>model complexity + partition entropy</td><td>well-separated clusters</td></tr></table>
<figure><img src="image_375.png" /></figure>

<figure><img src="image_376.png" /></figure>

<figure><img src="image_377.png" /></figure>

<figure><img src="image_378.png" /></figure>

<figure><img src="image_379.png" /></figure>

<figure><img src="image_380.png" /></figure>

<figure><img src="image_381.png" /></figure>

<figure><img src="image_382.png" /></figure>

BIC criterion: integrated likelihood

■ Posterior likelihood of m:

$$p(\mathbf{m}|\mathcal{D}) \propto p(\mathcal{D}|\mathbf{m}) \underbrace{p(\mathbf{m})}_{\text{prior on } \mathbf{m}}$$

■ Ideal model in a Bayesian context:

$$\hat{\mathbf{m}}^* \in \arg \max_{\mathbf{m} \in \mathcal{M}} p(\mathbf{m}|\mathcal{D})$$

■ Integrated likelihood: if $p(\mathbf{m}) = \text{cst}$, it is equivalent to maximize

$$p(\mathcal{D}|\mathbf{m}) = \int_{\Theta} p(\mathcal{D}; \theta, \mathbf{m}) \underbrace{p(\theta|\mathbf{m})}_{\text{prior on } \theta} d\theta$$

■ Difficulties:

■ Choose the prior $p(\theta|\mathbf{m})$

■ Evaluate the integral
<figure><img src="image_383.png" /></figure>

<figure><img src="image_384.png" /></figure>

<figure><img src="image_385.png" /></figure>

<figure><img src="image_386.png" /></figure>

<figure><img src="image_387.png" /></figure>

<figure><img src="image_388.png" /></figure>

BIC criterion: genesis

■ Laplace-Metropolis approximation: under standard regularity conditions, we have

$$ \ln p(D|\mathbf{m}) = \ell(\hat{\theta}; D) - \frac{\nu}{2} \ln(n) + O_p(1) $$

■ BIC criterion (Bayesian Information Criterion): retain **m** maximizing

$$ \text{BIC}_{\mathbf{m}} = \ell(\hat{\theta}_{\mathbf{m}}; D) - \frac{\nu_{\mathbf{m}}}{2} \ln(n) $$
<figure><img src="image_389.png" /></figure>

<figure><img src="image_390.png" /></figure>

<figure><img src="image_391.png" /></figure>

<figure><img src="image_392.png" /></figure>

<figure><img src="image_393.png" /></figure>

<figure><img src="image_394.png" /></figure>

<figure><img src="image_395.png" /></figure>

<figure><img src="image_396.png" /></figure>

BIC criterion: consistency

■ **Consistency:** BIC asymptotically selects the best

$$ \mathbf{m}^* = \arg \inf_{\mathbf{m} \in \mathcal{M}} \mathrm{KL}(\rho, \rho_{\theta^*}) $$

■ **Theoretical illustration of consistency:** $\mathbf{m}_1 \subseteq \mathbf{m}_2$, $\mathbf{m}_1$ being the true model, $\Delta\nu = \nu_2 - \nu_1$, $\Delta\ell = \ell(\hat{\theta}_2; \mathcal{D}) - \ell(\hat{\theta}_1; \mathcal{D})$, we have

$$ 2(\mathrm{BIC}_2 - \mathrm{BIC}_1) + \Delta\nu \ln(n) = 2\Delta\ell \xrightarrow{d} \chi_{\Delta\nu}^2 $$

With $\mu = \Delta\nu$ and $\sigma^2 = 2\Delta\nu$ the mean and the variance of $\chi_{\Delta\nu}^2$

$$ \mathrm{P}(\chi_{\Delta\nu}^2 > \Delta\nu \ln(n)) \le \mathrm{P}(|\chi_{\Delta\nu}^2 - \mu| > \Delta\nu \ln(n) - \mu) \le \frac{\sigma^2}{(\Delta\nu \ln(n) - \mu)^2} \xrightarrow{n \to \infty} 0 $$

by using the Chebyshev inequality. Thus, asymptotically, BIC will select $\mathbf{m}_1$

■ **Special case of K:** be careful on the $\chi^2$ approximation validity...
<figure><img src="image_397.png" /></figure>

<figure><img src="image_398.png" /></figure>

<figure><img src="image_399.png" /></figure>

<figure><img src="image_400.png" /></figure>

<figure><img src="image_401.png" /></figure>

<figure><img src="image_402.png" /></figure>

<figure><img src="image_403.png" /></figure>

<figure><img src="image_404.png" /></figure>

<figure><img src="image_405.png" /></figure>

<figure><img src="image_406.png" /></figure>

<figure><img src="image_407.png" /></figure>

<figure><img src="image_408.png" /></figure>

<figure><img src="image_409.png" /></figure>

<figure><img src="image_410.png" /></figure>

## ICL criterion: genesis

* Revisiting the fundamental decomposition: if z known, retain m maximizing

$$ \underbrace{\ln p(\mathbf{x}, \mathbf{z}|\mathbf{m})}_{\text{all data evidence}} = \underbrace{\ln p(\mathbf{x}|\mathbf{m})}_{\text{data x evidence}} + \underbrace{\ln p(\mathbf{z}|\mathbf{x}, \mathbf{m})}_{\text{partition z evidence}} $$

Thus models leading to overlapping groups are more penalized (low z evidence)

* ICL criterion (Integrated Classification Likelihood): replace z by $\hat{z}$

$$ \text{ICL} = \ln p(\mathbf{x}, \hat{z}|\mathbf{m}) $$

* BIC-like approximation of ICL:

$$ \ln p(\mathbf{x}, \mathbf{z}|\mathbf{m}) = \ln p(\mathbf{x}, \mathbf{z}; \hat{\theta}_{\mathbf{x},\mathbf{z}}|\mathbf{m}) - \frac{\nu}{2} \ln n + O_p(1) $$

In case of the right model m: $\hat{\theta}_{\mathbf{x},\mathbf{z}} \xrightarrow{\text{a.s.}} \theta^*$ and $\hat{\theta}_{\mathbf{x}} \xrightarrow{\text{a.s.}} \theta^*$. Thus, for n large enough, $\hat{\theta}_{\mathbf{x},\mathbf{z}} \approx \theta_{\mathbf{x}}$. Then, we take $\hat{z} = \text{MAP}(\hat{\theta}_{\mathbf{x}})$ (or also $\hat{z} = \mathbf{t}(\hat{\theta}_{\mathbf{x}})$). It gives

$$ \begin{align*} \text{ICLbic} &= \ln p(\mathbf{x}, \hat{z}; \hat{\theta}_{\mathbf{x}}) - \frac{\nu}{2} \ln n \\ &= \text{BIC} - \xi(\hat{\theta}_{\mathbf{x}}; \hat{z}) \\ &= \text{CL} - \frac{\nu}{2} \ln n \end{align*} $$
<figure><img src="image_411.png" /></figure>

<figure><img src="image_412.png" /></figure>

<figure><img src="image_413.png" /></figure>

<figure><img src="image_414.png" /></figure>

<figure><img src="image_415.png" /></figure>

<figure><img src="image_416.png" /></figure>

<figure><img src="image_417.png" /></figure>

<figure><img src="image_418.png" /></figure>

<figure><img src="image_419.png" /></figure>

<figure><img src="image_420.png" /></figure>

<figure><img src="image_421.png" /></figure>

<figure><img src="image_422.png" /></figure>

<figure><img src="image_423.png" /></figure>

<figure><img src="image_424.png" /></figure>

## ICL criterion: robustness to model misspecification

* A bivariate mixture of a uniform and a Gaussian cluster:

* non-Gaussian component: $\pi_1 = 0.5$, $p_1(x_1) = 0.25 I_{[-1,1]}(x^1) I_{[-1,1]}(x^2)$

* Gaussian component: $\pi_2 = 0.5$, $\mu_2 = (3.3, 0)'$, $\Sigma_2 = I$

* 50 simulated data sets of size $n = 200$

<figure><img src="image_425.png" /></figure>

<table><thead><tr><th>K</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>BIC</td><td>.</td><td><b>60</b></td><td>.</td><td>32</td><td>8</td></tr><tr><td>ICLbic</td><td>.</td><td><b>100</b></td><td>.</td><td>.</td><td>.</td></tr></tbody></table>
<figure><img src="image_426.png" /></figure>

<figure><img src="image_427.png" /></figure>

<figure><img src="image_428.png" /></figure>

<figure><img src="image_429.png" /></figure>

<figure><img src="image_430.png" /></figure>

<figure><img src="image_431.png" /></figure>

## ICL criterion: consistency?

* Assumption: true model with two groups and parameter $\theta_2^*$

* Theoretical result:

* Preliminaries: $\delta_n = n(\theta_2^* - \theta_2^{*p})'J(\theta_2^*)(\theta_2^* - \theta_2^{*p})$, $J(\theta_2^*)$ the Fisher matrix for a data unit calculated with the true parameter $\theta_2$ and $\theta_2^{*p}$ its projected value on the parameter subspace associated to the one component case, $\mu_n = [\text{E}[\chi_{\Delta\nu}^2(\delta_n)] = \Delta\nu + \delta_n$, $\sigma_n^2 = \text{Var}[\chi_{\Delta\nu}^2(\delta_n)] = 2(\Delta\nu + \delta_n)$

* Asymptotically: by Chebishev inequality, with $\mu_n - \Delta\nu \ln n - 2n \ln 2 > 0$

$$p(\text{choose wrong model}) = p(\text{ICLBic}_2 < \text{ICLBic}_1) \le \frac{\sigma_n^2}{(\mu_n - \Delta\nu \ln n - 2n \ln 2)^2}$$

Thus it goes towards 0 for well-separated groups

* Experimental result: 100 samples from a univariate Gaussian mixture

$$\pi_1 = \pi_2, \quad \mu_1 = 0, \quad \mu_2 = \Delta\mu, \quad \sigma_1^2 = \sigma_2^2 = 1$$

<table><thead><tr><th>Δμ</th><th>2.9</th><th>3.0</th><th>3.1</th><th>3.2</th><th>3.3</th></tr><tr><th>n</th><th>BIC</th><th>ICL</th><th>BIC</th><th>ICL</th><th>BIC</th><th>ICL</th><th>BIC</th><th>ICL</th></tr></thead><tbody><tr><td>100</td><td>94</td><td>23</td><td>96</td><td>31</td><td>97</td><td>44</td><td>95</td><td>45</td><td>97</td><td>60</td></tr><tr><td>400</td><td>100</td><td>9</td><td>100</td><td>21</td><td>100</td><td>48</td><td>100</td><td>70</td><td>100</td><td>85</td></tr><tr><td>700</td><td>100</td><td>8</td><td>100</td><td>15</td><td>100</td><td>39</td><td>100</td><td>72</td><td>100</td><td>96</td></tr><tr><td>1 000</td><td>100</td><td>6</td><td>100</td><td>16</td><td>100</td><td>56</td><td>100</td><td>75</td><td>100</td><td>91</td></tr></tbody></table>
<figure><img src="image_432.png" /></figure>

<figure><img src="image_433.png" /></figure>

<figure><img src="image_434.png" /></figure>

<figure><img src="image_435.png" /></figure>

<figure><img src="image_436.png" /></figure>

<figure><img src="image_437.png" /></figure>

Large *n*: BIC behaviour (1/2)

* The mixture density is wrong (as all models)

* Mixtures allow to estimate any distribution by increasing the number of components (high flexibility)

<figure><img src="image_438.png" /></figure>

<figure><img src="image_439.png" /></figure>

<figure><img src="image_440.png" /></figure>
<figure><img src="image_441.png" /></figure>

<figure><img src="image_442.png" /></figure>

<figure><img src="image_443.png" /></figure>

<figure><img src="image_444.png" /></figure>

<figure><img src="image_445.png" /></figure>

<figure><img src="image_446.png" /></figure>

## Large *n*: BIC behaviour (2/2)

Since BIC is consistent, as *n* grows, it adds components for improving the true density estimation

### Real example

<figure><img src="image_447.png" /></figure>

Marketing context :

Looking for classes of customers

<figure><img src="image_448.png" /></figure>

It was only d=2
in this example

**Reality is even worst:** n=10⁵ customers, d=77, mixed,
1 day computer for 20 classes, more than 40 classes!
<figure><img src="image_449.png" /></figure>

<figure><img src="image_450.png" /></figure>

<figure><img src="image_451.png" /></figure>

<figure><img src="image_452.png" /></figure>

<figure><img src="image_453.png" /></figure>

<figure><img src="image_454.png" /></figure>

Missing data: illustration with the cancer data set (1/2)

■ Strategy "mice<fnref n="5" /> + mixture": mixture on the dataset completed by mice

> data.imp=mice(data)

> data.comp.mice=complete(data.imp)

■ Strategy "full mixture": mixture on the observed (no completed) dataset

<footnote n="5">⁵http://cran.r-project.org/web/packages/mice/mice.pdf</footnote>
<figure><img src="image_455.png" /></figure>

<figure><img src="image_456.png" /></figure>

<figure><img src="image_457.png" /></figure>

<figure><img src="image_458.png" /></figure>

<figure><img src="image_459.png" /></figure>

<figure><img src="image_460.png" /></figure>

Missing data: illustration with the cancer data set (2/2)

<figure><img src="image_461.png" /></figure>

mice + mixture
$\hat{K} = 7$

<figure><img src="image_462.png" /></figure>

full mixture
$\hat{K} = 2$

... may lose some cluster information when imputation before clustering
<figure><img src="image_463.png" /></figure>

<figure><img src="image_464.png" /></figure>

<figure><img src="image_465.png" /></figure>

<figure><img src="image_466.png" /></figure>

<figure><img src="image_467.png" /></figure>

<figure><img src="image_468.png" /></figure>

CAUTION

Impossible to use BIC/ICL for kernel/spectral clustering (data set has changed...)
<figure><img src="image_469.png" /></figure>

<figure><img src="image_470.png" /></figure>

<figure><img src="image_471.png" /></figure>

<figure><img src="image_472.png" /></figure>

<figure><img src="image_473.png" /></figure>

<figure><img src="image_474.png" /></figure>

<figure><img src="image_475.png" /></figure>

<figure><img src="image_476.png" /></figure>

Reformulate K-means: elbow as a slope heuristics (1/3)

■ SH (Slope Heuristics) criterion: retain **m** maximizing

$$SH_m = \ell(\hat{\theta}_m; D) - 2\text{variance}_m$$

■ Estimating the penalty: optimal penalty is linear in $\nu_m$

$$2\text{variance}_m = \kappa\nu_m.$$

and also

$$2\widehat{\text{variance}}_m = 2\left\{ \ell(\hat{\theta}_m; D) - p(D) \right\} + 2\left\{ p(D) - \ell(\theta_m; D) \right\} \\ \approx \kappa\nu_m \\ \text{bias} \approx \text{cust for too complex models}$$

thus, for complex enough models, $\ell(\hat{\theta}_m; D)$ behaves linearly with $\nu_m$ and the corresponding slope is $\kappa/2$

■ CAPUSHE<fnref n="6" /> (CALibrated Penalty Using Slope HEuristics): $\kappa/2$ can be estimated by a linear regression of $\ell(\hat{\theta}_m; D)$ on $\frac{\kappa}{2}\nu_m$

<footnote n="6">⁶http://cran.r-project.org/web/packages/capushe/</footnote>
<figure><img src="image_477.png" /></figure>

<figure><img src="image_478.png" /></figure>

<figure><img src="image_479.png" /></figure>

<figure><img src="image_480.png" /></figure>

<figure><img src="image_481.png" /></figure>

<figure><img src="image_482.png" /></figure>

Reformulate K-means: elbow as a slope heuristics (2/3)

<figure><img src="image_483.png" /></figure>

<figure><img src="image_484.png" /></figure>
<figure><img src="image_485.png" /></figure>

<figure><img src="image_486.png" /></figure>

<figure><img src="image_487.png" /></figure>

<figure><img src="image_488.png" /></figure>

<figure><img src="image_489.png" /></figure>

<figure><img src="image_490.png" /></figure>

Reformulate K-means: elbow as a slope heuristics (2/3)

<figure><img src="image_491.png" /></figure>

<figure><img src="image_492.png" /></figure>
<figure><img src="image_493.png" /></figure>

<figure><img src="image_494.png" /></figure>

<figure><img src="image_495.png" /></figure>

<figure><img src="image_496.png" /></figure>

<figure><img src="image_497.png" /></figure>

<figure><img src="image_498.png" /></figure>

## Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
<figure><img src="image_499.png" /></figure>

<figure><img src="image_500.png" /></figure>

<figure><img src="image_501.png" /></figure>

<figure><img src="image_502.png" /></figure>

<figure><img src="image_503.png" /></figure>

<figure><img src="image_504.png" /></figure>

Outliers: Two possibilities

■ "After": exclude data outside the confidence area of clusters

<figure><img src="image_505.png" /></figure>

■ "During": model outliers as a particular cluster in the mixture

$$p(\mathbf{x}_1; \theta) = \sum_{k=1}^{K} \pi_k p(\mathbf{x}_1; \alpha_k) + \pi_0 p(\mathbf{x}_1; \alpha_0)$$
<figure><img src="image_506.png" /></figure>

<figure><img src="image_507.png" /></figure>

<figure><img src="image_508.png" /></figure>

<figure><img src="image_509.png" /></figure>

<figure><img src="image_510.png" /></figure>

<figure><img src="image_511.png" /></figure>

Outliers: “during” example with acoustic emission control

■ Data: $n = 2$ 061 event locations in a rectangle of $R^2$ representing the vessel

■ Model: Diagonal Gaussian mixture + uniform (noise)

■ Groups: sound locations = vessel defects

<figure><img src="image_512.png" /></figure>
<figure><img src="image_513.png" /></figure>

<figure><img src="image_514.png" /></figure>

<figure><img src="image_515.png" /></figure>

<figure><img src="image_516.png" /></figure>

<figure><img src="image_517.png" /></figure>

<figure><img src="image_518.png" /></figure>

<figure><img src="image_519.png" /></figure>

<figure><img src="image_520.png" /></figure>

<figure><img src="image_521.png" /></figure>

<figure><img src="image_522.png" /></figure>

# Units: changing the data units

## Principle of data units transformation u:

$$ 
\begin{array}{lclcl}
u: & \mathcal{X} = \mathcal{X}^{\text{id}} & \longrightarrow & \mathcal{X}^{\text{u}} \\
& \mathbf{x} = \mathbf{x}^{\text{id}} = \text{id}(\mathbf{x}) & \longrightarrow & \mathbf{x}^{\text{u}} = \mathbf{u}(\mathbf{x})
\end{array}
$$

* u is a bijective mapping to preserve the whole data set information quantity

* We denote by $u^{-1}$ the reciprocal of u, so $u^{-1} \circ u = \text{id}$

* Thus, id is only a particular unit u

* Often a meaningful restriction<fnref n="7" /> on u: it proceeds lines by lines and rows by rows

$$ \mathbf{u}(\mathbf{x}) = (\mathbf{u}(\mathbf{x}_1), \dots, \mathbf{u}(\mathbf{x}_n)) \quad \text{with} \quad \mathbf{u}(\mathbf{x}_i) = (\mathbf{u}_1(x_{i1}), \dots, \mathbf{u}_d(x_{id})) $$

* Advantage to respect the variable definition, transforming only its unit

* $\mathbf{u}(\mathbf{x}_i)$ means that u applied to the data set $\mathbf{x}_i$, restricted to the single individual i

* $\mathbf{u}_i$ corresponds to the specific (bijective) transformation unit associated to variable j

<footnote n="7">⁷Possibility to relax this restriction, including for instance linear transformations involved in PCA (principal component analysis). But the variable definition is no longer respected.</footnote>
<figure><img src="image_523.png" /></figure>

<figure><img src="image_524.png" /></figure>

<figure><img src="image_525.png" /></figure>

<figure><img src="image_526.png" /></figure>

<figure><img src="image_527.png" /></figure>

<figure><img src="image_528.png" /></figure>

# Units: revisiting units as a modelling component

* Explicitly exhibiting the “canonical” unit id in the model

$$p_m = \{ \cdot \in \mathcal{X} \mapsto p(\cdot; \theta) : \theta \in \Theta_m \} = \{ \cdot \in \mathcal{X}^{\text{id}} \mapsto p(\cdot; \theta) : \theta \in \Theta_m \} = p_m^{\text{id}}$$

* Thus the variable space and the probability measure are embedded

* As the standard probability theory: a couple (variable space, probability measure)!

* Changing id into u, while preserving m, is expected to produce a new modelling

$$p_m^u = \{ \cdot \in \mathcal{X}^u \mapsto p(\cdot; \theta) : \theta \in \Theta_m \}.$$

A model should be systematically defined by a couple (u,m), denoted by $p_m^u$
<figure><img src="image_529.png" /></figure>

<figure><img src="image_530.png" /></figure>

<figure><img src="image_531.png" /></figure>

<figure><img src="image_532.png" /></figure>

<figure><img src="image_533.png" /></figure>

<figure><img src="image_534.png" /></figure>

Units: interpretation and identifiability of $p_m^u$

* Standard probability theory (again): there exists a measure $u^{-1}(m)$ s.t.<fnref n="8" />

$$u^{-1}(m) \in \{m' \in \mathbb{M} : p_{m'}^{id} = p_m^u\}$$

* There exists two alternative interpretations of strictly the same model:

* $p_m^u$: data measured with unit $u$ arise from measure $m$;

* $p_{u^{-1}(m)}^u$: data measured with unit $id$ arise from measure $u^{-1}(m)$

* Two points of view:

**Statistician**

The model $p_m^u$ is not identifiable over the couple $(m, u)$

**Practitioner**

Freedom to choose the interpretation which is the most meaningful for him

<footnote n="8">⁸ This set is usually restricted to a single element</footnote>
<figure><img src="image_535.png" /></figure>

<figure><img src="image_536.png" /></figure>

<figure><img src="image_537.png" /></figure>

<figure><img src="image_538.png" /></figure>

<figure><img src="image_539.png" /></figure>

<figure><img src="image_540.png" /></figure>

Units: opportunity for designing new models

Great opportunity to build easily numerous new **meaningful models** $p_m^u!$

■ Just combine a standard model family {$m$} with a standard unit family {$u$}

■ New family can be huge! **Combinatorial problems** can occur...

■ Some model stability can exist in some (specific) cases: $m = u^{-1}(m)$
<figure><img src="image_541.png" /></figure>

<figure><img src="image_542.png" /></figure>

<figure><img src="image_543.png" /></figure>

<figure><img src="image_544.png" /></figure>

<figure><img src="image_545.png" /></figure>

<figure><img src="image_546.png" /></figure>

Units: model selection

As any model, possible to choose between $p_{m1}^{u1}$ and $p_{m2}^{u2}$

However, caution when using likelihood-based model selection criteria (as BIC)

■ Prohibited to compare $m_1$ in unit $u_1$ and $m_2$ in unit $u_2$

■ But allowed after transforming in identical unit id

■ Thus compare their equivalent expression: $p_{u_1^{-1}(m_1)}^{id}$ and $p_{u_2^{-1}(m_2)}^{id}$

■ Example for abs. continuous x and differentiable u, the density transform in id is:

$$p_{u^{-1}(m)}^{id} = \{ \cdot \in \mathcal{X}^{id} \mapsto p(\mathbf{u}(\cdot); \boldsymbol{\theta}) \times |\mathbf{J}^{\mathbf{u}}(\cdot)| : \boldsymbol{\theta} \in \Theta_{\mathbf{m}} \}$$

with $\mathbf{J}^{\mathbf{u}}(\cdot)$ the Jacobian associated to the transformation $\mathbf{u}$
<figure><img src="image_547.png" /></figure>

<figure><img src="image_548.png" /></figure>

<figure><img src="image_549.png" /></figure>

<figure><img src="image_550.png" /></figure>

<figure><img src="image_551.png" /></figure>

<figure><img src="image_552.png" /></figure>

Units: prostate cancer data (1/2)

■ **Individuals:** 506 patients with prostatic cancer grouped on clinical criteria into two Stages 3 and 4 of the disease

■ **Variables:** $d = 12$ pre-trial variates were measured on each patient, composed by

■ Eight continuous variables (age, weight, systolic blood pressure, diastolic blood pressure, serum haemoglobin, size of primary tumour "SZ", index of tumour stage and histologic grade, serum prostatic acid phosphatase "AP")

■ Two ordinal variables (performance rating, cardiovascular disease history)

■ Two categorical variables with various numbers of levels (electrocardiogram code, bone metastases)

■ **Some missing data:** 62 missing values ($\approx$ 1%)

■ Two historical units for performing the clustering task:

■ Raw units id: [McParland & Gormley, 2015]<fnref n="9" />

■ Transformed data u: since SZ and AP are skewed, [Jorgensen & Hunt, 1996]<fnref n="10" /> propose

$$u_{SZ} = \sqrt{\cdot} \text{ and } u_{AP} = \ln(\cdot)$$

<footnote n="9">⁹McParland, D. and Gormley, I. C. (2015). Model based clustering for mixed data: clustmd. arXiv preprint arXiv:1511.01720.</footnote>

<footnote n="10">¹⁰Jorgensen, M. and Hunt, L. (1996). Mixture model clustering of data sets with categorical and continuous variables. In Proceedings of the Conference ISIS, volume 96, pages 375–384.</footnote>
<figure><img src="image_553.png" /></figure>

<figure><img src="image_554.png" /></figure>

<figure><img src="image_555.png" /></figure>

<figure><img src="image_556.png" /></figure>

<figure><img src="image_557.png" /></figure>

<figure><img src="image_558.png" /></figure>

Units: prostate cancer data (2/2)

■ Model **m**: full mixed data $\mathbf{x} = (\mathbf{x}^{\text{cont}}, \mathbf{x}^{\text{cat}}, \mathbf{x}^{\text{ordi}}, \mathbf{x}^{\text{int}}, \mathbf{x}^{\text{rank}})$ (missing data are allowed also) are simply modeled by **inter conditional independence**

$$p(\mathbf{x}; \alpha_k) = p(\mathbf{x}^{\text{cont}}; \alpha_k^{\text{cont}}) \times p(\mathbf{x}^{\text{cat}}; \alpha_k^{\text{cat}}) \times p(\mathbf{x}^{\text{ordi}}; \alpha_k^{\text{ordi}}) \times \dots$$

In addition, for symmetry between types, **intra conditional independence** for each

■ Results:

■ New units $u_{SZ}$ and $u_{AP}$ are selected by ICL

■ New units allow to select two groups and provides a lower error rate

<figure><img src="image_559.png" /></figure>

<table><thead><tr><th colspan="2">clusters</th></tr><tr><th>1</th><th>2</th></tr></thead><tbody><tr><td>287</td><td>5</td></tr><tr><td>52</td><td>162</td></tr></tbody></table>

<table><thead><tr><th colspan="2">clusters</th></tr><tr><th>1</th><th>2</th></tr></thead><tbody><tr><td>270</td><td>22</td></tr><tr><td>23</td><td>191</td></tr></tbody></table>

Table: Raw units: 11% misclassified

Table: New units: 9% misclassified
<figure><img src="image_560.png" /></figure>

<figure><img src="image_561.png" /></figure>

<figure><img src="image_562.png" /></figure>

<figure><img src="image_563.png" /></figure>

<figure><img src="image_564.png" /></figure>

<figure><img src="image_565.png" /></figure>

Variable selection for Gaussians<fnref n="11" />

Definition

$$p(\mathbf{x}_1; \theta) = \underbrace{\left\{ \sum_{k=1}^{K} \pi_k p(\mathbf{x}_1^S; \mu_k, \Sigma_k) \right\}}_{\text{clustering variables}} \times \underbrace{\left\{ p(\mathbf{x}_1^U; \mathbf{a} + \mathbf{x}_1^R \mathbf{b}, \mathbf{C}) \right\}}_{\text{redundant variables}} \times \underbrace{\left\{ p(\mathbf{x}_1^W; \mathbf{u}, \mathbf{V}) \right\}}_{\text{independent variables}}$$

where

- all parts are Gaussians

- S: set of variables useful for clustering

- U: set of redundant clustering variables, expressed with $R \subseteq S$

- W: set of variables independent of clustering

Trick

Variable selection is recasted as a particular model selected by BIC

<footnote n="11">¹¹Raftery and Dean (2006), Maugis et al. (09a), Maugis et al. (09b)</footnote>
<figure><img src="image_566.png" /></figure>

<figure><img src="image_567.png" /></figure>

<figure><img src="image_568.png" /></figure>

<figure><img src="image_569.png" /></figure>

<figure><img src="image_570.png" /></figure>

<figure><img src="image_571.png" /></figure>

# Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
<figure><img src="image_572.png" /></figure>

<figure><img src="image_573.png" /></figure>

<figure><img src="image_574.png" /></figure>

<figure><img src="image_575.png" /></figure>

<figure><img src="image_576.png" /></figure>

<figure><img src="image_577.png" /></figure>

# (R)MixtComp package

## MixtComp

<figure><img src="image_578.png" /></figure>

launch

binder

Visit notebooks repo

MixtComp takes mixture model analysis one step further and deals with mixed, missing or uncertain data which are common in today's data sets. Mixed data concerned by MixtComp include continuous, categorical, integer and functional (as time series) ones. All of them can be combined, offering many possibilities as multivariate time series, etc.

A related publication for this software is still ongoing but you can find early information here. Here is the related package on the CRAN.

## Two ways to be used

* On the CRAN:
  https://cran.r-project.org/web/packages/RMixtComp/index.html

* As a note book (nothing to install):
  https://team.inria.fr/modal/software/notebooks/
<figure><img src="image_579.png" /></figure>

<figure><img src="image_580.png" /></figure>

<figure><img src="image_581.png" /></figure>

<figure><img src="image_582.png" /></figure>

<figure><img src="image_583.png" /></figure>

<figure><img src="image_584.png" /></figure>

Example 1: prostate cancer data<fnref n="12" />

■ **Individuals**: 506 patients with prostatic cancer grouped on clinical criteria into two Stages 3 and 4 of the disease

■ **Variables**: *d* = 12 pre-trial variates were measured on each patient, composed by eight continuous variables (age, weight, systolic blood pressure, diastolic blood pressure, serum haemoglobin, size of primary tumour, index of tumour stage and histolic grade, serum prostatic acid phosphatase) and four categorical variables with various numbers of levels (performance rating, cardiovascular disease history, electrocardiogram code, bone metastases)

■ **Some missing data**: 62 missing values (≈ 1%)

We forget the classes (Stages of the disease) for performing clustering

Questions

■ How many clusters?

■ Which partition?

<footnote n="12">¹²Byar DP, Green SB (1980): Bulletin Cancer, Paris 67:477-488</footnote>
<figure><img src="image_585.png" /></figure>

<figure><img src="image_586.png" /></figure>

<figure><img src="image_587.png" /></figure>

<figure><img src="image_588.png" /></figure>

<figure><img src="image_589.png" /></figure>

<figure><img src="image_590.png" /></figure>

# RMixtComp basic example

## Data

Load the CSV data file as dataframe.

```r
[1]: data <- read.table("mixtcomp-example.csv", sep = ";", header = TRUE)
head(data)
```

A data.frame: 6 x 12

<table><thead><tr><th>Age</th><th>Wt</th><th>PF</th><th>HX</th><th>SBP</th><th>DBP</th><th>EKG</th><th>HG</th><th>SZ</th><th>SG</th><th>AP</th><th>BM</th></tr><tr><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;chr&gt;</th><th>&lt;dbl&gt;</th><th>&lt;chr&gt;</th></tr></thead><tbody><tr><td>75</td><td>76</td><td>1</td><td>1</td><td>15</td><td>9</td><td>5</td><td>138</td><td>1.4142</td><td>8</td><td>1.0986</td><td>1</td></tr><tr><td>76</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>5.3852</td><td>9</td><td>2.4849</td><td>?</td></tr><tr><td>54</td><td>116</td><td>1</td><td>1</td><td>13</td><td>7</td><td>4</td><td>146</td><td>6.4807</td><td>?</td><td>1.9459</td><td>1</td></tr><tr><td>69</td><td>102</td><td>1</td><td>2</td><td>14</td><td>8</td><td>5</td><td>134</td><td>1.7321</td><td>9</td><td>1.0986</td><td>1</td></tr><tr><td>66</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>?</td><td>1</td><td>9</td><td>2.3979</td><td>?</td></tr><tr><td>75</td><td>94</td><td>2</td><td>2</td><td>14</td><td>7</td><td>2</td><td>176</td><td>2</td><td>8</td><td>2.1972</td><td>1</td></tr></tbody></table>
<figure><img src="image_591.png" /></figure>

<figure><img src="image_592.png" /></figure>

<figure><img src="image_593.png" /></figure>

<figure><img src="image_594.png" /></figure>

<figure><img src="image_595.png" /></figure>

<figure><img src="image_596.png" /></figure>

## Clustering with RMixtComp

Launch the RMixtComp package.

[2] library(RMixtComp)

Loading required package: RMixtCompUtilities

Define the distribution used for each variable.

[3] model <- list(Age = "Gaussian", Mt = "Gaussian", PF = "Multinomial", HX = "Multinomial", SBP = "Gaussian", DBP = "Gaussian", FKG = "Multinomial", HG = "Gaussian", SZ = "Gaussian", SG = "Gaussian", AP = "Gaussian", BM = "Multinomial")

Define the SEM algorithm's parameters

[4] algo = List(nBurnInIter = 50, nIter = 100, nbGibbsBurnInIter = 50, nbGibbsIter = 100, nIntPerClass = floor(nrow(data)/2), nStry = 5, confidenceLevel = 0.95, ratioStableCriterion = 0.99, nStableCriterion = 10)

Choose the desired number of classes and the number of runs for each given number of classes.

[5] nClass <= 1:8
nRun <= 3

```r
[1] run <- RMixtCompUtilities, model = 'nClass = nClass, criterion = "ICL", shus = 0, nrun = 0

=e= Run MixtComp in learn mode with 3 run(s) per number of classes and 1 core(s)
Data: 506 Individuals and 12 variables.
  - K = 2

  Start run 1 on thread number 83
  Run 1 DONE on thread number 83
  Start run 2 on thread number 83
  Run 2 DONE on thread number 83
  Start run 3 on thread number 83
  Run 3 DONE on thread number 83
  Run completed in 8.222s
  - K = 2

  Start run 1 on thread number 83
  Run 1 DONE on thread number 83
  Start run 2 on thread number 83
  Run 2 DONE on thread number 83
  Start run 3 on thread number 83
  Run 3 DONE on thread number 83
  Run completed in 8.525s
  - K = 3

  Start run 1 on thread number 83
  Run 1 DONE on thread number 83
  Start run 2 on thread number 83
  Run 2 DONE on thread number 83
  Start run 3 on thread number 83
  Run 3 DONE on thread number 83
  Run completed in 8.655s
```

```r
***
  - K = 8
```

```r
Start run 1 on thread number 83
Run 1 DONE on thread number 83
Start run 2 on thread number 83
Run 2 DONE on thread number 83
Start run 3 on thread number 83
Run 3 DONE on thread number 83
Run completed in 1.161s
Total runtime: 6.864s
Best model according to ICL: 2 clusters.
```
<figure><img src="image_597.png" /></figure>

<figure><img src="image_598.png" /></figure>

<figure><img src="image_599.png" /></figure>

<figure><img src="image_600.png" /></figure>

<figure><img src="image_601.png" /></figure>

<figure><img src="image_602.png" /></figure>

# Output's Analysis

## Criterion

Draw the criterion value (BIC and ICL) for each model that was built. The higher the value (close to 0) the better the model.

[9]: plotCrit(res, pkg = "plotly")

<figure><img src="image_603.png" /></figure>

Choose the number of classes to study in the following.

[10]: K <- 3
resK <- extractMixtCompObject(res, K)
<figure><img src="image_604.png" /></figure>

<figure><img src="image_605.png" /></figure>

<figure><img src="image_606.png" /></figure>

<figure><img src="image_607.png" /></figure>

<figure><img src="image_608.png" /></figure>

<figure><img src="image_609.png" /></figure>

Variables

Draw the discriminating level of each variable. A high value (close to one) means that the variable is highly discriminating. A low value (close to zero) means that the variable is poorly discriminating.

[11] plotDiscrimVar(resk, pkg = "plotly")

<figure><img src="image_610.png" /></figure>

Draw the similarity between every pair of variables. A high value (close to one) means that the two variables provide the same information for the clustering task (i.e. similar partitions). A low value (close to zero) means that the two variables provide some different information for the clustering task (a. different partitions).

[12] heatmapplot(resk, pkg = "plotly")

<figure><img src="image_611.png" /></figure>

Select a variable to draw its distribution.

[13] variable <- "SG" plotbeta <- plotplot(resk, variable, gr = TRUE, pkg = "plotly")

<figure><img src="image_612.png" /></figure>
Need to formalize

<figure><img src="image_613.png" /></figure>

Formaling estimation

<figure><img src="image_614.png" /></figure>

<figure><img src="image_615.png" /></figure>

Formaling selection

<figure><img src="image_616.png" /></figure>

<figure><img src="image_617.png" /></figure>

More advanced formalizing

<figure><img src="image_618.png" /></figure>

Experiments

<figure><img src="image_619.png" /></figure>

To go further

<figure><img src="image_620.png" /></figure>

Classes

Draw the proportion of individuals in each class.

<figure><img src="image_621.png" /></figure>

Draw the similarity level between each pair of classes. A high value (close to one) means that the 2 classes are strongly different (i.e. low overlapping). A low value (close to zero) means that the 2 classes are similar for the clustering task (i.e. high overlapping).

<figure><img src="image_622.png" /></figure>

Draw the discriminating level of each variable for the selected class.

<figure><img src="image_623.png" /></figure>

Select a variable to draw its distribution for the selected class.

<figure><img src="image_624.png" /></figure>
Need to formalize

<figure><img src="image_625.png" /></figure>

Formalizing estimation

<figure><img src="image_626.png" /></figure>

Formalizing selection

<figure><img src="image_627.png" /></figure>

More advanced formalizing

<figure><img src="image_628.png" /></figure>

Experiments

<figure><img src="image_629.png" /></figure>

To go further

<figure><img src="image_630.png" /></figure>

Probabilities

Draw the probability of assignment to a class for each individual. Individuals have been reordered in decreasing assignment probability.

[18]: heatmapTikSorted(resK, pkg = "plotly")

Probabilities of classification

<figure><img src="image_631.png" /></figure>

<figure><img src="image_632.png" /></figure>

Probability that observation 120 belongs to class 2 : 0
Need to formalize

<figure><img src="image_633.png" /></figure>

Formalizing estimation

<figure><img src="image_634.png" /></figure>

Formalizing selection

<figure><img src="image_635.png" /></figure>

More advanced formalizing

<figure><img src="image_636.png" /></figure>

Experiments

<figure><img src="image_637.png" /></figure>

To go further

<figure><img src="image_638.png" /></figure>

## Advanced

Visualize in a Gaussian-like way, and onto R2, results of Gaussian or non-Gaussian based clustering.

[19]: library(ClusVis)

[20]: logTik <- getTik(resK, log = TRUE)
prop <- getProportion(resK)
resVisu <- clusvisu(logTik, prop)

## Component Interpretation

[21]: plotDensityClusVisu(resVisu, addobs = FALSE)

<figure><img src="image_639.png" /></figure>

## Observation Scatter-plot

[22]: plotDensityClusVisu(resVisu, addobs = TRUE)

Difference between entropies: 0.22

<figure><img src="image_640.png" /></figure>
<figure><img src="image_641.png" /></figure>

<figure><img src="image_642.png" /></figure>

<figure><img src="image_643.png" /></figure>

<figure><img src="image_644.png" /></figure>

<figure><img src="image_645.png" /></figure>

<figure><img src="image_646.png" /></figure>

Example 2: Canadian weather data<fnref n="13" />

* This dataset is extracted from the fda package (Ramsay, Wickham, Graves, and Hooker 2018). It contains daily *temperature* and *precipitation* at 35 different locations in Canada averaged over 1960 to 1994.

* Data are contained in a list of 3 elements: `tempav`, a 365×35 matrix where one column corresponds to the daily temperature, `precav` a 365 × 35 matrix where one column corresponds to the daily precipitation and time, a vector containing the index of the day.

<footnote n="13">¹³Ramsay JO, Silverman BW (2005). Functional Data Analysis. Springer Series in Statistics, 2nd edition.
Springer. ISBN 038740080X.</footnote>
<figure><img src="image_647.png" /></figure>

<figure><img src="image_648.png" /></figure>

<figure><img src="image_649.png" /></figure>

<figure><img src="image_650.png" /></figure>

<figure><img src="image_651.png" /></figure>

<figure><img src="image_652.png" /></figure>

## Visualize the data

### Temperature

```R
plottemp <- plottemp(CanadianWeatherData, CanadianWeatherTemp, type = "1")
```

<figure><img src="image_653.png" /></figure>

### Precipitation

```R
plotprecip <- plotprecip(CanadianWeatherData, CanadianWeatherPrecip, type = "1")
```

<figure><img src="image_654.png" /></figure>

## Transform the data

The data must be converted in the RlmStComp format, which is for one functional

```
"time_value_1, time_value_2, ..., time_value_n"
```

```R
plot[18] data <- list(temper = apply(CanadianWeatherTemp, 2, function(x) createFunctional(CanadianWeatherTemp, x)))
```

Which gives for the temperatures of the first city:

```R
plot[11] data <- plotTemp[1]
```

St. Johns :

```
1: -3.6,2.3-3.1,3.3-3.4,4.4-5.2,8.6-4.5,7.5,8.3-3.1,9.4-10.5-11.4,8.12-5.2,13.5-5.14-5.4,15.5-4.4,16.4-6.17-5.9,18.5-19.4,8.2
```
<figure><img src="image_655.png" /></figure>

<figure><img src="image_656.png" /></figure>

<figure><img src="image_657.png" /></figure>

<figure><img src="image_658.png" /></figure>

<figure><img src="image_659.png" /></figure>

<figure><img src="image_660.png" /></figure>

<figure><img src="image_661.png" /></figure>

<figure><img src="image_662.png" /></figure>

<figure><img src="image_663.png" /></figure>

<figure><img src="image_664.png" /></figure>

## Clustering with RMixtComp

Define the distribution used for each variable and the associated hyperparameters. The functional model requires 2 hyperparameters: `nSub` is the number of subregressions into which the function will be decomposed; `nCoeff` is the number of polynomial coefficients of each subregression (2 = line).

In [12]:
`nSub` <- 4
`nCoeff` <- 2
`func <- list(type = 'Fum_CS', paramStr = paste0('nSub: ', nSub, ', nCoeff: ', nCoeff))`
`model <- list(tempav = func, precav = func)`

Define the SEM algorithm's parameters

In [13]:
`algo <- createAlgo()`

Choose the desired number of classes and the number of runs for each given number of classes.

In [14]:
`nClass <- 2*4`
`nRun <- 3`

In [15]:
`res <- mixtCompLearn(data, model, algo, nClass = nClass, criterion = 'ICL', nRun = nRun, nCore = 2)`

```r
# Run hierarchical MixtComp in learn mode with 3 run(s) per number of classes and 2 core(s)
Data: 35 individuals and 3 variables.
```

```r
-- K = 1
Run times: 32.274s
-- K = 2
```

Number of splits to perform: 1
Split a cluster in two
Compute criterion

Run times: 32.655s

```r
-- K = 3
```

Number of splits to perform: 2
Split a cluster in two
Split a cluster in two
Compute criterion

Run times: 96.574s

Number of splits to perform: 2
Split a cluster in two
Split a cluster in two
Compute criterion

Run times: 96.574s

```r
-- K = 4
```

Number of splits to perform: 2
Split a cluster in two
Split a cluster in two
Compute criterion

Run times: 86s

Total runtime: 247.804s

Best model according to ICL: 4 clusters.
<figure><img src="image_665.png" /></figure>

<figure><img src="image_666.png" /></figure>

<figure><img src="image_667.png" /></figure>

<figure><img src="image_668.png" /></figure>

<figure><img src="image_669.png" /></figure>

<figure><img src="image_670.png" /></figure>

# Output's Analysis

## Criterion

This chart represents the criterion value (choose between BIC and ICL) for each model that was built. The higher the value (close to 0) the better the model.

In [16]: plotCrit(res)

<figure><img src="image_671.png" /></figure>

Choose the number of classes to study in the following.

In [17]: K <- 3
resK <- extractMixtCompObject(res, K)
<figure><img src="image_672.png" /></figure>

<figure><img src="image_673.png" /></figure>

<figure><img src="image_674.png" /></figure>

<figure><img src="image_675.png" /></figure>

<figure><img src="image_676.png" /></figure>

<figure><img src="image_677.png" /></figure>

## Variables

This chart represents the discriminating level of each variable. A high value (close to one) means that the variable is highly discriminating. A low value (close to zero) means that the variable is poorly discriminating.

<figure><img src="image_678.png" /></figure>

This graph displays the similarity between every pair of variables. A high value (close to one) means that the two variables provide the same information for the clustering task (i.e. similar partitions). A low value (close to zero) means that the two variables provide some different information for the clustering task (i.e. different partitions).

<figure><img src="image_679.png" /></figure>

Draw the distribution of the two variables.

<figure><img src="image_680.png" /></figure>

<figure><img src="image_681.png" /></figure>

Mean curves and 95%-level confidence intervals per class

<figure><img src="image_682.png" /></figure>
<figure><img src="image_683.png" /></figure>

<figure><img src="image_684.png" /></figure>

<figure><img src="image_685.png" /></figure>

<figure><img src="image_686.png" /></figure>

<figure><img src="image_687.png" /></figure>

<figure><img src="image_688.png" /></figure>

## Classes

This chart shows the proportion of individuals in each class.

In [221]: plotFpropClassplot()

Proportion

<figure><img src="image_689.png" /></figure>

This chart represents the connectivity level between each pair of classes. A high value (blue to red) means that the 2 classes are strongly overlapping (i.e. low overlapping). A low value (red to blue) means that the 2 classes are similar for the clustering task (i.e. high overlapping).

In [221]: plotConnClassplot()

Connectivity between classes

<figure><img src="image_690.png" /></figure>

The graph displays the discriminating level of each variable for the selected class.

In [221]: plotDiscVarplot() variable = "class"

In [221]: plotDiscVarplot() class = class

Discriminant level of the class in rows

<figure><img src="image_691.png" /></figure>

This chart summarizes the distribution of the selected variable for the selected class.

In [221]: plotDetailThread() variable = class = class, get = TRUE

Mean values and 95% confidence intervals per class

<figure><img src="image_692.png" /></figure>
<figure><img src="image_693.png" /></figure>

<figure><img src="image_694.png" /></figure>

<figure><img src="image_695.png" /></figure>

<figure><img src="image_696.png" /></figure>

<figure><img src="image_697.png" /></figure>

<figure><img src="image_698.png" /></figure>

Probabilities

<figure><img src="image_699.png" /><figcaption>This chart shows the probability of assignment to a class for each individual. Individuals have been reordered in decreasing assignment probability.</figcaption></figure>

Advanced

<figure><img src="image_700.png" /><figcaption>Visualizes in a Gaussian-like way, and onto R2, results of Gaussian or non-Gaussian based clustering.</figcaption></figure>

Component interpretation

<figure><img src="image_701.png" /></figure>

<figure><img src="image_702.png" /></figure>
<figure><img src="image_703.png" /></figure>

<figure><img src="image_704.png" /></figure>

<figure><img src="image_705.png" /></figure>

<figure><img src="image_706.png" /></figure>

<figure><img src="image_707.png" /></figure>

<figure><img src="image_708.png" /></figure>

## Graphic visualization on a map with $K = 4$

<figure><img src="image_709.png" /><figcaption>Figure 6: Geographic positions of the stations. The colour indicates the class: black = 1, red = 2, green = 3 and blue = 4.</figcaption></figure>
<figure><img src="image_710.png" /></figure>

<figure><img src="image_711.png" /></figure>

<figure><img src="image_712.png" /></figure>

<figure><img src="image_713.png" /></figure>

<figure><img src="image_714.png" /></figure>

<figure><img src="image_715.png" /></figure>

## Outline

1. Need to formalize

2. Formalizing estimation

3. Formalizing selection

4. More advanced formalizing

5. Experiments

6. To go further
<figure><img src="image_716.png" /></figure>

<figure><img src="image_717.png" /></figure>

<figure><img src="image_718.png" /></figure>

<figure><img src="image_719.png" /></figure>

<figure><img src="image_720.png" /></figure>

<figure><img src="image_721.png" /></figure>

Some remaining questions

- More on dependent data (like times series)

- High-dimensional data

- Missing not at random data (MNAR)

- ...
<figure><img src="image_722.png" /></figure>

<figure><img src="image_723.png" /></figure>

<figure><img src="image_724.png" /></figure>

<figure><img src="image_725.png" /></figure>

<figure><img src="image_726.png" /></figure>

<figure><img src="image_727.png" /></figure>

Next lesson

Clustering : une vision unifiée pour une utilisation éclairée
**Traitement de la grande dimension & co-clustering**