Les méthodes de boosting

add.png

\bullet Présentation:

Le boosting est une famille d’outils d’analyse discriminante très puissants. L’idée est d’améliorer la règle de discrimination d’un classifieur sur plusieurs sous-parties de l’hyperplan formé par (Y, X ^1, \cdots, X ^P) et d’établir la classification finale au travers d’un vote global pondéré. Chaque partie du plan se voit alors attribuer une pondération en fonction de sa performance de discrimination et qui est mise à jour itération par itération afin d’adapter la frontière décisionnelle en fonction des zones les plus difficiles à classer.

Le classifieur utilisé est la plupart du temps la régression linéaire, logistique ou les arbres de décision. Ils portent le nom de « classifieurs fiables » à la fois car ils sont simples d’utilisation mais également car l’objectif du boosting est d’améliorer la règle de discrimination lorsqu’elle est peu performante.

Le boosting appartient à la catégorie des méthodes d’apprentissage statistique nécessitant un jeu de d’apprentissage pour la conception de la règle d’apprentissage et un jeu de test pour sa validation afin de limiter les risques de sur-apprentissage. Comme la plupart d’entre elles, le boosting a le principal défaut de concevoir une règle décisionnelle peu lisible d’où son surnom de boîte noire.

Parmi les méthodes de boosting, la plus célèbre reste l’algorithme AdaBoost conçu par Yoav Freun et Robert Shapire en 1997 et définit pour la discrimination d’une variable réponse Y binaire à partir d’un ensemble de variables continues explicatives \mathbf{X} = (X ^1, \cdots, X ^P).

D’autres algorithmes de boosting sont particulièrement répandus, nous pouvons ainsi citer: l’algorithme Gentle Adaboost, Real Adaboost, MART et RankBoost.

\bullet Les algorithmes de Boosting:

Hypothèse préliminaire: Y binaire ou continue, \mathbf{X} continue.

L’algorithme générale:

Avant d’aborder les algorithmes Adaboost et MART, nous présentons l’algorithme de boosting de base.

– Étape initiale: définir un classifieur faible f et un nombre T d’itérations d’entrainement pour l’optimisation des pondérations.

– Étape 1: Soit \mathbf{E}_1 = \mathbf{E} = [Y \mathbf{X}], la matrice de données d’origine et pour laquelle chaque observation a un poids similaire aux autres et égal à 1.

– Étape 2: Pour t \in \lbrace 2, \cdots, T \rbrace,

—- Nous appliquons notre classifieur faible: h_t = f(\mathbf{E}_1)

—- Nous mesurons le taux d’erreurs de classification de la règle décisionnelle conçue précédemment: \epsilon_t = Pr_{x ~ \mathbf{E}_t, y} I[h_t(x) \neq y]

—- Nous mettons à jour les pondérations et chaque observation mal classée se voit attribuer un poids plus important. Ainsi, \mathbf{E}_t devient \mathbf{E}_{t + 1}

– Étape 3: Nous déterminons la règle décisionnelle finale: H(x) = \mbox{Vote a partir des } h_t(x), \forall t \in [1,T]

L’algorithme AdaBoost

L’Algorithme AdaBoost a été pensé pour le cas où Y est binaire. Il se décline de la manière suivante:

– Étape initiale: initialisation des pondérations associées aux observations,

w_i ^1 = \frac{1}{N}, \forall i \in [1,N]

– Étape 1: \forall t \in [2,T]

—- Application du classifieur paramétré sur le jeu d’apprentissage en utilisant les pondérations w_i ^{t - 1}, \forall i \in [1,N], nous obtenons la règle décisionnelle f_t(x).

—- Calcul du taux d’erreurs de la règle décisionnelle définie précédemment,

\epsilon_t = Pr_{x ~ \mathbf{E}_t, y} I[f_t(x) \neq y]

—-Si \epsilon_t >0.5 alors le classifieur fait mieux qu’une règle aléatoire et nous conservons le modèle, sinon nous ajustons,

  • le coefficient de pondération:

\alpha_t = \frac{1}{2} ln(\frac{1 - \epsilon_t}{\epsilon_t})

  • les pondérations:

w_i ^t = w_i ^{t - 1} e ^{\alpha_t \cdot I_{(y_i \neq f_t (x_i))}}, \forall i \in [1,N]

– Étape 2: Nous obtenons le modèle final,

H(x) = sign[\sum_{t = 1} ^T \alpha_t f_t (x)]

L’algorithme Multiple Additive Regression Trees (MART)

L’Algorithme MART se base sur l’utilisation d’arbres de régression. Par conséquent, il est adapté au cas où Y est continue. Les différentes étapes d’estimation sont les suivantes:

– Étape initiale: Initialiser f_0 (x) = arg min_{\delta} \sum_{i = 1} ^N L(y_i, \delta), avec L fonction de perte à minimiser. Et détermination du nombre d’itérations T.

– Étape 1: Pour t \in [1, T],

— Pour i \in [1, N] calculer,

r_{i,t} = - [\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}]_{f = f_{t - 1}}

— Estimer le modèle par arbre de régression sur r_{i,t} et donnant les régions décisionnelles R_{j,t}, \forall j \in \lbrace 1, \cdots, J_t \rbrace

— Pour j \in \lbrace 1, \cdots, J_t \rbrace, calculer,

\delta_{j,t} = arg min_{\delta} \sum_{x_i \in R_{j,t}} L(y_i, f_{t - 1} (x_i) + \delta)

— Mettre à jour f_t (x) = f_{t - 1} (x) + \sum_{j = 1} ^{J_t} \delta_{j,t} I (x \in R_{j,t})

– Étape 2: Nous avons enfin la règle décisionnelle finale de forme: \hat{f} (x) = f_T (x)

\bullet Annexe théorique:

Cette partie de l’article présente une esquisse de la justification de l’usage de la fonction de perte exponentielle.

L’algorithme AdaBoost se base sur l’idée nouvelle d’une fonction de perte exponentielle dont le principal avantage est computationnel. Ainsi, nous montrons facilement que,

f ^* (x) = arg min_{f(x)} E_{Y|x} (e ^{-Y f(x)}) = \frac{1}{2} log \frac{P(Y = 1 | x)}{P(Y = -1 | x)}

Et équivalent à,

P(Y = 1 | x) = \frac{1}{1 + e ^{-2 f^* (x)}}

De cette manière, l’expansion additive produite par l’algorithme AdaBoost revient à estimer une partie du log-odds de P(Y = 1 | x). Ce qui justifie l’usage de son signe comme règle de classification décrite par la formule,

G(x) = sign (\sum_{m = 1} ^M \alpha_m G_m (x))

Un autre critère de minimisation de perte est la log-vraisemblance de la loi binomial négative (appelée également déviance ou entropie croisée), interprétant f comme une transformation logit. Nous avons,

p(x) = P(Y = 1 | x) = \frac{e ^{f(x)}}{e ^{-f(x)} + e ^{f(x)}}  = \frac{1}{1 + e ^{-2 f(x)}}

, avec Y ' = \frac{Y + 1}{2} \in \lbrace 0, 1 \rbrace. La fonction de perte est alors,

l(Y,p(x)) = Y ' log (p(x)) + (1 - Y ') log(1 - p(x))

Et de déviance,

-l (Y, f(x)) = log (1 + e ^{-2Y f(x)})

De la formule précédente, nous voyons que,

E_{Y|x} [-l(Y,f(x))] = E_{Y|x} [e ^{-Y f(x)}]

\bullet Application informatique:

Procédure SAS: ND.

Package et fonction Rhttps://www.rdocumentation.org/packages/JOUSBoost/versions/2.1.0/topics/adaboost

\bullet Bibliographie:

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– The top ten algorithms in data mining de Xindong Wu et Vipin Kumar.

– Probabilité, analyse des données et Statistique de Gilbert Saporta.

– Data mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry.

Le modèle additif généralisé

add.png

\bullet Présentation:

Le modèle additif généralisé, que nous pouvons retrouvons sous l’acronyme GAM, a été élaboré par Trevor Hastie et Rob Tibshirani en 1990.  Le modèle GAM se base sur la somme (doù le terme d’additif) de fonctions f_1, \cdots, f_P de transformation des différentes variables explicatives \mathbf{X} = (X ^1, \cdots, X^P).

La variable à expliquer Y peut être de tout format et les variables explicatives \mathbf{X} = (X ^1, \cdots, X ^P) doivent êtres continues.

Le modèle GAM nécessite un jeu d’apprentissage pour sa construction et un jeu de test pour sa validation car il est sujet à un possible fort sur-apprentissage.  En dépit de sa puissance avérée pour modéliser un phénomène, il reste peu conseillé si nous disposons d’un nombre P de variables explicatives trop important du fait d’un coût de calcul très élevé.

\bullet Le modèle additif généralisé:

Hypothèse préliminaire: \mathbf{X} continue.

Le modèle:

La forme générale du modèle est,

E[Y | X ^1, \dots, X ^P] = \alpha + \sum_{p = 1} ^P f_p (X ^p)

Avec f_p, \forall p \in [1, P], fonction de transformation associée à la variable X ^p ou au uplet de variables explicatives. Deux principaux algorithmes d’estimation des paramètres sont utilisés: l’algorithme Backfitting et l’algorithme général local de scoring,

L’Algorithme backfitting (Friedman et Stuezle, 1981):

– Etape d’initialisation: poser f_0 = E[Y] = \frac{1}{n} \sum_{i = 1} ^n Y_i, et f_p ^0 = 0, \forall p \in[1, P]

– Etape itérative t: pour p \in [1,P], déterminer,

  • R_p = Y - f_0 - \sum_{k = 1} ^{p - 1} f_k ^t (X ^k) - \sum_{k = p + 1} ^P f_k ^{t - 1} (X ^k)
  • f_p ^t = E[R_p | X ^p]

Tant que,

\frac{1}{n} || Y - s_0 - \sum_{p = 1} ^P s_p ^t (X ^p) || ^2 décroît ou satisfait le critère de convergence

A noter que pour p = 1, \sum_{k = 1} ^{1 - 1} s_k ^t (X ^k) = \sum_{k = 1} ^0 s_k ^t (X ^k) = 0 et pour p = P, \sum_{k = P + 1} ^P s_k ^t (X ^k) = 0.

De plus, l’estimation des différentes fonctions de lissage à l’instant t, f_p ^t, se fait selon la fonction lien retenue (modèle linéaire, logistique, etc).

Les transformations:

Le modèle GAM se base donc sur la transformation de variables ou de combinaisons de variables afin d’optimiser le modèle prédictif construit. A ce titre, trois principales fonctions existent: La fonction Cubic Smoothing Spline, la fonction de régression locale (LOESS) et la fonction Thin-Plate Smoothing Spline.

– La fonction Cubic Smoothing Spline de paramètre \lambda \in [0, + \infty[:

RSS(f,\lambda) = \sum_{i = 1} ^n [y_i - f(x_i)] ^2 + \lambda \int (f '' (t)) ^2 dt

– La fonction LOESS de paramètre le noyau K:

min_{\alpha(x_0), \beta(x_0)} \sum_{i = 1} ^n K_{\lambda} (x_0,x_i) [y_i - \alpha(x_0) - \beta(x_0) x_i] ^2

– La fonction Thin-Plate Smooting Spline:

f(x) = \beta_0 + \beta ^t x + \sum_{i = 1} ^n \alpha_i h_i (x), h_i (x) = \theta(|| x - x_p||), \theta(z) = z ^2 log z ^2

\bullet Annexe théorique:

Cette partie de l’article présente une esquisse de la démonstration des fonctions à estimer pour déterminer les valeurs de \mathbf{X} selon la projection f = (f_1, f_2, \cdots).

– Nous avons vu que la fonction Cubic Smoothing Spline de paramètre \lambda \in [0, + \infty[ était définie ainsi:

RSS(f,\lambda) = \sum_{i = 1} ^n [y_i - f(x_i)] ^2 + \lambda \int (f '' (t)) ^2 dt

Lorsque \lambda = 0 alors f peut être n’importe quel fonction d’interpolation de X. Si \lambda = \infty, alors f est la fonction des moindres carrés partiels et aucune derivée seconde ne peut être obtenue.

Pour obtenir une version simplifiée de la fonction Cubic Smoothing Spline, nous posons D_j(x) un sous-ensemble de D fonctions basiques et \theta un vecteur de paramètres à estimer. La fonction d’origine peut alors se réécrire:

RSS(\theta,\lambda) = (y - D \theta) ^t (y - D \theta) + \lambda \theta ^t \Omega_N \theta

, avec \lbrace D \rbrace_{i,j} = D_j (x_i) et \lbrace \Omega_N \rbrace_{j,k} = \int D_j '' (t) D_k '' (t) dt. La solution est alors,

\hat{\theta} = (D ^t D + \lambda \Omega_D) ^{-1} D ^t y

, qui s’apparente à la régression ridge généralisée. Après développement, nous obtenons ainsi la forme simplifiée de f,

\hat{f}(x) = \sum_{j = 1} ^D D_j (x) \hat{\theta}_j

– Nous avons vu que la fonction LOESS de paramètre le noyau K était définie ainsi:

min_{\alpha(x_0), \beta(x_0)} \sum_{i = 1} ^n K_{\lambda} (x_0,x_i) [y_i - \alpha(x_0) - \beta(x_0) x_i] ^2

La fonction d’estimation est de la forme,

\hat{f} (x_0) = \hat{\alpha} (x_0) + \hat{\beta} (x_0) x_0

Nous définissons alors la fonction du vecteur de valeurs: b(x) ^t = (1,x), \mathbf{B} la matrice de régression de taille N \times 2 composée des b(x_i) ^t et \mathbf{W}(x_0) la matrice diagonale de taille N \times N et dont l’élément i est K_{\lambda} (x_0, x_i). Nous avons alors,

\hat{f} (x_0) = b(x_0) ^t \mathbf{B} ^t \mathbf{W} (x_0) \mathbf{B} ^{-1} \mathbf{B} ^t \mathbf{W} (x_0) y = \sum_{i = 1} ^N l_i (x_0) y_i

Ce qui nous donne une expression explicite pour l’estimation de la LOESS.

– Nous avons vu que la fonction Thin-Plate Smooting Spline était définie ainsi:

f(x) = \beta_0 + \beta ^t x + \sum_{i = 1} ^n \alpha_i h_i (x)

, avec h_i (x) = \theta(|| x - x_p||) et \theta(z) = z ^2 log z ^2.

Soit \mathbf{X} \in R ^2 et les fonctions de base de la forme h_{1,k} (X ^1), \forall k \in 1, \cdots, M_1 des coordonnées de X ^1 et de même pour la fonction h_{2,k} (X ^2), \forall 1, \cdots, M_2. Nous définissons alors,

g_{j,k} (\mathbf{X}) = h_{1,j} (X ^1) h_{2,k} (X ^2), j = 1, \cdots, M_1 et k = 1, \cdots, M_2

, que nous pouvons généraliser avec la forme suivante,

g(\mathbf{X}) = \sum_{j = 1} ^{M_1} \sum_{k = 1} ^{M_2} \theta_{j,k} g_{j,k} (\mathbf{X})

Pour estimer cette fonction, l’idée est de résoudre le problème suivant,

min_J \sum_{i = 1} ^N \lbrace y_i - f(x_i) \rbrace ^2 + \lambda J[f]

, avec J fonction de pénalisation visant à stabiliser f.

De plus, si le paramètre \lambda \rightarrow 0 alors la solution s’approche d’une fonction d’interpolation et si \lambda \rightarrow \infty alors elle s’approche d’un plan des moindres carrés. Pour les valeurs intermédiaires de \lambda, la solution peut être représentée comme une fonction linéaire dont les coefficients sont obtenus par régression ridge.

La solution est alors de la forme,

f(x) = \beta_0 + \beta^t x + \sum_{j = 1} ^N \alpha_j h_j (x)

, où h_j(x) = \eta(||x - x_j||) et \eta(z) = z ^2 log z ^2.

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_gam_sect004.htm

Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/gam.html

\bullet Bibliographie:

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– Data Mining et statistique décisionnelle. L’Intelligence des données de Stéphane Tufféry.

https://www.math.univ-toulouse.fr/~besse/pub/Appren_stat.pdf

La classification ascendante hiérarchique des variables

add

\bullet Présentation:

La classification ascendante hiérarchique des variables, élaborée par Michael R. Anderberg en 1973, comprends un large panel d’outils de classification agissant directement sur un ensemble de variables (X ^1, \cdots, X ^P) continues.

Elle permet de s’affranchir de méthodologies un peu plus longues comme celles basées sur l’Analyse en Composantes Principales (ACP) ou encore celles basées sur l’Analyse des Correspondances Multiples qui requièrent, dans un premier temps, la construction des axes factoriels puis l’application d’une Classification Ascendante Hiérarchique pour arriver à un résultat similaire.

L’idée de la classification ascendante hiérarchique des variables est de lancer une série itérative d’ACP, via une rotation oblique et non orthogonale, sur les groupes et sous-groupes de variables conçus au fur et à mesure en fonction de leurs corrélations et, à chaque itération, de redistribuer les variables dans les différents groupes pour optimiser la classification finale. La principale force de cet outil, et qui en donne l’une des propriétés fondamentales, est que la rotation oblique amène la création de variables latentes très faiblement corrélées entre elles.

\bullet La classification de variables

Hypothèse préliminaire: (X ^1, \cdots, X ^P) continues.

L’Algorithme général, que nous retrouvons souvent dans la littérature sous l’appellation varclus, consiste à appliquer les étapes qui sont les suivantes,

Étape 0: Notons E l’ensemble de départ contenant toutes les variables explicatives X ^1, \cdots, X ^P et « ACP quartimax » la version de l’Analyse en Composantes Principales basée sur la rotation oblique.

Étape 1:

–  Exécuter une ACP quartimax sur E et sortir les deux valeurs propres les plus grandes si elles sont supérieures à 1. Si ce n’est pas le cas l’algorithme s’arrête.

– Pour chaque variable, récupérer le coefficient associé aux deux premières composantes et déterminer une première version du groupe d’attribution en fonction de la valeur absolue maximale. Nous obtenons alors G_1 ', G_2 ' les deux sous-groupes issus de notre première séparation de E.

– Réattribuer, si c’est possible, les variables misent dans G_2 ' (soit celui dont les coefficients étaient plus grands sur la seconde composante que sur la première) à G_1 ' (donc associé à la première composante). Pour cela il faut appliquer, pour chacune des variables parcourues une à une X ^p, p \in G_2 ', la procédure suivante,

  • Lancer une ACP quartimax sur G_1 ' pour obtenir la valeur propre \lambda_{G_1 '}.
  • Lancer une nouvelle ACP quartimax sur le groupe composé des variables de G_1 ' et de la variable en cours de test X ^p, nous retenons \lambda_p valeur propre associée.

– Nous laissons alors dans G_2 ' la variable X ^p dont la valeur propre \lambda_p était la moins contributive et réaffectons les autres variables au groupe G_1 '. Cette manipulation se justifie par le fait que l’algorithme varclus est de type classification ascendante hiérarchique. A noter, que si \lambda_p < \lambda_{G_1' } alors la variable associée reste dans G_2 '. Cette manipulation terminée, nous obtenons nos groupes G_1, G_2 finaux.

Étape i > 1: Réitérer l’étape une sur les sous-groupes G_1 et G_2 tant que la seconde valeur propre est supérieure à 1, sinon le groupe n’est pas divisé et l’algorithme s’arrête.

L’objectif de l’algorithme est la recherche, par étapes et divisons successives, de la ré-affectation des différentes variables aux groupes ou sous-groupes permettant de maximiser la variance expliquée.

La puissance de l’approche, comparée à celles basées sur les méthodologies ACP+CAH ou ACM+CAH, est que nous substituons la distance euclidienne par le coefficient de corrélation (de Pearson, Spearman ou Kendal) et gagnons en puissance si les liens entre nos variables sont exclusivement linéaire.

\bullet Annexe théorique:

Cette partie de l’article présente la démonstration de l’inertie totale du nuage de points de l’ACP.

Sa définition générale et mathématique en un point a est,

I_a = \sum_{i = 1} ^n p_i (e_i - a) ' \mathbf{M} (e_i - a)

Où,

\mathbf{M} = \begin{pmatrix} \frac{1}{s_1 ^2} & & & & & 0 \\ & \frac{1}{s_2 ^2} & & & & \\ & & . & & \\ & & & & . & \\ 0 & & & & & \frac{1}{s_P ^2} \\ \end{pmatrix}

En notant i_G l’inertie totale (au centre de gravité G), nous avons par la relation de Huyghens,

I_a = I_G + (G - a) ' \mathbf{M} (G - a) = I_G + || G - a || ^2

A noter que si G = 0,

I_G = \sum_{i = 1} ^n p_i e_i ' \mathbf{M} e_i

Nous avons alors,

2 I_g = \sum_{i = 1} ^n \sum_{j = 1} ^n p_i p_j (e_i - e_j)' \mathbf{M} (e_i - e_j) = \sum_{i = 1} ^n \sum_{j = 1} ^n p_i p_j || e_i - e_j || ^2

, soit la moyenne des carrés de toutes les distances entre les n individus.

Enfin, l’inertie totale peut également être définie comme la trace de la matrice \mathbf{M V}. Et comme p_i e_i ' \mathbf{M} e_i est un scalaire nous avons par propriété de commutativité sous la trace,

I_g = trace(\sum_{i = 1} ^n p_i e_i ' \mathbf{M} e_i)

= trace(\sum_{i = 1} ^n \mathbf{M} e_i ' p_i e_i)

= trace(\mathbf{M} \mathbf{X} ' \mathbf{D} \mathbf{X})

= trace(\mathbf{M} \mathbf{V})

\bullet Exemple:

Soit le jeu de données suivants,

add

Nous débutons l’algorithme par une ACP quartimax sur l’ensemble complet de variables explicatives. Nous obtenons deux valeurs propres supérieures à 1:

\lambda_1 = 2.072474087, \lambda_2 = 1.07269809

Nous pouvons donc lancer une première itération, ou première division. Lors de la manipulation précédente, nous obtenons également les coefficients suivants associés à nos deux premières composantes,

\begin{tabular} {|l|c|c|} \hline & Composante 1 & Composante 2 \\ \hline X1 & 0.79661 & 0.53644 \\ \hline X2 & -0.25751 & 0.64743 \\ \hline X3 & 0.08645 & 0.55907 \\ \hline X4 & 0.089118 & -0.02693 \\ \hline X5 & 0.29312 & 0.73244 \\ \hline \end{tabular}

Nous déterminons ainsi deux groupes,

– celui des variables dont la plus forte contribution est celle pour le premier axe factoriel et qui sont: X ^1, X ^4,

– et celui des variables dont la plus forte contribution est celle pour le second axe factoriel et qui sont: X ^2, X^3, X ^5.

Avant de passer à la seconde étape, nous devons filtrer la variable qui apporterait le moins d’information au premier groupe si nous l’y avions incluse. Nous lançons donc une nouvelle ACP quartimax sur le premier groupe X ^1, X ^4 et nous obtenons comme valeur propre:

\lambda_1 = 1.6229470

Maintenant, regardons celles obtenues après l’apport de chaque variable,

– Valeur propre obtenue suite à l’ACP quartimax sur X ^1, X ^4, X ^2:

\lambda_{X^1, X^2,X^4} = 1.6407263

– Valeur propre obtenue suite à l’ACP quartimax sur X ^1, X ^4, X ^3:

\lambda_{X1,X^3,X^4} = 1.7163887

– Valeur propre obtenue suite à l’ACP quartimax sur X ^1, X ^4, X ^5:

\lambda_{X^1,X^4,X^5} = 1.9245231

Étant donné que nous cherchons à classer nos variables une par une, il en ressort que la variable qui apporte le moins à la restitution de la variance est X ^2. Notre première division finalisée sera donc X ^2 pour le groupe 1 versus X^1, X^3, X^4, X^5 pour le groupe 2.

Maintenant, relançons une nouvelle ACP quartimax sur le sous-ensemble X^1, X^3, X^4, X^5 soit le groupe 2. Nous obtenons les deux valeurs propres suivantes,

\lambda_{2,1} = 2.034982 > 1, \lambda_{2,2} = 0.9987 < 1

Seule la première valeur propre \lambda_{2,1} est supérieur à 1, la seconde \lambda_{2,2} ne l’étant pas, nous stoppons là l’algorithme. Il ne nous reste plus qu’à calculer les parts de variance expliquée afin de déterminer la longueur des branches de liaisons.

La première valeur propre était \lambda_1 = 2.072474087. Il nous faut désormais celles associées aux deux groupes construits dont nous savons déjà que l’une d’elle est \lambda_{2,1} = 0.2034982. Mais il nous faut encore celle de l’autre groupe (qui n’est pas \lambda_{2,2} car associée à la seconde composante quand nous voulons celle de ce groupe sur la première composante). Mais étant donné que le groupe en question est celui basé sur l’unique variable X ^2, forcément sa valeur propre est 1. Nous en déduisons que,

\lambda_2 = 2.034982 + 1 = 3.034982

Enfin, pour déterminer la longueur des branches, il nous faut les parts de variances restituées suivantes,

br_1 = \frac{2.072474087}{2.072474087+3.034982} = \frac{2.072474087}{5.107456087} = 0.40577

br_2 = \frac{3.034982}{2.072474087+3.034982} = \frac{3.034982}{5.107456087} = 0.59423

Nous avons ainsi tous les éléments pour la construction du dendrogramme qui est le suivant,

add

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#varclus_toc.htm

Package et fonction R: https://www.rdocumentation.org/packages/ClustOfVar/versions/0.8/topics/ClustOfVar-package

\bullet Bibliographie:

– Cluster analysis for applications de Michael R. Anderberg.

– What are the sensory differences among coffees? Multi-panel analysis of variance and flash analysis de Pascal Schlich.

–  Data mining et statistique décisionnelle de Stéphane Tufféry.

– Classification de variables. Classification autour de variables latentes de Ricco Rakotomalala.

Le positionnement multidimensionnel

\bullet Présentation:

Le positionnement multidimensionnel ou échelonnement multidimensionnel, élaboré par Joseph Bernard Kruskal en 1964, est un outil de la famille de l’analyse exploratoire permettant d’étudier les similarités ou dissimilarités au sein d’une matrice \mathbf{X} = (X ^1, \cdots, X ^P) de P variables continues.

Le positionnement multidimensionnel, connu sous l’appellation MDS également, se concentre exclusivement sur la projection des observations dans un plan visuellement interprétable.

Elle peut être vue comme une généralisation voir une extension de l’Analyse en Composantes Principales (ACP) ou de l’Analyse des Correspondances Multiples (ACM) car basée sur le même procédé mais permettant la gestion de relations non linéaires. Cette propriété provient du fait qu’elle travaille directement sur la matrice de similarité/dissimilarité issue de \mathbf{X} qui peut être paramétrée selon un large panel de fonctions de lien.

\bullet Le positionnement multidimensionnel:

Hypothèse préliminaire: Aucune.

L’Algorithme MDS et assez simple voir classique d’utilisation et se déroule selon les étapes suivantes,

– Étape 1: Construction de la matrice de similarité/dissimilarité \mathbf{D_X} de taille n \times n selon une fonction lien fixée et qui peut être linéaire ou non linéaire.

– Étape 2: Élévation au carré de \mathbf{D_X}. A noter que nous aurons besoin de la version symétrique de cette matrice, qui l’est bien évidemment par construction mais dont la représentation se résume souvent aux éléments triangulaires inférieurs exclusivement.

– Etape 3: Sortie des valeurs propres \lambda et des vecteurs propres \mathbf{V} associés à \mathbf{D_X} ^2.

– Etape 4: Calcul des projections des observations en multipliant les facteurs par la racine carré de leur valeur propre associées,

\mathbf{Proj} = (V ^1 \sqrt{\lambda_1}, \cdots, V ^n \sqrt{\lambda_n})

L’Interprétation se fait comme pour une ACP ou une ACM. A savoir que plus les observations sont proches et plus elles décrivent un profil commun. Enfin, les observations ou groupes observations s’opposent en fonction de leur coordonnée sur les axes retenus.

\bullet Annexe théorique:

La partie suivante de l’article présente la notion de valeurs propres et de vecteurs propres.

Soit E un espace vectoriel sur un corps commutatif K et donc munit des opérateurs fondamentaux: l’addition, la soustraction, la multiplication et la division. Par propriété, les éléments de E sont alors des vecteurs tandis que ceux de K sont des scalaires. Et soit f un endomorphisme de forme,

f: E \rightarrow E

Une valeur propre est alors: Un scalaire noté \lambda associé àf tel que,

\exists x \ne 0, f(x) = \lambda x

Les valeurs propres d’une matrice carrée \mathbf{A} de taille P \times P sont les valeurs propres de l’endomorphisme de K ^P de la matrice \mathbf{A} dans la base canonique et dont la propriété fondamentale est que les coordonnées de tout vecteur de l’espace dans cette base sont données par les composantes mêmes qui les constituent.

Si E de dimension finie P, les valeurs propres de f sont:

– les racines du polynôme caractéristique commun de forme,

det(\mathbf{X} \cdot Id - f) = det (\mathbf{X} \cdot Id - \mathbf{A})

– toutes les valeurs spectrales de f ou \mathbf{A}.

Les valeurs propres \lambda se déterminent en recherchant le vecteur,

\lambda, det(\mathbf{A} - \lambda Id) = 0

Lors de la manipulation de l’objet ci-dessus, \lambda est identique d’une élément à l’autre de la diagonale, c’est lors de sa résolution que les différentes valeurs possibles de \lambda forment alors le vecteur \lambda_1, \dots, \lambda_P. Il s’agit là de la principale explication pour laquelle, si tous les vecteurs de \mathbf{A} sont indépendants et donc peuvent être vue comme distribués aléatoirement les uns par rapport aux autres, alors \lambda_1 \approx \lambda_2 \approx \cdots \lambda_P.

Un vecteur propre est alors: un vecteur x non nul de E, dit vecteur propre de f et associé à \lambda s’il est de la forme,

f(x) = \lambda x

Les vecteurs propres (associés à une valeur propre \lambda) d’une matrice carrée \mathbf{A} sont les vecteurs propres (associés à la valeur propre \lambda) de l’endomorphisme de K ^P représenté par \mathbf{A}. L’ensemble des vecteurs propres associés à une matrice porte le nom de spectre. Deux propriétés fondamentales sont à considérer,

– un même vecteur propre ne peut être associé à deux valeurs propres différentes,

– une famille de P vecteurs propres associés chacun à une unique valeur propre constitue une famille libre, autrement dit, la seule combinaison linéaire des P vecteurs propres qui donnent le vecteur nul est celle pour laquelle \lambda = 0.

Les vecteurs propres x se déterminent en recherchant le vecteur,

(x ^1, \dots, x ^P), \lambda_p, p \in [1, P],( \mathbf{A} - \lambda Id) \mathbf{X} = 0

La résolution se fait alors valeur propre par valeur propre, consituant ainsi de manière itérative la construction de nos vecteurs propres. La propriété d’affectation unique d’une valeur propre à un vecteur propre provient de cette méthode de calcul, puisque chaque valeur propre sera unique, nous en déduisons que le vecteur propre calculé le sera également.

D’Un point de vue pratique, une valeur propre peut être vue comme une mesure de l’information contenue dans le vecteur propre qui lui est associé. Plus elle est grande, par rapport aux autres valeurs propres, et plus le vecteur propre étire le nuage de point et donc restitue la variance de l’échantillon. A contrario, si les vecteurs de \mathbf{A} sont indépendants, les valeurs propres seront alors très proches et donc les vecteurs propres également, d’où l’absence d’information ou plutôt de discrimination induit par les axes.

\bullet Exemple:

Soit le jeu de données suivant,

add

Nous nous limitons à la recherche de deux axes factoriels afin de simplifier la lecture des résultats finaux.

Dans un premier temps, calculons la matrice de distance selon la norme euclidienne,

\mathbf{D_X} = \begin{pmatrix} Obs & 1 & 2 & 3 & \cdots & 18 & 19 & 20 \\ 1 & 0 & 2.399264 & 7.883224 & \cdots & 21.421600 & 23.647763 & 24.357353 \\ 2 & 2.399264 & 0 & 7.958839 & \cdots & 19.718820 & 21.653041 & 22.264297 \\ 3 & 7.883224 & 7.958839 & 0 & \cdots & 20.082169 & 21.165185 & 22.325087 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 18 & 21.421600 & 19.718820 & 20.082169 & \cdots & 0 & 8.294253 & 9.414303 \\ 19 & 23.647763 & 21.653041 & 21.165185 & \cdots & 8.294253 & 0& 2.783625 \\ 20 & 24.357353 & 22.264297 & 22.325087 & \cdots & 9.414303 & 2.783625 & 0 \\ \end{pmatrix}

, que nous élevons au carré,

\mathbf{D_X} ^2 = \begin{pmatrix} Obs & 1 & 2 & 3 & \cdots & 18 & 19 & 20 \\ 1 & 0 & 5.75647 & 62.14522 & \cdots & 458.88495 & 559.216688 & 593.280624 \\ 2 & 5.75647 & 0 & 63.34313 & \cdots & 388.83187 & 468.854192 & 495.69612 \\ 3 & 62.14522 & 63.34313 & 0 & \cdots & 388.83187 & 468.854192 & 495.69612 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 18 & 458.88495 & 388.83187 & 403.29352 & \cdots & 0 & 68.794629 & 88.629093 \\ 19 & 559.21669 & 468.85419 & 447.96507 & \cdots & 68.79463 & 0& 7.748568 \\ 20 & 593.28062 & 495.69891 & 498.40953 & \cdots & 88.62909 & 7.748568 & 0 \\ \end{pmatrix}

Maintenant, nous pouvons extraire directement les valeurs propres et vecteurs propres associés à nos deux axes d’intérêt,

\lambda = (1232.0742, 603.0851)

, et,

\mathbf{V} = \begin{pmatrix} -0.33694356 &-0.2318659 \\ -0.28077941 & -0.2464647 \\ -0.25183323 & -0.2319835 \\ -0.23399253 & -0.1517208 \\ -0.19099322 & -0.1269469 \\ -0.26124237 & 0.2349641 \\ -0.21332495 & 0.2355275 \\ -0.16415724 & 0.2476648 \\ -0.14203587 & 0.1697586 \\ -0.09004166 & 0.1638683 \\ 0.07253286 & 0.2232178 \\ 0.08083402 & 0.2673542 \\ 0.14760548 & 0.2670581 \\ 0.20451644 & 0.2378346 \\ 0.22376925 & 0.3204969 \\ 0.19139369 & -0.0885641 \\ 0.20316062 & -0.1760915 \\ 0.25231958 & -0.1902502 \\ 0.32202642 & -0.2158966 \\ 0.33651525 & -0.3021356 \end{pmatrix}

Nous réajustons \mathbf{V} par \sqrt{\lambda} en calculant V ^1 \times \sqrt{\lambda_1} et V ^2 \times \sqrt{\lambda_2} et obtenons,

\mathbf{V} ^* = \begin{pmatrix} -11.827027 & -5.694113 \\ -9.855614 & -6.052628 \\ -8.839577 & -5.697002 \\ -8.213352 & -3.725929 \\ -6.704037 & -3.117535 \\ -9.169846 & 5.770199 \\ -7.487901 & 5.770199 \\ -7.487901 & 5.784036 \\ -5.762069 & 6.082101 \\ -4.985589 & 4.168897 \\ -3.160545 & 4.024244 \\ 2.545970 & 5.481737 \\ 2.837348 & 6.565629 \\ 5.181087 & 6.558358 \\ 7.178714 & 5.840693 \\ 7.854505 & 7.870695 \\ 6.718093 & -2.174939 \\ 7.131123 & -4.324418 \\ 8.856648 & -4.672125 \\ 11.303422 & -5.301944 \\ 11.811993 & -7.419784 \\ \end{pmatrix}

Nous avons donc pour chaque observation les coordonnées sur les deux axes factoriels retenus. Nous pouvons alors tracer la carte des individus dans le repère construit,

add

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mds_sect004.htm

Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cmdscale.html

\bullet Bibliographie:

Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis de Joseph Bernard Kruskal.

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard.

L’Analyse d’Hill et Smith

add.png

\bullet Présentation:

Publié en 1976 par Mark Oliver Hill et Anthony John Edwin Smith, l’analyse d’Hill et Smith est une approche multivariée permettant d’étudier les relations multvariées d’un ensemble de variables mixtes continues, ordinales et/ou qualitatives distinctes \mathbf{X} = (X ^1, \cdots, X ^p).

L’objectif de l’analyse d’hill et Smith est de synthétiser l’information en réduisant le nombre de dimensions afin d’avoir une lecture visuelle et simple des interactions se produisant entre les diverses variables du jeu de données.

L’Outil repose sur une approche mixte ACP (Analyse en Composantes Principales), ACM (Analyse des Correspondances Multiples), puisant dans la force des deux approches pour déterminer les différentes relations entre les différentes variables malgré que leur format diverge.

\bullet L’Analyse d’Hill et Smith:

Hypothèse préliminaire: Aucune.

L’Algorithme:

L’Analyse d’Hill et Smith se base donc sur la possibilité de mixer les deux formats principaux de variables: continue et qualitative et de produire des axes factoriels synthétisant les corrélations et liaisons au sein d’un jeu de données. La méthodologie est une sorte de mélange d’ACM et d’ACP, ainsi la première permettra de transformer les variables qualitatives en variables continues et la seconde permettra de construire le lien entre les deux formats de variables.

L’algorithme se déroule ainsi,

– Étape 1: Séparer les variables continues, donnant la matrice \mathbf{A} de taille n \times P_c, et les variables qualitatives, donnant la matrice \mathbf{B} de taille n \times P_q.

– Étape 2: Standardiser \mathbf{A} selon le vecteur de pondération en ligne et de taille 1 \times n,

W_l = (\frac{1}{n}, \cdots, \frac{1}{n})

, à noter que W_l est fixé ainsi par défaut, mais peut être égal aux poids individuels s’ils existent.

– Étape 3: Déterminer la matrice disjonctive complète associée à \mathbf{B} et que nous noterons \mathbf{B} ^D et appliquer la transformation,

\frac{B_{i,j} ^D}{\sum_{i = 1} ^n B_{i,j} ^D} - 1, \forall i \in [1,n], j \in [1, P_q ^* = \sharp \lbrace \mbox{Modalites au sein de } \mathbf{B} \rbrace]

, nous obtenons alors \mathbf{B}_T ^D.

– Étape 4: Reformer la matrice complète (\mathbf{A}, \mathbf{B}_T ^D) et appliquer en premier la pondération en ligne \sqrt{W_l} puis la racine carré de la pondération en colonne,

W_c = (1, \cdots, 1, \frac{\sum_{i = 1} ^n (B_T ^D)_i ^1}{n}, \cdots, \frac{\sum_{i = 1} ^n (B_T ^D)_i ^{P_q ^*}}{n})

, de taille 1 \times (P_c + P_q ^*). Notons, que dans le vecteur W_c, les variables continues ont un poids unitaire contrairement aux variables qualitatives.

L’Application de ces deux pondérations successives nous donne la matrice (\mathbf{A}, \mathbf{B}_T ^D)_P.

– Étape 5: Faire le produit croisé,

(\mathbf{A}, \mathbf{B}_T ^D)_P ^t (\mathbf{A}, \mathbf{B}_T ^D)_P

– Étape 6: Sortir les valeurs propres \lambda et les vecteurs propres associés \mathbf{V}_{HS}.

– Étape 7: Appliquer à \mathbf{V}_{HS} la première pondération en colonne \frac{1}{\sqrt{W_c}} puis la seconde \sqrt{\lambda} pour obtenir les composantes factorielles \mathbf{CF}_{HS}.

Les parts d’information restituée sur chacun des axes factorielles sont décrites au travers du vecteur \lambda. Nos composantes factorielles \mathbf{CF}_{HS} nous donnent les coordonnées de nos différentes variables projetées sur le plan factoriel.

L’Interprétation:

– Pour les modalités des variables qualitatives, la lecture est la même que pour une ACM. La présence d’un groupe en terme de proximité de deux ou plusieurs modalités de variables qualitatives impliquent un lien entre elles et donc que les individus qui ont choisis l’une des modalités du groupe ont également choisis les autres modalités de ce même groupe.

– Pour les variables quantitatives, la lecture est la même que pour une ACP centrée-réduite. Ainsi, un groupe de variables en terme de proximité sur le cercle de corrélation implique une corrélation soit un lien linéaire croissant. Deux groupes de variables opposés sur le cercle de corrélation implique une anti-corrélation soit un lien linéaire décroissant. Deux groupes de variables dont l’angle formé est de 90° implique une indépendance et donc l’absence de corrélation ou d’anti-corrélation entre eux.

– Pour le mixte entre variable(s) continue(s) et variable(s) qualitative(s), il faudra mixer les deux interprétations. Ainsi, si une ou plus variables qualitatives sont au sein d’un groupe de variables continues, cela implique que plus les variables continues croient et plus les individus ont choisis cette ou ces modalité(s).

\bullet Annexe théorique:

Nous présentons ici une esquisse des démonstrations des six propositions justifiant l’aspect théorique de l’analyse de Hill et Smith.

Proposition 1: l’ACP, qui se base sur la matrice des corrélations, est une approche sur la recherche des variables les plus corrélées entre elles.

Démonstration: Soit \mathbf{A} une matrice centrée-réduite (standardisée) de taille n \times P. Nous avons,

 \rho ^2 = \sum_{p = 1} ^P \frac{(Y ^t A ^p) ^2}{\frac{(Y ^t Y)}{P}}

, avec Y vecteur arbitraire des scores et de taille 1 \times n.

Posons \lambda = P \rho ^2, nous obtenons alors les valeurs stationnaires de \rho ^2 lorsque,

\lambda Y = (\sum_{p = 1} ^P (A ^p (A ^p) ^t) Y

, soit la formule d’une composante principale.

Proposition 2: Les composantes principales sont la résultante de la somme de leur propre estimateur des moindres carrés.

Démonstration: L’estimateur des moindres carrés de Y dans la régression linéaire sur la variable A ^p est,

Y_p = (A ^p) ^t Y A ^p = A ^p (A ^p) ^t Y

En utilisant la proposition 1, nous constatons que la composante principale satisfait la relation,

\lambda Y = \sum_{p = 1} ^P Y_p

Proposition 3: Si les A ^p, p \in [1,P] varient entre eux sous contrainte linéaire, alors leur détermination et celle de Y revient à maximiser \rho ^2 sous formulation d’une problème de valeurs propres extraites depuis une matrice carré symétrique.

Démonstration: Soit Y fixe et égal, par exemple, à A ^p. Assumons que Y est standardisé et supposons que A ^p peut s’écrire selon la combinaison linéaire,

A ^p = \theta_1 \zeta_1 + \cdots + \theta_K \zeta_K

, avec (\zeta_1, \cdots, \zeta_K) orthonormal.

Nous avons alors,

\rho_p ^2 = \frac{((A ^p) ^t Y) ^2}{(A ^p) ^t A ^p} = \frac{(\sum_{k = 1} ^K \theta_k \zeta^t Y) ^2}{\sum_{k = 1} ^K \theta_k ^2}

, dont le coefficient de corrélation au carré est maximisé quand,

\theta_k \propto \zeta_k ^t Y

Notons alors \mathbf{C} = (\mathbf{A}, \zeta_1, \cdots, \zeta_K), le coefficient de corrélation au carré est alors maximisé quand,

\theta = \mathbf{C} ^t Y

Par la proposition 2, nous avons que pour \theta  fixé, la moyenne au carré des corrélations est maximisée lorsque \lambda Y est la somme des estimateurs des moindres carrés de Y. Par conséquent, le maximum global de \rho ^2 est obtenu quand,

\theta = \mathbf{C} ^t Y; \lambda Y = \mathbf{C} \theta

, de forme,

\lambda Y = \mathbf{C} \mathbf{C} ^t Y

Proposition 4: Pour \mathbf{B} matrice de variables qualitatives, l’approche revient à un modèle additif sur les modalités des variables, ramenant le problème à une analyse des correspondances.

Démonstration: Pour une classification à un facteur, l’estimation de Y est un simple moyenne des \theta selon la correspondance à Y. D’où les valeurs stationnaires obtenues pour \rho ^2 lorsque,

\theta_r = moyenne de Y pour la modalité r,

P \rho ^2 Y_i = somme des \theta_r, r \in [1,R] pour les caractéristiques modales de l’individu i.

D’où, le triplet (\rho,Y,\frac{\theta}{\rho}) est solution au problème d’analyse des correspondances définie par (0,I) la matrice d’incidence individuel (en ligne) et des catégories (en colonne).

Proposition 5: Les propositions 1, 2, 3, 4 implique que la combinaison de \mathbf{A} et \mathbf{B} et l’étude de leur lien multivariée revient à une analyse multivariée.

Démonstration: La proposition 2 peut être appliquée directement, utilisant un modèle de classification à un facteur pour variables qualitatives et un modèle linéaire pour variables continues.

Proposition 6: La proposition 2 fournit un algorithme itératif direct et pratique pour calculer les solutions de la proposition 5.

Démonstration: Soit,

Y ^* = \sum_{p = 1} ^P Y_p

Selon les notations de la proposition 3, nous pouvons réécrire le problème sous l’itération,

Y ^* = \mathbf{C} \mathbf{C} ^t Y

, dont la convergence peut être obtenue rapidement par la méthode décrite par Clint et Jennings ou celle de Hill.

\bullet Exemple:

Soit le jeu de données suivants,

add

Appliquons l’analyse d’Hill et Smith afin d’obtenir une représentation synthétique des liaisons entre nos neuf variables.

Dans un premier temps, nous séparons les variables continues des variables qualitatives et nous transformons les deux matrices indépendamment,

– pour la matrice restreinte aux variables continues (X ^1, X ^2, X ^3, X ^4, X ^5), nous définissons tout d’abord le vecteur des pondérations des observations,

W_l = (\frac{1}{20}, \cdots, \frac{1}{20}) et de taille 1 \times 20

, et obtenons après standardisation pondérée selon W_l,

\mathbf{A} = (X_s ^1, X_s ^2, X_s ^3, X ^4, X ^5) = \begin{pmatrix} -1.66667688 & 0.53699601 & -0.7559172 & -1.59438477 & -0.24569927 \\ -1.45850552 & 0.82019342 & -0.4080684 & -1.48730298 & 0.02327323 \\ -1.29588196 & -1.60161550 & -1.669232 & -1.37958370 & 0.20543383 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1.33207517 & 0.85113461 & -1.4779835 & -0.02322048 & 1.40414742 \\ 1.47493991 & 0.46699740 & 1.3348975 & 0.15243809 & 1.62667035 \\ 1.63455708 & 0.98718255 & 1.6922264 & -0.14491013 & 1.74313248 \\ \end{pmatrix}

– pour la matrice restreinte aux variables qualitatives (X ^6, X ^7, X ^8, X ^9), nous construisons d’abord la matrice disjonctive complète,

\mathbf{B} ^D = (X ^6 = "A", X ^6 = "B", X ^7 = "A", \cdots, X ^9 = "A", X ^9 = "B", X ^9 = "C")

= \begin{pmatrix} 0 & 1 & 1 & \cdots & 1 & 0 & 0 \\ 0 & 1 & 1& \cdots & 1 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 1 & 0 & \cdots & 0 & 0 & 1 \\ 0 & 1 & 0 & \cdots & 0 & 0 & 1 \\ 0 & 1 & 0 & \cdots & 0 & 0 & 1 \\ \end{pmatrix}

, nous déterminons le vecteur des transformations,

W_c = (\frac{\sharp \lbrace X ^6 = "A" \rbrace}{20}, \frac{\sharp \lbrace X ^6 = "B" \rbrace}{20}, \frac{\sharp \lbrace X ^7 = "A" \rbrace}{20}, \cdots, \frac{\sharp \lbrace X ^9 = "A" \rbrace}{20}, \frac{\sharp \lbrace X ^9 = "A" \rbrace}{20}, \frac{\sharp \lbrace X ^9 = "A" \rbrace}{20})

= (\frac{10}{20}, \frac{10}{20}, \frac{6}{20}, \cdots, \frac{7}{20}, \frac{6}{20}, \frac{7}{20})

, et enfin nous pouvons calculer la matrice standardiser des variables qualitatives,

\mathbf{B}_T ^D = \begin{pmatrix} \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{1}{\frac{6}{20}} - 1 & \cdots & \frac{1}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{0}{\frac{7}{20}} - 1 \\ \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{1}{\frac{6}{20}} - 1 & \cdots & \frac{1}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{0}{\frac{7}{20}} - 1 \\ \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{1}{\frac{6}{20}} - 1 & \cdots & \frac{1}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{0}{\frac{7}{20}} - 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{0}{\frac{6}{20}} - 1 & \cdots & \frac{0}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1 \\ \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{0}{\frac{6}{20}} - 1 & \cdots & \frac{0}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1 \\ \frac{0}{\frac{10}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1  & \frac{0}{\frac{6}{20}} - 1 & \cdots & \frac{0}{\frac{7}{20}} - 1 & \frac{0}{\frac{6}{20}} - 1 & \frac{1}{\frac{10}{20}} - 1 \\ \end{pmatrix}

= \begin{pmatrix} -1 & 1 & 2.333333 & \cdots & 1.857143 & -1 & -1 \\ -1 & 1 & 2.333333 & \cdots & 1.857143 & -1 & -1 \\ 2.333333 & -1 & -1 & \cdots & 1.857143 & -1 & -1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ -1 & 1 & -1 & \cdots & -1 & -1 & 1.857143 \\ -1 & 1 & -1 & \cdots & -1 & -1 & 1.857143 \\ -1 & 1 & -1 & \cdots & -1 & -1 & 1.857143 \\ \end{pmatrix}

Maintenant nous reformons notre matrice complète (\mathbf{A}, \mathbf{B}_T ^D) et nous lui appliquons la pondération en ligne \sqrt{W_l} = (\sqrt{\frac{1}{20}}, \cdots, \sqrt{\frac{1}{20}}) de taille 1 \times 20, nous obtenons,

(\mathbf{A}, \mathbf{B}_T ^D)_{\sqrt{W_c}} = \begin{pmatrix} -0.37268028 & 0.120075959 & -0.16902823 & \cdots & 0.4152698 & -0.2236068 & -0.2236068 \\ -0.32613175 & 0.183400823 & -0.09124687 & \cdots & 0.4152698 & -0.2236068 & -0.2236068 \\ -0.28976801 & -0.358132113 & -0.03732516 & \cdots & 0.4152698 & -0.2236068 & -0.2236068 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0.29786106 & 0.190319484 & -0.33048717 & \cdots & -0.2236068 & -0.2236068 & 0.4152698 \\ 0.32980659 & 0.104423793 & 0.29849215 & \cdots & -0.2236068 & -0.2236068 & 0.4152698 \\ 0.36549807 & 0.220740729 & 0.37839333 & \cdots & -0.2236068 & -0.2236068 & 0.4152698 \\ \end{pmatrix}

Nous appliquons une nouvelle pondération en colonne: \sqrt{W_c} = (\sqrt{1}, \sqrt{1}, \sqrt{1}, \sqrt{1}, \sqrt{1}, \sqrt{\frac{10}{20}}, \sqrt{\frac{10}{20}}, \sqrt{\frac{6}{20}}, \sqrt{\cdots, \frac{7}{20}}, \sqrt{\frac{6}{20}}, \sqrt{\frac{7}{20}}), de taille 1 \times 15 et dont les premiers éléments correspondent à des poids unitaires associés aux variables continues. Nous obtenons alors,

(\mathbf{A}, \mathbf{B})_{\sqrt{W_c}, \sqrt{W_l}} = (\mathbf{A}, \mathbf{B})_P = \begin{pmatrix} -0.37268028 & 0.120075959 & -0.16902823 & \cdots & 0.2456769 & -0.1224745 & -0.1322876 \\ -0.32613175 & 0.183400823 & -0.09124687 & \cdots & 0.2456769 & -0.1224745 & -0.1322876 \\ -0.28976801 & -0.358132113 & -0.03732516 & \cdots & 0.2456769 & -0.1224745 & -0.1322876 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0.29786106 & 0.190319484 & -0.33048717 & \cdots & -0.1322876 & -0.1224745 & 0.2456769 \\ 0.32980659 & 0.104423793 & 0.29849215 & \cdots & -0.1322876 & -0.1224745 & 0.2456769 \\ 0.36549807 & 0.220740729 & 0.37839333 & \cdots & -0.1322876 & -0.1224745 & 0.2456769 \\ \end{pmatrix}

Maintenant nous pouvons faire le produit croisé,

(\mathbf{A}, \mathbf{B})_P ^t \times (\mathbf{A}, \mathbf{B})_P  = \begin{pmatrix} 1 & 0.130828486 & 0.26407954 & \cdots & -0.66671271 & 0.0003475779 & 0.666390918 \\ 0.1308284856 & 1 & 0.10607370 & \cdots & -0.19052973 & 0.1994050956  & 0.005916488 \\ 0.2640795394 & 0.106073704 & 1 & \cdots & -0.08051511 & -0.3103825387 & 0.367873498 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ -0.6667127122 & -0.190529734 & -0.08051511 & \cdots & 0.65 & -0.3240370349 & -0.35 \\ 0.0003475779 & 0.199405096 & -0.31038254 & \cdots & -32403703 & 0.7 & -0.324037035 \\ 0.6663909176 & 0.005916488 & 0.36787350  & \cdots & -0.35 & -0.3240370349 & 0.65 \\ \end{pmatrix}

Nous pouvons ainsi extraire les vecteurs propres,

\mathbf{V}_{HS} = \begin{pmatrix} -0.46521443 & 0.06264258 & 0.05431917 & \cdots & 0 & 0 & 0 \\ -0.11428424 & -0.05132924 & -0.82637677 & \cdots &  0 & 0 & 0 \\ -0.18002968 & 0.18564043 & 0.01720829 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0.33623486 & 0.18304311 & 0.05532521 & \cdots & 0.2870376 & 0.1386134 & -0.4922145 \\ -0.02431857 & -0.42243290 & -0.18784324 & \cdots & 0.2657452 & 0.128331 & -0.4557021 \\ -0.31372024 & 0.20805376 & 0.11858383 & \cdots & 0.2870376 & 0.1386134 & -0.4922145 \\ \end{pmatrix}

, et les valeurs propres associées,

\lambda = (4.079046, 3.102799, 1.170648, 0.9732852, 0.7077521, 0.3976803, 0.2601413, 0.1877866, 0.08289118, 0.03216202, 0.0058086, 0, 0, 0, 0)

Afin d’obtenir nos composantes principales, nous devons encore recourir à quelques opérations. Dans un premier dernier temps, nous devons appliquer les pondérations,

\frac{1}{\sqrt{W_c}} = (\frac{1}{\sqrt{1}}, \frac{1}{\sqrt{1}}, \frac{1}{\sqrt{1}}, \frac{1}{\sqrt{1}}, \frac{1}{\sqrt{1}}, \frac{1}{\sqrt{\frac{10}{20}}}, \frac{1}{\sqrt{\frac{10}{20}}}, \frac{1}{\sqrt{\frac{6}{20}}}\frac{1}{\sqrt{\frac{7}{20}}}, \frac{1}{\sqrt{\frac{7}{20}}}, \frac{1}{\sqrt{\frac{7}{20}}}, \frac{1}{\sqrt{\frac{13}{20}}}, \frac{1}{\sqrt{\frac{7}{20}}}, \frac{1}{\sqrt{\frac{6}{20}}}, \frac{1}{\sqrt{\frac{7}{20}}})

, à \mathbf{V}_{HS},

(\mathbf{V}_{HS})_P = \begin{pmatrix} -0.46521443 & 0.06264258 & 0.05431917 & \cdots & 0 & 0 & 0 \\ -0.11428424 & -0.05132924 & -0.82637677 & \cdots & 0 & 0 & 0 \\ -0.18002968 & 0.18564043 & 0.01720829 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0.56834065 & 0.30939932 & 0.09351668 & \cdots & -0.4851821 & 0.232994 & -0.8319944 \\ -0.04439944 & -0.77125342 & -0.34295326 & \cdots & 0.4851821 & 0.2342994 & -0.8319944 \\ -0.53028398 & 0.35167504 & 0.20044326 & \cdots & 0.4851821 & 0.2342994 & -0.8319944 \\ \end{pmatrix}

Calculons maintenant le dernier vecteur des pondérations, basé sur la racine carré de \lambda,

\sqrt{\lambda} = (2.01966491, 1.76147627, 1.08196471, 0.98655217, 0.84128003, 0.63061902, 0.51004050, 0.43334352, 0.28790828, 0.17933773, 0.07621424)

, que nous appliquons dans un second dernier temps à la matrice (\mathbf{V}_{HS})_P, ce qui nous donne nos composantes factorielles,

CF_{HS} = \begin{pmatrix} -0.93957726 & 0.11034342 & 0.058771422 & \cdots & 0 & 0 & 0 \\ -0.09614505 & -0.3236920 & -0.421485623 & \cdots & 0 & 0 & 0 \\ -0.05183204 & 0.03329233 & 0.001311517 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots  \vdots & \vdots \\ & 0.47813364 & 0.19511310 & 0.047697294 & \cdots & 0.3059651 & 0.0195022 & -0.3605394 \\ -0.01278297 & -0.13831484 & -0.026137924 & \cdots & 0.08701146 & 0.001785695 & -1.680350 \\ -0.93408265 & 0.38049998 & 0.197747732 & \cdots & 0.5249499 & 0.2311485 & -0.6999402 \\ \end{pmatrix}

Si nous retenons les deux premières composantes factorielles, c’est:

100 \times (\frac{4.079046}{\sum \lambda} + \frac{3.102799}{\sum \lambda}) = 100 \times (0.3708223968 + 0.2820726045) = 65.2895 \%

, de l’information qui est restituée, ce qui nous permet d’obtenir un plan d’assez bonne qualité.

Ci-dessous la projection des variables,

add.png

En terme d’interprétation, nous constatons trois cas de figures:

– Les modalités X ^6 = "A", X ^7 = "B", X ^9 = "B" sont liées entre elles mais aussi avec la variable continue X ^4 mettant en évidence que plus un individu a une forte valeur pour cette caractéristique et plus il a de chance d’avoir répondu aux trois modalités évoquées.

– Les modalités X ^7 = "C", X ^8 = "B", X ^9 = "C" sont liées entre elles mais aussi avec les variables continues X ^3, X ^5 (même si la relation entre elles est faible vue leur distance avec le cercle de corrélation). Ainsi, plus un individu a de fortes valeurs pour ces deux informations et plus il a de chance d’avoir répondu aux trois modalités évoquées. A noter une relation à modérer avec la variable continue X ^1. Enfin, ce groupe est indépendant au premier groupe évoqué.

– Les modalités X ^6 = "B", X ^7 = "A", X ^9 = "A" sont liées entre elles mais ne présentent aucun lien avec les variables continues. La modalité X ^8 = "A", quand a elle, est anti-corrélé avec le second groupe évoqué.

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_prinqual_sect004.htm

Package et fonction R: pbil.univ-lyon1.fr/ADE-4/ade4-html/dudi.hillsmith.html

\bullet Bibliographie:

– Principal component analysis of taxonomic data with multi-state discrete characters de M. O. Hill et J. E. Smith

– La page web: https://pbil.univ-lyon1.fr/R/pdf/thema2B.pdf

Les méthodes de partitionnement

add

\bullet Présentation:

Les méthodes de partitionnement font parties des trois familles d’outils d’analyse non supervisée les plus répandues avec la classification ascendante hiérarchique (CAH) et les méthodes à estimation de densité.

Les méthodes de partitionnement sont particulièrement intéressantes dans la construction de clusters/groupes d’observations d’une matrice de données à P \geq 1 variable(s) continue(s), \mathbf{X} = (X ^1, \cdots, X ^P), à partir de la structure même des données sans apport informatif d’une variable auxiliaire. A noter que des travaux récents ont permis d’étendre cette famille d’outils au cas des variables qualitatives.

Les méthodes de partitionnement les plus connues sont les nuées dynamiques, la méthode des centres mobiles, des K-means, des K-medoids, des K-modes, des K-prototypes et des K-médianes. La littérature diverge assez souvent au sujet des nuées dynamiques, dans certains cas elles sont mentionnées comme l’algorithme commun à toutes les variantes de méthodes de partitionnement, dans d’autres elles sont assimilées aux K-means.

Les avantages de cet outil d’analyse sont, principalement, sa simplicité et sa rapidité d’exécution même sur un grand nombre d’observations. A contrario, nous pourrons lui reprocher une classification finale un peu trop approximative et le besoin de spécifier à priori un nombre de classes à obtenir en fin de procédure. Par conséquent, les méthodes de partitionnement sont généralement utilisées dans un objectif d’analyse préliminaire avant le déploiement d’une méthode plus robuste comme la classification hiérarchique ascendante ou encore les méthodes à estimation de densité.

\bullet L’Algorithme des nuées dynamiques:

Par abus de langage, nous partirons du postulat que l’algorithme commun des méthodes de partitionnement est celui des nuées dynamiques.

L’Algorithme des nuées dynamiques a été conçu par Hugo Steinhaus et Stuart Lloyd en 1957. Les nuées dynamiques reposent sur la convergence d’un algorithme itératif relativement trivial. Il nécessite de paramétrer le nombre de clusters K que nous souhaitons obtenir en fin de procédure. L’algorithme se décrit ainsi,

– Itération 1: Sélectionner K points aléatoirement au sein du jeu de données. Ensuite, séparer le jeu de données en K premiers groupes selon une méthode de partitionnement fixée. Nous obtenons ainsi une première partition des données.

– Jusqu’à convergence des K groupes: Recalculer les centres des K clusters selon le noyau sélectionné et appliquer le même repartitionnement qu’en étape une.

La convergence est obtenue lorsque les groupes ne varient plus.

\bullet Les variantes:

Globalement, toutes les méthodes de partitionnement se base sur l’algorithme des nuées dynamiques. Difficile de partir sur un langage commun étant donné que d’un pays à l’autre les appellations divergent. Parmi les variantes des nuées dynamiques nous avons,

– La méthode des centres mobiles: conçue par Edward Forgy en 1965,  c’est la distance avec le centre de gravité le plus proche de chaque partition qui est utilisée pour déterminer les K classes.

– La méthode des K-moyennes ou K-means: conçue par James McQueen en 1967, c’est la distance moyenne entre les points et le centre de gravité de chaque partition qui est utilisée pour déterminer les K classes.

– La méthode des K-medoids: conçue par Leonard Kaufman et Peter J. Rousseeuw en 1987, c’est la distance entre les observations sélectionnées via permutation et le désigné de la classe (medoïd) de chaque partition qui est utilisée pour déterminer les K classes.

– La méthode des K-médianes: décrite par Anil K. Jain et Richard C. Dubes en 1988, c’est la distance médiane entre les points et le centre de gravité de chaque partition qui est utilisée pour déterminer les K classes.

– La méthode des K-modes: décrite par Zhexue Huang en 1988 et qui est adaptée aux variables exclusivement qualitatives. C’est la similarité entre les observations de chaque partition qui est utilisée pour déterminer les K classes.

– La méthode des K-prototypes: conçue par Zhexue Huang en 1988, qui est adaptée aux mixtes de variables continues et qualitatives et combine l’approche par méthode des K-modes et des K-means. C’est la similarité et la distance moyenne au barycentre entre les observations de chaque partition qui est utilisée pour déterminer les K classes.

\bullet Exemple:

Soit le jeu de données ci-dessous,

add1

La projection des observations nous donne la carte des individus suivante,

add2.png

Appliquons l’algorithme des nuées dynamiques pour déterminer nos K = 2 clusters fixés au préalable.

Itération 1: La première itération consiste donc à sélectionner deux points au hasard et qui formeront nos deux premiers clusters. Nous tirons donc l’observation N°5 et l’observation N°2.

add

Cette première itération divise le plan en deux clusters, celui regroupant les observations N°3, 5, 6, 7, 8. Et celui regroupant les observations N° 1, 2, 4, 9, 10.

Itération 2: Nous calculons maintenant les centres de gravité relatifs aux deux groupes construit au cours de l’itération 1. Nous obtenons le barycentre du groupe vert, de coordonnées (6.3236,10.5678), et celui du groupe bleu, de coordonnées (9.0579,10.7572).

add

Cette seconde itération met à jour les deux clusters construit auparavant. Ainsi, pour le groupe vert, les observations 3, 6, 7, 8 sont conservées et l’observation 5 est envoyé dans le groupe bleu qui contient toujours les observations 1, 2, 4, 9, 10.

Itération 3: Nous calculons maintenant les centres de gravité relatifs aux deux groupes construits. Nous obtenons le barycentre du groupe vert, de coordonnées (3.3650,5.898), et celui du groupe bleu, de coordonnées (9.113,7.754).

add.png

La division de l’itération 3 confirme celle de l’itération 2, nous en déduisons que l’algorithme a convergé et que les deux clusters conçus sont,

– groupe 1: observations N°3, 6, 7, 8,

– groupe 2: observations N° 1, 2, 4, 5, 9, 10.

\bullet Application informatique:

– Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#cluster_toc.htm

– Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/hclust.html

\bullet Bibliographie:

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard

– Data mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani, Jerome Friedman

– Probabilités, analyse des données et statistique de Gilbert Saporta

Sur la division des corps matériels en parties de Hugo Steinhaus

– Least square quantization in PCM de Stuart Lloyd

– Cluster analysis of multivariate data: efficiency versus interpretability of classifications de  Edward Forgy

– Some Methods for classification and Analysis of Multivariate Observations de James McQueen

Clustering by means of Medoids, in Statistical Data Analysis Based on the L_1–Norm and Related Methods de Leonard Kaufman et Peter J. Rousseeuw

– Algorithms for clustering data de Anil K. Jain et Richard C. Dubes

– Clustering large data sets with mixed numeric and categorical values de Zhexue Huang

Les méthodes à estimation de densité

add.png

\bullet Présentation:

Les méthodes à estimation de densité font parties des trois familles d’outils d’analyse non supervisée les plus répandues avec la classification ascendante hiérarchique (CAH) et les méthodes de partitionnement.

Les méthodes à estimation de densité sont particulièrement efficaces dans la construction de clusters/groupes d’observations d’une matrice de données à P \geq 1 variable(s) continue(s), \mathbf{X} = (X ^1, \cdots, X ^P), à partir de la structure même des données sans apport informatif d’une variable auxiliaire.

Le clustering, ou l’action de créer des groupes d’observations, ou clusters dans un argo-anglicisme bien souvent utilisé pour simplifier le propos, s’exécute via un algorithme ascendant en fonction d’un paramètre préalablement fixé et en fonction de la méthode choisie.

Les trois méthodes à estimation de densité les plus connues sont les K-plus proches voisins, la méthode du noyau uniforme et la méthode hybride de Wong.

\bullet Les K-plus proches voisins:

Les K-plus proches voisins ont été élaborés par Anthony M. Wong et Tom Lane en 1983 et se base sur le principe de vote majoritaire en fonction de la proximité aux K voisins les plus proches pour déterminer à quel groupe appartient l’observation que nous cherchons à regrouper.

L’Algorithme nécessite au préalable de fixer le paramètre K dont l’influence sera majeure sur le résultat final. Il n’existe pas de méthodes pour optimiser ce paramètre.

Le choix du regroupement est donc calculé au travers de la formule:

Si d_{X_{i_1},X_{i_2}} \leq max(r_K (X_{i_1}), r_K (X_{i_2})) alors d_{X_{i_1},X_{i_2}} ^* = \frac{1}{2} (\frac{1}{f(X_{i_1})} + \frac{1}{f(X_{i_2})})

Avec,

r_K (X_i) la distance entre l’observation X_i et son Kème voisin le plus proche au sens de la distance euclidienne.

f(X_i) la fonction estimée de densité en X_i ou le rapport entre le nombre d’observations incluses dans la sphère, ou du cercle si P = 2, de rayon r_K (X_i) centrée en X_i et le volume, respectivement l’aire, de la sphère, respectivement du cercle.

Les K-plus proches voisins restent l’un des algorithmes les plus efficaces de la famille des méthodes à estimation de densité, principalement de part sa logique et sa facilité de déroulement. Une version pour analyse supervisée existe également et demeure tout aussi efficace.

add2

\bullet La méthode du noyau uniforme:

La méthode du noyau uniforme se base sur le principe de la création d’une sphère de rayon R (ou d’un cercle si P = 2) centrée en chacune des observations et du nombre d’observations qui y sont incluses pour déterminer à quel groupe appartient l’observation que nous cherchons à regrouper.

L’Algorithme nécessite au préalable de fixer le paramètre R dont l’influence sera majeure sur le résultat final. Il n’existe pas de méthodes pour optimiser ce paramètre.

Le choix du regroupement est donc calculé au travers de la formule:

Si d_{X_{i_1},X_{i_2}} \leq R alors d_{X_{i_1},X_{i_2}} ^* = \frac{1}{2} (\frac{1}{f(X_{i_1})} + \frac{1}{f(X_{i_2})})

, avec f(X_i) la fonction estimée de densité en X_i ou le rapport entre le nombre d’observations incluses dans la sphère, ou du cercle si P = 2, de rayon R centrée en X_i et le volume, respectivement l’aire, de la sphère, respectivement du cercle.

add2

\bullet La méthode hybride de Wong:

La méthode hybride de Wong, développée par Anthony M. Wong et C. Schaack en 1982, se base sur un algorithme en deux temps (d’où le terme d’hydride).

– Dans un premier temps il s’agit de lancer une partition selon la méthode des K-means de la famille des méthodes de partitionnement et de conserver la classification préliminaire ainsi obtenue. La méthode nécessite donc de fixer le paramètre K du nombre de clusters à obtenir.

– A partir de là, le second temps consiste à regrouper les différents clusters préliminaires. La recherche de deux clusters adjacents se fait au travers de la formule,

d ^2 (Cl_{k_1}, Cl_{k_2}) < d ^2 (Cl_{k_1}, Cl_{k_3}) + d ^2 (Cl_{k_2}, Cl_{k_3})

Dés lors la distance à minimiser pour déterminer les deux clusters adjacents à regrouper est,

d (Cl_{k_1}, Cl_{k_2}) ^* = \frac{(\sum_{i \in Cl_{k_1}} ||X_i - \overline{X_{Cl_{k_1}}}|| ^2 + \sum_{i \in Cl_{k_2}} ||X_i - \overline{X_{Cl_{k_2}}}|| ^2 + \frac{n_{k_1} + n_{k_2}}{4} d ^2(Cl_{k_1}, Cl_{k_2})) ^{\frac{\lambda}{2}}}{(n_{k_1} + n_{k_2}) ^{1 + \frac{\lambda}{2}}}

La méthode requiert également de fixer le paramètre de lissage \lambda.

L’Algorithme nécessite au préalable de fixer les paramètre K, \lambda dont les influences seront majeures sur le résultat final. Il n’existe pas de méthodes pour optimiser ces paramètres. La méthode hybride de Wong est considérée comme a éviter lorsque n < 100.

add2

\bullet Exemple:

Soit le jeu de données ci-dessous,

add1

La projection des observations nous donne la carte des individus suivante,

add2.png

Appliquons la méthode des 2-plus proches voisins afin d’établir notre classification.

Dans un premier temps nous calculons la matrice des distances sur \mathbf{X},

D_{\mathbf{X}}= \begin{pmatrix} & obs1 & obs2 & obs3 & obs4 & obs5 & obs6 & obs7 & obs8 & obs9 \\  obs2 & 1.730431 & & & & & & & & \\ obs3 & 7.306695 & 7.851494 & & & & & & & \\ obs4 & 0.991125 & 1.378891 & 8.214224 & & & & & & \\ obs5 & 2.229135 & 2.740852 & 5.190977 & 3.050761 & & & & & \\ obs6 & 10.860636 & 12.570271 & 10.627881 & 11.603027 & 10.847912 & & & & \\ obs7 & 8.594840 & 10.314997 & 9.309023 & 9.311560 & 8.746753 & 2.311944 & & & \\ obs8 & 6.403463 & 8.123649 & 9.287647 & 6.955009 & 7.149683 & 5.065955 & 2.830877 & & \\ obs9 & 5.463783 & 6.765099 & 11.353934 & 5.386607 & 7.317928 & 9.069772 & 6.941757 & 4.141981 & \\ obs10 & 4.224611 & 5.452226 & 10.553691 & 4.075979 & 6.198213 & 9.640028 & 7.401122 & 4.578377 & 1.327253 \\ \end{pmatrix}

Dans un second temps, calculons la distance moyenne avec les deux observations les plus proche de chaque observation prise individuellement,

D_{2-\mbox{proches voisins}}= \begin{pmatrix} Obs & 1er & 2nd & \overline{D} \\ obs1 & obs2 & obs4 & 1.360778 \\ obs2 & obs1 & obs4 & 1.554661 \\ obs3 & obs1 & obs5 & 6.248836 \\ obs4 & obs1 & obs2 & 1.360778 \\ obs5 & obs1 & obs2 & 2.484993 \\ obs6 & obs7 & obs8 & 3.688949 \\ obs7 & obs6 & obs8 & 3.688949 \\ obs8 & obs7 & obs9 & 3.486429 \\ obs9 & obs8 & obs10 & 2.734617 \\ obs10 & obs1 & obs9 & 2.775932 \\ \end{pmatrix}

– La distance minimale sur ces deux tables est tout d’abord celle entre les observations N°1 et N°4, formant ainsi le cluster Cl_{1,4}.

– Vient ensuite, la distance entre l’observation N°9 et l’observation N°10, formant le cluster Cl_{9,10}.

– Le troisième regroupement concerne l’observation N°2 avec le cluster Cl_{1,4} et formant désormais le cluster Cl_{1,2,4}.

– Le prochain regroupement est celui de l’observation N°6 et N°7, formant le cluster Cl_{6,7}. A noter que la distance entre les observations N°1 et N°5 est inférieure à celle citée ci-avant, mais étant donné que l’observation N°1 est déjà regroupée, il convient de prendre en compte sa distance pondérée par les 2 plus proches voisins (seconde table).

– Le cinquième regroupement met l’observation N°5 avec le cluster Cl_{1,2,4}, formant ainsi le cluster Cl_{1,2,4,5}.

– Maintenant c’est au tour de l’observation N°8 d’être regroupé eau cluster Cl_{6,7} même si sa distance avec ses deux plus proches voisins est plus petite, malheureusement elle implique deux clusters distincts Cl_{6,7} et Cl_{9,10} et reste plus proche du premier que du second. Nous obtenons ainsi le cluster Cl_{6,7,8}.

– Enfin, le dernier regroupement est celui de l’observation N°3 au cluster Cl_{1,2,4,5}, formant le cluster Cl_{1,2,3,4,5}.

Nous obtenons ainsi les trois groupes suivants,

add

\bullet Application informatique:

– Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#cluster_toc.htm

– Package et fonction R: http://finzi.psych.upenn.edu/R/library/spatstat/html/nndist.html

\bullet Bibliographie:

– Data Mining et statistique décisionnelle. L’Intelligence des données de Stéphane Tufféry

– A K-th Nearest Neighbor Clustering Procedure de Anthony M. Wong et Tom Lane

– A Hybrid Clustering Method for Identifying High-Density Clusters de Anthony M. Wong