I je suis intéressé par le calcul manuel de laire sous la courbe (AUC), ou la statistique c, pour un modèle de régression logistique binaire.
Par exemple, dans lensemble de données de validation, jai la vraie valeur de la variable dépendante, rétention (1 = conservé; 0 = non conservé), ainsi quun état de rétention prévu pour chaque observation générée par mon analyse de régression à laide dun modèle qui a été construit à laide de la formation set (cela va de 0 à 1).
Mes premières réflexions étaient didentifier le nombre « correct » de classifications de modèles et de diviser simplement le nombre dobservations « correctes » par le nombre dobservations totales à calculer la statistique c. Par « correct », si le véritable statut de rétention dune observation = 1 et le statut de rétention prévu est> 0,5, alors cest une classification « correcte ». De plus, si le statut de rétention réel dune observation = 0 et le statut de rétention prévu est < 0,5, il sagit également dune classification « correcte ». Je suppose quune « égalité » se produirait lorsque la valeur prédite = 0,5, mais ce phénomène ne se produit pas dans mon ensemble de données de validation. En revanche, les classifications « incorrectes » seraient si le véritable statut de rétention dune observation = 1 et le statut de rétention prévu est < 0,5 ou si le véritable statut de rétention dun résultat = 0 et létat de rétention prévu est> 0,5. Je connais TP, FP, FN, TN, mais je ne sais pas comment calculer la statistique c compte tenu de ces informations.
Réponse
Je recommanderais & larticle de 1982 de McNeil La signification et lutilisation de la zone sous une caractéristique de fonctionnement du récepteur (ROC ) courbe .
Exemple
Ils ont le tableau suivant de létat de la maladie et du résultat du test (correspondant, par exemple, au risque estimé à partir dun modèle logistique). Le premier chiffre à droite est le nombre de patients avec un statut de maladie vrai « normal » et le deuxième nombre est le nombre de patients avec un statut de maladie vrai « anormal »:
(1) Tout à fait normal: 33/3
(2) Probablement normal: 6/2
(3) Susceptible: 6/2
(4) Probablement anormal: 11/11
(5) Certainement 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 lorsquil est 5, « Certainement anormal , le patient est généralement anormal (vrai pour 33 des 35 patients), donc le prédicteur a du sens. Mais comment juger un patient avec un score de 2, 3 ou 4? Ce que nous définissons notre 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 estimée sensibilité et spécificité pour différents seuils. (Jécrirai simplement «sensibilité» et «spécificité» à partir de maintenant, laissant la nature estimée des valeurs être implicite.)
Si nous choisissons notre seuil de sorte que nous classions tout les patients comme anormaux, peu importe ce que disent les résultats de leurs tests (cest-à-dire que nous choisissons le seuil 1+), nous obtiendrons une sensibilité de 51/51 = 1. La spécificité sera 0/58 = 0. ça sonne si bien.
OK, alors choisissons une coupure moins stricte. Nous classons les patients comme anormaux uniquement sils ont un résultat de test é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 cela en choisissant différents seuils (3, 4, 5,> 5). (Dans ce dernier cas, nous ne classerons aucun patient comme anormal, même s’ils ont le score le plus élevé possible au test de 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)) )
Le résultat 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 ceci, 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)
Calcul manuel de la AUC
Nous pouvons très facilement calculer laire sous la courbe ROC, en utilisant la formule de laire dun 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
LAUC peut également être considérée comme une mesure de concordance.Si nous prenons toutes les paires possibles de patients dont lun est normal et lautre anormal, nous pouvons calculer à quelle fréquence cest lanormal qui a le résultat de test le plus élevé (le plus « anormal ) (si ils ont la même valeur, nous comptons cela comme une demi-victoire):
o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))
La réponse est à nouveau 0.8931711, la zone sous la courbe ROC. Ce sera toujours le cas.
Une vue graphique de la concordance
Comme la souligné Harrell dans sa réponse, cela a aussi une interprétation graphique. Tracons le score du test (estimation du risque) sur laxe y et létat réel de la maladie sur laxe x (ici avec quelques instabilités, pour montrer les points qui se chevauchent):
plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")
Trouvons maintenant une ligne entre chaque point à gauche (un patient «normal») et chaque point à droite (un patient «anormal»). La proportion de lignes avec une pente positive (cest-à-dire la proportion de paires concordantes ) est lindice de concordance (les lignes plates comptent pour «50% de concordance»).
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 nervosité et de transparence, nous pouvons obtenir un tracé 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 sont inclinées vers le haut, donc lindice de concordance sera élevé. On voit également la contribution à lindice de chaque type de couple dobservation. La plupart proviennent de patients normaux avec un score de risque de 1 et des patients anormaux avec un score 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 lindice de concordance réel basé sur 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, soit lAUC.
Le test de Wilcoxon – Mann – Whitney
Il existe un lien étroit entre la mesure de concordance et le Wilcoxon – Mann – Whitney test. En fait, ce dernier teste si la probabilité de concordance (cest-à-dire que cest le patient anormal dans une paire aléatoire normal – anormal qui aura le résultat de test le plus «anormal») est exactement de 0,5. Et sa statistique de test nest quune 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 nous le divisons par le nombre de paires possibles, nous obtenons un numéro familier:
w = wi$statistic w/(length(abnorm)*length(norm))
Oui, cest 0.8931711, laire sous la courbe ROC.
Des moyens plus simples de calculer lAUC (en R)
Mais facilitons-nous la vie. Il existe différents packages qui calculent automatiquement lAUC pour nous.
Le package Epi
Le package Epi
crée une belle courbe ROC avec divers statistiques (y compris lAUC) intégrées:
library(Epi) ROC(testres, truestat) # also try adding plot="sp"
Le package pROC
Jaime aussi le package pROC
, car il peut lisser lestimation ROC (et calculer une estimation AUC basée sur le ROC lissé):
(La ligne rouge est le ROC dorigine, et la ligne noire est le ROC lissé. Notez également le rapport hauteur / largeur par défaut 1: 1. Il est logique de lutiliser, car la sensibilité et la spécificité ont un 0–1 .)
LAUC estimée à partir du ROC lissé est de 0,9107, similaire à, mais légèrement plus grande que, lAUC du ROC non lissé (si vous regardez t le chiffre, vous pouvez facilement voir pourquoi il est plus grand). (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 paquet rms
de Harrell peut calculer diverses statistiques de concordance associées en utilisant la fonction rcorr.cens()
. Le C Index
dans sa sortie est lAUC:
> library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711
Le package caTools
Enfin, nous avons le package caTools
et sa fonction colAUC()
. Il présente quelques avantages par rapport aux autres packages (principalement la vitesse et la possibilité 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 celle que nous avons calculée à maintes reprises:
library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711
Derniers mots
Beaucoup de gens semblent penser que lAUC nous dit à quel point bien un test est. Et certaines personnes pensent que lASC est la probabilité que le test classifie correctement un patient. Il sagit de et non de . Comme vous pouvez le voir dans lexemple et les calculs ci-dessus, lAUC nous dit quelque chose sur une famille de tests, un test pour chaque seuil possible.
Et lAUC est calculée sur la base des seuils que lon nutiliserait jamais dans la pratique. Pourquoi devrions-nous nous soucier de la sensibilité et de la spécificité des valeurs seuils «absurdes»? Pourtant, cest sur quoi repose (partiellement) la CUA. (Bien sûr, si lAUC est très proche de 1, presque tous les tests possibles auront un grand pouvoir discriminant, et nous serions tous très heureux.)
La normale aléatoire –Une interprétation anormale des paires de lASC est agréable (et peut être étendue, par exemple aux modèles de survie, où nous voyons si cest la personne avec le risque (relatif) le plus élevé qui meurt le plus tôt). Mais on ne lutiliserait jamais dans la pratique. Cest un cas rare où lon sait que lon a une personne en bonne santé et une malade, ne sait pas quelle personne est la personne malade et doit décidez lequel dentre eux traiter. (Dans tous les cas, la décision est facile; traitez celui qui présente le risque estimé le plus élevé.)
Je pense donc quétudier la courbe ROC réelle sera plus utile que de simplement regarder la mesure récapitulative de lASC. 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 arriver quelque part.
Notez également que lAUC 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 regarde que les valeurs de risque relatives (ou les classements, si vous voulez, voir l’interprétation du test de Wilcoxon – Mann – Whitney), pas les valeurs absolues, ce que vous devriez soyez intéressé. Par exemple, si vous divisez chaque estimation de risque de votre modèle logistique par 2, vous obtiendrez exactement la même AUC (et ROC).
Lors de lévaluation dun modèle de risque, létalonnage est également très important. Pour examiner cela, vous examinerez tous les patients avec un score de risque denviron 0,7, par exemple, et vous verrez si environ 70% dentre eux étaient réellement malades. Faites cela pour chaque score de risque possible (éventuellement en utilisant une sorte de lissage / régression locale). Tracez les résultats et vous obtiendrez une mesure graphique de calibrage .
Si vous avez un modèle avec les deux un bon calibrage et une bonne discrimination, alors vous commencer à avoir un bon modèle. 🙂
Commentaires
Réponse
Jetez un œil à cette question: Comprendre la courbe ROC
Voici comment créer une courbe ROC (à partir de cette question):
Dessiner une courbe ROC
étant donné un ensemble de données traité par votre classificateur de classement
- exemples de test de classement sur score décroissant
- commencent par $ (0, 0) $
- pour chaque exemple $ x $ (dans le ordre décroissant)
- si $ x $ est positif, déplacez $ 1 / \ text {pos} $ vers le haut
- si $ x $ est négatif, déplacez $ 1 / \ text {neg} $ vers la droite
où $ \ text {pos} $ et $ \ text {neg} $ sont respectivement les fractions dexemples positifs et négatifs.
Vous pouvez utiliser cette idée pour calculer manuellement AUC ROC en utilisant lalgorithme 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 par un gif devrait illustrer cela processus plus clair
Commentaires
- Merci @Alexey Grigorev, cest un excellent visuel et il se révélera probablement utile à lavenir! +1
- Pourriez-vous expliquer un peu les » fractions dexemples positifs et négatifs « , voulez-vous dire le plus petite valeur unitaire de deux axes?
- @Allan Ruin:
pos
signifie ici le nombre de données positives. Disons que vous avez 20 points de données, dont 11 points valent 1. Ainsi, lors du dessin du graphique, nous avons un rectangle 11×9 (hauteur x largeur). Alexey Grigorev a fait une mise à léchelle, mais laissez-le simplement tel quel ‘ si vous le souhaitez. Maintenant, déplacez simplement 1 sur le graphique à chaque étape.
Réponse
Le message de Karl a beaucoup mais je nai pas encore vu au cours des 20 dernières années un exemple de courbe ROC qui a changé la façon de penser de quiconque dans la bonne direction. À mon humble avis, la seule valeur dune courbe ROC est que sa surface 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.
Pour ce qui est de calculer manuellement lindex $ c $, faites un tracé avec $ Y = 0,1 $ sur le x $ -axis et le prédicteur continu ou la probabilité prédite que $ Y = 1 $ sur laxe $ y $. Si vous connectez chaque point avec $ Y = 0 $ avec chaque point avec $ Y = 1 $, la proportion des lignes qui ont une pente positive est la probabilité de concordance.
Toutes les mesures qui ont un dénominateur de $ n $ dans ce paramètre sont des règles de notation de précision incorrectes et doivent être évités. Cela inclut la proportion classée correctement, la sensibilité et la spécificité.
Pour la fonction R Hmisc
rcorr.cens
, imprimez la fonction tout le résultat pour voir plus dinformations, en particulier une erreur standard.
Commentaires
- Merci, @Frank Harell, japprécie votre point de vue. Jutilise simplement la statistique c comme une probabilité de concordance, car je ne ‘ que comme des seuils. Merci encore!
Réponse
Voici une alternative à la manière naturelle de calculer lAUC en utilisant simplement le trapèze règle pour obtenir laire sous la courbe ROC.
LAUC est égale à la probabilité quune observation positive échantillonnée au hasard ait une probabilité prédite (dêtre positive) supérieure à une observation négative échantillonnée au hasard. Vous pouvez lutiliser pour calculer lAUC assez facilement dans nimporte quel langage de programmation en passant par toutes les combinaisons par paires dobservations positives et négatives. Vous pouvez également échantillonner des observations de manière aléatoire si la taille de léchantillon est trop grande. Si vous souhaitez calculer lAUC en utilisant un stylo et du papier, ce nest peut-être pas la meilleure approche à moins que vous nayez 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 package pROC
:
library(pROC) auc(y, probs) Area under the curve: 0.6287
Utilisation de léchantillonnage aléatoire:
mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896
Réponse
- Vous avez une valeur vraie pour les observations.
- Calculez la probabilité postérieure, puis classez les observations en fonction de cette probabilité.
- En supposant une probabilité de coupure de $ P $ et le nombre dobservations $ N $:
$$ \ frac {\ text {Somme des vrais rangs} -0,5PN (PN + 1)} { PN (N-PN)} $$
Commentaires
- @ user73455 … 1) Oui, jai la vraie valeur pour les observations. 2) La probabilité postérieure est-elle synonyme de probabilités prédites pour chacune des observations? 3) compris; cependant, quest-ce que » Somme des vrais rangs » et comment calculer cette valeur? Peut-être quun exemple vous aiderait à expliquer cette réponse plus en détail? Merci!
sens=c(sens,0); omspec=c(omspec,0)
, ne devrait pas ‘ t ce soitsens=c(0, sens); omspec=c(0, omspec)
? Il trace correctement avec le premier0
mais pas comme il est actuellement dans la réponse.sens
etomspec
avant de tracer).