#########################################################################################
# R code to simulate survival data and plot maximum likelihood estimates for lamba 	#
# 25 October 2005									#
# Cancer Clinical Trials Working Group							#
#########################################################################################

######################################
# Function to simulate survival data #
# n = sample size		     #
# lambda = hazard rate 		     #
# atime = accrual time (months)      #
# nsim = number of data sets	     #
######################################


simarray <- function(n, lambda, atime, nsim)  {

	simdat <- array(NA, dim=c(n, 5, nsim)) 

		for (i in 1:nsim) {
		s.time <- round(rexp(n,lambda))
		accrual.time <- round(runif(n,0,atime))
		fu.time <- 60-accrual.time
		event <- ifelse(fu.time < s.time,0,1)
		s.time <- ifelse(fu.time < s.time, fu.time, s.time)
		id <- as.integer(rank(accrual.time, ties.method = "first"))
		simdat[,,i] <- cbind(id, s.time, accrual.time, fu.time, event)
		}

	assign("simdat", simdat, env=.GlobalEnv)

}

###########################################
# Function to calculate the MLE of lambda #
# lmba: sequence of values for lambda	  #
###########################################

# lmba <- seq(0.0001, 0.08, by = 0.00001) 

getlik <- function(lmba, simdat) {

	likvals <- matrix(NA, ncol=dim(simdat)[3], nrow=length(lmba))
 
	likmat <- matrix(NA, ncol=dim(simdat)[3], nrow = length(lmba))
	mle <- likatmle <- rep(NA, dim(simdat)[3])

	for (i in 1:dim(simdat)[3]) {
	mle[i] <- 1/(sum(simdat[,5,i]*(simdat[,2,i]) + (1 - simdat[,5,i])*simdat[,4,i]))*sum(simdat[,5,i])
	likatmle[i] <- round(sum(-1*(simdat[,5,i]*(log(1/mle[i]) + (simdat[,2,i])*mle[i]) + (1 - simdat[,5,i])*mle[i]*simdat[,4,i])), 2)

	for (j in 1:length(lmba)) {
	likmat[j,i] <- round(sum(-1*(simdat[,5,i]*(log(1/lmba[j]) + (simdat[,2,i])*lmba[j]) + (1 - simdat[,5,i])*lmba[j]*simdat[,4,i])), 2)
	}

	likmat[,i] <- likmat[,i] - likatmle[i]
	}

	assign("likmat", likmat, env=.GlobalEnv)
	assign("mle", mle, env=.GlobalEnv)

}

################################################################
# Get 1/8 support intervals - point where lik function = 0.125 #
################################################################

getints <- function(likmat, mle, lmba) {

low <- upp <- rep(NA, length(mle))

for (i in 1:length(mle)) {

low[i] <- min(lmba[exp(likmat[,i]) > 0.12 & exp(likmat[,i]) < 0.13])
upp[i] <- max(lmba[exp(likmat[,i]) > 0.12 & exp(likmat[,i]) < 0.13])

}

likpars <- cbind(mle, low, upp)

assign("likpars", likpars, env=.GlobalEnv)

}


###################################################
# Format data for the number of people to look at #
# Uses data through end of follow-up		  #
# num: Number of patients to evaluate		  #
###################################################

getlookdat <- function(simdat, num) {

lookdat <- array(NA, dim=c(num, dim(simdat)[2], dim(simdat)[3]))

for (i in 1:dim(simdat)[3]) {
lookdat[,,i] <- matrix(simdat[,,i][simdat[,1,i]<=num], nrow=num, byrow=F)
}

return(lookdat)

}

##############################################################
# Plot mle and support intervals for each simulated data set #
# The first function (plotmle) plots data at the end of trial#
# The 2nd function (plotbounds) plots data at each look	     #
# all.looks: array of the data for every 5 pts. - see below  #
##############################################################

plotmle <- function(likpars, color, lsim) {
plot(range(likpars), range(0:dim(likpars)[1]), type="n", xlab=expression(paste(lambda)), ylab="Simulated Data Set")
points(likpars[,1], 1:dim(likpars)[1], pch = 19)
segments(likpars[,2], 1:dim(likpars)[1], likpars[,3], 1:dim(likpars)[1], col=color)
abline(v = lsim, lty=2, col="grey")
}

plotbounds <- function(all.looks, color, lsim) {
par(bg="black", fg="white", col.axis = "white", col.lab="white")
plot(range(all.looks[,2:3,]), range(0:dim(all.looks)[1]), type="n", xlab=expression(paste(lambda)), ylab="Simulated Data Set")
for (i in 1:dim(all.looks)[3]) {
segments(all.looks[,2,i], 1:dim(all.looks)[1], all.looks[,3,i], 1:dim(all.looks)[1], col=color[i])
}
abline(v = lsim, lty=2, col="grey")
}


#####################################################################################################################################

# Try an example 

# 100 data sets
simarray(40, 1/50, 36, 100)
lmba <- seq(0.0001, 0.08, by = 0.00001) 
getlik(lmba, simdat)
getints(likmat, mle, lmba)
plotmle(likpars, "black", 1/100)

# 50 data sets
simarray(40, 1/100, 36, 50)
lmba <- seq(0.0001, 0.08, by = 0.00001) 
getlik(lmba, simdat)
getints(likmat, mle, lmba)
plotmle(likpars, "black", 1/100)

# 50 data sets - longer accrual time
simarray(40, 1/100, 48, 50)
lmba <- seq(0.0001, 0.08, by = 0.00001) 
getlik(lmba, simdat)
getints(likmat, mle, lmba)
plotmle(likpars, "black", 1/100)

# 50 data sets - shorter accrual time
simarray(40, 1/100, 24, 50)
lmba <- seq(0.0001, 0.08, by = 0.00001) 
getlik(lmba, simdat)
getints(likmat, mle, lmba)
plotmle(likpars, "black", 1/100)

# 50 data sets - N=20
simarray(20, 1/100, 36, 50)
lmba <- seq(0.0001, 0.08, by = 0.00001) 
getlik(lmba, simdat)
getints(likmat, mle, lmba)
plotmle(likpars, "black", 1/100)

# Plot likelihood functions for the data sets #

plot(range(c(0, 0.03)), range(0:1), type="n", xlab=expression(paste(lambda)), ylab="Likelihood Function")
choices <- colors()
color = sample(choices,25, replace = F) 
for (i in 1:50) {
lines(lmba, exp(likmat[,i]), col=color[i])
}

# For each of 50 data sets, start looking at intervals of 5 people

simdat5 <- getlookdat(simdat, 5)
simdat10 <- getlookdat(simdat, 10)
simdat15 <- getlookdat(simdat, 15)
simdat20 <- getlookdat(simdat, 20)
simdat25 <- getlookdat(simdat, 25)
simdat30 <- getlookdat(simdat, 30)
simdat35 <- getlookdat(simdat, 35)

getlik(lmba, simdat5)
likpars5 <- getints(likmat, mle, lmba)
likpars5[,2][likpars5[,1]==0] <- min(lmba) # Replaces interval to cover whole line in case of no events #
likpars5[,3][likpars5[,1]==0] <- max(lmba)

getlik(lmba, simdat10)
likpars10 <- getints(likmat, mle, lmba)
likpars10[,2][likpars10[,1]==0] <- min(lmba)
likpars10[,3][likpars10[,1]==0] <- max(lmba)

getlik(lmba, simdat15)
likpars15 <- getints(likmat, mle, lmba)
likpars15[,2][likpars15[,1]==0] <- min(lmba)
likpars15[,3][likpars15[,1]==0] <- max(lmba)

getlik(lmba, simdat20)
likpars20 <- getints(likmat, mle, lmba)

getlik(lmba, simdat25)
likpars25 <- getints(likmat, mle, lmba)

getlik(lmba, simdat30)
likpars30 <- getints(likmat, mle, lmba)

getlik(lmba, simdat35)
likpars35 <- getints(likmat, mle, lmba)

getlik(lmba, simdat)

likpars <- getints(likmat, mle, lmba)

all.looks <- array(c(likpars5, likpars10, likpars15, likpars20), dim=c(50,3,4))

# Plot one data set first #
plot(range(c(likpars5[1,], likpars10[1,], likpars15[1,])), range(0:9), yaxt="n", type="n", xlab=expression(paste(lambda)), ylab="Number of patients used to estimate MLE")
points(likpars5[1,1], 1, pch=19)
points(likpars10[1,1], 2, pch=19)
points(likpars15[1,1], 3, pch=19)
points(likpars20[1,1], 4, pch=19)
points(likpars25[1,1], 5, pch=19)
points(likpars30[1,1], 6, pch=19)
points(likpars35[1,1], 7, pch=19)
points(likpars[1,1], 8, pch=19)

segments(likpars5[1,2], 1, likpars5[1,3], 1)
segments(likpars10[1,2], 2, likpars10[1,3], 2)
segments(likpars15[1,2], 3, likpars15[1,3], 3)
segments(likpars20[1,2], 4, likpars20[1,3], 4)
segments(likpars25[1,2], 5, likpars25[1,3], 5)
segments(likpars30[1,2], 6, likpars30[1,3], 6)
segments(likpars35[1,2], 7, likpars35[1,3], 7)
segments(likpars[1,2], 8, likpars[1,3], 8)
axis(2, at=c(1:8), c(5,10,15,20,25,30,35,40))
abline(v=0.01, lty=2, col="grey")


# Use plotbounds function to overlay interval estimates at every 5 patients

pdf("Support bounds plot 50 sims.pdf", height = 8, width=11)
plotbounds(all.looks, c("white", "red", "aquamarine", "yellow", "green"), 1/100)
legend(0.06, 13, col=c("black", "white", "red", "aquamarine", "yellow", "green"), lty=rep(1, 6), c("Estimate MLE after:", "5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), bty="n", cex=.85) 
dev.off()

plotbounds(all.looks, c("white", "red", "aquamarine", "yellow"), 1/100)
legend(0.06, 35, col=c("black", "white", "red", "aquamarine", "yellow"), lty=rep(1, 5), c("Estimate MLE after:", "5 patients", "10 patients", "15 patients", "20 patients"), bty="n", cex=.85)