Caribou!

Statistiques et Modélisation

Il est toujours assez étonnant de voir à quel point le statisticien n’est pas perçu à sa juste valeur par l’individu lambda de notre société… Qui n’a jamais entendu le fameux « ouarf les statistiques, on sait jamais d’où elles sortent ni comment ou par qui elles sont faites, il y a surement une arnaque quelque part! ». Mais pardonnons-leur leur incultisme à ce sujet!

Cette passionnante discipline qu’est celle des Statistiques (et de la Modélisation, mot qu’on oublie trop souvent de lui associer et qui fait la différence entre le simple statisticien et l’ingénieur en Statistiques) est en plein essor du fait du pouvoir de décision qu’elle détient au cœur d’un projet nécessitant une évaluation complète du risque pour passer d’une phase de développement à une autre. Ainsi cette discipline, en quantifiant l’information de manière rigoureuse et théoriquement solide, a permis le passage du choix binaire (oui/non) au choix flexible argumenté par la connaissance de la probabilité d’échec induite (nous pourrons désormais choisir « oui » si le risque associé est faible et « non » si ce dernier est trop important). Nous comprenons aisément pourquoi l’individu qui sait manier les Statistiques est alors celui détenant toutes les clés pour l’aboutissement d’un projet, en devenant même l’un des chaînons les plus importants. Mais la discipline ne se résume pas qu’à cela, notons aussi ses qualités indéniables et incontestables dans la construction de plans d’expérience et de sondage robustes assurant la fiabilité des résultats à venir sur des sous-échantillons représentatifs d’un contexte ou encore ses liens permanents avec des disciplines primordiales comme la théorie de la mesure, l’algèbre linéaire, l’analyse numérique, les séries de Fourier ou encore les algorithmes de la bio-informatique poussant toujours un peu plus l’élaboration théorique et l’apparition de nouvelles méthodologies. Les Statistiques sont également l’opportunité de découvrir un nombre extrêmement varié de domaines d’application tel que la biologie, l’écologie, l’agroalimentaire, le marketting, ect; car si d’entières différences sont flagrantes entre ces contextes multiples, le statisticien aura toujours la capacité de rajouter sa pierre à l’édifice du fait d’outils d’analyse adaptables à n’importe lequel d’entre eux.
Enfin, cette discipline ne serait pas ce qu’elle est sans le folklore qui l’entoure! Sa capacité à modéliser un phénomène et donc ses qualités prédictives ne peuvent qu’ajouter à son caractère séducteur, ainsi le statisticien n’est plus uniquement celui qui donne une photographie de ce qu’il observe mais il devient en quelque sorte un être surnaturel capable de prédire le futur, rien que cela! Alors finalement? Qui peut oser en dire autant?

L’objectif de ce site est de rassembler, de la manière la plus exhaustive qui soit, l’ensemble des tests statistiques, outils d’analyse exploratoire et méthodes d’analyse supervisée et non supervisée existantes. Il est également question de présenter la théorie se cachant derrière chaque outil ainsi que les algorithmes associés tout en présentant les syntaxes sous R et SAS pour les mettre en oeuvre. J’ai fait en sorte qu’il soit utile aussi bien au jeune diplômé qui recherche des informations pratiques et théoriques pour combler sa culture statheuse (oui, le ‘h’ est mis exprès par référence au mot ‘math’, il s’agit pour moi de la marque différenciant le simple statisticien utilisateur d’histogrammes sans savoir comment ça marche de celui capable de maîtriser la théorie et ainsi être plus performant dans l’usage des méthodes de Statistiques et Modélisation) qu’à la personne non affiliée à la discipline et qui se sent perdue lorsqu’elle se trouve face à son jeu de données et son désir compulsif de le faire parler.
L’ autre raison de la création de ce site est également dû aux multiples collaborations récentes que j’ai eu avec des étudiants en thèse venant de tous horizons ainsi que d’une réelle volonté de partager ma passion. Il est donc viser de créer un espace d’échanges dans le but d’élargir cette communication, proposant ainsi  mon expertise sur d’autres projets (aussi n’hésitez pas à me contacter même si je suis souvent over-pris!!!). Le site, qui n’en est encore qu’à un simple stade de fœtus je le reconnais (ouverture du site officiellement en: mi-2017, aussi mes excuses les plus sincères pour les pages incomplètes et les coquilles jusqu’à cette date-là…), sera ainsi en mue permanente tant l’ambition et les idées que j’ai pour lui sont légions, aussi ne soyez pas déboussolé si du jour au lendemain plus rien n’est à sa place!

Ci-dessous le planning du site:

add.png

Avant de conclure cette page de présentation, je tiens également à souligner le fait qu’il s’agit d’un site fait par passion pour ce domaine qu’est la Statistique, la Probabilité et la Modélisation et principalement motivé par les trois faits suivants: parfaire mes propres connaissances, centraliser les informations essentiels de chaque outil et offrir à l’intéressé le détail et l’application de ces différents outils (combien de fois ai-je trouvé de jolies formules théoriques servant à rien, « coucou les documents sur l’analyse discriminante de Fisher » pour ne citer qu’eux). Toutefois, comme tout site présent sur internet, il n’est pas exempt de tout reproche et ceci en dépit des efforts de rigueur effectués. Il sera toujours vivement conseillé de croiser les informations que vous pourrez en extraire avec d’autres sources afin de les valider (notamment celles indiquées dans la partie bibliographique des différents articles) et je dirais même: bien fou est celui qui se basera sur mon site pour en tirer des vérités absolues. Il est question ici de proposer un espace décrivant les différents outils de statistique et modélisation avec mes mots à moi et ma propre interprétation, rien de plus.

Enfin, car un bon site est un site capable de s’auto-hors-sujet avec panache! Sévissant sur cette formidable île qu’est Mayotte en tant qu’ingénieur ayant pour ambition de contribuer à son développement, c’est avec émotion que j’alimente ce site en promouvant le dernier, et de loin le plus féerique, département français (bon OK, à égalité avec la Corse…).

mayotte.jpg

Ce site est dédié à cette petite femme de gauche au coeur immense et cet ours polaire fada de l’OM, deux personnes exceptionnelles, qui me servent de parents.

Enfin je souhaite rendre hommage, dans un premier temps, aux trois personnes suivantes: Gilbert Saporta, Stéphane Tufféry et Ricco Rakotomalala, véritables sources d’inspiration lors de mes études et début de ma carrière professionnelle. Ne les connaissant pas personnellement, je continue de carresser délicatement le doux rêve de les rencontrer un jour (je m’imagine déjà complètement hystérique et perdu dans l’errection intellectuelle ce jour-là, je les plains d’avance).

Et dans un second temps, à tous ces célèbres mathématiciens (R. A. Fisher, K. Pearson, etc – liste à compléter rapidement) à qui nous devons notre discipline et qui ont réussi l’exploit de faire de l’or avec de l’oxygène. C’est grâce à leur travaux que notre civilisation ne cesse d’évoluer dans le bon sens (lorsque c’est le cas…) et rares sont ceux qui sont suffisamment intéressés et informés pour s’en rendre compte. Au nom de notre espèce, il convient de les remercier et de leur rendre hommage chaque jour pour leurs découvertes et ce qu’elles nous ont apporté et continuent de nous apporter quotidiennement. Merci.

Les méthodes de partitionnement (En cours de rédaction…)

La nuée dynamique

Origine: Hugo Steinhaus en 1957 puis étendue aux K-moyennes par J. B. McQueen en 1967 et finalement généralisé par Edwin Diday en 1971. Les nuées dynamiques appartiennent à la famille des méthodes de partitionnement des données. Ils sont une extension de la méthode des K-moyennes en proposant une plus grande plage de noyaux de partitionnement.

L’idée:

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. La première étape consiste à sélectionner C points aléatoirement au sein du jeu de données. Ensuite, nous séparons le jeu de données en C premiers groupes selon une méthode de partitionnement fixée. Nous recalculons alors les centres des C clusters selon le noyau sélectionné et appliquons le même repartitionnement qu’en étape une. La procédure continue jusqu’à ce que les centres ne bougent plus et que l’algorithme a convergé.

add1.png

Propriétés: Les nuées dynamiques manquent de précision et sont surtout utilisées pour un travail préliminaire. De plus, elles nécessitent le paramétrage du nombre de classes au préalable.

Les K-moyennes

Origine: Hugo Steinhaus en 1957 puis étendue aux K-moyennes par J. B. McQueen en 1967 et finalement généralisé par Edwin Diday en 1971. Les nuées dynamiques appartiennent à la famille des méthodes de partitionnement des données. Ils sont une extension de la méthode des K-moyennes en proposant une plus grande plage de noyaux de partitionnement.

L’idée:

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. La première étape consiste à sélectionner C points aléatoirement au sein du jeu de données. Ensuite, nous séparons le jeu de données en C premiers groupes selon une méthode de partitionnement fixée. Nous recalculons alors les centres des C clusters selon le noyau sélectionné et appliquons le même repartitionnement qu’en étape une. La procédure continue jusqu’à ce que les centres ne bougent plus et que l’algorithme a convergé.

add1.png

Propriétés: Les nuées dynamiques manquent de précision et sont surtout utilisées pour un travail préliminaire. De plus, elles nécessitent le paramétrage du nombre de classes au préalable.

Les méthodes à estimation de densité (En cours de rédaction…)

Les K-plus proches voisins

Origine: M. Anthony Wong et Tom Lane, 1983. Il s’agit de l’une des trois méthodes non paramétriques appartenant à la famille des estimations de densité.

L’idée:

La méthode des K-plus proches voisins 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. Le choix du regroupement est donc calculée au travers de la formule:

Si D_{Cl_{k_1}, Cl_{k_2}} ^* \leq max(r_k (X_{k_1}), r_k (X_{k_2})) alors D_{k_1,k_2} ^* = \frac{1}{2} (\frac{1}{f(X_{k_1})} + \frac{1}{f(X_{k_2}}})

add2

Propriétés: La méthode des K-plus proches voisins restent l’un des algorithmes les plus efficaces du clustering hiérarchique, principalement de part sa logique et sa facilité de déroulement. Une version pour analyse supervisée existe également et reste tout aussi efficace. Cependant, il reste très influencé par le paramètre K pas toujours évident à optimiser.

La méthode du noyau uniforme

Origine: Aucune trace…

L’idée:

La méthode du noyau uniforme se base sur le principe de la création d’une sphère de rayon R centrée en chacune des observations et du nombre d’observations voisines qui y sont incluses pour déterminer à quel groupe appartient l’observation que nous cherchons à regrouper. Le choix du regroupement est donc calculée au travers de la formule:

Si D_{Cl_{k_1}, Cl_{k_2}} ^* \leq R alors D_{k_1,k_2} ^* = \frac{1}{2} (\frac{1}{f(X_{k_1})} + \frac{1}{f(X_{k_2}}})

add2

Propriétés: Rest reste très influencé par le paramètre R pas toujours évident à optimiser.

La méthode hybride de Wong

Origine: M. Anthony Wong et C. Schaack, 1982. Il s’agit de l’une des trois méthodes non paramétriques appartenant à la famille des estimations de densité.

L’idée: La méthode hybride de Wong 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-moyennes et de conserver la classification préliminaire ainsi obtenue. A partir de là le second temps propose les regroupements des différents clusters préliminaires. Des clusters sont adjacents lorsque,

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

, la distance minimale implique que les deux clusters doivent être regroupés.

add2

Propriétés: La méthode hybride de Wong est a éviter lorsque n < 100.

La classification ascendante hiérarchique

add2

\bullet Présentation:

La classification ascendante hiérarchique, régulièrement mentionnée sous l’acronyme CAH, fait partie des trois familles d’outil d’analyse non supervisée les plus répandues avec les méthodes de partitionnement et les méthodes à noyaux de densité.

La classification ascendante hiérarchique est particulièrement efficace 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 définissant deux paramètres,

– une fonction de distance (également appelée fonction de proximité): distance euclidienne, de Manhattan, , de Minkowski, etc

– et une méthode spécifique de regroupement: distance maximale, méthode de Ward, méthode de Mc Quitty, etc

, à partir desquelles la typologie des observations de \mathbf{X} sera construite au travers de partitions emboîtées d’hétérogénéités croissantes d’où sa propriété hiérarchique. Il existe un très grand nombre de méthodes de regroupement et chacune d’elles dispose d’avantages et d’inconvénients.

A noter que la classification ascendante hiérarchique ne requiert aucune condition de validité, cependant l’application de méthodes de transformation des données peut avoir un effet plus ou moins accentué sur les clusters conçus en fin d’algorithme et sont même parfois recommandées avant application.

Enfin, le choix du nombre de clusters se fait en fonction de deux approches:

– soit nous définissons au préalable un nombre de groupes maximales à composer,

– soit en fonction de quatre outils décisionnels que nous présenterons: le R ^2, le pseudo-T ^2, le semi-partiel-R ^2 et le Cubic Clustering Criterion.

Le choix entre ces deux méthodes est purement subjectif et principalement lié au contexte car aucune n’est à considérer comme la meilleure.

\bullet L’algorithme de classification ascendante hiérarchique:

La CAH regroupe un ensemble de techniques construites sur le même principe: le regroupement successif des points par ordre de proximité croissante. L’Algorithme commun est le suivant:

Étape initiale: Considérer autant de classes qu’il y a d’observations et regrouper les deux observations qui sont les plus proches selon la fonction de distance choisie.

Étapes récursives suivantes: Tant qu’il reste une observation isolée (sans groupe) ou au moins deux groupes d’observations à fusionner, pour l’itération t,

—- Calculer la matrice de distance entre les différents groupes d’observations et les observations isolées, selon la fonction de distance choisie et lorsque nous sommes en présence de groupe d’observations, selon la méthode de regroupement choisie.

—- La distance minimale d_t au sein de la matrice de distance calculée à l’étape t indique le regroupement à faire. Il pourra alors s’agir aussi bien de la fusion de deux groupes d’observations comme d’un groupe observations et d’une observation isolée ou deux observations sans groupe. La distance d_t retenue sert enfin à construire l’arbre de classification ascendante hiérarchique.

—- Fin de l’itération t, nous passons à l’itération t+1 en reprenant les deux précédentes étapes.

Étape finale: A ce stade nous sommes dans la situation où il nous reste soit deux groupes d’observations, soit un groupe d’observations et une observation isolée, nous les regroupons pour conclure l’algorithme. Une fois toutes les observations classées, nous disposons de l’arbre de classification ascendante hiérarchique présentant la structure du clustering et permettant de déterminer à quel niveau couper l’arbre pour visualiser les différents classifications possibles en fonction des différentes classes emboîtées.

Bruno Falissard, dans son ouvrage (voir bibliographie), propose l’exemple suivant afin d’illustrer l’algorithme présenté ci-dessus. Soit les cinq observations que nous cherchons à regrouper selon une fonction de proximité triviale.

add1

Appliquons l’algorithme itératif de classification hiérarchique ascendante. Nous commençons donc en supposons que chaque observation est une classe, nous remarquons que la distance la plus petite entre les 4 \times 3 = 12 paires d’observations possibles est celle entre l’observation 2 et 3 (d_1).

add2

Nous les rassemblons donc et formons ainsi le premier regroupement de notre arbre de classification. Maintenant nous mettons à jour la matrice des distances en recalculant les différentes distances et en considérant le regroupement des observations 2,3 comme une entité unique (notée 6).

add3

Nous constatons que, cette fois-ci, la distance minimale est celle entre les observations 4 et 5 (d_2). Évidemment, cette distance est supérieure à d_1, d’où une paire de branche plus haute. Nous passons à l’itération suivante en mettant à jour l’information avec le cluster regroupant les observations 4, 5 (noté 7) à l’instar du cluster noté 6.

add4

La distance minimale est désormais celle entre le cluster noté 6 et l’observation isolée 1. Nous mettons à jour l’arbre de classification avec ce nouveau regroupement (noté 8) et procédons à l’étape finale.

add5.png

L’algorithme a fini de tourner et l’arbre construit permet de retrouver les différents groupes construits.

\bullet Les fonctions distances:

Comme expliquer en introduction, il faut définir au préalable deux paramètres afin d’utiliser une CAH. Nous présentons le premier: la fonction de distance à utiliser pour le calcul de la distance entre deux observations décrites au travers de P variables. Six principales fonctions existent, les voici:

– La distance euclidienne: d(X_{i_1}, X_{i_2}) = \sqrt{\sum_{j = 1} ^P (X_{i_1} ^j - X_{i_2} ^j) ^2}

– La distance maximale: d(X_{i_1}, X_{i_2}) = max_{j \in [1, P]} | X_{i_1} ^j - X_{i_2} ^j |

– La distance de Manhattan: d(X_{i_1}, X_{i_2}) = \sum_{j = 1} ^P | X_{i_1} ^j - X_{i_2} ^j |

– La distance de Canberra: d(X_{i_1}, X_{i_2}) = \sum_{j = 1} ^P \frac{| X_{i_1} ^j - X_{i_2} ^j |}{| X_{i_1} ^j | + | X_{i_2} ^j |}

– La distance de Minkowski (pour q, paramètre à fixer), qui est entre autre une généralisation de la distance euclidienne: d(X_{i_1}, X_{i_2}) = (\sum_{j = 1} ^P (X_{i_1} ^j - X_{i_2} ^j) ^q) ^{\frac{1}{q}}

\bullet Les méthodes de regroupement:

Nous abordons maintenant le second paramètre à déterminer: la méthode de regroupement. Plusieurs méthodes existent pour déterminer la proximité entre deux clusters ou un cluster et une observation isolée, en voici les principales.

La distance moyenne

Origine: Robert Reuven Sokal et Charles Duncan Michener, 1958. Nous pouvons la trouver sous l’appellation de « CAH du saut moyen » ou encore « average linkage ».

L’idée:

La distance moyenne se détermine en calculant toutes les distances possibles entre les observations respectives aux deux groupes étudiés et en déterminant la moyenne de ces distances. Les deux groupes dont la distance moyenne est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add

Propriétés: La distance moyenne tend à créer des clusters de faible variance identique. De plus, elle est moins sensible au bruit.

La distance médiane

Origine: John C. Gower, 1961. Nous pouvons la trouver sous l’appellation « median method ».

L’idée:

La distance des médianes se détermine en calculant toutes les distances possibles entre les observations respectives aux deux groupes étudiés et en déterminant la valeur médiane de ces distances. Les deux groupes dont la distance médiane est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add

Propriétés: La distance des médianes présentent l’avantage d’être particulièrement robuste aux outliers s’ils sont trop nombreux. Pour le reste, elle demeure l’une des moins populaires et très rarement utilisée lors d’une CAH.

La distance minimale

Origine: K. Florek, J. Lukaszewicz, J. Perkal, S. Zubrzycki, Louis L. McQuitty et P. H. A. Sneath, 1951-1957. Nous pouvons la trouver sous l’appellation de « CAH du saut minimum » ou encore « single linkage ».

L’idée:

La distance minimale consiste à calculer l’ensemble des distances possibles entre les observations respectives aux deux groupes considérés et à en retenir la plus petite. Les deux groupes dont la distance minimale est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La distance minimale est plus robuste que son homologue la distance maximale mais présente néanmoins une faiblesse (voir une force si c’est l’objectif visé). Soit la situation de deux classes éloignées et qui peuvent être reliées par un faible nombre d’observations proches les unes des autres et faisant la jonction entre les deux groupes. Dés lors, le regroupement basé sur la distance minimale ne considérera pas trois clusters (les deux groupes séparés et le lot d’observations entre eux) mais un unique cluster.

La distance maximale

Origine: Thorvald Julius Sørensen, 1948. Nous pouvons la trouver sous l’appellation de « CAH du saut maximum » ou critère du diamètre ou encore « complete linkage ».

L’idée:

La distance maximale consiste à calculer l’ensemble des distances possibles entre les observations respectives aux deux groupes considérés et à en retenir la plus grande.  Les deux groupes dont la distance maximale est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La distance maximale produit généralement des clusters de diamètres à peu près égaux, cependant elle est fortement influencée par la présence d’outliers. Logique par construction, étant donné le fait qu’elle se base sur la distance maximale entre deux clusters. A l’instar de la distance médiane, la distance maximale est très peu utilisée.

La distance des barycentres

Origine: Robert Reuven Sokal et Charles Duncan Michener, 1958. Nous pouvons la trouver sous l’appellation « centroid method ».

L’idée:

La distance des barycentres se détermine en calculant les barycentres des différents groupes à partir de la moyenne des coordonnées des observations respectives à ces clusters. Les deux groupes dont la distance des barycentres est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La distance des barycentres demeure plus robuste aux outliers que les autres méthodes, cependant elle est réputée comme étant moins performante que la méthode de Ward ou encore la distance moyenne.

La distance de Ward

Origine: Joe H. Ward, 1963. Nous pouvons la trouver sous l’appellation de méthode de Ward ou encore « Ward’s minimum-variance method ».

L’idée:

La distance de Ward se base sur la distance au carré entre les barycentres des deux groupes considérés et pondérés par leur effectif respectif. En notant Cl_{k_1}, Cl_{k_2} les deux clusters considérés et de barycentres respectifs G_{k_1}, G_{k_2}, la formule d’usage est alors,

D_{Cl_{k_1},Cl_{k_2}} = \frac{|| G_{k_1} - G_{k_2} || ^2}{\frac{1}{n_{k_1}} + \frac{1}{n_{k_2}}}

Les deux groupes dont la distance de Ward est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La distance de Ward est la méthode la plus souvent utilisée car elle a été conçue dans un objectif de maximisation de l’inertie inter-classe et de minimisation de l’inertie intra-classe. Elle aura tendance à produire des clusters de même taille et reste très sensible aux outliers. De plus, contrairement à la distance minimale, elle est très peu efficace face aux classes « étirées ».

La méthode du maximum de vraisemblance

Origine: Warren S. Sarle et M. J. Symons, 1983. Nous pouvons la trouver sous l’appellation de méthode EML (pour Expected Maximun-Likelihood).

L’idée:

La méthode procède au regroupement d’observations par maximisation de la vraisemblance à chaque itération sous les hypothèses de mélange de lois multi-normales, d’égalité des matrices de variances-covariances et de probabilités d’échantillonnage inégales. Le choix du regroupement est donc déterminé au travers de la formule:

D_{Cl_{k_1}, Cl_{k_2}} = n \cdot P \cdot log[1 + \frac{\mathbf{\Sigma_{k_1 \cup k_2}} - \mathbf{\Sigma_{k_1}} - \mathbf{\Sigma_{k_2}}}{\sum_{k = 1} ^K \mathbf{\Sigma_k}}] - \lambda [n_{k_1 \cup k_2} log (n_{k_1 \cup k_2}) - n_{k_1} log (n_{k_1}) - n_{k_2} log (n_{k_2})]

Où,

n désigne le nombre d’observations,

\mathbf{\Sigma_{k_1 \cup k_2}} la matrice de variances-covariances issues du regroupement des clusters Cl_{k_1}, Cl_{k_2},

n_{k_1 \cup k_2} l’effectif obtenu lorsque que nous regroupons les clusters Cl_{k_1}, Cl_{k_2},

– le paramètre \lambda, fixé à 2 par défaut, représente la pénalité à appliquer aux matrices de variances-covariances afin de s’affranchir d’éventuel biais dû aux outliers.

Les deux groupes dont la distance calculées à partir de la méthode EML est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La méthode EML est similaire à la méthode de Ward, réputée comme l’une des plus performantes méthodes de classification. Cependant, elle a tendance à créer des groupes de tailles disproportionnées.

La méthode du \beta-flexible

Origine: Godfrey N. Lance et William Thomas Williams, 1967.

L’idée:

La méthode du \beta-flexible se base directement sur la maximisation de l’inertie inter-classe et la minimisation de celle intra-classe au travers de la formule suivante:

D_{Cl_{k_1}, Cl_{k_2} \cup Cl_{k_3}} = (D_{Cl_{k_1}, Cl_{k_2}} + D_{Cl_{k_1}, Cl_{k_3}}) \frac{1 - \beta}{2} + \beta \cdot D_{Cl_{k_2}, Cl_{k_3}}

En général, le paramètre \beta est fixé à -0.25 par défaut. Les distances D_{Cl_{k_1}, Cl_{k_2}} et D_{Cl_{k_1}, Cl_{k_3}} correspondent à celles entre le groupe à fusionner aux deux autres groupes (maximisation de l’inertie inter-groupe) et la distance D_{Cl_{k_2}, Cl_{k_3}} est celle entre les deux groupes déjà construits (minimisation de l’inertie intra-groupe).

Les deux groupes dont la distance calculée à partir de la méthode du \beta-flexible est la plus petite sont ceux qui seront fusionnés pour l’itération en cours.

add1

Propriétés: La méthode du \beta-flexible est optimale dans une approche visant à obtenir les meilleurs critères de classification possibles. Néanmoins, il nécessite un choix sur le paramètre \beta dont l’incidence sur la classification finale reste majeure. Le problème est de taille puisqu’il n’y a pas vraiment d’approche pour déterminer le choix idéal de ce paramètre.

La distance de McQuitty

Origine: Louis L. Mc Quitty, 1966. 

L’idée: La distance de Mc Quitty consiste à calculer la moyenne des moyennes des moyennes des ….. des moyennes des distances entre les observations des deux groupes à fusionner.

Par exemple, soit le cluster regroupant les observations 1,4 au niveau t_1 de la classification puis l’observation 2 avec ce groupe au niveau t_2. Au niveau t_3, La distance entre une observation 5 et le groupe des observations 1, 2, 4 sera la distance entre 5 et 2 additionnée à la somme divisée par deux des distances entre 5 et 1 et 5 et 4, le tout divisé à son tour par deux.

La formule itérative et générale est alors,

D_{Cl_{k_1}, Cl_{k_2} \cup Cl_{k_3}} = \frac{D_{Cl_{k_1}, Cl_{k_2}} + D_{Cl_{k_1}, Cl_{k_3}}}{2}

add.png

Propriétés:  La distance de Mc Quitty se base sur une approche combinatoire des groupes plutôt que sur les observations isolées dans ces groupes.

\bullet Critères décisionnels pour le choix du nombre de clusters optimal:
 
Les algorithmes de classification ne fournissent pas directement le nombre de cluster optimisant la classification. Il faut pour cela s’appuyer sur le contexte ainsi que sur quelques indicateurs que nous présentons ici.
  
En définissant K_t le nombre de groupes à l’itération t du processus de classification, nous avons,
 
– Le R _t ^2, de formule,
 
R_t ^2 = \frac{\mbox{Inertie inter-classe}}{\mbox{Inertie totale}} = \frac{\sum_{i = 1} ^n || X_i - \overline{X} ||  ^2 - \sum_{k = 1} ^{K_t} \sum_{i_k \in Cl_k} || X_{i_k} - \overline{X_{k}} || ^2}{\sum_{i = 1} ^n || X_i - \overline{X} ||  ^2}
  
La formule induit que les observations qui seront encore non regroupées auront un poids nul sur le calcul du R_t ^2 puisque leur coordonnées sont confondues avec leur barycentre. Ainsi, pour le cas initial où une observation équivaut à un groupe, le R ^2 vaudra 1.
 
Le R_t ^2 mesure à quel proportion l’inertie inter-classe est maximal et (par contraposé, puisque égale à 1 - R_t ^2,à quel point l’inertie intra-classe est minimale). Plus le R_t ^2 est proche de 1 et plus cela implique que la classification est bonne. Généralement, le choix se porte sur la cassure de la courbe que nous pouvons tracer à partir des R_t ^2, t \in [1, t_{max}] (cf graphe extrait de l’ouvrage de Stéphane Tufféry et présent en bibliographie).
add
– Le semi-partiel-R_t ^2, de formule,
  
semi-partiel-R_t ^2 = R_{t - 1} ^2 - R_t ^2
  
A noter que lors de l’itération initiale (t = 1), R_1 = 1.
 
Le semi-partiel-R_t ^2 caractérise l’évolution de la chute du R ^2 au cours des différentes itérations de la procédure de classification et notamment la perte d’inertie observéeà chaque regroupement de deux clusters. En général, le choix se porte sur le nombre de classes dont le semi-partiel-R_{t-1} ^2 est faible et le semi-partiel-R_t ^2 est fort, ce qui est symbolisé par un « pic » suivi d’un « creux » lorsque nous traçons la courbe pour t \in [1, t_{max} (cf graphe extrait de l’ouvrage de Stéphane Tufféry et présent en bibliographie).
add
  
– Le pseudo-F_t, de formule,

pseudo-F_t = \frac{\frac{R_t ^2}{K_t - 1}}{\frac{1 - R_t ^2}{n - K_t}}

Le pseudo-F_t mesure le rapport entre l’inertie inter-classe et l’inertie intra-classe. En effet,

\frac{\frac{R_t ^2}{K_t - 1}}{\frac{1 - R_t ^2}{n - K_t}} = \frac{R_t ^2}{1 - R_t ^2} \cdot \frac{n - K_t}{K_t - 1} = \frac{\frac{\mbox{Inertie inter-classe}}{\mbox{Inertie totale}}}{\frac{\mbox{Inertie intra-classe}}{\mbox{Inertie totale}}} \cdot \delta = \frac{\mbox{Inertie inter-classe}}{\mbox{Inertie intra-classe}} \cdot \delta

Plus le rapport est grand et plus cela nous conforte dans le fait que l’inertie inter-classe est maximale et l’inertie intra-classe est minimale et donc que nous sommes en présence d’une bonne classification.

– Le pseudo-T_t ^2 pour le regroupement des clusters Cl_{k_1}, Cl_{k_2}, de formule,

pseudo-T_t ^2 = \frac{\sum_{i \in Cl_{k_1} \cup Cl_{k_2}} || X_i - \overline{X_{Cl_{Cl_1 \cup Cl_2}}} || ^2 - \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{\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}{n_{k_1} + n_{k_2} - 2}}

Le pseudo-T_t ^2 ne peut être calculé que pour le regroupement de deux groupes d’observations ou d’un groupe et d’une observation isolée. Il consiste à mesurer l’évolution du barycentre avant et après regroupement des deux classes. En effet, en notant M_k = \sum_{i \in Cl_k} || X_i - \overline{X_{Cl_k}} || ^2, nous avons,

\frac{M_{k_1 \cup k_2} - M_{k_1} - M_{k_2}}{\frac{M_{k_1} + M_{k_2}}{n_{k_1} + n_{k_2} - 2}} = (n_{k_1} + n_{k_2} - 1) \cdot (\frac{M_{k_1 \cup k_2}}{M_{k_1} + M_{k_2}} - 1)

Le pseudo-T_t ^2 vaut 0 si le rapport est égal à 1 soit que les barycentres de Cl_{Cl_{k_1} \cup Cl_{k_2}}, Cl_{k_1}, Cl_{k_2} sont confondus et donc que le regroupement ne permet pas de maximiser l’inertie inter-groupe.

A contratio, plus le pseudo-T_t ^2 est grand et plus la classification est de bonne qualité.

– Le Cubic Criterion Clustering, noté CCC_t:

CCC_t = K_t \cdot log [\frac{1 - E[R_t ^2]}{1 - R_t ^2}] \cdot \frac{\sqrt{\frac{n p ^*}{2}}}{(0.001 + E[R_t ^2]) ^{1.2}}

Pour expliquer les calculs de E[R_t ^2] et p ^*, il convient de détailler directement tout l’algorithme de calcul à appliquer et mélangeant les deux objets. D’abord, il faut savoir que E[R_t ^2] est, par définition, la valeur théorique du R_t ^2 dans le cas d’une distribution uniforme, et sans regroupement, des observations. Nous devons dans un premier temps, déterminer le vecteur\lambda des valeurs propres de la matrice de taille P \times P,

\mathbf{T} = \frac{\mathbf{X ^T} \cdot \mathbf{X}}{n - 1}

Nous devons ensuite calculer le vecteur des pondérations,

c = (\frac{\prod_{j = 1} ^P \sqrt{\lambda_j}}{K_t}) ^{\frac{1}{P}}

, que nous appliquons aux valeurs propres \lambda. Nous obtenons alors le vecteur,

u = \frac{\sqrt{\lambda}}{c}

Ce vecteur est la version initiale de celui à appliquer pour le calcul de E[R_t ^2]. En effet, nous devons appliquer une correction en fonction du nombre de valeurs du vecteur u supérieurs à 1. Cette correction se rapporte directement à P,

p ^* = min(\sharp \lbrace j \in [1, P]; u_j > 1 \rbrace, K - 1) 

Enfin, deux versions du calcul de E[R_t ^2] existent en fonction de la valeur de p ^*. Si p ^* \in ]0, P[, alors nous recalculons le vecteur u avec la nouvelle valeur de c,

c ^* = (\frac{\prod_{j = 1} ^{P ^*} \sqrt{\lambda_j}}{K}) ^{\frac{1}{p ^*}}

, et,

u ^* = \frac{\sqrt{\lambda}}{c ^*}

Nous avons alors,

E[R_t ^2] = 1 - [\frac{\sum_{j = 1} ^{p ^*} \frac{1}{n + u_j} + \sum_{j = p ^* + 1} ^P \frac{u_j ^2}{n + u_j}}{\sum_{j = 1} ^P u_j}] \cdot [\frac{(n - K_t) ^2}{n}] \cdot [1 + \frac{4}{n}]

Dans le cas où p ^* = 0p ^* = p, alors,

E[R_t ^2] = 1 - [\frac{\sum_{j = 1} ^P \frac{1}{n + u_j}}{\sum_{j = 1} ^P u_j}] \cdot [\frac{(n - K_t) ^2}{n}] \cdot [1 + \frac{4}{n}]

En ce qui concerne la lecture du CCC_t,

– Si CCC_t > 2: alors la classification peut être considérée comme bonne,

– Si 0 < CCC_t < 2: alors la classification peut être considérée sans pouvoir confirmer qu’elle soit de qualité,

– Si CCC_t < 0: c’est que nous sommes en présence d’outliers gênants (notamment si CCC_t < -30).

En général, un creux pour k < K classes suivi d’un pic pour k + 1 classes indique une bonne classification pour k+1 classes, surtout si nous avons une croissance ou une décroissance lente à partir de k+2 classes (cf graphe extrait de l’ouvrage de Stéphane Tufféry et présent en bibliographie).

add

Enfin, il convient de ne pas utiliser le Cubic Clustering Criterion pour la méthode de la distance minimale.

\bullet Exemple:

Soit le jeu de données ci-dessous,

add1

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

add2.png

Construisons la classification ascendante hiérarchique selon les deux paramètres suivants,

– la fonction de distance: distance euclidienne,

– la méthode de regroupement: la distance minimale.

Afin de trouver un compromis entre détail et rapidité des calculs, nous considérerons (., ., \cdots, .) comme le vecteur des éléments à sommer.

\bigstar L’itération initial,

Nous initialisons l’algorithme en calculant la matrice des distances (selon la fonction de distance sélectionnée) et obtenons,

\mathbf{D} = \begin{pmatrix} obs & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ 2 & 1.730431 & . & . & . & . & . & . & . & . \\ 3 & 7.306695 & 7.851494 & . & . & . & . & . & . & . \\ 4 & 0.991125 & 1.378891 & 8.214224 & . & . & . & . & . & . \\ 5 & 2.229135 & 2.70852 & 5.190977 & 3.050761 & . & . & . & . & . \\ 6 & 10.860636 & 12.570271 & 10.627881 & 11.603027 & 10.847912 & . & . & . & . \\ 7 & 8.594840 & 10.314997 & 9.309023 & 9.311560 & 8.746753 & 2.311944 & . & . & . \\ 8 & 6.403463 & 8.123649 & 9.287647 & 6.955009 & 7.149683 & 5.065955 & 2.830877 & . & . \\ 9 & 5.463783 & 6.765099 & 11.353934 & 5.386607 & 7.317928 & 9.069772 & 6.941757 & 4.141981 & . \\ 10 & 4.224611 & 5.452226 & 10.553691 & 4.075979 & 6.198213 & 9.640028 & 7.401122 & 4.578377 & 1.327253 \\ \end{pmatrix}

Afin de justifier le tableau \mathbf{D}, prenons le temps de montrer le calcul de la distance entre l’observation 1 et l’observation 2. Nous avons donc,
 
D(X_1, X_2) = \sqrt{(8.1472 - 9.0579) ^2 + (9.2858 - 10.7572) ^2}
 
= \sqrt{0.8293745 + 2.165018}
 
= \sqrt{2.994392}
 
= 1.730431
  
Revenons sur le processus itératif, la distance la plus petite au sein de \mathbf{D} est d(X_1, X_4) = 0.991125, correspondant au premier regroupement pour l’itération initiale de l’algorithme. Nous pouvons construire la première branche du dendogramme,

addbbbbbb.png

Calculons les différentes indicateurs pour l’itération initiale. Nous obtenons,

R_1 ^2 = \frac{\sum_{i = 1} ^9 (3.113284 ^2, 4.837503 ^2, \cdots, 5.482 ^2) - \sum_{k_1 = 1} ^9 (0.4955625 ^2, 0 ^2, 0 ^2, 0.4955625 ^2, 0 ^2, \cdots, 0 ^2)}{\sum_{i = 1} ^9 (3.113284 ^2, 4.837503 ^2, \cdots, 5.482 ^2)}

= \frac{245.9374 - 0.4911644}{245.9374}

= 0.9980029

semi-partiel-R_1 ^2 = 1 - R_1 ^2 = 1 - 0.9980029 = 0.0019971

pseudo-F_1 = \frac{\frac{0.9980029}{9 - 1}}{\frac{1 - 0.9980029}{10 - 9}} = \frac{0.1247504}{0.0019971} = 62.4654

Le pseudo-T_1 ^2 n’étant pas calculable pour cette première itération étant donné qu’il mesure l’évolution de l’inertie lorsque deux groupes de plusieurs observations sont fusionnés. De même pour le CCC_1, dont l’algorithme ne permet pas de retenir la moindre valeur propre.

\bigstar L’itération N° 2,

Nous passons désormais à la seconde itération dont la particularité sera de considérer les observations 1, 4 comme une unique entité et dont la distance sera déterminée en fonction de la méthode choisie. La figure ci-dessous présente ce premier regroupement (intitulé cluster A) et la méthode appliquée pour déterminer les distances à retenir entre les différentes observations et le cluster A.

add3

La distance la plus petite (1.327253, cf \mathbf{D}) est celle entre les observations 9 et 10. Nous procédons donc au regroupement et appellerons ce cluster le B. Nous pouvons donc mettre à jour le dendogramme.

addbbbbbb.png

Calculons les différentes indicateurs pour l’itération 2. Nous obtenons,

R_2 ^2 = \frac{254.9375 - \sum_{k_1 = 1} ^9 (0.4955625 ^2, 0 ^2, 0 ^2, 0.4955625 ^2, 0 ^2, \cdots, 0 ^2, 0.6636267^2,  0.6636267^2)}{254.9375}

= \frac{245.9374 - 1.371965}{245.9374}

= 0.9944215

semi-partiel-R_2 ^2 = R_1 ^2 - R_2 ^2 = 0.9980029 - 0.9944215 = 0.0035814

pseudo-F_2 = \frac{\frac{0.9944215}{8 - 1}}{\frac{1 - 0.9944215}{10 - 8}} = \frac{0.1420602}{0.00278925} = 50.93121

Le pseudo-T_2 ^2 et le CCC_2 n’étant à nouveau pas calculable pour cette seconde itération et pour les même raisons décrites en première.

\bigstar L’itération N° 3,

Passons ensuite à l’itération N° 3. La figure ci-dessous caractérise la procédure suivie pour déterminer le nouveau regroupement à effectuer en fonction des groupes A, B conçus.

add2

La distance minimale (1.378891, cf \mathbf{D}) est cette fois-ci celle entre l’observation 2 et le groupe A. Par conséquent nous l’insérons au groupe et pouvons mettre à jour le dendogramme,

addbbbbbb.png

, et nous pouvons calculer les nouveaux indicateurs pour cette itération.

R_3 ^2 = \frac{254.9375 - \sum_{k_1 = 1} ^9 (0.8200341 ^2, 0.9893395 ^2, 0 ^2, 0.5550719 ^2, 0 ^2, \cdots, 0 ^2, 0.6636267^2,  0.6636267^2)}{254.9375}

= \frac{245.9374 - 2.840154}{245.9374}

= 0.9884517

semi-partiel-R_3 ^2 = R_2 ^2 - R_3 ^2 = 0.9944215 - 0.9884517 = 0.0059698

pseudo-F_3 = \frac{\frac{0.9884517}{7 - 1}}{\frac{1 - 0.9884517}{10 - 7}} = \frac{0.1647419}{0.003849433} = 42.79649

– Bonne nouvelle, nous sommes ici dans le cas du regroupement de deux groupes (celui contenant l’observation 2 et le groupe A issu de la fusion des observations 1, 4). Nous avons donc,

pseudo-T_3 ^2 = \frac{(0.8200341 ^2 + 0.9893395 ^2 + 0.5550719 ^2) - (0.4955625 ^2 + 0.4955625 ^2) - 0 ^2}{\frac{ (0.4955625 ^2 + 0.4955625 ^2) + 0 ^2}{2 + 1 - 2}}

= \frac{1.468189}{0.4911644}

= 2.9892

– Enfin, le CCC_3 n’étant à nouveau pas calculable pour cette troisième itération et pour les même raisons décrites en première et seconde.

\bigstar etc etc,

Nous continuons cet enchainement de calcul jusqu’à l’avant-dernier regroupement. A ce stade nous avons la carte des individus suivantes avec l’application de la méthode visant à déterminer les distances à retenir entre les différentes observations des clusters C et D,

addbbb

La distance minimale (4.141981, cf \mathbf{D}) est celle entre le cluster C et le cluster D. Par conséquent, nous les regroupons pour cette avant-dernière itération et obtenons le dendogramme suivant,

addbbbbbb

Calculons désormais les différents indicateurs,

R_8 ^2

= \frac{254.9375 - \sum_{k_1 = 1} ^9 (3.298928 ^2, 5.019707 ^2, 0 ^2, 3.887283 ^2, 4.314452 ^2, 7.767051 ^2, 5.459712 ^2, 3.104740 ^2, 3.590493 ^2,  3.009369^2)}{254.9375}

= \frac{245.9374 - 191.5287}{245.9374}

= 0.2212299

semi-partiel-R_8 ^2 = R_7 ^2 - R_8 ^2 \approx 0.7504299 - 0.2212299 \approx 0.592

pseudo-F_8 = \frac{\frac{0.2212299}{2 - 1}}{\frac{1-0.2212299}{10 - 2}} = \frac{0.2212299}{0.09734626} = 2.272607

pseudo-T_8 ^2

= \frac{(3.298928 ^2 + 5.019707 ^2 + 3.887283 ^2 + 4.314452 ^2 + 7.767051 ^2 + 5.459712 ^2 + 3.104740 ^2 + 3.590493 ^2 +  3.009369^2) - (2.4495831 ^2 + 0.3422138 ^2 + 2.6248904 ^2) - (1.174442 ^2 + 2.566814 ^2 + 1.254978 ^2 + 3.301218 ^2 + 4.312358 ^2 + 3.054969 ^2)}{\frac{(2.4495831 ^2 + 0.3422138 ^2 + 2.6248904 ^2) + (1.174442 ^2 + 2.566814 ^2 + 1.254978 ^2 + 3.301218 ^2 + 4.312358 ^2 + 3.054969 ^2)}{3 + 6 - 2}}

= \frac{191.5287 - 13.00762 - 48.37013}{\frac{13.00762 + 48.37013}{7}}

= \frac{130.1509}{8.76825}

= 14.84343

– Enfin, le CCC_8 est calculable pour cette itération alors profitons-en!

Donc nous débutons par le calcul de,

\mathbf{T} = \begin{pmatrix} 8.1472 & 9.0579 & \cdots & 9.6489 \\ 9.2858 & 10.7572 &\cdots & 5.3371 \\ \end{pmatrix} \times \begin{pmatrix} 8.1472 & 9.2858 \\ 9.0579 & 10.7572 \\ \cdots & \cdots \\ 5.4688 & 3.4694 \\ 9.5751 & 4.0119 \\ 9.6489 & 5.3371 \\ \end{pmatrix} \times \frac{1}{10 - 1}

= \begin{pmatrix} 55.20539 & 5.85148 \\ 50.85148 & 67.13948 \\ \end{pmatrix}

Nous en déduisons le vecteur de valeurs propres suivants,

\lambda = (112.372813, 9.972056)

 Le calcul de,

c = (\frac{\prod_{j = 1} ^2 \sqrt{\lambda_j}}{2}) ^{\frac{1}{2}} = \sqrt{\frac{10.60064 \times 3.157856}{2}} = 4.09116

, nous permet de calculer le vecteur u,

u = (\frac{10.600604}{4.09116}, \frac{3.157856}{4.09116}) = (2.5910999, 0.7718731)

Maintenant, déterminons le paramètre,

p ^* = min(\sharp \lbrace j \in [1,2]; u_j > 1 \rbrace, 2 - 1) = (1,1) = 1 \in ]0, 2[

Ceci fait, mettons à jour le vecteur u. Nous avons donc,

c ^* = (\frac{\prod_{j = 1} ^1 \sqrt{\lambda_j}}{2}) ^{\frac{1}{1}} = \frac{\sqrt{112.372813}}{2} = 5.300302

, et donc,

u ^* = (\frac{\sqrt{112.372813}}{5.300302},\frac{\sqrt{9.972056}}{5.300302}) = (2,0.595788)

Nous avons ensuite,

E[R_8 ^2] = 1 - [\frac{\sum_{j = 1} ^1 \frac{1}{n + u_j} + \sum_{j = 2} ^2 \frac{u_j ^2}{n + u_j}}{\sum_{j = 1} ^2 u_j}] \cdot [\frac{(10 - 2) ^2}{10}] \cdot [1 + \frac{4}{10}]

= 1 - [\frac{\frac{1}{10 + 2} + \frac{0.5957882}{10 + 0.5957888}}{2 ^2 + 0.595788 ^2}] \times 6.5 \times 1.4

= 1 - (\frac{0.1168338}{4.354963}) \times 6.4 \times 1.4

= 1 - 0.02682773 \times 6.4 \times 1.4

= 1 - 0.3981331

= 0.6018668889

Maintenant que nous avons le p ^* et le E[R_8 ^2], nous pouvons déterminer le CCC_8 à cette itération,

CCC_8 = log[\frac{1 - 0.6018668889}{1 - 0.2212299}] \times \frac{\sqrt{\frac{10 \times 1}{2}}}{(0.001 + 0.6018668889) ^{1.2}}

= log(0.5112332) \times 4.104113

= -2.75357

\bigstar La dernière itération,

Nous concluons cette dernière étape en reliant les deux groupes restant, celui contenant l’observation 3 et celui contenant les observations 1, 2, 4, 5, 6, 7, 8, 9, 10. La distance entre les deux groupes selon la méthode choisie est celle entre l’observation 3 et l’observation 5 du groupe décrit ci-dessus (5.190977, cf \mathbf{D}).

Nous obtenons alors le dendogramme suivant,

addb

\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

– A Statistical Method for Evaluating Systematic Relationships de Robert  Reuven Sokal et Charles Duncan Michener

– A Method of Establishing Groups of Equal Amplitude in Plant Sociology Based on Similarity of Species Content and Its Application to Analyses of the Vegetation on Danish Commons de Thorvald Julius Sørensen

– Clustering Criteria and Multivariate Normal Mixtures de M. J. Symons

– Cubic Clustering Criterion de Warren S. Sarle

– A General Theory of Classificatory Sorting Strategies. I. Hierarchical Systems de Godfrey N. Lance et William Thomas Williams

– A Comparison of Some Methods of Cluster Analysis de John C. Gower

– Sur la Liaison et la Division des Points d’un Ensemble Fini de K. Florek, J. Lukaszewicz, J. Perkal et S. Zubrzycki

– Taksonomia Wroclawska de K. Florek, J. Lukaszewicz, J. Perkal et S. Zubrzycki

– Elementary Linkage Analysis for Isolating Orthogonal and Oblique Types and Typal Relevancies de Louis L. McQuitty

– The Application of Computers to Taxonomy de P. H. A. Sneath

– Hierarchical Grouping to Optimize an Objective Function de Joe H. Ward

– Similarity Analysis by Reciprocal Pairs for Discrete and Continuous Data de Louis L. McQuitty

La régression linéaire (En cours de rédaction…)

\bullet Présentation:

\bullet La régression linéaire:

– R2 de Nagelkerke:

Le R ^2 de Nagelkerke est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 de Nagelkerke  est,

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– R2 de Mc Fadden:

Le R ^2 de Mc Fadden est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 de Mc Fadden est,

R ^2 = 1 - \frac{ln \hat{L} (M_{Complet})}{ln \hat{L} (M_{Intercept})}

Où, \hat{L} estimateur de log-vraisemblance en fonction du modèle M considéré.

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– R2 de Efron:

Le R ^2 de Efron est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 de Efron est,

R ^2 = 1 - \frac{\sum_{i = 1} ^N (y_i - \hat{\pi_i}) ^2}{\sum_{i = 1} ^N (y_i - \overline{y}) ^2}

Où, \hat{\pi} probabilité de prédiction du modèle construit.

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– Le R2 de Cox-Snell:

Le R ^2 de Cox-Snell est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 de Cox-Snell est,

R ^2 = 1 - \lbrace \frac{L(M_{Intercept})}{L(M_{Complet})} \rbrace ^{\frac{2}{N}}

Où, L(M) est la probabilité conditionnel de la variable réponse relativement aux variables explicatives.

http://statisticalhorizons.com/r2logistic

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

https://books.google.fr/books?id=AyIYAAn4a2kC&pg=PA458&lpg=PA458&dq=le+R+2+de+cox-snell&source=bl&ots=0Fy5boRh_J&sig=IQ8RANtQY1pBkSYPSyLnCIqPYuY&hl=fr&sa=X&ved=0CD8Q6AEwA2oVChMI7JHK35juxwIVR1caCh23SQZi#v=onepage&q=le%20R%202%20de%20cox-snell&f=false

– Le Cp de Mallows:

Publié par Collin Lingwood Mallow, Le Cp de Mallows est une approche permettant de choisir le meilleur modèle construit par régression des moindres carrés.

En posant p variables et N observations, la formule du Cp de Mallows est:

C_p = \frac{SSE_p}{S ^2} - N + 2 \cdot P

Où,

SSE_p = \sum_{i = 1} ^N (Y_i - \hat{Y_i}) ^2

– Criterion for the Rejection of Doubtful Observations de Benjamin Peirce

– On Peirce’s criterion de Benjamin Peirce

– Le R2 de McKelvey & Zavoina:

Le R ^2 de McKelvey & Zavoina est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 de McKelvey & Zavoina est,

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– Le R2 count:

Le R ^2 count est un équivalent au R ^2 de la régression linéaire.

La formule du R ^2 count est,

http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– Le test RESET de Ramsay:

\bullet Annexe théorique:

\bullet Exemple:

\bullet Application informatique:

\bullet Bibliographie:

Les séries temporelles (En cours de rédaction…)

add.png

\bullet Présentation:

\bullet Les modèles:

\bullet Les tests:

– Le test de Box-Peirce

Présentation:

Le test de Box-Peirce porte également le nom de test du portemanteau. Il s’agit d’un estimateur de corrélation pour série temporelle.

Le test:

Soit \epsilon = (\epsilon_1, \cdots, \epsilon_n) une famille de variables aléatoires de loi à densité. La statistique du test de portemanteau est:

$\displaystyle \hat{\rho}_{\varepsilon}(h)=\frac{\frac{1}{n-h}\sum_{i=1}^{n-h}{\varepsilon}_i{\varepsilon}_{i+h}}{\frac{1}{n}\sum_{i=1}^{n}{\varepsilon}_i^2} $

La statistique de Box-Peirce, qui s’inspire de celle du portemanteau, est alors:

$\displaystyle Q_{BP}=n\sum_{h=1}^k\hat{\rho}_{\varepsilon}^2(h) $

La statistique de test de Ljung-Box suit une loi du \chi ^2 à k degrés de liberté.

– Le test de Ljung-Box

Présentation:

Le test de Ljung-Box porte également le nom de test du portemanteau. Il s’agit d’un estimateur de corrélation pour série temporelle.

Le test:

Soit \epsilon = (\epsilon_1, \cdots, \epsilon_n) une famille de variables aléatoires de loi à densité. La statistique du test de portemanteau est:

$\displaystyle \hat{\rho}_{\varepsilon}(h)=\frac{\frac{1}{n-h}\sum_{i=1}^{n-h}{\varepsilon}_i{\varepsilon}_{i+h}}{\frac{1}{n}\sum_{i=1}^{n}{\varepsilon}_i^2} $

La statistique de Ljung-Box, qui s’insipire de celle du portemanteau, est alors:

$\displaystyle Q_{LB}=n(n+2)\sum_{h=1}^k\frac{\hat{\rho}_{\varepsilon}^2(h)}{n-h}, $

La statistique de test de Ljung-Box suit une loi du \chi ^2 à k degrés de liberté.

– Le test de Durbin-Watson

Présentation:

Le test de Durbin-Watson permet de détecter la présence d’autocorrélation au sein des résidus ou erreurs de prédiction associés à une analyse par régression.

Le test:

Posons e_t le résidu associé à l’observation au temps t. La statistique de test de Durbin-Watson est:

d = {\sum_{t=2}^T (e_t - e_{t-1})^2 \over {\sum_{t=1}^T e_t^2}},

Avec T nombre d’observations.

Afin de déterminer si nous sommes en présence d’autocorrélation positive ou négative il faut comparer la statistique de test d aux valeurs critiques basses ou hautes.

– Le test de Chow

Présentation:

Le test de Chow permet de déterminer si deux séries de coefficients de deux régressions linéaires distinctes sont égaux. Il est notamment très prisé dans le cadre de séries temporelles pour étudier les points de rupture à un instant t précis.

Enfin, il est également utilisé afin d’évaluer l’influence des variables indépendantes sur les deux groupes ainsi construits.

Le test:

Soit les deux types de modèle suivants:

y_t = a_1 + b_1 \cdot x_{1,t} + c_1 \cdot x_{2,t} + \cdots + \epsilon

y_t = a_12 + b_2 \cdot x_{1,t} + c_2 \cdot x_{2,t} + \cdots + \epsilon

La statistique du test de Chow est alors,

F= \frac{\frac{S_c - (S_1 + S_2)}{k}}{\frac{S_1 + S_2}{N_1 + N_2 - 2 \cdot k}}

Avec S_c la somme des carrés des résidus estimés du modèle initial, et S_1, S_2 celles respectives aux modèles 1 et 2. Enfin, k le nombre total de paramètres à estimer. 

La statistique du test suit une loi de Fisher avec (k, N_1 + N_2 - 2k degrés de liberté. L’hypothèse nulle du test de Chow est: « a_1 = a_2, b_1 = b_2, c_1 = c_2, \dots« .

– Le test de Dickey-Fuller

Présentation:

Le test:

– Le test du bruit blanc

Présentation:

Le test:

– Le test de Zivot-Andrews

Présentation:

Le test:

– Le test KPSS

Présentation:

Le test:

– Le test de Schmidt-Phillips

Présentation:

Le test:

– Le test de Phillips-Perron

Présentation:

Le test:

\bullet Exemple:

\bullet Application informatique:

\bullet Bibliographie:

– Séries chronologiques de Jean-Jacques Droesbeke, Bernard Fichet, Philippe Tassi

http://www.math-info.univ-paris5.fr/~cabanal/Annales/Stat-Processus/Poly/Poly-0/node30.html

https://en.wikipedia.org/wiki/Durbin%E2%80%93Watson_statistic

Applied Regression Analysis in Econometrics de E. Doran Howard

Introduction to Econometrics de Christopher Dougherty

– Tests of Equality Between Setes of Coefficients in Two Linear Regression de Gregory C. Chow

 

Les tests pour la détection d’outliers

add.png

\bullet Introduction:

Les outliers ou observations aberrantes sont des valeurs d’une variable continue X douteuses au sens d’un écart fort avec le reste de la distribution de X. L’intérêt au problème des outliers débute en 1750 avec les travaux de Roger Joseph Boxcovitch, Pierre Simon de Laplace et Adrien-Marie Legendre. 

Deux approches sont à suggérer pour le traitement des outliers: les conserver en adaptant le choix de l’outil statistique ou alors leur suppression. La principale motivation à leur détection est qu’ils peuvent étirer la moyenne ou l’écart-type à tort apportant ainsi un biais au calcul des estimateurs d’intérêt.

La plupart du temps les outliers ont pour origine une erreur de saisie ou de mesure ce qui privilégierait leur suppression, néanmoins le manque de certitude sur leur origine fait qu’il n’y a pas vraiment de méthode universelle car un outlier peut également être une réelle valeur et simplement inattendue. C’est d’ailleurs sur ce point précis que la communauté scientifique s’oppose lors de leur traitement.

Toute une batterie de tests et de critères statistiques a été conçue afin de déterminer au sens statistique du terme si une observation est un outlier ou non, le plus populaire de tous restant le critère décisionnel de Tukey.

add1

\bullet Les différents tests:

– Le critère de Tukey

Qui et quand: John Wilder Tukey, 1977

Le critère:

L’Outil possède l’avantage de pouvoir déterminer un seuil à partir duquel les observations, en fonction de leur valeur, doivent être considérées comme aberrante relativement à la distribution de X. Permettant ainsi la détection d’un ou plusieurs outliers simultanément.

L’algorithme décisionnel du critère de Tukey est, \forall i \in [1,n],

si X_i < Q_1 - 1.5 IQR ou X_i > Q_3 - 1.5 IQR

, alors l’observation i est un outlier.

Nous noterons Q_1 le premier quartile (soit le quantile à 25\%), Q_3 le troisième quartile (soit celui à 75\%) et IQR l’interquartile. Leurs formules respectives sont,

Q_1 = X_{E[\frac{n+1}{4}]} + 0.75 (X_{E[\frac{n+1}{4}] + 1} - X_{E[\frac{n+1}{4}]}),

Q_3 = X_{3 E[\frac{n+1}{4}]} + 0.75 (X_{3 E[\frac{n+1}{4}] + 1} - X_{3 E[\frac{n+1}{4}]}),

IQR = Q_3 - Q_1.

E[.] désigne la partie entière de ..

– Le test Q de Dixon

Qui et quand: Wilfrid Joseph Dixon, 1951.

Le test:

Le test Q de Dixon se base principalement sur la distance entre X_i, l’observation suspectée d’être aberrante, et celle qui la précède afin de déterminer s’il s’agit d’un outlier ou non.

Supposons X trié, la formule du test Q de Dixon est,

– si X_i \geq mediane(X) alors Q = \frac{X_{i + 1} - X_i}{X_n - X_1},

– si X_i \leq mediane(X) alors Q = \frac{X_i - X_{i - 1}}{X_n - X_1}.

Nous reportons ensuite Q à la table de Dixon ci-dessous selon un niveau de confiance \alpha fixé,

add

Si Q > Q_{\alpha \%}, alors l’observation i peut être considérée comme un outlier.

– Le critère de Chauvenet

Qui et quand: William Chauvenet, 1863.

Le critère:

Le critère de Chauvenet se base sur une approche stochastique pour la détermination d’outliers et considère les observations les unes après les autres. La procédure pour départager tous les outliers consiste à supprimer les observations en partant de la valeur maximale ou minimale et de relancer le test à chaque nouvelle observation à vérifier.

Pour une observation X_i suspectée d’être un outlier, le critère de Chauvenet est,

si n \cdot erfc (\frac{|X_i - \overline{X}|}{\sigma_X}) < 0.5 alors X_i est un outlier

, où que erfc est la fonction d’erreur complémentaire.

– Le test de Grubbs

Qui et quand: Frank E. Grubbs, 1969.

Le test:

Le test de Grubbs, également connu sous le nom de test du maximum des résidus normalisés, recherche à déterminer si la valeur maximale ou minimale de X est un outlier. Lorsque nous suspectons un lot d’observations, il faudra se tourner vers sa généralisation: le test de Tietjen-Moore.

La formule du test de Grubbs est,

G_{min} = \frac{\overline{X} - min_{i \in 1, \cdots, n} (X)}{s} et G_{max} = \frac{max_{i \in 1, \cdots, n} (X) - \overline{X}}{s}

Le critère de décision, selon un seuil de confiance \alpha fixé et permettant de déterminer si la valeur maximale et/ou la valeur minimale sont des outliers est,

Si G_{min}, G_{max} > \frac{n - 1}{\sqrt{n}} \sqrt{\frac{t_{\frac{\alpha}{n},n-2} ^2}{n - 2 + t_{\frac{\alpha}{n},n - 2} ^2}} alors ce sont des outliers

, t_{\delta,df} représentant la fonction quantile associée à la distribution de Student pour le seuil de confiance \delta et pour df degrés de liberté.

– Le test de Tietjen-Moore

Qui et quand: Gary Tietjen et Roger Moore, 1972

Le test:

Le test de Tietjen-Moore est une généralisation de celui de Grubbs permettant de détecter un lot d’outliers et plus uniquement si la valeur maximale et minimale en sont. A noter que les valeurs testées doivent se succéder au niveau de la distribution. De plus, le test de Tietjen-Moore reste très dépendant du paramètre I à fixer et représentant le nombre d’observations à tester, ce qui peut nuire à la qualité de la détection s’il est mal paramétré.

Trois version du test de Tietjen-Moore existe en fonction de ce que nous cherchons à déterminer (les données seront supposées triées),

– Pour la ou les valeurs minimale(s) de la première à la I-ème:

L_{I,min} = \frac{\sum_{i = 1} ^{n-I} (X_i - \overline{X_{-[1,I]}}) ^2}{\sum_{i = 1} ^n (X_i - \overline{X}) ^2}

– Pour la ou les valeurs maximale(s) de la I-ème à la n-ème:

L_{I,max} = \frac{\sum_{i = I + 1} ^n (X_i - \overline{X_{-[I,n]}}) ^2}{\sum_{i = 1} ^n (X_i - \overline{X}) ^2}

, où X_{-[.,.]} représente le vecteur X privé des éléments inclus dans l’intervalle [.,.],

– Pour la ou les valeurs minimale(s) et maximale(s) (approche symétrique):

E_I = \frac{\sum_{i = 1} ^{n - I} (z_i - \overline{z_{-[I,n]}}) ^2}{\sum_{i = 1} ^n (z_i - \overline{z_i}) ^2}

, où z est le vecteur X ordonné en fonction du score résiduel \forall i \in [1,n], |X_i - \overline{X}|.

Afin de déterminer la valeur seuil à partir de laquelle nous pouvons tester notre hypothèse il faut simuler selon une loi normale centrée-réduite un nombre L de vecteurs de taille identique n à celle de X et appliquer le test de Tietjen-Moore pour chacun. La valeur seuil correspond au quantile à 5\% de la distribution des tests réalisés.

Si la valeur L_{I,min}, L_{I,max}, E_I est inférieur à la valeur seuil obtenu alors nous pouvons conclure que les observations considérées sont des outliers.

– Le test généralisé de la déviation extrême de Student

Qui et quand: Bernard Rosner, 1983

Le test:

Le test généralisé ESD (Extrem Studentized Deviation) permet la détection simultanée d’un ou plusieurs outliers et ne nécessite aucun paramétrage.

La première étape pour pouvoir calculer le test généralisé ESD est de ranger par ordre croissant X et déterminer le vecteur Z tel que,

Z_1 = \frac{max_{i \in [1, n]} (| X_i - \overline{X} |)}{\sigma_X}

\forall i \in [2,n],

— Suppression de l’observation qui maximise max_I(|X_i - \overline{X_I}|) et mise à jour de l’ensemble des observations à considérer I et dépourvu d’elle.

— Calcul de Z_i = \frac{max_{i \in I} (| X_i - \overline{X_I} |)}{\sigma_{X_I}}

, où X_I représente le vecteur X dépourvu des observations dans l’intervalle I indiqué.

En seconde étape, il faut pour chaque observation du vecteur Z, calculer la statistique de test \lambda,

\forall i \in [1,n], \lambda_i = \frac{(n - i) t_{1 - \frac{\alpha}{2 (n - i + 1)},n - i - 1}}{\sqrt{(n - i - 1 + t_{1 - \frac{\alpha}{2 (n - i + 1)},n - i - 1} ^2) (n - i + 1)}}

, où t_{\delta,df} représente la fonction quantile associée à la distribution de Student pour le seuil de confiance \delta et pour df degrés de liberté pour un seuil \alpha fixé.

Si Z_i > \lambda_i alors l’observation i peut être considérée comme un outlier.

– Le test \tau modifié de Thompson

Qui et quand: David J. Thompson, 1950.

Le test:

Le test \tau modifié de Thompson repose sur un algorithme permettant de vérifier directement quelles observations au sein de X peuvent être considérées comme des outliers.

Il faut dans un premier temps calculer,

\forall i \in [1,n], \delta_i = | X_i - \overline{X}|

Et pour l’observation maximisant \delta_i, pour un seuil de confiance \alpha fixé, la comparer à la valeur de référence,

\tau = \frac{t_{\frac{\alpha}{2},n-2} (n - 1)}{\sqrt{n (n - 2 + t_{\frac{\alpha}{2},n-2} ^2)}}

, où t_{\beta,df} représente la fonction quantile de Student pour le seuil de confiance \beta et pour df degrés de liberté.

Si \delta_i > \tau \sigma_X alors l’observation peut être considérée comme un outlier. 

La procédure doit être réitérée jusqu’à ce que \delta_i \leq \tau \sigma_X en supprimant à chaque nouvelle itération l’observation détectée comme un outlier et en mettant à jour la moyenne et l’écart-type.

– Le critère de Peirce

Présentation: Benjamin Peirce, 1877

Le test:

Le critère de Peirce présente de nombreux avantages: plus rigoureux que la majorité des autres méthodes de détection, il permet de déterminer directement l’ensemble des outliers au sein de X et ne repose pas sur un critère arbitraire pour le rejet d’une observation. De plus, il reste très simple à mettre en oeuvre.

L’algorithme de Peirce est le suivant,

1- Déterminer \overline{X}, \sigma_X,

2- A l’aide de la table de Peirce, calculer R pour un outlier suspecté sur n observations de X,

3- Calculer | X_i - \overline{X} |_{max} = \sigma_X R,

4- Rechercher les observations tel que \forall i \in [1,n], | X_i - \overline{X} | > |X_i - \overline{X} |_{max}, celles répondant à ce critère pourront être classées comme outlier,

5- Si au moins une observation répond à ce critère alors nous la supprimons de X et recalculons | X_i - \overline{X} |_{max} pour les même valeurs de n, \overline{X}, \sigma_X, seul R change et doit être relevé sur la table de Peirce pour deux outliers suspectés sur n observations. L’Opération est à renouveler à chaque fois que nous détectons un nouvel outlier.

6- Si plus aucune observation ne répond au critère, alors il faut mettre à jour n, \overline{X}, \sigma_X et reprendre depuis l’étape une à 5 tant qu’à l’étape 4 au moins une observation est détectée comme outlier.

Ci-dessous, la table de Peirce pour la détermination de R,

add

\bullet Exemple:

Soit l’échantillon suivant,

add2

Le boxplot associé à X est,

add3.png

– Le critère de Tukey

Nous cherchons donc à calculer le seuil, selon le critère de Tukey, à partir duquel une observation peut être considérée comme un outlier.

Nous avons donc,

Q_1 = 3.1547 + 0.75 \times (0.35896 - 3.1457)

= 3.1547 + 0.75 \times 0.4349

= 3.1547 + 0.326175

= 3.480875

Q_3 = 9.3500 + 0.25 \times (9.4578 - 9.3500)

= 9.3500 + 0.25 \times 0.1078

= 9.3500 + 0.02695

= 9.37695

IQR = 9.37695 - 3.480875 = 5.896075

Par conséquent,

S_1 = Q_1 - 1.5 \times IQR

= 3.480875 - 1.5 \times 5.896075

= 3.480875 - 8.44112

= - 5.363237

S_2 = Q_3 + 1.5 \times IQR

= 9.37695 + 1.5 \times 5.896075

= 9.37695 + 8.844112

= 18.22106

Nous n’avons pas d’observation pour laquelle sa valeur est inférieure à S_1, cependant nous avons une observation (X_{12}) dont la valeur est supérieure à S_2. Nous en déduisons que X_{12} est un outlier.

– Le test Q de Dixon

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Si nous regardons la table de Dixon, la valeur seuil pour n = 20 et au risque de 5\% est: 0.450. Nous avons donc,

Q = \frac{2.1475 - 1.9593}{19.1245 - 1.9593} = \frac{0.1882}{17.1652} = 0.1882 < 0.450,

Q = \frac{19.1245 - 9.8308}{19.1245 - 1.9593} = \frac{9.2937}{17.1652} = 0.5414268 > 0.450.

Nous en déduisons, qu’au risque de 5\%, la valeur minimale (1.9593) n’est pas un outlier contrairement à la valeur maximale (19.1245) qui en est un.

– Le critère de Chauvenet

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Dans un premier temps calculons la moyenne et l’écart-type,

\mu = 6.605515,

\sigma = 4.095388.

Nous pouvons maintenant calculer les éléments à soumettre au critère de Chauvenet. Nous avons alors,

\frac{1.9593 - 6.605515}{4.095388} = \frac{-4.646215}{4.095388} = -1.134499

\Rightarrow 20 \times erfc(|-1.134499|) = 20 \times 0.1086206525 = 2.1724130491

\Rightarrow 2.1724130491 > 0.5, nous en concluons que la valeur minimale n’est pas un outlier.

\frac{19.1245 - 6.605515}{4.095388} = \frac{12.51899}{4.095388} = 3.05685

\Rightarrow 20 \times erfc(|3.05685|) = 20 \times 0.0000153895 = 0.0003077894

\Rightarrow 0.0003077894 < 0.5, nous en concluons que la valeur maximale est un outlier.

– Le test de Grubbs

Nous cherchons à voir si la plus petite et la plus grande valeur de X sont des outliers. Dans un premier temps calculons la moyenne et l’écart-type,

\mu = 6.605515,

\sigma = 4.095388.

Nous avons donc,

G_{min} = \frac{6.605515 - 1.9593}{4.095388} = \frac{4.646215}{4.095388} = 1.134499

G_{max} = \frac{19.1245 - 6.605515}{4.095388} = \frac{12.51899}{4.095388} = 3.056851

Nous devons comparer ces deux valeurs au seuil,

\frac{20 - 1}{\sqrt{20}} \times \sqrt{\frac{t_{\frac{0.05}{20},20-2} ^2}{20 - 2 + t_{\frac{0.05}{20},20-2} ^2}} = \frac{19}{4.472136} \times \sqrt{\frac{t_{0.00125,18} ^2}{18 + t_{0.00125,18} ^2}}

= 4.248529 \times \sqrt{\frac{(-3.196574) ^2}{18 + (-3.196574) ^2}}

= 4.248529 \times \sqrt{\frac{10.21809}{28.21809}}

= 4.248529 \times 0.6017568

= 2.556381

Nous constatons que G_{min} < 2.556381 et G_{max} > 2.556381 et en concluons que seul la valeur maximale de X peut être considérée comme un outlier.

– Le test de Tietjen-Moore

Nous cherchons à déterminer si la valeur maximale et minimale sont des outliers, nous fixons donc I = 2. Nous appliquons donc la troisième version du test de Tietjen-Moore.

Dans un premier temps, nous avons les valeurs résiduelles,

\forall i \in[1,20], |X_i - \overline{X}| = (4.646215, 4.058315, \cdots, 1.482015, 1.340115)

Les valeurs maximales sont alors 4.646215, 12.518985, qui correspondante au minimum de X et à son maximum. Le vecteur z ordonné selon les valeurs résiduelles et contenant les valeurs de X est alors,

z = (7.5456, 5.2654, 5.2575, 5.1235, 8.1457, 8.9854, 4.1493, 4.1254

, 9.3500, 9.4578, 9.5965, 9.6160, 3.5896, 9.8308, 3.1547, 3.1386, 2.5472, 2.1475, 1.9593, 19.1245)

Nous avons donc,

E_2 = \frac{\sum_{i = 1} ^{20 - 2} (z_i - \overline{z_I}) ^2}{\sum_{i = 1} ^n (z_i - \overline{z}) ^2} = \frac{\sum_{i = 1} ^{18} (z_i - 6.168139) ^2}{\sum_{i = 1} ^n (z_i - 6.605515) ^2} = \frac{136.9162}{318.6719} = 0.4296463

Maintenant, simulons le seuil. Pour L = 100000 simulations, nous obtenons un quantile à 5\% égal à 0.4150498. Nous constatons que E_2 > 0.4150498, par conséquent nous ne pouvons conclure que la valeur minimale et maximale sont des outliers.

– Le test généralisé de la déviation extrême de Student

Nous cherchons à déterminer au sein de X si les trois observations aux valeurs les plus grandes et les plus petites peuvent être considérées comme des outliers.

Dans un premier temps déterminons les différents éléments du vecteur Z,

\overline{X} = 6.605515, sigma_X = 4.095388,

Z_1 = \frac{max_{i \in [1,20]} (| X_i -  6.605515|)}{4.095388} = \frac{12.51899}{4.095388} = 3.05685

L’observation pour laquelle max_{i \in [1,20]} (|X_i - \overline{X}|) est atteint est la valeur maximale de X, 19.1245. Elle sera donc supprimée pour le calcul de Z_2.

\overline{X_{-20}} = 5.946621, \sigma_{X_{-20}} = 5.946621,

Z_2 = \frac{max_{i \in [1,19]} (| X_i -  5.946621|)}{2.92212} = \frac{3.987321}{2.92212} = 1.36453

L’observation pour laquelle max_{i \in [1,19]} (|X_i - \overline{X_{-20}}|) est atteint est la valeur minimal de X_{-20}, 1.9593. Elle sera donc supprimé pour le calcul de Z_3.

\overline{X_{-\lbrace 1,20 \rbrace }} = 6.168139, \sigma_{X_{-\lbrace 1,20 \rbrace }} = 2.837938,

Z_3 = \frac{(max_{i \in [2,19]} (| X_i -  6.168139|)}{2.837938} = \frac{4.020639}{2.837938} = 1.416747

Maintenant, calculons la valeur seuil qui nous permettra de déterminer si nos observations peuvent être considérées comme des outliers,

\lambda_1 = \frac{(20 - 1) \times t_{(1 - \frac{0.05}{2 (20 - 1 + 1)},20-1-1)}}{\sqrt{(20 - 1 - 1 + t_{(1 - \frac{0.05}{2 (20 - 1 + 1)},20-1-1)} ^2) \times (20 - 1 + 1)}}

= \frac{19 \times t_{(0.99875,18)}}{\sqrt{(18 + t_{(0.99875,18)}) \times 20}}

= \frac{19 \times 3.510104}{\sqrt{(18 + 3.510104) \times 20}}

= \frac{66.69198}{24.62553}

= 2.708245

\lambda_2 = \frac{(20 - 2) \times t_{(1 - \frac{0.05}{2 (20 - 2 + 1)},20-2-1)}}{\sqrt{(20 - 2 - 1 + t_{(1 - \frac{0.05}{2 (20 - 2 + 1)},20-2-1)} ^2) \times (20 - 2 + 1)}}

= \frac{18 \times t_{(0.9986842,17)}}{\sqrt{(17 + t_{(0.9986842,17)}) \times 19}}

= \frac{18 \times 3.5193}{\sqrt{(17 + 3.5193) \times 19}}

= \frac{63.3474}{23.62888}

= 2.680931

\lambda_3 = \frac{(20 - 3) \times t_{(1 - \frac{0.05}{2 (20 - 3 + 1)},20-3-1)}}{\sqrt{(20 - 3 - 1 + t_{(1 - \frac{0.05}{2 (20 - 3 + 1)},20-3-1)} ^2) \times (20 - 3 + 1)}}

= \frac{17 \times t_{(0.9986111,17)}}{\sqrt{(17 + t_{(0.9986111,17)}) \times 18}}

= \frac{17 \times 3.530648}{\sqrt{(17 + 3.530648) \times 18}}

= \frac{60.02102}{22.63578}

= 2.651599

Nous constations que Z_1 > \lambda_1, Z_2 < \lambda_2, Z_3 < \lambda_3, par conséquent nous ne pouvons conclure que seul la valeur maximale de X peut être considérée comme un outlier.

– Le test \tau modifié de Thompson

Nous cherchons à déterminer les outliers au sein de de X. Nous avons \overline{X} = 6.605515. Calculons les \delta_i, \forall i \in [1,20],

\delta = (|1.9593 - 6.605515|, \cdots, |5.2654 - 6.605515|) = (4.6426215, \cdots, 1.340115)

L’Observation pour laquelle la valeur de \delta est maximale est celle qui correspond au maximum de X, soit \delta_{12} = 12.51899.

Maintenant, déterminons la valeur de référence pour un seuil de confiance de 5\%,

\tau = \frac{t_{\frac{0.05}{2},20-2} \times (20 - 1)}{\sqrt{20 \times (20 - 2 + t_{\frac{0.05}{2},20-2} ^2)}}

= \frac{t_{0.025,18} \times 19}{\sqrt{20 \times (18 + t_{0.025,18} ^2)}}

= \frac{-2.100922 \times 19}{\sqrt{20 \times (18 +4.413873)}}

= \frac{-39.91752}{\sqrt{448.2775}}

= \frac{-39.91752}{21.17256}

= -1.885342

Nous avons \sigma_X ^2 = 16.7722 et donc que,

\delta_{12} = 12.51899 > |-1.885342| \times 16.7722 = 31.61233

, nous en concluons que la valeur maximale de X est un outlier.

Lançons la seconde itération, nous supprimons donc la valeur maximale de X et obtenons,

\overline{X_{-12}} = 5.946621 et \sigma_{X_{-12}} ^2 = 8.538788

Après calcul de \delta, il en ressort que la valeur maximale correspond à la valeur minimale de X. Nous obtenons en valeur de référence et pour n = 19 observations,

\tau = -1.881106

Nous avons alors, \delta_1 = 3.9873211 < |-1.881106| \times 8.538788 = 16.06237. Nous en déduisons que la valeur minimale ne peut être considérée comme un outlier et stoppons l’algorithme à cette itération.

– Le critère de Peirce

Déterminons l’ensemble des outliers de X par le critère de Peirce. Nous débutons l’algorithme en calculons \overline{X} = 6.6055145, \sigma_X = 4.095388.

En nous reportant à la table de Peirce pour n = 20 et un outlier suspecté, nous obtenons R = 2.209. Par conséquent,

|X_i - \overline{X} |_{max} = 4.095388 \times 2.209 = 9.046712

Nous constatons que la seule observation remplissant le critère,

|X_i - \overline{X}| > 9.046712

, est la numéro douze, soit la valeur maximale de X.

Nous recherchons la présence d’un second outlier en mettant à jour la valeur de R = 1.914 depuis la table de Peirce pour n = 20 et deux outliers suspectés, nous obtenons,

|X_i - \overline{X} |_{max} = 4.095388 \times 1.914 = 7.838573

Aucune observation ne remplit le critère,

|X_i - \overline{X}| > 7.838573

Pour cette première itération nous en déduisons que X_{12} est un outlier et le supprimons du jeu de données.

Lançons une seconde itération en mettant à jour \overline{X_{-12}} = 5.946621, \sigma_{X_{-12}} = 2.92212. Nous obtenons R = 2.185 pour un outlier suspecté et n = 19, ce qui implique,

|X_i - \overline{X} |_{max} = 2.92212 \times 2.185 = 6.384832

Il se trouve qu’aucune observation ne remplit le critère,

|X_i - \overline{X_{-12}}| > 6.384832

Par conséquent, nous pouvons arrêter l’algorithme et conclure que seule la valeur maximale de X peut être considérée comme un outlier.

\bullet Application informatique:

Procédure SAS:

Package et fonction R:

\bullet Bibliographie: 

– Statistique. Dictionnaire encyclopédique de Yadolah Dodge

– Exploratory Data Analysis de John Wilder Tukey

– Simplified Statistics for Small Numbers of Observations de R. B. Dean et W. J. Dixon

– A Manual fo Spherical and Practical Astronomy de William A. Chauvenet

– Procedures for Detecting Outlying de F. E. Grubbs

– Rejecting Outliers in Factorial Designs de W. Stefansky

– Some Grubbs-Type Statistics for the Detection of Outliers de Gary Tietjen et Roger Moore

– Percentage Points for a Generalized ESD Many-Outlier Procedure de Barnard Rosner

– Outliers de John M. Cimbala

– Criterion for the Rejection of Doubtful Observations de Benjamin Peirce

L’Analyse de survie (En cours de rédaction…)

add.png

\bullet Introduction:

L’Analyse de survie est un outil statistique permettant de construire des indicateurs sur une durée de vie d’une variable aléatoire T afin de quantifier le temps de survie d’un instant à l’autre.

Universellement, nous notons S la fonction de survie à modéliser et de formule,

S(t) = P(T > t)

, avec t une variable temps.

La notion de données censurés aura également un impact important dans une analyse de survie. Une donnée censuré est une observation pour laquelle l’évènement que nous souhaitons mesuré ne s’est pas produit au cours de la période d’observation. Ces données portent également le nom de données tronquées.

De nombreux outils existent, les trois plus connus demeure la courbe de survie représentant le temps de survie, la régression de Cox pour le modéliser et l’estimateur de Kaplan-Meier pour founir un indicateur statistique.

\bullet La courbe de survie:

La courbe de survie se construit selon deux méthodes. S’il n’y a pas de données censurées alors il s’agit de la proportion de survivants en fonction du temps.

En présence de données censurées il convient de passer par l’estimateur de Kaplan-Meier afin de pouvoir construire la courbe de survie.

\bullet La régression de Cox:

Le modèle: La régression de Cox a été publiée en 1975 par David Cox dans l’objectif de modéliser des données de survies. Elle porte également le nom de modèle de Cox ou de modèle à risque proportionnel.

Le modèle de Cox modélisant la probabilité qu’un observation décède au temps t alors qu’il était en vie au temps t-1 s’écrit,

\lambda (t, X_1, \cdots, X_n) = \lambda_0 (t) e ^{\sum_{i = 1} ^n \beta_i X_i}

Avec \lambda_0 (t) le risque instantanée de décès lorsque (X_1, \cdots, X_n) sont nulles. Le modèle de Cox est semi-paramétrique du fait que nous n’ayons pas besoin d’estimer ce paramètre indépendant du temps t.

L’estimation: Les paramètres du modèle de Cox se détermine de manière classique par maximisation de la vraisemblance,

L(\beta) = \prod_{i = 1} ^m \frac{e ^{\beta ^t ^{\sum_{k \in D_i} X_k}}}{(\sum_{l \in R_i} e ^{\beta ^t X_l}) ^{d_i}}

Les coefficients ainsi déterminer permettent de donner le risque relatif entre la probabilité de décès et de survie au travers de la formule e ^{\beta_k}. A noter que cette interprétation est efficace sur des données indépendantes du temps, dans le cas inverse nous ne pouvons plus utiliser le risque relatif.

Les hypothèses de validité: Le modèle de Cox repose sur l’hypothèse forte des risques proportionnels et nécessite qu’elle soit vérifier avant application. Pour cela il s’agit de vérifier que chaque variable est indépendante du temps qu’il soit bénéfique, nocif ou nul.

Le risque relatif ajusté ou hazard ratio:

\bullet Estimateur de Kaplan-Meier:

Présentation: Publié en 1958 par Edward L. Kaplan et Paul Meier, l’estimateur de Kaplan-Meier, également connu sous le nom de l’estimateur produit-limite, est un estimateur pour estimer la fonction de survie.

L’estimateur:

Soit S(t) = P(T>t). Notons les durées observées jusqu’à chaque décès des membres de l’échantillon sont :

t_1 \leq t_2 \leq \cdots \leq t_N

À chaque n_i correspond un t_i, n_i étant le nombre de personnes « à risque » juste avant le temps, t_i et d_i le nombre de décès au temps t_i.

L’estimateur de Kaplan-Meier est l’estimation du maximum de vraisemblance non-paramétrique de <i>S</i>(<i>t</i>). C’est un produit de la forme:

\hat{S} (t) = \prod_{t_i < t} \frac{n_i - d_i}{n_i}

Lorsqu’il n’y aucune censure, n_i est le nombre de survivants juste avant le temps t_i. Lorsqu’il y a censure, n_i est le nombre de survivants moins le nombre de pertes (cas censurés). Ce sont seulement ces cas survivant qui continuent à être observés (qui n’ont pas encore été censurés) qui sont « à risque » de décès (observé).

\bullet Le test du log-rank:

Présentation: Publié par Nathan Mantel et David Cox, le test du log-rank (connu aussi sous le nom de test de Mantel-Cox) est une approche non paramétrique permettant de comparer les distributions de deux séries de survie X ^1, X ^2.

Le test:

Soit t \in[1, \cdots, T] les différents temps. Notons N_{1,t}, N_{2,t} le nombre de sujet à risques au temps t respectivement pour les séries X ^1, X ^2 et N_t = N_{1,t} + N_{2,t}.

Nous notons alors O_{1,t}, O_{2,t} le nombre d’évènements observés au temps t avec O_t = O_{1,t} + O_{2,t}.

Enfin nous avons,

E_{1,t} = \frac{O_t}{N_t} \cdot N_{1,t}

Et,

V_t = \frac{O_t \cdot \frac{N_{1,t}}{N_t} \cdot (1 - \frac{N_{1,t}}{N_t}) \cdot (N_t - O_t)}{N_t - 1}

\bullet Exemple:

Soit l’échantillon suivant,

add1.png

, où la variable « Temps » est le moment où nous observons une censure (« Survie = 1 ») ou non (« Survie = 0 »).

Estimateur de Kaplan-Meier:

Pour le groupe A, nous avons,

– Nous avons pour « Temps = 1 »,

\hat{S} (t = 1) = \frac{20 - 1}{20} = 0.95

– Puis pour « Temps = 1 à 2 »,

\hat{S} (t \in [1,2]) = \frac{19 - 0}{19} \times \hat{S} (t = 1) = 1 \times 0.95 = 0.95

– Puis pour « Temps = 1 à 3 »,

\hat{S} (t \in [1, 3]) = \frac{18 - 0}{18} \times \hat{S} (t \in [1,2]) = 1 \times 0.95 = 0.95

\cdots

– Puis pour « Temps = 1 à 18 »,

\hat{S} (t \in [1, 18]) = \frac{3 - 0}{3} \times \hat{S} (t \in [1,17]) = 1 \times 0.2724265 = 0.2724265

– Puis pour « Temps = 1 à 19 »,

\hat{S} (t \in [1, 19]) = \frac{2 - 0}{2} \times \hat{S} (t \in [1,18]) = 1 \times 0.2724265 = 0.2724265

– Puis pour « Temps = 1 à 20 »,

\hat{S} (t \in [1, 16]) = \frac{1 - 0}{1} \times \hat{S} (t \in [1,19]) = 1 \times 0.2724265 = 0.2724265

En procédant aux mêmes calcul pour le groupe B, nous obtenons alors les courbes de Kaplan-Meier suivantes,

add.png

Test du Log-Rank:

La courbe ci-dessus semble laisser penser qu’il n’y a pas de différence entre le groupe A et le groupe B. Vérifions le au sens statistique du terme,

Dans un premier temps, relevons les temps pour lesquelles nous n’avons pas de données censurées (survie = 0). Nous avons pour les groupes A et B:

T_A = (2,3,8,12,17,18,19,20) et T_B = (1,2,3,7,9,11,17,18,19)

Huit cas pour le groupe A et 9 cas pour le groupe B. Nous obtenons alors les temps fusionnés suivants,

T = (1,2,3,7,8,9,11,12,17,18,19,20)

Procédons itérativement,

– pour T = 1,

\begin{tabular}{|l|c|c|} \hline & <span style="color:#808080;">\sharp \lbrace T = 1 \rbrace</span> & \sharp \lbrace T > 1 \rbrace \\ \hline & Groupe A & 0 & 8 \\ \hline Groupe B & 1 & 8 \\ \hline \end{tabular}

\Rightarrow O_1 = O_{T = 1}, E_1 = E_{T = 1}, V_1 = V_{T = 1}

\Rightarrow O_1 = 1, E_1 = \frac{1 + 8}{8 + 9} = \frac{9}{17} = 0.5294118, V_1 = \frac{(0 + 8) \times (1 + 8) \times (0 + 1) \times ((8 + 8) - (0 + 1))}{(8 + 9) ^2 \times ((8 + 9) - 1} = \frac{8 \times 9 \times 1 \times 15}{17 ^2 \times 16} = \frac{1080}{4624} = 0.233564

\cdots

– pour T = 9,

\begin{tabular}{|l|c|c|} \hline & \sharp \lbrace T = 9 \rbrace & \sharp \lbrace T > 9 \rbrace \\ \hline & Groupe A & 0 & 5 \\ \hline Groupe B & 1 & 4 \\ \hline \end{tabular}

\Rightarrow O_6 = O_{T = 9}, E_6 = E_{T = 9}, V_6 = V_{T = 6}

\Rightarrow O_9 = 1, E_9 = \frac{1 + 4}{5 + 5} = \frac{5}{10} = 0.5, V_9 = \frac{(0 + 5) \times (1 + 4) \times (0 + 1) \times ((5 + 4) - (0 + 1))}{(5 + 5) ^2 \times ((5 + 5) - 1} = \frac{5 \times 5 \times 1 \times 8}{10 ^2 \times 9} = \frac{200}{900} = 0.222222222

\cdots

– pour T = 20,

\begin{tabular}{|l|c|c|} \hline & \sharp \lbrace T = 20 \rbrace & \sharp \lbrace T > 20 \rbrace \\ \hline & Groupe A & 1 & 0 \\ \hline Groupe B & 0 & 0 \\ \hline \end{tabular}

\Rightarrow O_{12} = O_{T = 20}, E_{12} = E_{T = 20}, V_{12} = V_{T = 20}

\Rightarrow O_{20} = 0, E_{20} = \frac{0 + 0}{1 + 0} = 0, V_{20} = \frac{(1 + 0) \times (0 + 0) \times (1 + 0) \times ((0 + 0) - (1 + 0))}{(1 + 0) ^2 \times ((1 + 0) - 1} = 0

En sommant les différents termes, nous avons,

O = \sum_{t = 1} ^{12} O_t =1 + 1 + 1 + 1 + 0 + 1 + 1 + 0 + 1 + 1 + 1 + 0 = 9

E = \sum_{t = 1} ^{12} E_t = 0.5294118 + 0.5 + 0.5 + 0.5 + 0.4545455 + 0.5 + 0.44444444 + 0.375 + 0.4285714 + 0.4 + 0.33333333 + 0 = 4.965306

V = \sum_{t = 1} ^{12} V_t = 0.2335640 + 0.4 + 0.3846154 + 0.2272727 + 0.2231405 + 0.222222 + 0.2160494 + 0.2008929 + 0.2448980 + 0.12000 - 0.22222222 + 0 = 2.250433

Par conséquent,

Z = \frac{}{\sqrt{}} =

, que nous reportons à la loi du \chi ^2 pour degrés de liberté et nous obtenons p = .

Régression de Cox:

\bullet Application informatique:

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

Package R:

Pour l’estimateur de Kaplan-Meier, https://cran.r-project.org/web/packages/epitools/epitools.pdf

Pour le reste, https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survfit.html

\bullet Bibliographie: 

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

– Nonparametric estimation from incomplete observations de Edward L. Kaplan et Paul Meier

– Evaluation of survival data and two new rank order de Nathan Mantel

– Asymptotically Efficient Rank Invariant Test Procedures de Richard Peto

– Linear Rank Tests in Survival Analysis de David Harrington

– The asymptotic properties of nonparametric tests for comparing survival distributions de D. Schoenfeld

– Dertimining the statistical significace of survivorship prediction models de H. P. Berty, H. Shi, J. Lyons-Weiler

– Analyse des durees de vie avec le logiciel R de Segolen Geray