library(survival)

#time.vec <- round(c(rexp(5,0.01), rexp(5,0.05)),1)
#cens.vec <- round(rexp(10,0.005),1)
#stat.vec <- time.vec <= cens.vec
#obs.vec  <- ifelse(stat.vec,time.vec,cens.vec)
#x.vec    <- c(rep(0,5),rep(1,5))



# the sampled values appearing in Exercise 4:

 obs.vec 	<- c(132.0,19.2,63.9,137.4,102.2,16.6,3.2,8.6,12.7,33.0)
 stat.vec 	<- c(1,1,1,0,0,1,1,1,1,1)
 x.vec	<- c(rep(0,5),rep(1,5))

#   Question 1:

Surv(obs.vec,stat.vec)
summary(coxph(Surv(obs.vec,stat.vec)~x.vec))


(a1    <- summary(survfit(Surv(obs.vec,stat.vec))))
(t.lik <- a1$time)
(k     <- length(t.lik))
(beta     <- coxph(Surv(obs.vec,stat.vec)~x.vec)$coef)
(se.beta  <- sqrt(coxph(Surv(obs.vec,stat.vec)~x.vec)$var))

(exp.beta <- ifelse(x.vec==0,1,exp(beta)))
ratio.0 <- rep(NA,k)
ratio.1 <- rep(NA,k)

for(i in 1:k)
{
	ratio.0[i] <- sum(t.lik[i] == obs.vec) 		 / sum(t.lik[i] <= obs.vec)
	ratio.1[i] <- sum(exp.beta[t.lik[i] == obs.vec]) / sum(exp.beta[t.lik[i] <= obs.vec])
}

sum(log(ratio.0))
sum(log(ratio.1))

coxph(Surv(obs.vec,stat.vec)~x.vec)$loglik
2*diff(coxph(Surv(obs.vec,stat.vec)~x.vec)$loglik)
summary(coxph(Surv(obs.vec,stat.vec)~x.vec))

survdiff(Surv(obs.vec,stat.vec)~x.vec,rho=0)  # log-rank test


# Question 2:  95% CI for exp(beta)

exp(beta) 
exp(beta - qnorm(0.975)*se.beta) 
exp(beta + qnorm(0.975)*se.beta) 

summary(coxph(Surv(obs.vec,stat.vec)~x.vec))


# Haashara: a larger sample size yields a more compelling CI:

	time.vec <- round(c(rexp(5000,0.01), rexp(5000,0.05)),1)
	cens.vec <- round(rexp(10000,0.005),1)
	stat.vec <- time.vec <= cens.vec
	obs.vec  <- ifelse(stat.vec,time.vec,cens.vec)
	x.vec    <- c(rep(0,5000),rep(1,5000))
	summary(coxph(Surv(obs.vec,stat.vec)~x.vec))


# Question 3: 

# ==>  refresh sampled obs:

	 obs.vec 	<- c(132.0,19.2,63.9,137.4,102.2,16.6,3.2,8.6,12.7,33.0)
	 stat.vec 	<- c(1,1,1,0,0,1,1,1,1,1)
	 x.vec	<- c(rep(0,5),rep(1,5))

par(mfcol=c(1,2))
plot(survfit(Surv(obs.vec,stat.vec)~x.vec),conf.int=F,main="KM fit")
plot( survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=0),conf.int=F),main="Cox PH model")
lines(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=1),conf.int=F))

t.vec <- seq(0,140,length=100)
lines(t.vec,exp(-0.05*t.vec),col=2)
lines(t.vec,exp(-0.01*t.vec),col=2)


# Question 4: 

summary(coxph(Surv(obs.vec,stat.vec)~x.vec))

summary(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=0) ))

# Q: how do we compute?

summary(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=1) ))

# A: recall S(t) = exp(-H(t)) &  H(t;x) = H(t;0)*exp(bx) 
#		==>  S(t;x) = exp(-H(t;x)) = (exp(-H(t;0)))^exp(bx)
 

a0 <- summary(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=0) ))
a1 <- summary(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec),data.frame(x.vec=1) ))

cbind(Time=a0$time,S0=a0$surv,S1=a1$surv,Ans=a0$surv^exp(beta))


# Haashara: R default 

(a2 <- summary(survfit(coxph(Surv(obs.vec,stat.vec)~x.vec) )))

cbind(Time=a2$time,R.default=a2$surv,Ans=a0$surv^exp(beta*mean(x.vec)))
mean(x.vec)