I zajímám se o ruční výpočet plochy pod křivkou (AUC) nebo c-statistiky pro binární logistický regresní model.

Například v validační datová sada, mám skutečnou hodnotu pro závislou proměnnou, retence (1 = zachována; 0 = nezachována), stejně jako předpokládaný stav uchování pro každé pozorování generované mojí regresní analýzou pomocí modelu, který byl vytvořen pomocí školení set (to bude v rozsahu od 0 do 1).

Moje počáteční myšlenky byly identifikovat „správný“ počet klasifikací modelu a jednoduše rozdělit počet „správných“ pozorování počtem celkových pozorování, která se mají vypočítat statistika c. „Správně“, pokud je skutečný retenční stav pozorování = 1 a předpokládaný retenční stav je> 0,5, pak jde o „správnou“ klasifikaci. Navíc pokud je skutečný stav uchování pozorování = 0 a předpokládaný stav uchování < 0,5, pak jde také o „správnou“ klasifikaci. Předpokládám, že ke „kravatě“ dojde, když je předpovězená hodnota = 0,5, ale tento jev se v mém souboru údajů o ověření nevyskytuje. Na druhou stranu by „nesprávnou“ klasifikací bylo, kdyby skutečný stav uchování pozorování = 1 a předpokládaný stav uchování byl < 0,5, nebo pokud by byl skutečný stav uchování výsledku = 0 a předpokládaný stav uchování je> 0,5. Jsem si vědom TP, FP, FN, TN, ale nevím, jak vypočítat statistiku c vzhledem k těmto informacím.

Odpověď

Doporučil bych Hanleyho & McNeilův dokument z roku 1982 Význam a využití oblasti pod provozní charakteristikou přijímače (ROC ) křivka .

Příklad

Mají následující tabulku stavu onemocnění a výsledek testu (odpovídající například odhadovanému riziku z logistického modelu). První číslo vpravo je počet pacientů se skutečným stavem choroby „normální“ a druhé číslo je počet pacientů se skutečným stavem choroby „abnormální“:

(1) Rozhodně normální: 33/3
(2) Pravděpodobně normální: 6/2
(3) Otázné: 6/2
(4) Pravděpodobně neobvyklé: 11/11
(5) Rozhodně abnormální: 2/33

Celkově tedy existuje 58 „normálních“ pacientů a „51“ abnormálních. Vidíme, že když je prediktor 1, „Rozhodně normální“, pacient je obvykle normální (platí pro 33 z 36 pacientů), a když je 5, „Určitě abnormální“, pacienti jsou obvykle abnormální (platí pro 33 z 35 pacientů), takže prediktor dává smysl. Jak bychom však měli soudit pacienta se skóre 2, 3 nebo 4? To, co jsme nastavili jako mezní hodnotu pro hodnocení pacientů jako abnormálních nebo normálních, abychom určili citlivost a specificitu výsledného testu.

Citlivost a specificita

Můžeme vypočítat odhad citlivost a specificita pro různé mezní hodnoty. (Od nynějška budu psát „citlivost“ a „specifičnost“, takže odhadovaná povaha hodnot bude implicitní.)

Pokud zvolíme mezní hodnotu, abychom mohli klasifikovat všechny pacienti jako abnormální, bez ohledu na to, co říkají jejich výsledky testů (tj. zvolíme mezní hodnotu 1+), získáme citlivost 51/51 = 1. Specifičnost bude 0/58 = 0. Ne zní to tak dobře.

Dobře, tak pojďme zvolit méně přísný limit. Pacienty klasifikujeme jako abnormální, pouze pokud mají výsledek testu 2 nebo vyšší. Poté nám chybí 3 abnormální pacienti a máme citlivost 48/51 = 0,94. Ale máme mnohem vyšší specificitu, 33/58 = 0,57.

Nyní můžeme pokračovat a zvolit různé mezní hodnoty (3, 4, 5,> 5). (V posledním případě nebudeme žádného pacienta klasifikovat jako abnormálního, i když bude mít nejvyšší možné skóre testu 5.)

ROC křivka

Pokud to uděláme pro všechny možné mezní hodnoty a vykreslíme citlivost proti 1 minus specifičnost, dostaneme ROC křivku. Můžeme použít následující R kód:

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

Výstup je:

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

Můžeme vypočítat různé statistiky:

 ( 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  

A pomocí toho můžeme vykreslit (odhadovanou) ROC křivku:

 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)  

křivka AUC

Ruční výpočet AUC

Můžeme velmi snadno vypočítat plochu pod křivkou ROC pomocí vzorce pro plochu lichoběžníku:

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

Výsledkem je 0,8931711.

Shoda opatření

AUC lze také chápat jako opatření shody.Pokud vezmeme všechny možné páry pacientů, kde jeden je normální a druhý je abnormální, můžeme vypočítat, jak často je to abnormální, který má nejvyšší (nejvíce „abnormálně vypadající“) výsledek testu (pokud mají stejnou hodnotu, počítáme to jako „poloviční vítězství“):

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

Odpověď je opět 0,8931711, což je oblast pod křivkou ROC. Bude tomu tak vždy.

Grafické zobrazení shody

Jak zdůraznil Harrell ve své odpovědi, má to také grafickou interpretaci. Vynesme skóre testu (odhad rizika) na ose y a skutečný stav onemocnění na ose x (zde s určitým chvěním, abychom ukázali překrývající se body):

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

Bodový graf skóre rizika proti skutečné nemoci status.

Pojďme nyní nakreslit čáru mezi každým bodem vlevo („normální“ pacient) a každým bodem vpravo („abnormální“ pacient). Podíl linií s kladným sklonem (tj. Podíl shodných párů) je index shody (ploché čáry se počítají jako „50% shoda“).

Je trochu obtížné vizualizovat skutečné řádky pro tento příklad, kvůli počtu vazeb (skóre se stejným rizikem), ale s určitým chvěním a průhledností můžeme získat rozumnou zápletku:

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

Bodový graf skóre rizika proti skutečnému stavu choroby, s čáry mezi všemi možnými pozorovacími páry.

Vidíme, že většina čar stoupá vzhůru, takže index shody bude vysoký. Vidíme také příspěvek do indexu od každého typu pozorovacího páru. Většina z nich pochází od normálních pacientů s rizikovým skóre 1 spárovaných s abnormálními pacienty s rizikovým skóre 5 (1–5 párů), ale dost také z 1–4 párů a 4–5 párů. A je velmi snadné vypočítat skutečný index shody na základě definice sklonu:

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

Odpověď je opět 0,8931711, tj. AUC.

Wilcoxon – Mann – Whitneyův test

Existuje úzké spojení mezi konkordančním měřítkem a Wilcoxon – Mann – Whitney test. Ve skutečnosti druhý testuje, zda je pravděpodobnost shody (tj. Že abnormální pacient v náhodném páru normální – abnormální bude mít nejvíce „abnormálně vypadající“ výsledek testu) je přesně 0,5. A jeho testovací statistika je pouze jednoduchá transformace odhadované pravděpodobnosti konkordance:

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

Statistika testu (W = 2642) počítá počet shodných párů. Vydělíme-li to počtem možných párů, dostaneme familiární číslo:

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

Ano, je to 0,8931711, oblast pod křivkou ROC.

Snadnější způsoby výpočtu AUC (v R)

Pojďme si však usnadnit život. Existují různé balíčky, které pro nás automaticky vypočítají AUC.

Balíček Epi

Balíček Epi vytváří pěknou ROC křivku s různými vložené statistiky (včetně AUC):

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

ROC křivka z balíčku Epi

Balíček pROC

Také se mi líbí balíček pROC, protože může vyhladit odhad ROC (a vypočítat odhad AUC na základě vyhlazeného ROC):

ROC křivka (nehladká a vyhlazená) z balíčku pROC

(Červená čára je původní ROC a černá čára je vyhlazený ROC. Všimněte si také výchozího poměru stran 1: 1. Má smysl to použít, protože jak citlivost, tak specifičnost má 0–1 rozsah.)

Odhadovaná AUC z vyhlazeného ROC je 0,9107, podobná, ale mírně větší než AUC z nehladeného ROC (pokud se podíváte na Na obrázku snadno zjistíte, proč je větší). (I když ve skutečnosti máme příliš málo možných odlišných hodnot výsledků testu, abychom mohli vypočítat plynulé AUC).

Balíček RMS

Balíček Harrell rms může vypočítat různé související statistiky konkordance pomocí funkce rcorr.cens(). C Index v jeho výstupu je AUC:

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

Balíček caTools

Konečně máme balíček caTools a jeho colAUC() funkci. Oproti jiným balíčkům má několik výhod (zejména rychlost a schopnost pracovat s vícerozměrnými daty – viz ?colAUC), které někdy mohou být užitečné.Ale samozřejmě dává stejnou odpověď, jakou jsme počítali znovu a znovu:

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

ROC křivka z balíčku caTools

Závěrečná slova

Zdá se, že mnoho lidí si myslí, že AUC nám říká, jak „dobrá“ test je. A někteří lidé si myslí, že AUC je pravděpodobnost, že test správně klasifikuje pacienta. Není to ne . Jak můžete vidět z výše uvedeného příkladu a výpočtů, AUC nám říká něco o rodině testů, jeden test pro každou možnou mez.

A AUC se počítá na základě cutoffs by člověk nikdy nepoužil v praxi. Proč bychom se měli starat o citlivost a specifičnost „nesmyslných“ hraničních hodnot? Na tom je (částečně) AUC stále založen. (Samozřejmě, pokud je AUC velmi blízko 1, téměř každý možný test bude mít velkou diskriminační sílu a všichni bychom byli velmi šťastní.)

„Náhodný normální –Abnormální párová interpretace AUC je hezká (a lze ji rozšířit, například na modely přežití, kde vidíme, zda jde o osobu s nejvyšším (relativním) nebezpečím, která zemře nejdříve). Ale člověk by to nikdy nepoužil v praxi. Je to vzácný případ, kdy jeden jednoho zdravého a jednoho nemocného, neví, která osoba je nemocný, a musí rozhodnout, s kým z nich zacházet. (V každém případě je rozhodnutí snadné; zacházejte s tím s nejvyšším odhadovaným rizikem.)

Takže si myslím, že studium skutečné křivky ROC bude užitečnější než pouhý pohled na souhrnné opatření AUC. A pokud použijete ROC společně s (odhady) nákladů falešných pozitiv a falešných negativů, spolu se základními sazbami toho, co studujete, můžete se někam dostat.

Pamatujte, že AUC měří pouze diskriminaci , nikoli kalibraci. To znamená, že měří, zda můžete rozlišovat mezi dvěma osobami (jednou nemocnou a jednou zdravou) na základě skóre rizika. Z tohoto důvodu se dívá pouze na relativní rizikové hodnoty (nebo hodnotíte, chcete-li, viz výklad testu Wilcoxon – Mann – Whitney), nikoli na absolutní hodnoty, které byste měli mít o ně zájem. Pokud například vydělíte každý odhad rizika z vašeho logistického modelu číslem 2, získáte přesně stejnou AUC (a ROC).

Při hodnocení modelu rizika kalibrace je také velmi důležitá. Chcete-li to prozkoumat, podíváte se na všechny pacienty s rizikovým skóre kolem, např. 0,7, a uvidíte, zda přibližně 70% z nich skutečně onemocnělo. Udělejte to pro každé možné skóre rizika (případně pomocí nějakého vyhlazení / lokální regrese). Umístěte výsledky do grafu a získáte grafickou míru kalibrace .

Pokud máte model s dobrou dobrou kalibrací a dobrou diskriminací, pak začít mít dobrý model. 🙂

Komentáře

  • Děkuji vám, @Karl Ove Hufthammer, toto je ta nejdůkladnější odpověď, jakou jsem kdy dostal. Zvláště oceňuji vaši sekci “ Závěrečná slova „. Skvělá práce! Ještě jednou děkujeme!
  • Děkuji vám za tuto podrobnou odpověď. Pracuji s datovou sadou, kde Epi :: ROC () v2.2.6 je přesvědčen, že AUC je 1,62 (ne to není mentalistická studie), ale podle ROC věřím mnohem více v 0,56, že výsledky výše uvedeného kódu in.
  • Myslím, že v sens=c(sens,0); omspec=c(omspec,0) je malá chyba, neměla by to být ‚ sens=c(0, sens); omspec=c(0, omspec)? Vykresluje správně s předními 0, ale ne tak, jak je aktuálně v odpovědi.
  • Ne, aktuální definice je, AFAICS, right, @steveb, a vede ke správnému grafu. Myslím, že to, co je možná matoucí, je, že křivka ROC je kreslena z zprava doleva (tj. Z pravého horního rohu do levého dolního rohu), nikoli z zleva do správně , jako většina pozemků. To je jen výsledek toho, jak jsem definoval proměnné; stejně dobře by se to dalo vykreslit zleva doprava (obrácením sens a omspec před vykreslením).

Odpověď

Podívejte se na tuto otázku: Porozumění křivce ROC

Zde je postup, jak vytvořit ROC křivku (z této otázky):

Kreslení ROC křivky

vzhledem k datové sadě zpracované vaším klasifikátor hodnocení

  • příklady testů na snížení skóre
  • začínají v $ (0, 0) $
  • pro každý příklad $ x $ (v sestupné pořadí)
    • pokud je $ x $ kladné, posuňte $ 1 / \ text {pos} $ nahoru
    • pokud je $ x $ záporné, posuňte $ 1 / \ text {neg} $ doprava

kde $ \ text {pos} $ a $ \ text {neg} $ jsou zlomky pozitivních a negativních příkladů.

Tuto myšlenku můžete použít k ručnímu výpočtu AUC ROC pomocí následujícího algoritmu:

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 

Tento pěkný obrázek ve formátu GIF by to měl ilustrovat proces jasnější

budování křivky

Komentáře

  • Díky @Alexey Grigoreve, toto je skvělý vizuál a pravděpodobně se v budoucnu osvědčí! +1
  • Mohl byste prosím trochu vysvětlit “ zlomky pozitivních a negativních příkladů „, máte na mysli nejmenší jednotková hodnota dvou os?
  • @ Allan Ruin: pos zde znamená počet kladných dat. Řekněme, že máte 20 datových bodů, ve kterých je 11 bodů 1. Takže při kreslení grafu máme obdélník 11×9 (výška x šířka). Alexey Grigorev udělal měřítko, ale pokud to chcete, nechte to ‚ s. Nyní stačí v každém kroku v grafu přesunout 1.

Odpovědět

Karlův příspěvek má mnoho vynikajících informací. Ale za posledních 20 let jsem ještě neviděl příklad ROC křivky, která změnila kohokoli myšlení dobrým směrem. Jedinou hodnotou ROC křivky podle mého skromného názoru je, že její plocha se shoduje s velmi užitečnou pravděpodobností shody. Samotná křivka ROC láká čtenáře k použití mezních hodnot, což je špatná statistická praxe.

Pokud jde o manuální výpočet $ c $ -indexu, vytvořte graf s $ Y = 0,1 $ na $ x $ -axis a spojitý prediktor nebo předpokládaná pravděpodobnost, že $ Y = 1 $ na $ y $ -axis. Pokud spojíte každý bod s $ Y = 0 $ s každým bodem s $ Y = 1 $, podíl linií, které mají kladný sklon, je pravděpodobnost shody.

Jakákoli opatření, která mají jmenovatele $ n $ v tomto nastavení jsou nesprávná pravidla bodování přesnosti a je třeba se jim vyhnout. Patří sem správně klasifikovaný podíl, citlivost a specificita.

Pro funkci R Hmisc balíčku rcorr.cens vytiskněte celý výsledek, abyste viděli více informací, zejména standardní chyba.

Komentáře

  • Děkuji, @Frank Harell, vážím si vaší perspektivy. Jednoduše používám c-statistiku jako pravděpodobnost konkordance, protože se mi ‚ nelíbí mezní hodnoty. Ještě jednou děkujeme!

Odpověď

Zde je alternativa k přirozenému způsobu výpočtu AUC jednoduše pomocí lichoběžníkového tvaru pravidlo, abyste dostali oblast pod ROC křivkou.

AUC se rovná pravděpodobnosti, že náhodně vzorkované pozitivní pozorování má předpovězenou pravděpodobnost (pozitivního) větší než náhodně vzorkované negativní pozorování. Tímto způsobem můžete poměrně snadno vypočítat AUC v jakémkoli programovacím jazyce procházením všech párových kombinací pozitivních a negativních pozorování. Pokud byla velikost vzorku příliš velká, můžete také náhodně vybrat pozorování. Pokud chcete vypočítat AUC pomocí pera a papíru, nemusí to být nejlepší přístup, pokud nemáte velmi malý vzorek / spoustu času. Například v 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 

Můžeme ověřit pomocí balíčku pROC:

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

Použití náhodného vzorkování:

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

Odpověď

  1. Máte skutečnou hodnotu pro pozorování.
  2. Vypočítejte zadní pravděpodobnost a poté seřazte pozorování podle této pravděpodobnosti.
  3. Za předpokladu pravděpodobnosti cut-off $ P $ a počtu pozorování $ N $:
    $$ \ frac {\ text {Součet skutečných hodnocení} -0,5PN (PN + 1)} { PN (N-PN)} $$

Komentáře

  • @ user73455 … 1) Ano, mám skutečnou hodnotu pro pozorování. 2) Je zadní pravděpodobnost synonymem pro předpokládané pravděpodobnosti pro každé z pozorování? 3) Rozumím; Co je však “ Součet skutečných řad “ a jak lze tuto hodnotu vypočítat? Možná by vám příklad pomohl důkladněji vysvětlit tuto odpověď? Děkuji!

Napsat komentář

Vaše e-mailová adresa nebude zveřejněna. Vyžadované informace jsou označeny *