#Version: 1.0  Build: 08/01/2006
#Version: 1.1  Build: 09/23/2008
# Note: Run functions file before running the program

# To begin, create options file and write it's location in below
#options.file="C:/RSource/TimeBMD/toxicoTet.(D)"
#Created Using R2.3.1

############################################################################
## The program consists of four parts                                     ## 
## Part 1:Exploratory Data Analysis                                       ##
## Part 2:Dose-Response Modeling                                          ##
## Part 3:Estimate of Benchmark Doses(BMD)                                ##
## Part 4:Bootstrap BMDL and figures                                      ##
############################################################################

#########################################################
#  Part 1:  Exloratory data analysis                    #
#########################################################

###Include all the necessary libraries
library(base)
library(nlme)
library(lattice)

#Change the working directory to option specified in options file (or leave as default)
#if(scan(options.file,skip=22,n=1,what=character(),allowEscapes=F, quiet=T)!=-9999) {
#	setwd(scan(options.file,skip=22,n=1,what=character(),allowEscapes=F, quiet=T))
#	}
setwd(rVars$info$DataDir)
if(rVars$info$DebugFile!="")
capture.output(cat("Set Work Dir (DataDir):",rVars$info$DataDir,"\n"), file = rVars$info$DebugFile, append = FALSE)

# loads dataset from information given in input file and call it dataset.user 
#dataset.user <- read.table(scan(options.file,skip=2,n=1,what=character(),allowEscapes=F, quiet=T),header=T)
dataset.user <- read.table(rVars$info$DataText,header=T)

# create data set of relevent variable (dataset) and standardize column headers
for (i in 1:length(dataset.user[1,])) {
#	if(names(dataset.user)[i]==scan(options.file,skip=4,n=4,what=character(), quiet=T)[1]) i.ID=i
#	if(names(dataset.user)[i]==scan(options.file,skip=4,n=4,what=character(), quiet=T)[2]) i.dose=i
#	if(names(dataset.user)[i]==scan(options.file,skip=4,n=4,what=character(), quiet=T)[3]) i.time=i
#	if(names(dataset.user)[i]==scan(options.file,skip=4,n=4,what=character(), quiet=T)[4]) i.response=i
	if(names(dataset.user)[i]==rVars$info$Header.ANID) i.ID=i
	if(names(dataset.user)[i]==rVars$info$Header.DOSE) i.dose=i
	if(names(dataset.user)[i]==rVars$info$Header.TIME) i.time=i
	if(names(dataset.user)[i]==rVars$info$Header.RESP) i.response=i
}
dataset <-  cbind(ID=dataset.user[i.ID],dose=dataset.user[i.dose],time=dataset.user[i.time],response=dataset.user[i.response])
rm(i.ID,i.dose,i.time,i.response)

# removes observations with missing responses (-9999) from dataset
j=1;dataset.tmp=dataset[1,] #initialize values
for (i in 1:length(dataset[,1])) {
	if(dataset[i,4]!=-9999) {
		dataset.tmp[j,]<-dataset[i,]
		j=j+1
	}
}
#match(tmp$ID, sort(tmp$ID))
dataset <- data.frame(ID=match(dataset.tmp[,1],sort(dataset.tmp[,1])), dose=dataset.tmp[,2], time=dataset.tmp[,3], response=dataset.tmp[,4])
rm(j,dataset.tmp)

# stores the original time values (time.orig) and transform the time values such that the minimum time is at 0.
time.orig<-dataset$time
min.time.orig<-min(time.orig)
dataset$time<-time.orig-min.time.orig
time.level.orig<-as.numeric(levels(as.factor(time.orig)))

attach(dataset, warn.conflicts=FALSE)

# read names of variables and units of measurement from options file
#dose.name <- scan(options.file,skip=5,n=4,what=character(), quiet=T)[2]
#	if(dose.name==-9999) dose.name<-scan(options.file,skip=4,n=4,what=character(), quiet=T)[2]
dose.name=rVars$info$Header.DOSE
if(rVars$info$DebugFile!="")
capture.output(cat("dose.name: ",dose.name,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#dose.units <- scan(options.file,skip=6,n=4,what=character(), quiet=T)[2]
#dose.units="mg/kg"
dose.units=-9999

#time.name <- scan(options.file,skip=5,n=4,what=character(), quiet=T)[3]
#	if(time.name==-9999) time.name<-scan(options.file,skip=4,n=4,what=character(), quiet=T)[3]
time.name=rVars$info$Header.TIME
if(rVars$info$DebugFile!="")
capture.output(cat("time.name: ",time.name,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#time.units <- scan(options.file,skip=6,n=4,what=character(), quiet=T)[3]
#time.units="week"
time.units=-9999

#response.name <- scan(options.file,skip=5,n=4,what=character(), quiet=T)[4]
#	if(response.name==-9999) response.name<-scan(options.file,skip=4,n=4,what=character(), quiet=T)[4]
response.name=rVars$info$Header.RESP
if(rVars$info$DebugFile!="")
capture.output(cat("response.name: ",response.name,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#response.units <- scan(options.file,skip=6,n=4,what=character(), quiet=T)[4]
#response.units="kg"
response.units=-9999

# other options read from file

# some charting options
#logscale=scan(options.file,skip=7,n=1,what=character(), quiet=T)[1]
#	if(logscale==-9999) logscale=FALSE; logscale<-as.logical(logscale)
logscale<-rVars$info$TALS=="Log"
if(rVars$info$DebugFile!="")
capture.output(cat("logscale: ",logscale,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#chart1.title <- scan(options.file,skip=8,n=1,what=character(), quiet=T)[1]
#	if(chart1.title==-9999) chart1.title<-paste(response.name,"vs.",time.name, "by", dose.name)   #default
chart1.title=rVars$info$CTIT
if(rVars$info$CTIT=="na") chart1.title=paste(rVars$info$Header.RESP,"vs.",rVars$info$Header.TIME, "by", dose.name)
if(rVars$info$DebugFile!="")
capture.output(cat("chart1.title: ",chart1.title,"\n"), file = rVars$info$DebugFile, append = TRUE)

# detects the different factor levels of the time variable
time.level<-as.numeric(levels(as.factor(time)))

# Creating a term ("logfix") to adjust log graph (to allow a 0 time intercept)
# Since time=0, we will plot log(time+logfix) instead of log(time)
     n  <- length(time.level)
     t1 <- time.level[1]
     t2 <- time.level[2]
     t3 <- time.level[3]
     logfix=t2/2 #initialize logfix= to the midpoint between the first and second time level (which is t2/2 since t1=0 always)
     if (2*t2-t3-t1!=0) logfix <- (-t2^2 + t1*t3)/(2*t2-t3-t1)
     rm(n,t1,t2,t3)


# creates scatterplot of original data
windows(rescale = "fit", restoreConsole = FALSE)
if (logscale) {
	observed_t <- xyplot(response~I(time+logfix) | as.factor(dose),data=dataset, groups=ID, panel=subj.panel, 
       scales=list(x=list(log=10, at=time.level+logfix, labels=c(time.level+min.time.orig))), 
       xlab=paste(time.name,if(time.units!=-9999)paste("(",time.units, "in log-scale )")), 
       ylab=paste(response.name,if(response.units!=-9999)paste("(",response.units, ")")), 
       main=chart1.title,panel.groups=paste(levels(as.factor(dose)),if(dose.units!=-9999)dose.units),
       layout=c(length(levels(as.factor(dose))),1,1),as.table=T, par.settings = list(par.main.text = list(cex = 0.9))) }   

if (!logscale) {
	observed_t <- xyplot(response~time | as.factor(dose),data=dataset, groups=ID, panel=subj.panel, 
	scales=list(x=list(at=time.level, labels=c(time.level+min.time.orig))),
       xlab=paste(time.name,if(time.units!=-9999)paste("(",time.units,")")), 
       ylab=paste(response.name,if(response.units!=-9999)paste("(",response.units, ")")), 
       main=chart1.title,panel.groups=paste(levels(as.factor(dose)),if(dose.units!=-9999)dose.units),
       layout=c(length(levels(as.factor(dose))),1,1),as.table=T,par.settings = list(par.main.text = list(cex = 0.9))) }   

# saves plot of original data
plot(observed_t)
savePlot(rVars$info$Observed_Trajectory,type="emf")

####################################################################
# Part 2: Dose-response modeling---ToxicoDiffusion model           #                          
####################################################################

# Use NLME in S-plus to fit response by toxico-diffusion model
# and obtain residuals and estimates of random effects
# Variance is assumed to be constant across dose levels

# Attached data set is already a grouped data frame; otherwise the following can be used to ceate a groupedDataframe:

dataset<-groupedData(response ~ time | ID, data=dataset, order.groups=FALSE,
		outer= ~dose,
	   labels=list(x=time.name, y=response.name),
	   units=list(x=time.units, y=response.units))
		
# define a vector of starting values (either from option file or by default)


# Read itial value function options from file
#ExpoTime<-scan(options.file,skip=23,n=1,what=numeric(), quiet=T)
#	if(ExpoTime==-9999) ExpoTime=0
ExpoTime=rVars$info$EXPT
if(rVars$info$DebugFile!="")
capture.output(cat("ExpoTime: ",ExpoTime,"\n"), file = rVars$info$DebugFile, append = TRUE)

#baselinePoly<-scan(options.file,skip=24,n=1,what=numeric(), quiet=T)+1
#	if(baselinePoly==-9998) baselinePoly=1
baselinePoly=rVars$info$BGDEG+1
if(rVars$info$DebugFile!="")
capture.output(cat("baselinePoly: ",baselinePoly,"\n"), file = rVars$info$DebugFile, append = TRUE)

# Read starting values from file
initialvalues<-rep(-9999,baselinePoly+3)

#A0<-scan(options.file,skip=20,n=1,what=numeric(), quiet=T)[1] #used to test if initial values are given in options file
A0 = rVars$info$SVA0

# From Geoffrey, needs some work
if(A0!=-9999) {
	initialvalues[1]=A0
#if (baselinePoly==3) initialvalues<-scan(options.file,skip=20,n=6,what=numeric(), quiet=T)
#if (baselinePoly==2) initialvalues<-scan(options.file,skip=20,n=5,what=numeric(), quiet=T)
#if (baselinePoly==1) initialvalues<-scan(options.file,skip=20,n=4,what=numeric(), quiet=T)
}
if(baselinePoly==3){
	if(rVars$info$SVA1!=-9999) initialvalues[2]=rVars$info$SVA1
	if(rVars$info$SVA2!=-9999) initialvalues[3]=rVars$info$SVA2
	if(rVars$info$SVB0!=-9999) initialvalues[4]=rVars$info$SVB0
	if(rVars$info$SVC0!=-9999) initialvalues[5]=rVars$info$SVC0
	if(rVars$info$SVK0!=-9999) initialvalues[6]=rVars$info$SVK0
}
if(baselinePoly==2){
	if(rVars$info$SVA1!=-9999) initialvalues[2]=rVars$info$SVA1
	if(rVars$info$SVB0!=-9999) initialvalues[3]=rVars$info$SVB0
	if(rVars$info$SVC0!=-9999) initialvalues[4]=rVars$info$SVC0
	if(rVars$info$SVK0!=-9999) initialvalues[5]=rVars$info$SVK0
}
if(baselinePoly==1){
	if(rVars$info$SVB0!=-9999) initialvalues[2]=rVars$info$SVB0
	if(rVars$info$SVC0!=-9999) initialvalues[3]=rVars$info$SVC0
	if(rVars$info$SVK0!=-9999) initialvalues[4]=rVars$info$SVK0
}
if(rVars$info$DebugFile!="")
capture.output(cat("initialvalues: ",initialvalues,"\n"), file = rVars$info$DebugFile, append = TRUE)

datasetm<--9999 		# initialization of object in which model will be stored
# use initial values provided by user and test if they work

dataset.expo<-data.frame(c(dataset,expo=ExpoTime))
attach(dataset.expo, warn.conflicts=FALSE)

#Problem1 here, GLN
if(A0!=-9999) {
	datasetm<-nlme(model= response~toxico.diffusion(A,B,C, K, time=time, ExpoTime=expo),
					 fixed= c(c(A~1,A~time,A~time+I(time^2))[baselinePoly], B~dose-1, C~dose-1, K~1),
					 random= A~1,
					 start= initialvalues, groups = ~ ID,
					 data= dataset.expo, na.action=na.omit, verbose=F)
	datasetm$initialValues<-initialvalues
	possible.init<-matrix(initialvalues,ncol=(3+baselinePoly))
}
if(rVars$info$DebugFile!="")
capture.output(datasetm, file = rVars$info$DebugFile, append = TRUE)

#End of Problem1 here, GLN

#Problem2 here, GLN
if((as.character(datasetm)==-9999)[1]) {
	#shortcut function to use nlme function with desired options
	d<-function(i) {datasetm<-try(nlme(model= response~toxico.diffusion(A,B,C, K, time, expo),
				fixed= c(c(A~1,A~time,A~time+I(time^2))[baselinePoly], B~dose-1, C~dose-1, K~1), random=A~1,
				start= c(fixed=matrix(possible.init,ncol=3+baselinePoly)[i,]), groups = ~ ID,
				data= dataset.expo, na.action=na.omit, verbose=F),silent=T)
				if (is.na(datasetm[3])) datasetm<--9999 #makes datasetm=-9999 if attempt failed
				datasetm}
	possible.init<-ToxicoDiffInitMatrix(dataset, ExpoTime, baselinePoly)
	possible.init.n<-length(possible.init[,1])
	for(i in 1:possible.init.n) if(as.character(datasetm)[1]==-9999) {
		datasetm<-d(i)
		if(baselinePoly==1) datasetm$initialValues<-c(A0=possible.init[i,1],B0=possible.init[i,2],
			C0=possible.init[i,3],K0=possible.init[i,4])
		if(baselinePoly==2) datasetm$initialValues<-c(A0=possible.init[i,1],A1=possible.init[i,2],
			B0=possible.init[i,3],C0=possible.init[i,4],K0=possible.init[i,5])
		if(baselinePoly==3) datasetm$initialValues<-c(A0=possible.init[i,1],A1=possible.init[i,2],
			A2=possible.init[i,3],B0=possible.init[i,4],C0=possible.init[i,5],
			K0=possible.init[i,6],)

	}
}
#End of Problem2 here, GLN	

# Plot the fitted panel graph with log-scaled time of the original scale 
# First obtain predicted response

dataset$Fitted[!is.na(response)]<-predict(datasetm)    # Retain NA for missing responses

if (logscale) {
	observed_fit <- xyplot(dataset$Fitted~I(time+logfix) | as.factor(dose),
	 data=dataset, groups=ID, panel=subj.panel, 
       scales=list(x=list(log=10, at=time.level+logfix, labels=c(time.level+min.time.orig))), 
       xlab=paste(time.name,if(time.units!=-9999)paste( "(",time.units, "in log-scale )")), 
       ylab=paste(response.name,if(response.units!=-9999)paste( "(",response.units, ")")),
       main = paste("Fitted Values of",response.name,"vs.",time.name,"by",dose.name), 
       strip=function(...) strip.default(...,strip.names = c(F, T), style = 1),
       layout=c(length(levels(as.factor(dose))),1,1),as.table=T, par.settings = list(par.main.text = list(cex = 0.9)))}

if (!logscale) {
	observed_fit <- xyplot(dataset$Fitted~time | as.factor(dose), 
	data=dataset, groups=ID, panel=subj.panel, 
	scales=list(x=list(at=time.level, labels=c(time.level+min.time.orig))),
       xlab=paste(time.name,if(time.units!=-9999)paste( "(",time.units, ")")), 
       ylab=paste(response.name,if(response.units!=-9999)paste( "(",response.units, ")")),
       main = paste("Fitted Values of",response.name,"vs.",time.name,"by",dose.name),
       strip=function(...) strip.default(...,strip.names = c(F, T), style = 1),
       layout=c(length(levels(as.factor(dose))),1,1),as.table=T, par.settings = list(par.main.text = list(cex = 0.9)))}

# saves graph
plot(observed_fit)
savePlot(rVars$info$Fitted_Trajectory,type="emf")


# Plot the residuals

# Computes standarzied residuals 
dataset$resid.sd[!is.na(response)]<-resid(datasetm)/datasetm$sigma

# plots residuals pooled across groups 
par(cex.main=.9) #makes the header slightly smaller
plot(dataset$Fitted,dataset$resid.sd, xlab=paste("Fitted values",if(response.units!=-9999)paste("(",response.units,")")),
        ylab="Standardized residuals",main=paste("Standardized Residuals vs. Fitted values of",response.name))
abline(0,0,lwd=2.5)

# saves graph
savePlot(rVars$info$Residual,type="emf")

# plot residuals by dose group
resbydosegrp <- xyplot(dataset$resid.sd~dataset$Fitted | as.factor(dose), 
       xlab=paste("Fitted values",if(response.units!=-9999)paste("(",response.units,")")),ylab="Standardized Residuals",
       main=paste("Standardized Residuals vs. Fitted Values of",response.name,"by",dose.name),
       panel.groups=paste(levels(as.factor(dose)), dose.units),
       layout=c(length(levels(as.factor(dose))),1,1),as.table=T,par.settings = list(par.main.text = list(cex = 0.9)))

# saves graph
plot(resbydosegrp)
savePlot(rVars$info$Residual_Grouped_By_Dose,type="emf")

#################################################
# Part 3: Estimate of Benchmark Dose(BMD)       #
#################################################


######## Preparing for Bootstrap ##########################################
# Create a data frame that include ID, dose, time, and residuals for      # 
# bootstrap. Also create another list that contains estimated parameters, #
# random effects, variance matrix of random effects, and variance matrix  #
# of residuals, and overall standard deviation.                           #
###########################################################################

### Creating the data frame for residuals 
### the program requires data organized in the order of ID, dose, time

ResData<-dataset[!is.na(dataset$response), 1:3]
ResData$Res<-ResData$BtrapRes<-as.data.frame(datasetm$resid)$ID

model.list<-NULL
model.list$ran<-as.matrix(datasetm$coef$random$ID) # extracting random effects 
ranVar<-as.matrix(datasetm$modelStruct$reStruct$ID)
model.list$Dmatr<-diag(as.numeric(ranVar),length(unique(dataset$ID)))      

model.list$Lambda<-diag(rep(1,length(ID)))
model.list$sigma<-datasetm$sigma
model.list$para<-datasetm$coef$fixed

# Scale dose and time to be in [0,1] for BMD computation and retained in 
# model.list 
# (The scaled dose and time are assumed by the BMD program) 
# so we need to scale time and dose as well as model parameters accordingly

# creates variables for later use
maxtime<-max(unique(ResData$time))
maxdose<-max(unique(ResData$dose))

# Scale dose and time to be in [0,1] for BMD computation and retained in 
# model.list 
# (The scaled dose and time are assumed by the BMD program) 
# so we need to scale time and dose as well as model parameters accordingly

N.para<-length(model.list$para)
model.list$para[N.para]<-datasetm$coef$fixed[N.para]*maxtime
model.list$para[c(N.para-2, N.para-1)]<-(datasetm$coef$fixed[c(N.para-2, N.para-1)])*maxtime*maxdose
model.list$dose<-ResData$dose/maxdose
model.list$time<-ResData$time/maxtime
model.list$ExpoTime<-ExpoTime/maxtime


# BMR/BMD/Bootstrap Options
#alpha.BMR= scan(options.file,skip=9,n=1,what=numeric(), quiet=T)
#	if(alpha.BMR==-9999) alpha.BMR=0.05                   #default alpha.BMR= 0.05
alpha.BMR=rVars$info$BMRL
if(rVars$info$DebugFile!="")
capture.output(cat("alpha.BMR: ",alpha.BMR,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#risk.type = scan(options.file,skip=10,n=1,what=character(), quiet=T)
#	if(risk.type==-9999) risk.type="extra"       #default risk.type="extra"
risk.type=rVars$info$RISK
if(rVars$info$RISK=="Extra") risk.type="extra"
if(rVars$info$RISK != "Extra")
	risk.type="additional"
if(rVars$info$DebugFile!="")
capture.output(cat("risk.type: ",risk.type,"\n"), file = rVars$info$DebugFile, append = TRUE)

#threshold = scan(options.file,skip=12,n=1,what=character(), quiet=T)
#	if(threshold==-9999) threshold=NULL
#	if(!is.null(threshold)) if(threshold=="F") threshold=NULL
#	if(!is.null(threshold)) if(threshold=="FALSE") threshold=NULL
#	if(!is.null(threshold)) threshold=as.numeric(threshold)

if(rVars$info$ADDE == "Background Rate")
{
	threshold=NULL
	threshold.upper=NULL
}

if(rVars$info$ADDE == "Cut Point")
{
	threshold = rVars$info$THOL
	threshold.upper=rVars$info$UPTL
}

if(rVars$info$DebugFile!="")
capture.output(cat("threshold: ",threshold,"\n"), file = rVars$info$DebugFile, append = TRUE)

#threshold.upper = scan(options.file,skip=13,n=1,what=character(), quiet=T)
	if(is.null(threshold)) threshold.upper=NULL
	if(!is.null(threshold.upper)) if(threshold.upper==-9999) threshold.upper=NULL
	if(!is.null(threshold.upper)) if(threshold.upper=="F") threshold.upper=NULL
	if(!is.null(threshold.upper)) if(threshold.upper=="FALSE") threshold.upper=NULL
	if(!is.null(threshold.upper)) threshold.upper=as.numeric(threshold.upper)
	
uppertail=rVars$info$UPTL
if(rVars$info$UPTL==-9999) uppertail=FALSE
if(rVars$info$DebugFile!="")
capture.output(cat("uppertail: ",uppertail,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#spontaneous.risk=scan(options.file,skip=11,n=1,what=numeric(), quiet=T)
#	if(spontaneous.risk==-9999) spontaneous.risk=0.05     #default...will be ignored if threshold!=F
spontaneous.risk = rVars$info$SPRL
if(rVars$info$DebugFile!="")
capture.output(cat("spontaneous.risk: ",spontaneous.risk,"\n"), file = rVars$info$DebugFile, append = TRUE)

BMDMethod=rVars$info$BMDM
if(rVars$info$DebugFile!="")
capture.output(cat("BMDMethod: ",BMDMethod,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#BMDMethod=scan(options.file,skip=14,n=1,what=character(), quiet=T)
#	if(BMDMethod==-9999) BMDMethod="Lowertail"            #default is "Lowertail"
	if(BMDMethod=="Bothside") if(is.null(threshold.upper)) threshold=NULL  
										#both thesholds must be given in twosided function if a threshold is used
	if(BMDMethod=="Bothside") if(!is.null(threshold)) threshold=c(threshold,threshold.upper)

uppertail=(BMDMethod=="Uppertail")  	                  #default is F (lowertail)
if(BMDMethod=="Uppertail") threshold=threshold.upper

#Climit.level=1-scan(options.file,skip=15,n=1,what=numeric(), quiet=T)
#	if(Climit.level==10000) Climit.level=0.05             #default is 0.05 (95% CI)
Climit.level=rVars$info$TCLVL
if(rVars$info$DebugFile!="")
capture.output(cat("Climit.level: ",Climit.level,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#CI.2sided=scan(options.file,skip=16,n=1,what=character(), quiet=T)
#	if(CI.2sided==-9999) CI.2sided=FALSE                  #default is for 1 sided (lower) interval (FALSE)
#	CI.2sided=as.logical(CI.2sided)
CI.2sided=as.logical(rVars$info$TSCI)	
if(rVars$info$DebugFile!="")
capture.output(cat("CI.2sided: ",CI.2sided,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#N.Bootstrap<-scan(options.file,skip=17,n=1, quiet=T)
#	if(N.Bootstrap==-9999) N.Bootstrap=2000               #number of bootstrap replications, default is 2000
N.Bootstrap=rVars$info$NBSR
if(rVars$info$DebugFile!="")
capture.output(cat("N.Bootstrap: ",N.Bootstrap,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#timepts<-scan(options.file,skip=18,n=1, quiet=T)
#	if(timepts==-9999) timepts=100                        #number of time points, default is 100
timepts=rVars$info$NTPT
if(rVars$info$DebugFile!="")
capture.output(cat("timepts: ",timepts,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
time.seq<-seq(1/timepts,1,1/timepts)

if(BMDMethod!="Bothside") {
	data.BMD.alphaBMR<-BMDNormOneSide(risk.type=risk.type,FUN="ToxicoDiffusion0",time=time.seq,
      	           param=model.list$para, spontaneous.risk=spontaneous.risk, risk.level=alpha.BMR, 
                       stdev=model.list$sigma, upper.tail=uppertail, threshold=threshold, 
                       controlData=dataset$response[dataset$time==0])}
if(BMDMethod=="Bothside") {
	data.BMD.alphaBMR<-BMDTwoSide(risk.type=risk.type,FUN="ToxicoDiffusion0",time=time.seq,
      	           param=model.list$para, spontaneous.risk=spontaneous.risk, risk.level=alpha.BMR, 
                       stdev=model.list$sigma, threshold=threshold, controlData=dataset$response[dataset$time==0],
                       interval=c(0,maxdose))}

#######################################
# Part 4: Bootstrap BMDL              # 
#######################################
#Call the functions: BtrapPARA, BtrapBMD, and BtrapBMDL  
#then, calculate the BMDs at 5% and 10% risk levels  
# and the BMDLs at 95% and 99% confidence level.    

#Call the function BtrapPARA to generate B(e.g. 1000) copies of bootstrap estimate of the model parameters 

data.BtrPARA<-BtrapPARA(DATA=ResData, model.list=model.list, B=N.Bootstrap, 
          dim=length(model.list$para), Zcol=1,Model=c("ToxicoDiffusion","ToxicoDiffusion1","ToxicoDiffusion2")[baselinePoly], BtrRandom="Yes", BtrMethod="CrossData")
# Call the function BtrapBMD to get "N.Bootstrap" (default=2000) bootstrapped BMD time-profile curves 
# of "timepts" (default=100) time points.

data.BtrBMD.alphaBMR<-BtrapBMD(BtrPARA=data.BtrPARA,B=N.Bootstrap,risk.type=risk.type,
      FUN="ToxicoDiffusion0",BMDMethod=BMDMethod,time=time.seq,
      spontaneous.risk=spontaneous.risk,risk.level=alpha.BMR,std=model.list$sigma, 
	threshold=threshold, controlData=dataset$response[dataset$time==0])

if(CI.2sided==F)
data.BtrBMDL.alphaBMRCL<-BtrapBMDL(BtrBMD=data.BtrBMD.alphaBMR,Climit.level=Climit.level)

if(CI.2sided==T)
data.BtrBMDL.alphaBMRCL<-BtrapBMDL(BtrBMD=data.BtrBMD.alphaBMR,Climit.level=Climit.level/2)
data.BtrBMDL.alphaBMRCU<-BtrapBMDL(BtrBMD=data.BtrBMD.alphaBMR,Climit.level=1-Climit.level/2)

# Dump all the results to a file named by BtrRESULTS (if options file says to do so).
# When you need the results next time, you may simply use data.restore S-Plus function to unload 

#dump.BtrapResults=scan(options.file,skip=21,n=1,what=character(), quiet=T)[1];
#	if(dump.BtrapResults==-9999) dump.BtrapResults=FALSE;
#	dump.BtrapResults<-as.logical(dump.BtrapResults)
dump.BtrapResults=rVars$info$SBSR
if(rVars$info$DebugFile!="")
capture.output(cat("dump.BtrapResults: ",dump.BtrapResults,"\n"), file = rVars$info$DebugFile, append = TRUE)

if(dump.BtrapResults)
dump(c("dataset","data.BMD.alphaBMR", "data.BtrPARA","data.BtrBMD.alphaBMR","data.BtrBMDL.alphaBMRCL")
,"BtrapResults")

# Plot of the Bootstrap Time-Profile Benchmark Doses 

# reads xaxis and yaxis options from options file
#xmin=scan(options.file,skip=19,n=4, quiet=T)[1]; if(xmin==-9999) xmin=min(time.orig)+ExpoTime
#xmax=scan(options.file,skip=19,n=4, quiet=T)[2]; if(xmax==-9999) xmax=max(time.orig)
#ymin=scan(options.file,skip=19,n=4, quiet=T)[3]; if(ymin==-9999) ymin=0
#ymax=scan(options.file,skip=19,n=4, quiet=T)[4]; if(ymax==-9999) ymax=max(data.BMD.alphaBMR*maxdose,median(dose, na.rm=TRUE),na.rm=T)
xmin=rVars$info$XAMI; if(xmin==-9999) xmin=min(time.orig)+ExpoTime
xmax=rVars$info$XAMX; if(xmax==-9999) xmax=max(time.orig)
ymin=rVars$info$YAMI; if(ymin==-9999) ymin=0
ymax=rVars$info$YAMX; if(ymax==-9999) ymax=max(data.BMD.alphaBMR*maxdose,median(dose, na.rm=TRUE),na.rm=T)
if(rVars$info$DebugFile!="")
capture.output(cat("\nX-Y Axis:",xmin,xmax,ymin,ymax,"\n"), file = rVars$info$DebugFile, append = TRUE)


# Creates bootstrap plot

if(logscale){
plot(c(xmin+logfix,xmax+logfix),c(ymin,ymax),log="x",xaxt="n",type="n",
xlab=paste(time.name,if(time.units!=-9999)paste("(",time.units,")")),
ylab=paste("Benchmark Doses",if(dose.units!=-9999)paste("(",dose.units,")")),
sub=paste("Dashed lines(s) is", (1-Climit.level)*100,"%","Confidence Band(s)"),
main=paste("BMD Time-Profile Based on",response.name,"\n (", risk.type, "Risk at",alpha.BMR*100,"% BMR Level)"))
axis(side=1, at=(time.level+min.time.orig+logfix), label=(time.level+min.time.orig))
for(i in 1:N.Bootstrap)
	{lines(maxtime*time.seq+min.time.orig+logfix,data.BtrBMD.alphaBMR[,i]*maxdose,col=8)}
lines(maxtime*time.seq+min.time.orig+logfix,data.BMD.alphaBMR*maxdose,lty=1,lwd=2)
lines(maxtime*time.seq+min.time.orig+logfix,data.BtrBMDL.alphaBMRCL*maxdose,lty=2,lwd=2)
if(CI.2sided==T) lines(maxtime*time.seq+min.time.orig+logfix,data.BtrBMDL.alphaBMRCU*maxdose,lty=2,lwd=2)
}

if(!logscale){
plot(c(xmin,xmax),c(ymin,ymax),
xlab=paste(time.name,if(time.units!=-9999)paste("(",time.units,")")),type="n",
ylab=paste("Benchmark Doses",if(dose.units!=-9999)paste("(",dose.units,")")),
sub=paste("Dashed lines(s) is", (1-Climit.level)*100,"%","Confidence Band(s)"),
main=paste("BMD Time-Profile for",response.name,"\n (", risk.type, "Risk at",alpha.BMR*100,"% BMR Level)"))
axis(side=1, at=(time.level+min.time.orig))
for(i in 1:N.Bootstrap)
	{lines(maxtime*time.seq+min.time.orig,data.BtrBMD.alphaBMR[,i]*maxdose,col=8)}
lines(maxtime*time.seq+min.time.orig,data.BMD.alphaBMR*maxdose,lty=1,lwd=2)
lines(maxtime*time.seq+min.time.orig,data.BtrBMDL.alphaBMRCL*maxdose,lty=2,lwd=2)
if(CI.2sided==T) lines(maxtime*time.seq+min.time.orig,data.BtrBMDL.alphaBMRCU*maxdose,lty=2,lwd=2)
}

# saves bootstrap plot
savePlot(rVars$info$Bootstrap,type="emf")


# saves dataset to file specified in options file (default is to not save)
# saves dataset (with residuals) as a tab deliminated ASCII text file.
# add original scaled time to output file
dataset2<-dataset
dataset2$time.orig<-dataset$time+min.time.orig

#The following was disabled for now, GLN 09/17/2008
#if(scan(options.file,skip=3,n=1,what=character(),quiet=T)[1]!=-9999)
#	write.table(dataset2, file=scan(options.file,skip=3,n=1,what=character(),quiet=T)[1],sep="\t")

#Find minimum BMD (and time, and confidence limits of that BMD)
minimumBMD<-min(data.BMD.alphaBMR*maxdose,na.rm=T)
if(rVars$info$DebugFile!="")
capture.output(cat("\ndata.BMD.alphaBMR:",data.BMD.alphaBMR,"\nmaxdose:",maxdose), file = rVars$info$DebugFile, append = TRUE)

minimumBMD.index<-seq(along = data.BMD.alphaBMR*maxdose)[data.BMD.alphaBMR*maxdose == minimumBMD]
minimumBMD.index<-minimumBMD.index[!is.na(minimumBMD.index)]
minimumBMD.time<-maxtime*time.seq[minimumBMD.index]+min.time.orig+ExpoTime
minimumBMD.ConfidenceLimit<-(data.BtrBMDL.alphaBMRCL*maxdose)[minimumBMD.index]
if(CI.2sided) minimumBMD.ConfidenceLimit.upper<-(data.BtrBMDL.alphaBMRCU*maxdose)[minimumBMD.index]

#creation of variables for output file

#adds names of variables to the initial value matrix
possible.initial.values<-data.frame(possible.init)
names(possible.initial.values)<-c("A0","B0","C0","K0")
if(length(possible.initial.values[1,])==5) names(possible.initial.values)<-c("A0","A1","B0","C0","K0")
if(length(possible.initial.values[1,])==6) names(possible.initial.values)<-c("A0","A1","A2","B0","C0","K0")
if(rVars$info$DebugFile!="")
capture.output(possible.initial.values, file = rVars$info$DebugFile, append = TRUE)

#model form
functional.form<-sprintf("%s%s",c("A","A0+A1*time","A0+A1*time+A2*time^2")[baselinePoly],"+B*dose*time*exp(-K*time)/(1+C*dose*time*exp(-K*time))")

dose.level.units<-dose.level<-levels(as.factor(dose))
minimumBMD.units<-round(minimumBMD,6)
minimumBMD.ConfidenceLimit.units<-round(minimumBMD.ConfidenceLimit,6)
if(CI.2sided) minimumBMD.ConfidenceLimit.upper.units<-round(minimumBMD.ConfidenceLimit.upper,6)
if(dose.units!=-9999) {
	dose.level.units<-paste(dose.level,dose.units)
	minimumBMD.units<-paste(minimumBMD.units,dose.units)
	minimumBMD.ConfidenceLimit.units<-paste(minimumBMD.ConfidenceLimit.units,dose.units)
	if(CI.2sided) minimumBMD.ConfidenceLimit.upper.units<-paste(minimumBMD.ConfidenceLimit.upper.units,dose.units)
	}
time.level.units<-time.level.orig
ExpoTime.units<-ExpoTime
minimumBMD.time.units<-minimumBMD.time
if(time.units!=-9999) {
	time.level.units<-paste(time.level,time.units)
	ExpoTime.units<-paste(ExpoTime,time.units)
	minimumBMD.time.units<-paste(minimumBMD.time.units,time.units)
	}
if(!is.null(threshold)) {
	threshold.units<-threshold
	if(response.units!=-9999) threshold.units<-paste(threshold.units,response.units)
	}
sample.size<-datasetm$dims$N #get sample size
randomeffects.stdDev<-data.frame((t(VarCorr(datasetm)[,2]))) #get random effects information
	names(randomeffects.stdDev)<-c("A", "Residual") #put in attractive column names
	row.names(randomeffects.stdDev)<-c("StdDev:") #put in an attractive row name
	
#Variables to be read from option file.  Default (if -9999) is to not be included in output file.
#study.name<-scan(options.file,skip=25,n=1,what=character(), quiet=T)
#	if(study.name==-9999) study.name<-paste(response.name,"Data")
study.name=rVars$info$USERNOTES
if(rVars$info$DebugFile!="")
capture.output(cat("study.name: ",study.name,"\n"), file = rVars$info$DebugFile, append = TRUE)
	
#chemical<-scan(options.file,skip=26,n=1,what=character(), quiet=T)
chemical=rVars$info$CHNA
if(rVars$info$DebugFile!="")
capture.output(cat("chemical: ",chemical,"\n"), file = rVars$info$DebugFile, append = TRUE)

#exposure<-scan(options.file,skip=27,n=1,what=character(), quiet=T)
exposure=rVars$info$EXTY
if(rVars$info$DebugFile!="")
capture.output(cat("exposure: ",exposure,"\n"), file = rVars$info$DebugFile, append = TRUE)

#species<-scan(options.file,skip=28,n=1,what=character(), quiet=T)
species=rVars$info$SPEC
if(rVars$info$DebugFile!="")
capture.output(cat("species: ",species,"\n"), file = rVars$info$DebugFile, append = TRUE)

#sex<-scan(options.file,skip=29,n=1,what=character(), quiet=T)
sex=rVars$info$GEND
if(rVars$info$DebugFile!="")
capture.output(cat("sex: ",sex,"\n"), file = rVars$info$DebugFile, append = TRUE)

outFile = rVars$info$OutFile

# saves model results to summary output file (outFile)
#write(paste("Analysis of",study.name),file="outFile")
capture.output(cat("Analysis of ",study.name,"\n"), file = outFile, append = FALSE)

#write("\nSTUDY DESCRIPTION",file=outFile,append=TRUE)
cat("STUDY DESCRIPTION","\n", file = outFile, append = TRUE)

if(chemical!="na")cat("Chemical:\t",chemical,"\n",file=outFile,append=TRUE)

cat("Dose Levels:\t",dose.level.units,"\n",file=outFile,append=TRUE)
cat("Test Times:\t",time.level.units,"\n",file=outFile,append=TRUE)

if(exposure!="na") cat("Exposure:\t",exposure,"\n",file=outFile,append=TRUE)

cat("Exposure Time:\t",ExpoTime.units,"\n",file=outFile,append=TRUE)

cat("Sample Size:\t",sample.size,"\n",file=outFile,append=TRUE)

if(species!="na") cat("Species:\t",species,"\n",file=outFile,append=TRUE)

if(sex!="na") cat("Sex:\t\t",sex,"\n",file=outFile,append=TRUE)

cat("\nDOSE-RESPONSE MODELING\n",file=outFile,append=TRUE)
cat(functional.form,"\n\n",file=outFile,append=TRUE)

#if(rVars$info$DebugFile!="")
#capture.output(datasetm, file = rVars$info$DebugFile, append = TRUE)

#sink(file=outFile,append=TRUE)
	capture.output(data.frame(AIC=summary(datasetm)$AIC,BIC=summary(datasetm)$BIC,logLik=datasetm$logLik,row.names=c(" ")), file = outFile, append = TRUE)

	cat("\nRandom effects:",file=outFile,append=TRUE)

	cat("\nFormula: A ~ 1 | ID\n",file=outFile,append=TRUE)

	capture.output(randomeffects.stdDev, file = outFile, append=TRUE)

	cat("\nFixed effects:\n",file=outFile,append=TRUE)
	capture.output(summary(datasetm)$tTable, file=outFile,append=TRUE)

	cat("\nCorrelation:\n",file=outFile,append=TRUE)

	capture.output(summary(datasetm)$corFixed, file=outFile,append=TRUE)

	cat("\nStandardized Within-Group Residuals:\n",file=outFile,append=TRUE)
	capture.output(summary(datasetm)$residuals, file=outFile,append=TRUE)
#sink()

cat("\nInitial Values:\t",datasetm$initialValues,"\n\n",file=outFile,append=TRUE)

cat("Possible Initial Values\n",file=outFile,append=TRUE)

#sink(file=outFile,append=TRUE)
	capture.output(possible.initial.values, file=outFile, append=TRUE)
#sink()

cat("\nBENCHMARK DOSE ESTIMATION",file=outFile,append=TRUE)
cat("\nRisk Type:\t\t\t",risk.type,file=outFile,append=TRUE)

if(is.null(threshold)) cat("\nSpontaneous Risk Level:\t\t",spontaneous.risk*100,"%",file=outFile,append=TRUE)
#cat("\nSpontaneous Risk Level:\t\t",spontaneous.risk*100,"%",file=outFile,append=TRUE)

if(!is.null(threshold)) cat("\nThreshold:\t\t\t",threshold.units,file=outFile,append=TRUE)

cat("\nArea of Adverse Effects:\t",BMDMethod,file=outFile,append=TRUE)
cat("\nBMR Level:\t\t\t",alpha.BMR*100,"%\n",file=outFile,append=TRUE)

cat("\nBOOTSTRAP ESTIMATION OF BMDL\n",file=outFile,append=TRUE)
cat("\nBootstrap Replications:\t\t",N.Bootstrap,file=outFile,append=TRUE)
cat("\nMinimum BMD:\t\t\t",minimumBMD.units,file=outFile,append=TRUE)
cat("\nAt Test Time:\t\t\t",minimumBMD.time.units,file=outFile,append=TRUE)
cat("\nConf. Level:\t\t\t",100*(1-Climit.level),"%",file=outFile,append=TRUE)
cat("\nBMDL:\t\t\t\t",minimumBMD.ConfidenceLimit.units,file=outFile,append=TRUE)
if(CI.2sided) cat("\nBMDU:\t\t\t\t",minimumBMD.ConfidenceLimit.upper.units,file=outFile,append=TRUE)


if(as.character(datasetm)[1]==-9999) 
cat("Toxico Diffusion model could not be fit due to inability to find a valid starting value using initial value algorithm. " ,
	"Please enter new starting value manually in the options file.", 
	"If an initial value was entered manually, please try another initial value or consider fitting a different model.",
	file=outFile, append=TRUE)