Je suis intéressé par le calcul de l'aire sous la courbe (AUC), ou la statistique C, à la main pour un modèle de régression logistique binaire.

Par exemple, dans le jeu de données de validation, j'ai la valeur vraie pour la variable dépendante, rétention (1 = retenue; 0 = non conservée), ainsi qu'un statut de rétention prévu pour chaque observation générée par mon analyse de régression à l'aide d'un modèle construit à l’aide du kit d’entraînement (celui-ci va de 0 à 1).

Mes idées initiales étaient d'identifier le nombre "correct" de classifications de modèles et de diviser simplement le nombre d'observations "correctes" par le nombre total d'observations pour calculer la statistique c. Par "correct", si le statut de rétention réel d'une observation = 1 et le statut de rétention prévu est> 0,5, il s'agit d'une classification "correcte". De plus, si le statut de rétention réel d'une observation = 0 et le statut de rétention prévu est <0,5, il s'agit également d'une classification "correcte". Je suppose qu'un "lien" se produirait lorsque la valeur prédite = 0,5, mais ce phénomène ne se produit pas dans mon jeu de données de validation. D'autre part, les classifications "incorrectes" seraient si le statut de rétention réel d'une observation = 1 et le statut de rétention prédit est <0. 5 ou si le statut de rétention réel pour un résultat = 0 et si le statut de rétention prévu est> 0,5. Je connais TP, FP, FN, TN, mais je ne sais pas comment calculer la statistique c étant donné ces informations.

Je recommanderais l'article de 1982 de Hanley & McNeil, intitulé « La signification et l'utilisation de la courbe de la zone sous une caractéristique de fonctionnement du récepteur ».

Exemple

Ils présentent le tableau suivant de l'état de la maladie et du résultat du test (correspondant par exemple au risque estimé à partir d'un modèle logistique). Le premier chiffre à droite est le nombre de patients avec le statut réel de la maladie «normal» et le deuxième nombre est le nombre de patients avec le statut réel de la maladie «anormal»:

(1) Absolument normal: 33/3
(2) Probablement normal: 6/2
(3) Incertain: 6/2
(4) Probablement anormal: 11/11
(5) Absolument anormal: 2/33

Il y a donc au total 58 patients «normaux» et «51» anormaux. Nous voyons que lorsque le prédicteur est 1, «tout à fait normal», le patient est généralement normal (vrai pour 33 des 36 patients) et lorsqu'il est 5, «tout à fait anormal», le patient est habituellement anormal (vrai pour 33 des 35 patients), le prédicteur a donc du sens. Mais comment devrions-nous juger un patient avec un score de 2, 3 ou 4? Ce que nous définissons comme seuil pour juger un patient comme anormal ou normal détermine la sensibilité et la spécificité du test résultant.

Sensibilité et spécificité

Nous pouvons calculer la sensibilité et la spécificité estimées pour différentes valeurs limites. (J'écrirai simplement 'sensibilité' et 'spécificité' à partir de maintenant, laissant la nature estimée des valeurs être implicite.)

Si nous choisissons notre seuil pour classer tous les patients comme anormaux, peu importe les résultats de leurs tests (c’est-à-dire que nous choisissons le seuil 1+), nous obtiendrons une sensibilité de 51/51 = 1. La spécificité sera 0 / 58 = 0. Cela ne semble pas si bon.

OK, choisissons donc un seuil moins strict. Nous classons les patients comme anormaux uniquement si leur résultat de test est égal ou supérieur à 2. Nous manquons alors 3 patients anormaux et avons une sensibilité de 48/51 = 0,94. Mais nous avons une spécificité beaucoup plus grande, de 33/58 = 0,57.

Nous pouvons maintenant continuer en choisissant diverses valeurs de coupure (3, 4, 5,> 5). (Dans le dernier cas, nous ne classerons aucun patient comme anormal, même si son score au test est le plus élevé possible: 5).

La courbe ROC

Si nous faisons cela pour tous les seuils possibles et que nous traçons la sensibilité par rapport à 1 moins la spécificité, nous obtenons la courbe ROC. Nous pouvons utiliser le code R suivant:

# Data
norm     = rep(1:5, times=c(33,6,6,11,2))
abnorm   = rep(1:5, times=c(3,2,2,11,33))
testres  = c(abnorm,norm)
truestat = c(rep(1,length(abnorm)), rep(0,length(norm)))

# Summary table (Table I in the paper)
( tab=as.matrix(table(truestat, testres)) )

La sortie est:

        testres
truestat  1  2  3  4  5
       0 33  6  6 11  2
       1  3  2  2 11 33

Nous pouvons calculer diverses statistiques:

( tot=colSums(tab) )                            # Number of patients w/ each test result
( truepos=unname(rev(cumsum(rev(tab[2,])))) )   # Number of true positives
( falsepos=unname(rev(cumsum(rev(tab[1,])))) )  # Number of false positives
( totpos=sum(tab[2,]) )                         # The total number of positives (one number)
( totneg=sum(tab[1,]) )                         # The total number of negatives (one number)
(sens=truepos/totpos)                           # Sensitivity (fraction true positives)
(omspec=falsepos/totneg)                        # 1 − specificity (false positives)
sens=c(sens,0); omspec=c(omspec,0)              # Numbers when we classify all as normal

Et en utilisant cela, nous pouvons tracer la courbe ROC (estimée):

plot(omspec, sens, type="b", xlim=c(0,1), ylim=c(0,1), lwd=2,
     xlab="1 − specificity", ylab="Sensitivity") # perhaps with xaxs="i"
grid()
abline(0,1, col="red", lty=2)

Calculer manuellement la AUC

Nous pouvons très facilement calculer l'aire sous la courbe ROC, en utilisant la formule de l'aire d'un trapèze:

height = (sens[-1]+sens[-length(sens)])/2
width = -diff(omspec) # = diff(rev(omspec))
sum(height*width)

Le résultat est 0.8931711.

Une mesure de concordance

La CUA peut également être considérée comme une mesure de concordance. Si nous prenons toutes les paires possibles de patients où l’une est normale et l’autre est anormale, nous pouvons calculer la fréquence à laquelle le résultat anormal du test est le plus élevé (le plus «anormal») (si elles ont la même valeur, nous comptez comme une demi-victoire):

o = outer(abnorm, norm, "-")
mean((o>0) + .5*(o==0))

La réponse est à nouveau 0,8931711, l'aire sous la courbe ROC. Ce sera toujours le cas.

Une vue graphique de la concordance

Comme Harrell l’a souligné dans sa réponse, il existe également une interprétation graphique. Représentons le score du test (estimation du risque) sur l’ axe y et le statut réel de la maladie sur l’ axe x (ici, avec quelques tremblements, pour montrer les points qui se chevauchent):

plot(jitter(truestat,.2), jitter(testres,.8), las=1,
     xlab="True disease status", ylab="Test score")

Tracons maintenant une ligne entre chaque point à gauche (un patient «normal») et chaque point à droite (un patient «anormal»). La proportion de lignes présentant une pente positive (c'est-à-dire la proportion de paires concordantes ) est l'indice de concordance (les lignes plates comptent pour une "concordance de 50%").

Il est un peu difficile de visualiser les lignes réelles pour cet exemple, en raison du nombre de liens (score de risque égal), mais avec un peu de frémissement et de transparence, nous pouvons obtenir un graphique raisonnable:

d = cbind(x_norm=0, x_abnorm=1, expand.grid(y_norm=norm, y_abnorm=abnorm))
library(ggplot2)
ggplot(d, aes(x=x_norm, xend=x_abnorm, y=y_norm, yend=y_abnorm)) +
  geom_segment(colour="#ff000006",
               position=position_jitter(width=0, height=.1)) +
  xlab("True disease status") + ylab("Test\nscore") +
  theme_light()  + theme(axis.title.y=element_text(angle=0))

Nous voyons que la plupart des lignes ont une pente ascendante, donc l'indice de concordance sera élevé. Nous voyons également la contribution à l'indice de chaque type de paire d'observateurs. La majeure partie provient de patients normaux présentant un indice de risque de 1 associé à des patients anormaux présentant un indice de risque de 5 (1 à 5 paires), mais une grande partie provient également de 1 à 4 paires et de 4 à 5 paires. Et il est très facile de calculer l’indice de concordance réel en fonction de la définition de la pente:

d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm))
mean((d$slope > 0) + .5*(d$slope==0))

La réponse est à nouveau 0.8931711, c’est-à-dire l’ASC.

Le test de Wilcoxon – Mann – Whitney

Il existe un lien étroit entre la mesure de concordance et le test de Wilcoxon – Mann – Whitney. En réalité, ce dernier teste si la probabilité de concordance (c'est-à-dire que c'est le patient anormal dans une paire aléatoire normale anormale qui aura le résultat de test le plus «anormal») est exactement 0,5. Et sa statistique de test n’est qu’une simple transformation de la probabilité de concordance estimée:

> ( wi = wilcox.test(abnorm,norm) )
    Wilcoxon rank sum test with continuity correction

data:  abnorm and norm
W = 2642, p-value = 1.944e-13
alternative hypothesis: true location shift is not equal to 0

La statistique de test ( W = 2642) compte le nombre de paires concordantes. Si on le divise par le nombre de paires possibles, on obtient un nombre familier:

w = wi$statistic
w/(length(abnorm)*length(norm))

Oui, il s’agit de 0.8931711, l’aire sous la courbe ROC.

Des moyens plus simples pour calculer l'AUC (en R)

Mais rendons la vie plus facile pour nous-mêmes. Différents forfaits calculent automatiquement l'AUC pour nous.

Le forfait Epi

Le Epipaquet crée une belle courbe ROC avec diverses statistiques (y compris la AUC) intégrée:

library(Epi)
ROC(testres, truestat) # also try adding plot="sp"

Le paquet pROC

J'aime aussi le pROCpaquet, car il peut lisser l'estimation du RDC (et calculer une estimation de l'ASC basée sur le ROC lissé):

(La ligne rouge correspond au ROC d'origine et la ligne noire au ROC lissé. Notez également le format d'image par défaut 1: 1. Il est judicieux de l'utiliser, car la sensibilité et la spécificité ont une plage de 0 à 1.)

L'ASC estimée à partir du RDC lissé est de 0,9107, similaire à l' ASC du ROC non lissé , mais légèrement plus grande (si vous regardez la figure, vous pouvez facilement voir pourquoi elle est plus grande). (Bien que nous ayons vraiment trop peu de valeurs de résultat de test distinctes possibles pour calculer une AUC lisse).

Le paquetage rms

Le rmspackage de Harrell peut calculer diverses statistiques de concordance connexes à l'aide de la rcorr.cens()fonction. Le C Indexdans sa sortie est l'AUC:

> library(rms)
> rcorr.cens(testres,truestat)[1]
  C Index 
0.8931711

Le package caTools

Enfin, nous avons le caToolspaquet et sa colAUC()fonction. Il présente quelques avantages par rapport aux autres packages (principalement la vitesse et la capacité de travailler avec des données multidimensionnelles - voir ?colAUC) qui peuvent parfois être utiles. Mais bien sûr, cela donne la même réponse que nous avons calculée encore et encore:

library(caTools)
colAUC(testres, truestat, plotROC=TRUE)
             [,1]
0 vs. 1 0.8931711

Mots finaux

Beaucoup de gens semblent penser que l'AUC nous dit à quel point un test est «bon». Et certaines personnes pensent que l'AUC est la probabilité que le test classe correctement un patient. Ce n'est pas . Comme vous pouvez le constater à l'aide de l'exemple et des calculs ci-dessus, l'AUC nous dit quelque chose à propos d'une famille de tests, un test pour chaque seuil possible.

Et la AUC est calculée sur la base de seuils inutilisés dans la pratique. Pourquoi devrions-nous nous préoccuper de la sensibilité et de la spécificité des valeurs de seuil "sans sens"? Pourtant, c'est sur quoi l'AUC est (partiellement) basée. (Bien sûr, si l'AUC est très proche de 1, presque tous les tests possibles auront un grand pouvoir discriminant, et nous en serions tous très heureux.)

L’interprétation des paires «aléatoire-normal-anormal» de la séquence AUC est agréable (et peut être étendue, par exemple, aux modèles de survie, où nous voyons si c’est la personne présentant le risque le plus élevé (relatif) qui décède le plus tôt). Mais on ne l'utilisera jamais dans la pratique. C'est un cas rare où l'on sait qu'on a une personne en bonne santé et une personne malade, on ne sait pas quelle personne est la personne malade et il faut décider laquelle d'entre elles doit traiter. (Dans tous les cas, la décision est facile; traitez celui qui présente le risque estimé le plus élevé.)

Je pense donc que l’étude de la courbe ROC réelle sera plus utile que la simple analyse de la mesure de synthèse de la CUA. Et si vous utilisez le ROC avec (les estimations des) coûts des faux positifs et des faux négatifs, ainsi que les taux de base de ce que vous étudiez, vous pouvez aller quelque part.

Notez également que l'AUC ne mesure que la discrimination , pas l'étalonnage. Autrement dit, il mesure si vous pouvez faire la distinction entre deux personnes (une malade et une en bonne santé), en fonction du score de risque. Pour cela, il ne considère que les valeurs de risque relatives (ou les rangs, si vous préférez, voir l'interprétation du test de Wilcoxon-Mann-Whitney), et non les valeurs absolues, qui devraient vous intéresser. Par exemple, si vous divisez chaque risque Estimez votre modèle logistique par 2, vous obtiendrez exactement la même AUC (et le ROC).

Lors de l'évaluation d'un modèle de risque, l' étalonnage est également très important. Pour examiner cela, vous examinerez tous les patients avec un score de risque voisin de 0,7 (par exemple) et verrez si environ 70% d'entre eux étaient réellement malades. Faites ceci pour chaque score de risque possible (en utilisant éventuellement une sorte de régression de lissage / locale). Tracez les résultats et vous obtiendrez une mesure graphique de la calibration .

Si vous avez un modèle avec à la fois un bon calibrage et une bonne discrimination, vous commencez à avoir un bon modèle. :)

Regardez cette question: Comprendre la courbe ROC

Voici comment construire une courbe ROC (à partir de cette question):

Dessin de la courbe ROC

étant donné un ensemble de données traité par votre classificateur de classement

• classer les exemples de tests sur un score décroissant
• commencer dans$(0, 0)$
• pour chaque exemple (dans l'ordre décroissant) $x$
  • si est positif, déplacez vers le haut1 /$x$$1/\text{pos}$
  • si est négatif, déplacez droite1 /$x$$1/\text{neg}$

où et sont les fractions des exemples positifs et négatifs respectivement.$\text{pos}$$\text{neg}$

Vous pouvez utiliser cette idée pour calculer manuellement le ROC AUC en utilisant l'algorithme suivant:

auc = 0.0
height = 0.0

for each training example x_i, y_i
  if y_i = 1.0:
    height = height + tpr
  else 
    auc = auc + height * fpr

return auc

Cette belle image animée gif devrait illustrer ce processus plus clairement

Le message de Karl contient beaucoup d'excellentes informations. Mais, au cours des 20 dernières années, je n’ai pas encore vu un exemple de courbe ROC qui aurait changé l’attitude de quiconque dans la bonne direction. À mon humble avis, la seule valeur d'une courbe ROC est que son aire correspond à une probabilité de concordance très utile. La courbe ROC elle-même incite le lecteur à utiliser des seuils, ce qui est une mauvaise pratique statistique.

En ce qui concerne le calcul manuel de l’ index- , faites un tracé avec sur l’ axe et le prédicteur continu ou la probabilité prédite que sur l’ axe- . Si vous connectez chaque point avec avec chaque point avec , la proportion des lignes ayant une pente positive est la probabilité de concordance.Y = 0 , 1 x Y = 1 y Y = 0 Y = $c$$Y=0,1$$x$$Y=1$$y$$Y=0$$Y=1$

Toute mesure ayant un dénominateur dans ce paramètre constitue une règle de notation de l'exactitude inappropriée et doit être évitée. Cela inclut la proportion correctement classée, la sensibilité et la spécificité.$n$

Pour la fonction de Hmiscpackage R rcorr.cens, imprimez le résultat complet pour obtenir plus d'informations, notamment une erreur standard.

Voici une alternative à la manière naturelle de calculer l'ASC en utilisant simplement la règle trapézoïdale pour obtenir l'aire sous la courbe ROC.

L'ASC est égale à la probabilité qu'une observation positive échantillonnée de manière aléatoire ait une probabilité prédite (d'être positive) supérieure à une observation négative de ce type échantillonnée de manière aléatoire. Vous pouvez l’utiliser pour calculer l’ASC assez facilement dans n’importe quel langage de programmation en passant en revue toutes les combinaisons par paires d’observations positives et négatives. Vous pouvez également échantillonner des observations au hasard si la taille de l'échantillon est trop grande. Si vous souhaitez calculer l'ASC à l'aide d'un stylo et du papier, cette approche n'est peut-être pas la meilleure, sauf si vous avez un très petit échantillon / beaucoup de temps. Par exemple dans R:

n <- 100L

x1 <- rnorm(n, 2.0, 0.5)
x2 <- rnorm(n, -1.0, 2)
y <- rbinom(n, 1L, plogis(-0.4 + 0.5 * x1 + 0.1 * x2))

mod <- glm(y ~ x1 + x2, "binomial")

probs <- predict(mod, type = "response")

combinations <- expand.grid(positiveProbs = probs[y == 1L], 
        negativeProbs = probs[y == 0L])

mean(combinations$positiveProbs > combinations$negativeProbs)
[1] 0.628723

Nous pouvons vérifier en utilisant le pROCpackage:

library(pROC)
auc(y, probs)
Area under the curve: 0.6287

En utilisant un échantillonnage aléatoire:

mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE))
[1] 0.62896

1. Vous avez une valeur réelle pour les observations.
2. Calculer la probabilité a posteriori, puis classer les observations en fonction de cette probabilité.
3. En supposant que la probabilité de coupure de et le nombre d'observations :N Somme des vrais rangs - 0.5 P N ( P N + 1 $P$$N$
   $$\frac{\text{Sum of true ranks}-0.5PN(PN+1)}{PN(N-PN)}$$