Obiective:

-introducerea unor variatelor continue in studiul timpului de supravietuire;

-analiza efectului acestor variabile asupra timpului de supravietuire;

Modelul Cox

A fost introdus de catre David Cox, prin lucrarea “Regression models and life tables, (1972, J.R. Stat. Soc. B, 34:187–220)”, citata de mai mult de 55000 ori.

Instrumentul matematic folosit este functia hazard \(h(t,X)\); vom obtine rezultate de forma: riscul de deces al categoriei A este de k-ori fata de riscul de deces al categoria B;

Modelul Cox presupune urmatoarea forma pentru functia hazard \[h(t,\textbf{X})=h_0(t,\alpha) e^{\beta^T\textbf{X}}, \] unde \(\beta^T=(\beta_1,\ldots, \beta_n)\) este vectorul coeficientilor, iar \(\textbf{X}=(X_1,X_2,\ldots, X_n)^T\) este vectorul variabilelor inregistrate in studiu.

Observatii:

-termenul \(h_0(t,\alpha)\) depinde doar de timp, nu si de covariate; din acest motiv se numeste functie hazard de baza( toate variabilele \(X_i\) sunt =0).

-termenul \(e^{\beta^T\textbf{X}}=e^{\beta_1X_1+\beta_2X_2+\ldots+\beta_n X_n}\) depinde de covariate, nu de timp.

-se pot obtine informatii despre timpul de supravietuire fara a cunoste forma exacta a lui \(h_0(t,\alpha);\) din acest motiv acest model se numeste semi-parametric.

Cum putem obtine informatii din modelul Cox?

Prin determinarea unui raport: HR-proportional hazard ratio.

De exemplu, daca facem raportul dintre functia hazard, cazul in care variabila \(X_1=a\), iar toate celelalte sunt fixate la o valoare medie, si dintre functia hazard, cazul in care variabila \(X_1=b\), iar toate celelalte sunt fixate la o valoare medie, observam ca se obtine valoarea \(HR=\frac{h(t,X_1=a)}{h(t,X_1=b)}=e^{\beta_1(a-b)}\).

Deci, risul de aparitie al evenimentului urmarit pentru subiectul cu covariata \(X_1=a\) este de \(e^{\beta_1(a-b)}\) fata de cel al subiectului cu covariata \(X_1=b\) (la un moment de timp fixat).

Pentru a decide daca riscul de apartie al evenimentului creste sau scade intre cele doua clase, vom face bbservatie ca:

-daca \(HR>1\), atunci riscul de apartie al evenimentului creste;

-daca \(HR<1\), deci riscul de apartie al evenimentului scade.

Observatie: Modelul Cox este un model, poate fi gresit, fiind construit pe urmatoarele presupuneri (PHA=proportional hazard assumptions):

-nu exista interactiune intre timp si variabilele din studiu (nu depind de timp);

-functia \(\log h(t,X)\) este liniara in covariate.

In continuare vom studia modul in care diverse tipuri de variabile (continue, de categorie) intervin in modelul Cox si cum putem face interpretari legate de cresterea sau scaderea riscului de aparitie al evenimentului prin analiza raportului HR, pentu fiecare caz.

Vom exemplifica teoria folosind setul de date uis.txt. Acest set de date contine 628 de inregistrari si 12 variabile si urmareste timpul pana cand consumatorii de droguri vor reincepe consumul lor. Mai multe informatii la https://rdrr.io/github/dave-harrington/eventtimedata/man/uis.html

Vom apela setul de uis (salvat pe calculator)

library(survival)

q<- read.table("uis.txt",header=TRUE, na.strings=".")
head(q)
##   ID AGE BECKTOTA HERCOC IVHX NDRUGTX RACE TREAT SITE LOS TIME CENSOR
## 1  1  39     9.00      4    3       1    0     1    0 123  188      1
## 2  2  33    34.00      4    2       8    0     1    0  25   26      1
## 3  3  33    10.00      2    3       3    0     1    0   7  207      1
## 4  4  32    20.00      4    3       1    0     0    0  66  144      1
## 5  5  24     5.00      2    1       5    1     1    0 173  551      0
## 6  6  30    32.55      3    3       1    0     1    0  16   32      1

Ne propunem sa avem diverse modele de tip Cox, introducand pe rand variabile (continue sau de tip categorie) si studiind valoarea HR.

Modelul Cox cu o singura variabila continua

Abordare teoretica

Se considera setul de date in care:

-variabila \(X_1\) este continua;

-parametrul ei este \(\beta_1\);

-functia hazard \(h(t,X)=h_0(t,\alpha)e^{\beta_1 X_1}\);

-avem \(HR=\frac{h(t,X_1=a)}{h(t,X_1=b)}=e^{\beta_1(a-b)}\);

-de exemplu, daca \(X_1\)-variabila varsta, iar \(a=b+1\), atunci \(HR=e^{\beta_1}\) ne da riscul de a se produce evenimentul cand la varsta se adauga un an;

Abordare in R

Vom analiza influenta variabile AGE-varsta la momentul inrolarii in studiu- asupra riscului de apartiei a evenimentului, prin construirea unui model Cox, ce contine doar aceasta variabila.

m1=coxph(formula=Surv(q$TIME,q$CENSOR)~q$AGE)
m1
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$AGE)
## 
##            coef exp(coef)  se(coef)      z      p
## q$AGE -0.012878  0.987204  0.007188 -1.792 0.0732
## 
## Likelihood ratio test=3.25  on 1 df, p=0.07155
## n= 623, number of events= 504 
##    (5 observations deleted due to missingness)

Am obtinut urmatoare forma a functiei de hazard \[h(t, AGE)=h_0(t)exp(-0.013*AGE)\] Se observa ca in acest caz \(\exp(\beta_1)=0.98\approx 1\), deci nu putem spune ca la o diferenta de un an, riscul ca o persoana sa reincepa sa consume droguri creste sau scade, fiind o valoare apropriata de 1. Din output-ul de la R se observa ca avem variabila AGE are un p-value=0.07>0.05, deci putem considera variabila age nesemnificativa.

In continuare ne propunem sa construim alt model Cox, in care sa introducem doua variabile continue.

Modelul Cox cu doua variabile continue

Abordare teoretica

Consideram variabilele continue \(X_1, X_2\), avand in modelul Cox parametrii sunt \(\beta_1,\beta_2\) si \(\beta_{12}\)-daca exista interactiune intre variabile (se influenteaza reciproc).

In cazul modelului fara interactiune intre \(X_1\) si \(X_2\), functia hazard este \[h(t,X)=h_0(t,\alpha)e^{\beta_1 X_1+\beta_2X_2}\]

In acest caz, \(HR=\frac{h(t,X_1=a,X_2-fixat)}{h(t,X_1=b, X_2-fixat)}=e^{\beta_1(a-b)}\)

De exemplu, daca \(a=b+1\), atunci \(HR=e^{\beta_1}\) ne da riscul de a se produce evenimentul cand variabilei \(X_1\) i se adauga o unitate de masura;

In cazul in care se presupune ca exista interactiune intre variabilele \(X_1\) si \(X_2\), atunci functia hazard se scrie \[h(t,X)=h_0(t,\alpha)e^{\beta_1 X_1+\beta_2X_2+\beta_{12}X_1X_2}\]

In acest caz interpretarile sunt mai greu de facut.

Abordare in R

Folosind setul de date uis.txt, vom construi modele Cox folosind variabilele: Becktota=scorul Beck de depresie (cu valori in intervalul [0,54]) si NDRUG=numarul de tratamente antidrog anterioare (cu valori in intervalul [0,40]).

Construim modelul fara interactiune

m2=coxph(formula=Surv(q$TIME,q$CENSOR)~q$BECKTOTA+q$NDRUGTX)
m2
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$BECKTOTA + q$NDRUGTX)
## 
##                coef exp(coef) se(coef)     z        p
## q$BECKTOTA 0.009553  1.009599 0.004818 1.983   0.0474
## q$NDRUGTX  0.030238  1.030699 0.007649 3.953 7.72e-05
## 
## Likelihood ratio test=17.58  on 2 df, p=0.0001524
## n= 580, number of events= 468 
##    (48 observations deleted due to missingness)

Din ouutput-ul, se obtine ca \[h(t,Beck, Ndrug)=h_0(t)exp(0.01*Beck+0.03*Ndrug) \]

Variabila scorul BECK are urmatoarea interpretare:

valori 0-9 nu exista depresie

valori 10-18 depresie moderata

valori 19-29 depresie moderat-severa

valori 30-54 depresie severa

In functie de ce dorim sa analizam, putem face diverse comparatii: de exemplu, putem determina riscul de aparitie al evenimentului urmarit (reinceperea consumului de droguri) de o persoana cu scor Beck 40=depresie-severa fata de o persoana cu scor Beck 15=depresie medie, Obtinem \(HR=exp(0.01*(40-15))=1.27>1\). Deci, riscul creste de 1.27 ori. Se observa ca variabila scor Beck are un p-value apropiat de 0.05.

In cazul variabilei NDRUG-numar de tratemente antidrog anterioare, putem obtine, de exemplu, ca o persoana care a avut 5 tratamente abtidrog anterioare are \(e^{0.03*5}\approx 1.16\) sanse sa reapuce de consum, in raport cu o persoana care nu a avut tratamente antidrog. Intr-adevar, \[HR=\frac{e^{0.03*5}}{e^{0.03*0}}=e^{0.03*5}=1.16>1\] Sau, putem compara o persoana cu 20 tratemente anterioare cu o persoana cu ) tratamenet. Avem

\[HR=\frac{e^{0.03*20}}{e^{0.03*0}}=e^{0.03*20}=1.8>1\] Se observa ca variabila NDRG este semnificativa, avand un p-value asociat foarte mic.

Construim modelul cu interactiune

m2i=coxph(formula=Surv(q$TIME,q$CENSOR)~q$BECKTOTA+q$NDRUGTX+q$BECKTOTA*q$NDRUGTX)
m2i
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$BECKTOTA + q$NDRUGTX + 
##     q$BECKTOTA * q$NDRUGTX)
## 
##                           coef exp(coef)  se(coef)     z     p
## q$BECKTOTA           0.0042471 1.0042561 0.0063580 0.668 0.504
## q$NDRUGTX            0.0078984 1.0079297 0.0192870 0.410 0.682
## q$BECKTOTA:q$NDRUGTX 0.0012341 1.0012349 0.0009519 1.296 0.195
## 
## Likelihood ratio test=19.24  on 3 df, p=0.000244
## n= 580, number of events= 468 
##    (48 observations deleted due to missingness)

In acest caz, functia hazard este

\[h(T,X)=h_0(t)exp(0.004*Beck+0.0079*Ndrug+0.0012*Beck*Ndrug)\] dar nici una dintre variabilele modelului nu mai sunt semnificative.

Modelul Cox cu o singura variabila binara

Abordare teoretica

Consideram variabila binara \(X_1\in\{0,1\}\) si notam cu \(\beta_1\) parametrul ei din modelul Cox asociat.

In acest caz functia hazard \(h(t,X_1)=h_0(t,\alpha)e^{\beta_1 X_1}\).

Mai exact, \(h(t,X_1=0)=h_0(t)\) si \(h(t,X_1=1)=h_0(t)e^{\beta_1}\)

Raportul de hazard \(HR=\frac{h(t,X_1=1)}{h(t,X_1=0)}=e^{\beta_1}\)

Abordare in R

Vom construi un model Cox prin care sa studiem cum influenteaza variabila binara RACE riscul de aparitie al evenimentului. Variabila RACE are doua nivele, white=0, alti=1.

m3=coxph(formula=Surv(q$TIME,q$CENSOR)~q$RACE)
m3
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$RACE)
## 
##           coef exp(coef) se(coef)      z      p
## q$RACE -0.2846    0.7523   0.1059 -2.687 0.0072
## 
## Likelihood ratio test=7.59  on 1 df, p=0.005881
## n= 622, number of events= 504 
##    (6 observations deleted due to missingness)

In acest caz functi hazard este \[h(t,Race)=h_0(t)e^{-0.29*Race}\]

Se observa ca raportul \(HR=\frac{h(t,X_1=1)}{h(t,X_1=0)}=e^{\beta_1}=0.7\) este mai mic decat 1, deci riscul de reconsum de droguri al categoriei 1(=altii) fata de cei din categoria 0(=white) este de 0.7 ori, deci scade.

In plus, se obesrva ca p-value este foarte mic, deci variabila Race este semnificativa.

Modelul Cox cu o singura variabila cu K nivele

Abordare teoretica

Consideram o variabila \(X_1\) de tip categorie, avand urmatoarele \(K\) valori \(\{c_0,c_1,\ldots c_{K-1}\}\). Ne propunem sa cream \(K-1\) variabile dummy, notate cu \(C_k,k=1,\ldots K-1\)

Functia hazard \(h(t,X_1)=h_0(t,\alpha)e^{\beta_1 C_1+\ldots \beta_{K-1}C_{K-1}}\), adica:

\(h(t,c_0)=h_0(t)\)

\(h(t,c_1)=h_0(t)e^{\beta_1}\)

\(\ldots\)

\(h(t,c_{K-1})=h_0(t)e^{\beta_{K-1}}\)

Obtinem:

\(HR=\frac{h(t,c_k)}{h(t,c_0)}=e^{\beta_k},\quad k\neq 0\)

\(HR=\frac{h(t,c_k)}{h(t,c_j)}=e^{\beta_k-\beta_j},\quad k,j\neq 0\)

Abordare in R

In setul de date uis.txt variabila HERCOC este o variabila cu 4 nivele:

\(c_0\)-consum de heroina si cocaina

\(c_1\)-consum doar de heroina

\(c_2\)-consum de cocaina

\(c_3\)- consum de alte droguri (nu cocaina, nu heroina)

Pentru a analiza efectul variabilei Hercoc asupra riscului reinceperii consumului de droguri, se vor transforma aceste variabile in 3 variabile dummy, in care categoria \(c_0\) este de referinta:

C1<-(q$HERCOC=="2")*1
C2<-(q$HERCOC=="3")*1
C3<-(q$HERCOC=="4")*1
m4=coxph(formula=Surv(q$TIME,q$CENSOR)~C1+C2+C3)
m4
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ C1 + C2 + C3)
## 
##        coef exp(coef) se(coef)      z      p
## C1  0.07759   1.08068  0.14450  0.537 0.5913
## C2 -0.25489   0.77501  0.13508 -1.887 0.0592
## C3 -0.16220   0.85027  0.13017 -1.246 0.2127
## 
## Likelihood ratio test=7.79  on 3 df, p=0.0505
## n= 610, number of events= 493 
##    (18 observations deleted due to missingness)

Avem:

\(h(t,c_0)=h_0(t)\)

\(h(t,c_1)=1.08*h_0(t)\)

\(h(t,c_2)=0.78*h_0(t)\)

\(h(t,c_3)=0.85*h_0(t)\)

Se observa ca variabilele \(C_1,C_2,C_3\) sunt nesemnificative, deci se considera ca variabila Hercoc nu influenteaza reconsumul de droguri.

Modelul Cox cu o variabila binara si o variabila continua

Consideram variabila \(X_1\in\{0,1\}\) si \(X_2\in\mathbb{R}\).

Fie parametrii \(\beta_1,\beta_2\) pentru modelul fara interactiune si \(\beta_1,\beta_2,\beta_{12}\) pentru modelul cu interactiune.

Functia hazard pentru modelul fara interactiune

\(h(t,X_1,X_2)=h_0(t,\alpha)e^{\beta_1 X_1+ \beta_{2}X_2}\), adica \(h(t,X_1=0,X_2)=h_0(t)e^{\beta_2X_2}\)

\(h(t,X_1=1,X_2)=h_0(t)e^{\beta_1+\beta_2X_2}.\)

Raportul de hazard este

\(HR=\frac{h(t,X_1=1)}{h(t,X_1=0)}=e^{\beta_1}.\)

Functia hazard pentru modelul cu interactiune \[h(t,X_1,X_2)=h_0(t,\alpha)e^{\beta_1 X_1+ \beta_{2}X_2+\beta_{12}X_1X_2},\]

\(h(t,X_1=0,X_2)=h_0(t)e^{\beta_2X_2}\)

\(h(t,X_1=1,X_2)=h_0(t)e^{\beta_1+(\beta_2+\beta_{12})X_2}.\)

Avem \[HR=\frac{h(t,X_1=0,X_2=a+1)}{h(t,X_1=0,X_2=1)}=e^{\beta_2},\] \[HR=\frac{h(t,X_1=1,X_2=a+1)}{h(t,X_1=1,X_2=1)}=e^{\beta_2+\beta_{12}}. \]

Abordarea in R

Ne propunem sa studiem cum este influentata reinceperea consumului de droguri de variabilele Race si Ndrug-numarul de tratamente antidrog anterioare.

Modelul fara interactiune ce implica aceste doua variabile este

m5=coxph(formula=Surv(q$TIME,q$CENSOR)~q$RACE+q$NDRUGTX)
m5
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$RACE + q$NDRUGTX)
## 
##                coef exp(coef)  se(coef)      z        p
## q$RACE    -0.256819  0.773508  0.107428 -2.391 0.016820
## q$NDRUGTX  0.026618  1.026976  0.007616  3.495 0.000474
## 
## Likelihood ratio test=18.28  on 2 df, p=0.0001074
## n= 607, number of events= 493 
##    (21 observations deleted due to missingness)

iar functia hazard asociata este

\(h(t,Race, Ndrug)=h_0(t)exp(-0.26*Race+0.027*Ndrug)\)

Observam ca modelul cu interactiune aduce variabila \(RACE*NDRUGTX\) nesemnificativa, cu p-value=0.95.

m5i=coxph(formula=Surv(q$TIME,q$CENSOR)~q$RACE+q$NDRUGTX+q$RACE*q$NDRUGTX)
m5i
## Call:
## coxph(formula = Surv(q$TIME, q$CENSOR) ~ q$RACE + q$NDRUGTX + 
##     q$RACE * q$NDRUGTX)
## 
##                       coef exp(coef)  se(coef)      z        p
## q$RACE           -0.251486  0.777644  0.145241 -1.732 0.083361
## q$NDRUGTX         0.026761  1.027122  0.008044  3.327 0.000879
## q$RACE:q$NDRUGTX -0.001351  0.998650  0.024813 -0.054 0.956576
## 
## Likelihood ratio test=18.28  on 3 df, p=0.0003848
## n= 607, number of events= 493 
##    (21 observations deleted due to missingness)

Astfel, se va prefera modelul ce contine variabilele Race si Ndrug, fara interactiune.