## Exercise 19.1 Leukhaemia survival analysis patients with acute myelogenous leukaemia
## original data from Feigl, P. & Zelen, M. (1965) Estimation of exponential survival probabilities with concomitant information. Biometrics 21, 826–838.
## the vector 'ag' indicates presence or absence in patients of Auer rods and/or significant granulation of the leukaemic cells in the bone marrow at the time of diagnosis

### preliminaries

library(MASS)
library(survival)
data(leuk)
leuk
#      wbc      ag time
#1    2300 present   65
#2     750 present  156
#3    4300 present  100
#4    2600 present  134
#5    6000 present   16
#6   10500 present  108
#7   10000 present  121
#8   17000 present    4
#9    5400 present   39
#10   7000 present  143
#11   9400 present   56
#12  32000 present   26
#13  35000 present   22
#14 100000 present    1
#15 100000 present    1
#16  52000 present    5
#17 100000 present   65
#18   4400  absent   56
#19   3000  absent   65
#20   4000  absent   17
#21   1500  absent    7
#22   9000  absent   16
#23   5300  absent   22
#24  10000  absent    3
#25  19000  absent    4
#26  27000  absent    2
#27  28000  absent    3
#28  31000  absent    8
#29  26000  absent    4
#30  21000  absent    3
#31  79000  absent   30
#32 100000  absent    4
#33 100000  absent   43


###  part 1
with(leuk, plot(survfit(Surv(time) ~ ag),xlab="Time after diagnosis (weeks)",ylab="Survival probability",lty = c(4,1), col = c("black","blue"),lwd=2))
legend(100,0.9,c("without Auer rods","with Auer rods"),lty = c(4,1), col = c("black","blue"),lwd=2,bty=FALSE)


### part 2
cox<-with(leuk,coxph(Surv(time) ~ ag+log(wbc)))
summary(cox)
#Call:
#coxph(formula = Surv(time) ~ ag + log(wbc))
#
#  n= 33, number of events= 33 
#
#             coef exp(coef) se(coef)      z Pr(>|z|)   
#agpresent -1.0691    0.3433   0.4293 -2.490  0.01276 * 
#log(wbc)   0.3677    1.4444   0.1360  2.703  0.00687 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#          exp(coef) exp(-coef) lower .95 upper .95
#agpresent    0.3433     2.9126     0.148    0.7964
#log(wbc)     1.4444     0.6923     1.106    1.8857
#
#Concordance= 0.726  (se = 0.047 )
#Likelihood ratio test= 15.64  on 2 df,   p=4e-04
#Wald test            = 15.06  on 2 df,   p=5e-04
#Score (logrank) test = 16.49  on 2 df,   p=3e-04