Interesuje mnie ręczne obliczanie pola pod krzywą (AUC) lub statystyki c dla binarnego modelu regresji logistycznej.

Na przykład w zbiorze danych sprawdzania poprawności mam prawdziwą wartość zmiennej zależnej retencji (1 = zachowane; 0 = nie zachowane), a także przewidywany status retencji dla każdej obserwacji wygenerowanej przez moją analizę regresji przy użyciu modelu, który był zbudowany przy użyciu zestawu treningowego (będzie wynosił od 0 do 1).

Moje początkowe myśli polegały na zidentyfikowaniu „poprawnej” liczby klasyfikacji modeli i po prostu podzieleniu liczby „poprawnych” obserwacji przez liczbę wszystkich obserwacji w celu obliczenia statystyki c. Przez „poprawne”, jeśli prawdziwy status retencji obserwacji = 1, a przewidywany status retencji wynosi> 0,5, to jest to „poprawna” klasyfikacja. Dodatkowo, jeśli prawdziwy status retencji obserwacji = 0, a przewidywany status retencji wynosi <0,5, to jest to również „poprawna” klasyfikacja. Zakładam, że „remis” wystąpiłby, gdy przewidywana wartość = 0,5, ale zjawisko to nie występuje w moim zbiorze danych walidacyjnych. Z drugiej strony „niepoprawne” klasyfikacje miałyby miejsce, gdyby prawdziwy status retencji obserwacji = 1, a przewidywany status retencji wynosi <0. 5 lub jeśli prawdziwy status retencji dla wyniku = 0, a przewidywany status retencji wynosi> 0,5. Jestem świadomy TP, FP, FN, TN, ale nie wiem, jak obliczyć statystykę c, biorąc pod uwagę te informacje.

Poleciłbym artykuł Hanleya i McNeila z 1982 r. „ Znaczenie i wykorzystanie obszaru pod krzywą charakterystyki odbiornika (ROC) ”.

Przykład

Posiadają poniższą tabelę statusu choroby i wyniku testu (odpowiadającą na przykład oszacowanemu ryzyku z modelu logistycznego). Pierwsza liczba po prawej to liczba pacjentów z prawdziwym statusem choroby „normalnym”, a druga liczba to liczba pacjentów z prawdziwym statusem choroby „nienormalnym”:

(1) Zdecydowanie normalny: 33/3
(2) Prawdopodobnie normalny: 6/2
(3) Wątpliwy: 6/2
(4) Prawdopodobnie nienormalny: 11/11
(5) Zdecydowanie nienormalny: 2/33

Tak więc 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 gdy wynosi 5, „Zdecydowanie nieprawidłowy” pacjenci są zwykle nienormalni (prawda dla 33 z 35 pacjentów), więc predyktor ma sens. Ale jak powinniśmy oceniać pacjenta z wynikiem 2, 3 lub 4? To, co ustaliliśmy jako granicę dla oceny pacjentów jako nienormalne lub normalne, określa determinację czułości i swoistości wynikowego testu.

Czułość i swoistość

Możemy obliczyć szacunkową czułość i swoistość dla różnych wartości odcięcia. (Po prostu napiszę teraz „czułość” i „specyficzność”, pozwalając na domniemany charakter wartości.)

Jeśli wybierzemy naszą wartość graniczną, abyśmy sklasyfikowali wszystkich pacjentów jako nienormalnych, bez względu na to, co mówią ich wyniki testu (tj. Wybieramy wartość graniczną 1+), uzyskamy czułość 51/51 = 1. Swoistość wyniesie 0 / 58 = 0. Nie brzmi tak dobrze.

OK, więc wybierzmy mniej ścisłą granicę. Klasyfikujemy pacjentów jako nienormalnych tylko wtedy, gdy mają wynik testu 2 lub wyższy. Następnie tęsknimy za 3 nienormalnymi pacjentami i mamy czułość 48/51 = 0,94. Ale mamy znacznie zwiększoną specyficzność, wynoszącą 33/58 = 0,57.

Możemy teraz to kontynuować, wybierając różne wartości odcięcia (3, 4, 5,> 5). (W ostatnim przypadku nie będziemy klasyfikować żadnych pacjentów jako nienormalnych, nawet jeśli mają najwyższy możliwy wynik testu wynoszący 5).

Krzywa ROC

Jeśli zrobimy to dla wszystkich możliwych wartości odcięcia, a wykreślymy czułość względem 1 minus swoistości, 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)) )

Dane wyjściowe 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

Za pomocą 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ć powierzchnię pod krzywą ROC, stosując wzór na 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, u których jedna jest normalna, a druga jest nienormalna, możemy obliczyć, jak często jest to ta nienormalna, która ma najwyższy (najbardziej „nietypowy”) wynik testu (jeśli mają tę samą wartość, policz to jako „połowę zwycięstwa”):

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

Odpowiedź to ponownie 0,8931711, obszar pod krzywą ROC. Tak będzie zawsze.

Graficzny widok zgodności

Jak zauważył Harrell w swojej odpowiedzi, ma to również interpretację graficzną. Narysujmy wynik testu (oszacowanie ryzyka) na osi y i prawdziwy stan choroby na osi x (tutaj z pewnym drżeniem, 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 („nienormalny” pacjent). Proporcja linii o dodatnim nachyleniu (tj. Proporcja par zgodnych ) jest wskaźnikiem 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ę powiązań (wynik równego ryzyka), ale przy pewnym drżeniu 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 wskaźnik zgodności będzie wysoki. Widzimy także wkład do indeksu z każdego rodzaju pary obserwacji. Większość pochodzi od normalnych pacjentów z wynikiem ryzyka 1 w połączeniu z nietypowymi pacjentami z wynikiem ryzyka 5 (1–5 par), ale całkiem sporo pochodzi również od 1–4 par i 4–5 par. I bardzo łatwo jest obliczyć rzeczywisty wskaźnik zgodności na podstawie definicji nachylenia:

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

Odpowiedź to ponownie 0,8931711, tj. AUC.

Test Wilcoxona – Manna – Whitneya

Istnieje ścisły związek między miarą zgodności a testem Wilcoxona – Manna – Whitneya. W rzeczywistości to drugie sprawdza, czy prawdopodobieństwo zgodności (tj. Czy jest to nieprawidłowy pacjent w losowej parze normalna-nienormalna, która będzie miała najbardziej „nienormalnie” wynik testu) wynosi dokładnie 0,5. A jego statystyka testowa jest po prostu transformacją 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 znaną liczbę:

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

Tak, to 0,8931711, obszar pod krzywą ROC.

Łatwiejsze sposoby obliczania AUC (w R)

Ale ułatwmy sobie życie. Istnieją różne pakiety, które automatycznie obliczają dla nas AUC.

Pakiet Epi

EpiPakiet tworzy piękny krzywej ROC z różnych statystyk (w tym AUC) osadzone:

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

Pakiet pROC

Podoba mi się również ten pROCpakiet, ponieważ może wygładzić oszacowanie ROC (i obliczyć oszacowanie AUC na podstawie wygładzonego ROC):

(Czerwona linia jest oryginalnym ROC, a czarna linia jest wygładzonym ROC. Zwróć również uwagę na domyślny współczynnik kształtu 1: 1. Sensowne jest użycie tego, ponieważ zarówno czułość, jak i swoistość ma zakres 0–1.)

Szacowana wartość AUC dla wygładzonego ROC wynosi 0,9107, podobnie jak, ale nieco większa, niż wartość AUC dla nie wygładzonego ROC (jeśli spojrzysz na rysunek, możesz łatwo zobaczyć, dlaczego jest większy). (Chociaż naprawdę mamy zbyt mało możliwych wyraźnych wartości wyników testu, aby obliczyć płynne AUC).

Pakiet rms

rmsPakiet Harrella umożliwia obliczanie różnych powiązanych statystyk zgodności za pomocą rcorr.cens()funkcji. C IndexW jego wyjściu jest AUC:

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

Pakiet caTools

Wreszcie mamy caToolspakiet i jego colAUC()funkcję. Ma kilka zalet w porównaniu z innymi pakietami (głównie szybkość i możliwość pracy z danymi wielowymiarowymi - patrz ?colAUC), które czasem mogą być pomocne. Ale oczywiście daje tę samą odpowiedź, którą obliczaliśmy w kółko:

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

Ostatnie słowa

Wiele osób uważa, że ​​AUC mówi nam, jak „dobry” jest test. Niektórzy uważają, że AUC to prawdopodobieństwo, że test prawidłowo sklasyfikuje pacjenta. To nie . Jak widać z powyższego przykładu i obliczeń, AUC mówi nam coś o rodzinie testów, po jednym teście dla każdej możliwej wartości granicznej.

AUC oblicza się na podstawie wartości granicznych, których nigdy nie zastosowalibyśmy w praktyce. Dlaczego powinniśmy dbać o wrażliwość i swoistość „nonsensownych” wartości odcięcia? Nadal na tym opiera się (częściowo) AUC. (Oczywiście, jeśli AUC jest bardzo bliskie 1, prawie każdy możliwy test będzie miał wielką moc dyskryminacyjną i wszyscy bylibyśmy bardzo szczęśliwi.)

Interpretacja AUC „losowa normalna – nienormalna” jest dobra (i może być rozszerzona, na przykład na modele przeżycia, gdzie widzimy, czy to osoba o najwyższym (względnym) ryzyku, która umrze najwcześniej). Ale nigdy nie użyłby tego w praktyce. To rzadki przypadek, gdy wiadomo, że ma się jedną osobę zdrową i jedną chorą, nie wie, która osoba jest chora i musi zdecydować, który z nich leczyć. (W każdym razie decyzja jest łatwa; traktuj tę o najwyższym szacowanym ryzyku).

Myślę więc, że badanie rzeczywistej krzywej ROC będzie bardziej przydatne niż tylko spojrzenie na miarę podsumowującą AUC. A jeśli użyjesz ROC razem z (szacunkowymi) kosztami fałszywie dodatnich i fałszywych negatywów, wraz z podstawowymi stawkami tego, co studiujesz, możesz gdzieś dostać.

Należy również pamiętać, że AUC mierzy tylko dyskryminację , a nie kalibrację. Oznacza to, że mierzy, czy możesz rozróżnić dwie osoby (jedną chorą i jedną zdrową), na podstawie oceny ryzyka. W tym celu uwzględnia jedynie względne wartości ryzyka (lub rangi, jeśli chcesz, por. Interpretację testu Wilcoxona – Manna – Whitneya), a nie wartości bezwzględne, którymi powinieneś być zainteresowany. Na przykład, jeśli podzielisz każde ryzyko oszacuj na podstawie modelu logistycznego o 2, otrzymasz dokładnie takie same AUC (i ROC).

Przy ocenie modelu ryzyka bardzo ważna jest również kalibracja . Aby to zbadać, przyjrzysz się wszystkim pacjentom z wynikiem ryzyka wynoszącym około, np. 0,7, i zobaczysz, czy około 70% z nich faktycznie chorowało. Zrób to dla każdego możliwego wyniku ryzyka (ewentualnie używając pewnego rodzaju wygładzania / regresji lokalnej). Wykreśl wyniki, a otrzymasz graficzną miarę kalibracji .

Jeśli masz model z obu dobrej kalibracji i dobrej dyskryminacji, potem zaczynasz mieć dobry model. :)

Spójrz na to pytanie: Zrozumienie krzywej ROC

Oto jak zbudować krzywą ROC (z tego pytania):

Rysowanie krzywej ROC

biorąc pod uwagę zestaw danych przetworzony przez klasyfikator rankingu

• przykłady testów rangowych dotyczących malejącego wyniku
• $(0, 0)$
• $x$
  • $x$$1/\text{pos}$
  • $x$$1/\text{neg}$

$\text{pos}$$\text{neg}$

Możesz użyć tego pomysłu do ręcznego obliczania AUC ROC przy użyciu 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 animowany obraz gif powinien lepiej zilustrować ten proces

Post Karla zawiera wiele doskonałych informacji. Ale nie widziałem jeszcze w ciągu ostatnich 20 lat przykładu krzywej ROC, która zmieniłaby czyjeś myślenie w dobrym kierunku. Jedyną wartością krzywej ROC jest, moim skromnym zdaniem, to, że jej powierzchnia jest równa bardzo przydatnemu prawdopodobieństwu zgodności. Sama krzywa ROC kusi czytelnika do stosowania wartości odcięcia, co jest złą praktyką statystyczną.

$c$$Y=0,1$$x$$Y=1$$y$$Y=0$$Y=1$

$n$

W przypadku funkcji Hmiscpakietu R rcorr.censwydrukuj cały wynik, aby wyświetlić więcej informacji, zwłaszcza błąd standardowy.

Oto alternatywa dla naturalnego sposobu obliczania AUC po prostu za pomocą reguły trapezoidalnej, aby uzyskać pole pod krzywą ROC.

AUC jest równe prawdopodobieństwu, że losowo dobrana obserwacja dodatnia ma przewidywane prawdopodobieństwo (bycia dodatnim) większa niż losowo dobrana obserwacja ujemna. Możesz tego użyć do dość łatwego obliczenia AUC w dowolnym języku programowania, przechodząc przez wszystkie kombinacje pozytywnych i negatywnych obserwacji. Można również losowo próbkować obserwacje, jeśli wielkość próby jest zbyt duża. Jeśli chcesz obliczyć AUC za pomocą pióra 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ą pROCpakietu:

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

Korzystanie z losowego próbkowania:

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

1. Masz prawdziwą wartość dla obserwacji.
2. Oblicz prawdopodobieństwo tylne, a następnie uszereguj obserwacje według tego prawdopodobieństwa.
3. $P$$N$
   $$\frac{\text{Sum of true ranks}-0.5PN(PN+1)}{PN(N-PN)}$$