I. Ich bin daran interessiert, die Fläche unter der Kurve (AUC) oder die c-Statistik von Hand für ein binäres logistisches Regressionsmodell zu berechnen.
Zum Beispiel in Im Validierungsdatensatz habe ich den wahren Wert für die abhängige Variable Retention (1 = beibehalten; 0 = nicht beibehalten) sowie einen vorhergesagten Aufbewahrungsstatus für jede Beobachtung, die durch meine Regressionsanalyse unter Verwendung eines Modells generiert wurde, das mithilfe des Trainings erstellt wurde set (dies reicht von 0 bis 1).
Meine ersten Gedanken waren, die „richtige“ Anzahl von Modellklassifikationen zu identifizieren und einfach die Anzahl der „richtigen“ Beobachtungen durch die Anzahl der zu berechnenden Gesamtbeobachtungen zu dividieren die c-Statistik. Wenn der wahre Retentionsstatus einer Beobachtung = 1 ist und der vorhergesagte Retentionsstatus> 0,5 ist, ist dies eine „korrekte“ Klassifizierung. Wenn der wahre Aufbewahrungsstatus einer Beobachtung = 0 und der vorhergesagte Aufbewahrungsstatus < 0,5 ist, ist dies ebenfalls eine „korrekte“ Klassifizierung. Ich gehe davon aus, dass ein „Gleichstand“ auftreten würde, wenn der vorhergesagte Wert = 0,5 ist, aber dieses Phänomen tritt in meinem Validierungsdatensatz nicht auf. Andererseits wären „falsche“ Klassifizierungen, wenn der wahre Aufbewahrungsstatus einer Beobachtung = 1 und der vorhergesagte Aufbewahrungsstatus < 0,5 ist oder wenn der wahre Aufbewahrungsstatus für ein Ergebnis = 0 und der vorhergesagte Aufbewahrungsstatus ist> 0,5. Mir sind TP, FP, FN, TN bekannt, aber ich weiß nicht, wie die c-Statistik anhand dieser Informationen berechnet werden soll.
Antwort
Ich würde Hanleys & McNeils 1982er Artikel Die Bedeutung und Verwendung des Bereichs unter einer Empfängerbetriebscharakteristik (ROC) empfehlen ) Kurve .
Beispiel
Sie haben die folgende Tabelle mit dem Krankheitsstatus und dem Testergebnis (entsprechend beispielsweise dem geschätzten Risiko aus einem logistischen Modell). Die erste Zahl rechts ist die Anzahl der Patienten mit dem Status „wahr“ „normal“ und die zweite Zahl ist die Anzahl der Patienten mit dem Status „wahr“ „abnormal“:
(1) Definitiv normal: 33/3
(2) Wahrscheinlich normal: 6/2
(3) Fraglich: 6/2
(4) Wahrscheinlich abnormal: 11/11
(5) Definitiv abnormal: 2/33
Es gibt also insgesamt 58 „normale“ Patienten und „51“ abnormale. Wir sehen, dass wenn der Prädiktor 1 ist, „definitiv normal“, der Patient normalerweise normal ist (wahr für 33 der 36 Patienten), und wenn es 5 ist, „definitiv abnormal“, ist der Patient normalerweise abnormal (wahr für 33 der 36 Patienten) 35 Patienten), also macht der Prädiktor Sinn. Aber wie sollen wir einen Patienten mit einer Punktzahl von 2, 3 oder 4 beurteilen? Was wir unseren Grenzwert für die Beurteilung eines Patienten als abnormal oder normal festgelegt haben, um die Sensitivität und Spezifität des resultierenden Tests zu bestimmen.
Sensitivität und Spezifität
Wir können die geschätzte
Wenn wir unseren Cutoff so wählen, dass wir all die Patienten als abnormal, unabhängig davon, was ihre Testergebnisse aussagen (dh wir wählen den Cutoff 1+), erhalten wir eine Sensitivität von 51/51 = 1. Die Spezifität ist 0/58 = 0. Nicht klingt so gut.
OK, also wählen wir einen weniger strengen Cutoff. Wir klassifizieren Patienten nur dann als abnormal, wenn sie ein Testergebnis von 2 oder höher haben. Wir vermissen dann 3 abnormale Patienten und haben eine Empfindlichkeit von 48/51 = 0,94. Wir haben jedoch eine stark erhöhte Spezifität von 33/58 = 0,57.
Wir können dies nun fortsetzen und verschiedene Grenzwerte auswählen (3, 4, 5,> 5). (Im letzten Fall werden keine Patienten als abnormal eingestuft, selbst wenn sie die höchstmögliche Testnote von 5 haben.)
Die ROC-Kurve
Wenn wir dies für alle möglichen Grenzwerte tun und die Empfindlichkeit gegen 1 minus die Spezifität darstellen, erhalten wir die ROC-Kurve. Wir können den folgenden R-Code verwenden:
# 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)) )
Die Ausgabe lautet:
testres truestat 1 2 3 4 5 0 33 6 6 11 2 1 3 2 2 11 33
Wir können verschiedene Statistiken berechnen:
( 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
Und damit können wir die (geschätzte) ROC-Kurve zeichnen:
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)
Manuelles Berechnen der AUC
Wir können die Fläche unter der ROC-Kurve sehr einfach anhand der Formel für die Fläche eines Trapezes berechnen:
height = (sens[-1]+sens[-length(sens)])/2 width = -diff(omspec) # = diff(rev(omspec)) sum(height*width)
Das Ergebnis ist 0,8931711.
Ein Konkordanzmaß
Die AUC kann auch als Konkordanzmaß angesehen werden.Wenn wir alle möglichen Paare von Patienten nehmen, bei denen einer normal und der andere abnormal ist, können wir berechnen, wie häufig der abnormale das höchste (am meisten „abnormal aussehende“) Testergebnis aufweist (wenn Sie haben den gleichen Wert, wir zählen dies als einen halben Sieg):
o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))
Die Antwort lautet erneut 0,8931711, der Bereich unter der ROC-Kurve. Dies wird immer der Fall sein.
Eine grafische Ansicht der Konkordanz
Wie Harrell in seiner Antwort hervorhob, hat dies auch eine grafische Interpretation. Zeichnen wir den Testergebnis (Risikoschätzung) auf der y -Achse und den tatsächlichen Krankheitsstatus auf der x -Achse (hier mit etwas Jitter, um überlappende Punkte anzuzeigen):
plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")
Zeichnen wir nun eine Linie zwischen jedem Punkt links (ein „normaler“ Patient) und jedem Punkt rechts (ein „abnormaler“ Patient). Der Anteil von Linien mit einer positiven Steigung (d. H. Der Anteil von konkordanten Paaren) ist der Konkordanzindex (flache Linien zählen als „50% Konkordanz“).
Aufgrund der Anzahl der Bindungen (gleiches Risiko) ist es etwas schwierig, die tatsächlichen Linien für dieses Beispiel zu visualisieren, aber mit etwas Jittering und Transparenz können wir eine vernünftige Darstellung erhalten:
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))
Wir sehen, dass die meisten Linien nach oben geneigt sind, sodass der Konkordanzindex hoch ist. Wir sehen auch den Beitrag jeder Art von Beobachtungspaar zum Index. Das meiste davon kommt von normalen Patienten mit einem Risikowert von 1 gepaart mit abnormalen Patienten mit einem Risikowert von 5 (1–5 Paare), aber ziemlich viel kommt auch von 1–4 Paaren und 4–5 Paaren. Und es ist sehr einfach, den tatsächlichen Konkordanzindex basierend auf der Steigungsdefinition zu berechnen:
d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm)) mean((d$slope > 0) + .5*(d$slope==0))
Die Antwort lautet erneut 0,8931711, dh die AUC.
Der Wilcoxon-Mann-Whitney-Test
Es besteht ein enger Zusammenhang zwischen dem Konkordanzmaß und dem Wilcoxon-Mann-Whitney-Test Prüfung. Tatsächlich testet letzteres, ob die Wahrscheinlichkeit einer Übereinstimmung (d. H. Dass es der abnormale Patient in einem zufälligen normal-abnormalen Paar ist, das das am meisten „abnormal aussehende“ Testergebnis hat) genau 0,5 beträgt. Die Teststatistik ist nur eine einfache Transformation der geschätzten Konkordanzwahrscheinlichkeit:
> ( 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
Die Teststatistik (W = 2642
) zählt die Anzahl der übereinstimmenden Paare. Wenn wir es durch die Anzahl der möglichen Paare teilen, erhalten wir eine vertraute Zahl:
w = wi$statistic w/(length(abnorm)*length(norm))
Ja, es ist 0,8931711, der Bereich unter der ROC-Kurve.
Einfachere Methoden zur Berechnung der AUC (in R)
Aber machen wir uns das Leben leichter. Es gibt verschiedene Pakete, die die AUC für uns automatisch berechnen.
Das Epi-Paket
Das Epi
-Paket erstellt eine schöne ROC-Kurve mit verschiedenen Statistiken (einschließlich der AUC) eingebettet:
library(Epi) ROC(testres, truestat) # also try adding plot="sp"
Das pROC-Paket
Ich mag auch das pROC
-Paket, da es kann Glätten Sie die ROC-Schätzung (und berechnen Sie eine AUC-Schätzung basierend auf dem geglätteten ROC):
(Die rote Linie ist der ursprüngliche ROC und die schwarze Linie ist der geglättete ROC. Beachten Sie auch das Standard-Seitenverhältnis von 1: 1. Es ist sinnvoll, dies zu verwenden, da sowohl die Empfindlichkeit als auch die Spezifität 0–1 haben Bereich.)
Die geschätzte AUC aus dem geglätteten ROC beträgt 0,9107, ähnlich, aber etwas größer als die AUC aus dem ungeglätteten ROC (wenn Sie a In der Abbildung können Sie leicht erkennen, warum sie größer ist. (Obwohl wir wirklich zu wenige mögliche unterschiedliche Testergebniswerte haben, um eine glatte AUC zu berechnen).
Das Effektivpaket
Harrells rms
-Paket Mit der Funktion rcorr.cens()
können verschiedene zugehörige Konkordanzstatistiken berechnet werden. Die C Index
in ihrer Ausgabe ist die AUC:
> library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711
Das Paket caTools
Schließlich haben wir das Paket caTools
und dessen Funktion colAUC()
. Es hat einige Vorteile gegenüber anderen Paketen (hauptsächlich Geschwindigkeit und die Fähigkeit, mit mehrdimensionalen Daten zu arbeiten – siehe ?colAUC
), die manchmal hilfreich sein können.Aber natürlich gibt es die gleiche Antwort, die wir immer wieder berechnet haben:
library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711
Letzte Worte
Viele Leute scheinen zu glauben, dass die AUC uns sagt, wie gut ein Test ist. Und einige Leute denken, dass die AUC die Wahrscheinlichkeit ist, dass der Test einen Patienten korrekt klassifiziert. Es ist nicht . Wie Sie dem obigen Beispiel und den Berechnungen entnehmen können, sagt die AUC etwas über eine Familie von Tests aus, einen Test für jeden möglichen Grenzwert.
Und die AUC wird basierend auf berechnet Cutoffs würde man in der Praxis niemals verwenden. Warum sollten wir uns um die Sensitivität und Spezifität von „unsinnigen“ Grenzwerten kümmern? Darauf basiert die AUC (teilweise). (Wenn die AUC sehr nahe bei 1 liegt, hat fast jeder mögliche Test eine große Unterscheidungskraft, und wir würden uns alle sehr freuen.)
Die zufällige Normalität –Anormale Paarinterpretation der AUC ist nett (und kann zum Beispiel auf Überlebensmodelle ausgedehnt werden, bei denen wir sehen, ob es die Person mit der höchsten (relativen) Gefahr ist, die am frühesten stirbt). Aber man würde es niemals in der Praxis anwenden. Es ist ein seltener Fall, in dem man weiß , dass man eine gesunde und eine kranke Person hat, nicht weiß, welche Person die kranke ist und muss Entscheiden Sie, welche von ihnen behandelt werden sollen. (In jedem Fall ist die Entscheidung einfach; behandeln Sie die Entscheidung mit dem höchsten geschätzten Risiko.)
Daher denke ich, dass das Studium der tatsächlichen ROC-Kurve nützlicher ist als nur das Betrachten die AUC-Zusammenfassungsmaßnahme. Und wenn Sie den ROC zusammen mit (Schätzungen der) Kosten für falsch positive und falsch negative Ergebnisse zusammen mit den Basisraten für das, was Sie studieren, verwenden, können Sie irgendwohin gelangen.
Beachten Sie auch, dass die AUC nur die Diskriminierung misst, nicht die Kalibrierung. Das heißt, es wird anhand der Risikobewertung gemessen, ob Sie zwischen zwei Personen (einer kranken und einer gesunden) unterscheiden können. Hierzu werden nur relative Risikowerte betrachtet (oder Ränge, wenn Sie so wollen, vgl. Die Wilcoxon-Mann-Whitney-Testinterpretation), nicht die absoluten, die Sie sollten interessiert sein. Wenn Sie beispielsweise jede Risikoschätzung aus Ihrem Logistikmodell durch 2 teilen, erhalten Sie genau die gleiche AUC (und ROC).
Bei der Bewertung eines Risikomodells Die Kalibrierung ist ebenfalls sehr wichtig. Um dies zu untersuchen, werden Sie alle Patienten mit einem Risikowert von etwa 0,7 untersuchen und feststellen, ob ungefähr 70% davon tatsächlich krank waren. Tun Sie dies für jede mögliche Risikobewertung (möglicherweise unter Verwendung einer Art Glättung / lokaler Regression). Zeichnen Sie die Ergebnisse, und Sie erhalten ein grafisches Maß für die Kalibrierung .
Wenn Sie ein Modell mit sowohl guter Kalibrierung als auch guter Unterscheidung haben, sind Sie es fange an, ein gutes Modell zu haben. 🙂
Kommentare
- Vielen Dank, @Karl Ove Hufthammer, dies ist die gründlichste Antwort, die ich je erhalten habe. Ich schätze besonders Ihren Abschnitt “ Final Words „. Ausgezeichnete Arbeit! Nochmals vielen Dank!
- Vielen Dank für diese detaillierte Antwort. Ich arbeite mit einem Datensatz, in dem Epi :: ROC () v2.2.6 davon überzeugt ist, dass die AUC 1,62 beträgt (nein, es ist keine mentalistische Studie), aber laut ROC glaube ich viel mehr an 0,56, dass der obige Code resultiert in.
- Ich denke, es gibt einen kleinen Fehler in
sens=c(sens,0); omspec=c(omspec,0)
, sollte ‚ nichtsens=c(0, sens); omspec=c(0, omspec)
? Es wird korrekt mit dem führenden0
dargestellt, jedoch nicht so, wie es derzeit in der Antwort angegeben ist. - Nein, die aktuelle Definition lautet AFAICS, korrekt, @steveb, und führt zu einer korrekten Darstellung. Ich denke, was vielleicht verwirrend ist, ist, dass die ROC-Kurve von rechts nach links (dh von der oberen rechten Ecke zur unteren linken Ecke) gezeichnet wird, nicht von links nach links richtig , wie die meisten Handlungen. Dies ist nur das Ergebnis meiner Definition der Variablen. man hätte es genauso gut von links nach rechts zeichnen können (indem man sowohl
sens
als auchomspec
vor dem Zeichnen umkehrt).
Antwort
Sehen Sie sich diese Frage an: ROC-Kurve verstehen
So erstellen Sie eine ROC-Kurve (aus dieser Frage):
Zeichnen der ROC-Kurve
anhand eines von Ihrem verarbeiteten Datensatzes Ranking-Klassifikator
- Rank-Testbeispiele für abnehmende Punktzahl
- beginnen in $ (0, 0) $
- für jedes Beispiel $ x $ (in der absteigende Reihenfolge)
- Wenn $ x $ positiv ist, verschieben Sie $ 1 / \ text {pos} $ nach oben
- Wenn $ x $ negativ ist, verschieben Sie $ 1 / \ text {neg} $ nach rechts
wobei $ \ text {pos} $ und $ \ text {neg} $ die Bruchteile positiver bzw. negativer Beispiele sind.
Mit dieser Idee können Sie den AUC-ROC mithilfe des folgenden Algorithmus manuell berechnen:
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
Dieses schöne gif-animierte Bild sollte dies veranschaulichen Prozess klarer
Kommentare
- Danke @Alexey Grigorev, dies ist ein großartiges Bild und es wird sich wahrscheinlich in Zukunft als nützlich erweisen! +1
- Könnten Sie bitte etwas über “ Bruchteile positiver und negativer Beispiele “ erklären kleinster Einheitswert von zwei Achsen?
- @Allan Ruin:
pos
bedeutet hier die Anzahl der positiven Daten. Nehmen wir an, Sie haben 20 Datenpunkte, von denen 11 Punkte 1 sind. Wenn Sie also das Diagramm zeichnen, haben wir ein Rechteck 11×9 (Höhe x Breite). Alexey Grigorev hat zwar skaliert, aber lassen Sie es einfach so, wie es ‚ ist, wenn Sie möchten. Verschieben Sie jetzt bei jedem Schritt einfach 1 im Diagramm.
Antwort
Karls Beitrag hat viel zu bieten Aber ich habe in den letzten 20 Jahren noch kein Beispiel für eine ROC-Kurve gesehen, die das Denken eines Menschen in eine gute Richtung verändert hat. Der einzige Wert einer ROC-Kurve ist meiner bescheidenen Meinung nach, dass ihre Fläche einer sehr nützlichen Konkordanzwahrscheinlichkeit entspricht. Die ROC-Kurve selbst verleitet den Leser dazu, Grenzwerte zu verwenden, was eine schlechte statistische Praxis ist.
Wenn Sie den $ c $ -Index manuell berechnen, erstellen Sie einen Plot mit $ Y = 0,1 $ auf dem $ x $ -Achse und der kontinuierliche Prädiktor oder die vorhergesagte Wahrscheinlichkeit, dass $ Y = 1 $ auf der $ y $ -Achse ist. Wenn Sie jeden Punkt mit $ Y = 0 $ mit jedem Punkt mit $ Y = 1 $ verbinden, ist der Anteil der Linien mit einer positiven Steigung die Konkordanzwahrscheinlichkeit.
Alle Kennzahlen mit einem Nenner von $ n $ in dieser Einstellung sind falsche Regeln für die Genauigkeitsbewertung und sollten vermieden werden. Dies umfasst korrekt klassifizierte Proportionen, Sensitivität und Spezifität.
Drucken Sie für die Funktion R Hmisc
Paket rcorr.cens
die Funktion Das gesamte Ergebnis, um weitere Informationen zu erhalten, insbesondere einen Standardfehler.
Kommentare
- Vielen Dank, @Frank Harell, ich schätze Ihre Perspektive. Ich benutze einfach die c-Statistik als Konkordanzwahrscheinlichkeit, da ich ‚ keine Cutoffs mag. Nochmals vielen Dank!
Antwort
Hier ist eine Alternative zur natürlichen Methode der AUC-Berechnung durch einfaches Verwenden des Trapezes Regel, um die Fläche unter der ROC-Kurve zu erhalten.
Die AUC ist gleich der Wahrscheinlichkeit, dass eine zufällig ausgewählte positive Beobachtung eine vorhergesagte Wahrscheinlichkeit (positiv zu sein) hat, die größer ist als eine zufällig ausgewählte negative Beobachtung. Sie können dies verwenden, um die AUC in jeder Programmiersprache ganz einfach zu berechnen, indem Sie alle paarweisen Kombinationen von positiven und negativen Beobachtungen durchgehen. Sie können Beobachtungen auch zufällig auswählen, wenn die Stichprobengröße zu groß ist. Wenn Sie die AUC mit Stift und Papier berechnen möchten, ist dies möglicherweise nicht der beste Ansatz, es sei denn, Sie haben eine sehr kleine Stichprobe / viel Zeit. Zum Beispiel in 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
Wir können dies mit dem Paket pROC
überprüfen:
library(pROC) auc(y, probs) Area under the curve: 0.6287
Zufallsstichprobe verwenden:
mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896
Antwort
- Sie haben einen echten Wert für Beobachtungen.
- Berechnen Sie die hintere Wahrscheinlichkeit und ordnen Sie die Beobachtungen nach dieser Wahrscheinlichkeit.
- Unter der Annahme einer Grenzwahrscheinlichkeit von $ P $ und der Anzahl der Beobachtungen $ N $:
$$ \ frac {\ text {Summe der wahren Ränge} -0,5PN (PN + 1)} { PN (N-PN)} $$
Kommentare
- @ user73455 … 1) Ja, ich habe den wahren Wert für Beobachtungen. 2) Ist die hintere Wahrscheinlichkeit gleichbedeutend mit den vorhergesagten Wahrscheinlichkeiten für jede der Beobachtungen? 3) verstanden; Was ist jedoch “ Summe der wahren Ränge “ und wie berechnet man diesen Wert? Vielleicht hilft Ihnen ein Beispiel, diese Antwort genauer zu erklären? Vielen Dank!