###################################################
### chunk number 1: preliminaries
###################################################
#line 26 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Load the spsurvey package
library(spsurvey)



###################################################
### chunk number 2: data
###################################################
#line 40 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Read the data file and determine the number of rows in the file
SC_estuaries <- read.delim("SC_estuaries.tab")
nr <- nrow(SC_estuaries)



###################################################
### chunk number 3: data
###################################################
#line 48 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Display the initial six lines in the data file
head(SC_estuaries)



###################################################
### chunk number 4: data
###################################################
#line 55 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Display the final six lines in the data file
tail(SC_estuaries)



###################################################
### chunk number 5: figure
###################################################
#line 78 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Plot estuaries in Indiana
stratum <- levels(SC_estuaries$stratum)
cols <- rainbow(2)
plot(SC_estuaries$xcoord, SC_estuaries$ycoord, type="n", xlab="x-coordinate",
     ylab="y-coordinate")
for(i in 1:2) {
   ind <- SC_estuaries$stratum == stratum[i]
   points(SC_estuaries$xcoord[ind], SC_estuaries$ycoord[ind], pch=20,
          col=cols[i])
   }
legend(x="bottomright", inset=0.05, legend=stratum, pch=20, cex=0.8, col=cols)
title("Plot of South Carolina Estuary Sites Color-Coded by Stratum")



###################################################
### chunk number 6: Statuseval
###################################################
#line 105 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Use the table and addmargins functions to create a table displaying the count
# for each code of the status variable
cat("\nA table displaying the number of values for each level of the status
variable follows:\n")
addmargins(table(SC_estuaries$status))



###################################################
### chunk number 7: Statuseval
###################################################
#line 121 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
#
# Conduct an analysis of site status evaluation variables
#

# Create the sites data frame, which identifies sites to use in the analysis
# Note that all sites will be used to estimate number of estuaries in each category
sites <- data.frame(siteID=SC_estuaries$siteID,
                    Use=rep(TRUE, nr))



###################################################
### chunk number 8: Statuseval
###################################################
#line 139 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Create the subpop data frame, which defines populations and subpopulations for
# which estimates are desired
subpop <- data.frame(siteID=SC_estuaries$siteID,
                     All_Estuaries=rep("All Estuaries", nr), 
							       Estuary_Type=SC_estuaries$stratum)



###################################################
### chunk number 9: Statuseval
###################################################
#line 155 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Create the design data frame, which identifies the stratum code, weight,
#    x-coordinate, and y-coordinate for each site ID
design <- data.frame(siteID=SC_estuaries$siteID,
                     wgt=SC_estuaries$wgt,
                     xcoord=SC_estuaries$xcoord,
                     ycoord=SC_estuaries$ycoord)



###################################################
### chunk number 10: Statuseval
###################################################
#line 170 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Create the data.cat data frame, which specifies the variables to use in the
# analysis
data.cat <- data.frame(siteID=SC_estuaries$siteID,
                       Status=SC_estuaries$status)



###################################################
### chunk number 11: Statuseval
###################################################
#line 179 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Calculate extent estimates for the site status evaluation variables
Extent_Estimates <- cat.analysis(sites, subpop, design, data.cat)



###################################################
### chunk number 12: Statuseval
###################################################
#line 186 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Print the extent estimates
print(Extent_Estimates)



###################################################
### chunk number 13: Statuseval
###################################################
#line 197 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write results as a comma-separated value (csv) file
write.table(Extent_Estimates, file="Extent_Estimates.csv", sep=",",
            row.names=FALSE)



###################################################
### chunk number 14: Conditioneval
###################################################
#line 211 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Use the table and addmargins functions to create a table displaying the count
# for each code of the IBI status variable
cat("\nA table displaying the number of values for each level of the IBI status
variable follows:\n")
addmargins(table(SC_estuaries$IBI_status))



###################################################
### chunk number 15: Conditioneval
###################################################
#line 223 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Use the table and addmargins functions to create a table displaying the count
# for each code of the WQ status variable
cat("\nA table displaying the number of values for each level of the WQ status variable follows:\n")
addmargins(table(SC_estuaries$WQ_status))



###################################################
### chunk number 16: Conditioneval
###################################################
#line 239 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
#
# Conduct an analysis of estuary condition variables
#

# Create the sites data frame
# Note that only sampled sites are used
sites <- data.frame(siteID=SC_estuaries$siteID,
                    Use=SC_estuaries$status == "Sampled")

# Note that the existing subpop and design data frames can be reused



###################################################
### chunk number 17: Conditioneval
###################################################
#line 259 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Create the data.cat data frame, which specifies the variables to use in the
# analysis
data.cat <- data.frame(siteID=SC_estuaries$siteID,
                       IBI_Status=SC_estuaries$IBI_status,
                       WQ_Status=SC_estuaries$WQ_status)



###################################################
### chunk number 18: Conditioneval
###################################################
#line 269 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Calculate estimates for the categorical variables
Condition_Estimates <- cat.analysis(sites, subpop, design, data.cat)



###################################################
### chunk number 19: Conditioneval
###################################################
#line 276 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Print the condition estimates for all basins combined
print(Condition_Estimates[c(1:3, 16:19),])



###################################################
### chunk number 20: Conditioneval
###################################################
#line 287 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write results as a csv file
write.table(Condition_Estimates, file="Condition_Estimates.csv", sep=",", 
            row.names=FALSE)



###################################################
### chunk number 21: Conditionevalpop
###################################################
#line 303 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
#
# Conduct an analysis of estuary condition variables correcting for population
# size
#

# Note that the existing sites, subpop, design, and data.cont data frames can be
# reused

# Assign frame size values
framesize <- c("Open Water"=628.509298, "Tidal Creek"=105.829522)



###################################################
### chunk number 22: Conditionevalpop
###################################################
#line 323 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Calculate estimates for the estuary condition variables
Condition_Estimates_popsize <- cat.analysis(sites, subpop, design, data.cat,
   popsize=list(All_Estuaries=sum(framesize),
                Estuary_Type=as.list(framesize)))



###################################################
### chunk number 23: Conditionevalpop
###################################################
#line 332 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Print the estuary condition estimates for all sites combined
print(Condition_Estimates_popsize[1:12,])



###################################################
### chunk number 24: Conditionevalpop
###################################################
#line 343 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write results as a csv file
write.table(Condition_Estimates_popsize, file="Condition_Estimates_popsize.csv",
            sep=",", row.names=FALSE)



###################################################
### chunk number 25: Quanteval
###################################################
#line 357 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Use the summary function to summarize the data structure of the IBI score
# variable
cat("\nSummarize the data structure of the IBI score variable:\n")
summary(SC_estuaries$IBI_score)



###################################################
### chunk number 26: Quanteval
###################################################
#line 368 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Use the summary function to summarize the data structure of the WQ score
# variable
cat("\nSummarize the data structure of the WQ score variable:\n")
summary(SC_estuaries$WQ_score)



###################################################
### chunk number 27: Quanteval
###################################################
#line 384 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
#
# Conduct an analysis of estuary condition variables
#

# Note that the existing sites, subpop, and design data frames can be reused

# Create the data.cat data frame, which specifies the variables to use in the
# analysis
data.cont <- data.frame(siteID=SC_estuaries$siteID,
                        IBI_Score=SC_estuaries$IBI_score,
                        WQ_Score=SC_estuaries$WQ_score)



###################################################
### chunk number 28: Quanteval
###################################################
#line 406 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Calculate CDF estimates for the quantitative variables
CDF_Estimates <- cont.analysis(sites, subpop, design, data.cont,
   popsize=list(All_Estuaries=sum(framesize),
                Estuary_Type=as.list(framesize)))



###################################################
### chunk number 29: Quanteval
###################################################
#line 422 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write CDF estimates as a csv file
write.table(CDF_Estimates$CDF, file="CDF_Estimates.csv", sep=",",
            row.names=FALSE)



###################################################
### chunk number 30: Quanteval
###################################################
#line 432 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
cont.cdfplot("CDF_Estimates.pdf", CDF_Estimates$CDF)



###################################################
### chunk number 31: Quanteval
###################################################
#line 438 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Print the percentile estimates for IBI score for all sites combined
print(CDF_Estimates$Pct[1:10,])



###################################################
### chunk number 32: Quanteval
###################################################
#line 450 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write percentile estimates as a csv file
write.table(CDF_Estimates$Pct, file="Percentile_Estimates.csv", sep=",",
            row.names=FALSE)



###################################################
### chunk number 33: Quanteval
###################################################
#line 463 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Test for statistical difference between CDFs for strata
CDF_Tests <- cont.cdftest(sites, subpop[,c(1,3)], design, data.cont,
   popsize=list(Estuary_Type=as.list(framesize)))



###################################################
### chunk number 34: Quanteval
###################################################
#line 471 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Print results of the statistical tests for difference between strata CDFs for
# IBI score and WQ score
print(CDF_Tests)



###################################################
### chunk number 35: Quanteval
###################################################
#line 479 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Area_Analysis.Rnw"
# Write CDF test results as a csv file
write.table(CDF_Tests, file="CDF_Tests.csv", sep=",", row.names=FALSE)



