I är intresserad av att beräkna arean under kurvan (AUC), eller c-statistiken, för hand för en binär logistisk regressionsmodell.

Till exempel i valideringsdataset, jag har det verkliga värdet för den beroende variabeln, retention (1 = retained; 0 = not retained), samt en förutsagd retentionstatus för varje observation som genereras av min regressionsanalys med hjälp av en modell som byggdes med hjälp av träningen set (detta kommer att sträcka sig från 0 till 1).

Mina första tankar var att identifiera det ”korrekta” antalet modellklassificeringar och helt enkelt dela antalet ”korrekta” observationer med antalet totala observationer att beräkna c-statistiken. Med ”korrekt”, om den verkliga retentionsstatusen för en observation = 1 och den förutsagda retentionsstatusen är> 0,5 så är det en ”korrekt” klassificering. Dessutom, om den verkliga retentionsstatusen för en observation = 0 och den förutsagda retentionsstatusen är < 0,5 så är det också en ”korrekt” klassificering. Jag antar att en ”tie” skulle inträffa när det förutspådda värdet = 0,5, men det fenomenet förekommer inte i min valideringsdataset. Å andra sidan skulle ”felaktiga” klassificeringar vara om den verkliga retentionsstatusen för en observation = 1 och den förutsagda retentionsstatusen är < 0,5 eller om den verkliga retentionsstatusen för ett resultat = 0 och den förväntade retentionsstatusen är> 0,5. Jag är medveten om TP, FP, FN, TN, men är inte medveten om hur man beräknar c-statistiken med den här informationen.

Svar

Jag skulle rekommendera Hanleys & McNeils 1982-papper Betydelsen och användningen av området under en mottagaregenskaper (ROC ) kurva .

Exempel

De har följande tabell över sjukdomsstatus och testresultat (motsvarar till exempel den uppskattade risken från en logistisk modell). Den första siffran till höger är antalet patienter med sann sjukdomsstatus normal och den andra siffran är antalet patienter med sann sjukdomsstatus onormal:

(1) Definitivt normalt: 33/3
(2) Förmodligen normalt: 6/2
(3) Tvivelaktigt: 6/2
(4) Troligtvis onormalt: 11/11
(5) Definitivt onormalt: 2/33

Så det finns totalt 58 normala patienter och 51 onormala. Vi ser att när prediktorn är 1, ”Definitivt normal”, är patienten vanligtvis normal (sant för 33 av de 36 patienterna), och när det är 5, ”Definitivt onormalt” är patienterna vanligtvis onormala (sant för 33 av 35 patienter), så prediktorn är vettig. Men hur ska vi bedöma en patient med poängen 2, 3 eller 4? Vad vi ställer in vår gräns för att bedöma en patient som onormal eller normal för att avgöra känsligheten och specificiteten för det resulterande testet.

Känslighet och specificitet

Vi kan beräkna uppskattat känslighet och specificitet för olika cutoffs. (Jag skriver bara ”känslighet” och ”specificitet” från och med nu, så att värderingarnas uppskattade natur blir underförstådda.)

Om vi väljer vår gräns så att vi klassificerar alla patienterna som onormala, oavsett vad deras testresultat säger (dvs, vi väljer avgränsningen 1+), får vi en känslighet på 51/51 = 1. Specificiteten blir 0/58 = 0. Inte låter så bra.

OK, så låt oss välja en mindre strikt avskärning. Vi klassificerar endast patienter som onormala om de har ett testresultat på 2 eller högre. Vi saknar sedan 3 onormala patienter och har en känslighet på 48/51 = 0,94. Men vi har en mycket ökad specificitet, 33/58 = 0,57.

Vi kan nu fortsätta detta genom att välja olika avbrott (3, 4, 5,> 5). (I det sista fallet klassificerar vi inte några patienter som onormala, även om de har högsta möjliga testpoäng på 5.)

ROC-kurvan

Om vi gör detta för alla möjliga cutoffs och plot känsligheten mot 1 minus specificiteten får vi ROC-kurvan. Vi kan använda följande R-kod:

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

Utgången är:

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

Vi kan beräkna olika statistik:

 ( 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  

Och med hjälp av detta kan vi plotta den (uppskattade) ROC-kurvan:

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

Manuell beräkning av AUC

Vi kan mycket enkelt beräkna arean under ROC-kurvan med hjälp av formeln för en trapetsformad area:

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

Resultatet är 0.8931711.

Ett överensstämmelsemått

AUC kan också ses som ett överensstämmelsemått.Om vi tar alla möjliga par av patienter där den ena är normal och den andra är onormal, kan vi beräkna hur ofta det är den onormala som har det högsta (mest ”onormala utseende”) testresultatet (om de har samma värde, vi räknar med att detta är en halv seger):

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

Svaret är återigen 0,8931711, området under ROC-kurvan. Detta kommer alltid att vara fallet.

En grafisk bild av överensstämmelse

Som Harrell påpekade i sitt svar har detta också en grafisk tolkning. Låt oss plotta testpoäng (riskuppskattning) på y -axen och sann sjukdomsstatus på x -axen (här med lite skakningar, för att visa överlappande punkter):

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

Spridningsdiagram över riskpoäng mot sann sjukdom status.

Låt oss nu rita en linje mellan varje punkt till vänster (en ”normal” patient) och varje punkt till höger (en ”onormal” patient). Andelen rader med en positiv lutning (dvs. andelen överensstämmande par) är överensstämningsindexet (platta linjer räknas som 50% överensstämmelse).

Det är lite svårt att visualisera de faktiska linjerna för detta exempel på grund av antalet band (lika riskpoäng), men med viss skakning och transparens kan vi få en rimlig plot:

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

Scatterdiagram över riskpoäng mot sann sjukdomsstatus, med linjer mellan alla möjliga observationspar.

Vi ser att de flesta av linjerna lutar uppåt, så konkordansindexet kommer att vara högt. Vi ser också bidraget till indexet från varje typ av observationspar. Det mesta kommer från normala patienter med en riskpoäng på 1 parat med onormala patienter med en riskpoäng på 5 (1–5 par), men en hel del kommer också från 1–4 par och 4–5 par. Och det är väldigt enkelt att beräkna det faktiska konkordansindexet baserat på lutningsdefinitionen:

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

Svaret är återigen 0,8931711, dvs AUC.

Wilcoxon – Mann – Whitney-testet

Det finns en nära koppling mellan överensstämmelsemåttet och Wilcoxon – Mann – Whitney testa. Egentligen testar den senare om sannolikheten för överensstämmelse (dvs. att det är den onormala patienten i ett slumpmässigt normalt – onormalt par som får det mest ”onormala utseende” testresultatet) är exakt 0,5. Och dess teststatistik är bara en enkel omvandling av den uppskattade överensstämmelsesannolikheten:

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

Teststatistiken (W = 2642) räknar antalet överensstämmande par. Om vi delar det med antalet möjliga par får vi ett känt nummer:

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

Ja, det är 0,8931711, området under ROC-kurvan.

Enklare sätt att beräkna AUC (i R)

Men låt oss göra livet lättare för oss själva. Det finns olika paket som automatiskt beräknar AUC för oss.

Epi-paketet

Epi -paketet skapar en fin ROC-kurva med olika statistik (inklusive AUC) inbäddad:

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

ROC-kurva från Epi-paketet

pROC-paketet

Jag gillar också paketet pROC eftersom det kan jämna ut ROC-uppskattningen (och beräkna en AUC-uppskattning baserat på den utjämnade ROC):

ROC-kurva (ojämn och utjämnad) från pROC-paketet

(Den röda linjen är den ursprungliga ROC, och den svarta linjen är den utjämnade ROC. Observera också standardformatet 1: 1. Det är vettigt att använda detta, eftersom både känslighet och specificitet har 0–1 intervall.)

Den beräknade AUC från utjämnade ROC är 0,9107, liknar, men något större än, AUC från den ojämna ROC (om du ser en t figuren kan du enkelt se varför den är större). (Även om vi verkligen har för få möjliga distinkta testresultatvärden för att beräkna en jämn AUC).

rms-paketet

Harrells rms -paket kan beräkna olika relaterade överensstämmelsestatistik med funktionen rcorr.cens(). C Index i dess utgång är AUC:

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

CaTools-paketet

Slutligen har vi caTools -paketet och dess colAUC() -funktion. Det har några fördelar jämfört med andra paket (främst hastighet och förmågan att arbeta med flerdimensionell data – se ?colAUC) som kan ibland vara till hjälp.Men givetvis ger det samma svar som vi har beräknat om och om igen:

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

ROC-kurva från caTools-paketet

Slutord

Många verkar tro att AUC säger till oss hur ”bra” ett test är. Och vissa tror att AUC är sannolikheten för att testet kommer att klassificera en patient korrekt. Det är inte . Som du kan se från ovanstående exempel och beräkningar berättar AUC för oss något om en familj av test, ett test för varje möjlig cutoff.

Och AUC beräknas baserat på cutoffs man aldrig skulle använda i praktiken. Varför bör vi bry oss om känsligheten och specificiteten hos meningslösa gränsvärden? Ändå är det vad AUC bygger (delvis) på. (Naturligtvis, om AUC är väldigt nära 1, kommer nästan alla möjliga test ha stor diskriminerande kraft, och vi skulle alla vara väldigt glada.)

Den slumpmässiga normalen –Abnormalt ”partolkning av AUC är trevligt (och kan utvidgas till exempel till överlevnadsmodeller, där vi ser om det är personen med den högsta (relativa risken som dör tidigast). Men man skulle aldrig använda det i praktiken. Det är ett sällsynt fall där man vet man har en frisk och en sjuk person, inte vet vilken person som är sjuk och måste bestäm vilken av dem som ska behandlas. (Hur som helst är beslutet enkelt; behandla den med den högsta uppskattade risken.)

Så jag tror att studera den faktiska ROC-kurvan kommer att vara mer användbar än att bara titta på AUC: s sammanfattande åtgärd. Och om du använder ROC tillsammans med (uppskattningar av) kostnader för falska positiva och falska negativ, tillsammans med basfrekvenser för vad du studerar, kan du komma någonstans.

Observera också att AUC bara mäter diskriminering , inte kalibrering. Det vill säga, det mäter om du kan skilja mellan två personer (en sjuk och en frisk) baserat på riskpoängen. För detta tittar det bara på relativa riskvärden (eller rankas, om du vill, jfr Wilcoxon – Mann – Whitney-tolkningen), inte de absoluta, som du ska vara intresserad av. Om du till exempel delar upp varje riskuppskattning från din logistikmodell med 2 får du exakt samma AUC (och ROC).

När du utvärderar en riskmodell, kalibrering är också mycket viktigt. För att undersöka detta kommer du att titta på alla patienter med en riskpoäng på runt, t.ex. 0,7, och se om cirka 70% av dessa faktiskt var sjuka. Gör detta för varje möjlig riskpoäng (möjligen med någon form av utjämning / lokal regression). Rita upp resultaten så får du ett grafiskt mått på kalibrering .

Om du har en modell med både bra kalibrering och god diskriminering, då börja ha en bra modell. 🙂

Kommentarer

  • Tack, @Karl Ove Hufthammer, detta är det mest grundliga svaret jag någonsin har fått. Jag uppskattar särskilt ditt ” Slutord ” avsnitt. Utmärkt arbete! Tack igen!
  • Tack så mycket för det detaljerade svaret. Jag arbetar med en dataset där Epi :: ROC () v2.2.6 är övertygad om att AUC är 1.62 (nej det är inte en mentaliststudie), men enligt ROC tror jag mycket mer på 0.56 som ovanstående kod resulterar in.
  • Jag tror att det finns ett litet fel i sens=c(sens,0); omspec=c(omspec,0), borde ’ t detta vara sens=c(0, sens); omspec=c(0, omspec)? Det plottas korrekt med det ledande 0 men inte som det för närvarande är i svaret.
  • Nej, den nuvarande definitionen är, AFAICS, korrekt, @steveb, och resulterar i en korrekt plot. Jag tror att det som kanske är förvirrande är att ROC-kurvan dras från höger till vänster (dvs. från det övre högra hörnet till det nedre vänstra hörnet), inte från vänster till det rätt , som de flesta tomter är. Det är bara resultatet av hur jag definierade variablerna; man kunde lika gärna ha plottat det från vänster till höger (genom att vända både sens och omspec innan man plottade).

Svar

Ta en titt på den här frågan: Förstå ROC-kurvan

Så här bygger du en ROC-kurva (från den frågan):

Ritar ROC-kurva

givet en datamängd som bearbetats av din rankingklassificering

  • exempel på rankningstest på minskande poäng
  • startar i $ (0, 0) $
  • för varje exempel $ x $ (i minskande ordning)
    • om $ x $ är positivt, flytta $ 1 / \ text {pos} $ upp
    • om $ x $ är negativt, flytta $ 1 / \ text {neg} $ höger

där $ \ text {pos} $ och $ \ text {neg} $ är fraktionerna av positiva respektive negativa exempel.

Du kan använda denna idé för att manuellt beräkna AUC ROC med följande algoritm:

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 

Denna fina gif-animerade bild ska illustrera detta processrensare

bygga kurvan

Kommentarer

  • Tack @Alexey Grigorev, det här är en fantastisk bild och det kommer sannolikt att vara användbart i framtiden! +1
  • Kan du förklara lite om ” fraktioner av positiva och negativa exempel ”, menar du minsta enhetsvärdet på två axlar?
  • @Allan Ruin: pos betyder här antalet positiva data. Låt oss säga att du har 20 datapunkter, där 11 poäng är 1. Så när vi ritar diagrammet har vi en rektangel 11×9 (höjd x bredd). Alexey Grigorev gjorde skalning men lät det bara vara som det ’ är om du vill. Nu är det bara att flytta 1 i diagrammet vid varje steg.

Svar

Karls inlägg har mycket men jag har ännu inte sett de senaste 20 åren ett exempel på en ROC-kurva som förändrade någons tänkande i en bra riktning. Det enda värdet av en ROC-kurva enligt min ödmjuka uppfattning är att dess område råkar motsvara en mycket användbar överensstämmelsessannolikhet. Själva ROC-kurvan frestar läsaren att använda cutoffs, vilket är dålig statistisk praxis.

När det gäller manuell beräkning av $ c $ -index, gör en plot med $ Y = 0,1 $ på $ x $ -axis och den kontinuerliga prediktorn eller förutsagd sannolikhet att $ Y = 1 $ på $ y $ -axeln. Om du förbinder varje punkt med $ Y = 0 $ med varje punkt med $ Y = 1 $, är andelen linjer som har en positiv lutning sannolikheten för överensstämmelse.

Alla mått som har en nämnare på $ n $ i denna inställning är felaktiga regler för noggrannhet och bör undvikas. Detta inkluderar proportion klassificerad korrekt, känslighet och specificitet.

För R Hmisc -paketet rcorr.cens -funktionen, skriv ut hela resultatet för att se mer information, särskilt ett standardfel.

Kommentarer

  • Tack, @Frank Harell, jag uppskattar ditt perspektiv. Jag använder helt enkelt c-statistiken som en överensstämmelsesannolikhet, eftersom jag inte ’ inte gillar cutoffs. Tack igen!

Svar

Här är ett alternativ till det naturliga sättet att beräkna AUC genom att helt enkelt använda den trapetsformade regel för att få området under ROC-kurvan.

AUC är lika med sannolikheten att en slumpmässigt samplad positiv observation har en förutsagd sannolikhet (att vara positiv) större än en slumpmässigt samplad negativ observation. Du kan använda detta för att beräkna AUC ganska enkelt i vilket programmeringsspråk som helst genom att gå igenom alla parvisa kombinationer av positiva och negativa observationer. Du kan också slumpmässigt prova observationer om urvalsstorleken var för stor. Om du vill beräkna AUC med penna och papper kanske det här inte är den bästa metoden om du inte har ett mycket litet urval / mycket tid. Till exempel i 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 

Vi kan verifiera med pROC -paketet:

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

Med slumpmässig sampling:

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

Svar

  1. Du har sant värde för observationer.
  2. Beräkna posterior sannolikhet och rangordna sedan observationer efter denna sannolikhet.
  3. Förutsatt att gränsen är sannolik för $ P $ och antal observationer $ N $:
    $$ \ frac {\ text {Summan av sanna rankningar} -0.5PN (PN + 1)} { PN (N-PN)} $$

Kommentarer

  • @ user73455 … 1) Ja, jag har det verkliga värdet för observationer. 2) Är posterior sannolikhet synonym med förutsagda sannolikheter för var och en av observationerna? 3) Förstått; vad är dock ” Summan av sanna rankningar ” och hur beräknar man detta värde? Kanske skulle ett exempel hjälpa dig att förklara detta svar mer ingående? Tack!

Lämna ett svar

Din e-postadress kommer inte publiceras. Obligatoriska fält är märkta *