Cox-model: example from prof Tuan

This one is guided by using material of Prof Nguyen Van Tuan.

R Markdown

This is an example how to proceed for the survival analysis inluding Kaplan-Meierand cox-rank test in R

Step 1: loading the packages

library(survival)

step 2: input data

group = c(1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1,1,0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,0)
time =c(6,6,6,7,10,13,16,22,23,6,9,10,11,17,19,20,25,32,32,34,35,1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)
status = c(0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,0)
wbc =c(2.31,4.06,3.28,4.43,2.96,2.88,3.60,2.32,2.57,3.20,2.80,2.70,2.60,2.16,2.05,2.01,1.78,2.20,2.53,1.47,1.45,2.80,5.00,4.91,4.48,4.01,4.36,2.42,3.49,3.97,3.52,3.05,2.32,3.26,3.49,2.12,1.50,3.06,2.30,2.95,2.73,1.97)

Step 3: produce the result for survival

dat=data.frame(group,time,status,wbc)
baseline = Surv(time, status==0)
km = survfit(baseline ~ 1)
summary(km)

## Call: survfit(formula = baseline ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     42       2    0.952  0.0329       0.8901        1.000
##     2     40       2    0.905  0.0453       0.8202        0.998
##     3     38       1    0.881  0.0500       0.7883        0.985
##     4     37       2    0.833  0.0575       0.7279        0.954
##     5     35       2    0.786  0.0633       0.6709        0.920
##     6     33       3    0.714  0.0697       0.5899        0.865
##     7     29       1    0.690  0.0715       0.5628        0.845
##     8     28       4    0.591  0.0764       0.4588        0.762
##    10     23       1    0.565  0.0773       0.4325        0.739
##    11     21       2    0.512  0.0788       0.3783        0.692
##    12     18       2    0.455  0.0796       0.3227        0.641
##    13     16       1    0.426  0.0795       0.2958        0.615
##    15     15       1    0.398  0.0791       0.2694        0.588
##    16     14       1    0.369  0.0784       0.2437        0.560
##    17     13       1    0.341  0.0774       0.2186        0.532
##    22      9       2    0.265  0.0765       0.1507        0.467
##    23      7       2    0.189  0.0710       0.0909        0.395

Cox-model (predictors of survival )

cox = coxph(Surv(time, status==0) ~ group + wbc)
summary(cox)

## Call:
## coxph(formula = Surv(time, status == 0) ~ group + wbc)
## 
##   n= 42, number of events= 30 
## 
##          coef exp(coef) se(coef)      z Pr(>|z|)    
## group -1.3861    0.2501   0.4248 -3.263   0.0011 ** 
## wbc    1.6909    5.4243   0.3359  5.034  4.8e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## group    0.2501     3.9991    0.1088    0.5749
## wbc      5.4243     0.1844    2.8082   10.4776
## 
## Concordance= 0.852  (se = 0.04 )
## Likelihood ratio test= 46.71  on 2 df,   p=7e-11
## Wald test            = 33.6  on 2 df,   p=5e-08
## Score (logrank) test = 46.07  on 2 df,   p=1e-10

Including Plots

You can also embed plots, for example:

plot(km, xlab="Time to death", ylab="Prob of survival")

Kaplan – Meier plot

```r
plot(km,
xlab="Time to death", ylab="Prob of survival")

plot(km, fun="event",
xlab="Time to death", ylab="Prob of survival")

#If we want to compare by group, then


```r
time = c(6, 6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22,
23, 25, 32, 32, 34, 35, 1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8,
11, 11, 12, 12, 15, 17, 22, 23)
status = c(0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1,
1, 1, 1, 1, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,0)
group = c(1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1,1,
0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,0)
dat = data.frame(time, status, group)
library(survival)
km = survfit(Surv(time, status==0) ~ group); km

## Call: survfit(formula = Surv(time, status == 0) ~ group)
## 
##          n events median 0.95LCL 0.95UCL
## group=0 21     21      8       4      12
## group=1 21      9     23      16      NA

plot(km, lty=c(1,4), lwd=2, xlab="Weeks", ylab="S(t)")
legend("topright", c("Placebo", "Treatment"), lty=c(1,4),
lwd=2)

km

## Call: survfit(formula = Surv(time, status == 0) ~ group)
## 
##          n events median 0.95LCL 0.95UCL
## group=0 21     21      8       4      12
## group=1 21      9     23      16      NA

km.diff = survdiff(Surv(time, status==0) ~ group)
km.diff

## Call:
## survdiff(formula = Surv(time, status == 0) ~ group)
## 
##          N Observed Expected (O-E)^2/E (O-E)^2/V
## group=0 21       21     10.7      9.77      16.8
## group=1 21        9     19.3      5.46      16.8
## 
##  Chisq= 16.8  on 1 degrees of freedom, p= 4e-05