Breviar teoretic*

Notam cu \(T\) variabila aleatoare ce inregistreaza timpul pana la evenimentul urmarit, iar cu \(X_1\ldots X_m\) alte variabile inregistrate in studiu (de exemplu, varsta, gen, etc). Ne propunem a studia modul cum aceste variabile influenteza timpul de supravietuire, prin definirea unei relatii matematice.

Un model parametric de regresie pentru \(T\) este o relatie matematica de forma \[\ln T=\beta_0+\beta_1 X_1+\ldots +\beta_m X_m +\sigma\varepsilon,\]

unde parametrii \(\beta_i\) se numesc coeficienti de regresie, iar \(\sigma\) este o constanta reala si \(\varepsilon\) -eroare aleatoare.

Exista doua modele parametrice de regresie folosite mai des:

-modelul exponential;

-modelul Weibull.

Modelul exponential

Consideram urmatoarele ipoteze:

-\(\sigma=1\);

-variabila aleatoare (eroarea) \(\varepsilon\) urmeaza o distributie cu valori-extreme, adica \(f_{\varepsilon}(x)=e^{x-e^x}, -\infty<x<\infty\)

-variabila aleatoare \(T\) urmeaza o distributie exponentiala \(f(t)=\lambda e^{-\lambda t}, t\geq 0\)

Vom estima parametrul \(\lambda\) prin metoda verosimilitatii maxime aplicata pentru acest model, avand n observatii de forma \((t_i, \delta_i, x_{1i},\ldots,x_{mi}),\) unde \(t_i\) urmeaza o distributie exponentiala de parametru \(\lambda\).

Se obtine \[\hat{\lambda}=\exp\{-(\hat{\beta_0}+\hat{\beta_1}x_1+\ldots+\hat{\beta_m}x_m)\}.\] In consecinta, \(\hat{S}(t)=\exp\{-\hat{\lambda}t\},\quad t\geq 0\)

Observatii:

-timpul mediu de supravietuire \[E(T)=\frac{1}{\hat{\lambda}}=\exp\{\hat{\beta_0}+\hat{\beta_1}x_1+\ldots +\hat{\beta_m}x_m\}\]

-pentru o variabila numerica \(X_i\), cantitatea \(e^{\beta_i}\) reprezinta modificarea relativa in timpul mediu de supravietuire adusa de fiecare crestere cu o unitate a covariatei \(X_i\).

-pentru o variabila categorica \(X\) cu \(l\) nivele, modelul de regresie va contine \(l-1\) variabile dummy (\(y_i=1,\) pentru \(x=i\) si \(0\), altfel ). Avem: \[ E(T)=\frac{1}{\hat{\lambda}}=\exp\{\hat{\beta_0}+\hat{\beta_1}y_1+\ldots +\hat{\beta_{l-1}}y_{l-1}+\text{alti termeni }\}\] Astfel, cantitatea \(\exp\{\beta_i-\beta_j\}\) reprezinta raportul dintre timpul de supravietuire cu variabila \(X\) la nivelul \(i\) si aceeasi variabila, la nivelul \(j\).

Modelul Weibull

Consideram urmatoarele ipoteze:

-\(\sigma\) se estimeaza din date;

-\(\varepsilon\) este o variabila aleatoare urmeaz o distribuie cu valori-extreme \(f_{\varepsilon}(x)=e^{x-e^x}, -\infty<x\infty\)

-\(T\) este o variabila aleatoare urmeaza o distributie Weibull \(f(t)=\alpha \lambda t^{\alpha-1} e^{-\lambda t^{\alpha}}, t\geq 0\)

Parametrii \(\alpha\) si \(\lambda\) se vor estia prin metoda MML. Avem urmatoarele rezultate:

\[\hat{\alpha}=\frac{1}{\hat{\sigma}}\]

\[\hat{\lambda}=\exp\{-(\hat{\beta_0}+\hat{\beta_1}x_1+\ldots +\hat{\beta_m}x_m)/{\hat{\sigma}}\}\] Astfel, functia de supravietuire are expresia: \[\hat{S}(t)=\exp\{-\hat{\lambda}t^{\hat{\alpha}}\}, t\geq 0\]

Cel mai potrivit model…

In general ne intereseaza sa alegem “Cel mai bun model dintre doua modele, exponential si Weibull.

Vom testa urmatoare ipoteza:

\(H_0:\) modelul exponenial este mai potrivit

\(H_a\): modelul Weibull este mai potrivit

Se calculeza statistica \(S=-2[\ln L(\beta_0,\beta_1,\ldots, \beta_m)-\ln L(\beta_0,\beta_1,\ldots, \beta_m,\sigma)]\), unde \(L\) este functia de verosimilitate maxima a modelului exponential si, respectiv Weibul calculata pt parametrii modelului. Decizia se va lua folosind faptul ca \(S\sim \chi^2(1)\).

Aplicatie in R (Exemplu 3.7.4/55)

Un studiu pe 5 ani urmareste evolutia a 20 de pacienti in urma unei interventii chirurgicale de inlocuire de valva la inima. Variabila \(T\) inregistreaza timpul (in ani) dintre operatie si un eveniment. Sunt inregistrate urmatoarele variabile:

\(x_1\)-varsta pacientului la momentul operatiei

\(x_2\) clasa NYHA a pacientului, adica

  -clasa I: pacienti fara limita in activitatea fizica 
  -clasa II: pacienti cu recomandare de activitatea fizica moderata;
  -clasa III: pacienti cu recomandare de  odihna
  -clasa IV: pacienti carora le este recomandat repausul la pat

Sa se asocieze un model parametric setului de date si sa se decida care dintre modelele exponential si Weibull este mai potrivit pentru a descrie setul de date.

Introducem datele problemei

library(survival)
age<-c(56,28,68,54,65,68,59,67,38,52,67,62,52,72,59,54,46,62,73,47)
NYHA<-c(2,3,1,3,4,3,2,3,1,1,3,2,2,2,1,3,2,4,2,1)
Time<-c(2.7,2.8,0.7,2.9,0.2,0.9,3,0.6,3.5,2,0.8,0.1,4,0.3,1.3,4.7,1.8,0.2,0.3,4.6)
Status<-c(0,1,0,1,1,0,0,1,0,1,1,1,1,0,0,0,1,1,0,0)

Definim setul de date valve

valve<-data.frame(cbind(Subject=1:20,age,NYHA,Time,Status))
head(valve)
##   Subject age NYHA Time Status
## 1       1  56    2  2.7      0
## 2       2  28    3  2.8      1
## 3       3  68    1  0.7      0
## 4       4  54    3  2.9      1
## 5       5  65    4  0.2      1
## 6       6  68    3  0.9      0

Construim modelul parametric Weibull (fara a implica variabilele din model)

wei=survreg(Surv(Time,Status)~ 1,dist="w")
wei
## Call:
## survreg(formula = Surv(Time, Status) ~ 1, dist = "w")
## 
## Coefficients:
## (Intercept) 
##    1.351018 
## 
## Scale= 1.088501 
## 
## Loglik(model)= -23.1   Loglik(intercept only)= -23.1
## n= 20

Notam variabila age cu \(Y\) (pentru a simplifica notatia), iar variabilei NYHA ii asociem 3 variabile dummy (daca inregistrarea este in NYHA=i, atunci variabiala \(Y_i\) are valoarea 1 si in rest are valoarea 0).

Y <- age
NYHA
##  [1] 2 3 1 3 4 3 2 3 1 1 3 2 2 2 1 3 2 4 2 1
Y1 <- (NYHA=="1")*1
Y1
##  [1] 0 0 1 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 1
Y2 <- (NYHA=="2")*1
Y2
##  [1] 1 0 0 0 0 0 1 0 0 0 0 1 1 1 0 0 1 0 1 0
Y3 <- (NYHA=="3")*1
Y3
##  [1] 0 1 0 1 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0

Vom construi modele parametrice de regresie (care includ variabile \(Y\) si \(Y_1,Y_2\) si \(Y_3\)) astfel:

wei1=survreg(Surv(Time,Status)~Y+Y1+Y2+Y3,dist="w")
wei1
## Call:
## survreg(formula = Surv(Time, Status) ~ Y + Y1 + Y2 + Y3, dist = "w")
## 
## Coefficients:
## (Intercept)           Y          Y1          Y2          Y3 
## -0.20058220 -0.02217401  3.31305544  2.69709754  2.48016432 
## 
## Scale= 0.6858893 
## 
## Loglik(model)= -17.1   Loglik(intercept only)= -23.1
##  Chisq= 12 on 4 degrees of freedom, p= 0.0173 
## n= 20
exp1=survreg(Surv(Time,Status)~Y+Y1+Y2+Y3,dist="exponential")
exp1
## Call:
## survreg(formula = Surv(Time, Status) ~ Y + Y1 + Y2 + Y3, dist = "exponential")
## 
## Coefficients:
## (Intercept)           Y          Y1          Y2          Y3 
##  -0.5846296  -0.0161341   3.8569922   2.8764140   2.5788070 
## 
## Scale fixed at 1 
## 
## Loglik(model)= -18   Loglik(intercept only)= -23.2
##  Chisq= 10.38 on 4 degrees of freedom, p= 0.0345 
## n= 20

Vom testa ipoteza ca modelul exponential este mai potrivit

teststat=2*(wei1$logl[2]-exp1$logl[2])
pchisq(teststat,df=1,lower.tail=F)
## [1] 0.1889525

Rezultatul obtinut, \(p-value=0.1889\), ne indica faptul ca nu putem respinge ipoteza nula, deci modelul exponential este mai potrivit. Acest model il vom imbunatatii prin analiza variabilelor implicate in model:

summary(exp1)
## 
## Call:
## survreg(formula = Surv(Time, Status) ~ Y + Y1 + Y2 + Y3, dist = "exponential")
##               Value Std. Error     z      p
## (Intercept) -0.5846     2.4089 -0.24 0.8082
## Y           -0.0161     0.0362 -0.45 0.6562
## Y1           3.8570     1.3351  2.89 0.0039
## Y2           2.8764     0.9594  3.00 0.0027
## Y3           2.5788     0.9433  2.73 0.0063
## 
## Scale fixed at 1 
## 
## Exponential distribution
## Loglik(model)= -18   Loglik(intercept only)= -23.2
##  Chisq= 10.38 on 4 degrees of freedom, p= 0.034 
## Number of Newton-Raphson Iterations: 5 
## n= 20

Rezultatele ne sugereaza faptul ca variabila \(Y\)-age nu este semnificativa, deci va fi eliminata din model. Vom constui astfel un nou model (care sa nu contina aceasta variabila).

exp2=survreg(Surv(Time,Status)~Y1+Y2+Y3,dist="exponential")
summary(exp2)
## 
## Call:
## survreg(formula = Surv(Time, Status) ~ Y1 + Y2 + Y3, dist = "exponential")
##              Value Std. Error     z       p
## (Intercept) -1.609      0.707 -2.28 0.02284
## Y1           4.103      1.225  3.35 0.00081
## Y2           3.012      0.913  3.30 0.00097
## Y3           2.765      0.866  3.19 0.00141
## 
## Scale fixed at 1 
## 
## Exponential distribution
## Loglik(model)= -18.1   Loglik(intercept only)= -23.2
##  Chisq= 10.17 on 3 degrees of freedom, p= 0.017 
## Number of Newton-Raphson Iterations: 5 
## n= 20

In acest model toate variabilele sunt semnificative, iar coeficientii modelului sunt: \(\beta_0=-1.609, \beta=4.103, \beta_2=3.012, \beta_3=2,765\)

Deci, \[\hat{\lambda}=\exp\{-(\hat{\beta_0}+\hat{\beta_1}Y_1+ +\hat{\beta_2}Y_2+\hat{\beta_3}Y_3)\}=\exp\{1.609-4.103Y_1-3.012Y_2-2.765Y_3\}\]

Astfel, \(\hat{S}(t)=\exp\{-\hat{\lambda}t\}, t\geq 0\).

Variabila \(T\sim Exp(\lambda)\), atunci timpul mediu de supravietuire este \[E(T)=1/{\hat{\lambda}}=\exp\{\hat{\beta_0}+ \hat{\beta_{1}}Y_1+\hat{\beta_{2}}Y_2+\hat{\beta_{3}}Y_3\}\] Rezultate:

\(\frac{E(T_{Y_2=1})}{E(T_{Y_1=1})}=\frac{e^\hat{\beta_2}}{e^{\hat{\beta_1}}}=e^{\hat{\beta_2}-\hat{\beta_1}}=e^{3.012-4.103}\approx 0.336\);

Deci, timpul mediu de supravietuire pentru subiectii din clasa II (\(Y_2=1\)) este \(33.6\%\) din timpul mediu de supravietuire pentru subiectii din clasa I (\(Y_1=1\))

\(\frac{E(T_{Y_3=1})}{E(T_{Y_1=1})}=\frac{e^\hat{\beta_3}}{e^{\hat{\beta_1}}}=e^{\hat{\beta_3}-\hat{\beta_1}}=e^{2.765-4.103}\approx 0.2624\);

Deci, timpul mediu de supravietuire pentru subiectii din clasa III (\(Y_3=1\)) este \(26.24\%\) din timpul mediu de supravietuire pentru subiectii din clasa I (\(Y_1=1\));

\(\frac{E(T_{Y_1=Y_2=Y_3=0})}{E(T_{Y_1=1})}=\frac{e^0}{e^{\hat{\beta_1}}}=e^{0-4.103}\approx 0.165\);

Deci, timpul mediu de supravietuire pentru subiectii din clasa IV (\(Y_1=Y_2=Y_3=0\)) este \(1.65\%\) din timpul mediu de supravietuire pentru subiectii din clasa I (\(Y_1=1\));

Analog si pentru alte comparatii…

*Olga Korosteleva, Clinical Statistics: Introducing Clinical Trials, Survival Analysis, and Longitudinal Data Analysis.