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)  

curba AUC

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

Scatter diagramă de scor de risc împotriva adevăratei boli status.

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

Scatter diagramă a scorului de risc împotriva stării adevărate a bolii, cu linii între toate perechile de observație posibile.

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"  

Curba ROC din pachetul Epi

Pachetul pROC

Îmi place și pachetul pROC, deoarece poate neteziți estimarea ROC (și calculați o estimare AUC pe baza ROC netezită):

Curba ROC (netezită și netezită) din pachetul pROC

(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  

Curba ROC din pachetul caTools

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

  • Mulțumesc, @Karl Ove Hufthammer, acesta este cel mai amănunțit răspuns pe care l-am primit vreodată. Apreciez în special secțiunea dvs. ” Cuvinte finale „. Muncă excelentă! Vă mulțumim din nou!
  • Vă mulțumesc foarte mult pentru acest răspuns detaliat. Lucrez cu un set de date în care Epi :: ROC () v2.2.6 este convins că ASC este 1,62 (nu, nu este un studiu mentalist), dar conform ROC, cred mult mai mult în 0,56 că rezultatul codului de mai sus în.
  • Cred că există o mică eroare în sens=c(sens,0); omspec=c(omspec,0), nu ar trebui să ‘ să fie sens=c(0, sens); omspec=c(0, omspec)? Se complotează corect cu 0 principal, dar nu așa cum este în prezent în răspuns.
  • Nu, definiția actuală este, AFAICS, corectă, @steveb, și are ca rezultat un complot corect. Cred că ceea ce este probabil confuz este că curba ROC este trasată de la dreapta la stânga (adică de la colțul din dreapta sus la colțul din stânga jos), nu de la stânga la dreapta , la fel ca majoritatea parcelelor. Acesta este doar rezultatul modului în care am definit variabilele; s-ar fi putut la fel de bine să-l complotăm de la stânga la dreapta (inversând atât sens cât și omspec înainte de complot).

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

construirea curbei

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

  1. Aveți adevărată valoare pentru observații.
  2. Calculați probabilitatea posterioară și apoi clasificați observațiile după această probabilitate.
  3. 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!

Lasă un răspuns

Adresa ta de email nu va fi publicată. Câmpurile obligatorii sunt marcate cu *