#********************************************************************************
# Input from "MasterEstimator.Input Parameters.Oct13.2011.R"
$Macros
#**************************************************************************************
require(survival)


#Change Data Formats ####################################################################################

####Cas: Remove Erroneous NA rows, Create DWP 100%
Cas[,turbine.col]	 <- as.character(Cas[,turbine.col])
if (is.na(propsrch)) {
	Cas <- Cas[!is.na(Cas[,turbine.col]),]
	propsrch            <- "Hundred"
	Cas[,propsrch]      <- 1 
	} 
if (!is.na(propsrch)) { 
	Cas[,propsrch]	 <- as.numeric(as.character(Cas[,propsrch])) 
	if(sum(is.na(Cas[,turbine.col])*is.na(Cas[,propsrch]))>=1)  {Cas <-  Cas[!is.na(Cas[,turbine.col]),]}
	}

####Cas
Cas.names	    <- c(names(Cas)) #HERE
Cas[,date.col]      <- as.Date(as.character(Cas[,date.col]), "%m/%d/%Y")
Cas[,prev.date.col] <- as.Date(as.character(Cas[,prev.date.col]), "%m/%d/%Y")
Cas[,"Interval"]    <- as.numeric(difftime(Cas[,date.col], Cas[,prev.date.col], units = c("days")))
Cas[,date.col]      <- as.factor(Cas[,date.col])        #to preserve when converted to matrix
Cas[,prev.date.col] <- as.factor(Cas[,prev.date.col])
INT.calc	    <- min(as.numeric(names(table(Cas[,"Interval"]))[table(Cas[,"Interval"])==max(table(Cas[,"Interval"]))]))
Turbines            <- unique(Cas[,turbine.col])
how.many.reqs       <- length(summary.list)


####SE: From Factor to Character, Remove Erroneous NA rows, SE 0&1 values
SE[,response.col.SE] <- as.character(SE[,response.col.SE])
SE <- SE[!is.na(SE[,response.col.SE]),]
if (length(unique(SE[,response.col.SE]))<=2){
SE[,response.col.SE] <- as.numeric(SE[,response.col.SE] == max(SE[,response.col.SE])) }

####CP: From Factor/Character to Numeric, Remove Erroneous NA rows
CP[, event.col]      <- as.numeric(as.character(CP[,event.col]))  
CP[, Left.col]       <- as.numeric(as.character(CP[,Left.col]))   
CP[, Right.col]      <- as.numeric(as.character(CP[,Right.col]))  
if(sum(is.na(CP[,event.col])*is.na(CP[,Left.col]))>=1){CP <- CP[!is.na(CP[,event.col]),]}

####SE and CP: Covariate Lists 
coef.list.SE <- SE.model.covariates
	if(!is.na(SE.model.covariates)) {coef.list.SE <- unique(unlist(strsplit(strsplit(SE.model.covariates," \\+ ")[[1]]," \\: ")))}
coef.list.CP <- CP.model.covariates
	if(!is.na(CP.model.covariates)) {coef.list.CP <- unique(unlist(strsplit(strsplit(CP.model.covariates," \\+ ")[[1]]," \\: ")))}


if(!identical(coef.list.SE,NA)){
		how.many.SEcovs<-length(coef.list.SE)
		if.SEcovs<-TRUE}
if(!identical(coef.list.CP,NA)){
		how.many.CPcovs<-length(coef.list.CP)
		if.CPcovs<-TRUE}

if(identical(coef.list.SE,NA)){
		if.SEcovs<-FALSE
		coef.list.SE<-"DUMMY.SE"
		SE.model.covariates<-"DUMMY.SE"
		how.many.SEcovs<-length(coef.list.SE)
		SE$DUMMY.SE<-1
		Cas$DUMMY.SE<-1
		}
if(identical(coef.list.CP,NA)){
		if.CPcovs<-FALSE
		coef.list.CP<-"DUMMY.CP"
		CP.model.covariates<-"DUMMY.CP"
		how.many.CPcovs<-length(coef.list.CP)
		CP$DUMMY.CP<-1
		Cas$DUMMY.CP<-1
		}

#Create Error Checks ###############################################################################

Data.Error.Check.Col.Names<-function(Cas, CP, SE, coef.list.SE, coef.list.CP){
	

	# Check to make sure that the interval between searches >0
	check.int <- sum(Cas[,"Interval"]<=0)
	if(check.int >0)
		return(list(Err = 1, 
		message ="Fatality data date of previous search >= search date.  Check dates."))

	# Check to make sure that response variable has only 2 levels, and can be recoded to 0,1
	check.levels.SE.resp <- length(unique(SE[,response.col.SE]))
	if(check.levels.SE.resp > 2)
		return(list(Err = 1, 
		message ="SE response variable has >2 values. Check spelling and case sensitivity."))

	# Check to make sure that explanatory variables selected for use in adjusting
	# observed casualties exist in the searcher efficiency and carcass persistence data
	check.covariates.match.SE<-sum(coef.list.SE%in%names(Cas)==0)
	if(check.covariates.match.SE>0)
		return(list(Err = 1, 
		message ="Fatality data does not contain explanatory variables indicated by searcher efficiency data. Check for spelling errors and case sensitivity."))

	check.covariates.match.CP<-sum(coef.list.CP%in%names(Cas)==0)
	if(check.covariates.match.CP>0)
		return(list(Err = 1, 
		message ="Fatality data does not contain explanatory variables indicated by carcass persistence data. Check for spelling errors and case sensitivity."))
    		
	#check that there are no missing values in fatality data for chosen SEcovariate 
	check.miss.SE <- 0
	if(!identical(coef.list.SE,NA)){
	check.miss.SE1 <-  sum(is.na(Cas[,coef.list.SE]))
	check.miss.SE2 <-  nrow(Cas[which(Cas[,coef.list.SE]==""),])
	check.miss.SE <- check.miss.SE1+check.miss.SE2}
	if(check.miss.SE>0)
		return(list(Err=1,
		message="Covariate variable has missing value in fatality data"))
	
	#check that there are no missing values in fatality data for chosen CPcovariate 
	check.miss.CP <- 0
	if(!identical(coef.list.CP,NA)){
	check.miss.CP1 <-  sum(is.na(Cas[,coef.list.CP]))
	check.miss.CP2 <-  nrow(Cas[which(Cas[,coef.list.CP]==""),])
	check.miss.CP <- check.miss.CP1+check.miss.CP2}
	if(check.miss.CP>0)
		return(list(Err=1,
		message="Covariate variable has missing value in fatality data"))

	#### run these two lines after all the separate error checks
	Err <- list(Err = 0)
  	Err
} # ends Data.Error.Check.Col.Names


Data.Error.Check.Suff.Data<-function(Cas, CP, SE, coef.list.SE, coef.list.CP){

	## check to make sure that the number of turbines in the study is less than the number of turbines at the site 
	## check number of turbines in the study if greater than = number found in Fatality data
	## check for missing Cas turbine IDs
	## check for missing CP Censoring Type values
	## check for missing Cas DWP values
	## check that DWP is between 0 and 1
	## check that there are not missing values for summary groups
	## check CP left is greater than 0
	## check CP left is less than or equal to CP right
	## check CP left is less than CP right when event code=3 (interval censoring)
	## check CP event codes are 0,1,2,3
	## check to make sure that CP and SE data have at least 10 observations for all combos of
	## covariates, for now this includes all interaction wise combos
	## check to make sure that SE do not have only 0 (not found) data for some combo of covariates
	## Check that if CP has 100% right censored data the values are reasonable relative the interval length


    #Check Nturbines in study is less than Nturbines at site
	if(Num.Turbines.in.Study>Num.Total.Turbines.Site)
		return(list(Err = 1, 
		message =" Total Number of Turbines in the Study should be less than or equal to the Total Number of Turbines at the Site "))
    
    #Check Nturbines in study is greater than how.many.turbines (unique in Cas)
	if(Num.Turbines.in.Study<length(Turbines))
		return(list(Err = 1, 
		message =" Total Number of Turbines in the Study should be greater than or equal to the number of values found in the Fatality data"))    

    #check for missing values in fatality Turbine ID 
	check.miss.ID <- sum(is.na(Cas[,turbine.col]))
	if(check.miss.ID>0)
		return(list(Err = 1,
		message="Missing turbine ID in fatality data"))

    #check for missing values in CP event 
	check.miss.event <- sum(is.na(CP[,event.col]))
	if(check.miss.event>0)
		return(list(Err = 1,
		message="Missing CP Censoring Type value"))

   #check DWP Missing
	if (!identical(propsrch, "Hundred")) {
	check.DWP <- sum(is.na(Cas[,propsrch]))
	if(check.DWP > 0)
		return(list(Err = 1,
		message="Missing DWP value in fatality data"))}

   #check that DWP are between 0 and 1
	if (!identical(propsrch, "Hundred")) {
	if(max(Cas[,propsrch]) > 1)
		return(list(Err = 1,
		message="Fatality DWP values should be between 0 and 1."))
	if(min(Cas[,propsrch]) <= 0)
   		return(list(Err = 1,
		message="Fatality DWP values should be between 0 and 1."))}

    #check that there are no missing values in summary groups
	check.miss.cov <- 0
	if(!identical(summary.list,NA)){
	check.miss.cov1 <-  sum(is.na(Cas[,summary.list]))
	check.miss.cov2 <-  nrow(Cas[which(Cas[,summary.list]==""),])
	check.miss.cov <- check.miss.cov1+check.miss.cov2}
	if(check.miss.cov>0)
		return(list(Err=1,
		message="Fatality summary variable has missing values"))

    #check that CP left values are greater than 0
	check.Left <- sum(is.na(CP[,Left.col]))
	if(check.Left > 0)
		return(list(Err = 1,
		message="Missing CP left column values"))
	if(min(CP[,Left.col])<=0)
		return(list(Err=1,
		message="CP left censored times must be greater than 0"))


    #check that CP right values are always greater than CP left values
	check.Right <- sum(is.na(CP[,Right.col]))
	if(check.Right > 0)
		return(list(Err = 1,
		message="Missing CP right column values"))	
	if(sum(CP[,Left.col] > CP[,Right.col]) >0)
		return(list(Err=1,
		message="CP left censored times must be less than or equal to right censored times"))

    #check that CP right values are strictly greater than CP left values for Event 3	
	if(sum(CP[which(CP[,event.col]==3),Left.col] == CP[which(CP[,event.col]==3),Right.col]) >0)
		return(list(Err=1,
		message="CP left censored times must be less than right censored times for CP trials in which Event=3"))
  
    #check CP event codes are between 0 and 3
	if(max(CP[,event.col]) > 3)
		return(list(Err = 1,
		message="CP Censoring Type values should be 0, 1, 2 or 3."))
	if(min(CP[,event.col]) < 0)
		return(list(Err = 1,
		message="CP Censoring Type values should be 0, 1, 2 or 3."))

    # Check to see if at least 10 of all SE combos 
	SE.combos1<-table(SE[,coef.list.SE]) #will note difference in spacing (should already be gone)
	SE.combos<-as.data.frame(SE.combos1, row.names = NULL,responseName = "Count", stringsAsFactors = TRUE)
	SE.Less.10<-sum(SE.combos$Count<10)
	if(SE.Less.10>0)
		return(list(Err = 1, 
		message ="Not enough searcher efficiency trials at all combinations of covariates (at least 10 of each needed)"))

    # Check to see if at least 10 of all CP combos
	CP.combos1<-table(CP[,coef.list.CP])
	CP.combos<-as.data.frame(CP.combos1, row.names = NULL,responseName = "Count", stringsAsFactors = TRUE)
	CP.Less.10<-sum(CP.combos$Count<10)
	if(CP.Less.10>0)
		return(list(Err = 1, 
		message ="Not enough carcass persistence trials at all combinations of covariates (at least 10 of each needed)"))
    
   # Revised code needs needs to alert user if Fat data contains extra covariate level, otherwise it will just exclude it from total 
   # Check to see if all SE combos needed for Fat are present 
	Cas.SE.combos1 <- table(Cas[,coef.list.SE])
	Cas.SE.combos<-as.data.frame(Cas.SE.combos1, row.names = NULL,responseName = "Count", stringsAsFactors = TRUE)
	if(how.many.SEcovs>1){Cas.SE.Miss<-sum(is.na(merge(Cas.SE.combos,SE.combos, by=coef.list.SE, all.x=T)$Count.y))}
	if(how.many.SEcovs==1){Cas.SE.Miss<-sum(is.na(merge(Cas.SE.combos,SE.combos, by="Var1", all.x=T)$Count.y))}
	if(Cas.SE.Miss>0)
		return(list(Err = 1, 
		message ="Missing searcher efficiency trials for all combinations of SE covariates in Fatality data"))

   # Check to see if all CP combos needed for Fat are present 
	Cas.CP.combos1 <- table(Cas[,coef.list.CP])
	Cas.CP.combos<-as.data.frame(Cas.CP.combos1, row.names = NULL,responseName = "Count", stringsAsFactors = TRUE)
	if(how.many.CPcovs>1){Cas.CP.Miss<-sum(is.na(merge(Cas.CP.combos,CP.combos, by=coef.list.CP, all.x=T)$Count.y))}
	if(how.many.CPcovs==1){Cas.CP.Miss<-sum(is.na(merge(Cas.CP.combos,CP.combos, by="Var1", all.x=T)$Count.y))}
	if(Cas.CP.Miss>0)
		return(list(Err = 1, 
		message ="Missing carcass persistence trials for all combinations of CP covariates in Fatality data"))


    # Check to see if all observations in any SE combos are 0
	SE.zeroes1<-table(SE[SE[,response.col.SE]==1,coef.list.SE])
	SE.zeroes <-as.data.frame(SE.zeroes1, row.names = NULL,responseName = "Count", stringsAsFactors = TRUE)
	if(how.many.SEcovs>1){SE.All.0 <- sum(merge(Cas.SE.combos[which(Cas.SE.combos$Count != 0),], SE.zeroes, by=coef.list.SE, all.x=T)$Count.y == 0)}
	if(how.many.SEcovs==1){SE.All.0<-sum(merge(Cas.SE.combos, SE.zeroes, by="Var1", all.x=T)$Count.y == 0)}		
	if(SE.All.0>0 | is.na(SE.All.0) ==T)
		return(list(Err = 1, 
		message ="Some combinations of covariates have all SE=0. Cannot estimate fatality"))

    # Check that if CP has 100% right censored data the values are reasonable relative the interval length
    #1 there must be at % CP trial longer than the max search interval that is > % interval in Fat data
    #2 or >10 CP trials longer than each search interval
    #3 if covariates, each subset (with right censored data only) must pass the above conditions 
	CP.zeroes <- 0
	if(sum(CP[,event.col])==0 & if.CPcovs==F){
		tab <- as.data.frame(table(Cas[,"Interval"]))
		tab$FatPerc <- table(Cas[,"Interval"])/length(Cas[,"Interval"])
		for (i in 1:nrow(tab)){
		tab[i,"CPCumCount"] <- sum(CP[,Right.col] > as.numeric(as.character(tab[i,1])))
		}
		tab$CPCumPerc <- tab$CPCumCount/length(CP[,Right.col])
		CP.zeroes <- sum(tab$CPCumCount<10 & tab$FatPerc>tab$CPCumPerc)}
	if(sum(CP[,event.col])==0 & if.CPcovs==T){  
		#Check all subsets, data is completely right censored
		tab <- as.data.frame(table(Cas[,"Interval"],Cas[,"CP.Cov.Ind"]))
		tab$FatPerc <- tab$Freq/rep(aggregate(tab$Freq,by=list(tab$Var2),FUN="sum")$x,each=(nrow(tab)/length(unique(tab$Var2))))
		for (i in 1:nrow(tab)){
		tab[i,"CPCumCount"] <- sum(CP[CP$CP.Cov.Ind == tab$Var2[i],Right.col] > as.numeric(as.character(tab$Var1[i])))
		}
		tab <- merge(tab,aggregate(as.data.frame(table(CP[,Right.col],CP[,"CP.Cov.Ind"]))$Freq,by=list(as.data.frame(table(CP[,Right.col],CP[,"CP.Cov.Ind"]))$Var2),FUN="sum"),by.x="Var2",by.y="Group.1")
		tab$CPCumPerc <- tab$CPCumCount/tab$x
		CP.zeroes <- sum(tab$CPCumCount<10 & tab$FatPerc>tab$CPCumPerc)}
	if(sum(CP[,event.col])!=0 & if.CPcovs==T & sum(colSums(table(CP[,event.col],CP[,"CP.Cov.Ind"],exclude = c(0)))==0)>0){  
		#Check certain subsets that have only right censored
		tab <- as.data.frame(table(Cas[,"Interval"],Cas[,"CP.Cov.Ind"]))
		exclude <- as.data.frame(colSums(table(CP[,event.col],CP[,"CP.Cov.Ind"],exclude = c(0)))); names(exclude) <- c("ind");exclude$names <- names(colSums(table(CP[,event.col],CP[,"CP.Cov.Ind"],exclude = c(0))))
		tab <- merge(tab,exclude, by.y="names", by.x="Var2")
		tab$FatPerc <- tab$Freq/rep(aggregate(tab$Freq,by=list(tab$Var2),FUN="sum")$x,each=(nrow(tab)/length(unique(tab$Var2))))
		for (i in 1:nrow(tab)){
		tab[i,"CPCumCount"] <- sum(CP[CP$CP.Cov.Ind == tab$Var2[i],Right.col] > as.numeric(as.character(tab$Var1[i])))
		}
		tab <- merge(tab,aggregate(as.data.frame(table(CP[,Right.col],CP[,"CP.Cov.Ind"]))$Freq,by=list(as.data.frame(table(CP[,Right.col],CP[,"CP.Cov.Ind"]))$Var2),FUN="sum"),by.x="Var2",by.y="Group.1")
		tab$CPCumPerc <- tab$CPCumCount/tab$x
		CP.zeroes <- sum(tab[tab$ind==0,]$CPCumCount<10 & tab[tab$ind==0,]$FatPerc>tab[tab$ind==0,]$CPCumPerc)}
	if(CP.zeroes>0)
		return(list(Err = 1, 
		message ="The CP data is completely right censored or completely right censored for a covariate combination, there is insufficient data to estimate CP for each search interval in the fatality data. Cannot estimate fatality"))		


	#### run these two lines after all the separate error checks
	Err <- list(Err = 0)
  	Err
} # ends Data.Error.Check.Suff.Data

### First Check
#If RunIT is converted to a character string, error message will be displayed
	RunIT <- 0 
	Err<-Data.Error.Check.Col.Names(Cas, CP, SE, coef.list.SE, coef.list.CP)
	if(Err$Err == 1) {RunIT <- print(Err$message)}

if(mode(RunIT)!="character") {

################################################
#SE CP Covariate Data Checks (For this to run col.now must be in Cas and SE, First Check Passed)

#Formating of covariate columns, remove spacing differences
for(i in 1:how.many.SEcovs){
	col.now<-coef.list.SE[i]
	Cas[,col.now]<-as.character(Cas[,col.now])
for (j in 1:dim(Cas)[1]){
	Cas[j,col.now]<-paste(strsplit(as.matrix(Cas[j,col.now]), " "),collapse=".")}
	SE[,col.now]<-as.character(SE[,col.now])
for (j in 1:dim(SE)[1]){
	SE[j,col.now]<- paste(strsplit(as.matrix(SE[j,col.now]), " "),collapse=".")}
} # ends loop

for(i in 1:how.many.CPcovs){
	col.now<-coef.list.CP[i]
	Cas[,col.now]<-as.character(Cas[,col.now])
for (j in 1:dim(Cas)[1]){
	Cas[j,col.now]<-paste(strsplit(as.matrix(Cas[j,col.now]), " "),collapse=".")}
	CP[,col.now]<-as.character(CP[,col.now])
for (j in 1:dim(CP)[1]){
	CP[j,col.now]<- paste(strsplit(as.matrix(CP[j,col.now]), " "),collapse=".")}
} # ends loop

##Make sure coef.list.CP is ordered by large to small for CP plots
plotdims<-data.frame(rep(0,how.many.CPcovs),rep(0,how.many.CPcovs))
names(plotdims) <- c("dims","vars")
for (i in 1:how.many.CPcovs){
	plotdims[i,1] <- nrow(unique(CP[coef.list.CP][i]))  
	plotdims[i,2] <- names(unique(CP[coef.list.CP])[i])}
plotdims <- plotdims[order(-plotdims$dims),]
coef.list.CP <- unlist(plotdims$vars) 			#Reordered maxdims=firstlisted

if (how.many.CPcovs>=3){	                        #Can't plot more than 3 factors
CP$CP.Cov.IndPlotSub<-NA
for (i in 1:dim(CP)[1]){CP$CP.Cov.IndPlotSub[i]<-paste(strsplit(as.matrix(CP[i,coef.list.CP[1:3]]), " "),collapse=".")}
CP$CP.Cov.IndPlotSub<-as.factor(CP$CP.Cov.IndPlotSub)}


##Covariate Combinations in Cas Data
Cas$SE.Cov.Ind<-NA
Cas$CP.Cov.Ind<-NA
for (i in 1:dim(Cas)[1]){
Cas$SE.Cov.Ind[i]<-paste(strsplit(as.matrix(Cas[i,coef.list.SE]), " "),collapse=".")
Cas$CP.Cov.Ind[i]<-paste(strsplit(as.matrix(Cas[i,coef.list.CP]), " "),collapse=".")
	}
Cas$SE.Cov.Ind<-as.factor(Cas$SE.Cov.Ind)
Cas$CP.Cov.Ind<-as.factor(Cas$CP.Cov.Ind)

SE$SE.Cov.Ind<-NA
CP$CP.Cov.Ind<-NA	
for (i in 1:dim(SE)[1]){
	SE$SE.Cov.Ind[i]<-paste(strsplit(as.matrix(SE[i,coef.list.SE]), " "),collapse=".")}
	
for (i in 1:dim(CP)[1]){
	CP$CP.Cov.Ind[i]<-paste(strsplit(as.matrix(CP[i,coef.list.CP]), " "),collapse=".")}

#####List unique covariate combos 
SE.p.m.d.EST<-unique(SE[,c("SE.Cov.Ind",coef.list.SE)])
SE.p.m.d.EST<-SE.p.m.d.EST[order(SE.p.m.d.EST$SE.Cov.Ind) ,]
CP.p.m.d.EST<-unique(CP[,c("CP.Cov.Ind",coef.list.CP)])
CP.p.m.d.EST<-CP.p.m.d.EST[order(CP.p.m.d.EST$CP.Cov.Ind) ,]

### Second Check 
	Err<-Data.Error.Check.Suff.Data(Cas, CP, SE, coef.list.SE, coef.list.CP)
	if(Err$Err == 1) {RunIT <- print(Err$message)}


if(mode(RunIT)!="character") {

#Bootstrap Code###########################################################################

## This is a basic bootstrap resampling function.
Get.Boot<-function(data){
	datasize<-dim(data)[1]
	boot.row<-sample(seq(1:datasize),size=datasize, replace=TRUE)
	boot.sample <- rbind.data.frame(data[boot.row,])
	names(boot.sample)<-names(data)
	boot.sample
} # ends Get.Boot function

#### an alternate set of Get.Boot functions, this time, boot strap within all combos of covariates
Get.Boot.CovariateLevels<-function(data, covariates){
	### Need to make sure ahead of time that SE and CP data ONLY include values of covariates that Cas does,
	### this check needs to be done ONCE in master function so time isn't taken to re-do it every bootstrap iteration
	boot.by.combos<-by(data,data[,covariates], Get.Boot, simplify=TRUE)
	boot.sample<-do.call(rbind, lapply(boot.by.combos,data.frame))
	boot.sample
} # ends Get.Boot.CovariateLevels function


## small function for use by Get.Casualty.Boot, this
## gets the rows of data corresponding to the bootrapped turbines
which.row<-function(data, turbine.col,crnt.turbine){
	which(data[,turbine.col]==crnt.turbine)
}

## Bootstrap function for turbines of Fatality Data
Get.Casualty.Boot<-function(data, turbine.col){
	### Bootstrap the turbines, not the observed fatalities
	turbine.boot<-sort(sample(Turbines, size=length(Turbines), replace=TRUE))
	# which rows = these turbine numbers
	boot.row<-unlist(sapply(turbine.boot, which.row, turbine.col=turbine.col, data=data, simplify=TRUE))
	boot.sample<-rbind.data.frame(data[boot.row,])
	names(boot.sample)<-names(data)
	boot.sample
} # ends Get.Casualty.Boot

# this function adds boot.ID's, relative to iteration number, to each boostrapped estimate
boot.ID<-function(data, iteration.number){
	data$boot.ID<-rep(iteration.number, times=dim(data)[1])
	data
} # ends boot.ID


#Get Bootstrapped Data ###############################################################

# Get Bootstrap samples of Fatality Data (LIST1 is original)
	Fatality.obs2<-list(Cas)
	Fatality.obs<-list()
	Fatality.obs[1:num.boot.samples]<-Fatality.obs2 
	Fatality.Boot<-c(Fatality.obs2,lapply(Fatality.obs, Get.Casualty.Boot, turbine.col=turbine.col))
	names(Fatality.Boot)<-rep("Fatality.Data",(num.boot.samples+1))

# Get Bootstrap samples of Persistence Trials Data
	Persist.obs2<-list(CP)
	Persist.obs<-list()
	Persist.obs[1:num.boot.samples]<- Persist.obs2
	Persist.Boot<- c(Persist.obs2, lapply(Persist.obs[1:num.boot.samples], Get.Boot.CovariateLevels, covariates=c("CP.Cov.Ind")))
	names(Persist.Boot)<-rep("Persist.Trials",(num.boot.samples+1))

#Get Bootstrap samples of Searcher Efficiency Trials Data
	SearchEF.obs2<-list(SE)
	SearchEF.obs<-list()
	SearchEF.obs[1:num.boot.samples]<-SearchEF.obs2
	SearchEF.Boot<-c(SearchEF.obs2, lapply(SearchEF.obs, Get.Boot.CovariateLevels, covariates="SE.Cov.Ind"))  
	names(SearchEF.Boot)<-rep("SearchEF.Trials",(num.boot.samples+1))

# Create Boot ID
for (i in 1:(num.boot.samples+1)){Fatality.Boot[[i]]<-boot.ID(Fatality.Boot[[i]], i-1)}
for (i in 1:(num.boot.samples+1)){Persist.Boot[[i]]<-boot.ID(Persist.Boot[[i]], i-1)}
for (i in 1:(num.boot.samples+1)){SearchEF.Boot[[i]]<-boot.ID(SearchEF.Boot[[i]], i-1)}

#Calculate SE #############################################################################

## This function fits a logistic regression for the searcher effciency data,
## it outputs an object that is a fitted glm, for later use in predicting SE
## of the bootstrapped Fatality Data

SE.p.m.d.EST2 <- data.frame(SE.p.m.d.EST)

Search.Eff.logit<-function(SE.data, SE.model.covariates, response.col.SE, if.SEcovs, SE.p.m.d.EST){
	if(if.SEcovs==TRUE){fmla<-as.formula(paste(response.col.SE,SE.model.covariates, sep="~"))}
	if(if.SEcovs==FALSE){fmla<-as.formula(paste(response.col.SE,1, sep="~"))}
	Srch.EF<-glm(fmla, family=binomial(link = "logit"),data=SE.data)
	

	if(if.SEcovs==TRUE){
		
		SE.pred<-predict(Srch.EF, SE.p.m.d.EST, type="response") 	  
		SE.to.merge <- data.frame(SE.p.m.d.EST$SE.Cov.Ind, SE.pred)	  
		names(SE.to.merge)<-c("SE.Cov.Ind","SE.pred")

										}
	if(if.SEcovs==FALSE){
		SE.pred<-predict(Srch.EF, type="response")[1]
		SE.to.merge<-data.frame(SE.p.m.d.EST$SE.Cov.Ind, SE.pred)		
		names(SE.to.merge)<-c("SE.Cov.Ind","SE.pred")
						}
	if(SE.data$boot.ID[1] == 0){
		list(SE.to.merge, extractAIC(Srch.EF))}
	else {SE.to.merge}

} ## ends Search.Eff.logit

#Calculate CP and RHAT#################################################################################

### Functions for computing effective intervals and r.hat
## effintt.funct = Function to calculate effective interval for any of 4 possible distributions
## Takes data from merged bootstrapped CP, SE and Carcasses
## Input: Interval(i)     = interval preceding carcass find (int)
##        CP.pred(c)      = estimated 1/lambda = average carcass persistence time (t-bar-hat) for exponential, but just 1/lambda for weibull  
##        shape           = shape parameter from fit of weibull, loglogistic or lognormal
effint.fnct <- function(c, shape){
     switch(distn,
            exponential = c*(-log(0.01)),
            weibull     = c*(-log(0.01))^shape,
            loglogistic = c*99^shape, 
            lognormal   = c*exp(qnorm(0.99))^shape  
             )
   }

prob.GTi <- function(i, c, shape){
     switch(distn,
            exponential = exp(-(i/c)),
            weibull     = exp(-(i/c)^(1/shape)),       
            loglogistic = 1/(1+(i/c)^(1/shape)),  
            lognormal   = 1-pnorm(log((i/c)^(1/shape))) 
             )
   }

r.hat   <- function(i,c,shape){integrate(prob.GTi, lower=0, upper=i, c=c ,shape=shape)$value/i}
V.r.hat <- Vectorize(r.hat)

## This function fits a survial model with specified distribution

CP.p.m.d.EST2 <- data.frame(CP.p.m.d.EST)

#CP plots#########################################################################

setwd(dirname(write.path.SECP)) #Start Plotting, Save as jpeg in same folder as SECP
SURV<- Surv(time=Persist.Boot[[1]][,Left.col],time2=Persist.Boot[[1]][,Right.col], event=Persist.Boot[[1]][,event.col],type=c('interval'))
	
	if(if.CPcovs==FALSE){
		jpeg('CPplot.jpeg')
		par(mfrow=c(1,1),oma=c(1,0,0,0), mar=c(6, 4, 4, 2) + 0.1, xpd=F)
		plot(survfit(as.formula(paste("SURV",1, sep="~")), data=Persist.Boot[[1]]), main="Carcass Persistence", ylab="Percent Persisting", xlab="Time")
		
		lines(predict(survreg(as.formula(paste("SURV",1, sep="~")), data=Persist.Boot[[1]], dist=distn),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="blue")
		lines(predict(survreg(as.formula(paste("SURV",1, sep="~")), data=Persist.Boot[[1]], dist="exponential"),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="red")
		mtext(text=bquote(phantom("------ ") * "exponential " * phantom("------ ") * .(distn)), side=1, col="black", outer=TRUE)
		mtext(text=bquote(phantom("------ ") * phantom("exponential ") * "------ " * phantom(.(distn))), side=1, col="blue", outer=TRUE)
		mtext(text=bquote("------ " * phantom("exponential ------ ") * phantom(.(distn))), side=1, col="red", outer=TRUE)
		dev.off()}

	if(if.CPcovs==TRUE){	

		if (how.many.CPcovs<=2){

			if (how.many.CPcovs==1){
				jpeg('CPplot.jpeg', width = 480+100*plotdims$dims[1], height = 480)
				par(mfcol=c(1,plotdims$dims[1]),oma=c(3,3,2,0),mar=c(4, 4, 4, 2) + 0.1, xpd=F)} 
			if (how.many.CPcovs==2){
				jpeg('CPplot.jpeg', width = 480+100*plotdims$dims[1], height = 480+100*plotdims$dims[2])
				par(mfcol=c(plotdims$dims[2],plotdims$dims[1]),oma=c(3,3,2,0),mar=c(4, 4, 4, 2) + 0.1, xpd=F)} 
		
			
			f1 <- function(dat) {
					if((class(try(plot(survfit(as.formula(paste("SURV",1, sep="~")), data=dat), ylim=c(0,1), xlim=c(0,max(CP[,Left.col])), main=dat[,"Persist.Boot..1......CP.Cov.Ind.."][1]), silent=TRUE)) == "try-error") == TRUE){
  					plot(1,1, xlab="", ylab="", type="n", ylim=c(0,1), xlim=c(0,max(CP[,Left.col])))
					text(.5*max(CP[,Left.col]), .5, "Insufficient Data")}
					else{
					lines(predict(survreg(as.formula(paste("SURV",1, sep="~")), data=dat, dist=distn),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="blue")
					lines(predict(survreg(as.formula(paste("SURV",1, sep="~")), data=dat, dist="exponential"),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="red")}}
			by(data.frame(SURV,Persist.Boot[[1]][,"CP.Cov.Ind"]),data.frame(SURV,Persist.Boot[[1]][,"CP.Cov.Ind"])[,"Persist.Boot..1......CP.Cov.Ind.."], f1, simplify= T)
			mtext(text=bquote(phantom("------ ") * "exponential " * phantom("------ ") * .(distn)), side=1, col="black", outer=TRUE)
			mtext(text=bquote(phantom("------ ") * phantom("exponential ") * "------ " * phantom(.(distn))), side=1, col="blue", outer=TRUE)
			mtext(text=bquote("------ " * phantom("exponential ------ ") * phantom(.(distn))), side=1, col="red", outer=TRUE)
				
			
			mtext("Carcass Persistence", line=0, adj=.5, cex=1, outer=TRUE)
			mtext("Percent Persisting", line=0, side=2, adj=.5, cex=1, outer=TRUE)
			dev.off()}

		if (how.many.CPcovs>=3){
				
			for (i in 1:plotdims$dims[3]){
				plotname <- paste('CPplot',unique(CP.p.m.d.EST2[,4])[i],'.jpeg', sep='')
				jpeg(plotname, width = 480+100*plotdims$dims[1], height = 480+100*plotdims$dims[2])
				par(mfcol=c(plotdims$dims[2],plotdims$dims[1]),oma=c(3,3,2,0),mar=c(4, 4, 4, 2) + 0.1, xpd=F)
				Persist.BootSub <-  Persist.Boot[[1]][which(Persist.Boot[[1]][,plotdims$vars[3]]==unique(CP.p.m.d.EST2[,4])[i]),]
				SURVSub <- Surv(time=Persist.BootSub[,Left.col],time2=Persist.BootSub[,Right.col], event=Persist.BootSub[,event.col],type=c('interval'))				


			f1 <- function(dat) {
					if((class(try(plot(survfit(as.formula(paste("SURVSub",1, sep="~")), data=dat), ylim=c(0,1), xlim=c(0,max(CP[,Left.col])), main=dat[,"Persist.BootSub....CP.Cov.IndPlotSub.."][1]), silent=TRUE)) == "try-error") == TRUE){
  					plot(1,1, xlab="", ylab="", type="n", ylim=c(0,1), xlim=c(0,max(CP[,Left.col])), main=dat[,"Persist.BootSub....CP.Cov.IndPlotSub.."][1])
					text(.5*max(CP[,Left.col]), .5, "Insufficient Data")}
					else{
					lines(predict(survreg(as.formula(paste("SURVSub",1, sep="~")), data=dat, dist=distn),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="blue")
					lines(predict(survreg(as.formula(paste("SURVSub",1, sep="~")), data=dat, dist="exponential"),type="quantile",p=seq(.01,.99,by=.01))[1,],seq(.99,.01,by=-.01),col="red")}}
				
			by(data.frame(SURVSub,Persist.BootSub[,"CP.Cov.IndPlotSub"]),data.frame(SURVSub,Persist.BootSub[,"CP.Cov.IndPlotSub"])[,"Persist.BootSub....CP.Cov.IndPlotSub.."], f1, simplify= T)				
			mtext(text=bquote(phantom("------ ") * "exponential " * phantom("------ ") * .(distn)), side=1, col="black", outer=TRUE)
			mtext(text=bquote(phantom("------ ") * phantom("exponential ") * "------ " * phantom(.(distn))), side=1, col="blue", outer=TRUE)
			mtext(text=bquote("------ " * phantom("exponential ------ ") * phantom(.(distn))), side=1, col="red", outer=TRUE)
					
				
			mtext("Carcass Persistence", line=0, adj=.5, cex=1, outer=TRUE)
			mtext("Percent Persisting", line=1, side=2, adj=.5, cex=1, outer=TRUE)
			dev.off()}}}

#End Plotting


Fit.Surv<-function(CP.data, covariates, distn, event.col, Left.col, Right.col, if.CPcovs, CP.p.m.d.EST){
# Create CP estimates. Note CP.est$scale is default output from survreg, but scale represents what is usually known as shape


	# First: make Surv object for R to understand
	SURV<- Surv(time=CP.data[,Left.col],time2=CP.data[,Right.col], event=CP.data[,event.col],type=c('interval'))
	if(if.CPcovs==TRUE){fmla<-as.formula(paste("SURV",covariates, sep="~"))}
	if(if.CPcovs==FALSE){fmla<-as.formula(paste("SURV",1, sep="~"))}
	CP.est<-survreg(fmla, data=CP.data, dist=distn)

	if(if.CPcovs==TRUE){
		CP.pred<-predict(CP.est, CP.p.m.d.EST)
		CP.to.merge<-data.frame(CP.p.m.d.EST$CP.Cov.Ind, CP.pred, CP.est$scale)
		names(CP.to.merge)<-c("CP.Cov.Ind", "CP.pred", "shape") 										}

	if(if.CPcovs==FALSE){
		CP.pred<-predict(CP.est)[1]
		CP.to.merge<-data.frame(CP.p.m.d.EST$CP.Cov.Ind, CP.pred, CP.est$scale)
		names(CP.to.merge)<-c("CP.Cov.Ind", "CP.pred", "shape")
		}

	if(CP.data$boot.ID[1] == 0){
		list(CP.to.merge, extractAIC(CP.est))}
	else {CP.to.merge}
		

} ## ends Fit.Surv


#Combine####################################################################################################

SE.est <- lapply(SearchEF.Boot, Search.Eff.logit, SE.model.covariates=SE.model.covariates, response.col=response.col.SE, if.SEcovs=if.SEcovs,
			SE.p.m.d.EST=SE.p.m.d.EST2)

CP.est <- lapply(Persist.Boot, Fit.Surv, covariates=CP.model.covariates, distn=distn,
	 event.col=event.col, Left.col=Left.col, Right.col=Right.col, if.CPcovs, CP.p.m.d.EST=CP.p.m.d.EST2)

Fatality.Boot[[1]] <- merge(Fatality.Boot[[1]], CP.est$Persist.Trials[[1]] ) 		       #Raw estimates
Fatality.Boot[[1]] <- merge(Fatality.Boot[[1]], SE.est$SearchEF.Trials[[1]] ) 		       #Raw estimates
for (i in 2:(num.boot.samples+1)){Fatality.Boot[[i]]<- merge(Fatality.Boot[[i]], CP.est[[i]])} #Boot estimates
for (i in 2:(num.boot.samples+1)){Fatality.Boot[[i]]<- merge(Fatality.Boot[[i]], SE.est[[i]])} #Boot estimates

f = function(x) function(i) unlist(lapply(x, `[[`, i), use.names=FALSE)
Cas <- as.data.frame(Map(f(Fatality.Boot), names(Fatality.Boot[[1]])))

Cas$Effint  <- pmin(effint.fnct(Cas$CP.pred, shape=Cas$shape) ,Cas$Interval)
Cas$propeff <- Cas$Effint/Cas$Interval											 
Cas$rhat    <- V.r.hat(Cas$Effint, Cas$CP.pred, Cas$shape)

#Calculate AdjFat##############################################################################################

# set min SE.pred value to .05
Cas$AdjFat<-1/(as.numeric(Cas[,propsrch])*Cas$propeff*pmax(0.05,Cas$SE.pred)*Cas$rhat)

### Get AICc values

AIC.values<-data.frame(matrix(c(NA,NA),ncol=2))
names(AIC.values)<-c("SE.AIC","CP.AIC")
	AICSE <- SE.est$SearchEF.Trials[[2]] # vector c(edf, -2*logL) where AIC = -2*logL + k*edf and k=2 by default
	AICCP <- CP.est$Persist.Trials[[2]] #    edf is equiv. df = K in B&A 
	AIC.values[1,1]<-round(AICSE[2] + ((2*AICSE[1]*(AICSE[1]+1))/(dim(SE)[1]-AICSE[1]-1)),2)
	AIC.values[1,2]<-round(AICCP[2] + ((2*AICCP[1]*(AICCP[1]+1))/(dim(CP)[1]-AICCP[1]-1)),2)


#Create Output: BOOT ############################################################################

Max_Excel<-1048576
COLUMN.names.to.keep.boot <- c(Cas.names, "SE.Cov.Ind","CP.Cov.Ind", "SE.pred","CP.pred", "rhat", "shape", "propeff", "AdjFat", "boot.ID")
Boot <- Cas[,COLUMN.names.to.keep.boot]

if (nrow(Boot)<Max_Excel){
write.table(Boot, file=write.path.AD, row.names=FALSE, col.names=TRUE, sep=",")
}

if (nrow(Boot)>=Max_Excel){
cutoff1 <- (Boot[Max_Excel,]$boot.ID)-1
write.table(Boot[which(Boot$boot.ID <= cutoff1),], file=write.path.AD, row.names=FALSE, col.names=TRUE, sep=",")
Boot <- Boot[which(Boot$boot.ID > cutoff1),]

iter <- 0
write.path.ADOver <- paste(c(strsplit(write.path.AD, "")[[1]][1:(length(strsplit(write.path.AD, "")[[1]])-4)], "OverflowBoot1" , ".csv"),  collapse="")

repeat {
iter <- iter + 1
cutoff <- Boot$boot.ID[Max_Excel]-1

if (is.na(cutoff) == FALSE){
write.table(Boot[which(Boot$boot.ID <= cutoff),], file=write.path.ADOver, row.names=FALSE, col.names=TRUE, sep=",")
Boot <- Boot[which(Boot$boot.ID > cutoff),]
write.path.ADOver <- paste(c(strsplit(write.path.ADOver, "")[[1]][1:(length(strsplit(write.path.ADOver, "")[[1]])-5)], (iter+1) , ".csv"),  collapse="")}

if (is.na(cutoff)==TRUE) {
write.table(Boot, file=write.path.ADOver, row.names=FALSE, col.names=TRUE, sep=",")
break}
}}

#Create Ouptut: SECP #############################################################################

#SE Heading ######

Header.a.3<-matrix(paste("User-defined alpha level of", paste(alpha,";",sep=""), "all reported confidence intervals are", paste((1-alpha)*100, "%",sep=""), "confidence intervals based on", num.boot.samples, "bootstraps."),ncol=1)
Header.a.3<-rbind(Header.a.3, matrix(" ", ncol=1))

write.table(Header.a.3, file = write.path.SECP, sep = ",",
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")

if(!identical(coef.list.SE,"DUMMY.SE")){ 
	SE.Title1<-paste("Searcher Efficiency Estimates; AICc for model =",SE.model.covariates, sep=" ")
	SE.Ests.Title<-paste(SE.Title1,AIC.values$SE.AIC,sep=": ")
	}
if(identical(coef.list.SE,"DUMMY.SE")){  
	SE.Ests.Title<-paste("Searcher Efficiency Estimates; AICc for mean-only model",AIC.values$SE.AIC,sep=": ")
	}

#SE.all.Final ##### 

SE.BOOT <- do.call(rbind, lapply(SE.est[-1],data.frame))

if(identical(coef.list.SE,"DUMMY.SE")){ 
SE.all.Final <- data.frame(matrix(as.numeric(),ncol=6))
names(SE.all.Final)<-c("","Found", "Placed", "SE", "Lwr", "Upr")
SE.all.Final[1,1]<-"Mean"
SE.all.Final$Found <- sum(SE[,response.col.SE])                 
SE.all.Final$Placed <- length(SE[,response.col.SE])		   
SE.all.Final$SE <- round(SE.est$SearchEF.Trials[[1]]$SE.pred,2)
SE.all.Final$Lwr <- round(quantile(SE.BOOT$SE.pred, alpha/2),2)
SE.all.Final$Upr <- round(quantile(SE.BOOT$SE.pred, 1-alpha/2),2)
}

if(!identical(coef.list.SE,"DUMMY.SE")){ 
SE.found <- by(SE[,response.col.SE], SE$SE.Cov.Ind, function(dat) {sum(dat==1)}, simplify=T)
SE.all.Final<-data.frame(SE.p.m.d.EST[,-1], matrix(as.numeric(SE.found),ncol=1))
names(SE.all.Final) <- c(coef.list.SE, "Found")
SE.all.Final$Placed <- matrix(as.numeric(by(SE[,response.col.SE], SE$SE.Cov.Ind, function(dat) {length(dat)}, simplify=T)),ncol=1)
SE.all.Final$SE <- round(SE.est$SearchEF.Trials[[1]]$SE.pred,2)
SE.all.Final$Lwr <- matrix(as.numeric(by(SE.BOOT$SE.pred, SE.BOOT$SE.Cov.Ind, function(dat) {round(quantile(dat, alpha/2),2)}, simplify=T)), ncol=1)
SE.all.Final$Upr <- matrix(as.numeric(by(SE.BOOT$SE.pred, SE.BOOT$SE.Cov.Ind, function(dat) {round(quantile(dat, 1-alpha/2),2)}, simplify=T)), ncol=1)
}

add.line<-matrix(rep(" ",ncol(SE.all.Final)), ncol=ncol(SE.all.Final))
SE.all.Final<-rbind(as.matrix(SE.all.Final),add.line)

#Write Table #####

Header3 <- data.frame(SE.Ests.Title)

write.table(Header3, file = write.path.SECP, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")
write.table(SE.all.Final, file = write.path.SECP, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = TRUE,eol = "\n")

#CP Heading ####

if(!identical(coef.list.CP,"DUMMY.CP")){
	CP.Title1<-paste("Carcass Persistence Estimates; model =",CP.model.covariates, sep=" ")
	CP.Title2<-paste(CP.Title1, "and distribution =", distn, "with shape =", round(CP.est$Persist.Trials[[1]]$shape[1],2), sep=" ")
	CP.Ests.Title<-paste(CP.Title2, ", AICc =", AIC.values$CP.AIC,sep=" ")

	}
if(identical(coef.list.CP,"DUMMY.CP")){
	CP.Title1<-paste("Carcass Persistence Estimates; mean-only model and distribution =", distn, "with shape =", round(CP.est$Persist.Trials[[1]]$shape,2), sep=" ")
	CP.Ests.Title<-paste(CP.Title1, ", AICc =", AIC.values$CP.AIC,sep=" ")
}
Header4a<-data.frame(CP.Ests.Title)

#Create CP.all.Final and INT.calc ############

CP.BOOT <- do.call(rbind, lapply(CP.est[-1],data.frame))
CP.est$Persist.Trials[[1]]$Interval <- INT.calc #Mode from original CP data
CP.BOOT$Interval <- INT.calc

CP.BOOT$Effint  <- pmin(effint.fnct(CP.BOOT$CP.pred, shape=CP.BOOT$shape) , CP.BOOT$Interval)
CP.BOOT$propeff <- CP.BOOT$Effint/CP.BOOT$Interval											 
CP.BOOT$rhat    <- V.r.hat(CP.BOOT$Effint, CP.BOOT$CP.pred, CP.BOOT$shape)

CP.est$Persist.Trials[[1]]$Effint  <- pmin(effint.fnct(CP.est$Persist.Trials[[1]]$CP.pred, shape=CP.est$Persist.Trials[[1]]$shape) , CP.est$Persist.Trials[[1]]$Interval)
CP.est$Persist.Trials[[1]]$propeff <- CP.est$Persist.Trials[[1]]$Effint/CP.est$Persist.Trials[[1]]$Interval											 
CP.est$Persist.Trials[[1]]$rhat    <- V.r.hat(CP.est$Persist.Trials[[1]]$Effint, CP.est$Persist.Trials[[1]]$CP.pred, CP.est$Persist.Trials[[1]]$shape)

if(identical(coef.list.CP,"DUMMY.CP")){ 
CP.all.Final <- data.frame(matrix(as.numeric(),ncol=8))
names(CP.all.Final)<-c("", "Placed", "CP", "Lwr", "Upr", "r", "r.Lwr", "r.Upr")
CP.all.Final[1,1]<-"Mean"             
CP.all.Final$Placed <- length(CP[,event.col])		   
CP.all.Final$CP <- round(CP.est$Persist.Trials[[1]]$CP.pred,2)
CP.all.Final$Lwr <- round(quantile(CP.BOOT$CP.pred, alpha/2),2)
CP.all.Final$Upr <- round(quantile(CP.BOOT$CP.pred, 1-alpha/2),2)
CP.all.Final$r <- round(CP.est$Persist.Trials[[1]]$rhat,2)
CP.all.Final$r.Lwr <- round(quantile(CP.BOOT$rhat, alpha/2),2)
CP.all.Final$r.Upr <- round(quantile(CP.BOOT$rhat, 1-alpha/2),2)
}

if(!identical(coef.list.CP,"DUMMY.CP")){ 

CP.placed <- by(CP$CP.Cov.Ind, CP$CP.Cov.Ind, function(dat) {length(dat)}, simplify=T)
CP.all.Final<-data.frame(CP.p.m.d.EST[,-1], matrix(as.numeric(CP.placed),ncol=1))
names(CP.all.Final) <- c(coef.list.CP, "Placed")
CP.all.Final$CP <- round(CP.est$Persist.Trials[[1]]$CP.pred,2)
CP.all.Final$Lwr <- matrix(as.numeric(by(CP.BOOT$CP.pred, CP.BOOT$CP.Cov.Ind, function(dat) {round(quantile(dat, alpha/2),2)}, simplify=T)), ncol=1)
CP.all.Final$Upr <- matrix(as.numeric(by(CP.BOOT$CP.pred, CP.BOOT$CP.Cov.Ind, function(dat) {round(quantile(dat, 1-alpha/2),2)}, simplify=T)), ncol=1)
CP.all.Final$r <- round(CP.est$Persist.Trials[[1]]$rhat,2)
CP.all.Final$r.Lwr <- matrix(as.numeric(by(CP.BOOT$rhat, CP.BOOT$CP.Cov.Ind, function(dat) {round(quantile(dat, alpha/2),2)}, simplify=T)), ncol=1)
CP.all.Final$r.Upr <- matrix(as.numeric(by(CP.BOOT$rhat, CP.BOOT$CP.Cov.Ind, function(dat) {round(quantile(dat, 1-alpha/2),2)}, simplify=T)), ncol=1)
}


#Write Tables #######

Header4b<-data.frame(paste("Estimates and Confidence Intervals for r based on an interval of",INT.calc, "days.",sep= " "))

write.table(Header4a, file = write.path.SECP, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")
write.table(Header4b, file = write.path.SECP, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")
write.table(CP.all.Final, file = write.path.SECP, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = TRUE,eol = "\n")


#Create Output: FAT ###############################################################################################################

### Create Output.Heading, which has low boot # warning(if necessary), # turbines in study area vs # sampled1

NTiS<-matrix(NA, ncol=1)
if(Num.Turbines.in.Study<Num.Total.Turbines.Site) NTiS[1,]<-paste("This study included",Num.Turbines.in.Study, "turbines, and these estimates have been expanded to reference the", Num.Total.Turbines.Site, "at the study site.")
if(Num.Turbines.in.Study==Num.Total.Turbines.Site) NTiS[1,]<-paste("This study included all",Num.Turbines.in.Study, "turbines located at the study site.")

if(num.boot.samples<1000){
	bootstrap.N.warning<-matrix(paste("Results based on ", num.boot.samples, " bootstrap resamples: STANDARD ERROR NOT CONSIDERED VALID"), ncol=1)
	Output.Heading<-rbind(bootstrap.N.warning, NTiS)}# ends if num.boot.samples<1000
if(num.boot.samples>=1000){
	Output.Heading<-NTiS}# ends if num.boot.samples>=1000

Output.Heading<-rbind(Output.Heading, paste("User-defined alpha level of", paste(alpha,";",sep=""), "all reported confidence intervals are", paste((1-alpha)*100, "%",sep=""), "confidence intervals based on", num.boot.samples, "bootstraps."))
Output.Heading<-rbind(Output.Heading," ")

### Create Fat Table (AllReqs.Final)

AllReqs.Final <- data.frame(matrix(as.numeric(),ncol=9))
names(AllReqs.Final)<-c("Group","Level","Number.Found","PerTurb.Est","PerTurb.Lwr","PerTurb.Upr","SiteTot.Est","SiteTot.Lwr","SiteTot.Upr")
AllReqs.Final[1,1]<-"Overall"

AllReqs.Final$Number.Found <- nrow(Cas[which(Cas$boot.ID == 0), ])
AllReqs.Final$Level <- ""
AllReqs.Final$PerTurb.Est <- round(sum(Cas[which(Cas$boot.ID == 0), ]$AdjFat)/Num.Turbines.in.Study,2)			
AllReqs.Final$SiteTot.Est <- ceiling(Num.Total.Turbines.Site*AllReqs.Final$PerTurb.Est)
AllReqs.Final$PerTurb.Lwr <- round(quantile(with(Cas[which(Cas$boot.ID != 0),], aggregate(AdjFat, by=list(boot.ID), sum))[,2]/Num.Turbines.in.Study,alpha/2),2)
AllReqs.Final$PerTurb.Upr <- round(quantile(with(Cas[which(Cas$boot.ID != 0),], aggregate(AdjFat, by=list(boot.ID), sum))[,2]/Num.Turbines.in.Study,1-alpha/2),2)
AllReqs.Final$SiteTot.Lwr <- floor(Num.Total.Turbines.Site*AllReqs.Final$PerTurb.Lwr) 				####FIXED To Floor
AllReqs.Final$SiteTot.Upr <- ceiling(Num.Total.Turbines.Site*AllReqs.Final$PerTurb.Upr)

add.line <- data.frame(matrix(c(" "," "," "," "," "," "," "," ", " "), ncol=9))
names(add.line)<-names(AllReqs.Final)
AllReqs.Final <- rbind(AllReqs.Final, add.line)

if(!identical(summary.list,NA)){ 

for (i in 1:how.many.reqs){
num.levels <- length(unique(Cas[,summary.list[i]]))
Cas <- Cas[order(Cas[,summary.list[i]]),]
List1 <- data.frame(matrix(as.numeric(),ncol=9))
names(List1)<-c("Group","Level","Number.Found","PerTurb.Est","PerTurb.Lwr","PerTurb.Upr","SiteTot.Est","SiteTot.Lwr","SiteTot.Upr")

for (j in 1:num.levels) {
List1[j,"Group"] <- rep(summary.list[i])
List1[j,"Level"]  <- as.character(unique(Cas[,summary.list[i]])[j])
List1[j,"Number.Found"] <- nrow(Cas[which(Cas$boot.ID == 0 & Cas[,summary.list[i]] == List1[j,"Level"]), ])
List1[j,"PerTurb.Est"] <- round(sum(Cas[which(Cas$boot.ID == 0 & Cas[,summary.list[i]] == List1[j,"Level"]), ]$AdjFat)/Num.Turbines.in.Study,2)			
List1[j,"SiteTot.Est"] <- ceiling(Num.Total.Turbines.Site*List1[j,"PerTurb.Est"])
List1[j,"PerTurb.Lwr"] <- round(quantile(with(Cas[which(Cas$boot.ID != 0 & Cas[,summary.list[i]] == List1[j,"Level"]),], aggregate(AdjFat, by=list(boot.ID), sum))[,2]/Num.Turbines.in.Study,alpha/2),2)
List1[j,"PerTurb.Upr"] <- round(quantile(with(Cas[which(Cas$boot.ID != 0 & Cas[,summary.list[i]] == List1[j,"Level"]),], aggregate(AdjFat, by=list(boot.ID), sum))[,2]/Num.Turbines.in.Study,1-alpha/2),2)
List1[j,"SiteTot.Lwr"] <- floor(Num.Total.Turbines.Site*List1[j,"PerTurb.Lwr"]) 				
List1[j,"SiteTot.Upr"] <- ceiling(Num.Total.Turbines.Site*List1[j,"PerTurb.Upr"]) }

add.line <- data.frame(matrix(c(" "," "," "," "," "," "," "," ", " "), ncol=9))
names(add.line)<-names(List1)
List1 <- rbind(List1, add.line)
AllReqs.Final <- rbind(AllReqs.Final, List1)}}

### Write Table 

if (sum(as.numeric(AllReqs.Final$Number.Found)>as.numeric(AllReqs.Final$SiteTot.Lwr),na.rm=T)>0){
	LowerBound.warning<-matrix(paste("Estimates for Site Total Lower Bound are less than Number Found at Site, revise interval accordingly"), ncol=1)
	Output.Heading1<-rbind(Output.Heading,LowerBound.warning)
} else {Output.Heading1<-Output.Heading}


if(min(as.numeric(AllReqs.Final$Number.Found), na.rm=TRUE)<5){Output.Heading1<-rbind(Output.Heading1,"The number of observed casualties for some requested groups is less than 5, use caution when interpreting estimates."," ")}  

write.table(Output.Heading1, file = write.path.Final, sep = ",",
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")
Header1 <- data.frame("Requested Summaries")
write.table(Header1, file = write.path.Final, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = FALSE,eol = "\n")
write.table(AllReqs.Final, file = write.path.Final, sep = ",",append=TRUE,
            na = "NA", dec = ".", row.names = FALSE, col.names = TRUE,eol = "\n")

}} # ends if mode(RunIT)!="character")

#Allows the error messages to print
if(mode(RunIT)=="character") {
	Final.Output<-RunIT
	write.table(Final.Output, file=write.path.SECP,row.names=FALSE, col.names=FALSE)
	write.table(Final.Output, file=write.path.AD,row.names=FALSE, col.names=FALSE)
	write.table(Final.Output, file=write.path.Final,row.names=FALSE, col.names=FALSE)

} # ends if mode(RunIT)=="character"