Jak wykonać analizę ROC w R za pomocą modelu Coxa

10

Stworzyłem kilka modeli regresji Coxa i chciałbym zobaczyć, jak dobrze działają te modele, i pomyślałem, że być może krzywa ROC lub statystyka c mogą być przydatne podobnie jak w przypadku tych artykułów:

JN Armitage och JH van der Meulen, „Identyfikacja chorób współistniejących u pacjentów chirurgicznych przy użyciu danych administracyjnych z Royal College of Surgeons Charlson Score”, British Journal of Surgery, vol. 97, num. 5, ss. 772–781, maj 2010 r.

Armitage zastosował regresję logistyczną, ale zastanawiam się, czy można użyć modelu z pakietu przetrwania, survivalROC daje wskazówkę, że jest to możliwe, ale nie mogę wymyślić, jak to zrobić, aby działało z regresją Coxa.

Byłbym wdzięczny, gdyby ktoś pokazał mi, jak przeprowadzić analizę ROC na tym przykładzie:

library(survival)
data(veteran)

attach(veteran)
surv <- Surv(time, status)
fit <- coxph(surv ~ trt + age + prior, data=veteran)
summary(fit)

Jeśli to możliwe, doceniłbym zarówno wyjściową statystykę c, jak i ładny wykres

Dzięki!

Aktualizacja

Dziękuję bardzo za odpowiedzi. @Dwin: Chciałbym mieć pewność, że zrozumiałem to przed wybraniem odpowiedzi.

Obliczenia, jakie rozumiem, zgodnie z sugestią DWina:

library(survival)
library(rms)
data(veteran)

fit.cph <- cph(surv ~ trt + age + prior, data=veteran, x=TRUE, y=TRUE, surv=TRUE)

# Summary fails!?
#summary(fit.cph)

# Get the Dxy
v <- validate(fit.cph, dxy=TRUE, B=100)
# Is this the correct value?
Dxy = v[rownames(v)=="Dxy", colnames(v)=="index.corrected"]

# The c-statistic according to the Dxy=2(c-0.5)
Dxy/2+0.5

Nie jestem zaznajomiony z funkcją sprawdzania poprawności i ładowaniem, ale po przyjrzeniu się prof. Odpowiedź Franka Harrela tutaj na R-help Pomyślałem, że to prawdopodobnie sposób na zdobycie Dxy. Pomoc w sprawdzaniu poprawności stanów:

... Korelacja rang Somxa do obliczenia przy każdym ponownym próbkowaniu (zajmuje to nieco więcej czasu niż statystyki oparte na prawdopodobieństwie). Wartości odpowiadające wierszowi Dxy są równe 2 * (C - 0,5), gdzie C jest wskaźnikiem C lub prawdopodobieństwem zgodności.

Chyba najbardziej wprawiają mnie w zakłopotanie kolumny. Doszedłem do wniosku, że należy użyć poprawionej wartości, ale tak naprawdę nie zrozumiałem poprawnego wyniku:

      index.orig training    test optimism index.corrected   n
Dxy      -0.0137  -0.0715 -0.0071  -0.0644          0.0507 100
R2        0.0079   0.0278  0.0037   0.0242         -0.0162 100
Slope     1.0000   1.0000  0.2939   0.7061          0.2939 100
...

W pytaniu R-help zrozumiałem, że powinienem mieć „surv = TRUE” w cph, jeśli mam warstwy, ale nie jestem pewien, jaki jest cel parametru „u = 60” w funkcji sprawdzania poprawności. Byłbym wdzięczny, gdybyś pomógł mi je zrozumieć i sprawdzić, czy nie popełniłem żadnych błędów.

Max Gordon
źródło
2
Prawdopodobnie rzuciłbym okiem na pakiet rms i jego cph()komendę.
chl
2
index.correctednależy podkreślić. Są to szacunki prawdopodobnych przyszłych wyników. u=60nie jest potrzebne, validateponieważ nie masz warstw. Jeśli posiadasz warstwy, krzywe przeżycia mogą się krzyżować i musisz określić konkretny punkt czasowy, aby uzyskać uogólniony obszar ROC.
Frank Harrell,

Odpowiedzi:

2

@chl wskazał konkretną odpowiedź na twoje pytanie. Funkcja pakietu „rms” cphwytworzy Somers-D, który można w prosty sposób przekształcić w indeks c. Jednak Harrell (który wprowadził indeks c do praktyki biostatystycznej) uważa, że ​​jest to nierozsądne jako ogólna strategia oceny środków prognostycznych, ponieważ ma niską zdolność do dyskryminacji między alternatywami. Zamiast polegać na literaturze chirurgicznej w zakresie wskazówek metodycznych, mądrzej byłoby poszukać zgromadzonej mądrości w tekście Harrella, „Strategie modelowania regresji” lub „Modelach prognoz klinicznych” Steyerberga.

DWin
źródło
4
DxyC
Dziękuję za odpowiedzi, moja sytuacja jest taka, że ​​mam trzy różne wyniki, które chcę porównać i zobaczyć, jak się zachowują. Nie miałem czasu zajrzeć do części Somers-D i wrócę, kiedy będę miał czas (rzuciłem okiem i nie znalazłem nic przydatnego). Zamówiłem także książkę @FrankHarrell, „Regression Modeling Strategies”, ISBN 13: 978-0387952321 i mam nadzieję, że pomoże mi w dokonaniu wyboru.
Max Gordon
2
Ponieważ Dxy = 2 * (c-0,5), obliczenie c dla Dxy powinno być trywialne.
DWin
3

χ2

Frank Harrell
źródło
+1 za prowadzenie mnie we właściwym kierunku. Właśnie skończyłem robić statystyki C, a bardziej szczegółowy wynik, na który patrzę, miał statystykę C wynoszącą 0,4365081, podczas gdy drugi miał 0,4414625 (chyba powinienem liczyć 0,5-Dxy / 2 w moim przypadku). Sporo czasu poświęciłem na obliczenie mojej 140 000 próbki; Musiałem obniżyć bootstrapy do 10 i nie jestem pewien, jaki to ma wpływ. Z niecierpliwością czekam na przeczytanie twojej książki (jest w mailu) i mam nadzieję, że pomoże mi to lepiej zrozumieć metodologię i porównać statystyki C ze wskaźnikiem adekwatności.
Max Gordon,
Dobrze. Nie jest łatwo stwierdzić, czy 0,44 vs. 0,43 znaczy wiele bez patrzenia na rozkłady przewidywanych wartości.
Frank Harrell,
Rozumiem, że ciężko jest komentować takie liczby. Spróbuję zajrzeć do dystrybucji. Moją główną interpretacją wyniku jest to, że mój model wyjaśnia bardzo mało i chociaż istnieje niewielka różnica, prawdopodobnie nie jest ona bardzo ważna. Byłoby interesujące, czego można się spodziewać w warunkach przeżycia - osiągnięcie wartości 0,8, tak jak to zrobili w analizie, do której nawiązałem w moim pytaniu, wydaje się dość odległe ... ale z drugiej strony moim przetrwaniem jest przeżycie wszczepionej protezy i nie przeżycie pacjenta. Wykorzystali także regresję logistyczną, która może zmienić oszacowanie.
Max Gordon,
Regresja logistyczna nie działałaby, jeśli czas jest ważny lub czas obserwacji różni się w zależności od pacjenta. Wracając do pierwotnego pytania, przewidywane ryzyko będzie miało wąski rozkład, jeśli model wyjaśni bardzo niewielką zmienność.
Frank Harrell,
Właśnie dostałem twoją książkę ... Miałem szybką blokadę w części poświęconej przetrwaniu, ale kiedy wypróbowałem twoje studium przypadku w rozdziale 20, ale dostaję błąd w części przypisującej (w, sz): „zmienna sz nie ma atrybut names () ”. Postępowałem zgodnie z rozdziałem. 8: załadowałem ramkę danych z getHdata (prostata) (nie mogłem znaleźć strony w książce), zrobiłem w <- transcan (~ sz + sg + ap + sbp + dbp + age + wt + hg + ekg + pf + bm + HX, przypisana = T = T przekształcone, imcat = "drzewo", data = prostaty), ale nie mogę znaleźć nic na nazywanie ...
Max Gordon