# library(survival)


		#########################
		#				#
		#	 Question 4		#
		#				#
		#########################
			
	
t  <- c(5 ,5 ,8 ,8 ,12,16,23,27,30,33,43)
st <- c(1 ,1 ,1 ,0 ,1 ,0 ,1 ,1 ,1 ,1 ,1 ) 

surv.x <- Surv(t,st) 

plot(survfit(surv.x~1)) 			
summary(survfit(surv.x~1,conf.type="none")) 

summary(survfit(surv.x~1,conf.type="plain")) 
s.14  <- (9/11)*(8/9)*(6/7)				#  K-M estimate for S(14)
sd.14 <- sqrt( 2/(11*9) + 1/(9*8) + 1/(7*6))	# sd for log(S(14))
s.14 * sd.14						# Greenwood sd for S(14)
s.14 + qnorm(.975)*.150	  				# ul for S(14)
s.14 - qnorm(.975)*.150	  				# ll for S(14)

summary(survfit(surv.x~1,conf.type="log")) 	
exp(log(s.14) + qnorm(.975)*sd.14)			# ul for S(14)
exp(log(s.14) - qnorm(.975)*sd.14)			# ll for S(14)


summary(survfit(surv.x~1,conf.type="log-log")) 	
sd.lglg <-  - sd.14 / log(s.14)			# sd for log(-log(S(14))
s.lglg  <- log(-log(s.14))				# est of log(-log(S(14))
exp(-exp(s.lglg - qnorm(.975)*sd.lglg))		# ul for S(14)
exp(-exp(s.lglg + qnorm(.975)*sd.lglg)) 		# ll for S(14)

print(survfit(surv.x~1),show.rmean=T)	# K-M based inference
summary(survfit(surv.x~1)) 	

a <- summary(survfit(surv.x~1))
t.i    <- c(0,a$time,161)
surv.i <- c(1,a$surv)
sum(diff(t.i) * surv.i)	# Restricted mean estimation





		#########################
		#				#
		#	 Question 5		#
		#				#
		#########################



# Survival estimation of censored exponential distribution 

n <- 4000

par(mfcol=c(1,1))
time <- rexp(n,1/10)
cens <- rexp(n,1/30)
stat <- time <= cens
x <- apply(cbind(time,cens),1,min)

plot(survfit(Surv(x,stat)~1))
lines(sort(time),1-pexp(sort(time),1/10),col=3)

# Integrated hazard estimation 

a <- summary(survfit(Surv(x,stat)~1))

H.est <- cumsum(a$n.event / a$n.risk)
H.sd  <- sqrt(cumsum(a$n.event / a$n.risk^2))

plot(sort(a$time),H.est,type="b")
lines(sort(a$time),H.est + qnorm(.975)*H.sd)
lines(sort(a$time),H.est - qnorm(.975)*H.sd)
abline(0,1/10,col=3)

#	Examine mean and median estimation


print(survfit(Surv(x,stat)~1),show.rmean=T)
-log(0.5) *10	#   exp(-lambda * t) = 0.5 delta  ==>  median =   - log(0.5) / lambda



# 	An example of an invalid censoring scheme

n        <- 4000
p.cens   <- 0.5
stat     <- c(rep(1,n*(1-p.cens)),rep(0,n*p.cens))

x <- rexp(n,1/10)
plot(survfit(Surv(x,stat)~1))
lines(sort(x),1-pexp(sort(x),1/10),col=3)

# Integrated hazard estimation 

a <- summary(survfit(Surv(x,stat)~1))

H.est <- cumsum(a$n.event / a$n.risk)
H.sd  <- sqrt(cumsum(a$n.event / a$n.risk^2))

plot(sort(a$time),H.est,type="b")
lines(sort(a$time),H.est + qnorm(.975)*H.sd)
lines(sort(a$time),H.est - qnorm(.975)*H.sd)
abline(0,1/10,col=3)

H.low <- lowess(sort(a$time),H.est,f=1/5)
lines(H.low,col=4)
approx(H.low,xout=20)
points(approx(H.low,xout=20),pch=3,col=2,cex=3)

print(survfit(Surv(x,stat)~1),show.rmean=T)




# An example of a similar censoring scheme for which the survival is exp(.1)

n        <- 40000
p.obs    <- runif(n)
obs      <- rbinom(n,1,p.obs)
x 	   <- rexp(n,0.1  / p.obs)
plot(survfit(Surv(x,obs)))
lines(sort(x),1-pexp(sort(x),1/10),col=3)