# File: estimation.txt
# Purpose: Estimation for Estuaries Example with Benthic IBI Status
# Programmer:  Tony Olsen
# Date:  March 2005
#
#Example updated for R 20.1 and psurvey.analysis 2.6
#
# Read in data file (designstatus) from weight adjustment for Open Water & Tidal Creeks into DesignStatus
# Note that designstatus file contains the data results as well as design status information

DesignStatus <- read.csv('designstatus.csv')
head(DesignStatus)

# Do estimates for IBI and WQ status
# select which sites to use - only TS-target sampled sites
# Set up sites, subpop, dsgn, and data data frames for estimation functions

# sites identifies which sites to use in estimation

# sites must have siteID and column indicating if site will be used
sites <- data.frame(siteID=DesignStatus$Site.id, 
        Use=DesignStatus$EE.Status=='TS')

# Define subpopulations for which estimates are desired
subpop <- data.frame(siteID=DesignStatus$Site.id, 
				EstuaryType=DesignStatus$Type)

# Specify design information
dsgn <- data.frame(siteID=DesignStatus$Site.id,
						stratum=DesignStatus$Type, 
						wgt=DesignStatus$final.wt,
						xcoord=DesignStatus$xmarinus,
						ycoord=DesignStatus$ymarinus)

# Specify data dataframe
data <- DesignStatus[,c('Site.id','IBI.Status','WQ.Status')]
names(data) <- c('siteID','IBIStatus','WQStatus')


# Calculate the estimates: note no longer know the size of the target stream length
IBIWQSum <- cat.analysis(sites,subpop,dsgn,data)

# Write results out as csv file to read in spreadsheet
write.table(IBIWQSum, file = "ow_ibiwqsum.csv", sep = ",",col.names=NA)


# Complete continuous indicator estimates 
# Can use same site, subpop, and dsgn data.frames
# Set up the data data frame
data <- DesignStatus[,c('Site.id','EE.IBI','WQ.IS')]
dimnames(data)[[2]] <- c('SiteID','IBI Score','WQ Score')

ScoresSum <- cont.analysis(sites,subpop,dsgn,data,
			pctval=c(5, 10, 25, 50, 75, 90, 95),
			conf=95 ,vartype='Local')

# Write results out as csv file to read in spreadsheet
write.table(ScoresSum$CDF, file = "scoressumcdf.csv", sep = ",",col.names=NA)
write.table(ScoresSum$Pct, file = "scoressumpct.csv", sep = ",",col.names=NA)
write.table(ScoresSum$Tot, file = "scoressumtot.csv", sep = ",",col.names=NA)



# Can use R to Plot CDF estimates using function: cdfplot.fcn
#	File: cdfplot.fcn.r
#	Programmer: Tony Olsen
#	Date: Sept 23, 2002
# Input:
#	cdfest - dataframe with
#			x-value for cdf plot
#			y-value (cdf estimate)
#			lower confidence bound
#			upper confidence bound
#
# Output:
#	plot of cdf with confidence bounds	
cdfplot.fcn <- function(cdfest,prop=T,...){

if(prop == T) {
	plot(cdfest[,1],cdfest[,2]  ,
		type='l',ylim=c(0,100),...)
	tvalue <- cdfest[,2]>=5 & cdfest[,2]<=95
}
else {
	plot(cdfest[,1],cdfest[,2],type='l',...)
	tvalue <- cdfest[,2]>=0.05*max(cdfest[,2]) & 
			cdfest[,2]<=0.95*max(cdfest[,2])
}
value <- cdfest[,1][tvalue]
upper <- cdfest[,4][tvalue]
lower <- cdfest[,3][tvalue]
lines(value,lower,lty=2)
lines(value,upper,lty=2)

legend(x=min(cdfest[,1]),y=max(cdfest[,2]),
	legend=c('CDF estimate','95% Confidence Limits'),
	lty=c(1,2), bty='n', cex=.7)
}

pdf('results.pdf',width=8,height=11.5)
op <- par(mfrow=c(2,2))

# select Open water cdf estimates
tsttype <- ScoresSum$CDF$Subpopulation == 'Open Water'
# Plot IBI Score CDFs for percent and km2 area
tst <- tsttype & ScoresSum$CDF$Indicator=='IBI Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.P','LCB95Pct.P','UCB95Pct.P')]
xlab <- 'IBI Score'
ylab <- 'Percent Open Water Area km2'
cdfplot.fcn(cdf,prop=T,xlab=xlab,ylab=ylab)
title('Open Water Area: IBI Score CDF Estimate')

tst <- tsttype & ScoresSum$CDF$Indicator=='IBI Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.U','LCB95Pct.U','UCB95Pct.U')]
xlab <- 'IBI Score'
ylab <- 'Open Water Area km2'
cdfplot.fcn(cdf,prop=F,xlab=xlab,ylab=ylab)
title('Open Water Area: IBI Score CDF Este')

# Plot WQ Score CDFs for percent and km2 area
tst <- tsttype & ScoresSum$CDF$Indicator=='WQ Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.P','LCB95Pct.P','UCB95Pct.P')]
xlab <- 'WQ Score'
ylab <- 'Percent Open Water Area km2'
cdfplot.fcn(cdf,prop=T,xlab=xlab,ylab=ylab)
title('Open Water Area: WQ Score CDF Est')

tst <- tsttype & ScoresSum$CDF$Indicator=='WQ Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.U','LCB95Pct.U','UCB95Pct.U')]
xlab <- 'WQ Score'
ylab <- 'Open Water Area km2'
cdfplot.fcn(cdf,prop=F,xlab=xlab,ylab=ylab)
title('Open Water Area: WQ Score CDF Est')

# select Tidal Creek cdf estimates
tsttype <- ScoresSum$CDF$Subpopulation == 'Tidal Creek'
# Plot IBI Score CDFs for percent and km2 area
tst <- tsttype & ScoresSum$CDF$Indicator=='IBI Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.P','LCB95Pct.P','UCB95Pct.P')]
xlab <- 'IBI Score'
ylab <- 'Percent Tidal Creek Area km2'
cdfplot.fcn(cdf,prop=T,xlab=xlab,ylab=ylab)
title('Tidal Creek Area: IBI Score CDF Estimate')

tst <- tsttype & ScoresSum$CDF$Indicator=='IBI Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.U','LCB95Pct.U','UCB95Pct.U')]
xlab <- 'IBI Score'
ylab <- 'Tidal Creek Area km2'
cdfplot.fcn(cdf,prop=F,xlab=xlab,ylab=ylab)
title('Tidal Creek Area: IBI Score CDF Este')

# Plot WQ Score CDFs for percent and km2 area
tst <- tsttype & ScoresSum$CDF$Indicator=='WQ Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.P','LCB95Pct.P','UCB95Pct.P')]
xlab <- 'WQ Score'
ylab <- 'Percent Tidal Creek Area km2'
cdfplot.fcn(cdf,prop=T,xlab=xlab,ylab=ylab)
title('Tidal Creek Area: WQ Score CDF Est')

tst <- tsttype & ScoresSum$CDF$Indicator=='WQ Score'
cdf <- ScoresSum$CDF[tst,c('Value','Estimate.U','LCB95Pct.U','UCB95Pct.U')]
xlab <- 'WQ Score'
ylab <- 'Tidal Creek Area km2'
cdfplot.fcn(cdf,prop=F,xlab=xlab,ylab=ylab)
title('Tidal Creek Area: WQ Score CDF Est')


par(op)
graphics.off()