rm(list=ls())

library(survival)
attach(veteran)

pdf("Survival1.pdf")
tt = time[trt == 1]
ss = status[trt == 1]
out = survfit(Surv(tt,ss) ~ 1)
plot(out$time,out$surv,xlim=c(0,450),type="s",lwd=3,col="purple",
         xlab="",ylab="",cex.axis=1.4)
lines(out$time,out$lower,type="s",lwd=3,col="white")   
lines(out$time,out$upper,type="s",lwd=3,col="white")

polygon(c(rev(out$time),out$time),c(rev(out$upper),out$lower),
         density=NULL,angle=90,border=NA,col="thistle") 
lines(out$time,out$surv,type="s",lwd=3,col="purple")
title(ylab="Survival", line=2.5, cex.lab=1.8)
title(xlab="Time",line=3, cex.lab=2)

dev.off()

pdf("Survival2.pdf")
out = survfit(Surv(time,status) ~ trt)  ##compare treated and untreated
n1 = out$strata[1]
n2 = out$strata[2]
n = n1 + n2
plot(out$time[1:n1],out$surv[1:n1],type="s",lwd=3,xlim=c(0,1000),
               col="red",ylab="",xlab="",cex.axis=1.4)
lines(out$time[(n1+1):n],out$surv[(n1+1):n],type="s",lwd=3,col="purple")
title(ylab="Survival", line=2.5, cex.lab=1.8)
title(xlab="Time",line=3, cex.lab=2)
dev.off()

test = survdiff(Surv(time,status) ~ trt)
survdiff(formula = Surv(time, status) ~ trt)

# print(test)

# Call:
# survdiff(formula = Surv(time, status) ~ 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 


out = coxph(Surv(time,status) ~ trt + karno + diagtime + age + prior)

# print(summary(out))

## Call:
## coxph(formula = Surv(time, status) ~ trt + karno + diagtime + 
##     age + prior)
## 
##   n= 137, number of events= 128 
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## trt       0.193053  1.212947  0.186446  1.035    0.300    
## karno    -0.034084  0.966490  0.005341 -6.381 1.76e-10 ***
## diagtime  0.001723  1.001725  0.009003  0.191    0.848    
## age      -0.003883  0.996125  0.009247 -0.420    0.675    
## prior    -0.007764  0.992266  0.022152 -0.350    0.726    
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1