Tag: analiza przeżycia

Przecinające się krzywe przeżycia

Spotkałem się ostatnio z ciekawym problemem.
Mamy dwie grupy pacjentów na dwóch różnych schematach leczenia i chcemy porównać ich dalsze losy, a konkretnie krzywe niepowodzenia leczenia (prawdopodobieństwo zgonu/wznowy).
Na pierwszy rzut oka klasyczna analiza przeżycia, test log-rank i po sprawie.

Ale szybko okazuje się, że krzywe przeżycia się przecinają, co więcej oczekiwać można tego po wcześniejszej rozmowie z lekarzem. Jeden schemat leczenia jest bardziej agresywny więc może prowadzić do gorszych rokowań krótkookresowych, ale lepszych w dalszej perspektywie.

Klasyczny test dla krzywych przeżycia oparty jest o odległość pomiędzy krzywymi, mierzoną jest jako ważona suma kwadratów odległości w poszczególnych punktach. Ale gdy krzywe się przecinają to taki test ma niską moc i nie ma sensu go stosować.

A więc co robić?
Ciekawe studium symulacyjne porównujące różne podejścia do testowania przecinających się krzywych zostało opublikowane dwa lata temu w Plos One (!).
Okazuje się, że dobrze sprawdza się rodzina testów Renyi, która jest oparta o supremum ważonych odległości pomiędzy krzywymi przeżycia.
W R te testy są zaimplementowane w pakiecie survMisc w funkcji comp. Jest to znacznie mocniejszy test dla przecinających się krzywych.

A przy okazji, okazuje się, że zmianę w hazardach w rozpatrywanym problemie dobrze ilustrują reszty Schonefelda. Poniższy wykres pokazuje, że hazard w jednej grupie jest znacznie wyższy do 12 miesiąca, a później gorsze losy czekają pacjentów drugiej grupy.

Oba wykresy wykonane pakietem survminer.

Opisy osi usunąłem ponieważ wyniki tych analiz jeszcze nie są opublikowane, ale też nazwy nie mają większego znaczenia.

Który telefon się wcześniej zepsuje?

Uczestniczyłem ostatnio w ciekawym (od strony metodologicznej) projekcie (patrz też komentarz na dole). A ponieważ przy okazji używałem kilku ciekawych narzędzi, więc poniżej krótko o nich opowiem.

Problem:

Przypuśćmy, że mamy taki problem. Mamy telefony dwóch marek, np A i B. Załóżmy też, że w tych telefonach może zepsuć się bateria lub ekran, innych uszkodzeń nie rozważamy. Jeżeli się cokolwiek zepsuje, telefon jest wymieniany na nowy, nie ma sensu go naprawiać.
Chcemy teraz sprawdzić czy tempo psucia się ekranu w telefonach marki A jest inne niż w telefonach marki B. Podobnie chcemy porównać działanie baterii.
Przypuśćmy, że w badaniu obserwujemy dużo (tysiące) telefonów każdej z marek. Badanie prowadzimy przez 3 lata, ale niektóre telefony uczestniczą w badaniu od roku, inne od 2 lub 3 lat.
Pytanie, jaką metoda porównać te marki?

Podejście brutalne 1:

Najbardziej brutalne podejście to policzenie jaki procent ekranów zepsuł się w marce A a jaki w marce B i porównanie tych procentów np. dokładnym testem Fishera (w R funkcja fisher.test {stats}, więcej o teście tutaj) lub regresją logistyczną.
Ale…
To nie jest dobre podejście, ponieważ telefony nie były obserwowane przez tak samo długi czas. Mogło się zdarzyć, że w jednej grupie pozyskano więcej telefonów na początku, 3 lata temu, więc ich czas obserwacji był dłuższy. A może w innej grupie pozyskano telefony niedawno, np. rok temu, czas obserwacji jest krótszy i dlatego mniej ich się zepsuło.

Podejście brutalne 2:

Skoro czas obserwacji jest różny, to porównajmy liczbę dni do czasu zepsucia się ekranu lub baterii. Możemy to zrobić np. testem Wilcoxona (w R funkcja wilcox.test {stats} więcej o teście np. tutaj).
Ale…
To nie jest dobre podejście, ponieważ nie wszystkie telefony się zepsuły, nie dla wszystkich znam czas do zepsucia. Być może w jednej grupie średni czas do uszkodzenia jest niższy, ale też mniej telefonów się zepsuło. Nie chcemy też odrzucać informacji o tych wszystkich telefonach, które się nie zepsuły.

Mniej brutalne podejście 3:

Skoro mamy różne czasy obserwacji i nie dla wszystkich obserwacji wiemy jaki był czas uszkodzenia ekranu lub baterii, to zastosujmy metody analizy przeżycia. Przedstawmy czas działania ekranu/baterii za pomocą krzywych przeżycia (np. Kaplan-Meier) i sprawdźmy czy marki różnią się używając np. testu log rank (w R funkcja survdiff {survival}, więcej o teście).
Ale…
Obserwujemy dwa zdarzenia, uszkodzenie ekranu lub uszkodzenie baterii. Jeżeli wystąpi jedno uszkodzenie to nie widzimy drugiego (telefon oddajemy i nie wiemy kiedy nastąpiłaby druga awaria). Ale analizując uszkodzenia baterii, nie możemy traktować uszkodzeń ekranu jako obserwacje cenzorowane ponieważ nie musi być to cenzorowanie nieinformatywne (brak spełnionych założeń). Być może uszkodzony ekran wpływa na sposób pracy baterii, a być może wadliwe ekrany i baterie produkowane są w jednej fabryce, itp.

Podejście najmniej brutalne 4:

Skoro różne telefony obserwujemy przez różny czas, nie dla wszystkich znamy całkowity czas działania, czyli mamy cenzorowanie, ale obserwujemy dwa (lub więcej) elementów, które mogą ulec uszkodzeniu, to możemy modelować je łącznie używając modelu konkurujących ryzyk.

W R można to zrobić używając pakietu mstate (artykuł z JSS) czy też od jakiegoś czasu można to zrobić pakietem survival (teraz funkcja Surv{survival} może przyjąć argument type=mstate).
Do porównania krzywych można wykorzystać test Gray’s log-rank test. A w pakiecie ggplotit jest nawet przeciążona funkcja ggplotit() rysujące skumulowane funkcje zdarzeń dla obiektów klasy survfitms lub Cuminc.

btw:
Projekt na którym pracowałem dotyczył biologii molekularnej i analizy przeżycia pacjentów, ale dla przykładu i wnioskowania to akurat bez znaczenia. Podobna metodologię można zastosować do analizy usterek w samochodach czy kliknięć internautów na stronach www.

Fundacja Nauki Polskiej, analiza przeżycia, wiek habilitacji i apel w sprawie danych

Zdarzało mi się na ramach tego bloga czepiać wykresów, że czegoś nie widać lub że widać coś czego nie ma. Dziś będę czepiał się wypowiedzi, które sugerują że pewien wniosek został wysnuty na podstawie wnioskowania statystycznego, ale gdy się zastanowić nad tym co to za wnioskowanie to pojawia się więcej pytań niż odpowiedzi.

Takie problemy są powszechne w gazetach skierowanych do szerokiego grona odbiorców. W tym przypadku jednak rzecz dotyczy zdania z ramki na 9 stronie rocznego raportu działania Fundacji na rzecz Nauki Polskiej (a więc największej w Polsce pozarządowej organizacji wspierającej Polską Naukę). FNP to organizacja od kórej można wymagać więcej, a mam też nadzieję, że wybaczy mi czepialstwo.

Chodzi o zdanie
,,Z badania karier laureatów programu START przeznaczonego dla najmłodszych uczonych, który fundacja realizuje od 1993r., wynika, że uzyskuja oni habilitację średnio o 9-10 lat wcześniej niż osoby niekorzystające z tego programu stypendialnego”.

Ok, co jest nie tak z tym zdaniem?

1. brak odnośnika do danych lub raportu na którym się opierano. Od innych organizacji bym tego nie oczekiwał ale FNP powinno wytyczać kierunki, a tym samym mogłoby udostępniać dane na podstawie których wnioskują.

2. Nie jest jasne jaką relację to zdanie ma pokazać. Czy to, że ci młodzi uczeni szybciej zrobili habilitację dzięki stypendium FNP, czy też czy fundacja umiejętnie odnajduje osoby które szybko zrobią habilitację (zgodnie z dewizą fundacji ,,wspierać najlepszych, aby mogli stać się jeszcze lepsi”)

3. Nie jest jasne kim są ci ,,niekorzystający z tego programu”.
Czy to rówieśnicy osób korzystajacych z programu, o podobnym potencjalne. Taka grupa kontrolna? Raczej nie. Bezsensowne byłoby losowe nieprzydzielanie stypendiów tylko po do by zbadać efekt programu.
Czy to rówieśnicy osób korzystających z programu, którzy nie aplikowali lub aplikowli ale nie otrzymali stypendium.
Czy tez wszyscy naukowcy bez wzgledu na wiek. To ostatnie rozwiązanie byłoby niedobre. Kiedys habilitacje i doktoraty robilo sie dłużej, nie bylo cztero czy trzyletnich studiow doktoranckich ale doktoraty i habilitacje robili asystenci/adiunkci laczac prace naukowa z innymi obowiazkami bez presji ze po ośmiu latach zatrudnienia będą wyrzuceni jeżeli habilitacji nie zrobią.

4. Najlepszą grupa kontrolną byliby rówieśnicy, ale wiele wskazuje że tak nie było.
Jeżeli program realizowany jest od 1993 roku a raport dotyczy roku 2010 to najstarsi stypendyści są 17 lat po otrzymaniu stypendium. Przyjmijmy uproszczenie, ze co roku podobna liczba osob otrzymuje stypendium START, wiec stypendyści są średnio 8.5 roku po otrzymaniu sypendium.
Stypendium START jest dla osob mlodych (do 30 roku zycia) najczesciej swieżo po doktoracie.
Trudno odgadnąć w ile lat robi sie szybko habilitaję, ale mysle ze srednio 6 lat w grupie stypendystow to bylby dobry wynik.
W grupie niestypendystow musialoby to być wiec 15-16 lat lub więcej po doktoracie. Ale program stypendialny nie jest tak długo prowadzony by mieć rówieśników robiących habilitację w takim wieku.

Ok, czyli wybór grupy kontrolnej jest niejasny, teraz zastanówmy się co z metodologią.
Do porównania obu grup uzyto średnich liczby lat do habilitacji.
Ale liczba lat do habilitacji to zmienna cenzurowana. W badanej grupie z pewnoscia sa osoby ktore jeszcze nie zrobiły habilitacji i pracuja w nauce oraz osoby ktore zrezygnowały ze ścieżki naukowej i nie beda robily habilitacji.
Liczenie średniej z tylko tych osób które zrobiły habilitacje jest błędem, poniewaz gubi informację jaka frakcja osób zrobiła habilitacje. Nawet pomijając te problemy to dla wielu rozkładów średnia nie jest dobrym miernikiem czegokolwiek.

Ok, sposób porównywani grup pozostawia wiele do zyczenia, ale takich porównań będzie coraz więcej, więc warto się zastanowić jak je robić. Np. czy czas do habilitacji różni sie i jak pomiędzy róznymi jednostkami naukowymi.

Odpowiednie byłyby narzędzia z analizy przyżycia, np. krzywa Kaplana Meiera pokazujące jaka frakcja osób zrobila habilitacje do k-tego roku. Lub funkcja intensywnosci / funkcja hazardu pokazujaca jaka jest częstość robienia habilitacji w k-tym roku.
Z krzywych Kaplana Meiera mozna by zobaczyć w jakim wieku najczęściej robiona jest habilitacja.

Ok, ale aby zrobić taka analizą trzeba mieć dane. Najlepiej w postaci wektora danych dla osob z informacja w jakim wieku dana osoba zrobila doktorat / habilitacje / czy otrzymala stypendium FNP i np gdzie teraz pracuje.

Ale skad takie dane wziac?

Tu prosba do Was, ktokolwiek widział, ktokolwiek wie o miejscu z ktorego mozna takie dane dostac prosze o informacje.

Próbowałem serwisu http://nauka-polska.pl, niestety jest tam limit na liczbe zapytań do bazy danych które mozna wykonac z jednego adresu IP w jednostce czasu. Ten limit można by obejsc ale sam fakt ze go ustawiono oznacza, ze twórcom nie spodobaloby sie twórcom gdyby ich dane ściągnąć i upublicznić. Ale może są jakies inne źródła publicznie dostepnych danych?