In acest material ne propunem sa studiem urmatoarele:
-metode de estimare a parametrilor modelului Cox, adica a coeficientilor \(\beta_i\) din expresia functiei de hazard a modelului Cox;
-verficarea conditiilor PHA (de proportionalitate a functiei de hazard);
-modalitati de scalare a variabilelor in modelul Cox pentru a determina cea mai buna forma sub care acestea apar in modelul Cox.
Pentru a exemplifica aceste concepte, vom folosi setul de date whas500. Acest set de date contine 500 de observatii si 22 de variabile, studiindu-se timpul de supravietuire pentru persoane cu infarct miocardic acut (mai multe info la https://web.archive.org/web/20170517071528/http://www.umass.edu/statdata/statdata/data/whas500.txt)
Vom apela setul de whas (salvat pe calculator)
library(survival)
library(splines)w <- read.csv(file = 'whas500.csv')
head(w)## id age gender hr sysbp diasbp bmi cvd afb sho chf av3 miord mitype year
## 1 1 83 0 89 152 78 25.54051 1 1 0 0 0 1 0 1
## 2 2 49 0 84 120 60 24.02398 1 0 0 0 0 0 1 1
## 3 3 70 1 83 147 88 22.14290 0 0 0 0 0 0 1 1
## 4 4 70 0 65 123 76 26.63187 1 0 0 1 0 0 1 1
## 5 5 70 0 63 135 85 24.41255 1 0 0 0 0 0 1 1
## 6 6 70 0 76 83 54 23.24236 1 0 0 0 1 0 0 1
## admitdate disdate fdate los dstat lenfol fstat
## 1 1/13/1997 1/18/1997 12/31/2002 5 0 2178 0
## 2 1/19/1997 1/24/1997 12/31/2002 5 0 2172 0
## 3 1/1/1997 1/6/1997 12/31/2002 5 0 2190 0
## 4 2/17/1997 2/27/1997 12/11/1997 10 0 297 1
## 5 3/1/1997 3/7/1997 12/31/2002 6 0 2131 0
## 6 3/11/1997 3/12/1997 3/12/1997 1 1 1 1Estimarea parametrilor modelului Cox se face folosind metoda Breslow sau metoda Efron ( mai multe detalii se pot gasi la pagina 11 a materialul https://cv.upt.ro/pluginfile.php/322403/mod_resource/content/1/notes3.pdf)
Vom construi un model Cox care sa aiba doar variabila bmi (body max index). Mentionam ca variabila lenfol=inregistreaza timpul pana la evenimentul urmarit, iar
fstat=inregistreaza tipul observatiei (cenzurat sau nu).
Primul model Cox cu o singura variabila, pe care il propunem, este:
phm.bmi_b = coxph(Surv(w$lenfol,w$fstat)~w$bmi,method = "breslow")
phm.bmi_b## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi, method = "breslow")
##
## coef exp(coef) se(coef) z p
## w$bmi -0.09827 0.90641 0.01474 -6.665 2.64e-11
##
## Likelihood ratio test=48.16 on 1 df, p=3.92e-12
## n= 500, number of events= 215phm.bmi_e = coxph(Surv(w$lenfol,w$fstat)~w$bmi,method = "efron")
phm.bmi_e## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi, method = "efron")
##
## coef exp(coef) se(coef) z p
## w$bmi -0.09847 0.90622 0.01475 -6.676 2.46e-11
##
## Likelihood ratio test=48.33 on 1 df, p=3.602e-12
## n= 500, number of events= 215Se observa ca cele doua metode dau rezultate asemanatoare. In R este implementata by defaut metoda Efron. Din outout, observam ca variabila bmi este semnificativa.
phm.bmi = coxph(Surv(w$lenfol,w$fstat)~w$bmi)
phm.bmi## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi)
##
## coef exp(coef) se(coef) z p
## w$bmi -0.09847 0.90622 0.01475 -6.676 2.46e-11
##
## Likelihood ratio test=48.33 on 1 df, p=3.602e-12
## n= 500, number of events= 215Vom construi un alt model Cox ce contine doar variabila binara gender (0 = Masculin, 1 = Feminin)
phm.gender = coxph(Surv(w$lenfol,w$fstat)~w$gender,method = "breslow")
phm.gender## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$gender, method = "breslow")
##
## coef exp(coef) se(coef) z p
## w$gender 0.3813 1.4641 0.1376 2.771 0.00559
##
## Likelihood ratio test=7.59 on 1 df, p=0.005876
## n= 500, number of events= 215Din output, se observa ca variabila gender este semnificativa.
Vom construi un alt model Cox, cu variabila binara cvd, ce inregistreaza istoricul bolilor cardiovasculare ale pacientilor (0 = Nu, 1 = Da ).
phm.cvd = coxph(Surv(w$lenfol,w$fstat)~w$cvd,method = "breslow")
phm.cvd## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$cvd, method = "breslow")
##
## coef exp(coef) se(coef) z p
## w$cvd 0.2828 1.3269 0.1680 1.684 0.0922
##
## Likelihood ratio test=2.99 on 1 df, p=0.0838
## n= 500, number of events= 215Din nou se observa ca variabila cvd nu este semnificativa, avand un p-value=0.092>0.05.
Ne propunem sa construim un model Cox cu variabilele bmi, gender si age (varsta).
m1=coxph(Surv(w$lenfol,w$fstat)~w$bmi + w$age + w$gender)
m1## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi + w$age + w$gender)
##
## coef exp(coef) se(coef) z p
## w$bmi -0.042084 0.958789 0.015427 -2.728 0.00637
## w$age 0.060791 1.062677 0.006503 9.348 < 2e-16
## w$gender -0.093065 0.911135 0.141085 -0.660 0.50949
##
## Likelihood ratio test=150 on 3 df, p=< 2.2e-16
## n= 500, number of events= 215In acest model, variabila gender devine nesemnificativa in acest model.
Un posibil model cu interactiune intre variabile este
mi1=coxph(Surv(w$lenfol,w$fstat)~w$bmi*w$gender)
mi1## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi * w$gender)
##
## coef exp(coef) se(coef) z p
## w$bmi -0.11716 0.88944 0.02346 -4.993 5.93e-07
## w$gender -0.80152 0.44865 0.76538 -1.047 0.295
## w$bmi:w$gender 0.04094 1.04179 0.03021 1.355 0.175
##
## Likelihood ratio test=52.56 on 3 df, p=2.28e-11
## n= 500, number of events= 215sau
mi2=coxph(Surv(w$lenfol,w$fstat)~w$bmi*w$age)
mi2## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi * w$age)
##
## coef exp(coef) se(coef) z p
## w$bmi -0.0308178 0.9696523 0.0916772 -0.336 0.7368
## w$age 0.0635256 1.0655867 0.0305271 2.081 0.0374
## w$bmi:w$age -0.0001345 0.9998655 0.0011717 -0.115 0.9086
##
## Likelihood ratio test=149.6 on 3 df, p=< 2.2e-16
## n= 500, number of events= 215Sau, model fara interactiune cu variabilele age si bmi:
mi=coxph(Surv(w$lenfol,w$fstat)~w$bmi+w$age)
mi## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi + w$age)
##
## coef exp(coef) se(coef) z p
## w$bmi -0.041194 0.959642 0.015288 -2.695 0.00705
## w$age 0.060103 1.061946 0.006433 9.344 < 2e-16
##
## Likelihood ratio test=149.6 on 2 df, p=< 2.2e-16
## n= 500, number of events= 215In modelul mi variabile age si bmi sunt semnificative, avand coeficientii \(\beta_{age}=0.06, \beta_{bmi}=-0.04\).
Vom construi modelul Cox
\[h(t,age, bmi)=h_0(t)exp(0.06*age-0.04*bmi)\]
Modelarea Cox se bazeaza pe presupunerea ca functiile de hazard sunt proportionale. Mai exact, se presupune ca functiile de hazard pentru oricare doua persoane la orice moment de timp sunt proportionale (daca o persoana are riscul de deces la un moment de timp de k ori fata de alt individ, atunci la un moment de timp ulterior, acest risc este la fel). Acest lucru inseamna ca variabile nu depind de timp.
Presupunerea PHA trebuie testata.
Exista urmatoarele metodele pentru a verifica PHA:
-grafic, reprezentand pentru fiecare grup functia log-log (S(t,X)); daca cele doua grafice sunt paralele, atunci PHA este verificata;
-folosind testarea residurilor(a se vedea https://cv.upt.ro/pluginfile.php/330195/mod_resource/content/3/notes4.pdf)
-extinderea modelului prin includerea dependentei de timp a variabilelor modelului.
In R se poate studia problema verificarii PHA folosind comanda cox.zph pentru fiecare model considerat si alegerea modelului in care outputul furnizeaza p-value cu valori mari (datorita ipotezelor de testare folosite).
library(splines)phm.bmi1 = coxph(Surv(w$lenfol,w$fstat)~w$bmi)
cox.zph(phm.bmi1,global=F)## chisq df p
## w$bmi 2.14 1 0.14phm.bmi_age_sysbp = coxph(Surv(w$lenfol,w$fstat)~w$bmi+w$age+w$sysbp)
cox.zph(phm.bmi_age_sysbp, global=F)## chisq df p
## w$bmi 0.987 1 0.321
## w$age 0.182 1 0.669
## w$sysbp 6.294 1 0.012Obiectivul urmarit este acela de a studia modul in care o variabila este inclusa in model. Mai exact, ne propunem sa determinam cel mai potrivit model Cox care descrie setul de date, implicand o variabila \(X\) in model in diverse scalari, cum ar fi, de exemplu, \(X^2, \log(X), \sqrt{X}\) etc.
Dintre modelele propuse vom alege modelul cu cel mai mic AIC(Akaike Information Criterion).
Iata mai jos, pentru setul de date Whas500 cum scalam variabila bmi si cum vom determina, in R, AIC-ul fiecarui model.
Primul model Cox este cel in care variabila bmi este nescalata. Al doilea model propus este cel in care variabila bmi apare la puterile 1 si 2. In al treilea model, se propune log(bmi), iar in modelul 4 apare radical(bmi) . In modelul 5 variabila bmi apare scalata in mai multe forme, anume \(bmi^2, 1/bmi^2, 1/bmi^3\).
phm.bmi1 = coxph(Surv(w$lenfol,w$fstat)~w$bmi)
phm.bmi2=coxph(Surv(w$lenfol,w$fstat)~w$bmi+I(w$bmi^2))
phm.bmi3=coxph(Surv(w$lenfol,w$fstat)~log(w$bmi^2))
phm.bmi4= coxph(Surv(w$lenfol,w$fstat)~sqrt(w$bmi^2))
phm.bmi5 = coxph(Surv(w$lenfol,w$fstat)~w$bmi^2+I(w$bmi^2)+I(w$bmi^3))
aic1=extractAIC(phm.bmi1)[2]
aic1## [1] 2408.311aic2=extractAIC(phm.bmi2)[2]
aic2## [1] 2403.903aic3=extractAIC(phm.bmi3)[2]aic3## [1] 2404.616aic4=extractAIC(phm.bmi4)[2]
aic4## [1] 2408.311aic5=extractAIC(phm.bmi5)[2]
aic5## [1] 2404.57Se observa ca modelul 2 are cel mai mic AIC.
phm.bmi2=coxph(Surv(w$lenfol,w$fstat)~w$bmi+I(w$bmi^2))
phm.bmi2## Call:
## coxph(formula = Surv(w$lenfol, w$fstat) ~ w$bmi + I(w$bmi^2))
##
## coef exp(coef) se(coef) z p
## w$bmi -0.328297 0.720149 0.086538 -3.794 0.000148
## I(w$bmi^2) 0.004391 1.004400 0.001617 2.715 0.006627
##
## Likelihood ratio test=54.74 on 2 df, p=1.299e-12
## n= 500, number of events= 215Deci, \(h(t,bmi)=h_0(t)exp(-0.3*bmi+0.004* (1/bmi^2))\)
Concluzie:
Exista nenumarate moduri de a scala variabilele implicate intr-un model. Practica sugereaza se se incerce initial cu scalarea patratica \(X^2\) sau cubica \(X^3\). Apoi, sa se incerce functii de baza ale variabilei \(X\), cum ar fi \(\log(X)\) sau \(\sqrt{X}\).
In cazurile in care exista un numar mare de variabile, este indicat ca sa se scalaze cate o variabila pe rand si sa se observe daca, modificand variabila, se obtine un model mai bun. Daca nu, sa se ramana la modelul liniar (al functiei \(\log (h(T,X))\)).