L’Aire sous la courbe ROC

add\bullet Présentation:

Conçu en 1941, le calcul de l’aire sous la courbe ROC (Receiver Operating characteristic) est une approche  non paramétrique permettant de mesurer l’association entre une variable continue X et une variable binaire Y.

Soit X|_{g_1} et X|_{g_2} les sous-échantillons de X restreint aux groupes 1 et 2 de Y, la courbe ROC permet ainsi de détermination la spécificité et la sensibilité de la distribution issues de ces deux sous-échantillons. Elle trouve notamment tout son intérêt quand X est un vecteur décisionnel issue d’une modèle prédictif afin de mettre en évidence les performances du modèle construit.

Le calcul de l’aire sous la courbe ROC (AUC, Area Under the Curbe) a été étendu à Y avec L \geq 2 modalités. Nous parlons alors de VUS (Volume Under the Surface ROC) dans le cas trichotomique et de HUM (Hypervolume Under the Manifold) dans le cas polychotomique strictement supérieur à 3.

\bullet L’indice:

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

La courbe ROC:

La courbe ROC permet de représenter l’évolution du taux de classification de Y à partir des valeurs de X. En faisant varier un cutoff c \in [min(X); max(X)] nous pouvons déterminer les deux scores suivant (dont les définitions sont parfois inversées si nous sommes sur des conventions françaises ou anglaises):

– la spécificité qui est la probabilité de prédire correctement \overbrace{Y} = \mbox{classe } 1 sachant que Y = \mbox{classe } 1 pour le cutoff c fixé

– la spécificité qui est la probabilité de prédire correctement \overbrace{Y} = \mbox{classe } 2 sachant que Y = \mbox{classe } 2 pour le cutoff c fixé

, et tracer la courbe passant par les points de coordonnées (Sensibilite_c, 1 - Specificite_c)

Ces taux de bonne classification sont basés sur la théorie de Bayes et sont très prisés dans la littérature. Nous sommes donc dans un plan à deux dimensions, la diagonale passant par les points de coordonnées (0,0) et (1,1) est la droite de référence d’une classification aléatoire. Plus la courbe colle le côté haut-gauche du plan plus la classification est de bonne qualité.

Ci-dessous un exemple de courbe ROC:

add

L’ aire sous la courbe ROC:

Le calcul de l’aire sous la courbe ROC dans le cas binaire consiste en celui des paires concordantes/discordantes. Il s’agit d’une mesure empirique qui offre un très bon indicateur de la répartition de X par rapport aux deux groupes de la variable Y. En effet, dans le cadre d’une liaison entre X et Y, nous nous attendons à ce que les valeurs de X pour le groupe 2 de Y soient généralement plus fortes que celles pour le groupe 1.

L’AUC tient son nom d’indicateur empirique du fait qu’il ne travaille pas sur la forme continue des variables mais plutôt sous un aspect complexe de leur rang.

Rappelons, pour X|_{g_1} et X|_{g_2} distributions de X pour les deux groupes de Y et i_1, i_2 deux observations des groupes 1 et 2 du jeu de données observés, une paire est dites:

– concordante, si (X|_{g_2})_{i_2} > (X|_{g_1})_{i_1}

– discordante, si (X|_{g_2})_{i_2} < (X|_{g_1})_{i_1}

– ex-aequo, si (X|_{g_1})_{i_1} = (X|_{g_2})_{i_2}

Le calcul de l’AUC se fait alors par la formule:

AUC = \frac{1}{n_1 \cdot n_2} [\sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} 1_{(X|_{g_1})_{i_1} > (X|_{g_2})_{i_2}} + \frac{1}{2} \sum_{i_1 = 1} ^n \sum_{i_2 = 1} ^n 1_{(X|_{g_1})_{i_1} = (X|_{g_2})_{i_2}}]

L’AUC varie entre [\frac{1}{2} - 1]. Il se peut qu’il soit inférieur à \frac{1}{2}, dans ce cas il faut rectifier sa valeur par 1 - AUC. Cet effet s’explique par le fait que nous avons supposé que le groupe de référence admet en général des valeurs fortes par rapport à l’autre groupe alors qu’en réalité les données observées tendent vers le contraire (par exemple, nous nous attendions à ce que X dans le groupe 2 ait généralement des valeurs supérieures à celles du groupe 1 alors que c’est l’inverse).

Lorsque l’AUC est autour de \frac{1}{2} cela indique que les distributions de X|_{g_1} et X|_{g_2} sont confondues et que nous sommes dans un cas similaire à un tirage aléatoire de X par rapport à Y (la courbe ROC correspondante sera proche de la diagonale). Plus l’AUC s’approche de 1 et plus les deux distributions sont séparées (la courbe tend à ressembler à un carré de côtés de longueurs 1). Ainsi une AUC comprise entre,

]0.5-0.6] implique que X|_{g_1} et X|_{g_2} sont quasi-identiquement distribuées

]0.6-0.7] implique que X|_{g_1} et X|_{g_2} sont légèrement bien ordonnées

]0.7-0.8] implique que X|_{g_1} et X|_{g_2} sont globalement bien ordonnées

[0.8-0.9] implique que X|_{g_1} et X|_{g_2} sont bien distinctement ordonnées

]0.9-1] implique que X|_{g_1} et X|_{g_2} sont quasi-parfaitement séparées

Tendance de l’indicateur:

La force de l’AUC réside dans le fait qu’il offre une mesure relative à chaque observations par rapport à toutes les autres. Par exemple, pour deux variables X ^1, X ^2 différentes comprenant chacune une unique observation attendue dans le groupe 1 mais finalement observées dans le 2, si cette observation est plus proche du groupe attendue pour X ^1 que pour X ^2 alors la première variable aura une meilleure AUC que la seconde.

L’hypervolume sous hyper-courbe:

L’AUC a été étendue au cas multiclasse (comprendre L > 2), devenant plus généralement le HUM (Hypervolume Under the Manifold). Son calcul est similaire à l’AUC, il faut désormais ne plus voir les couples concordants/discordants mais la cohérence des L-uplets présents. Le HUM varie ainsi entre [\frac{1}{L!} - 1].

Un exemple de HUM dans le cas à 3 classes (nous parlons alors de VUS, Volume Under the Surface):

add

Sa formule générale devient alors:

HUM = \frac{1}{n_1 \times \cdots \times n_L} [\sum_{i_1 = 1} ^{n_1} \cdots \sum_{i_L = 1} ^{n_L} 1_{(X|_{g_1})_{i_1} > \cdots > (X|_{g_L})_{i_L}}]

Où,

1_{(X|_{g_1})_{i_1} > \cdots > (X|_{g_L})_{i_L}} = 1, si chacune des L-1 conditions de l’indicatrice sont respectées

1_{(X|_{g_1})_{i_1} > \cdots > (X|_{g_L})_{i_L}} = 0, si au moins une seule condition parmi L de l’indicatrice n’est pas respectée

1_{(X|_{g_1})_{i_1} > \cdots > (X|_{g_L})_{i_L}} = \frac{1}{(K+1)!}, si l’ordre des groupes est respecté mais que pour certains cas nous sommes en présence de K ex-aequos

\bullet La première contrainte à cette généralisation de l’AUC est que plus le nombre de classes augmente, plus le nombre de comparaisons augmente et devient ainsi rapidement coûteux. L’idée proposée est de regrouper les classes  qui varient entre elles (par exemple nous nous attendons à ce que  les premières classes varient beaucoup entre elles mais très peu avec les dernières) et de calculer des HUM sur des sous-groupes de classes.

\bullet La seconde est que plus nous rajoutons de classes et plus l’indicateur HUM devient pessimiste (en effet nous voyons bien que plus nous rajoutons de conditions dans les indicatrices et moins nous avons de chance de toutes les satisfaire), ainsi un excellent classifieur tendra vers un classifieur médiocre rendant la mesure HUM obsolète au fur et mesure que le nombre de groupes augmente.

\bullet Le troisième étant la détermination de l’ordre des groupes afin d’exploiter l’indicateur HUM de manière optimale. Plusieurs critères sont utilisables comme la moyenne (peu robuste) ou la médiane (plus réaliste de part son aspect empirique qui permet de cibler 50% de la population de chacun des groupes et ainsi la tendance générale d’un groupe par rapport aux autres pour les ordonner). La façon optimale de procéder étant de calculer les HUM de l’ensemble des différentes combinaisons de classes et retenir la valeur la plus forte.

Intervalle de confiance et test de Delong:

L’intervalle de confiance et le test de Delong permettent de construire un test statistique concluant si deux valeurs du HUM sont statistiquement différentes l’une de l’autre.

L’intervalle de confiance se calcul par bootstrap. Soit B le nombre de tirages aléatoires avec remises sur l’échantillon de base, nous pouvons alors définir l’ensemble des indicateurs de performance HUM calculés sur chacun des B sous-échantillons: \lbrace HUM_b: b = 1 \cdots B \rbrace.

L’erreur standard (se) bootstrap du HUM est:

se_B(HUM) = \sqrt{\frac{1}{B-1} \sum_{b = 1} ^B (HUM_b - HUM)^2}

Enfin, au risque \alpha fixé, l’intervalle de confiance sera:

[\overbrace{HUM} - z_{\frac{\alpha}{2}} \cdot se_B(HUM) ; \overbrace{HUM} + z_{\frac{\alpha}{2}} \cdot se_B(HUM)]

Le test de Delong permet, à la base, de comparer deux hyper-surfaces ROC S_1 et S_2 entre elles, sa formule est:

z = \frac{HUM_{S_1} - HUM_{S_2}}{\sqrt{var(HUM_{S_1}) + var(HUM_{S_2})}} \sim N(0,1)

En supposant que l’une des deux hyper-surfaces ROC ait pour valeur de HUM la référence aléatoire \Theta_L, nous pouvons nous servir de ce cas particulier afin de juger si l’indicateur de performance HUM est différent d’une répartition aléatoire au sens statistique. La formule devient alors:

z = \frac{HUM - \Theta_L}{\sqrt{var(HUM)}} \sim N(0,1)

\bullet Annexe théorique:

Nous noterons,

TC_l = P[ \mbox{classe predite = }l | \mbox{classe reelle = }l]

, soit la probabilité de classer correctement une observation de la classe l, l \in [ 1, \cdots, L].

En posant L-1 cutoffs (c_1, \ldots, c_{L-1}), tel que c_1 < c_2 < \cdots < c_{L-1}, que nous faisons varier sous contraintes, et en posant l’algorithme décisionnel suivant pour une observation i, i \in [ 1, \cdots, n ], de valeur X_i:

– Si X_i \leq c_1 \Longrightarrow i \in classe 1

– Si c_1 < X_i \leq c_2 \Longrightarrow i \in classe 2

– Si c_2 < X_i \leq c_3 \Longrightarrow i \in classe 3

\vdots

– Si X_i > c_{L-1} \Longrightarrow i \in classe L

Nous pouvons calculer le L-uplets de performances (TC_1, \ldots, TC_L) à partir de l’algorithme décisionnel ci-dessus. L’indicateur HUM se base alors sur l’évolution des taux de bonne classification relativement entre eux, nous pouvons alors définir une fonction:

f^{(L-1)}:(TC_1, \ldots, TC_{L-1}) \longmapsto TC_L

[0,1] ^{L-1} \longmapsto [0,1]

L’hyper-surface ROC est alors tracée à partir de ce L-upletUn indicateur de performance est représenté par l’aire sous l’hyper-surface ROC:

HUM = \int_0 ^1 \int_0 ^{f^1 (TC_1)} \ldots \int_0 ^{f^{L-2}(TC_1, \cdots, TC_{L-2})} f^{(L-1)}(TC_1, \ldots, TC_{L-1}) \partial(TC_{L-1}) \ldots \partial(TC_1)

En posant X|_{g_l} les valeurs de X dans la classe l et G_l(.) la fonction de répartition de X|_{g_l}, nous pouvons réecrire le L-uplet:

(TC_1, \ldots, TC_L) = (G_1(c1), \ldots, 1-G_L(c_{L-1}))

Les cutoffs (c_1, \cdots, c_{L-1}) peuvent donc s’écrire:

\left\{ \begin{array}{ll} c_1 = G_1 ^{-1} (TC_1) \\ c_2 = G_2 ^{-1} \left[ TC_2 + G_2(G_1 ^{-1}(TC_1)) ] = G_2 ^{-1} [ TC_2 + G_2(c_1) \right] \\ c_3 = G_3 ^{-1} \left[ TC_3 + G_3(G_2 ^{-1}(TC_2 + G_2(G_1 ^{-1}(TC_1)))) \right] = G_3 ^{-1} \left[ TC_3 + G_3(c_2) \right] \\ \vdots \\ c_{L-1} = 1 - G_L \left[ G_{L-1} ^{-1}(TC_{L-1} + G_{L-1}(G_{L-2}^{-1}(TC_{L-2} + G_{L-2}(c_{L-3})) \right] \end{array} \right.

Alors nous pouvons écrire:

HUM = P[(X|_1)_{i_1} < (X|_2)_{i_2} < \ldots < (X|_L)_{i_L}]

= \int_0 ^1 \int_0 ^{f^1(TC_1)} \ldots \int_0 ^{f^{(L-2)} (TC_1, \ldots, TC_{L-2})} f^{(L-1)} (TC_1, \ldots, TC_{L-1}) \partial TC_{L-1} \ldots \partial TC_1

= \int_0 ^1 \ldots \int_0 ^{f^{L-2} (TC_1, \ldots, TC_{L-2})} 1 - G_L [ G_{L-1} ^{-1}(TC_{L-1} + G_{L-1}(G_{L-2}^{-1}(TC_{L-2} + G_{L-2}(c_{L-3}))] \partial TC_{L-1} \ldots \partial TC_1

Une écriture du HUM, autre que l’écriture sous forme d’intégrale, et plus propice au calcul est la suivante:

\overbrace{HUM} = \frac{1}{n_1 \times \ldots \times n_L} \cdot \sum_{i_1 = 1} ^{n_1} \sum_{i_2 = 1} ^{n_2} \cdots \sum_{i_L = 1} ^{n_L} CR((X|_{g_1})_{i_1}, (X|_{g_2})_{i_2}, \ldots, (X|_{g_L})_{i_L})

Où, CR =

+ 1 \mbox{ si } (X|_{g_1})_{i_1} < (X|_{g_2})_{i_1} < \ldots < (X|_{g_L})_{i_1}

\vdots

+ \frac{1}{k!} pour k classes ex-aequos et concordantes parmi L par rapport aux autres

\vdots

0 sinon

La démonstration du \overbrace{HUM} se base sur le fait que HUM coïncide avec le pourcentage de bonne classification (PC). Nous pouvons rappeler que:

PC = \int \cdots \int g ^1 (X|_{g_1}) \cdots g ^L (X|_{g_L}) \times CR (X|_{g_1}, \cdots, X|_{g_L}) \partial X|_{g_1} \cdots \partial X|_{g_L}

Or, une autre manière de le formuler est le taux de L-uplets concordants d’où:

\overbrace{HUM} = \widehat{PC}

Le HUM varie alors dans [ \frac{1}{L!} - 1 ]. La valeur seuil synonyme d’aléa pour L classes notée D_L (1) est calculée par:

D_L (1) = \frac{1}{2\cdot(L-1)!} \int_{-\infty} ^{1-} x^{L-1} \cdot sign(x) + \sum_{k=1}^L (-1)^k \cdot C_L ^k \cdot (x - k) ^{L-1} \cdot sign(x - k) dx

= \frac{1}{L!}

\bullet Exemple:

Soit les variables X, Y ci-dessous,

add

Nous présentons la représentation dotplots qui permet de comparer  X|_{g_1}, X|_{g_2}, les distributions restreintes de X aux deux groupes de Y.

add

La courbe ROC associée est,

add

Une façon « propre » de résumer le calcul de l’AUC est d’ériger le tableau de taille n_1 \times n_2 contenant le score de concordance des différentes observations i_1 \in [1, \cdots, n_1] du groupe 1 par rapport aux i_2 \in [1, \cdots, n_2] du groupe 2. Nous obtenons alors le tableau T, pour n_1 = n_2 = 10:

add

Par exemple, l’observation N° 1 qui appartient au groupe 1 est attendue comme ayant une valeur X_1 < X_j, \forall j, j variant dans l’ensemble des numéros d’identifiant d’observation du groupe 2X_1 = 9.2858 et X|_{g_2} = (11.0759, 9.054, 10.5308, 11.7792, 9.934, 16.1622, 17.7943, 18.3112, 19.5285 ,20.1656)

Par conséquent, le vecteur des différentes 1_{(X|_{g1})_1 < X|_{g2}},  \forall j est (1,0,1,1,1,1,1,1,1,1) puisque qu’hormis avec l’observation N° 12 avec laquelle l’observation N° 1 a une valeur supérieure, pour toutes les autres les paires sont concordantes.

Une fois ce tableau obtenu (notons que nous n’avons aucun ex-aequo sinon quoi le tableau aurait contenu des \frac{1}{2}), il suffit de sommer chacune des cellules entre elles et de diviser le produit des effectifs des deux groupes, ainsi:

AUC = \frac{1}{n_1 \times n_2} \times (\sum_{i = 1} ^{n_1} \sum_{j = 1} ^{n_2} T_{i,j}) = \frac{1}{100} \times 88 = 0.88

Nous en déduisons que la liaison entre X et Y est forte avec un tel score.

\bullet Application informatique:

Les logiciels SAS et R ne permettent que de procéder au calcul du HUM mais uniquement de l’AUC. Nous proposons ici une macro SAS et une fonction R offrant la possibilité de calculer le HUM quelque soit le nombre de classes. La version R du code est beaucoup plus rapide que celle de SAS du fait que SAS, sans le module IML, n’est pas vraiment optimisé pour la comparaison de k-uplets concordants comparé à R qui fait la différence de part son aspect algorithmique plus poussé. Les deux programmes calculent le HUM pour les différentes permutations de classes possibles afin de connaître l’ordre optimal de Y sachant la distribution de X. A noter que ces fonctions n’ont pas subies la vague de tests nécessaires pour assurer de leur fiabilité, aussi il est fortement recommandé de vérifier les résultats depuis les tables k-uplets retournées.

Macro SAS:

%macro HUM(DATA=,Y=);

         /* Ce programme calcule l’hypervolume sous la courbe ROC (HUM). La structure attendue du jeu de données DATA est: 2
colonnes – 1 variable réponse Y (format numérique ou caractère) / 1 variable explicative (format numérique). Il
fournit en sortie les différents scores de HUM pour les différentes combinaisons de classe de la variable réponse
(versus). Remerciement à Olivier Decourt pour son aide non négligeable et aux forumeurs Suistrop et Jerome_pdv2 du
forum developpez.net pour leur intérêt constructif */

/* 1) Options de simplification de la log */
options nonotes spool;

/* 2) Verificaton N°1 si bonne structure des données */
proc contents data = &DATA. out = contents noprint;
run;

proc sql noprint; /* Décompte du nombre de variables du jeu de données */
select count(name) into: nb_vars from contents;
quit;

%if &nb_vars. ^= 2 %then %do; /* Si on a 1 ou plus de deux variables alors le programme ne peut s’exécuter car pas adapté */
%put La structure n est pas conforme à celle attendue, vérifier qu une variable réponse et une variable explicative (uniquement) sont bien présentes, puis relancez.;
%end;

%if &nb_vars. = 2 %then %do; /* Si on a bien 1 variable réponse et 1 variable explicative alors le programme peut s’exécuter */

 /* 3) Informations sur les classes */
proc freq data = &DATA.;
tables &Y. / out = recap_classe;
ods exclude all;
run;

proc sql noprint;
select distinct(&Y.) into :liste_classe_O separated by ‘ ‘ from recap_classe; /* Nom des différentes classes */
select count(&Y.) into :nb_classe from recap_classe; /* Nombre totale de classes */
quit;

data _null_;
set contents;
if name ^= « &Y. »;
call symputx(« x »,name); /* Nom de la variable x */
call symputx(« format_x »,type); /* Vérification N°2 si bonne structure des données */
run;

%if &format_x. ^= 1 %then %do; /* Si le format n’est pas sur 1 c’est qu’il ne s’agit pas d’une variable en format numérique, le programme n’est dont pas adapté */
%put La structure n est pas conforme à celle attendue, vérifier que la variable explicative est bien en format numérique, puis relancez.;
%end;

%if &format_x. = 1 %then %do; /* Si le format est sur 1 c’est qu’il s’agit bien d’une variable au format numérique alors le programme peut s’exécuter */

data _null_;
set contents;
if name = « &Y. »;
call symputx(« y_long »,length); /* Taille de Y */
run;

 /* 4) Convertion des labels des classes en numérique afin d’homogénéiser le processus */
data norm;
set &DATA.;
%do it_norm = 1 %to &nb_classe.;
if &Y. = « %scan(&liste_classe_O.,&it_norm.) » then y_num = &it_norm.; /* Conversion des classes en format numérique */
%end;
if &X. ne .;/* Suppression des données manquates */
run;

proc sql noprint;
select distinct(y_num) into :liste_classe separated by ‘ ‘ from norm; /* Liste des numéros de classe */
quit;

 /* 5) Récupération des effectifs par classe */
%do it_eff = 1 %TO &nb_classe.;

proc sql noprint;
select count(*) into :eff from norm (where = (y_num = %scan(&liste_classe.,&it_eff.))); /* Effectif par classe */
quit;

 /* Création de la macro-variable eff contenant les effectifs par classe */
%if &it_eff. = 1 %then %do;
%let effectif = [ &eff. ( %scan(&liste_classe_O.,&it_eff.) );
%end;
%if (&it_eff. ^= 1 and &it_eff. ^= &nb_classe.) %then %do;
%let effectif = &effectif. x &eff. ( %scan(&liste_classe_O.,&it_eff.) );
%end;
%if &it_eff. = &nb_classe. %then %do;
%let effectif = &effectif. x &eff. ( %scan(&liste_classe_O.,&it_eff.) ) ];
%end;

%end;

/* 6) Table des combinaisons des classes */
data hum;
array x [&nb_classe.] (&liste_classe.);
nfact = fact(dim(x));
do i = 1 to nfact;
call allperm(i, of x[*]);
output;
end;
call symputx(« nb_combi »,nfact);
drop nfact;
run;

 /* 7) Calcul du HUM */
proc transpose data = hum (drop = i) out = versus;
run;

%do it_combi = 1 %TO &nb_combi.; /* Parcours de toutes les combinaisons de classe disponibles */

proc sql noprint;
select col&it_combi. into :liste_classe_versus separated by ‘ ‘ from versus; /* Liste des numéros de classe du versus en cours */
quit;

%do it_ordre = 1 %to &nb_classe.; /* Codage des classes afin de respecter l’ordre du versus étudié */

data norm;
set norm;
if y_num = « %scan(&liste_classe_versus.,&it_ordre.) » then y_ordre = &it_ordre.;
run;

%end;

proc sort data = norm;
by y_ordre;
run;

data format; /* Mise sous format des données afin de pouvoir simplifier le système de comparaison */
set norm end = fin;
by y_ordre;
if first.y_ordre then do;
__y_ordre + 1;
start = 0;
end;
start + 1;
type = « N »;
fmtName = cats(« x »,__y_ordre, »_ »); /* Création de la catégorie de chaque valeur numérique */
label = &x.;
if last.y_ordre then call symputx(cats(« n »,__y_ordre), start); /* Nombre d’observations par classe */
run;

proc format cntlin = format;
run;

data kuplets (keep = c d t tot);
array x {&nb_classe.};
array i {&nb_classe.};
format %do z = 1 %to &nb_classe.; /* Utilisation des formats */
i&z. x&z._.
%end;
;
%do z = 1 %to &nb_classe.; /* Adaptation au nombre d’observations par classe */
do i&z. = 1 to &&n&z..;
%end;
do j = 1 to &nb_classe.; /* Initialisation du vecteur de valeurs à comparer */
x[j] = vvalue(i[j]);
end;
concord = 0;
exaequo = 0;
discord = 0;
do j = 1 to (&nb_classe. – 1);
x_j = vvalue(x[j]); /* On isole la valeur en cours à comparer au reste du vecteur */
x[j] = .; /* On met à vide la cellule en cours et qui a déjà été isoler afin de la soustraire à la phase de comparaison */
min = min(of x[*]); /* On cherche la plus petite valeur du vecteur restant à comparer */
if x_j > min then do; /* Si la valeur en cours est plus grande que la valeur minimale du reste du vecteur à comparer alors c’est discordant car on s’attend à une variable up-régulée, on peut stopper directement les comparaisons */
discord = 1;
goto fin;
end;
else if x_j = min then exaequo = exaequo + 1; /* Si la valeur en cours est égale à la valeur minimale c’est qu’on a un ex-aequo */
end;
fin :
if discord = 1 then do; /* Si on a discordance, alors on se fiche du nombre d’ex-aequos */
exaequo = 0;
concord = 0;
end;
if discord = 0 and exaequo = 0 then concord = 1; /* Si on a ni discordance ni ex-aequos alors ça implique qu’on a un vecteur totalement concordant */
if discord = 0 and exaequo > 0 then do; /* Si on a aucune discordance mais des ex-aequos alors on applique la correction associée en fonction du nombre d’ex-aequos */
exaequo = 1 / fact(exaequo + 1);
concord = 0;
end;
c + concord; /* On incrémente le nombre de vecteurs concordants retenus */
d + discord; /* On incrémente le nombre de vecteurs disconcordants retenus */
t + exaequo; /* On incrémente le nombre de vecteurs ex-aequos retenus */
tot + 1;
%do z = 1 %to &nb_classe.;
end;
%end;
run;

data kuplets;
set kuplets;
bilan = (c + t) / tot; /* Calcul du score HUM */
call symputx(« hum »,bilan);
run;

data hum;
set hum;
if i = &it_combi. then hum = &hum.;
run;

%end;

/* 8) Optimisation */
proc sql noprint;
create table hum2 as select *, max(hum) as top from hum; /* Calcul du max de HUM */
quit;

data hum; /* Conservation uniquement de la ligne maximisant le HUM en fonction du versus de combinaisons de classes */
set hum2;
if hum = top then selection = « * »;
run;

/* 9) Edition finale de la table */
data hum;
set hum;
%do it_cl = 1 %TO &nb_classe.;
%do it_label = 1 %TO &nb_classe.;
if x&it_cl. = &it_label. then xc&it_cl. = « %scan(&liste_classe_O.,&it_label.) »; /* Remise des labels originaux des classes */
%end;
%end;
ref = 1 / &nb_classe.; /* Insertion de la référence aléatoire */
effectif = compress(« &effectif. »); /* Insertion des effectifs */
%do it_label = 1 %TO &nb_classe.;
rename xc&it_label. = Classe_&it_label.; /* Rename des colonnes de la table pour plus de lisibilité */
%end;
drop top i;
%do it_cl = 1 %TO &nb_classe.;
drop x&it_cl.;
%end;
run;

%end;

%end;

/* 10) Nettoyage de la work */
proc datasets lib = work nolist;
delete contents norm recap_classe versus hum2 kuplets format;
run;

/* 11) Reset des options */
options notes nospool;

%mend;

Fonction R:

HUM = function(DATA) {

# Ce programme calcule l’hypervolume sous la courbe ROC (HUM). La structure attendue du jeu de données DATA est: 2
# colonnes – la première doit contenir la variable réponse Y (format numérique ou caractère) / la seconde doit contenir
# la variable explicative (format numérique). Il fournit en sortie les différents scores de HUM pour les différentes
# combinaisons de classe de la variable réponse (versus). Remerciement au forumeur Sengar pour son aide non négligeable
# sur le forum developpez.net

library(combinat) # Chargement du package permettant d’effectuer les permutations de classes de Y

TABLE = na.omit(DATA) # Suppression des données manquantes
biblio = summary(as.factor(TABLE[,1])) # Récupération des différents classes
nb_class = length(biblio)
Class = names(biblio)

permut = permn(Class) # Calcul des différentes permutations de classes
nb_permut = factorial(nb_class) # Calcul de la référence
tab_cas = matrix(0,nb_permut,nb_class)
for (k in 1:factorial(nb_class)) {tab_cas[k,] = permut[[k]]} # Création de la table des permutations des classes

nobs_tot = dim(TABLE)[1]
# Récupération du nombre d’effectif par classe
for (k in 1:nb_class) {
nobs = dim(TABLE[which(TABLE[,1] == Class[k]),])[1]
if (k == 1) {effectif = paste(« [« ,nobs, »(« ,round(100*nobs/nobs_tot), »% ) »)}
if ((k > 1) && (k < nb_class)) {effectif = paste(effectif, »+ »,nobs, »(« ,round(100*nobs/nobs_tot), »% ) »)}
if (k == nb_class) {effectif = paste(effectif, »+ »,nobs, »(« ,round(100*nobs/nobs_tot), »% ) ] »)}
}

tab_HUM = cbind(tab_cas,effectif,1/nb_class,matrix(0,nb_permut,2)) # Création de la table finale contenant le score de HUM pour chaque permutation de classes

# Calcul du HUM pour les différentes permutations de classes
for (k in 1:nb_permut) {
# Création de la table ordonnée en fonction de l’ordre de la permutation étudiée
for (l in 1:nb_class) {
if (l == 1) {TABLE_c = cbind(l,TABLE[which(TABLE[,1] == tab_HUM[k,l]),2])}
if (l > 1) {TABLE_c = rbind(TABLE_c,cbind(l,TABLE[which(TABLE[,1] == tab_HUM[k,l]),2]))}
}
l = list()
# Création de la table des différents k-uplets
for (c in 1:nb_class) {l[[c]] = TABLE_c[,2][TABLE_c[,1] == c]}
kuplet = expand.grid(l)
X = TABLE_c[,2]
kuplet$res = apply(kuplet,1,function(X) {as.logical(min(order(X) == 1:length(X)))}) # Comparaison des différents k-uplet
kuplet$score = 0 # Variable contenant le score finale directement issue de kuplet$res et corrigé en fonction des ex-aequos non pris en compte encore
ind.verif = which(kuplet$res == TRUE) # Récupération des cas où le k-uplet est ordonné, à noter qu’en cas d’ex-aequos le k-uplet est supposé ordonné à tort pour le moment
kuplet$score[ind.verif] = 1 # Au départ, tous les k-uplets ordonnés sont supposés avoir le score maximal
nb.verif = length(ind.verif)
# Correction du score en fonction des ex-aequos
for (i_verif in 1:nb.verif) {
exaequos = 0
for (j in 2:nb_class) {if (kuplet[ind.verif[i_verif],j-1] == kuplet[ind.verif[i_verif],j]) {exaequos = exaequos + 1}} # Si le terme j-1 est égal au terme j alors on a bien un ex-aequos et on incrémente le score
kuplet$score[ind.verif[i_verif]] = 1/factorial(exaequos+1) # Calcul du score final
}
tab_HUM[k,nb_class + 3] = sum(kuplet$score)/dim(kuplet)[1] # Calcul du HUM après avoir corrigé en fonction des ex-aequos
}
# Détermination de la permutation maximisant le HUM
MAX = which((tab_HUM[,nb_class + 3] == max(tab_HUM[,nb_class + 3])) == TRUE)
tab_HUM[MAX,nb_class + 4] = « * »
tab_HUM = as.data.frame(tab_HUM)
# Labellisation de la table de synthèse des différents HUM en fonction des différentes permutations de classes
names(tab_HUM) = c(rep(« Classe »,nb_class), »Effectif », »Reference », »HUM », »BEST »)

}

\bullet Bibliographie: 

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

– Signal detection theory and psychophysics de David M. Green et John A. Swets

– ROC analysis with multiple classes and multiple tests: methodology and its application in microarray studies, J. Li et J. P. Fine

– Three-way ROCs, D. Mossman

– Ordered multiple-class ROC analysis with continuous measurements, C. T. Nakas et C. T. Yiannoutsos

– The meaning and use of the volume under a three class ROC surface (VUS), X. Hee et E. C. Frey

– Volume under the ROC surface for multi-class problems, C. Ferri, J. Hernandez-Orallo et M. A. Salido

– Efficient multiclass ROC approximation by decomposition via confusion matrix perturbation analysis, T. C. W. Langrebe et R. P. W. Duin

– The hypervolume under the ROC hypersurface of « Near-Guessing » and « Near-Perfect » observers in N-Class classification tasks, D. C. Edwards, C. E. Metz et R. M. Nishikawa