#Censoring Calendar Time

###########################################################3
# setwd("X://temp//Rdata")
library(survival)
library(Hmisc)
library(survival)

#####################################################
# TIME AXIS PLOT
#####################################################
####################################




########################################################
# LIFE TABLE 
########################################################
# KMsurv package and library
# lifetab(tis, ninit, nlost, nevent)
# tis     : end points of time intervals, 
#           length is 1 greater than nlost and nevent
# ninit   : number initially entering the study.  
# nlost   : number lost follow or withdrawn
# nevent  : number the event
tis<- c(0, 1, 2, 3, 4,  5,  6,  7,  8,  9, NA)
nsubs <- c(913, 505, 335, 228, 169, 122,  76,  56,  43, 32)
# ninit <- nsubs[1]
nlos <- c(19,  3,  4,  3,  5,  1,  0,  2, 0, 0)
nwdr <- c(77, 71, 58, 27, 35, 36, 17, 10, 8, 0)
nlost<-nlos+nwdr
nevent <- c(312, 96,  45, 29,   7,   9,   3,  1,  3, 32)
# lifetab(tis, ninit, nlost, nevent)
lifetab(tis, nsubs[1], nlost, nevent)


     
############################################################
# Kaplan-Meier Estimator
############################################################


###################################
#### SIMPLE DATA
###################################
exkm.time <- c(1,3,4,5,5,6,7,7,7,8)
exkm.censor<-c(0,1,0,1,1,0,1,1,0,0) 
data.frame(exkm.time,exkm.censor)
survfit(Surv(exkm.time,exkm.censor))
summary(survfit(Surv(exkm.time,exkm.censor),
	 type=c("kaplan-meier"), 
	 error=c("greenwood"),conf.type=c("log")))
plot(survfit(Surv(exkm.time,exkm.censor)))
plot(survfit(Surv(exkm.time,exkm.censor),  
     conf.type="none"), bty="l", 
     xlab="Time", ylab="Survival Rate")







########################################
### AML-6MP GROUP DATA
########################################
AML01.time<- c(10,7,32,23,22,6,16,34,32,25,11,19,
		19,6,17,35,6,13,9,6,10)
AML01.censor<-c(1,1, 0, 1, 1,1, 1, 0, 0, 0, 0, 0,
 		 0,1, 0, 0,0, 1,0,0, 0)


data.frame(AML01.time,AML01.censor)
survfit(Surv(AML01.time,AML01.censor))
summary(survfit(Surv(AML01.time,AML01.censor),
	 type="kaplan-meier", 
	 error="greenwood",conf.type="plain"))



par(mfrow = c(1,2))
plot(survfit(Surv(AML01.time,AML01.censor), 
	type="kaplan-meier", 
	 error="greenwood", conf.type="none"), 
        bty="l",xlab="Time", ylab="Survival Rate")
plot(survfit(Surv(AML01.time,AML01.censor),
	 type="kaplan-meier", 
	 error="greenwood",conf.type="plain"),
	bty="l", xlab="Time", ylab="Survival Rate")
par(mfrow = c(1,1)) 



#########################################
### BMJast DATA
#########################################
ast.time<-c(6,  13, 21, 30, 31, 37, 38,  47,  49,  50,       
            63, 79, 80, 82, 82, 86, 98, 149, 202, 219)
ast.censor<-c(1,1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,1)   

ast<-data.frame(ast.time,ast.censor)
summary(survfit(Surv(ast.time, ast.censor)), data=ast,
		type=c("kaplan-meier"),                      
		error=c("greenwood"), conf.type=c("plain"))


plot(survfit(Surv(ast.time,ast.censor),
     conf.type="none"), bty="l", 
     xlab="Time", ylab="Survival Rate")

## NA HAZARD ###
temp <- summary(survfit(Surv(ast.time, ast.censor), data=ast,type="fh"))
cbind(temp$time, -log(temp$surv))