Analiza przeżycia
Zbigniew Karpinski
16 Sep 2015
Kilka słów tytułem wprowadzenia
Metodę określaną w biologii i epidemiologii mianem „analizy przeżycia” nazywa się w naukach społecznych „analizą historii zdarzeń”. Dzisiejszą notatkę poświęcam szacowaniu modeli historii zdarzeń w R.
Dane wykorzystane w tej notatce zaczerpnąłem z ciekawego artykułu poświęconego związkom między analizą przeżycia a regresją logistyczną. Zanim jednak przejdę do części „właściwej” kilka słów wprowadzenia. Analiza historii zdarzeń to ogólna nazwa całej klasy narzędzi statystycznych zajmujących się opisem i modelowaniem czasu przed wystąpieniem zdarzenia. Typowym przykładem problemu, do którego stosuje się analizę historii zdarzeń, jest czas przeżycia pacjentów po operacji lub terapii. W przykładzie tym mamy więc grupę pacjentów cierpiących na pewną chorobę (np. raka płuc), którzy poddali się leczeniu (np. chemioterapii lub chemioterapii połączonej z zabiegiem chirurgicznym). Interesuje nas czas, jaki upłynie przed śmiercią pacjenta lub remisją choroby. Dla wszystkich pacjentów objętych ryzykiem śmierci lub nawrotu choroby możemy określić stopę ryzyka, \(h(t)\): \[h(t) = \lim_{\Delta t \rightarrow 0} \frac{P\left(t \leqslant T < t + \Delta t \vert T \geqslant t\right)}{\Delta t},\] gdzie \(T\) oznacza rzeczywisty czas wystąpienia zdarzenia (śmierci lub nawrotu choroby w naszym przykładzie), zaś \(t\) to moment, dla którego szacujemy stopę ryzyka. Mówiąc inaczej, mianownik w wyrażeniu powyżej określa prawdopodobieństwo warunkowe zajścia zdarzenia w pewnym niewielkim przedziale czasowym \(\Delta t\), o ile nie wystąpiło ono przed momentem czasowym \(t\).
Inną wielkością, która odgrywa istotną rolę w analizie historii zdarzeń, jest funkcja przeżycia, którą definiuje się w terminach prawdopodobieństwa przeżycia co najmniej do czasu \(t\): \[S(t) = P(T\geqslant t)\] Mówiąc nieco inaczej, dla każdego momentu czasowego \(t\) funkcja przeżycia określa odsetek badanej zbiorowości, dla którego analizowane zdarzenia jeszcze nie zaszło.
Jeśli cały okres objęty analizą podzielimy na mniejsze przedziały (np. dni lub tygodnie lub miesiące), dla każdego przedziału możemy określić trzy następujące wielkości:
- \(r(t)\), czyli liczbę osób objętych ryzykiem wystąpienia zdarzenia na początku przedziału \(t\),
- \(d(t)\), czyli liczbę przypadków wystąpienia zdarzenia w przedziale \(t\),
- \(w(t)\), czyli liczbę przypadków uciętych (censored), tzn. osób, które z różnych powodów wycofały się z badania lub osób, dla których analizowane zdarzenie nie wystąpiło przed zakończeniem obserwacji.
Stopę ryzyka szacuje się za pomocą tych trzech wielkości w następujący sposób: \[\hat{h}(t) = \frac{d(t)}{r(t) - w(t)}\] Z kolei oszacowanie funkcji przeżycia \(S(t)\) ma postać: \[\hat{S}(t) = \left[1 - \hat{h}(1)\right] \cdot \left[1 - \hat{h}(2)\right] \cdots \left[1 - \hat{h}(t - 1)\right]\]
Oczywiście, analiza przeżywalności po zabiegu czy terapii to nie jedyne zastosowania analizy historii zdarzeń. Przykłady użycia analizy historii zdarzeń w naukach społecznych obejmują:
- analizę czasu do znalezienia pierwszej pracy po ukończeniu nauki w grupie absolwentów wyższych uczelni
- analizę czasu przed zawarciem małżeństwa lub wejściem w związek partnerski w grupie osób rozpoczynających samodzielne życie
- analizę czasu przed urodzeniem pierwszego dziecka od momentu rozpoczęcia aktywności seksualnej
- analizę czasu przed powrotem do nałogu w grupie pacjentów, którzy opuszczają ośrodek leczenia uzależnienień,
- analizę czasu przed rezygnacją z dalszego udziału w badaniu w grupie uczestników badania panelowego
- …
O innych zastosowaniach i własnościach analizy historii zdarzeń można poczytać w literaturze źródłowej, np. w bardzo przystępnie napisanej książce Paula D. Allison; autor ten przygotował zresztą wiele różnych opracowań dotyczących analizy historii zdarzeń, możne je wyszukać tutaj. Sporo ciekawych materiałów na temat analizy historii zdarzeń można także znaleźć na stronie brytyjskiej statystyczki Fiony Steele. Materiałów na temat samej metody i jej zastosowań — od prostych wprowadzeń dla początkujących badaczy po zaawansowane opracowania statystyczne dotyczące arkanów szacowania stopy przeżycia — jest zresztą mnóstwo.
Dane
Jak wspomniałem, dane, które wykorzystuję w dzisiejszej notatce, pochodzą z cytowanego artykułu Efrona. Mówiąc dokładnie, Efron wykorzystał w swoim artykule dane dotyczące 96 pacjentów chorych na nowotwory szyi i głowy, których poddano dwóch rodzajom terapii: radioterapii lub radioterapii wspieranej chemioterapią. Badaczy interesowała różnica między tymi dwoma rodzajami terapii pod względem przeżywalności pacjentów. Mówiąc inaczej, zdarzeniem, które analizowano w opisywanym badaniu, była śmierć pacjenta lub nawrót choroby. Dla każdego pacjenta określono czas do wystąpienia tego zdarzenia.
Źródłowe dane dotyczące przeżywalności pacjentów znajdują się na stronach 4 i 5 artykułu Efrona, pod tabelami 1 i 2. Dane te to po prostu zbiór wartości liczbowych określających, dla każdego z pacjentów osobno, liczbę dni przed wystąpieniem zdarzenia. Wartości te pokazano w okienku poniżej:
## [1] "7" "34" "42" "63" "64" "74+" "83" "84"
## [9] "91" "108" "112" "129" "133" "133" "139" "140"
## [17] "140" "146" "149" "154" "157" "160" "160" "165"
## [25] "173" "176" "185+" "218" "225" "241" "248" "273"
## [33] "277" "279+" "297" "319+" "405" "417" "420" "440"
## [41] "523" "523+" "583" "594" "1101" "1116+" "1146" "1226+"
## [49] "1349+" "1412+" "1417" " 37" "84" "92" "94" "110"
## [57] "112" "119" "127" "130" "133" "140" "146" "155"
## [65] "159" "169+" "173" "179" "194" "195" "209" "249"
## [73] "281" "319" "339" "432" "469" "519" "528+" "547+"
## [81] "613+" "633" "725" "759+" "817" "1092+" "1245+" "1331+"
## [89] "1557" "1642+" "1771+" "1776" "1897+" "2023+" "2146+" "2297+"Jak widzimy, poszczególne wartości są otoczone cudzysłowami. Oznacza to, że program traktuje te wartości jako łańcuchy tekstowe, a nie jako wartości liczbowe. Wynika to stąd, że w niektórych przypadkach wartości kończą się znakiem +, który wskazuje na obserwację uciętą. Innymi słowy, elementy wektora pokazanego powyżej zawierają w sobie dwa rodzaje informacji: (a) o liczbie dni przed wystąpieniem zdarzenia i (b) o tym, czy przypadek jest ucięty. Przed przystąpieniem do analizy musimy więc rozbić wektor na dwa wektory „składowe”. Zaczniemy od wyodrębnienia wartości liczbowych informujących o czasie przed wystąpieniem zdarzenia. W kolejnych poleceniach zakładamy, że wektor pokazany w ostatnim okienku zapisaliśmy pod nazwą x.
Aby przekształcić wektor x do postaci liczbowej, wystarczy wyeliminować znak + z poszczególnych elementów tego wektora, albo zastąpić go znakiem „pustym”. Wykorzystamy do tego celu funkcję str_replace() z pakietu stringr przygotowanego z myślą o łatwiejszym przekształcaniu i przetwarzaniu łańcuchów tekstowych.
str_replace(string = x, pattern = "\\+", replacement = "")## [1] "7" "34" "42" "63" "64" "74" "83" "84" "91" "108"
## [11] "112" "129" "133" "133" "139" "140" "140" "146" "149" "154"
## [21] "157" "160" "160" "165" "173" "176" "185" "218" "225" "241"
## [31] "248" "273" "277" "279" "297" "319" "405" "417" "420" "440"
## [41] "523" "523" "583" "594" "1101" "1116" "1146" "1226" "1349" "1412"
## [51] "1417" " 37" "84" "92" "94" "110" "112" "119" "127" "130"
## [61] "133" "140" "146" "155" "159" "169" "173" "179" "194" "195"
## [71] "209" "249" "281" "319" "339" "432" "469" "519" "528" "547"
## [81] "613" "633" "725" "759" "817" "1092" "1245" "1331" "1557" "1642"
## [91] "1771" "1776" "1897" "2023" "2146" "2297"Podwójny ukośnik jest w ostatnim poleceniu potrzebny, ażeby program potraktował znak + jako znak, a nie jako operator arytmetyczny. Jak widzimy, polecenie zwraca nam wektor, którego elementy w dalszym ciągu są traktowane przez program jako tekst, ale nie zawierają już kłopotliwego plusa. Aby uzyskać wartości liczbowe, wstawmy ostatnie polecenie jako argument dla funkcji as.numeric():
as.numeric(str_replace(string = x, pattern = "\\+", replacement = ""))## [1] 7 34 42 63 64 74 83 84 91 108 112 129 133 133
## [15] 139 140 140 146 149 154 157 160 160 165 173 176 185 218
## [29] 225 241 248 273 277 279 297 319 405 417 420 440 523 523
## [43] 583 594 1101 1116 1146 1226 1349 1412 1417 37 84 92 94 110
## [57] 112 119 127 130 133 140 146 155 159 169 173 179 194 195
## [71] 209 249 281 319 339 432 469 519 528 547 613 633 725 759
## [85] 817 1092 1245 1331 1557 1642 1771 1776 1897 2023 2146 2297Uzyskaliśmy zatem zbiór wartości liczbowych mierzących czas przed wystąpieniem zdarzenia (lub przed ostatnią obserwacją w przypadku obserwacji uciętych). Aby uzyskać zmienną identyfikującą obserwacje ucięte, skorzystamy z funkcji str_detect() z pakietu stringr. Funkcja ta sprawdza, czy zadany znak występuje w danym ciągu znaków i zwraca wartość logiczną TRUE, jeśli tak. W naszym przykładzie, jeśli element wektora x zawiera znak +, element ten wskazuje na obserwację uciętą:
str_detect(string = x, pattern = "\\+")## [1] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [23] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [34] TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [45] FALSE TRUE FALSE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE
## [56] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [78] FALSE TRUE TRUE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE
## [89] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUEOstatnia rzecz. Jak pamiętamy, dane w naszym przykładzie dotyczą pacjentów poddanych dwóm rodzajom terapii: (a) radioterapii (nazwijmy ją „terapią A”) oraz (b) radioterapii wspieranej chemioterapią („terapia B”). Pierwszych 51 przypadku objęto terapią A, pozostałe 45 — terapią B. Zmienną określającą rodzaj terapii tworzymy następująco:
rep(LETTERS[1:2], times = c(51, 45))## [1] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"
## [18] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"
## [35] "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A" "A"
## [52] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B"
## [69] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B"
## [86] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B" "B"Wszystkie trzy utworzone przez nas zmienne połączymy teraz w jednej ramce danych, którą nazwiemy df:
df <- data.frame(
months = as.numeric(str_replace(string = x, pattern = "\\+", replacement = ""))/30.348,
status = !str_detect(string = x, pattern = "\\+"),
treatment = rep(LETTERS[1:2], times = c(51, 45))
)Kilka słów tytułem wyjaśnienia. Informację o czasie do wystąpienia zdarzenia zapisywano w badaniu opisywanym przez Efrona w dniach. Efron podzielił liczbę dni przez 30,3, czyli średnią długość miesiąca, aby czas mierzony w dniach zamienić na czas mierzony w miesiącach. Zmienna status określa, czy analizowane zdarzenie (a więc śmierć lub nawrót choroby) wystąpiło czy nie. Zmienną tę uzyskaliśmy przez negację logiczną utworzonej wcześniej zmiennej określającej przypadki ucięte. Mówiąc inaczej, jeśli dana jednostka jest przypadkiem uciętym, oznacza to, że zdarzenie w przypadku tej jednostki nie wystąpiło. A zatem „bycie przypadkiem uciętym” implikuje brak zdarzenia.
Podgląd pierwszych kilku wierszy uzyskanej ramki danych pokazano niżej:
## months status treatment
## 1 0.231 TRUE A
## 2 1.120 TRUE A
## 3 1.384 TRUE A
## 4 2.076 TRUE A
## 5 2.109 TRUE A
## 6 2.438 FALSE A
## 7 2.735 TRUE A
## 8 2.768 TRUE A
## 9 2.999 TRUE A
## 10 3.559 TRUE ASzacowanie funkcji przeżycia
Do szacowania funkcji przeżycia i modelowania historii zdarzeń wykorzystuje się funkcje w pakiecie survival. Pierwszym krokiem w analizie przeżycia jest utworzenie „zmiennej zależnej”, czyli tego, co chcemy modelować. Pamiętajmy, że przy szacowaniu funkcji przeżycia sama tylko informacja o czasie do wystąpienia zdarzenia nie wystarczy; potrzebna jest nam jeszcze informacja o obserwacjach uciętych. Za pomocą tych dwóch informacji oraz funkcji Surv() tworzy się specjalny obiekt, który następnie wykorzystuje się jako „zmienną zależną” w modelach:
library(survival)
Surv(time = df$months, event = df$status)## [1] 0.231 1.120 1.384 2.076 2.109 2.438+ 2.735 2.768
## [9] 2.999 3.559 3.691 4.251 4.382 4.382 4.580 4.613
## [17] 4.613 4.811 4.910 5.074 5.173 5.272 5.272 5.437
## [25] 5.701 5.799 6.096+ 7.183 7.414 7.941 8.172 8.996
## [33] 9.127 9.193+ 9.786 10.511+ 13.345 13.741 13.839 14.498
## [41] 17.233 17.233+ 19.210 19.573 36.279 36.773+ 37.762 40.398+
## [49] 44.451+ 46.527+ 46.692 1.219 2.768 3.032 3.097 3.625
## [57] 3.691 3.921 4.185 4.284 4.382 4.613 4.811 5.107
## [65] 5.239 5.569+ 5.701 5.898 6.393 6.425 6.887 8.205
## [73] 9.259 10.511 11.170 14.235 15.454 17.102 17.398+ 18.024+
## [81] 20.199+ 20.858 23.890 25.010+ 26.921 35.983+ 41.024+ 43.858+
## [89] 51.305 54.106+ 58.356+ 58.521 62.508+ 66.660+ 70.713+ 75.689+W poleceniu powyżej, podaliśmy dwa argumenty dla funkcji Surv: time, który określa czas do wystąpienia zdarzenia; w naszym przypadku argumentowi temu odpowiada zmienna months z naszej ramki danych df. Z kolei argument event określa, czy zdarzenie wystąpiło; temu argumentowi odpowiada zmienna status. W najprostszym przypadku te dwa argumenty wystarczą. W przyszłości opiszę reguły tworzenia obiektów type survival dla bardziej złożonych przypadków.
Aby oszacować funkcję przeżycia metodą Kaplana-Meiera (opisaną w cytowanym artykule Efrona), użyjemy funkcji survfit:
fit_1 <- survfit(Surv(time = months, event = status) ~ treatment, data = df)Jak widzimy, funkcja survfit ma składnię dość podobną do składni funkcji szacujących modele liniowe. Pierwszym argumentem jest formuła, drugim — nazwa zbioru danych. Po lewej stronie formuły znajduje się nasza „zmienna zależna”, czyli obiekt zawierający informacje potrzebne do oszacowania stopy ryzyka dla każdego momentu czasowego, po prawej stronie zaś — zmienne wyjaśniające. W naszym przypadku jest to tylko jedna zmienna, a mianowicie — rodzaj terapii. Użycie tej zmiennej w formule oznacza, że oszacowania funkcji przeżycia zostaną wykonane osobno dla obu gru pacjentów.
Wynik oszacowania pokazano w okienku poniżej:
## Call: survfit(formula = Surv(time = months, event = status) ~ treatment,
## data = df)
##
## treatment=A
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.231 51 1 0.980 0.0194 0.9431 1.000
## 1.120 50 1 0.961 0.0272 0.9090 1.000
## 1.384 49 1 0.941 0.0329 0.8788 1.000
## 2.076 48 1 0.922 0.0376 0.8507 0.998
## 2.109 47 1 0.902 0.0416 0.8239 0.987
## 2.735 45 1 0.882 0.0453 0.7975 0.975
## 2.768 44 1 0.862 0.0485 0.7719 0.962
## 2.999 43 1 0.842 0.0513 0.7470 0.949
## 3.559 42 1 0.822 0.0539 0.7227 0.934
## 3.691 41 1 0.802 0.0562 0.6989 0.920
## 4.251 40 1 0.782 0.0582 0.6755 0.905
## 4.382 39 2 0.742 0.0618 0.6299 0.873
## 4.580 37 1 0.722 0.0633 0.6076 0.857
## 4.613 36 2 0.681 0.0658 0.5640 0.823
## 4.811 34 1 0.661 0.0668 0.5426 0.806
## 4.910 33 1 0.641 0.0678 0.5214 0.789
## 5.074 32 1 0.621 0.0685 0.5005 0.771
## 5.173 31 1 0.601 0.0692 0.4799 0.753
## 5.272 30 2 0.561 0.0701 0.4393 0.717
## 5.437 28 1 0.541 0.0704 0.4193 0.698
## 5.701 27 1 0.521 0.0706 0.3996 0.680
## 5.799 26 1 0.501 0.0707 0.3800 0.661
## 7.183 24 1 0.480 0.0708 0.3597 0.641
## 7.414 23 1 0.459 0.0707 0.3397 0.621
## 7.941 22 1 0.438 0.0705 0.3199 0.601
## 8.172 21 1 0.418 0.0702 0.3004 0.580
## 8.996 20 1 0.397 0.0697 0.2811 0.560
## 9.127 19 1 0.376 0.0691 0.2621 0.539
## 9.786 17 1 0.354 0.0685 0.2420 0.517
## 13.345 15 1 0.330 0.0678 0.2207 0.494
## 13.741 14 1 0.307 0.0670 0.1998 0.470
## 13.839 13 1 0.283 0.0658 0.1793 0.446
## 14.498 12 1 0.259 0.0644 0.1594 0.422
## 17.233 11 1 0.236 0.0627 0.1400 0.397
## 19.210 9 1 0.210 0.0610 0.1185 0.371
## 19.573 8 1 0.183 0.0587 0.0979 0.344
## 36.279 7 1 0.157 0.0559 0.0783 0.316
## 37.762 5 1 0.126 0.0528 0.0552 0.286
## 46.692 1 1 0.000 NaN NA NA
##
## treatment=B
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1.22 45 1 0.978 0.0220 0.936 1.000
## 2.77 44 1 0.956 0.0307 0.897 1.000
## 3.03 43 1 0.933 0.0372 0.863 1.000
## 3.10 42 1 0.911 0.0424 0.832 0.998
## 3.62 41 1 0.889 0.0468 0.802 0.986
## 3.69 40 1 0.867 0.0507 0.773 0.972
## 3.92 39 1 0.844 0.0540 0.745 0.957
## 4.18 38 1 0.822 0.0570 0.718 0.942
## 4.28 37 1 0.800 0.0596 0.691 0.926
## 4.38 36 1 0.778 0.0620 0.665 0.909
## 4.61 35 1 0.756 0.0641 0.640 0.892
## 4.81 34 1 0.733 0.0659 0.615 0.875
## 5.11 33 1 0.711 0.0676 0.590 0.857
## 5.24 32 1 0.689 0.0690 0.566 0.838
## 5.70 30 1 0.666 0.0704 0.541 0.819
## 5.90 29 1 0.643 0.0716 0.517 0.800
## 6.39 28 1 0.620 0.0727 0.493 0.780
## 6.43 27 1 0.597 0.0735 0.469 0.760
## 6.89 26 1 0.574 0.0742 0.446 0.740
## 8.20 25 1 0.551 0.0747 0.423 0.719
## 9.26 24 1 0.528 0.0750 0.400 0.698
## 10.51 23 1 0.505 0.0752 0.377 0.676
## 11.17 22 1 0.482 0.0752 0.355 0.655
## 14.23 21 1 0.459 0.0750 0.333 0.633
## 15.45 20 1 0.436 0.0747 0.312 0.610
## 17.10 19 1 0.413 0.0742 0.291 0.588
## 20.86 15 1 0.386 0.0742 0.265 0.562
## 23.89 14 1 0.358 0.0739 0.239 0.537
## 26.92 12 1 0.328 0.0735 0.212 0.509
## 51.30 8 1 0.287 0.0749 0.172 0.479
## 58.52 5 1 0.230 0.0789 0.117 0.451Okienko zawiera dwie tabele, osobną dla każdej grupy pacjentów. Wartości w kolumnie z nagłówkiem survival to oszacowania funkcji przeżycia \(S(t)\).
Tabele te są jednak mało czytelne. Znacznie lepiej wygląda wykres sporządzony na podstawie tych tabel:
plot(fit_1, main = "Oszacowania funkcji przeżycia w dwóch grupach pacjentów\nw zależności od rodzaju terapii",
xlab = "Miesiące", ylab = "Prawdopodobieństwo przeżycia", col = c("navyblue",
"firebrick"), xlim = c(0, 80), las = 1)
grid()
legend("top", lty = 1, legend = c("Terapia A", "Terapia B"), col = c("navyblue",
"firebrick"), bty = "n", ncol = 2)
Widoczne na wykresie plusy identyfikują przypadki ucięte. Jeśli nie chcemy ich pokazywać, wystarczy dodać argument mark.time = FALSE do funkcji plot():
plot(fit_1, main = "Oszacowania funkcji przeżycia w dwóch grupach pacjentów\nw zależności od rodzaju terapii",
xlab = "Miesiące", ylab = "Prawdopodobieństwo przeżycia", col = c("navyblue",
"firebrick"), xlim = c(0, 80), las = 1, mark.time = FALSE)
grid()
legend("top", lty = 1, legend = c("Terapia A", "Terapia B"), col = c("navyblue",
"firebrick"), bty = "n", ncol = 2)