Breviar teoretic*

Vom aborda tematica estimarii functiei de supravietuire asociata unui set de date printr-o noua perspectiva, si anume a modelelor parametrice. In acest situatie se considera ca datele de timp sunt distribuite dupa o anumita distributie de probabilitate (exponentiala sau Weibull) si se vor estima parametrii acestei distributii, folosind metoda estimatorului de verosimilitate maxima.
Metoda folosita poarta denumirea de metoda parametrica pentru estimarea S(t) si consta, mai exact, in impunerea unei anumite forme algebrice a funcţiei S(t), formă ce depinde de unul sau mai mulţi, care vor fi estimaţi.

Distributiile de probabilitate exponentiala si Weibull
Distributiile de probabilitate exponentiala si Weibull

Recapitulare distributiile Exponentiala si Weibull

Distributia de probabilitate de tip exponential de parametru \(\lambda\) (notatie \(T\sim Exp(\lambda)\)): \[f(t)=\lambda e^{-\lambda t}\rightarrow S(t)=e^{-\lambda t}, \quad \lambda>0, t\geq 0\].

Distributia de probabilitat Weibull de parametrii \(\lambda\) si \(\alpha\) ( (notatie \(T\sim Weibull(\alpha, \lambda)\)): \[f(t)=\alpha\lambda t^{\alpha-1} e^{-\lambda t^{\alpha}}\rightarrow S(t)=e^{-\lambda t^{\alpha}}, \quad \alpha, \lambda>0, t\geq 0\]

Observatie: Distributia de probabilitat Weibull pentru \(\alpha=1\) este distributia exponentiala.

Estimarea parametriilor unui model

Estimarea parametrilor se face folosind metoda verosimilitii maxime. Vom folosi urmatoarele notatii:

  • \(T_i\)-v.a. ce inregistreaza timpul de suprav al persoanei \(i\), avand distributia de probabilitate \(f_i(t)\), iar functia de repartitie fiind notata cu \(F_i\)

  • \(C_i\)-v.a. ce inregistreaza momentul de cenzurare al persoanei \(i\), cu distributia de probabilitate \(g_i(t)\), iar functia de repartitie fiind notata \(G_i(t)\)

  • presupunem ca variabilele \(T_i\) si \(C_i\) sunt independente;

  • timpul de supravietuire observat este \(\min (T_i, C_i)\);

  • contributia la functia de suprav a subiectului \(i\) cu timpul de suprav \(t_i\) si statutul \(\delta_i=1\) este \(f_i(t_i)(1-G_i(t_i))\)

  • contributia la functia de suprav a subiectului \(i\) cu timpul de suprav \(t_i\) si \(\delta_i=0\) este \((1-F_i(t_i))g_i(t_i)\)

Estimarea parametriilor prin metoda MML se face parcurgand urmatorii pasi:

-definim functia de verosimilitate maxima \(L(\lambda, t_i)=\prod\limits_{i=1}^n [f_i(t_i)(1-G_i(t_i))]^{\delta_i}[(1-F_i(t_i))g_i(t_i)]^{1-\delta_i}\)

-\(L(\lambda, t_i)\approx (f_i(t_i))^{\delta_i}(1- F_i(t_i))^{1-\delta_i}\)

-se determina punctul de maxim al functiei \(l(\lambda, t_i)=log L(\lambda, t_i)\)

-se determina estimatorul de verosimilitate \(\hat{\lambda}\) si astfel vom avea expresia \(\hat{S}(t)=1-\hat{F}(t), t\geq 0\)

Estimarea parametrilor modelului exponential

Rezultate teoretice

Consideram cazul in care timpul de supravietuire \(T\) este distribuit exponential de parametru \(\lambda\): \[f(t)=\lambda e^{-\lambda t}, F(t)=1- e^{-\lambda t}\] Vom aplica metoda determinarii estimatorului de verosimilitate maxima in acest caz. Avem:

-\(ln L(\lambda, t_i)=ln \lambda \sum\limits_{i=1}^n\delta_i-\lambda \sum\limits_{i=1}^n t_i\)

-\(\hat{\lambda}=\frac{\sum\limits_{i=1}^n\delta_i}{\sum\limits_{i=1}^n t_i}=\frac{\text{ numarul obs cu eveniment }}{\text{ numarul total ani-pacienti-}}\)

-\(\hat{S}(t)=e^{-\hat{\lambda} t},\quad t\geq 0\)

Aplicatie

Se considera urmatoarele observatii: 1,1,2,3+,3,5,8,10,16+,18. Ne propunem asocierea unui model exponential, deci estimarea functiei de supravietuire asociata acestui model.

Solutie analitica

Prin folosrea metodei prezentate anterior, obtinem

\(\hat{\lambda}=\frac{8}{1+2+\ldots 67}=0.1194\)

Deci, \(\hat{S}(t)=e^{-0.1194 t},\quad t\geq 0\).

Solutie in R
library(survival)
time<-c(1,1,2,3,3,5,8,10,16,18)
status<-c(1,1,1,0,1,1,1,1,0,1)

Vom crea sirul de inregistrari dupa forma specifica din AS:

my.surv<-Surv(time,status)
my.surv
##  [1]  1   1   2   3+  3   5   8  10  16+ 18

Vom determina probabilitatile de suparvietuire prin metoda Kapan-Meier

my.fit<-survfit(my.surv~ 1)
my.fit
## Call: survfit(formula = my.surv ~ 1)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 10      8      5       2      NA
summary(my.fit)
## Call: survfit(formula = my.surv ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     10       2     0.80   0.126       0.5868        1.000
##     2      8       1     0.70   0.145       0.4665        1.000
##     3      7       1     0.60   0.155       0.3617        0.995
##     5      5       1     0.48   0.164       0.2458        0.938
##     8      4       1     0.36   0.161       0.1498        0.865
##    10      3       1     0.24   0.145       0.0732        0.786
##    18      1       1     0.00     NaN           NA           NA

Vom reprezenta grafic functia de suparvietuire si vom compara aceasta curba cu curba functie de supravietuire obtinuta prin metoda parametrica \(\hat{S}(t)=e^{-0.1194 t},\quad t\geq 0\). Comanda par(new=TRUE) ne foloseste pentru a reprezenta in acelasi sistem de axe cele doua curbe.

plot(my.fit, xlab="timp", ylab="probabilitatea de supravietuire", main="Curba Kaplan-Meier")
par(new=TRUE)
curve(exp(-0.1194*x),0,18,xlab="",ylab="", yaxt="n")

Graficul de mai sus ne da o prima intuitie despre cat de bun este modelul exponential, asezat peste curba de supravietuire determinata prin metoda Kaplan-Meier.

Parametru \(\lambda\) se paote estima si folosind R, prin comanda

modexp=survreg(Surv(time,status)~1,dist="exp")
summary(modexp)
## 
## Call:
## survreg(formula = Surv(time, status) ~ 1, dist = "exp")
##             Value Std. Error    z       p
## (Intercept) 2.125      0.354 6.01 1.8e-09
## 
## Scale fixed at 1 
## 
## Exponential distribution
## Loglik(model)= -25   Loglik(intercept only)= -25
## Number of Newton-Raphson Iterations: 4 
## n= 10

Din output se observa ca se obtine Intercept=2.13, deci parametrul \(\hat{\lambda}=exp(-2.13)=0.119\)

exp(-2.13)
## [1] 0.1188373

Deci, ca regula generala, din output vom obtine valoarea estimata a parametrului modelului exponential ca fiind \[\hat{\lambda}=exp(-intercept)\]

Estimarea parametrilor modelului Weibull

Rezultate teoretice

Vom aplica metoda prezentata mai sus pentru a estima parametrii distribuiei Weibull data prin densitatea de probabilitate \(f(t)=\alpha\lambda t^{\alpha-1} e^{-\lambda t^{\alpha}}\), cu functia de repartitie \(F(t)=1- e^{-\lambda t^{\alpha}}\).

In acest caz functia de verosimilitate maxima \(L(\lambda,\alpha)\) are o expresie mai complicata, depinde de doi parametrii \(\alpha, \lambda\) si, insasi metoda de estimare a acestor doi parametrii este mai dificila.

Derivarea partiala in functie de parametrii $$ i \(\alpha\) a funcia \(l(\lambda, \alpha)=ln(L)\) ne conduce la rezolvarea unui sistem, pentru care doar solutia numerica este posibila (sistemul este neliniar). Deci, involuntar, va trebui sa folosim calculatorul pentru estimarea parametriilor modelului Weibull. Dupa ce se determina valoarea numerica a acestor parametrii, vom avea ca functia de suparvietuire are expresia

\[\hat{S}(t)=e^{-\hat{\lambda} t^{\hat{\alpha}}},\quad t\geq 0\] #### Aplicatie

Se considera observaiile: 24+,16+,8,19,10,10,8+,5,17,20,10. Ne propunem asocierea unui model de tip Weibul acestui set de date, deci estimarea parametriilor funciei de supravieuire pentru acest model.

Solutie analitica

Prin metoda de versosimilitate maxima si rezolvarea numerica a sistemului asociat pentru determinarea parametriilor \(\lambda\) si \(\alpha\) avem:

\(\hat{\alpha}=2.2451,\hat{\lambda}=0.0015\) Deci, \[\hat{S}(t)=e^{-0.0015 t^{2.2451}},\quad t\geq 0\]

Solutie in R

Introducem datele problemei

T<-c(24,16,8,19,10,8,5,17,20,10)
Status<-c(0,0,1,1,1,0,1,1,1,1)
My.surv<-Surv(T,Status)
My.surv
##  [1] 24+ 16+  8  19  10   8+  5  17  20  10

Vom determina probabilitatile de supravituire prin metoda 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
plot(My.fit, xlab="timp", ylab="probabilitatea de supravietuire", main="Curba Kaplan-Meier")
par(new=TRUE)
curve(exp(-0.0015*x^(2.2451)),0,20,xlab="",ylab="", yaxt="n")

Pentru a obtine din R estimatorii pentru parametrii \(\lambda\) si \(\alpha\) putem folosi comanda

wei=survreg(My.surv~1,dist="w")
summary(wei)
## 
## Call:
## survreg(formula = My.surv ~ 1, dist = "w")
##              Value Std. Error     z      p
## (Intercept)  2.882      0.168 17.12 <2e-16
## Log(scale)  -0.809      0.309 -2.62 0.0088
## 
## Scale= 0.445 
## 
## Weibull distribution
## Loglik(model)= -25.4   Loglik(intercept only)= -25.4
## Number of Newton-Raphson Iterations: 6 
## n= 10

Din output, avem ca

\(scale= 0.445, \alpha=\frac{1}{scale}=2.247, \lambda=exp(-intercept * \alpha)=\exp(-2.882* 2.247)\)

alpha=1/0.445
alpha
## [1] 2.247191
lambda=exp(-2.882*2.247)
lambda
## [1] 0.001540183

Deci, regula generala pentru a estima parametrii distribuitie Weibull din output R este

\(\alpha=\frac{1}{scale}, \lambda=exp(-intercept * \alpha)\).