####### Fitting Cox Proportional Hazards Models #######

library(survival)
library(MASS)
data(motors)
attach(motors)

##### Example I: Motor data
motor.cox <- coxph(Surv(time, cens) ~ temp, motors)
summary(motor.cox)
Call:
coxph(formula = Surv(time, cens) ~ temp, data = motors)

  n= 40 
       coef exp(coef) se(coef)    z       p
temp 0.0919      1.10   0.0274 3.36 0.00079

     exp(coef) exp(-coef) lower .95 upper .95
temp      1.10      0.912      1.04      1.16

Rsquare= 0.472   (max possible= 0.936 )
Likelihood ratio test= 25.6  on 1 df,   p=4.3e-07
Wald test            = 11.3  on 1 df,   p=0.000786
Score (logrank) test = 22.7  on 1 df,   p=1.86e-06
# This Cox model can be used to compare several groups. 
# We see that the hazards (survival) differ significantly for 
# different temperature groups (but we don't know which 
# two temperatures lead to the difference). 

# predict at temperature 200
plot(survfit(motor.cox, newdata = data.frame(temp=200),
          conf.type = "log-log"))
# some info for temperature=130
summary( survfit(motor.cox, newdata = data.frame(temp=130)) )
Call: survfit.coxph(object = motor.cox, newdata = data.frame(temp = 130))

 time n.risk n.event survival  std.err lower 95% CI upper 95% CI
  408     40       4    1.000 0.000254        0.999            1
  504     36       3    1.000 0.000498        0.999            1
 1344     28       2    0.999 0.001910        0.995            1
 1440     26       1    0.998 0.002697        0.993            1
 1764     20       1    0.996 0.005325        0.986            1
 2772     19       1    0.994 0.007920        0.978            1
 3444     18       1    0.991 0.010673        0.971            1
 3542     17       1    0.988 0.013667        0.962            1
 3780     16       1    0.985 0.016976        0.952            1
 4860     15       1    0.981 0.020692        0.941            1
 5196     14       1    0.977 0.024941        0.929            1


###### Example II: VA data 

data(VA)
attach(VA)
# Cox model with treatment only, ignoring other covariates
va.cox1 <- coxph(Surv(stime) ~ treat)
> summary(va.cox1)
Call:
coxph(formula = Surv(stime) ~ treat)

  n= 137 
         coef exp(coef) se(coef)      z    p
treat2 0.0166      1.02    0.175 0.0949 0.92

       exp(coef) exp(-coef) lower .95 upper .95
treat2      1.02      0.984     0.722      1.43

Rsquare= 0   (max possible= 1 )
Likelihood ratio test= 0.01  on 1 df,   p=0.924
Wald test            = 0.01  on 1 df,   p=0.924
Score (logrank) test = 0.01  on 1 df,   p=0.924
# This Cox model can also be used to compare 
# the treatment groups, ignoring other covariates. 
# plot
plot(survfit(va.cox1)) 
# compare with log-rank test
survdiff(Surv(stime, status)~treat)
Call:
survdiff(formula = Surv(stime, status) ~ treat)

         N Observed Expected (O-E)^2/E (O-E)^2/V
treat=1 69       64     64.5   0.00388   0.00823
treat=2 68       64     63.5   0.00394   0.00823

 Chisq= 0  on 1 degrees of freedom, p= 0.928 
# We see that, when comparing two treatments, Cox model 
# with only treatment gives similar result as the log-rank test.
# So Cox models can be used to compare different groups. 
# Note that the above Cox model does NOT adjust for covariates, 
# but we can adjust for covariates (see below). 

# Cox model with treatment, ADJUSTED for covariates
va.cox2 <- coxph(Surv(stime) ~ treat + age + Karn + diag.time + cell)
Call:
coxph(formula = Surv(stime) ~ treat + age + Karn + diag.time + 
    cell)

  n= 137 
              coef exp(coef) se(coef)       z       p
treat2     0.30233      1.35  0.19999  1.5118 1.3e-01
age       -0.00957      0.99  0.00899 -1.0641 2.9e-01
Karn      -0.03060      0.97  0.00529 -5.7801 7.5e-09
diag.time  0.00026      1.00  0.00814  0.0319 9.7e-01
cell2      0.81511      2.26  0.25925  3.1441 1.7e-03
cell3      1.10804      3.03  0.28487  3.8897 1.0e-04
cell4      0.28574      1.33  0.27136  1.0530 2.9e-01

          exp(coef) exp(-coef) lower .95 upper .95
treat2         1.35      0.739     0.914      2.00
age            0.99      1.010     0.973      1.01
Karn           0.97      1.031     0.960      0.98
diag.time      1.00      1.000     0.984      1.02
cell2          2.26      0.443     1.359      3.76
cell3          3.03      0.330     1.733      5.29
cell4          1.33      0.751     0.782      2.27

Rsquare= 0.351   (max possible= 1 )
Likelihood ratio test= 59.3  on 7 df,   p=2.08e-10
Wald test            = 59.7  on 7 df,   p=1.76e-10
Score (logrank) test = 63.4  on 7 df,   p=3.13e-11
# We see that the treatment effect differs when covariates 
# are adjusted. There seems to be a stronger evidence of 
# treatment effect after covariates are adjusted. 
# This is the advantage of Cox model - it can adjust for 
# covariates when comparing several groups. The log-rank 
# test, on the other hand, has the advantage of assigning 
# different weights to survival differences appeared at 
# different time periods. 
# The validity of the proportional hazards assumption 
# will be checked later.