I Sunt interesat să calculez aria sub curba (AUC) sau statistica c, manual, pentru un model de regresie logistică binară.
De exemplu, în setul de date de validare, am adevărata valoare pentru variabila dependentă, retenție (1 = reținută; 0 = ne reținută), precum și o stare de reținere previzionată pentru fiecare observație generată de analiza mea de regresie utilizând un model care a fost construit folosind instruirea set (aceasta va varia de la 0 la 1).
Gândurile mele inițiale au fost să identific numărul „corect” al clasificărilor modelului și pur și simplu să împărțim numărul de observații „corecte” la numărul de observații totale de calculat statistica c. Prin „corect”, dacă starea de reținere adevărată a unei observații = 1 și starea de reținere prezisă este> 0,5 atunci aceasta este o clasificare „corectă”. În plus, dacă starea de reținere adevărată a unei observații = 0 și starea de reținere prezisă este < 0,5, atunci aceasta este, de asemenea, o clasificare „corectă”. Presupun că o „egalitate” ar apărea atunci când valoarea prezisă = 0,5, dar acel fenomen nu apare în setul meu de date de validare. Pe de altă parte, clasificările „incorecte” ar fi dacă starea de păstrare adevărată a unei observații = 1 și starea de păstrare anticipată este < 0,5 sau dacă starea de păstrare adevărată pentru un rezultat = 0 și starea de păstrare prevăzută este> 0,5. Sunt conștient de TP, FP, FN, TN, dar nu știu cum să calculez statistica c, având în vedere aceste informații.
Răspuns
Aș recomanda lui Hanley & lucrarea lui McNeil din 1982 Semnificația și utilizarea zonei sub o caracteristică de funcționare a receptorului (ROC) ) curba .
Exemplu
Au următorul tabel cu starea bolii și rezultatul testului (corespunzător, de exemplu, riscului estimat dintr-un model logistic). Primul număr din dreapta este numărul de pacienți cu starea de boală adevărată „normală” și al doilea număr este numărul de pacienți cu starea de boală adevărată „anormală”:
(1) Cu siguranță normal: 33/3
(2) Probabil normal: 6/2
(3) Întrebabil: 6/2
(4) Probabil anormal: 11/11
(5) Cu siguranță anormală: 2/33
Deci există în total 58 de pacienți „normali” și „51” anormali. Vedem că atunci când predictorul este 1, „Cu siguranță normal”, pacientul este de obicei normal (adevărat pentru 33 din cei 36 de pacienți), iar când este 5, „Cu siguranță anormal”, pacienții sunt de obicei anormali (adevărat pentru 33 dintre cei 35 de pacienți), astfel încât predictorul are sens. Dar cum ar trebui să judecăm un pacient cu un scor de 2, 3 sau 4? Determinarea sensibilității și specificității testului rezultat este ceea ce am stabilit pentru evaluarea unui pacient ca fiind anormal sau normal.
Sensibilitate și specificitate
Putem calcula estimatul sensibilitate și specificitate pentru diferite limite. (Voi scrie „sensibilitate” și „specificitate” de acum înainte, lăsând natura estimată a valorilor să fie implicită.)
Dacă ne alegem limita, astfel încât să clasificăm toate pacienții ca anormali, indiferent de ceea ce spun rezultatele testelor (de exemplu, alegem limita 1+), vom obține o sensibilitate de 51/51 = 1. Specificitatea va fi 0/58 = 0. Nu sună atât de bine.
OK, deci să alegem o limită mai puțin strictă. Clasificăm pacienții ca anormali numai dacă au un rezultat al testului de 2 sau mai mare. Ne lipsesc apoi 3 pacienți anormali și avem o sensibilitate de 48/51 = 0,94. Dar avem o specificitate mult crescută, de 33/58 = 0,57.
Acum putem continua acest lucru, alegând diferite limite (3, 4, 5,> 5). (În ultimul caz, nu vom clasifica niciunul dintre pacienți ca anormali, chiar dacă au cel mai mare scor de test posibil de 5.)
Curba ROC
Dacă facem acest lucru pentru toate limitele posibile, și trasăm sensibilitatea în raport cu 1 minus specificitatea, obținem curba ROC. Putem utiliza următorul cod R:
# 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)) )
Rezultatul este:
testres truestat 1 2 3 4 5 0 33 6 6 11 2 1 3 2 2 11 33
Putem calcula diverse statistici:
( 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
Și folosind aceasta, putem trasa curba (estimată) ROC:
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ând manual AUC
Putem calcula foarte ușor aria de sub curba ROC, folosind formula pentru aria unui trapez:
height = (sens[-1]+sens[-length(sens)])/2 width = -diff(omspec) # = diff(rev(omspec)) sum(height*width)
Rezultatul este 0.8931711.
O măsură de concordanță
ASC poate fi văzută și ca o măsură de concordanță.Dacă luăm toate perechile posibile de pacienți în care unul este normal și celălalt este anormal, putem calcula cât de des este cel anormal care are cel mai mare rezultat al testului (cel mai „cu aspect anormal”) (dacă au aceeași valoare, considerăm că aceasta este „jumătate de victorie”):
o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))
Răspunsul este din nou 0.8931711, zona de sub curba ROC. Acesta va fi întotdeauna cazul.
O vedere grafică a concordanței
După cum a subliniat Harrell în răspunsul său, aceasta are și o interpretare grafică. Să trasăm scorul testului (estimarea riscului) pe axa y și starea adevărată a bolii pe axa x (aici cu o anumită perturbare, pentru a arăta punctele suprapuse):
plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")
Să trasăm acum o linie între fiecare punct din stânga (un pacient „normal”) și fiecare punct din dreapta (un pacient „anormal”). Proporția liniilor cu o pantă pozitivă (de exemplu, proporția de perechi concordante ) este indicele de concordanță (liniile plate sunt „50% concordanță”).
Este un pic dificil să vizualizezi liniile reale pentru acest exemplu, datorită numărului de egalități (scor de risc egal), dar cu o anumită perturbare și transparență putem obține un complot rezonabil:
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))
Vedem că majoritatea liniilor înclină în sus, deci indicele de concordanță va fi ridicat. Vedem, de asemenea, contribuția la index din fiecare tip de pereche de observație. Cea mai mare parte provine de la pacienți normali cu un scor de risc de 1 asociat cu pacienți anormali cu un scor de risc de 5 (1-5 perechi), dar destul de multe provin și de la 1-4 perechi și 4-5 perechi. Și este foarte ușor să calculați indicele de concordanță real pe baza definiției pantei:
d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm)) mean((d$slope > 0) + .5*(d$slope==0))
Răspunsul este din nou 0.8931711, adică ASC.
Testul Wilcoxon – Mann – Whitney
Există o legătură strânsă între măsura de concordanță și Wilcoxon – Mann – Whitney Test. De fapt, acesta din urmă testează dacă probabilitatea de concordanță (adică faptul că este pacientul anormal dintr-o pereche aleatorie normală-anormală care va avea rezultatul testului cel mai „cu aspect anormal”) este exact 0,5. Și statistica sa de testare este doar o simplă transformare a probabilității de concordanță estimată:
> ( 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
Statistica testului (W = 2642
) contează numărul de perechi concordante. Dacă îl împărțim la numărul de perechi posibile, obținem un număr familar:
w = wi$statistic w/(length(abnorm)*length(norm))
Da, este 0.8931711, zona de sub curba ROC.
Moduri mai ușoare de a calcula ASC (în R)
Dar hai să ne ușurăm viața. Există diverse pachete care calculează automat ASC pentru noi.
Pachetul Epi
Pachetul Epi
creează o curbă ROC frumoasă cu diverse statistici (inclusiv AUC) încorporate:
library(Epi) ROC(testres, truestat) # also try adding plot="sp"
Pachetul pROC
Îmi place și pachetul pROC
, deoarece poate neteziți estimarea ROC (și calculați o estimare AUC pe baza ROC netezită):
(Linia roșie este ROC originală, iar linia neagră este ROC netezită. Rețineți, de asemenea, raportul de aspect implicit 1: 1. Este logic să folosiți acest lucru, deoarece atât sensibilitatea, cât și specificitatea au 0-1 .)
AUC estimat din ROC netezit este de 0,9107, similar, dar puțin mai mare decât AUC din ROC netezit (dacă arătați un cifra, puteți vedea cu ușurință de ce este mai mare). (Deși avem într-adevăr prea puține valori distincte posibile ale rezultatelor testului pentru a calcula o AUC netedă).
Pachetul rms
Pachetul rms
al lui Harrell poate calcula diferite statistici de concordanță conexe utilizând funcția rcorr.cens()
. C Index
din ieșirea sa este AUC:
> library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711
Pachetul caTools
În sfârșit, avem pachetul caTools
și funcția sa colAUC()
. Are câteva avantaje față de alte pachete (în principal viteza și capacitatea de a lucra cu date multidimensionale – consultați ?colAUC
) care pot uneori să vă fie de ajutor.Dar, desigur, oferă același răspuns pe care l-am calculat mereu:
library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711
Cuvinte finale
Mulți oameni par să creadă că ASC ne spune cât de „bine” un test este. Și unii oameni cred că ASC este probabilitatea ca testul să clasifice corect un pacient. Nu este nu . După cum puteți vedea din exemplul și calculele de mai sus, ASC ne spune ceva despre o familie de teste, câte un test pentru fiecare posibilă limită.
Și ASC se calculează pe baza tăieturi pe care nu le-ar folosi niciodată în practică. De ce ar trebui să ne pese de sensibilitatea și specificitatea valorilor limită „fără sens”? Totuși, pe asta se bazează (parțial) ASC. (Bineînțeles, dacă ASC este foarte aproape de 1, aproape fiecare test posibil va avea o mare putere discriminatorie și am fi cu toții foarte fericiți.)
„Normalul aleatoriu” –Interpretarea anormală a perechilor ASC este plăcută (și poate fi extinsă, de exemplu, la modelele de supraviețuire, unde vedem dacă este persoana cu cel mai mare risc (relativ) care moare cel mai devreme). Dar nu s-ar folosi niciodată în practică. Este un caz rar în care cineva știe are unul sănătos și unul bolnav, nu știe care este cel bolnav și trebuie decide care dintre ele să trateze. (În orice caz, decizia este ușoară; tratați-o pe cea cu cel mai mare risc estimat.)
Deci, cred că studierea curentei curbei ROC va fi mai utilă decât simpla examinare măsura sumară AUC. Și dacă utilizați ROC împreună cu (estimări ale) costurilor falselor pozitive și falselor negative, împreună cu ratele de bază ale a ceea ce studiați, puteți ajunge undeva.
Rețineți, de asemenea, că ASC măsoară numai discriminarea , nu calibrarea. Adică, măsoară dacă puteți discrimina între două persoane (una bolnavă și una sănătoasă), pe baza scorului de risc. Pentru aceasta, se uită doar la valorile de risc relative (sau se clasează, dacă doriți, cf. interpretarea testului Wilcoxon – Mann – Whitney), nu la cele absolute, pe care ar trebui să le De exemplu, dacă împărțiți fiecare estimare a riscului de modelul logistic la 2, veți obține exact aceleași ASC (și ROC).
Atunci când evaluați un model de risc, calibrarea este, de asemenea, foarte importantă. Pentru a examina acest lucru, veți privi toți pacienții cu un scor de risc de aproximativ, de exemplu, 0,7 și veți vedea dacă aproximativ 70% dintre aceștia erau de fapt bolnavi. Faceți acest lucru pentru fiecare scor de risc posibil (eventual folosind un fel de netezire / regresie locală). Graficează rezultatele și vei obține o măsură grafică a calibrării .
Dacă ai un model cu ambele calibrare bună și discriminare bună, atunci tu începe să ai un model bun. 🙂
Comentarii
Răspuns
Aruncați o privire la această întrebare: Înțelegerea curbei ROC
Iată cum să construiți o curbă ROC (din acea întrebare):
Desenarea curbei ROC
cu un set de date procesat de dvs. clasificator de clasificare
- clasificați exemplele de testare pe scorul descrescător
- începeți cu $ (0, 0) $
- pentru fiecare exemplu $ x $ (în ordinea descrescătoare)
- dacă $ x $ este pozitiv, mutați $ 1 / \ text {pos} $ în sus
- dacă $ x $ este negativ, mutați $ 1 / \ text {neg} $ dreapta
unde $ \ text {pos} $ și $ \ text {neg} $ sunt fracțiile de exemple pozitive și respectiv negative.
Puteți utiliza această idee pentru a calcula manual AUC ROC utilizând următorul 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
Această frumoasă imagine animată de gif ar trebui să ilustreze acest lucru proces mai clar
Comentarii
- Mulțumesc @Alexey Grigorev, acesta este un aspect vizual extraordinar și se va dovedi probabil util în viitor! +1
- Vă rog să explicați puțin despre ” fracții de exemple pozitive și negative „, vă referiți la cea mai mică valoare unitate a două axe?
- @Allan Ruin:
pos
aici înseamnă numărul de date pozitive. Să spunem că aveți 20 de puncte de date, în care 11 puncte sunt 1. Deci, când desenăm graficul, avem un dreptunghi 11×9 (înălțime x lățime). Alexey Grigorev a scalat, dar a lăsat-o așa cum dorește ‘. Acum, trebuie doar să mutați 1 pe grafic la fiecare pas.
Răspuns
Postarea lui Karl are multe de informații excelente, dar nu am văzut încă în ultimii 20 de ani un exemplu de curbă ROC care să schimbe gândirea oricui într-o direcție bună. Singura valoare a unei curbe ROC în umila mea părere este că aria sa se întâmplă să fie egală cu o probabilitate de concordanță foarte utilă. Curba ROC în sine îl tentează pe cititor să folosească limite, ceea ce este o practică statistică proastă.
În ceea ce privește calcularea manuală a indexului $ c $, faceți un grafic cu $ Y = 0,1 $ pe $ x $ -axa și predictorul continuu sau probabilitatea prezisă ca $ Y = 1 $ pe $ y $ -axa. Dacă conectați fiecare punct cu $ Y = 0 $ cu fiecare punct cu $ Y = 1 $, proporția liniilor care au o pantă pozitivă este probabilitatea de concordanță.
Orice măsură care are un numitor de $ n $ în această setare sunt reguli incorecte de notare a preciziei și ar trebui evitate. Aceasta include proporția clasificată corect, sensibilitatea și specificitatea.
Pentru funcția R Hmisc
pachet rcorr.cens
, tipăriți întregul rezultat pentru a vedea mai multe informații, în special o eroare standard.
Comentarii
- Mulțumesc, @Frank Harell, îți apreciez perspectiva. Pur și simplu folosesc statistica c ca probabilitate de concordanță, deoarece nu îmi plac ‘. Vă mulțumim din nou!
Răspuns
Iată o alternativă la modul natural de calcul al ASC prin simpla utilizare a trapezoidalului regula pentru a obține zona de sub curba ROC.
ASC este egală cu probabilitatea ca o observație pozitivă eșantionată să aibă o probabilitate estimată (de a fi pozitivă) mai mare decât o observație negativă eșantionată aleatoriu. Puteți utiliza acest lucru pentru a calcula AUC destul de ușor în orice limbaj de programare parcurgând toate combinațiile perechi de observații pozitive și negative. De asemenea, ați putea eșantiona observații în mod aleatoriu dacă dimensiunea eșantionului era prea mare. Dacă doriți să calculați AUC folosind stilou și hârtie, este posibil să nu fie cea mai bună abordare decât dacă aveți un eșantion foarte mic / mult timp. De exemplu, în 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
Putem verifica folosind pachetul pROC
:
library(pROC) auc(y, probs) Area under the curve: 0.6287
Utilizarea eșantionării aleatorii:
mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896
Răspuns
- Aveți adevărată valoare pentru observații.
- Calculați probabilitatea posterioară și apoi clasificați observațiile după această probabilitate.
- Presupunând o probabilitate de limită de $ P $ și numărul de observații $ N $:
$$ \ frac {\ text {Suma de ranguri adevărate} -0.5PN (PN + 1)} { PN (N-PN)} $$
Comentarii
- @ user73455 … 1) Da, am adevărata valoare pentru observații. 2) Probabilitatea posterioară este sinonimă cu probabilitățile prezise pentru fiecare dintre observații? 3) Înțeles; totuși, ce este ” Suma de ranguri adevărate ” și cum se calculează această valoare? Poate că un exemplu te-ar ajuta să explici mai bine acest răspuns? Vă mulțumim!
sens=c(sens,0); omspec=c(omspec,0)
, nu ar trebui să ‘ să fiesens=c(0, sens); omspec=c(0, omspec)
? Se complotează corect cu0
principal, dar nu așa cum este în prezent în răspuns.sens
cât șiomspec
înainte de complot).