I jestem zainteresowany ręcznym obliczeniem obszaru pod krzywą (AUC) lub statystyką c dla binarnego modelu regresji logistycznej.
Na przykład w zestaw danych walidacyjnych, mam prawdziwą wartość zmiennej zależnej, retencję (1 = zachowane; 0 = nie zachowane), a także przewidywany stan retencji dla każdej obserwacji wygenerowanej przez moją analizę regresji przy użyciu modelu, który został zbudowany przy użyciu treningu zestaw (zakres od 0 do 1).
Moje początkowe przemyślenia dotyczyły określenia „prawidłowej” liczby klasyfikacji modeli i po prostu podzielenia liczby „poprawnych” obserwacji przez liczbę wszystkich obserwacji do obliczenia statystyka c. Przez „prawidłowy”, jeśli rzeczywisty stan zachowania obserwacji = 1, a przewidywany stan retencji wynosi> 0,5, oznacza to „prawidłową” klasyfikację. Ponadto, jeśli prawdziwy stan zachowania obserwacji = 0, a przewidywany stan przechowywania wynosi < 0,5, to również jest to „poprawna” klasyfikacja. Zakładam, że „remis” wystąpiłby, gdy przewidywana wartość = 0,5, ale to zjawisko nie występuje w moim zbiorze danych do walidacji. Z drugiej strony „niepoprawne” klasyfikacje miałyby miejsce, gdyby rzeczywisty stan zachowania obserwacji = 1, a przewidywany stan przechowywania to < 0,5 lub jeśli prawdziwy stan zachowania wyniku = 0, a przewidywany stan przechowywania wynosi> 0,5. Znam TP, FP, FN, TN, ale nie wiem, jak obliczyć statystykę c na podstawie tych informacji.
Odpowiedź
Polecam & artykuł McNeila z 1982 roku Znaczenie i wykorzystanie obszaru pod charakterystyką operacyjną odbiornika (ROC ) krzywa .
Przykład
Mają poniższą tabelę statusu choroby i wyniku testu (odpowiadającą, na przykład, oszacowanemu ryzyku z modelu logistycznego). Pierwsza liczba po prawej stronie to liczba pacjentów z prawdziwym stanem chorobowym „normalny”, a druga liczba to liczba pacjentów z prawdziwym stanem choroby „nieprawidłowy”:
(1) Zdecydowanie normalne: 33/3
(2) Prawdopodobnie normalne: 6/2
(3) Wątpliwe: 6/2
(4) Prawdopodobnie nienormalne: 11/11
(5) Zdecydowanie nienormalny: 2/33
W sumie jest 58 „normalnych” pacjentów i „51” nienormalnych. Widzimy, że gdy predyktorem jest 1, „ zdecydowanie normalny , pacjent jest zwykle normalny (prawda dla 33 z 36 pacjentów), a kiedy wynosi 5, „ zdecydowanie nieprawidłowy , pacjenci są zwykle anormalni (prawda dla 33 z 36 pacjentów). 35 pacjentów), więc predyktor ma sens. Ale jak powinniśmy oceniać pacjenta z wynikiem 2, 3 lub 4? To, co ustaliliśmy dla oceny pacjenta jako nieprawidłowego lub normalnego, określa czułość i swoistość wynikowego testu.
Czułość i swoistość
Możemy obliczyć oszacowaną czułość i swoistość dla różnych wartości odcięcia. (Od teraz będę po prostu pisać „wrażliwość” i „specyficzność”, pozwalając, aby szacunkowa natura wartości była niejawna).
Jeśli wybierzemy wartość graniczną, aby sklasyfikować wszystkie pacjenci jako anormalni, bez względu na to, co mówią ich wyniki badań (tj. wybierzemy wartość odcięcia 1+), uzyskamy czułość 51/51 = 1. Specyficzność wyniesie 0/58 = 0. Nie brzmią tak dobrze.
OK, więc wybierzmy mniej ścisłe ograniczenie. Klasyfikujemy pacjentów jako anormalnych tylko wtedy, gdy mają wynik testu 2 lub wyższy. Następnie pomijamy 3 nieprawidłowych pacjentów i mamy czułość 48/51 = 0,94. Ale mamy znacznie zwiększoną specyficzność, 33/58 = 0,57.
Możemy teraz kontynuować, wybierając różne wartości odcięcia (3, 4, 5,> 5). (W ostatnim przypadku nie zaklasyfikujemy żadnego pacjenta jako anormalnego, nawet jeśli ma on najwyższy możliwy wynik testu 5)
Krzywa ROC
Jeśli zrobimy to dla wszystkich możliwych wartości odcięcia i wykreślimy czułość względem 1 minus swoistość, otrzymamy krzywą ROC. Możemy użyć następującego kodu 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)) )
Wynik to:
testres truestat 1 2 3 4 5 0 33 6 6 11 2 1 3 2 2 11 33
Możemy obliczyć różne statystyki:
( 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 używając tego, możemy wykreślić (szacunkową) krzywą 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)
Ręczne obliczanie AUC
Możemy bardzo łatwo obliczyć pole pod krzywą ROC, używając wzoru na pole powierzchni trapezu:
height = (sens[-1]+sens[-length(sens)])/2 width = -diff(omspec) # = diff(rev(omspec)) sum(height*width)
Wynik to 0,8931711.
Miara zgodności
AUC można również postrzegać jako miarę zgodności.Jeśli weźmiemy wszystkie możliwe pary pacjentów, z których jeden jest normalny, a drugi nieprawidłowy, możemy obliczyć, jak często to ten nienormalny ma najwyższy (najbardziej nienormalnie wyglądający) wynik testu (jeśli mają tę samą wartość, uważamy to za „połowę zwycięstwa”):
o = outer(abnorm, norm, "-") mean((o>0) + .5*(o==0))
Odpowiedź to znowu 0,8931711, obszar pod krzywą ROC. Tak będzie zawsze.
Graficzne przedstawienie zgodności
Jak zauważył Harrell w swojej odpowiedzi, ma to również graficzną interpretację. Wykreślmy wynik testu (oszacowanie ryzyka) na osi y i prawdziwy stan choroby na osi x (tutaj z pewnymi wahaniami, aby pokazać nakładające się punkty):
plot(jitter(truestat,.2), jitter(testres,.8), las=1, xlab="True disease status", ylab="Test score")
Narysujmy teraz linię między każdym punktem po lewej stronie („normalny” pacjent) a każdym punktem po prawej stronie („nienormalny” pacjent). Odsetek prostych o dodatnim nachyleniu (tj. Proporcja zgodnych par) to wskaźnik zgodności (linie płaskie liczą się jako „zgodność 50%”).
Trochę trudno jest wyobrazić sobie rzeczywiste linie dla tego przykładu ze względu na liczbę remisów (równy wynik ryzyka), ale przy pewnych fluktuacjach i przejrzystości możemy uzyskać rozsądny wykres:
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))
Widzimy, że większość linii jest nachylona w górę, więc indeks zgodności będzie wysoki. Widzimy również udział w indeksie każdego typu pary obserwacji. Większość z nich pochodzi od normalnych pacjentów z wynikiem ryzyka 1 w parze z pacjentami nienormalnymi z wynikiem ryzyka 5 (1–5 par), ale całkiem sporo pochodzi również od 1–4 par i 4–5 par. Obliczenie aktualnego wskaźnika zgodności na podstawie definicji nachylenia jest bardzo łatwe:
d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm)) mean((d$slope > 0) + .5*(d$slope==0))
Odpowiedź to znowu 0,8931711, tj. AUC.
Test Wilcoxona – Manna – Whitneya
Istnieje ścisły związek między miarą zgodności a miarą Wilcoxona – Manna – Whitneya. test. W rzeczywistości to ostatnie sprawdza, czy prawdopodobieństwo zgodności (tj. Czy to nieprawidłowy pacjent w losowej parze normalnie-nieprawidłowej ma najbardziej „nienormalnie wyglądający” wynik testu) wynosi dokładnie 0,5. Statystyka testowa to po prostu prosta transformacja szacowanego prawdopodobieństwa zgodności:
> ( 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
Statystyka testowa (W = 2642
) zlicza liczbę zgodnych par. Jeśli podzielimy to przez liczbę możliwych par, otrzymamy podobną liczbę:
w = wi$statistic w/(length(abnorm)*length(norm))
Tak, to 0,8931711, obszar pod krzywą ROC.
Łatwiejsze sposoby obliczania AUC (w R)
Ułatwmy sobie życie. Istnieją różne pakiety, które automatycznie obliczają AUC dla nas.
Pakiet Epi
Pakiet Epi
tworzy ładną krzywą ROC z różnymi osadzone statystyki (w tym AUC):
library(Epi) ROC(testres, truestat) # also try adding plot="sp"
Pakiet pROC
Podoba mi się również pakiet pROC
, ponieważ może wygładź oszacowanie ROC (i oblicz oszacowanie AUC na podstawie wygładzonego ROC):
(Czerwona linia to oryginalna ROC, a czarna to wygładzona ROC. Zwróć także uwagę na domyślny współczynnik proporcji 1: 1. Warto z tego skorzystać, ponieważ zarówno czułość, jak i swoistość mają wartość 0–1 zakres.)
Szacunkowa wartość AUC z wygładzonego ROC wynosi 0,9107, podobna do, ale nieco większa niż AUC z niewygładzonego ROC (jeśli spojrzysz na na figurze, możesz łatwo zobaczyć, dlaczego jest większa). (Chociaż naprawdę mamy zbyt mało możliwych odrębnych wartości wyników testu, aby obliczyć gładkie AUC).
Pakiet rms
rms
pakiet Harrella może obliczyć różne powiązane statystyki zgodności za pomocą funkcji rcorr.cens()
. C Index
na wyjściu to AUC:
> library(rms) > rcorr.cens(testres,truestat)[1] C Index 0.8931711
Pakiet caTools
Wreszcie mamy pakiet caTools
i jego funkcję colAUC()
. Ma kilka zalet w porównaniu z innymi pakietami (głównie szybkość i możliwość pracy z danymi wielowymiarowymi – patrz ?colAUC
), które mogą czasami być pomocne.Ale oczywiście daje tę samą odpowiedź, którą obliczaliśmy wielokrotnie:
library(caTools) colAUC(testres, truestat, plotROC=TRUE) [,1] 0 vs. 1 0.8931711
Ostatnie słowa
Wiele osób wydaje się myśleć, że AUC mówi nam, jak „dobrze” test jest. Niektórzy uważają, że AUC to prawdopodobieństwo, że test poprawnie sklasyfikuje pacjenta. To , a nie . Jak widać na powyższym przykładzie i obliczeniach, AUC mówi nam coś o rodzinie testów, po jednym teście dla każdego możliwego odcięcia.
AUC jest obliczane na podstawie odcięcia, których nigdy by się nie użyło w praktyce. Dlaczego powinniśmy dbać o wrażliwość i specyficzność „bezsensownych” wartości odcięcia? Jednak na tym (częściowo) opiera się wartość AUC. (Oczywiście, jeśli AUC jest bardzo bliskie 1, prawie każdy możliwy test będzie miał wielką moc rozróżniającą i wszyscy bylibyśmy bardzo szczęśliwi.)
– interpretacja AUC w przypadku nieprawidłowej pary jest dobra (i można ją rozszerzyć, na przykład, na modele przeżycia, w których sprawdzamy, czy to osoba o najwyższym (względnym) zagrożeniu umiera najwcześniej). Ale w praktyce nigdy by tego nie używał. To rzadki przypadek, gdy ktoś wie ma jedną zdrową i jedną osobę, nie wie, która osoba jest chora i musi zdecyduj, które z nich leczyć. (W każdym razie decyzja jest łatwa; potraktuj tę z najwyższym szacowanym ryzykiem.)
Myślę więc, że zbadanie rzeczywistej krzywej ROC będzie bardziej przydatne niż samo patrzenie sumaryczny pomiar AUC. A jeśli użyjesz ROC razem z (szacunkami) kosztów fałszywie pozytywnych i fałszywie negatywnych, wraz ze stawkami bazowymi tego, czego się uczysz, możesz gdzieś dostać.
> Należy również pamiętać, że AUC mierzy tylko dyskryminację , a nie kalibrację. Oznacza to, że mierzy, czy można rozróżnić dwie osoby (jedną chorą i jedną zdrową) na podstawie oceny ryzyka. W tym celu bierze pod uwagę jedynie względne wartości ryzyka (lub rangi, jeśli wolisz, por. Interpretacja testu Wilcoxona – Manna – Whitneya), a nie wartości bezwzględne, które powinieneś zainteresuj się. Na przykład, jeśli podzielisz każde oszacowanie ryzyka z modelu logistycznego przez 2, otrzymasz dokładnie takie same AUC (i ROC).
Podczas oceny modelu ryzyka kalibracja jest również bardzo ważna. Aby to zbadać, przyjrzysz się wszystkim pacjentom z wynikiem ryzyka około, np. 0,7 i zobaczysz, czy około 70% z nich faktycznie było chorych. Zrób to dla każdego możliwego wyniku ryzyka (prawdopodobnie używając jakiegoś rodzaju wygładzania / lokalnej regresji). Wykreśl wyniki, a otrzymasz graficzną miarę kalibracji .
Jeśli masz model z zarówno dobrą kalibracją, jak i dobrą dyskryminacją, to zacznij mieć dobry model. 🙂
Komentarze
Odpowiedź
Spójrz na to pytanie: Zrozumienie krzywej ROC
Oto jak zbudować krzywą ROC (na podstawie tego pytania):
Rysowanie krzywej ROC
biorąc pod uwagę zestaw danych przetwarzanych przez klasyfikator rankingowy
- przykłady testów rangowych przy malejącym wyniku
- zacznij od $ (0, 0) $
- dla każdego przykładu $ x $ (w malejąco)
- jeśli $ x $ jest dodatnie, przesuń $ 1 / \ text {pos} $ w górę
- jeśli $ x $ jest ujemne, przesuń $ 1 / \ text {neg} $ w prawo
gdzie $ \ text {pos} $ i $ \ text {neg} $ to odpowiednio ułamki przykładów pozytywnych i negatywnych.
Możesz użyć tego pomysłu do ręcznego obliczenia AUC ROC za pomocą następującego algorytmu:
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
Ten ładny obrazek animowany gif powinien to zilustrować jaśniejszy proces
Komentarze
- Dzięki @Alexey Grigorev, to świetna grafika i prawdopodobnie okaże się przydatna w przyszłości! +1
- Proszę wyjaśnić trochę o ” ułamkach pozytywnych i negatywnych przykładów „, czy masz na myśli najmniejsza wartość jednostkowa dwóch osi?
- @Allan Ruin:
pos
oznacza tutaj liczbę dodatnich danych. Powiedzmy, że masz 20 punktów danych, w których 11 punktów to 1. Tak więc podczas rysowania wykresu mamy prostokąt 11×9 (wysokość x szerokość). Alexey Grigorev skalował, ale po prostu pozwól mu ' s, jeśli chcesz. Teraz po prostu przesuń 1 na wykresie w każdym kroku.
Odpowiedź
Post Karla zawiera dużo doskonałych informacji, ale w ciągu ostatnich 20 lat nie widziałem jeszcze przykładu krzywej ROC, która zmieniłaby czyjekolwiek myślenie w dobrym kierunku. Moim skromnym zdaniem jedyną wartością krzywej ROC jest to, że jej pole jest równe bardzo użytecznemu prawdopodobieństwu zgodności. Sama krzywa ROC kusi czytelnika, aby użył wartości odcięcia, co jest złą praktyką statystyczną.
Jeśli chodzi o ręczne obliczanie $ c $ -index, utwórz wykres z $ Y = 0,1 $ na $ x $ -osi i ciągły predyktor lub przewidywane prawdopodobieństwo, że $ Y = 1 $ na osi $ y $. Jeśli połączysz każdy punkt z $ Y = 0 $ z każdym punktem z $ Y = 1 $, proporcja prostych o dodatnim nachyleniu będzie równa prawdopodobieństwu zgodności.
Wszelkie miary, które mają mianownik $ n $ w tym ustawieniu to niewłaściwe zasady oceny dokładności i należy ich unikać. Obejmuje to odpowiednio sklasyfikowane proporcje, czułość i swoistość.
Dla funkcji R Hmisc
pakiet rcorr.cens
wydrukuj cały wynik, aby zobaczyć więcej informacji, zwłaszcza błąd standardowy.
Komentarze
- Dziękuję, @Frank Harell, doceniam twoją perspektywę. Po prostu używam statystyki c jako prawdopodobieństwa zgodności, ponieważ ' nie lubię wartości odcięcia. Jeszcze raz dziękujemy!
Odpowiedź
Oto alternatywa dla naturalnego sposobu obliczania AUC przy użyciu prostego trapezu reguła, aby uzyskać obszar pod krzywą ROC.
AUC jest równe prawdopodobieństwu, że losowo próbkowana pozytywna obserwacja ma przewidywane prawdopodobieństwo (bycia dodatnim) większe niż losowo próbkowana negatywna obserwacja. Możesz użyć tego do dość łatwego obliczenia AUC w dowolnym języku programowania, przechodząc przez wszystkie parami kombinacje pozytywnych i negatywnych obserwacji. Możesz również losowo próbkować obserwacje, jeśli wielkość próby była zbyt duża. Jeśli chcesz obliczyć AUC za pomocą długopisu i papieru, może to nie być najlepsze podejście, chyba że masz bardzo małą próbkę / dużo czasu. Na przykład w 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
Możemy zweryfikować za pomocą pakietu pROC
:
library(pROC) auc(y, probs) Area under the curve: 0.6287
Korzystanie z próbkowania losowego:
mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE)) [1] 0.62896
Odpowiedź
- Masz prawdziwą wartość obserwacji.
- Oblicz prawdopodobieństwo późniejsze, a następnie uszereguj obserwacje według tego prawdopodobieństwa.
- Zakładając prawdopodobieństwo odcięcia $ P $ i liczbę obserwacji $ N $:
$$ \ frac {\ text {Suma prawdziwych rang} -0,5PN (PN + 1)} { PN (N-PN)} $$
Komentarze
- @ user73455 … 1) Tak, mam prawdziwą wartość do obserwacji. 2) Czy prawdopodobieństwo późniejsze jest synonimem prawdopodobieństwa przewidywanego dla każdej z obserwacji? 3) Zrozumiano; jednak czym jest ” Suma prawdziwych rang ” i jak obliczyć tę wartość? Może przykład pomoże ci dokładniej wyjaśnić tę odpowiedź? Dziękuję!
sens=c(sens,0); omspec=c(omspec,0)
, nie powinno ' t to jestsens=c(0, sens); omspec=c(0, omspec)
? Drukuje poprawnie z początkowym0
, ale nie tak, jak jest obecnie w odpowiedzi.sens
, jak iomspec
przed wydrukowaniem).