#########################################################################################
# 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	  #
###########################################

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)
	}

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

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

	return(likvals)

}

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

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

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)

likpars <- exp(-likpars*time)

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, cuttime) {

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

for (i in 1:dim(simdat)[3]) {
timdiff <- simdat[,4,i][simdat[,1,i]<=num] - min(simdat[,4,i][simdat[,1,i]<=num])
fu.timen <- timdiff + cuttime
eventn <- ifelse(fu.timen < simdat[,2,i][simdat[,1,i]<=num],0,1)
s.timen <- ifelse(fu.timen < simdat[,2,i][simdat[,1,i]<=num], fu.timen, simdat[,2,i][simdat[,1,i]<=num])

lookdat[,,i] <- cbind(simdat[,1,i][simdat[,1,i]<=num], s.timen, simdat[,3,i][simdat[,1,i]<=num], fu.timen, eventn)
}

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="Survival Probability at 6 Months", 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, xlab) {
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=xlab, 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 

# 50 data sets
simarray(40, 0.0594, 36, 50)
lmba <- seq(0.00001, 0.3, by = 0.00001) 


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

# 50 data sets - shorter accrual time
simarray(40, 0.0594, 24, 50)
lmba <- seq(0.00001, 0.3, 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.3, 1)), range(0:1), type="n", xlab="Survival Probability at 6 Months", ylab="Likelihood Function")
choices <- colors()
color = sample(choices,25, replace = F) 
for (i in 1:50) {
lines(exp(-1*lmba*6), exp(likmat40[,i]), col=color[i])
}

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

# 3 month cutoff #
simdat5 <- getlookdat(simdat, 5, 3)
simdat10 <- getlookdat(simdat, 10, 3)
simdat15 <- getlookdat(simdat, 15, 3)
simdat20 <- getlookdat(simdat, 20, 3)
simdat25 <- getlookdat(simdat, 25, 3)
simdat30 <- getlookdat(simdat, 30, 3)
simdat35 <- getlookdat(simdat, 35, 3)
simdat40 <- getlookdat(simdat, 40, 3)

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

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

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

likvals20 <- getlik(lmba, simdat20)
likmat20 <- likmat
likpars20 <- getints(likmat20, mle, lmba, 3)
likpars20[,2][likpars20[,1]==1] <- 0
likpars20[,3][likpars20[,1]==1] <- 1

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

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

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

likvals40 <- getlik(lmba, simdat40)
likmat40 <- likmat
likpars40 <- getints(likmat40, mle, lmba, 3)

pars5.3 <- likpars5
vals5.3 <- likvals5
pars10.3 <- likpars10
vals10.3 <- likvals10
pars15.3 <- likpars15
vals15.3 <- likvals15
pars20.3 <- likpars20
vals20.3 <- likvals20
pars25.3 <- likpars25
vals25.3 <- likvals25
pars30.3 <- likpars30
vals30.3 <- likvals30
pars35.3 <- likpars35
vals35.3 <- likvals35
pars40.3 <- likpars40
vals40.3 <- likvals40

# 6 month cutoff #

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

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

likvals10 <- getlik(lmba, simdat10)
likmat10 <- likmat
likpars10 <- getints(likmat10, mle, lmba, 6)
likpars10[,2][likpars10[,1]==1] <- 0
likpars10[,3][likpars10[,1]==1] <- 1

likvals15 <- getlik(lmba, simdat15)
likmat15 <- likmat
likpars15 <- getints(likmat15, mle, lmba, 6)
likpars15[,2][likpars15[,1]==1] <- 0
likpars15[,3][likpars15[,1]==1] <- 1

likvals20 <- getlik(lmba, simdat20)
likmat20 <- likmat
likpars20 <- getints(likmat20, mle, lmba,6)
likpars20[,2][likpars20[,1]==1] <- 0
likpars20[,3][likpars20[,1]==1] <- 1

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

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

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

likvals40 <- getlik(lmba, simdat40)
likmat40 <- likmat
likpars40 <- getints(likmat40, mle, lmba, 6)

pars5.6 <- likpars5
vals5.6 <- likvals5
pars10.6 <- likpars10
vals10.6 <- likvals10
pars15.6 <- likpars15
vals15.6 <- likvals15
pars20.6 <- likpars20
vals20.6 <- likvals20
pars25.6 <- likpars25
vals25.6 <- likvals25
pars30.6 <- likpars30
vals30.6 <- likvals30
pars35.6 <- likpars35
vals35.6 <- likvals35
pars40.6 <- likpars40
vals40.6 <- likvals40

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

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

likvals10 <- getlik(lmba, simdat10)
likmat10 <- likmat
likpars10 <- getints(likmat10, mle, lmba, 9)
likpars10[,2][likpars10[,1]==1] <- 0
likpars10[,3][likpars10[,1]==1] <- 1

likvals15 <- getlik(lmba, simdat15)
likmat15 <- likmat
likpars15 <- getints(likmat15, mle, lmba, 9)
likpars15[,2][likpars15[,1]==1] <- 0
likpars15[,3][likpars15[,1]==1] <- 1

likvals20 <- getlik(lmba, simdat20)
likmat20 <- likmat
likpars20 <- getints(likmat20, mle, lmba, 9)
likpars20[,2][likpars20[,1]==1] <- 0
likpars20[,3][likpars20[,1]==1] <- 1

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

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

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

likvals40 <- getlik(lmba, simdat40)
likmat40 <- likmat
likpars40 <- getints(likmat40, mle, lmba, 9)

pars5.9 <- likpars5
vals5.9 <- likvals5
pars10.9 <- likpars10
vals10.9 <- likvals10
pars15.9 <- likpars15
vals15.9 <- likvals15
pars20.9 <- likpars20
vals20.9 <- likvals20
pars25.9 <- likpars25
vals25.9 <- likvals25
pars30.9 <- likpars30
vals30.9 <- likvals30
pars35.9 <- likpars35
vals35.9 <- likvals35
pars40.9 <- likpars40
vals40.9 <- likvals40

# Plot one data set first #
plot(range(c(likpars5[1,], likpars10[1,], likpars15[1,])), range(0:9), yaxt="n", type="n", xlab="Survival Probability at 6 Months", 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(likpars40[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=exp(-0.0594*6), lty=2, col="grey")

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

all.looks3 <- array(c(pars5.3, pars10.3, pars15.3, pars20.3, pars40.3), dim=c(50,3,5))
all.looks6 <- array(c(pars5.6, pars10.6, pars15.6, pars20.6, pars40.6), dim=c(50,3,5))
all.looks9 <- array(c(pars5.9, pars10.9, pars15.9, pars20.9, pars40.9), dim=c(50,3,5))

plotbounds(all.looks3, c("white", "red", "blue", "yellow", "green"), exp(-0.0594*3), "Survival Probability at 3 months")
legend(-0.05, 12, col=c("black", "white", "red", "blue", "yellow", "green"), lty=rep(1, 6), c("Estimate MLE after:", "5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), bty="n", cex=.85, fill="black") 

plotbounds(all.looks6, c("white", "red", "blue", "yellow", "green"), exp(-0.0594*6), "Survival Probability at 6 months")
legend(-0.05, 12, col=c("black", "white", "red", "blue", "yellow", "green"), lty=rep(1, 6), c("Estimate MLE after:", "5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), bty="n", cex=.85, fill="black") 

plotbounds(all.looks9, c("white", "red", "blue", "yellow", "green"), exp(-0.0594*9), "Survival Probability at 9 months")
#legend(-0.05, 12, col=c("black", "white", "red", "blue", "yellow", "green"), lty=rep(1, 6), c("Estimate MLE after:", "5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), bty="n", cex=.85, fill="black") 


# Calculate likelihood ratios for survival probabilities = 0.70 and 0.45 #
# N = 40 has 90% power to detect 70% v 45% 6 month survival #
# N = 40 has 90% power to detect 84% v 62% 3 month survival #
# N = 40 has 90% power to detect 59% v 34% 9 month survival #

# 0.70 => lambda = 0.0594 => lmba value 5940
# 0.45 => lambda = 0.13308 => lmba value 13308

# 0.84 => lambda = 0.836775 => lmba value 5940
# 0.62 => lambda = 0.6281351 => lmba 15500

# 0.59 => lambda = 0.5859036 => lmba value 5940
# 0.34 => lambda = NOT IN THE VALUES OF LAMBDA. ARGH!

lr5.6 <- lr10.6 <- lr15.6 <- lr20.6 <- lr25.6 <- lr30.6 <- lr35.6 <- lr40.6 <- rep(NA, 50)
lr5.3 <- lr10.3 <- lr15.3 <- lr20.3 <- lr25.3 <- lr30.3 <- lr35.3 <- lr40.3 <- rep(NA, 50)
lr5.9 <- lr10.9 <- lr15.9 <- lr20.9 <- lr25.9 <- lr30.9 <- lr35.9 <- lr40.9 <- rep(NA, 50)

for (i in 1:50) {
lr5.6[i] <- exp(vals5.6[5940,i])/exp(vals5.6[13308,i])
lr10.6[i] <- exp(vals10.6[5940,i])/exp(vals10.6[13308,i])
lr15.6[i] <- exp(vals15.6[5940,i])/exp(vals15.6[13308,i])
lr20.6[i] <- exp(vals20.6[5940,i])/exp(vals20.6[13308,i])
lr25.6[i] <- exp(vals25.6[5940,i])/exp(vals25.6[13308,i])
lr30.6[i] <- exp(vals30.6[5940,i])/exp(vals30.6[13308,i])
lr35.6[i] <- exp(vals35.6[5940,i])/exp(vals35.6[13308,i])
lr40.6[i] <- exp(vals40.6[5940,i])/exp(vals40.6[13308,i])

lr5.9[i] <- exp(vals5.9[5940,i])/exp(vals5.9[13308,i])
lr10.9[i] <- exp(vals10.9[5940,i])/exp(vals10.9[13308,i])
lr15.9[i] <- exp(vals15.9[5940,i])/exp(vals15.9[13308,i])
lr20.9[i] <- exp(vals20.9[5940,i])/exp(vals20.9[13308,i])
lr25.9[i] <- exp(vals25.9[5940,i])/exp(vals25.9[13308,i])
lr30.9[i] <- exp(vals30.9[5940,i])/exp(vals30.9[13308,i])
lr35.9[i] <- exp(vals35.9[5940,i])/exp(vals35.9[13308,i])
lr40.9[i] <- exp(vals40.9[5940,i])/exp(vals40.9[13308,i])

lr5.3[i] <- exp(vals5.3[5940,i])/exp(vals5.3[15500,i])
lr10.3[i] <- exp(vals10.3[5940,i])/exp(vals10.3[15500,i])
lr15.3[i] <- exp(vals15.3[5940,i])/exp(vals15.3[15500,i])
lr20.3[i] <- exp(vals20.3[5940,i])/exp(vals20.3[15500,i])
lr25.3[i] <- exp(vals25.3[5940,i])/exp(vals25.3[15500,i])
lr30.3[i] <- exp(vals30.3[5940,i])/exp(vals30.3[15500,i])
lr35.3[i] <- exp(vals35.3[5940,i])/exp(vals35.3[15500,i])
lr40.3[i] <- exp(vals40.3[5940,i])/exp(vals40.3[15500,i])

}

# Scatterplot of LR for each data set #

# 6 month intervals #
par(bg="black", fg="white", col.axis = "white", col.lab="white")
plot(range(1:50), range(log(lr25.6)), xlab="Simulated Data Set", ylab="Log Scale Likelihood Ratio", ylim=c(0,35))
points(1:50, log(lr5.6), pch=19, col="white") 
abline(h = log(32), lty=2, col="white")
points(1:50, log(lr10.6), pch=19, col="red") 
points(1:50, log(lr15.6), pch=19, col="blue") 
points(1:50, log(lr20.6), pch=19, col="yellow") 
points(1:50, log(lr40.6), pch=19, col="green")
legend(0,20,c("5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), pch=rep(19,5), col=c("white", "red", "blue", "yellow", "green"), cex=.7, adj=0, x.intersp=.5, text.width=2.75)

# 3 month intervals #
par(bg="black", fg="white", col.axis = "white", col.lab="white")
plot(range(1:50), range(log(lr25.3)), xlab="Simulated Data Set", ylab="Log Scale Likelihood Ratio", ylim=c(0,35))
points(1:50, log(lr5.3), pch=19, col="white") 
abline(h = log(32), lty=2, col="white")
points(1:50, log(lr10.3), pch=19, col="red") 
points(1:50, log(lr15.3), pch=19, col="blue") 
points(1:50, log(lr20.3), pch=19, col="yellow") 
points(1:50, log(lr40.3), pch=19, col="green")
legend(0,30,c("5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), pch=rep(19,5), col=c("white", "red", "blue", "yellow", "green"), cex=.7, adj=0, x.intersp=.5, text.width=2.75)

# 9 month intervals #
par(bg="black", fg="white", col.axis = "white", col.lab="white")
plot(range(1:50), range(log(lr25.9)), xlab="Simulated Data Set", ylab="Log Scale Likelihood Ratio", ylim=c(0,20))
points(1:50, log(lr5.9), pch=19, col="white") 
abline(h = log(32), lty=2, col="white")
points(1:50, log(lr10.9), pch=19, col="red") 
points(1:50, log(lr15.9), pch=19, col="blue") 
points(1:50, log(lr20.9), pch=19, col="yellow") 
points(1:50, log(lr40.9), pch=19, col="green")
legend(0,20,c("5 patients", "10 patients", "15 patients", "20 patients", "40 patients"), pch=rep(19,5), col=c("white", "red", "blue", "yellow", "green"), cex=.7, adj=0, x.intersp=.5, text.width=2.75)

boxplot(lr5.6, lr10.6, lr15.6, lr20.6, lr40.6, names=c("N=5", "N=10", "N=15", "N=20", "N=40"), ylab="Log Scale Likelihood Ratio", ylim=c(0, 100))


# Use Freq. approach and calculate CI for each data set, using data once last patient reaches cutoff time #

bound6 <- bound3 <- bound9 <- matrix(1000, nrow=50, ncol=3)

# Cutoff time = 6 #

freqdat <- getlookdat(simdat, 40, 6)

for (i in 1:50) {

if(is.element(6, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound6[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 6]
bound6[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 6]
}

if(is.element(5, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound6[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 5]
bound6[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 5]
}

if(bound6[i,2]==1000) {
bound6[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 4]
bound6[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 4]
}

bound6[i,1] <- 0.5

}


freq.look6 <- array(c(pars40.6, bound6), dim=c(50,3,2))
plotbounds(freq.look6, c("red", "blue"), exp(-0.0594*6), "Survival Probability at 6 months")
legend(0.365, 8, col=c("black", "red", "blue"), lty=rep(1, 3), c("Estimation Method:", "Likelihood", "Frequentist (Kaplan-Meier)"), bty="n", cex=.85, fill="black") 

# Cutoff time = 3 #

freqdat <- getlookdat(simdat, 40, 3)

for (i in 1:50) {

if(is.element(3, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound3[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 3]
bound3[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 3]
}

if(is.element(2, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound3[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 2]
bound3[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 2]
}

if(bound3[i,2]==1000) {
bound3[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 1]
bound3[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 1]
}

bound3[i,1] <- 0.5

}

freq.look3 <- array(c(pars40.3, bound3), dim=c(50,3,2))
plotbounds(freq.look3, c("red", "blue"), exp(-0.0594*3), "Survival Probability at 3 months")
legend(0.62, 6, col=c("black", "red", "blue"), lty=rep(1, 3), c("Estimation Method:", "Likelihood", "Frequentist (Kaplan-Meier)"), bty="n", cex=.85, fill="black") 


# Cutoff time = 9 #

freqdat <- getlookdat(simdat, 40, 9)

for (i in 1:50) {

if(is.element(9, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound9[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 9]
bound9[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 9]
}

if(is.element(8, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound9[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 8]
bound9[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 8]
}

if(is.element(7, summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time) == "TRUE") {
bound9[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 7]
bound9[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 7]
}

if(bound9[i,2]==1000) {
bound9[i,2] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$lower[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 6]
bound9[i,3] <- summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$upper[summary(survfit(Surv(freqdat[,2,i], freqdat[,5,i])))$time == 6]
}


bound9[i,1] <- 0.5

}

freq.look9 <- array(c(pars40.9, bound9), dim=c(50,3,2))
plotbounds(freq.look9, c("red", "blue"), exp(-0.0594*9), "Survival Probability at 9 months")
legend(0.32, 6, col=c("black", "red", "blue"), lty=rep(1, 3), c("Estimation Method:", "Likelihood", "Frequentist (Kaplan-Meier)"), bty="n", cex=.85, fill="black")