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