In acest document vom prezenta rezolvarea probemei 3, pag 125, cartea Survival Analysis: A self learning book, Kleibaum.

Setul de date Anderson descrie inregistrarile legate de timpul pana la remisie a 42 de pacienti bolnavi de leucemie. Acest setde date il gasiti pe CV-UPT, la cursul 3. Mai multe informatii despre acest set de date se gasesc disponibile la

https://cardiomoon.github.io/autoReg/reference/anderson.html

Rezolvare in R

Vom apela setul de date anderson (salvat pe calculator)

library(survival)
q<- read.table("anderson2.txt",header=TRUE)
head(q)
##   st status sex log_WBC Rx
## 1 35      0   1    1.45  0
## 2 34      0   1    1.47  0
## 3 32      0   1    2.20  0
## 4 32      0   1    2.53  0
## 5 25      0   1    1.78  0
## 6 23      1   1    2.57  0
length(q$st)
## [1] 42

Modele Cox construite pentru setul de date Anderson

Se considera diferite modele Cox asociate setului de date Anderson.

Vom studia pentru inceput daca, in cazul modelului \(m_1,\) termenii de interactiune sunt semnificativi.

m1=coxph(Surv(q$st,q$status)~q$Rx+q$sex+q$log_WBC+q$Rx*q$sex+ q$Rx*q$log_WBC)
m1
## Call:
## coxph(formula = Surv(q$st, q$status) ~ q$Rx + q$sex + q$log_WBC + 
##     q$Rx * q$sex + q$Rx * q$log_WBC)
## 
##                    coef exp(coef) se(coef)      z        p
## q$Rx            0.76667   2.15258  1.84081  0.416   0.6771
## q$sex          -1.04804   0.35062  0.76099 -1.377   0.1684
## q$log_WBC       1.75101   5.76041  0.44420  3.942 8.08e-05
## q$Rx:q$sex      2.06385   7.87624  0.91614  2.253   0.0243
## q$Rx:q$log_WBC -0.08956   0.91433  0.53693 -0.167   0.8675
## 
## Likelihood ratio test=52.75  on 5 df, p=3.787e-10
## n= 42, number of events= 30

Se observa ca pentru variabila q\(Rx*q\)sex valoarea p-value =0.02<0.05, deci aceasta variabila este semnificativa, iar variabila q\(Rx*q\)\(log_WBC\) nu este semnificativa.

Se considera un nou model \(m_2\). Care dintre modelele \(m_1\) si \(m_2\) este de preferat?

m2=coxph(Surv(q$st,q$status)~q$Rx+q$sex+q$log_WBC+q$Rx*q$sex)
m2
## Call:
## coxph(formula = Surv(q$st, q$status) ~ q$Rx + q$sex + q$log_WBC + 
##     q$Rx * q$sex)
## 
##               coef exp(coef) se(coef)      z        p
## q$Rx        0.4752    1.6084   0.5607  0.847   0.3967
## q$sex      -1.0838    0.3383   0.7310 -1.483   0.1382
## q$log_WBC   1.7028    5.4895   0.3380  5.037 4.72e-07
## q$Rx:q$sex  2.1032    8.1925   0.8858  2.374   0.0176
## 
## Likelihood ratio test=52.72  on 4 df, p=9.747e-11
## n= 42, number of events= 30

Pentru a determina care dintre cele doua model este de preferat, vom calcula AIC fiecarui model:

library(MASS)
aic1=extractAIC(m1)[2]
aic1
## [1] 143.6197
aic2=extractAIC(m2)[2]
aic2
## [1] 141.6474

Se observa ca modelul \(m_2\) are AIC mai mic, deci este de preferat.

Folosindu-se modelul \(m_2\) sa se determine efectul tratamentului, i.e. al variabilei Rx (1=placebo, 0=tratement), variabilele celelalte fiind ajustate, adica considerate la o valoare medie.

Avem, \(HR=\frac{h(t,Rx=1)}{h(t,Rx=0)}\).

Observam ca acest raport de hazar are doua valori diferite, conform celor doua categorii ale variabilei sex.

Mai exact: \(HR(sex=female)=\exp(0.47+2.1)=\exp(2.57)\), iar \(HR(sex=male)=\exp(0.47)\).

exp(2.57)
## [1] 13.06582
exp(0.47)
## [1] 1.599994

Se considera modelul \(m_3\). Sa se studieze daca acest model este mai potrivit fata de modelul \(m2\) pentru a descrie setul de date.

m3=coxph(Surv(q$st,q$status)~q$Rx+q$sex+q$Rx*q$sex)
m3
## Call:
## coxph(formula = Surv(q$st, q$status) ~ q$Rx + q$sex + q$Rx * 
##     q$sex)
## 
##               coef exp(coef) se(coef)      z      p
## q$Rx        0.6243    1.8670   0.5411  1.154 0.2486
## q$sex      -1.1080    0.3302   0.7058 -1.570 0.1164
## q$Rx:q$sex  1.9570    7.0778   0.8145  2.403 0.0163
## 
## Likelihood ratio test=22.47  on 3 df, p=5.199e-05
## n= 42, number of events= 30
aic2=extractAIC(m2)[2]
aic2
## [1] 141.6474
aic3=extractAIC(m3)[2]
aic3
## [1] 169.8952

Se observa ca modelul \(m_3\) are un AIC mai mare, deci modelul se inrautateste.

Se considera modelul \(m_4\) de mai jos. Sa se studieze care din modelele \(m_1,m_2,m_2\) si \(m_4\) este mai potrivit pentru a descrie setul de date.

m4=coxph(Surv(q$st,q$status)~q$Rx+q$sex+q$log_WBC)
m4
## Call:
## coxph(formula = Surv(q$st, q$status) ~ q$Rx + q$sex + q$log_WBC)
## 
##             coef exp(coef) se(coef)     z        p
## q$Rx      1.5036    4.4978   0.4615 3.258  0.00112
## q$sex     0.3147    1.3698   0.4545 0.692  0.48872
## q$log_WBC 1.6819    5.3760   0.3366 4.997 5.82e-07
## 
## Likelihood ratio test=47.19  on 3 df, p=3.171e-10
## n= 42, number of events= 30
aic4=extractAIC(m4)[2]
aic4
## [1] 145.181

Conform criteriului celui mai mic AIC, observam ca modelul \(m_2\) este de preferat.