1. Breviar teoretic

Reamintim ca functia de supravietuire este definita prin formula \[S(t_i):=P(T>t_i). \] Ne propunem sa estimam valorile acestei functii pentru un set de intregistrari date \(t_i,i=1,2,\ldots, k\).

Formula de estimare (prin metoda Kaplan-Maier) a probabilitatii de supravietuire peste momentul de timp \(t\) este: \[\hat{S}(t_i)=\prod\limits_{j=1}^i(1-\frac{d_j}{n_j})\] unde:

Vom initializa \(t_0=0\), \(d_0=1\).

Pentru obtinerea acestei formule se foloseste regula produsului \(P(A\cap B)=P(A|B)P(B)\). Se obtine:

Deci, \(S(t_i)=P(T>t_i)=P(T>t_i|T>t_{i-1})P(T>t_{i-1})=P(T>t_i|T>t_{i-1})S(t_{i-1}),\, i=1,2\ldots k\)

Notand \(\pi_i=:P(T>t_i|T>t_{i-1})\), se obtine \[S(t_i)=\prod\limits_{j=1}^i \pi_j\]

Se aplica metoda verosimilitii maxime pentru estimarea \(\pi_i\). Se obtine \[\hat{\pi}_i=1-\frac{d_i}{n_i}\] Deci, \[\hat{S}(t_i)=\prod\limits_{j=1}^i \hat{\pi}_j=\prod\limits_{j=1}^i(1-\frac{d_j}{n_j})\]

2. Aplicatie

O companie realizeaza un studiu clinic pentru a testa eficacitatea unui nou model de valvă pentru inimă. Timpii de supravieţuire (în luni) pentru 10 pacienţi cu implant de valvă sunt înregistraţi: 24+ 16+ 8 19 10 8+ 5 17 20 10. Să se estimeze probabilitatileia de supravietuire.

3. Rezolvare analitica

Pentru inceput vom ordonam crescator timpii de supravietuire inregistrati: \[5<8<10<16<17<19<20<24\] Pentru fiecare moment de timp \(t_i\) vom centraliza \(n_i\) numarul de persoane la risc (persoanele in viata pana la momentul \(t_i\)-exclus acest moment- si nu au fost cenzurate, precum si \(d_i\)-numarul de persoane pentru care a avut loc evenimentul la momentul \(t_i\) si \(c_i\)-numarul de persoane cenzurate la \(t_i\). Vom calcula pentru fiecare moment \(t_i\) rata \(1-\frac{d_i}{n_i}\) si estimatorul \(\hat{S}(t_i)=\prod\limits_{j=1}^i(1-\frac{d_j}{n_j})\).

De exemplu, pentru \(t_0=0\), folosind initializarile \(t_0=0\), \(d_0=1\) si faptul ca la inceputul studiului toate persoanele din studiu sunt la risc, avem ca \(n_0=10\). In plus, \(c_0=0\) (nu sunt persoane cenzurate la \(t_0\). Astfel, \(S(t_0)=P(T>0)=1-\frac{0}{10}=1,\) deci \(S(t_0)=0\).

Mai departe, pentru \(t_1=5\), avem \(n_1=10\) si \(d_1=1,\) iar \(c_1=0\). Deci, \(1-\frac{d_1}{n_1}=1-\frac{1}{10}=0.9\). Astfel \(S(t_1)=1-\frac{d_1}{n_1} S(t_0)=1\times 0.9=0.9\).

In continuare, pentru \(t_2=8\), avem \(n_2=9\) si \(d_2=1,\) iar \(c_1=1\). Deci, \(1-\frac{d_2}{n_2}=1-\frac{1}{9}=0.89\). Astfel \(S(t_2)=1-\frac{d_2}{n_2} S(t_1)=0.9\times 0.89=0.8\).

La momentul \(t_3=10\), unde avem \(n_3=7\) (din persoanele aflate la risc anterior, adica 9, s-a scazut o persoana cu eveniment si o persoana cenzurata). Avem ca \(d_3=2,\) iar \(c_1=0\). Deci, \(1-\frac{d_3}{n_3}=1-\frac{2}{7}=0.71\). Astfel \(S(t_3)=1-\frac{d_3}{n_3} S(t_2)=0.8\times 0.71=0.57\).

Mai jos sunt prezentate calculele analitice pentru fiecare moment \(t_i\). Se observa ca pentru \(t_i=16\), probabilitatea de supravietuire nu se modifica, pentru ca, in realitate, la acest momnent de timp nu s-a produs un eveniment (inregistrarea 16 apare cu semn plus, deci consemneaza o inregistrare cenzurata).

as.table(matrix(c(0,10,0,0,1,1,5,10,1,0,0.9,0.9,8,9,1,1,0.89,0.8,10,7,2,0,0.71,0.57,16,5,0,1,1,0.57,17,4,1,0,0.75,0.43,19,3,1,0,0.67,0.29,20,2,1,0,0.5,0.15,24,1,0,1,1,0.15), nrow=9, byrow=TRUE, dimnames=list(c(),
 c("t_i", "n_i", "d_i", "(c_i la t_i)","(1-d_i/t_i)","Estimator S"))))
##     t_i   n_i   d_i (c_i la t_i) (1-d_i/t_i) Estimator S
## A  0.00 10.00  0.00         0.00        1.00        1.00
## B  5.00 10.00  1.00         0.00        0.90        0.90
## C  8.00  9.00  1.00         1.00        0.89        0.80
## D 10.00  7.00  2.00         0.00        0.71        0.57
## E 16.00  5.00  0.00         1.00        1.00        0.57
## F 17.00  4.00  1.00         0.00        0.75        0.43
## G 19.00  3.00  1.00         0.00        0.67        0.29
## H 20.00  2.00  1.00         0.00        0.50        0.15
## I 24.00  1.00  0.00         1.00        1.00        0.15

In cele ce urmeaza ne propunem sa obtinem probabilitatile de supravietuire peste un moment \(t_i\) folosind programul R.

4. Rezolvare in R

Se incarca pachetul folosit pentru Analiza Supravietuirii.

library(survival)
## Warning: package 'survival' was built under R version 4.0.5

Se Introduc datele problemei.

Se incarca timpii observati pana la aparitia evenimentului

time<-c(24,16,8,19,10,8,5,17,20,10)
time
##  [1] 24 16  8 19 10  8  5 17 20 10

Se incarca tipul de inregistrare pentru timp: cenzurat cu 0 (apare cu + in text) si 1 pentru data necenzurata

status<-c(0,0,1,1,1,0,1,1,1,1) 
status
##  [1] 0 0 1 1 1 0 1 1 1 1

Vom structura datele intr-un tabel, folosind urmatoarea comanda:

df<-data.frame(cbind(indiv=1:10,time,status))
df

Se folosesc comenzi din R specifice Analizei Supravietuirii

Se creeaza un vector de date folosit in Analiza Suprav in care data cenzurata apare cu +

my.surv<-Surv(time,status) 
my.surv
##  [1] 24+ 16+  8  19  10   8+  5  17  20  10

Calculul probabilitatilor de supravietuire

Urmatoarea comanda calculeaza probabilitatile de supravietuire conform formulei Kaplan Meier

my.fit<-survfit(my.surv~1) 
my.fit
## Call: survfit(formula = my.surv ~ 1)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 10      7     17      10      NA
summary(my.fit)
## Call: survfit(formula = my.surv ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5     10       1    0.900  0.0949       0.7320        1.000
##     8      9       1    0.800  0.1265       0.5868        1.000
##    10      7       2    0.571  0.1638       0.3258        1.000
##    17      4       1    0.429  0.1743       0.1931        0.951
##    19      3       1    0.286  0.1647       0.0923        0.884
##    20      2       1    0.143  0.1303       0.0239        0.854

Reprezentare grafica a curbei de supravietuire

Vom reprezenta grafic curba de supravietuire Kaplan-Meier (aici cu + sunt marcate momentele de timp in care are loc cenzurarea)

plot(survfit(my.surv~1),col="red",xlim = c(0, 25), ylim = c(0, 1),mark.time=TRUE, xlab="timp",ylab="probabilitati de supravietuire", conf.int=FALSE) 

title('Curba Kaplan-Maier')

Adugam curbei de supravietuire intervalele de incredere pentru probabilitatile de supravietuire determinate

plot(survfit(my.surv~1),col="red",xlim = c(0, 25), ylim = c(0, 1),mark.time=TRUE, xlab="timp",ylab="probabilitati de supravietuire", conf.int=TRUE) 

title('Curba Kaplan-Maier')

Estimarea timpului x de supravietuire

Folosind comanda de mai jos vom estima care este probabilitatea de supravietuire peste momentul \(t_m\), adica valoarea functiei \(\hat{S}(t_m):=P(T>t_m)\). De exemplu, vom estima \(\hat{S}(7)\), adica probabilitatea de supravietuire peste 7 zile (o saptamana):

summary(survfit(Surv(time, status) ~ 1), times = 7)
## Call: survfit(formula = Surv(time, status) ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     7      9       1      0.9  0.0949        0.732            1

Valoarea obtinuta este in acord cu faptul ca \(P(T>7)=1-\frac{1}{10}=0.9\); s-a calculat ca inainte de momentul 7, o singura persoana a decedat, din cei 10 aflati la risc.

Ne propunem sa calculam estimatorul \(\hat{S}(10)\).

summary(survfit(Surv(time, status) ~ 1), times = 10)
## Call: survfit(formula = Surv(time, status) ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    10      7       4    0.571   0.164        0.326            1

Folosind metoda Kaplan-Maier, se obtine ca \(\hat{S}(10)=0.571\).

O abordare naiva a calcului acestei probabilitati de supravietuire (fara a tine cont de faptul ca inainte de momentul t=10 au existat si observatii cenzurate) ne-ar conduce la rezultatul ca probabilitatea de supravietuire peste 10 zile este \((1-2/10)=0.8\) In realitate, se produce o supraestimare a probabilitatii de supravietuire, daca sunt ignorate informatii legate de cenzurarea inregistrarii ( si nu este adaptat la numitor numarul de persoane aflate la risc).

Estimarea timpului mediu de supravietuire

Un indicator statistic de interes in Analiza Supravituirii este timpul median de supravietuire (mediana). Mediana este momentul de timp \(t_m\) pentru care probabilitatea de supravietuire este 0.5 (adica \(S(t_m)=0.5\)).

Timpul median de supravietuire se poate obtine folosind comanda

survfit(Surv(time, status) ~ 1)
## Call: survfit(formula = Surv(time, status) ~ 1)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 10      7     17      10      NA

In consecinta \(t_m=17\) zile.

plot(survfit(my.surv~1),col="red",xlim = c(0, 25), ylim = c(0, 1),mark.time=TRUE, xlab="timp",ylab="probabilitati de supravietuire", conf.int=TRUE) 

title('Curba Kaplan-Maier')
abline(h=0.5,v=17)

Ne propunem sa comparam timpul median de supravietuire obtinut \(t_m=17\), cu timpul median al inregistrailor, din care vom elimina pe cele cenzurate. Pentru acest lucru, din setul de date vom elimina inregistrile cenzurate si vom calcula timpul median al momentelor de timp in care are loc evenimentul urmarit.

sub1=df[status == "1", ]  
sub1
summary(sub1$time)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00    9.00   10.00   12.71   18.00   20.00

Se observa ca in acest mod se obtine un timp median \(t_M=10\), deci prin acest procedeu are loc o subestimare a timpului median de supravietuire. Putem spune ca prin ignorarea observatiilor cenzurate, s-a pierdut informatie.