- Analiza prze¿ycia
- Cenzurowanie danych
- Funkcja prze¿ycia
- Metody estymacji funkcji prze¿ycia
- Testy statystyczne w analizie prze¿ycia
- Zadanie punktowane
Podczas zajêæ bêd± wykorzystywane biblioteki survival oraz survminer.
Analiza prze¿ycia
Analiza prze¿ycia to badanie statystyczne, które przedmiotem jest czas jaki up³ynie do wyst±pienia po raz pierwszy pewnego zdarzenia. Tym zdarzeniem mo¿e byæ ¶mieræ pacjenta, kolejne ocielenie siê krowy, czy migracja do innego miasta. Analiza prze¿ycia jest równie¿ wykorzystywana w naukach aktuarialnych przy modelowaniu ¶mierci ubezpieczonego lub momentu wyst±pienia szkody. Podstawowym pytaniem, na które próbuje odpowiedzieæ analiza prze¿ycia jest to, czy obserwacje w badanych grupach s± opisywane przez taki sam czas wyst±pienia zdarzenia.
Cenzurowanie danych
Charakterystyczne dla analizy prze¿ycia jest to, ¿e dane s± cenzurowane (ang. censored), czyli inaczej uciête. Je¶li badany osobnik nie dozna³ oczekiwanego zdarzenia do koñca trwania eksperymentu lub wycofa³ siê z niego przed jego zakoñczeniem, to taki osobnik jest cenzurowany prawostronnie. Innymi s³owy, w przypadku obserwacji cenzurowanych prawostronnie nie wiemy kiedy i czy w ogóle dosz³o do interesuj±cego nas zdarzenia. Cenzurowanie lewostronnie ma miejsce wtedy kiedy nie znamy dok³adnego czasu punktu pocz±tkowego, od którego mierzymy czas do wyst±pienia oczekiwanego zdarzenia.
W pakiecie R do analizy prze¿ycia wykorzystuje siê pakiet survival. Podstawow± funkcj± w pakiecie jest Surv, która zwraca Survival Object. Zapis \(t+\) oznacza, ¿e u danego osobnika nie wyst±pi³o zdarzenie do czasu \(t\).
library(survival)
obs.t <- c(7,6,6,5,2,4)
obs.cens <- c(0,1,0,0,1,1)
obs.surv <- Surv(obs.t,obs.cens)
obs.surv## [1] 7+ 6 6+ 5+ 2 4Funkcja prze¿ycia
Kluczowym pojêciem w analizie prze¿ycia jest funkcja prze¿ycia \(S(t) = \mathbb{P}(T>t)\), która okre¶la prawdopodobieñstwo, ¿e zdarzenie nie zostanie zaobserwowane do czasu \(t\). Je¶li czas do wyst±pienia zdarzenia \(T\) jest zmienn± losow± z rozk³adu \(f\) o dystrybuancie \(F\), to funkcja prze¿ycia dana jest wzorem \(S(t) = 1 - F(t)\).
Funkcja hazardu okre¶la warunkowe prawdopodobieñstwo zaobserwowania zdarzenia w jednostce czasu \(\Delta t\), je¶li nie wyst±pi³o ono do chwili \(t\), przypadaj±ce na jednostkê czasu \(\Delta t\)
\(\qquad h(t) = \lim\limits_{\Delta t \to 0} \frac{\mathbb{P}(t \leqslant T < t + \Delta t | T \geqslant t)}{\Delta t}\).
Skumulowany hazard to suma hazardu do chwili \(t\)
\(\qquad \Lambda(t) = \int\limits_{0}^t h(u) du\).
Z powy¿szych wzorów mo¿na wyprowadziæ zale¿no¶ci pomiêdzy funkcj± prze¿ycia a funkcj± hazardu i skumulowanym hazardem:
\(\qquad h(t) = \lim\limits_{\Delta t \to 0} \frac{\mathbb{P}(t \leqslant T < t + \Delta t | T \geqslant t)}{\Delta t} = \frac{\lim\limits_{\Delta t \to 0} \frac{\mathbb{P}(t \leqslant T < t + \Delta t)}{\Delta t}}{\mathbb{P}(T \geqslant t)} = \frac{f(t)}{S(t)} = -\frac{S'(t)}{S(t)} = -\frac{d\log{S(t)}}{dt}\),
\(\qquad \Lambda(t) = \int\limits_{0}^t h(u) du = \log{S(0)} - \log{S(t)} = -\log{S(t)}\),
\(\qquad S(t) = \exp(-\Lambda(t))\).
Metody estymacji funkcji prze¿ycia
W praktyce zale¿y nam na okre¶leniu funkcji prze¿ycia dla danych eksperymentalnych aby na jej podstawie móc obliczaæ hazard oraz porównywaæ funkcje prze¿ycia dla ró¿nych grup. Wyró¿niamy metody nieparametryczne oraz parametryczne estymacji funkcji prze¿ycia. Najpopularniejszymi metodami nieparametrycznymi s± estymatory Kaplana-Meiera oraz Flemingtona-Harringtona.
Estymator Kaplana-Meiera
Metoda ta zak³ada dyskretny czas oraz uwzglêdnia cenzurowanie. Dana jest wzorem \(\hat{S}(t) = P(t) \cdot P(t-1)\), gdzie \(P(t) = 1-\frac{d(t)}{r(t)}\) jest nazywane prawdopodobieñstwem prze¿ycia w chwili \(t\), \(d(t)\) jest liczb± wyst±pieñ zdarzenia w chwili \(t\), natomiast \(r(t)\) oznacza liczbê obserwacji objêtych ryzykiem. Poni¿sza tabela przedstawia przyk³ad liczenia estymatora Kaplana-Meiera.
| Czas | Obserwacje | Zdarzenia | Cenzurowane | Ryzyka | Pr. prze¿ycia | Funkcja prze¿ycia |
|---|---|---|---|---|---|---|
| \(t_j\) | \(n_j\) | \(d_j\) | \(c_j = n_j - n_{j+1} - d_j\quad\) | \(r_j = n_j - c_j\quad\) | \(P_j = (r_j-d_j)/r_j\qquad\) | \(S_j = P_j \cdot P_{j-1}\) |
| 0 | 41 | 3 | \(41 - 36 - 3 = 2\) | \(41 - 2 = 39\) | \((39-3)/39 = 0.92\) | \(0.92 \cdot 1.00 = 0.92\) |
| 1 | 36 | 2 | \(36 - 33 - 2 = 1\) | \(36 - 1 = 35\) | \((35-2)/35 = 0.94\) | \(0.94 \cdot 0.92 = 0.86\) |
| 2 | 33 | 1 | \(33 - 30 - 1 = 2\) | \(33 - 2 = 31\) | \((31-1)/31 = 0.97\) | \(0.97 \cdot 0.86 = 0.83\) |
| 3 | 30 | 2 | \(30 - 26 - 2 = 2\) | \(30 - 2 = 28\) | \((28-2)/28 = 0.93\) | \(0.93 \cdot 0.83 = 0.77\) |
Estymator Flemingtona-Harringtona
Ten estymator wykorzystuje zale¿no¶æ pomiêdzy funkcj± prze¿ycia i skumulowan± funkcj± hazardu \(\hat{S}(t) = \exp(-\hat{\Lambda}(t))\), gdzie \(\hat{\Lambda}(t)\) jest estymatorem skumulowanego hazardu Nelsona-Aalena
\(\qquad \hat{\Lambda}(t) = \sum\limits_{\tau<t} \frac{d(\tau)}{r(\tau)},\quad d(t)\) i \(r(t)\) oznaczaj± to samo co w przypadku estymatora Kaplana-Meiera.
Metody parametryczne
W przypadku kiedy dysponujemy informacj± z jakiej rodziny \(f(t)\) rozk³adów pochodz± nasze obserwacje (czasy do wyst±pienia zdarzenia) to mo¿emy wyznaczyæ dystrybuantê \(F(t)\), a st±d \(h(t)\), \(\Lambda(t)\) oraz \(S(t)\). Innymi s³owy dopasowuj±c rozk³ad \(f(t)\) do danych eksperymentalnych mo¿emy ³atwo wyznaczyæ funkcjê prze¿ycia. Najczê¶ciej wykorzystywanymi rodzinami rozk³adów w analizie prze¿ycia s±: wyk³adniczy, Weibulla, Gompertza i log-logistyczny.
- Rozk³ad wyk³adniczy ma funkcjê gêsto¶ci prawdopodobieñstwa postaci \(f(t) = \lambda \exp(-\lambda t)\), dla \(t \geqslant 0, \lambda > 0\). Dostajemy wiêc:
\(\qquad h(t) = \lambda\),
\(\qquad \Lambda(t) = \lambda t\),
\(\qquad S(t) = \exp(-\lambda t)\)
- Rozk³ad Weibulla jest dwuparametryczny i jego gêsto¶æ zapisuje siê jako \(f(t)=k\lambda(\lambda t)^{k-1} \exp(-(\lambda t)^k)\), dla \(t \geqslant 0, \lambda,k > 0\), co daje:
\(\qquad h(t) = k\lambda(\lambda t)^{k-1}\),
\(\qquad \Lambda(t) = (\lambda t)^k\),
\(\qquad S(t) = \exp(-(\lambda t)^k)\)
Parametr \(\lambda\) jest parametrem skali, natomiast \(k\) jest parametrem kszta³tu. Rozk³ad wyk³adniczy jest szczególnym przypadkiem rozk³adu Weibulla dla \(k=1\). Gdy \(k=1\) to hazard jest sta³y, je¶li \(k>1\) to rosn±cy, a dla \(k<1\) malej±cy w czasie.
Testy statystyczne w analizie prze¿ycia
W pakiecie R estymacjê funkcji prze¿ycia przy u¿yciu estymatora Kaplana-Meiera mo¿na wykonaæ wywo³uj±c funkcjê survfit. W dalszych przyk³adach bêdziemy siê pos³ugiwaæ zbiorem ovarian z pakietu survival zawieraj±cym dane dotycz±ce dwóch sposobów leczenia raka jajnika.
rak <- ovarian
rak.surv <- Surv(time = ovarian$futime, event = ovarian$fustat)
rak.surv## [1] 59 115 156 421+ 431 448+ 464 475 477+ 563 638
## [12] 744+ 769+ 770+ 803+ 855+ 1040+ 1106+ 1129+ 1206+ 1227+ 268
## [23] 329 353 365 377+rak.fit <- survfit(rak.surv~rx, data=rak)
summary(rak.fit)## Call: survfit(formula = rak.surv ~ rx, data = rak)
##
## rx=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 59 13 1 0.923 0.0739 0.789 1.000
## 115 12 1 0.846 0.1001 0.671 1.000
## 156 11 1 0.769 0.1169 0.571 1.000
## 268 10 1 0.692 0.1280 0.482 0.995
## 329 9 1 0.615 0.1349 0.400 0.946
## 431 8 1 0.538 0.1383 0.326 0.891
## 638 5 1 0.431 0.1467 0.221 0.840
##
## rx=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 353 13 1 0.923 0.0739 0.789 1.000
## 365 12 1 0.846 0.1001 0.671 1.000
## 464 9 1 0.752 0.1256 0.542 1.000
## 475 8 1 0.658 0.1407 0.433 1.000
## 563 7 1 0.564 0.1488 0.336 0.946Przy pomocy funkcji ggsurvplot mo¿na wykonaæ wykresy funkcji prze¿ycia, a ustawiaj±c parametr pval = TRUE otrzymamy p-warto¶æ dla testu statystycznego, którego hipoteza zerowa g³osi, ¿e funkcje prze¿ycia w ró¿nych grupach nie ró¿ni± siê statystycznie od siebie.
library(survminer)
ggsurvplot(rak.fit, data=rak, pval = TRUE)
Wysoka p-warto¶æ wskazuje na to, ¿e nie ma podstaw do odrzucenia hipotezy zerowej.
Zadanie punktowane
Dokoñcz analizê prze¿ycia dla danych ovarian. Zbadaj czy istnieje zale¿no¶æ miêdzy funkcj± prze¿ycia a czynnikami age, resid.ds i ecog.ps. Czynnik age nale¿y uprzednio zdyskretyzowaæ (wystarcz± dwie grupy, mo¿na zrobiæ histogram aby wybraæ punkt podzia³u).
