# library(survival)
# library(gregmisc)

#########################
#
#	Question 1:
#
#########################


time  <- c(2,7,13,15,20,25,34)
group <- c(0,0,1,1,1,0,0)

(x   <- group == 1)

(w.1 <- (1:7) 	      )
(w.2 <- cumsum(1/(7:1)) )


#	observed values

(obs.1   <-   sum(x*w.1))
(obs.2   <-   sum(x*w.2)) 

# theoretical rank mean and variance

mean(w.1) * sum(x)		
mean(w.2) * sum(x)		

(sum(w.1^2) -  7 * mean(w.1)^2) /(7-1) * ( sum(x^2) - 7 * mean(x)^2 )  
(sum(w.2^2) -  7 * mean(w.2)^2) /(7-1) * ( sum(x^2) - 7 * mean(x)^2 )  

# 	Normal approx p-values for log-rank test

(mn.2 <- mean(w.2) * sum(x)		)
(sd.2 <- sqrt((sum(w.2^2) -  7 * mean(w.2)^2) /(7-1) * ( sum(x^2) - 7 * mean(x)^2 ) ) )

((obs.2 - mn.2) / sd.2)
1-pnorm((sum(x*w.2) - mn.2) / sd.2)


# simulation-based rank mean and variance and p-value

s.1 <- rep(1,10000)
s.2 <- rep(1,10000)

for(i in 1:10000)
{
	s.1[i] <- sum(sample(x)*w.1)
	s.2[i] <- sum(sample(x)*w.2)
}


par(mfcol=c(1,2))
hist(s.1) ; hist(s.2)

table(s.1) ; 
table(s.2)  ; table(round(s.2,3))

mean(s.1) ; mean(s.2) 
var(s.1)  ; var(s.2)

#  permutation test p-values 

mean(s.1 >= obs.1) ; mean(s.1 > obs.1)
mean(s.2 >= obs.2) ; mean(s.2 > obs.2)



#########################
#
#	Question 2:
#
#########################


exp.05 <- rexp(200,.05)
exp.01 <- rexp(200,.01)
exp.20 <- rexp(200,.20)

time.0 <- rexp(200,.05)
time.1 <- c(exp.05,ifelse(exp.05 < 20,exp.05,exp.20+20))
time.2 <- c(exp.05,ifelse(exp.01 < 20,exp.01,exp.05+20))

x    <- c(rep(1,200),rep(2,200))

c.0 <- rexp(200,0.01)

stat.1 <- time.1 <= c.0
stat.2 <- time.2 <= c.0

obs.1 <- ifelse(stat.1,time.1,c.0)
obs.2 <- ifelse(stat.2,time.2,c.0)

surv.1 <- Surv(obs.1,stat.1)
surv.2 <- Surv(obs.2,stat.2)

par(mfcol=c(1,2))
plot(survfit(surv.1~x),col=c(1,2))
plot(survfit(surv.2~x),col=c(1,2))

survdiff(surv.1~x,rho=1)
survdiff(surv.1~x,rho=-1)

survdiff(surv.2~x,rho=1)
survdiff(surv.2~x,rho=-1)

# cumulative hazard plots

ind.1 <- summary(survfit(surv.1~x))$strata == "x=1"
ind.2 <- summary(survfit(surv.1~x))$strata == "x=2"
y    <- summary(survfit(surv.1~x))$surv
time <- summary(survfit(surv.1~x))$time

par(mfcol=c(1,2))
plot(time[ind.1],-log(y[ind.1]),xlim=c(0,40),xlab="time",ylab="Cum hazard",type="l")
lines(time[ind.2],-log(y[ind.2]),col=2)

ind.1 <- summary(survfit(surv.2~x))$strata == "x=1"
ind.2 <- summary(survfit(surv.2~x))$strata == "x=2"
y    <- summary(survfit(surv.2~x))$surv
time <- summary(survfit(surv.2~x))$time

plot(time[ind.1],-log(y[ind.1]),xlim=c(0,40),xlab="time",ylab="Cum hazard",type="l")
lines(time[ind.2],-log(y[ind.2]),col=2)



#########################
#
#	Question 3:
#
#########################


# Stratified 2 group nonparametric


time   <- c(3,8,18,20,21,22,24,25,10,13,15,16,17,19,25,32)
status <- c(1,1,1,0,1,1,1,1,0,0,1,1,0,1,1,1)
group  <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2)
sex    <- c(1,1,2,2,2,2,2,2,1,1,1,1,1,1,2,2)


Surv(time[group==1 & sex==1],status[group==1 & sex==1])
Surv(time[group==2 & sex==1],status[group==2 & sex==1])
Surv(time[group==1 & sex==2],status[group==1 & sex==2])
Surv(time[group==2 & sex==2],status[group==2 & sex==2])

survdiff(Surv(time,status)~ group,rho=0)
survdiff(Surv(time,status)~ group + strata(sex),rho=0)


a.1    <- summary(survfit(Surv(time[sex==1],status[sex==1])))
time.1 <- a.1$time
n.1    <- a.1$n.risk
d.1    <- a.1$n.event
d.11   <- c(1,1,0,0,0)
n.11   <- c(2,1,0,0,0)
e.11   <- n.11*d.1/n.1
v.11   <- n.11*(n.1-n.11)*d.1*(n.1-d.1)/(n.1^2*(n.1-1))
(tbl.1 <- cbind(time.1,n.1,n.11,d.1,d.11,e.11,v.11))

a.2    <- summary(survfit(Surv(time[sex==2],status[sex==2])))
time.2 <- a.2$time
n.2    <- a.2$n.risk
d.2    <- a.2$n.event
d.21   <- c(1,1,1,1,1,0)
n.21   <- c(6,4,3,2,1,0)
e.21   <- n.21*d.2/n.2
v.21   <- n.21*(n.2-n.21)*d.2*(n.2-d.2)/(n.2^2*(n.2-1))
(tbl.2 <- cbind(time.2,n.2,n.21,d.2,d.21,e.21,v.21))


sum(c(tbl.1[,5]-tbl.1[,6],tbl.2[,5]-tbl.2[,6]))^2 / (sum(tbl.1[1:4,7])+sum(tbl.2[1:5,7]))