I ben geïnteresseerd in het handmatig berekenen van de oppervlakte onder de curve (AUC), of de c-statistiek, voor een binair logistisch regressiemodel.

Bijvoorbeeld in de validatiedataset, heb ik de werkelijke waarde voor de afhankelijke variabele, retentie (1 = behouden; 0 = niet behouden), evenals een voorspelde retentiestatus voor elke observatie gegenereerd door mijn regressieanalyse met behulp van een model dat is gebouwd met behulp van de training set (dit zal variëren van 0 tot 1).

Mijn eerste gedachten waren om het juiste aantal modelclassificaties te identificeren en het aantal juiste waarnemingen te delen door het totale aantal te berekenen waarnemingen de c-statistiek. Met “correct”, als de ware retentiestatus van een waarneming = 1 en de voorspelde retentiestatus> 0,5 is, dan is dat een “correcte” classificatie. Bovendien, als de ware retentiestatus van een observatie = 0 en de voorspelde retentiestatus < 0,5 is, dan is dat ook een “correcte” classificatie. Ik neem aan dat er een “gelijkspel” zou optreden als de voorspelde waarde = 0,5, maar dat fenomeen komt niet voor in mijn validatiedataset. Aan de andere kant zouden onjuiste classificaties zijn als de ware bewaarstatus van een waarneming = 1 en de voorspelde bewaarstatus < 0,5 is of als de ware bewaarstatus voor een uitkomst = 0 en de voorspelde retentiestatus is> 0,5. Ik ken TP, FP, FN, TN, maar ik weet niet hoe ik de c-statistiek moet berekenen op basis van deze informatie.

Answer

Ik zou Hanleys & McNeils 1982 paper De betekenis en het gebruik van het gebied onder een ontvangerkarakteristiek (ROC ) curve .

Voorbeeld

Ze hebben de volgende tabel met ziektestatus en testresultaten (overeenkomend met bijvoorbeeld het geschatte risico van een logistiek model). Het eerste getal aan de rechterkant is het aantal patiënten met waar ziektestatus normaal en het tweede nummer is het aantal patiënten met waar ziektestatus abnormaal:

(1) Absoluut normaal: 33/3
(2) Waarschijnlijk normaal: 6/2
(3) Twijfelachtig: 6/2
(4) Waarschijnlijk abnormaal: 11/11
(5) Absoluut abnormaal: 2/33

Er zijn dus in totaal 58 normale patiënten en 51 abnormale. We zien dat wanneer de voorspeller 1 is, Absoluut normaal, de patiënt meestal normaal is (waar voor 33 van de 36 patiënten), en als het 5 is Absoluut abnormaal, de patiënt meestal abnormaal is (waar voor 33 van de 36 patiënten). 35 patiënten), dus de voorspeller is logisch. Maar hoe moeten we een patiënt beoordelen met een score van 2, 3 of 4? Wat we als grenswaarde hebben ingesteld om een patiënt als abnormaal of normaal te beoordelen, bepaalt de gevoeligheid en specificiteit van de resulterende test.

Gevoeligheid en specificiteit

We kunnen de geschatte gevoeligheid en specificiteit voor verschillende cutoffs. (Ik zal vanaf nu gewoon gevoeligheid en specificiteit schrijven, waarbij ik de geschatte aard van de waarden impliciet laat zijn.)

Als we onze cutoff kiezen zodat we alle de patiënten als abnormaal, ongeacht wat hun testresultaten zeggen (dwz we kiezen de cutoff 1+), we krijgen een gevoeligheid van 51/51 = 1. De specificiteit zal 0/58 = 0 zijn. klinken zo goed.

OK, dus laten we een minder strikte cutoff kiezen. We classificeren patiënten alleen als abnormaal als ze een testresultaat van 2 of hoger hebben. We missen dan 3 abnormale patiënten en hebben een gevoeligheid van 48/51 = 0,94. Maar we hebben een veel verhoogde specificiteit, van 33/58 = 0,57.

We kunnen nu doorgaan met het kiezen van verschillende cutoffs (3, 4, 5,> 5). (In het laatste geval zullen we geen patiënten als abnormaal classificeren, zelfs niet als ze de hoogst mogelijke testscore van 5 hebben.)

De ROC-curve

Als we dit doen voor alle mogelijke cutoffs, en de gevoeligheid uitzetten tegen 1 minus de specificiteit, krijgen we de ROC-curve. We kunnen de volgende R-code gebruiken:

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

De output is:

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

We kunnen verschillende statistieken berekenen:

 ( 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  

En hiermee kunnen we de (geschatte) ROC-curve plotten:

 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)  

AUC-curve

Handmatig de AUC

We kunnen heel gemakkelijk de oppervlakte onder de ROC-curve berekenen met behulp van de formule voor de oppervlakte van een trapezium:

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

Het resultaat is 0.8931711.

Een concordantiemaatstaf

De AUC kan ook worden gezien als een concordantiemaatstaf.Als we alle mogelijke paren patiënten nemen waarbij de ene normaal is en de ander abnormaal, kunnen we berekenen hoe vaak de abnormale patiënt het hoogste (meest abnormaal ogende) testresultaat heeft (als ze hebben dezelfde waarde, we rekenen dat als een halve overwinning):

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

Het antwoord is wederom 0.8931711, het gebied onder de ROC-curve. Dit zal altijd het geval zijn.

Een grafische weergave van concordantie

Zoals aangegeven door Harrell in zijn antwoord, heeft dit ook een grafische interpretatie. Laten we de testscore (risicoschatting) uitzetten op de y -as en de werkelijke ziektestatus op de x -as (hier met wat trillingen, om overlappende punten te tonen):

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

Scatterplot van risicoscore tegen echte ziekte status.

Laten we nu een lijn trekken tussen elk punt aan de linkerkant (een normale patiënt) en elk punt aan de rechterkant (een abnormale patiënt). De proportie lijnen met een positieve helling (d.w.z. de proportie concordante paren) is de concordantie-index (vlakke lijnen tellen als ‘50% concordantie ’).

Het is een beetje moeilijk om de feitelijke lijnen voor dit voorbeeld te visualiseren, vanwege het aantal gelijkspel (gelijke risicoscore), maar met wat trillingen en transparantie kunnen we een redelijk plot krijgen:

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

Scatterplot van risicoscore ten opzichte van de werkelijke ziektestatus, met lijnen tussen alle mogelijke observatieparen.

We zien dat de meeste lijnen omhoog lopen, dus de concordantie-index zal hoog zijn. We zien ook de bijdrage aan de index van elk type observatiepaar. Het meeste is afkomstig van normale patiënten met een risicoscore van 1 in combinatie met abnormale patiënten met een risicoscore van 5 (1–5 paren), maar er is ook veel afkomstig van 1–4 paren en 4–5 paren. En het is heel gemakkelijk om de werkelijke concordantie-index te berekenen op basis van de hellingsdefinitie:

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

Het antwoord is opnieuw 0.8931711, dwz de AUC.

De Wilcoxon – Mann – Whitney-test

Er is een nauw verband tussen de concordantiemaat en de Wilcoxon – Mann – Whitney test. De laatste test eigenlijk of de kans op overeenstemming (d.w.z. dat het de abnormale patiënt in een willekeurig normaal-abnormaal paar is dat het meest ‘abnormaal uitziende’ testresultaat zal hebben) precies 0,5 is. En de teststatistiek is slechts een eenvoudige transformatie van de geschatte concordantiekans:

 > ( 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  

De teststatistiek (W = 2642) telt het aantal concordante paren. Als we het delen door het aantal mogelijke paren, krijgen we een bekend getal:

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

Ja, het is 0,8931711, het gebied onder de ROC-curve.

Gemakkelijkere manieren om de AUC (in R) te berekenen

Maar laten we het onszelf gemakkelijker maken. Er zijn verschillende pakketten die de AUC automatisch voor ons berekenen.

Het Epi-pakket

Het Epi -pakket creëert een mooie ROC-curve met verschillende statistieken (inclusief de AUC) ingesloten:

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

ROC-curve van het Epi-pakket

Het pROC-pakket

Ik hou ook van het pROC -pakket, omdat het kan de ROC-schatting afvlakken (en een AUC-schatting berekenen op basis van de afgevlakte ROC):

ROC-curve (ongestroomlijnd en afgevlakt) uit het pROC-pakket

(De rode lijn is de originele ROC en de zwarte lijn is de afgevlakte ROC. Let ook op de standaard beeldverhouding van 1: 1. Het is logisch om deze te gebruiken, aangezien zowel de gevoeligheid als de specificiteit een 0-1 bereik.)

De geschatte AUC van de afgevlakte ROC is 0,9107, vergelijkbaar met, maar iets groter dan, de AUC van de niet-afgevlakte ROC (als je een In de figuur kun je gemakkelijk zien waarom het groter is). (Hoewel we echt te weinig mogelijke verschillende testresultaatwaarden hebben om een soepele AUC te berekenen).

Het rms-pakket

Harrells rms pakket kan verschillende gerelateerde concordantiestatistieken berekenen met de functie rcorr.cens(). De C Index in zijn uitvoer is de AUC:

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

Het pakket caTools

Ten slotte hebben we het caTools pakket en zijn colAUC() functie. Het heeft een paar voordelen ten opzichte van andere pakketten (voornamelijk snelheid en de mogelijkheid om met multidimensionale gegevens te werken – zie ?colAUC) die soms kunnen helpen.Maar het geeft natuurlijk hetzelfde antwoord als we keer op keer hebben berekend:

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

ROC-curve uit het caTools-pakket

Laatste woorden

Veel mensen lijken te denken dat de AUC ons vertelt hoe goed een test is. En sommige mensen denken dat de AUC de kans is dat de test een patiënt correct classificeert. Het is niet . Zoals je kunt zien aan de hand van het bovenstaande voorbeeld en de berekeningen, vertelt de AUC ons iets over een familie van tests, één test voor elke mogelijke cutoff.

En de AUC wordt berekend op basis van cutoffs die je in de praktijk nooit zou gebruiken. Waarom zouden we ons druk maken over de gevoeligheid en specificiteit van ‘onzinnige’ grenswaarden? Toch is dat waar de AUC (gedeeltelijk) op is gebaseerd. (Natuurlijk, als de AUC zeer dicht bij 1 ligt, zal bijna elke mogelijke test een groot onderscheidingsvermogen hebben, en we zouden allemaal erg blij zijn.)

De random normal –Abnormale paarinterpretatie van de AUC is mooi (en kan worden uitgebreid naar bijvoorbeeld overlevingsmodellen, waar we zien of de persoon met het hoogste (relatieve) risico het vroegst sterft). Maar je zou het in de praktijk nooit gebruiken. Het komt zelden voor dat men weet dat men één gezonde en één zieke persoon heeft, niet weet welke persoon de zieke is, en beslissen welke van hen te behandelen. (In elk geval is de beslissing eenvoudig; behandel degene met het hoogste geschatte risico.)

Dus ik denk dat het bestuderen van de werkelijke ROC-curve nuttiger zal zijn dan alleen kijken naar de AUC-samenvattende maatregel. En als je het ROC gebruikt samen met (schattingen van de) kosten van fout-positieven en fout-negatieven, samen met basistarieven van wat je studeert, kun je ergens komen.

Merk ook op dat de AUC alleen discriminatie meet, niet kalibratie. Dat wil zeggen, het meet of je onderscheid kunt maken tussen twee personen (een zieke en een gezonde) op basis van de risicoscore. Hiervoor kijkt het alleen naar relatieve risicowaarden (of rangschikt, als je wilt, cf. de Wilcoxon-Mann-Whitney-testinterpretatie), niet de absolute, die je zou moeten geïnteresseerd zijn. Als u bijvoorbeeld elke risicoschatting van uw logistieke model door 2 deelt, krijgt u exact dezelfde AUC (en ROC).

Bij het evalueren van een risicomodel, kalibratie is ook erg belangrijk. Om dit te onderzoeken, kijk je naar alle patiënten met een risicoscore van rond, bijvoorbeeld 0,7, en kijk je of ongeveer 70% van hen daadwerkelijk ziek was. Doe dit voor elke mogelijke risicoscore (mogelijk met behulp van een soort afvlakking / lokale regressie). Zet de resultaten uit en je krijgt een grafische maat voor kalibratie .

Als je een model hebt met beide goede kalibratie en goede discriminatie, dan begin een goed model te hebben. 🙂

Reacties

  • Bedankt, @Karl Ove Hufthammer, dit is het meest grondige antwoord dat ik ooit heb gekregen. Ik waardeer vooral uw ” Laatste woorden ” sectie. Goed werk! Nogmaals bedankt!
  • Hartelijk dank voor dit gedetailleerde antwoord. Ik werk met een dataset waarin Epi :: ROC () v2.2.6 ervan overtuigd is dat de AUC 1,62 is (nee het is geen mentalistische studie), maar volgens het ROC geloof ik veel meer in de 0,56 die de bovenstaande code oplevert in.
  • Ik denk dat er een kleine fout is in sens=c(sens,0); omspec=c(omspec,0), zou niet ‘ moeten zijn sens=c(0, sens); omspec=c(0, omspec)? Het plot correct met de leidende 0 maar niet zoals het momenteel in het antwoord staat.
  • Nee, de huidige definitie is, AFAICS, correct, @steveb, en resulteert in een correct plot. Ik denk dat het misschien verwarrend is dat de ROC-curve van rechts naar links wordt getrokken (dwz van de rechterbovenhoek naar de linkerbenedenhoek), niet van links naar de rechts , zoals de meeste plots zijn. Dat is slechts het resultaat van hoe ik de variabelen heb gedefinieerd; men had het net zo goed van links naar rechts kunnen plotten (door zowel sens als omspec om te draaien voordat u gaat plotten).

Antwoord

Bekijk deze vraag eens: ROC-curve begrijpen

Hier is hoe je een ROC-curve bouwt (op basis van die vraag):

ROC-curve tekenen

gegeven een dataset verwerkt door uw ranking classifier

  • testvoorbeelden rangschikken op afnemende score
  • beginnen in $ (0, 0) $
  • voor elk voorbeeld $ x $ (in de aflopende volgorde)
    • als $ x $ positief is, verplaats $ 1 / \ text {pos} $ omhoog
    • als $ x $ negatief is, verplaats $ 1 / \ text {neg} $ rechts

waarbij $ \ text {pos} $ en $ \ text {neg} $ respectievelijk de fracties zijn van positieve en negatieve voorbeelden.

U kunt dit idee gebruiken voor het handmatig berekenen van de AUC ROC met behulp van het volgende algoritme:

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 

Deze leuke gif-geanimeerde afbeelding zou dit moeten illustreren proces duidelijker

de curve opbouwen

Reacties

  • Bedankt @Alexey Grigorev, dit is een geweldige visual en het zal waarschijnlijk in de toekomst nuttig blijken! +1
  • Kan iets uitleggen over ” fracties van positieve en negatieve voorbeelden “, bedoelt u de kleinste eenheidswaarde van twee assen?
  • @Allan Ruin: pos betekent hier het aantal positieve gegevens. Stel dat je 20 gegevenspunten hebt, waarbij 11 punten 1 zijn. Dus bij het tekenen van de grafiek hebben we een rechthoek 11×9 (hoogte x breedte). Alexey Grigorev heeft wel geschaald, maar laat het gewoon zoals het ‘ s is als je wilt. Verplaats nu bij elke stap 1 op de kaart.

Antwoord

Karls bericht bevat veel van uitstekende informatie. Maar ik heb de afgelopen 20 jaar nog geen voorbeeld gezien van een ROC-curve die iemands denken in een goede richting veranderde. De enige waarde van een ROC-curve is naar mijn bescheiden mening dat het gebied toevallig gelijk is aan een zeer nuttige concordantiekans. De ROC-curve zelf verleidt de lezer om cutoffs te gebruiken, wat een slechte statistische praktijk is.

Wat betreft het handmatig berekenen van de $ c $ -index, maak een plot met $ Y = 0,1 $ op de $ x $ -as en de continue voorspeller of voorspelde kans dat $ Y = 1 $ op de $ y $ -as. Als je elk punt met $ Y = 0 $ verbindt met elk punt met $ Y = 1 $, is het deel van de lijnen dat een positieve helling heeft de concordantiekans.

Alle maten met een noemer van $ n $ in deze instelling zijn onjuiste regels voor nauwkeurigheidsscore en moeten worden vermeden. Dit omvat de correct geclassificeerde proporties, gevoeligheid en specificiteit.

Voor het R Hmisc pakket rcorr.cens functie, drukt u de volledige resultaat om meer informatie te zien, vooral een standaardfout.

Opmerkingen

  • Dank je, @Frank Harell, ik waardeer je perspectief. Ik gebruik gewoon de c-statistiek als een concordantiekans, aangezien ik ‘ niet van afkappunten houd. Nogmaals bedankt!

Antwoord

Hier is een alternatief voor de natuurlijke manier om de AUC te berekenen door simpelweg de trapeziumvormige regel om het gebied onder de ROC-curve te krijgen.

De AUC is gelijk aan de kans dat een willekeurig bemonsterde positieve waarneming een grotere voorspelde waarschijnlijkheid (om positief te zijn) heeft dan een willekeurig bemonsterde negatieve waarneming. U kunt dit gebruiken om de AUC vrij gemakkelijk in elke programmeertaal te berekenen door alle paarsgewijze combinaties van positieve en negatieve observaties te doorlopen. U kunt ook steekproeven nemen van waarnemingen als de steekproefomvang te groot was. Als u de AUC wilt berekenen met pen en papier, is dit misschien niet de beste aanpak, tenzij u een zeer kleine steekproef / veel tijd heeft. Bijvoorbeeld 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 

We kunnen verifiëren met behulp van het pROC pakket:

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

Willekeurige steekproeven gebruiken:

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

Antwoord

  1. Je hebt echte waarde voor observaties.
  2. Bereken de posterieure waarschijnlijkheid en rangschik de waarnemingen vervolgens op basis van deze waarschijnlijkheid.
  3. Uitgaande van de afkapkans van $ P $ en het aantal waarnemingen $ N $:
    $$ \ frac {\ text {Som of true ranks} -0,5PN (PN + 1)} { PN (N-PN)} $$

Reacties

  • @ user73455 … 1) Ja, ik heb de echte waarde voor observaties. 2) Is posterieure waarschijnlijkheid synoniem met voorspelde waarschijnlijkheden voor elk van de waarnemingen? 3) Begrepen; wat is echter ” Som van echte rangen ” en hoe berekent men deze waarde? Misschien helpt een voorbeeld u dit antwoord grondiger uit te leggen? Bedankt!

Geef een reactie

Het e-mailadres wordt niet gepubliceerd. Vereiste velden zijn gemarkeerd met *