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: fin 2016, 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

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.

La classification ascendante hiérarchique (En cours de rédaction…)

add2

\bullet Présentation:

La classification ascendante hiérarchique, régulièrement mentionnée sous l’acronyme CAH, est probablement la méthode d’analyse non supervisée la plus répandue. Elle 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 clusters dans un argo-anglicisme bien souvent utilisé pour faire simple, s’exécute par algorithme ascendant en définissant une fonction de proximité/distance à partir de laquelle la typologie des observations de \mathbf{X} sera construite au travers de partitions emboîtées d’hétérogénéités croissantes. Il existe un très grand nombre de fonctions de proximité ou de distance et chacune d’elles disposent de ses avantages et inconvénients.

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

Enfin, le choix du nombre de clusters se fait en fonction de deux approches: soit nous définissons un nombre de groupes maximales soit en fonction de quatre outils décisionnels que nous présenterons. Le choix entre ces deux méthodes est purement subjectif et principalement lié au contexte car aucune ne doit être considérée comme la meilleur.

\bullet Les trois familles de classification ascendante hiérarchique:

– La méthode classique

La CAH regroupe un ensemble de technique construites sur le mêm principe: le regroupement successif des points par ordre de proximité décroissante. L’Algorithme qui tourne quelque soit la fonction distance utilisé est le suivant:

  1. Considérer autant de classes qu’il y a d’observations
  2. Tant qu’il reste une observation à classer, à l’itération t

—- Calculer la matrice de distance entre les différents groupes d’observations et observations sans groupe et selon la fonction de proximité choisie.

—- La distance minimale d 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 une observation sans groupe ou deux observations sans groupe. La distance d retenue sert enfin à au dessin de l’arbre de classification d’ascendance hiérarchique.

—- Fin de l’itération t, nous passons à l’itération t+1 en reprenant toutes l’étape deux.

3. Une fois toutes les observations classées, nous disposons de l’arbre de classification hiérarchique présentant la structure du clustering et permettant de déterminer à quel niveau couper l’arbre pour visualiser à quel groupe appartiennent les différentes observations.

B. 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 5 \times 4 = 20 paires d’observations possibles est celle entre l’observation 2 et 3.

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.

add3

Nous constatons que cette fois-ci, la distance minimale est celle entre les observations 4 et 5. Evidemment, cette distance est supérieur à celle entre les observations 2, 3, d’où une paire de branche plus haute. Nous passons à l’itération 4 en mettant à jour avec le cluster regroupant les observations [late]4, 5[/latex] ) l’instar du cluster regroupant les observations 2, 3.

add4

La distance minimale est cette fois-ci celle entre le cluster 6 regrupant les observations 2, 3 et l’observation isolée 1. Nous mettons à jour l’arbre de classification avec ce nouveau regroupement et procédons à l’étape 5 qui sera la dernière.

add5.png

L’algorithme a finit de tourner et l’arbre dessiné permet de retrouver les différents clusters construits.

– La méthode de partitionnement

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.

– La méthode de noyau de densité

\bullet Les fonctions distances:

Le second paramètre à déterminer est la fonction à utiliser pour le calcul de la distance entre deux observations pour 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 l paramètre à fixer): d(X_{i_1}, X_{i_2}) = (\sum_{j = 1} ^P (X_{i_1} ^j - X_{i_2} ^j) ^l) ^{\frac{1}{l}}

\bullet Les procédures:

La distance moyenne

Origine: Robert Reuven Sokal et Charles Duncan Michener, 1958. Nous pouvons la trouvons 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 retenus.

add

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

La distance maximale

Origine: Thorvald Julius Sørensen, 1948. Nous pouvons la trouvons 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 retenus.

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, ce qui reste assez logique étant donné le fait qu’elle se base sur la distance maximale entre deux clusters. En général la distance maximale est très peu utilisée.

La distance minimale

Origine: K. Florek, J. Lukaszewicz, J. Perkal, S. Zubrzycki, L. L. McQuitty et P. H. A. Sneath, 1951-1957. Nous pouvons la trouvons 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 retenus.

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é) sur un jeu de données dont deux classes éloignées peuvent être reliées par un nombre faible d’observations proches entre elles et pouvant faire la jonction entre les deux groupes. Dés lors, le regroupement basé sur la distance minimale ne considérera pas 3 clusters (les deux groupes séparés et le lot d’observations entre eux) mais un unique cluster.

La distance des barycentres

Origine: Robert Reuven Sokal et Charles Duncan Michener, 1958. Nous pouvons la trouvons 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 sur chaque dimension des observations respectives à ces clusters. Les deux groupes dont la distance des barycentres est la plus petite sont ceux qui seront retenus.

add1

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

La distance des médianes

Origine: J. C. Gower, 1961. Nous pouvons la trouvons sous l’appellation median method.

L’idée:

La distance des médianes se détermine en calculant les « médianes » des différents groupes à partir de la médiane des coordonnées sur chaque dimension des observations respectives à ces clusters. Les deux groupes dont la distance des médianes est la plus petite sont ceux qui seront retenus.

add1

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 reste l’une des moins populaires et très peu souvent utilisée lors d’une Classification Ascendante Hiérarchique.

La distance de Ward

Origine: J. H. Ward, 1963. Nous pouvons la trouvons 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 respectifs des deux groupes considérés pondéré par leur effectif. En notant Cl_{k_1}, Cl_{k_2} les deux clusters considérés de barycentre respectif G_{k_1}, G_{k_2}, la formule d’usage est,

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 retenus.

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: W. S.Sarle et M. J. Symons, 1983. Nous pouvons la trouver sous l’appellation de méthode de EML (pour Expected Maximun-Likelihood).

L’idée:

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

D_{Cl_{k_1}, Cl_{k_2}} = n P log(1 + \frac{\Sigma - Sigma_{k_1} - Sigma_{k_2}}{\sum_{k = 1} ^K \Sigma_k}) - \lambda (n ^* log (n ^*) - n_{k_1} log (n_{k_1}) - n_{k_2} log (n_{k_2}))

Où, n désigne le nombre d’observations, P le nombre de dimensions (variables), \Sigma la matrice de variances-covariances issues du regroupement des clusters Cl_{k_1}, Cl_{k_2} et n ^* l’effectif obtenu lorsque que nous regroupons les clusters Cl_{k_1}, Cl_{k_2}. Enfin, le paramètre \lambda, fixé à 2 par défaut, représente la pénalité appliquée 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 retenus.

add1

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

La méthode du \beta-flexible

Origine: G. N. Lance et W. T. 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 à appliquer:

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

En général, le paramètre \beta est fixé à -0.25 par défaut. Les distances D_{Cl_{k_1}, Cl_{k_3}} et D_{Cl_{k_1}, Cl_{k_4}} correspondent à celles entre le groupe à fusionner aux deux autres groupes (maximisation de l’inertie inter-groupe) et la distance D_{Cl_{k_3}, Cl_{k_4}} 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 retenus.

add1

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

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}}})

A noter qu’une variable se basant directement sur le diamètre R de la sphère existe sous le nom de Uniform-kernel method, le paramètre n’est alors plus le nombre K mais R.

– 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. Cependant, il reste très influencé par le paramètre K pas toujours évident à optimiser.

La méthode du noyau uniforme

Origine:

L’idée:

Propriétés:

La méthode hybride de Wong

Origine:

L’idée:

Propriétés:

La distance de McQuitty

Origine:

L’idée:

Propriétés:

\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 pouvant être déduit.
 
– Le R _{it}^2, de formule,
R_{it} ^2 = \frac{\mbox{Inertie inter-classe}}{\mbox{Inertie totale}} = \frac{\sum_{i = 1} ^n || X_i - \overline{X} ||  ^2 - \sum_{k_{it} = 1} ^{K_{it}} \sum_{i_{k_{it}} \in Cl_{k_it}}} || X_{i_{k_{it}}} - \overline{X_{k_{it}}} || ^2}{\sum_{i = 1} ^n || X_i - \overline{X} ||  ^2}
, où K_{it} désigne le nombre de clusters à l’itération it de la procédure de classification. La formule induit que les observations qui seront encore non regroupée auront un poids nul sur le calcul du R_{it} ^2 puisque leur coordonnées sont confondus avec leur barycentre. Ainsi pour le cas où une observation équivaut à un groupe, le R ^2 vaut un.
– Le semi partiel R_{it} ^2, de formule,
sp-R_{it} ^2 = R_{it - 1} ^2 - R_{it} ^2
Le semi-partiel R_{it} ^2 caractérise l’évolution de la chute du R ^2 au cours des différentes itérations de la procédure de classification. A noter que lors de l’itération initiale (it = 1), R_1 = 1.
– Pseudo F, de formule,

pseudo-F = \frac{\frac{\sum_{i = 1} ^n || X_i - \overline{X} ||  ^2 - \sum_{k_{it} = 1} ^{K_{it}} \sum_{i_{k_{it}} \in Cl_{k_it}}} || X_{i_{k_{it}}} - \overline{X_{k_{it}}} || ^2}{K_{it} - 1}}{\frac{\sum_{k_{it} = 1} ^{K_{it}} \sum_{i_{k_{it}} \in Cl_{k_it}}} || X_{i_{k_{it}}} - \overline{X_{k_{it}}} || ^2}{n - K_{it}}}

 

– Pseudo t ^2 pour le regroupement des clusters G_1, G_2, de formule,

pseudo-t ^2 = \frac{\sum_{i \in Cl_{G_1} \union Cl_{G_2}} || X_i - \overline{X_{Cl_{G_1 \union G_2}} || ^2 - \sum_{i \in Cl_{G_1}} || X_i - \overline{X_{Cl_{G_1}} || ^2 - \sum_{i \in Cl_{G_2}} || X_i - \overline{X_{Cl_{G_2}} || ^2}{\frac{\sum_{i \in Cl_{G_1}} || X_i - \overline{X_{Cl_{G_1}} || ^2 + \sum_{i \in Cl_{G_2}} || X_i - \overline{X_{Cl_{G_2}} || ^2}{n_{G_1} + n_{G_2} - 2}}

 

– Cubic Criterion Clustering:

CCC = K log [\frac{1 - E[R ^2]}{1 - R ^2}] \cdot \frac{\sqrt{\frac{n p ^*}{2}}}{(0.001 + E[R ^2]) ^{1.2}}

Pour expliquer les calculs de E[R ^2] et p ^*, il convient de détailler directement toute l’algorithme de calcul à appliquer et mélangeant les deux. D’abord, il faut savoir que E[R ^2] représente la valeur théorique du R ^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 la pondération,

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

Que nous appliquons aux valeurs propres, nous obtenons alors le vecteur,

u = \frac{\lambda}{c}

Ce vecteur est la version initiale de celui à appliquer pour le calcul de E[R ^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 ^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 ^*} \lambda_j}{K}) ^{\frac{1}{p ^*}}</span></p> <p style="text-align:justify;"><span style="color:#808080;">, et,</span></p> <p style="text-align:center;"><span style="color:#808080;">[latex]u = \frac{\lambda}{c ^*}

Nous avons alors,

E[R ^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) ^2}{n}] \cdot [1 + \frac{4}{n}]

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

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

CCC > 2 : bonne classification, 0< CCC<2 : classification peut être OK mais à vérifier, CCC< 0 : présence d’outliers gênants (surtout si CCC < -30).

Nous traçons le CCC versus le nombre de classes. Un creux pour k classes suivi d’un pic pour k+1 classes indique une bonne classification en k+1 classes (surtout si on a une croissance ou décroissance douce à partir de k+2 classes). Rq : ne pas utiliser CCC et F en single linkage.

\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 la méthode de la distance minimale.

Nous initialisons l'algorithme en calculer la matrice des distances (selon la norme euclidienne) et obtenons,

\mathbf{D_1} = \begin{pmatrix} obs & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ 1 & 1.730431 & . & . & . & . & . & . & . & . \\ 2 & 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}

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} ^{10} (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} ^{10} (3.113284 ^2, 4.837503 ^2, \cdots, 5.482 ^2)}

= \frac{245.9374 - 0.4911644}{245.9374}

= 0.9980029

- sp-R_1 ^2 = 1 - R_1 ^2 = 1 - 0.9980029 = 0.0019971

- pseudo-F = \frac{\frac{\sum_{i = 1} ^{10} (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)}{9 - 1}}{\sum_{k_1 = 1} ^9 (0.4955625 ^2, 0 ^2, 0 ^2, 0.4955625 ^2, 0 ^2, \cdots, 0 ^2)}{10 - 9}}

= \frac{\frac{245.9374 - 0.4911644}{8}}{\frac{0.4911644}{1}}

= \frac{30.68078}{0.4911644}

= 62.4654

Le pseudo-t ^2 n'étant calculable pour cette première itération étant donné qu'il mesure l'évolution de l'inertie lorsque deux groupes de plusieurs observations sont regroupés. De même pour le CCC.

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

add3

Nous obtenons la matrice des distances suivantes,

\mathbf{D_2} = \begin{pmatrix} obs & A & 2 & 3 & 5 & 6 &7 & 8 & 9 \\ 2 & 1.998598 & . & . & . & . & . & . & . \\ 3 & 7.5475473 & 7.914920 & . & . & . & . & . & . \\ 5 & 3.210474 & 4.063529 & 5.562935 & . & . & . & . & . \\ 6 & 11.774134 & 13.191350 & 11.043181 & 10.893907 & . & . & . & . \\ 7 & 9.782901 & 11.462948 & 10.132024 & 8.972496 & 2.5189 45 & . & . & . \\ 8 & 8.023225 & 10.099191 & 10.548004 & 7.753578 & 5.446458 & 3.002310 & . & . \\ 9 & 7.349526 & 9.734812 & 12.841799 & 8.339789 & 9.553050 & 7.224126 & 4.260986 & / \\ 10 & 7.253523 & 9.681258 & 2.664138 & 7.963532 & 10.436960 & 7.986026 & 4.996152 & 1.66180 \\ end{pmatrix}

La distance la plus petite 1.327253 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

- sp-R_2 ^2 = R_1 ^2 - R_2 ^2 = 0.9980029 - 0.9944215 = 0.0035814

- pseudo-F = \frac{245.9374 - \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)}{8 - 1}}{\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)}{10 - 8}}

= \frac{\frac{245.9374 - 1.371965}{7}}{\frac{1.371965}{2}}

= \frac{34.93792}{0.6859825}

= 50.93121

Le pseudo-t ^2 et le CCC n'étant calculable pour cette seconde itération pour les même raisons décrites en première.

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çu.

add2

Nous obtenons la nouvelle matrice de distance ci-dessous,

\mathbf{D_3} = \begin{pmatrix} obs & B & A & 2 & 3 & 5 & 6 & 7 \\ A & 7.253523 & . & . & . & . & . & . \\ 2 & 9.681258 & 1.998598 & . & . & . & . & . \\ 3 & 12.664138 & 7.575473 & 7.914920 & . & . & . & . \\ 5 & 7.253523 & 3.210474 & 4.063529 & 5.562935 & . & . & . \\ 6 & 9.553050 & 11.774134 & 13.191350 & 11.043181 & 10.893907 & . & . \\ 7 & 7.224126 & 9.782901 & 11.462948 & 10.132024 & 8.972496 & 2.518945 & . \\ 8 & 4.260986 & 8.023225 & 10.099191 & 10.548004 & 7.753578 & 5.446458 & 3.002310 \\ end{pmatrix}

La distance minimale (1.378891) 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

- sp-R_3 ^2 = R_2 ^2 - R_3 ^2 = 0.9944215 - 0.9884517 = 0.0059698

- pseudo-F = \frac{245.9374 - \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)}{7 - 1}}{\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)}{10 - 7}}

= \frac{\frac{245.9374 - 2.840154}{6}}{\frac{2.840154}{3}}

= \frac{40.51621}{0.946718}

= 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 ^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 n'étant calculable pour cette troisième itération pour les même raisons décrites en première et seconde.

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 procédure selon la distance minimale,

addbbb

Calculons la nouvelle matrice des distances propre à l'itération en cours,

\mathbf{D_3} = \begin{pmatrix} obs & D & C \\ C & 3 & 4.141981 & . \\ 3 & 5.190977 & 9.287647 \\ end{pmatrix}

La distance minimale (4.141981) 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

- sp-R_8 ^2 = R_7 ^2 - R_8 ^2 \approx 0.7504299 - 0.2212299 \approx 0.592

- pseudo-F = \frac{245.9374 - \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)}{2 - 1}}{\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)}{10 - 2}}

= \frac{\frac{245.9374 - 191.5287}{1}}{\frac{191.5287}{8}}

= \frac{54.4087}{23.94109}

= 2.272607

- pseudo-t ^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 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 & 1.2699 & 9.1338 & 6.3236 & 0.9754 & 2.7850 & 5.4688 & 9.5751 & 9.6489 \\ 9.2858 & 10.7572 & 11.7537 & 9.3804 & 10.5678 & 1.1299 & 2.5688 & 3.4694 & 4.0119 & 5.3371 \\ \end{pmatrix} \times \begin{pmatrix} 8.1472 & 9.2858 \\ 9.0579 & 10.7572 \\ 1.2699 & 11.7537 \\ 9.1338 & 9.3804 \\ 6.3236 & 10.5678 \\ 0.9754 & 1.1299 \\ 2.7850 & 2.5688 \\ 5.4688 & 3.4694 \\ 9.5751 & 4.0119 \\ 9.6489 & 53371 \\ end{pmatrix} \times \frac{1}{10 - 1} = \begin{pmatrix}

= \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

, et nous permet de calculer le vecteur u,

[atex]u = (\frac{10.600604}{4.09116}, \frac{3.157856}{4.09116}) = (2.5910999, 0.7718731)[/latex]

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 donc,

E[R ^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 ^2], nous pouvons déterminer le CCC à cette itération,

CCC = 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

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écrite ci-dessus (d = 5.190977). 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

– A kth Nearest Neighbor Clustering Procedure de M. Anthony Wong et Tom Lane

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

– Cubic Clustering Criterion de W. S. Sarle

– A General Theory of Classificatory Sorting Strategies. I. Hierarchical Systems de G. N. Lance, G. N. et W. T. Williams

– A Comparison of Some Methods of Cluster Analysis de J. 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 L. L. McQuitty

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

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

– Some Methods for classification and Analysis of Multivariate Observations de J. B. MacQueen

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

– Une nouvelle méthode en classification automatique et reconnaissance des formes: la méthode des nuées dynamiques de Edwin Diday

– Le document en ligne: http://support.sas.com/documentation/onlinedoc/v82/techreport_a108.pdf

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é.

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

– 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é.

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

– 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.

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

– 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« .

Applied Regression Analysis in Econometrics de E. Doran Howard

Introduction to Econometrics de Christopher Dougherty:

 Tests of Equality Between Sets of Coefficients in Two Linear Regressions de Gregory C. Chow

– Le test de Mc Leod – Li

Présentation:

Le test de Mc Leod-Li permet d’étudier l’hétéroscédasticité conditionnelle d’un modèle ARCH.

Le test:

– Le test de Dickey-Fuller

Présentation:

Le test:

– Le test du bruit blanc

Présentation:

Le test:

– Le test de Franses

Présentation:

Le test:

– Le test de Zivot-Andrews

Présentation:

Le test:

– Le test de CUSUM

Présentation:

Le test:

– Le test de Gregory-Hansen

Présentation:

Le test:

– Le test de Hyllerberg-Engle-Granger-Yoo

Présentation:

Le test:

– Le test de Osborn-Chui-Smith-Birchenhall

Présentation:

Le test:

– Le test de la racine carré unité

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:

– Le test de Durbin

Présentation:

Le test:

– Le test de Geary

Présentation:

Le test:

– Le test de Moran

Présentation:

Le test:

\bullet Exemple:

\bullet Application informatique:

\bullet Bibliographie:

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:

\bullet Les différents tests:

– Le test de Crowley-Jones

Présentation: Benjamin Peirce, 1877

Le test:

(http://www.biecek.pl/statystykaMedyczna/2532042.pdf)

– Estimateur de Kaplan-Meier

Présentation: Également connu sous le nom de l’estimateur produit-limite, est un estimateur pour estimer la fonction de survie d’après des données de durée de vie.

Le test:

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 \le t_2 \le t_3 \le \cdots \le t_N.

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

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

\hat S(t) = \prod\limits_{t_i<t} \frac{n_i-d_i}{n_i}.

Lorsqu’il n’y aucune censure, ni est le nombre de survivants juste avant le temps ti.
Lorsqu’il y a censure, ni 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é).

– https://fr.wikipedia.org/wiki/Estimateur_de_Kaplan-Meier

– Le test de Gehan-Breslow-Wilcoxon

Présentation:

Le test:

– 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}

La statistique de test de Mantel-Cox est:

Z = \frac{\sum_{t = 1} ^T (O_{1,t} - E_{1,t}}{\sqrt{\sum_{t = 1} ^T V_t}}

La statistique de test suit une loi normale centrée-réduite.

– 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

– Le test de Fleming-Harrington

Présentation:

Le test:

– Le test de Peto-Peto-Prentice

Présentation:

Le test:

http://www.stata.com/manuals13/stststest.pdf

– Le test de Tarone–Ware

Présentation:

Le test:

http://www.stata.com/manuals13/stststest.pdf

– Le trend test

Présentation:

Le test:

http://www.stata.com/manuals13/stststest.pdf

– Le test de Gehan

Présentation:

Le test:

– Hazard Ratio

Présentation:

Le test:

\bullet Exemple:

\bullet Application informatique:

\bullet Bibliographie: 

Le test de Jonckheere

add

\bullet Présentation:

Publié en 1954 par Aimable Robert Jonckheere, le test de Jonckheere est une approche non paramétrique permettant de tester si les variables continues appariées X ^1, \cdots, X ^T sont appariées.

\bullet Le test:

Hypothèse préliminaire: Variables appariées continues pour T \geq 2.

La formule de la statistique de test de Jonckheere est,

J = \sum_{i = 1} ^n U_i

Où,

U_i = \sum_{t_1 = 1} ^{T-1} \sum_{t_2 = t_1 + 1} ^T I_{(X_i ^{t_1} < X_i ^{t_2})}

, est la statistique de Mann-Whitney adaptée aux échantillons appariés avec I(.) indicatrice valant notamment \frac{1}{2} si nous sommes en présence de deux termes ex-aequos.

A noter que le lien linéaire entre le U de Mann-Whitney et le test de Wilcoxon implique que la statistique de test de Jonckheere-Terpstra peut également se réécrire en fonction du W de Wilcoxon.

Si n est suffisamment grand, alors nous travaillons sur la statistique de test,

Z = | \frac{J - E[J]}{\sqrt{V(J)}} |

Avec,

E[J] = \frac{1}{4} n T (T - 1)

V(J) = \frac{n T (T - 1) (2 T + 5)}{72}

Ci-dessous la table de la loi normale centrée-réduite.addTendance pour le rejet de H_0:

Plus la statistique de test Z est grande et plus nous avons de chance de rejeter H_0. Ce qui revient à dire que,

Z \rightarrow \infty \Rightarrow J \rightarrow \infty

J est la somme des U_i de Mann-Withney qui explose si pour chaque observation nous pouvons constater que les valeurs sont toutes inférieures ou supérieures d’un temps successif à l’autre. Soit qu’au fur et à mesure que T croît, les valeurs de X augmentent ou diminuent marquant là un lien entre les deux informations.

\bullet Tendance lorsque n \longrightarrow \infty:

– Nous avons fait tourner cinq simulations afin d’étudier la robustesse du test de Jonckheere. Nous générons des échantillons de taille 10, 10^2, 10^3, 10 ^4 puis 10 ^5 de telle manière à ce que les différents sous-échantillons soient bien distincts et nous étudions si le fait d’augmenter le paramètre n a une influence sur le fait d’accepter H_0 à tort.

add

Nous constatons que nous rejetons à chaque fois l’hypothèse H_0, ce qui est en adéquation avec les hypothèses fixées.

– Nous avons fait tourner cinq simulations afin d’étudier la robustesse du test de Jonckheere. Nous générons des échantillons de taille 10, 10^2, 10^3, 10 ^4 puis 10 ^5 de telle manière à ce que les différents sous-échantillons soient légèrement confondus et nous étudions si le fait d’augmenter le paramètre n a une influence sur le fait de rejeter H_0 à tort.

add2

Nous constatons que passé un échantillon de taille N = 10 000 nous rejetons à tort H_0.

Nos simulations montrent bien que le test de Jonckheere est influencé par la taille d’échantillon.

\bullet Annexe théorique:

Cette partie présente la démonstration liant le U de Mann-Whitney et le test de Wilcoxon (W).

La statistique de Mann-Whitney est: U = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} et plus particulièrement:

– U_{X|_{g1}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{((X|_{g1})_{i_1} = (X|_{g2})_{i_1})})

U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{((X|_{g1})_{i_1} < (X|_{g2})_{i_2})} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

Il faut voir le calcul de la somme des rangs de manière inverse, c’est-à-dire qu’il faut comprendre que \forall l, (R|_{g1})_l est égal au \sharp \lbrace X|_{g2} < X|_{g1} \mbox{ a partir de l'indice } l \rbrace , remarquons que nous appliquons un décalage de l’indice l à chaque calcul de  (R|_{g1})_l. Dés lors nous avons que:

\sum_{l = 1} ^{n_1} (R|_{g1})_l = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_2}}) + \sum_{l = 1} ^n l
= U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2}

Or, U_{X|_{g1}} + U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}}) + \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

= \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}}) = n_1 \times n_2

Puisque 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}} = 1.

Par conséquent, \sum_{l = 1} ^{n_1} (R|_{g1})_l = U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} = n_1 \cdot n_2 - U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} et comme \sum_{l = 1} ^{n_1} (R|_{g1})_l est la somme des rangs provenant du test de Wilcoxon, le lien entre les deux statistiques est fait.

\bullet Exemple:

Soit les cinq échantillons appariés X ^1, X ^2, X ^3, X ^4, X ^5 suivants.

add

Ci-dessous, les boxplots associés à (X ^1, X ^2, X ^3, X ^4, X ^5).

add

Les boxplots montrent que X varie énormément en fonction du temps t \in \lbrace 1, 2, 3, 4, 5 \rbrace. Prouvons-le statistiquement.

Tout d’abord, commençons par déterminer la matrice du nombre de valeurs supérieurs à chaque temps t \in [1,5], nous obtenons,

U = \begin{pmatrix} 4 & 1 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 3 & 3 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 4 & 1 & 1 & 1 & 0 \\ 1 & 3 & 1 & 0 & 0 \\ 1 & 2 & 1 & 0 & 0 \\ 1 & 1 & 2 & 0 & 0 \\ 2 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 \\ 1 & 2 & 2 & 0 & 0 \\ 1 & 1 & 2 & 0 & 0 \\ 1 & 2 & 2 & 0 & 0 \\ 1 & 3 & 2 & 0 & 0 \\ 1 & 3 & 2 & 0 & 0 \\ 1 & 3 & 2 & 1 & 0 \\ 1 & 2 & 2 & 1 & 0 \\ 1 & 2 & 2 & 1 & 0 \\ 1 & 3 & 2 & 1 & 0 \\ 1 & 3 & 1 & 1 & 0 \\ \end{pmatrix}

Maintenant, nous pouvons calculer la statistique de test,

J = \sum_{i = 1} ^n \sum_{t = 1} ^T U_{i,t} = 4 + 1 + 1 + 1 + 0 + \cdots + 1 + 3 + 1 + 1 + 0 = 114

Reste à déterminer l’espérance et la variance de J. Nous avons,

E[J] = \frac{20 \times 5 \times (5 - 1)}{4} = \frac{400}{4} = 100

Et,

V(J) = \frac{20 \times 5 \times (5 - 1) \times (2 \times 5 + 5)}{72} = \frac{6000}{72} = 83.33333

Nous avons donc,

Z = | \frac{114 - 100}{\sqrt{83.33333}} | = | \frac{14}{9.128709} | = | 1.533623 | = 1.533623

Si nous reportons cette valeur de la statistique de test à la table de la loi normale centrée-réduite, nous obtenons une p-valeur de 0.06256121 > 5 \%. Nous en concluons que nous ne pouvons rejeter H_0 et qu’il n’y a d’influence du temps T sur X au sens statistique du terme.

\bullet Application informatique:

Procédure SAS:

%macro Jonckheere_test (TAB=);

/* La macro exécute le test de Jonckheere pour données appariées sur la table de données TAB contenant les variables continues */

/* Options pour nettoyer la log et les sorties */
options nonotes spool;
ods noresults;

/* Récupération du nombre d’observations n */
proc sql noprint;
select count(*) into: n from &TAB.;
quit;

proc contents data = &TAB. out = vars;
run;

/* Récupération du nombre de temps T et de la liste des variables appariées */
proc sql noprint;
select distinct(name) into: list_T separated by  »  » from vars;
select count(*) into: T from vars;
quit;

/* Initialisation de la statistique de test de Jonckheere */
%let J = 0;

%let Tmu = %eval(&T. – 1);

%do i = 1 %to &n.;

%do t1 = 1 %to &Tmu.;

%let tt1 = %eval(&t1. + 1);
%let c = 0;

/* Récupération de la valeur de référence à comparer à toutes les autres */
data _null_;
set &TAB.;
if _n_ = &i.;
call symputx(« X_ref »,%scan(&list_T.,&t1.));
run;

%do t2 = &tt1. %to &T.;

/* Récupération des valeurs à comparer à la valeur de référence */
data _null_;
set &TAB.;
if _n_ = &i.;
call symputx(« X_comp »,%scan(&list_T.,&t2.));
run;

/* Si inférieur alors on incrémente de 1 selon la définition de la fonction indicatrice */
%if &X_ref. < &X_comp. %then %do;
%let c = %sysevalf(&c. + 1);
%end;

/* Si égal alors on incrémente de 0.5 selon la définition de la fonction indicatrice */
%if &X_ref. = &X_comp. %then %do;
%let c = %sysevalf(&c. + 0.5);
%end;

%end;

/* Mise à jour de la statistique de test de Jonckheere */
%let J = %sysevalf(&J. + &c.);

%end;

%end;

/* Création de la table de résultats */
data Résultats;
Effectif = &n.;
Temps = &T.;
J = &J.;
Esperance = &n.*&T.*(&T.-1)/4;
Variance = &n.*&T.*(&T.-1)*(2*&T.+5)/72;
Z = abs((J-Esperance)/sqrt(Variance));
p_valeur = 1-probnorm(Z);
run;

/* Suppression des tables inutiles */
proc datasets lib = work nolist;
delete vars;
run;

/* Reset options */
options notes nospool;
ods results;

%mend;

Package et fonction R:

Jonckheere.test = function(TAB) {

# La fonction procède au test de Jonckheere pour données appariées sur la matrice TAB qui doit être composée de variables continues

# Suppression des données manquantes et récupération du nombre d’observations n et du nombre de temps T
TAB = na.omit(TAB)
n = dim(TAB)[1]
T = dim(TAB)[2]

# Création de la matrice synthétisant le nombre de valeur supérieur ou égale
TT = matrix(0,n,T)
for (i in 1:n) {
for(t1 in 1:(T-1)) {
c = 0
for (t2 in (t1+1):T) {
# Si la valeur est supérieur alors on incrémente de 1 selon la définition de l’indicatrice
if (TAB[i,t1]<TAB[i,t2]) {c = c + 1}
# Si la valeur est égale alors on incrémente de 0.5 selon la définition de l’indicatrice
if (TAB[i,t1]==TAB[i,t2]) {c = c + 0.5}
}
TT[i,t1] = c
}
}

# Calcul de la statistique de test de Jonckheere
J = sum(TT)
# Calcul de l’espérance
E = n*T*(T-1)/4
# Calcul de la variance
V = n*T*(T-1)*(2*T+5)/72

# Calcul de la statistique de test centrée-réduite
Z = abs((J-E)/sqrt(V))
# Calcul de la p-valeur selon la loi normale centrée-réduite
p = pnorm(Z,0,1,lower.tail=FALSE)

# Mise en forme et retour des résultats
return(data.frame(Pop = n, Temps = T, J = J, Esp = E, Var = V, Z = Z, p_valeur = p))

}

\bullet Bibliographie:

– Comparaison de populations. Tests non paramétriques de Ricco Rakotomalala

Le test de Fligner-Policello

add

\bullet Présentation:

Publié en 1981 par Michael A. Fligner et George E. Policello II, le test de Fligner-Policello, également appelé test des rangs robustes de Fligner-Policello, est une approche non paramétrique permettant de tester si X|_{Y = 1}, X|_{Y = 2}, les sous-échantillons d’une variable continue X restreinte aux deux groupes de Y, ont même médiane.

A noter que la manière dont a été pensé le test impose que passer une certaine taille d’échantillon, le test de Fligner-Policello est particulièrement consommateur en temps de calcul.

\bullet Le test:

Hypothèse préliminaire: X continue et Y binaire.

La formule de la statistique de test de Fligner-Policello est,

Z = \frac{\sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} - \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1}}{2 \sqrt{V_1 + V_2 + \overline{P_1} \overline{P_2}}}

Où,

\forall i_1 \in [1,n_1], (P_1)_{i_1} = \sum_{i_2 = 1} ^{n_2} [I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}} + \frac{1}{2} I_{(X|_{Y = 2})_{i_2} = (X|_{Y = 1})_{i_1}}]

\forall i_2 \in [1,n_2], (P_2)_{i_2} = \sum_{i_1 = 1} ^{n_1} [I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}]

Il s’agit en réalité des statistiques de test U_{X|_{g_1}}, U_{X|_{g_2}} de Mann-Whitney. Le lien linéaire existant entre elle et le W de Wilcoxon implique que le test de Fligner-Policello peut se réécrire en fonction de cette dernière. Enfin,

\forall k \in [1,2], \overline{P_k} = \frac{\sum_{i_k = 1} ^{n_k} (P_k)_{i_k}}{n_k}

\forall k \in [1,2], V_k = \sum_{i_k = 1} ^{n_k} [(P_k)_{i_k} - \overline{P_k}] ^2

La statistique de test suit une loi normale centrée-réduite et l’hypothèse H_0 est: « Les médianes sont égales / \theta_1 = \theta_2« .

Ci-dessous la table de la loi normale centrée-réduite.

add

Tendance pour le rejet de H_0:

Plus la statistique de test Z est grande et plus nous avons de chance de rejeter H_0. Ce qui implique,

\sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} - \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1} \rightarrow \infty

\Rightarrow \sum_{i_2 = 1} ^{n_2} (P_2)_{i_2} >>>> \sum_{i_1 = 1} ^{n_1} (P_1)_{i_1}

\Rightarrow \sum_{i_2 = 1} ^{n_2} \sum_{i_1 = 1} ^{n_1} [I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}] >>>> \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} [I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}} + \frac{1}{2} I_{(X|_{Y = 1})_{i_1} = (X|_{Y = 2})_{i_2}}]

\Rightarrow \forall i_1, i_2, \sum_{i_2 = 1} ^{n_2} \sum_{i_1 = 1} ^{n_1} I_{(X|_{Y = 1})_{i_1} < (X|_{Y = 2})_{i_2}} >>>> \sum_{i_1 = 1} ^{n_1} \sum_{i_1 = 2} ^{n_2} I_{(X|_{Y = 2})_{i_2} < (X|_{Y = 1})_{i_1}}

\Rightarrow \forall i_1, i_2, (X|_{Y = 2})_{i_2} > (X|_{Y = 1})_{i_1}

Soit que les valeurs de X sont ordonnées en fonction des groupes de Y et donc que les médianes sont différentes.

\bullet Tendance lorsque n \longrightarrow \infty:

Nous proposons ici de vérifier si le test de Fligner-Policello est sensible aux grands échantillons ou non.

Le tableau ci-dessous présente l’évolution des p-valeurs associées aux statistiques de test calculées sur plusieurs simulations dans le cas où les médianes sont différentes d’un groupe à l’autre.

add2

Globalement, quelque soit la taille de l’échantillon, le test statistique rejette H_0, ce qui est en accords avec nos hypothèses.

Procédons à la même expérience mais cette fois-ci dans un cas où les différentes médianes ne devraient pas être statistiquement différentes. Le tableau ci-dessous présente ces résultats.

add2

Jusqu’à N = 100 nous restons cohérent avec nos hypothèses, cependant nous voyons que pour un échantillon de taille N = 1000, le test rejette H_0 à tort.

Nous en déduisons que le test de Fligner-Policello est influencé par la taille de l’échantillon.

\bullet Annexe théorique:

Nous présentons ici une esquisse de la démonstration du lien entre le U de Mann-Whitney et le W de Wilcoxon.

La statistique de Mann-Whitney est: U = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} et plus particulièrement:

– U_{X|_{g1}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{((X|_{g1})_{i_1} = (X|_{g2})_{i_1})})

U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{((X|_{g1})_{i_1} < (X|_{g2})_{i_2})} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

Il faut voir le calcul de la somme des rangs de manière inverse, c’est-à-dire qu’il faut comprendre que \forall l, (R|_{g1})_l est égal au \sharp \lbrace X|_{g2} < X|_{g1} \mbox{ a partir de l'indice } l \rbrace , remarquons que nous appliquons un décalage de l’indice l à chaque calcul de  (R|_{g1})_l. Dés lors nous avons que:

\sum_{l = 1} ^{n_1} (R|_{g1})_l = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_2}}) + \sum_{l = 1} ^n l
= U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2}

Or, U_{X|_{g1}} + U_{X|_{g2}} = \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}}) + \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} < (X|_{g2})_{i_2}} + \frac{1}{2} \times 1_{(X|_{g1})_{i_1} = (X|_{g2})_{i_1}})

= \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} (1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}}) = n_1 \times n_2

Puisque 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} > (X|_{g2})_{i_2}} + 1_{(X|_{g1})_{i_1} = X|_{g2})_{i_2}} = 1.

Par conséquent, \sum_{l = 1} ^{n_1} (R|_{g1})_l = U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} = n_1 \cdot n_2 - U_{X|_{g2}} + \frac{n \cdot (n + 1)}{2} et comme \sum_{l = 1} ^{n_1} (R|_{g1})_l est la somme des rangs provenant du test de Wilcoxon, le lien entre les deux statistiques est fait.

\bullet Exemple:

Soit la variable aléatoire X distribuée selon deux groupes d’une variable Y:

add

Ci-dessous, les densités associées aux distributions de X selon nos deux groupes. Nous pourrons remarquer que nos deux sous-échantillons ont la même distribution et donc que les médianes semblent égales.

addDans un premier temps, déterminons les vecteurs,

P_1 = P((X|_{Y = 1}))

= (\sum_{i_2 = 1} ^{10} [I_{(X|_{Y = 2})_{i_2} < 8.1472} + \frac{1}{2} \times I_{(X|_{Y = 2})_{i_2} = 8.1472}], \cdots, \sum_{i_2 = 1} ^{10} [I_{(X|_{Y = 2})_{i_2} < 9.6489} + \frac{1}{2} \times I_{(X|_{Y = 2})_{i_2} = 9.6489}])

= (6,6,0,6,4,0,2,4,8,9)

\Rightarrow \sum_{i_1 = 1} ^{10} (P_1)_{i_1} = 6 + \cdots + 9 = 45

P_2 = P((X|_{Y = 2}))

= (\sum_{i_1 = 1} ^{10} [I_{(X|_{Y = 1})_{i_1} < 1.5761} + \frac{1}{2} \times I_{(X|_{Y = 1})_{i_1} = 1.5761}], \cdots, \sum_{i_1 = 1} ^{10} [I_{(X|_{Y = 1})_{i_1} < 9.5949} + \frac{1}{2} \times I_{(X|_{Y = 1})_{i_1} = 9.5949}])

= (2,10,8,3,5,2,3,8,5,9)

\Rightarrow \sum_{i_2 = 1} ^{10} (P_2)_{i_2} = 2 + \cdots + 9 = 55

A partir de ces deux objets, nous pouvons déterminer les éléments manquants,

\overline{P_1} = \frac{45}{10} = 4.5

\overline{P_2} = \frac{55}{10} = 5.5

V_1 = V_{X|_{Y = 1}}

= \sum_{i_1 = 1} ^{10} [(P_1)_{i_1} - 4.5] ^2

= (6 - 4.5) ^2 + \cdots + (9 - 4.5) ^2

= 86.5

V_2 = V_{X|_{Y = 2}}

= \sum_{i_2 = 1} ^{10} [(P_2)_{i_2} - 5.5] ^2

= (2 - 5.5) ^2 + \cdots + (9 - 5.5) ^2

= 82.5

Il ne nous reste plus qu’à calculer la statistique de test,

Z = \frac{55 - 45}{2 \times \sqrt{86.5+82.5 + 4.5 \times 5.5}} = \frac{10}{2 \times \sqrt{193.75}} = \frac{10}{27.83882} = 0.3592106

Si nous reportons la valeur de la statistique de test à la table de la loi normale centrée-réduite, nous obtenons une p-valeur de 0.7465 >>> 5 \%. Nous en concluons que nous ne pouvons rejeter H_0 et donc que les médianes sont égales.

\bullet Application informatique:

Procédure SAS: http://support.sas.com/documentation/cdl/en/statug/65328/HTML/default/viewer.htm#statug_npar1way_details20.htm

Package et fonction R: http://finzi.psych.upenn.edu/R/library/RVAideMemoire/html/fp.test.html

\bullet Bibliographie:

– Robust Rank Procedures for the Behrens-Fisher Problem de M. A. Fligner et G. E. Policello