In acest curs ne propunem sa imbunatatim modelul Cox, printr-o selectie a variabilelor modelului. Vom descrie aceasta selectie folosind setul de date Veteran, unul din seturile de date din recuzita R.

Iata cateva informatii legate de setul de date Veteran (Veterans’ Administration Lung Cancer stud.

Descriere: studiul clinic al unui tratament de cancer de plamani ( Randomised trial of two treatment regimens for lung cancer. This is a standard survival analysis data set Usage veteran data(cancer, package=“survival”) Variabile: trt: 1=standard 2=test; celltype: 1=squamous, 2=smallcell, 3=adeno, 4=large; time: survival time; status: censoring status; karno: Karnofsky performance score (100=good); diagtime: months from diagnosis to randomisation; age: in years; prior: prior therapy 0=no, 10=yes.

Source: D Kalbfleisch and RL Prentice (1980), The Statistical Analysis of Failure Time Data. Wiley, New York

library(survival,MASS)
attach(veteran)
head(veteran)
##   trt celltype time status karno diagtime age prior
## 1   1 squamous   72      1    60        7  69     0
## 2   1 squamous  411      1    70        5  64    10
## 3   1 squamous  228      1    60        3  38     0
## 4   1 squamous  126      1    60        9  63    10
## 5   1 squamous  118      1    70       11  65    10
## 6   1 squamous   10      1    20        5  49     0
length(trt)
## [1] 137

Abordare de tip Kaplan-Meier

S=Surv(time, status)
summary(S)
##       time           status      
##  Min.   :  1.0   Min.   :0.0000  
##  1st Qu.: 25.0   1st Qu.:1.0000  
##  Median : 80.0   Median :1.0000  
##  Mean   :121.6   Mean   :0.9343  
##  3rd Qu.:144.0   3rd Qu.:1.0000  
##  Max.   :999.0   Max.   :1.0000
KM<-survfit(S~trt)
KM
## Call: survfit(formula = S ~ trt)
## 
##        n events median 0.95LCL 0.95UCL
## trt=1 69     64  103.0      59     132
## trt=2 68     64   52.5      44      95
summary(KM)
## Call: survfit(formula = S ~ trt)
## 
##                 trt=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     3     69       1   0.9855  0.0144      0.95771        1.000
##     4     68       1   0.9710  0.0202      0.93223        1.000
##     7     67       1   0.9565  0.0246      0.90959        1.000
##     8     66       2   0.9275  0.0312      0.86834        0.991
##    10     64       2   0.8986  0.0363      0.83006        0.973
##    11     62       1   0.8841  0.0385      0.81165        0.963
##    12     61       2   0.8551  0.0424      0.77592        0.942
##    13     59       1   0.8406  0.0441      0.75849        0.932
##    16     58       1   0.8261  0.0456      0.74132        0.921
##    18     57       2   0.7971  0.0484      0.70764        0.898
##    20     55       1   0.7826  0.0497      0.69109        0.886
##    21     54       1   0.7681  0.0508      0.67472        0.874
##    22     53       1   0.7536  0.0519      0.65851        0.862
##    27     51       1   0.7388  0.0529      0.64208        0.850
##    30     50       1   0.7241  0.0539      0.62580        0.838
##    31     49       1   0.7093  0.0548      0.60967        0.825
##    35     48       1   0.6945  0.0556      0.59368        0.812
##    42     47       1   0.6797  0.0563      0.57782        0.800
##    51     46       1   0.6650  0.0570      0.56209        0.787
##    52     45       1   0.6502  0.0576      0.54649        0.774
##    54     44       2   0.6206  0.0587      0.51565        0.747
##    56     42       1   0.6059  0.0591      0.50040        0.734
##    59     41       1   0.5911  0.0595      0.48526        0.720
##    63     40       1   0.5763  0.0598      0.47023        0.706
##    72     39       1   0.5615  0.0601      0.45530        0.693
##    82     38       1   0.5467  0.0603      0.44049        0.679
##    92     37       1   0.5320  0.0604      0.42577        0.665
##    95     36       1   0.5172  0.0605      0.41116        0.651
##   100     34       1   0.5020  0.0606      0.39615        0.636
##   103     32       1   0.4863  0.0607      0.38070        0.621
##   105     31       1   0.4706  0.0608      0.36537        0.606
##   110     30       1   0.4549  0.0607      0.35018        0.591
##   117     29       2   0.4235  0.0605      0.32017        0.560
##   118     27       1   0.4079  0.0602      0.30537        0.545
##   122     26       1   0.3922  0.0599      0.29069        0.529
##   126     24       1   0.3758  0.0596      0.27542        0.513
##   132     23       1   0.3595  0.0592      0.26031        0.496
##   139     22       1   0.3432  0.0587      0.24535        0.480
##   143     21       1   0.3268  0.0582      0.23057        0.463
##   144     20       1   0.3105  0.0575      0.21595        0.446
##   151     19       1   0.2941  0.0568      0.20151        0.429
##   153     18       1   0.2778  0.0559      0.18725        0.412
##   156     17       1   0.2614  0.0550      0.17317        0.395
##   162     16       2   0.2288  0.0527      0.14563        0.359
##   177     14       1   0.2124  0.0514      0.13218        0.341
##   200     12       1   0.1947  0.0501      0.11761        0.322
##   216     11       1   0.1770  0.0486      0.10340        0.303
##   228     10       1   0.1593  0.0468      0.08956        0.283
##   250      9       1   0.1416  0.0448      0.07614        0.263
##   260      8       1   0.1239  0.0426      0.06318        0.243
##   278      7       1   0.1062  0.0400      0.05076        0.222
##   287      6       1   0.0885  0.0371      0.03896        0.201
##   314      5       1   0.0708  0.0336      0.02793        0.180
##   384      4       1   0.0531  0.0295      0.01788        0.158
##   392      3       1   0.0354  0.0244      0.00917        0.137
##   411      2       1   0.0177  0.0175      0.00256        0.123
##   553      1       1   0.0000     NaN           NA           NA
## 
##                 trt=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     68       2   0.9706  0.0205      0.93125        1.000
##     2     66       1   0.9559  0.0249      0.90830        1.000
##     7     65       2   0.9265  0.0317      0.86647        0.991
##     8     63       2   0.8971  0.0369      0.82766        0.972
##    13     61       1   0.8824  0.0391      0.80900        0.962
##    15     60       2   0.8529  0.0429      0.77278        0.941
##    18     58       1   0.8382  0.0447      0.75513        0.930
##    19     57       2   0.8088  0.0477      0.72056        0.908
##    20     55       1   0.7941  0.0490      0.70360        0.896
##    21     54       1   0.7794  0.0503      0.68684        0.884
##    24     53       2   0.7500  0.0525      0.65383        0.860
##    25     51       3   0.7059  0.0553      0.60548        0.823
##    29     48       1   0.6912  0.0560      0.58964        0.810
##    30     47       1   0.6765  0.0567      0.57394        0.797
##    31     46       1   0.6618  0.0574      0.55835        0.784
##    33     45       1   0.6471  0.0580      0.54289        0.771
##    36     44       1   0.6324  0.0585      0.52754        0.758
##    43     43       1   0.6176  0.0589      0.51230        0.745
##    44     42       1   0.6029  0.0593      0.49717        0.731
##    45     41       1   0.5882  0.0597      0.48216        0.718
##    48     40       1   0.5735  0.0600      0.46724        0.704
##    49     39       1   0.5588  0.0602      0.45244        0.690
##    51     38       2   0.5294  0.0605      0.42313        0.662
##    52     36       2   0.5000  0.0606      0.39423        0.634
##    53     34       1   0.4853  0.0606      0.37993        0.620
##    61     33       1   0.4706  0.0605      0.36573        0.606
##    73     32       1   0.4559  0.0604      0.35163        0.591
##    80     31       2   0.4265  0.0600      0.32373        0.562
##    84     28       1   0.4112  0.0597      0.30935        0.547
##    87     27       1   0.3960  0.0594      0.29509        0.531
##    90     25       1   0.3802  0.0591      0.28028        0.516
##    95     24       1   0.3643  0.0587      0.26560        0.500
##    99     23       2   0.3326  0.0578      0.23670        0.467
##   111     20       2   0.2994  0.0566      0.20673        0.434
##   112     18       1   0.2827  0.0558      0.19203        0.416
##   133     17       1   0.2661  0.0550      0.17754        0.399
##   140     16       1   0.2495  0.0540      0.16326        0.381
##   164     15       1   0.2329  0.0529      0.14920        0.363
##   186     14       1   0.2162  0.0517      0.13538        0.345
##   201     13       1   0.1996  0.0503      0.12181        0.327
##   231     12       1   0.1830  0.0488      0.10851        0.308
##   242     10       1   0.1647  0.0472      0.09389        0.289
##   283      9       1   0.1464  0.0454      0.07973        0.269
##   340      8       1   0.1281  0.0432      0.06609        0.248
##   357      7       1   0.1098  0.0407      0.05304        0.227
##   378      6       1   0.0915  0.0378      0.04067        0.206
##   389      5       1   0.0732  0.0344      0.02912        0.184
##   467      4       1   0.0549  0.0303      0.01861        0.162
##   587      3       1   0.0366  0.0251      0.00953        0.140
##   991      2       1   0.0183  0.0180      0.00265        0.126
##   999      1       1   0.0000     NaN           NA           NA
plot(KM,xlab='timp',ylab=' survival prob')
abline(h=0.5)
abline(v=103)
abline(v=52.5)

Se observa ce timpii mediani de supravietuire in cele doua grupuri (grupul trt=1 are timpul median 103, iar grupul trt=2 are timpul median 52.5 sap)

survdiff(S~trt)
## Call:
## survdiff(formula = S ~ trt)
## 
##        N Observed Expected (O-E)^2/E (O-E)^2/V
## trt=1 69       64     64.5   0.00388   0.00823
## trt=2 68       64     63.5   0.00394   0.00823
## 
##  Chisq= 0  on 1 degrees of freedom, p= 0.9

Se observa ca exista diferenta semnificativa intre cele doua grupuri.

Abordare folosind metodele parametrice

exp0=survreg(S~1, dist="exponential")
 wei0=survreg(S~1, dist="w")
 teststat=2*(wei0$log[2]-exp0$log[2])
 pchisq(teststat,df=1,lower.tail=F)
## [1] 0.01234947

Se obtine ca Weibul este mai potrivit, astfel se pot determina parametrii acestui model.

Selectia variabilelor din modelul Cox printr-o metoda automata bazata pe determinarea modelului cu cel mai mic AIC

Metoda urmatoare va initializa un model Cox fara variabile, notat aici cu phm_0 si va construi, adaugand/eliminand pe rand cate o variabila alte modele (adaugand variabile pana la scalari de gradul 2). Pentru fiecare dintre aceste model construite automat, va calcula AIC-ul si va selecta mereu modelul cu cel mai mic AIC, pe care va incerca la urmatorul pas sa-l imbunatateasca si astfel, sa scada AIC. Daca acest lucru nu mai este posibil, procedeul automat se va opri si astfel vom obtine “cel mai potrivit model”.

Iata cum se aplica aceasta metoda automata de selectie a variabilelor pentru setul de date veteran:

require(MASS)
## Loading required package: MASS
Scope=list(upper=~(trt+celltype+karno+diagtime+age+prior)^2,lowe=~1)
phm_0=coxph(Surv(time,status)~1)
phm_f=stepAIC(phm_0,Scope,direction="both")
## Start:  AIC=1010.9
## Surv(time, status) ~ 1
## 
##            Df     AIC
## + karno     1  970.87
## + celltype  3  992.05
## <none>        1010.90
## + diagtime  1 1011.99
## + age       1 1012.27
## + prior     1 1012.38
## + trt       1 1012.89
## 
## Step:  AIC=970.87
## Surv(time, status) ~ karno
## 
##            Df     AIC
## + celltype  3  959.53
## <none>         970.87
## + trt       1  971.93
## + age       1  972.80
## + prior     1  972.82
## + diagtime  1  972.85
## - karno     1 1010.90
## 
## Step:  AIC=959.53
## Surv(time, status) ~ karno + celltype
## 
##                  Df    AIC
## <none>              959.53
## + trt             1 959.83
## + age             1 961.09
## + prior           1 961.27
## + diagtime        1 961.35
## + celltype:karno  3 962.51
## - celltype        3 970.87
## - karno           1 992.05

Se obesrva cum, prin semnul + se specifica adaugarea unei variabile, iar prin semnul -, eliminarea ei; in dreptul fiecarei adaugari/eliminari a unei variabile, este calculat AIC-ul modelului. In cazul setului de date veteran, procedura automata ne arata ca modelul optim este cel care contine variabilel karno si celltype si nicio alta adaugare/eliminare de variabila nu produce un model cu un AIC mai mic.

Deci, modelul Cox optim aici este

m1=coxph(formula=Surv(time,status)~karno+celltype)
m1
## Call:
## coxph(formula = Surv(time, status) ~ karno + celltype)
## 
##                        coef exp(coef)  se(coef)      z        p
## karno             -0.031057  0.969421  0.005177 -5.999 1.99e-09
## celltypesmallcell  0.715334  2.044870  0.252686  2.831  0.00464
## celltypeadeno      1.157733  3.182711  0.292937  3.952 7.74e-05
## celltypelarge      0.325645  1.384923  0.276680  1.177  0.23921
## 
## Likelihood ratio test=59.37  on 4 df, p=3.931e-12
## n= 137, number of events= 128

Acest procedeu automat de selectie a variabilelor se poate face si in directie opusa celui anterior, si anume: se pleaca de la un model initial phm_0 ce contine variabilele modelului cu scalari inclusiv de ordin 2 si, prin eliminarea variabilelor, se incearca obtinerea unui model cu AIC mai mic. Desigur ca este o procedura mai complicata, ce presupune mai multi pasi. Iata cum arata aceasta procedura pe setul de date veteran:

phm_0=coxph(Surv(time,status)~(trt+celltype+karno+diagtime+age+prior)^2)
phm_f=stepAIC(phm_0,Scope,direction="both")
## Start:  AIC=967.53
## Surv(time, status) ~ (trt + celltype + karno + diagtime + age + 
##     prior)^2
## 
##                     Df    AIC
## - celltype:karno     3 964.80
## - celltype:prior     3 964.84
## - celltype:age       3 964.92
## - trt:age            1 965.53
## - diagtime:age       1 965.53
## - trt:celltype       3 965.96
## - karno:diagtime     1 966.17
## - age:prior          1 966.48
## - diagtime:prior     1 967.00
## <none>                 967.53
## - trt:diagtime       1 967.71
## - trt:karno          1 969.79
## - trt:prior          1 971.78
## - karno:age          1 971.91
## - karno:prior        1 972.93
## - celltype:diagtime  3 976.25
## 
## Step:  AIC=964.8
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:celltype + trt:karno + trt:diagtime + trt:age + 
##     trt:prior + celltype:diagtime + celltype:age + celltype:prior + 
##     karno:diagtime + karno:age + karno:prior + diagtime:age + 
##     diagtime:prior + age:prior
## 
##                     Df    AIC
## - celltype:age       3 961.72
## - trt:celltype       3 962.28
## - celltype:prior     3 962.57
## - diagtime:age       1 962.82
## - trt:age            1 962.88
## - karno:diagtime     1 963.37
## - age:prior          1 963.70
## - diagtime:prior     1 964.32
## - trt:diagtime       1 964.76
## <none>                 964.80
## - karno:age          1 967.51
## + celltype:karno     3 967.53
## - trt:karno          1 967.82
## - karno:prior        1 969.70
## - trt:prior          1 969.85
## - celltype:diagtime  3 972.90
## 
## Step:  AIC=961.72
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:celltype + trt:karno + trt:diagtime + trt:age + 
##     trt:prior + celltype:diagtime + celltype:prior + karno:diagtime + 
##     karno:age + karno:prior + diagtime:age + diagtime:prior + 
##     age:prior
## 
##                     Df    AIC
## - celltype:prior     3 959.00
## - trt:age            1 959.72
## - diagtime:age       1 959.74
## - karno:diagtime     1 960.19
## - trt:celltype       3 960.21
## - age:prior          1 960.27
## - diagtime:prior     1 961.13
## <none>                 961.72
## - trt:diagtime       1 962.37
## - karno:age          1 963.06
## - trt:karno          1 964.61
## + celltype:age       3 964.80
## + celltype:karno     3 964.92
## - karno:prior        1 967.72
## - trt:prior          1 967.85
## - celltype:diagtime  3 969.69
## 
## Step:  AIC=959
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:celltype + trt:karno + trt:diagtime + trt:age + 
##     trt:prior + celltype:diagtime + karno:diagtime + karno:age + 
##     karno:prior + diagtime:age + diagtime:prior + age:prior
## 
##                     Df    AIC
## - trt:celltype       3 956.95
## - trt:age            1 957.00
## - diagtime:age       1 957.00
## - age:prior          1 957.08
## - karno:diagtime     1 957.31
## - diagtime:prior     1 957.88
## <none>                 959.00
## - trt:diagtime       1 959.38
## - karno:age          1 960.24
## - trt:karno          1 960.87
## + celltype:prior     3 961.72
## + celltype:karno     3 962.18
## + celltype:age       3 962.57
## - trt:prior          1 963.54
## - karno:prior        1 965.39
## - celltype:diagtime  3 966.55
## 
## Step:  AIC=956.95
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:age + trt:prior + 
##     celltype:diagtime + karno:diagtime + karno:age + karno:prior + 
##     diagtime:age + diagtime:prior + age:prior
## 
##                     Df    AIC
## - age:prior          1 954.97
## - diagtime:age       1 955.00
## - trt:age            1 955.11
## - karno:diagtime     1 955.31
## - diagtime:prior     1 956.11
## <none>                 956.95
## - trt:diagtime       1 957.08
## + trt:celltype       3 959.00
## + celltype:age       3 959.69
## - trt:karno          1 960.05
## + celltype:prior     3 960.21
## - karno:age          1 960.44
## - trt:prior          1 960.94
## + celltype:karno     3 961.40
## - karno:prior        1 963.99
## - celltype:diagtime  3 964.80
## 
## Step:  AIC=954.97
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:age + trt:prior + 
##     celltype:diagtime + karno:diagtime + karno:age + karno:prior + 
##     diagtime:age + diagtime:prior
## 
##                     Df    AIC
## - diagtime:age       1 953.00
## - trt:age            1 953.14
## - karno:diagtime     1 953.33
## - diagtime:prior     1 954.16
## <none>                 954.97
## - trt:diagtime       1 955.08
## + age:prior          1 956.95
## + trt:celltype       3 957.08
## + celltype:age       3 957.69
## - trt:karno          1 958.06
## + celltype:prior     3 958.29
## - karno:age          1 958.63
## - trt:prior          1 958.99
## + celltype:karno     3 959.40
## - karno:prior        1 962.13
## - celltype:diagtime  3 962.80
## 
## Step:  AIC=953
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:age + trt:prior + 
##     celltype:diagtime + karno:diagtime + karno:age + karno:prior + 
##     diagtime:prior
## 
##                     Df    AIC
## - trt:age            1 951.22
## - karno:diagtime     1 951.38
## - diagtime:prior     1 952.23
## <none>                 953.00
## - trt:diagtime       1 953.16
## + diagtime:age       1 954.97
## + age:prior          1 955.00
## + trt:celltype       3 955.13
## - trt:karno          1 956.08
## + celltype:prior     3 956.31
## + celltype:age       3 956.33
## - karno:age          1 956.69
## - trt:prior          1 957.02
## + celltype:karno     3 957.40
## - karno:prior        1 960.18
## - celltype:diagtime  3 960.80
## 
## Step:  AIC=951.22
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:prior + celltype:diagtime + 
##     karno:diagtime + karno:age + karno:prior + diagtime:prior
## 
##                     Df    AIC
## - karno:diagtime     1 949.53
## - diagtime:prior     1 950.34
## <none>                 951.22
## - trt:diagtime       1 951.28
## + trt:age            1 953.00
## + trt:celltype       3 953.14
## + diagtime:age       1 953.14
## + age:prior          1 953.22
## - trt:karno          1 954.34
## + celltype:prior     3 954.41
## - karno:age          1 955.04
## + celltype:age       3 955.10
## - trt:prior          1 955.39
## + celltype:karno     3 955.41
## - karno:prior        1 958.35
## - celltype:diagtime  3 958.90
## 
## Step:  AIC=949.53
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:prior + celltype:diagtime + 
##     karno:age + karno:prior + diagtime:prior
## 
##                     Df    AIC
## - diagtime:prior     1 948.49
## - trt:diagtime       1 949.38
## <none>                 949.53
## + karno:diagtime     1 951.22
## + trt:age            1 951.38
## + diagtime:age       1 951.43
## + trt:celltype       3 951.49
## + age:prior          1 951.53
## - trt:karno          1 952.52
## + celltype:prior     3 953.11
## - karno:age          1 953.22
## + celltype:age       3 953.45
## - trt:prior          1 953.46
## + celltype:karno     3 953.99
## - celltype:diagtime  3 956.99
## - karno:prior        1 957.02
## 
## Step:  AIC=948.49
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:diagtime + trt:prior + celltype:diagtime + 
##     karno:age + karno:prior
## 
##                     Df    AIC
## - trt:diagtime       1 948.07
## <none>                 948.49
## + diagtime:prior     1 949.53
## + trt:celltype       3 950.14
## + karno:diagtime     1 950.34
## + diagtime:age       1 950.37
## + trt:age            1 950.41
## + age:prior          1 950.48
## - trt:karno          1 950.75
## - karno:age          1 951.51
## - trt:prior          1 952.23
## + celltype:prior     3 952.25
## + celltype:age       3 952.28
## + celltype:karno     3 953.08
## - celltype:diagtime  3 956.45
## - karno:prior        1 956.61
## 
## Step:  AIC=948.07
## Surv(time, status) ~ trt + celltype + karno + diagtime + age + 
##     prior + trt:karno + trt:prior + celltype:diagtime + karno:age + 
##     karno:prior
## 
##                     Df    AIC
## <none>                 948.07
## + trt:diagtime       1 948.49
## + diagtime:prior     1 949.38
## + karno:diagtime     1 949.88
## + diagtime:age       1 949.94
## + trt:age            1 949.99
## + age:prior          1 950.07
## + trt:celltype       3 950.22
## - trt:prior          1 950.37
## - karno:age          1 950.54
## - trt:karno          1 950.61
## + celltype:age       3 951.25
## + celltype:prior     3 951.83
## + celltype:karno     3 952.79
## - celltype:diagtime  3 954.47
## - karno:prior        1 955.66

In final, modelul Cox obtinut ca fiind optima este cel care contine variabilele trt, celltype, karno, diagtime, age, prior si termenii de tip interactiune \(trt*prior, celltype*diagtime, karno*age, karno*prior\).

m2=coxph(formula=Surv(time,status)~trt+celltype+karno+diagtime+age+prior+ trt*prior+celltype*diagtime+karno*age+karno*prior)
m2 
## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno + 
##     diagtime + age + prior + trt * prior + celltype * diagtime + 
##     karno * age + karno * prior)
## 
##                                  coef  exp(coef)   se(coef)      z        p
## trt                         0.6348091  1.8866619  0.2535320  2.504 0.012285
## celltypesmallcell           1.6723044  5.3244232  0.3724251  4.490 7.11e-06
## celltypeadeno               1.8499980  6.3598066  0.4403527  4.201 2.66e-05
## celltypelarge               0.8433017  2.3240276  0.4669327  1.806 0.070911
## karno                      -0.0768772  0.9260036  0.0325129 -2.365 0.018054
## diagtime                    0.0537421  1.0552125  0.0161643  3.325 0.000885
## age                        -0.0576226  0.9440061  0.0300640 -1.917 0.055280
## prior                       0.3514223  1.4210873  0.0932280  3.769 0.000164
## trt:prior                  -0.0783404  0.9246496  0.0433779 -1.806 0.070918
## celltypesmallcell:diagtime -0.0731561  0.9294557  0.0205542 -3.559 0.000372
## celltypeadeno:diagtime     -0.0535890  0.9478216  0.0428530 -1.251 0.211106
## celltypelarge:diagtime     -0.0450452  0.9559543  0.0416947 -1.080 0.279983
## karno:age                   0.0008366  1.0008370  0.0005220  1.603 0.108980
## karno:prior                -0.0038007  0.9962065  0.0012871 -2.953 0.003148
## 
## Likelihood ratio test=88.29  on 14 df, p=7.972e-13
## n= 137, number of events= 128

Se observa ca se obtine un model mult mai complex, in care nu toate variabilele modelului sunt semnificative, deci si acest model poate fi imbunatatit.

Selectia variabilelor folosind metoda Step-Wise Selection

Vom nota variabilele din setul de date veteran, pentru a le folosi mai usor in cele ce urmeaza.

x1=trt
x2=celltype
x3=karno
x4=diagtime
x5=age
x6=prior
S=Surv(veteran$time,veteran$status)
S
##   [1]  72  411  228  126  118   10   82  110  314  100+  42    8  144   25+  11 
##  [16]  30  384    4   54   13  123+  97+ 153   59  117   16  151   22   56   21 
##  [31]  18  139   20   31   52  287   18   51  122   27   54    7   63  392   10 
##  [46]   8   92   35  117  132   12  162    3   95  177  162  216  553  278   12 
##  [61] 260  200  156  182+ 143  105  103  250  100  999  112   87+ 231+ 242  991 
##  [76] 111    1  587  389   33   25  357  467  201    1   30   44  283   15   25 
##  [91] 103+  21   13   87    2   20    7   24   99    8   99   61   25   95   80 
## [106]  51   29   24   18   83+  31   51   90   52   73    8   36   48    7  140 
## [121] 186   84   19   45   80   52  164   19   53   15   43  340  133  111  231 
## [136] 378   49

Vom defini un prag de intrare \(\text{pi}\) si un prag de eliminare/iesire \(\text{pis}\) al variabilelor in model (se va folosi logaritmul acestor praguri). De obicei, decizi de a pastra o variabila in model este luata pe baza p-valorii, care este comparata cu nivelul de semnificatie 0.05. Aici valorile pe baza carora se iau deciziile, \(\text{pi}\) si \(\text{pis}\), au valori mai mari, tocmai pentru a nu elimina variabile din model, care pot sa fie importnate, din punct de vedere statistic.

 pi=0.15
logpi=log(pi)
logpi
## [1] -1.89712
pis=0.2
logpis=log(pis)
logpis
## [1] -1.609438

Vom construi modelul Cox fara variabile, notat \(m_0\) si vom sublinia valoare functie de verosimilitate pentru acest model:

m0=coxph(S~1) 
m0
## Call:  coxph(formula = S ~ 1)
## 
## Null model
##   log likelihood= -505.4491 
##   n= 137
ll0=m0$logl[1]
ll0
## [1] -505.4491

Vom construi modele Cox introducand pe rand cate o variabila \(x_i, i=1,\ldots, 6\) si vom face teste de diferentiere intre aceste modele construite si modelul initial \(m_0\).

m1=coxph(S~x1)
m2=coxph(S~x2)
m3=coxph(S~x3)
m4=coxph(S~x4)
m5=coxph(S~x5)
m6=coxph(S~x6)
ll1=m1$logl[2]
ll2=m2$logl[2]
ll3=m3$logl[2]
ll4=m4$logl[2]
ll5=m5$logl[2]
l65=m6$logl[2]
ll6=m6$logl[2]
pchisq(2*(ll1-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.08145637
pchisq(2*(ll2-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -14.29334
pchisq(2*(ll3-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -23.13312
pchisq(2*(ll4-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -1.079618
pchisq(2*(ll5-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.8461987
pchisq(2*(ll6-ll0),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.7519108

Vom selecta modelul cu cel mai mic log p-value, si anume modelul cu \(m_3\). Acest log p-value=-23.13 este mai mic si decat -1.9 (limita de intrare), deci continuam modelul cu x3, care va fi definit ca un nou model initial, notat cu \(m_{00}\), la care vom adauga pe rand cate o variabila si din noi, vom face teste de diferentiere dintre aceste modele construite si modelul initial \(m_{00}\).

m00=coxph(S~x3)
m01=coxph(S~x3+x1)
m02=coxph(S~x3+x2)
m04=coxph(S~x3+x4)
m05=coxph(S~x3+x5)
m06=coxph(S~x3+x6)
ll3=m00$logl[2]
ll00=m00$logl[2]
ll01=m01$logl[2]
ll02=m02$logl[2]
ll04=m04$logl[2]
ll05=m05$logl[2]
ll06=m06$logl[2]
pchisq(2*(ll01-ll00),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -1.098312
pchisq(2*(ll02-ll00),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -10.37362
pchisq(2*(ll04-ll00),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.1075621
pchisq(2*(ll05-ll00),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.2322606
pchisq(2*(ll06-ll00),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.1918498

Se observa ca cel mai maic log p-value este \(-10.37\) (la modelul \(m_2\)), deci variabila \(x_2\) se va include in model, unde \(x_2\)=celltype. Din nou pragul de intrare este verificat.

Continuam procedura de selectie a variabilelor, repatand pasii anteriori: definirea unui nou model initial, adaugarea de noi variabile, teste de diferentiere intre modele.

m000=coxph(S~x3+x2)
m001=coxph(S~x3+x2+x1)
m004=coxph(S~x3+x2+x4)
m005=coxph(S~x3+x2+x5)
m006=coxph(S~x3+x2+x6)
ll000=m000$logl[2]
ll001=m001$logl[2]
ll004=m004$logl[2]
ll005=m005$logl[2]
ll006=m006$logl[2]
pchisq(2*(ll001-ll000),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -1.647005
pchisq(2*(ll004-ll000),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.394469
pchisq(2*(ll005-ll000),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.6772833
pchisq(2*(ll006-ll000),lower.tail=FALSE,log.p=TRUE,df=1)
## [1] -0.4896083

Se obesrva ca modelul cu \(x_1\)=trt are cel mai mic log-p value, dar nu este mai mic decat pragul de intrare (\(p{i}=-1.9\)), deci nu selectam acesta variabila.

Deci, vom studia modelul Cox cu variabilele \(x_2\) si \(x_3\).

m000=coxph(S~x3+x2)
m000
## Call:
## coxph(formula = S ~ x3 + x2)
## 
##                  coef exp(coef)  se(coef)      z        p
## x3          -0.031057  0.969421  0.005177 -5.999 1.99e-09
## x2smallcell  0.715334  2.044870  0.252686  2.831  0.00464
## x2adeno      1.157733  3.182711  0.292937  3.952 7.74e-05
## x2large      0.325645  1.384923  0.276680  1.177  0.23921
## 
## Likelihood ratio test=59.37  on 4 df, p=3.931e-12
## n= 137, number of events= 128