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



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



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



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



###################################################
### chunk number 5: figure
###################################################
#line 79 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Plot streams in Indiana
strahler <- levels(IN_streams$Strahler_order)
cols <- rainbow(4)
plot(IN_streams$xcoord, IN_streams$ycoord, type="n", xlab="x-coordinate",
     ylab="y-coordinate")
for(i in 1:4) {
   ind <- IN_streams$Strahler_order == strahler[i]
   points(IN_streams$xcoord[ind], IN_streams$ycoord[ind], pch=20, col=cols[i])
   }
legend(x="topright", inset=0.05, legend=paste("Order", strahler), pch=20, cex=1,
       col=cols)
title("Plot of Indiana Stream Sites Color-Coded by Strahler Order")



###################################################
### chunk number 6: Statuseval
###################################################
#line 106 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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(IN_streams$status))



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



###################################################
### chunk number 8: Statuseval
###################################################
#line 134 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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 streams in each category
sites <- data.frame(siteID=IN_streams$siteID,
                    Use=rep(TRUE, nr))



###################################################
### chunk number 9: Statuseval
###################################################
#line 152 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Create the subpop data frame, which defines populations and subpopulations for
# which estimates are desired
subpop <- data.frame(siteID=IN_streams$siteID,
                     Upper_Wabash=rep("Upper Wabash", nr), 
							       Strahler_Order=IN_streams$Strahler_order)



###################################################
### chunk number 10: Statuseval
###################################################
#line 168 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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=IN_streams$siteID,
                     wgt=IN_streams$wgt,
                     xcoord=IN_streams$xcoord,
                     ycoord=IN_streams$ycoord)



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



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



###################################################
### chunk number 13: Statuseval
###################################################
#line 201 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Print the extent estimates for all basins combined
print(Extent_Estimates[c(1:8, 32:34),])



###################################################
### chunk number 14: Statuseval
###################################################
#line 212 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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 15: Conditioneval
###################################################
#line 226 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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(IN_streams$IBI_status))



###################################################
### chunk number 16: Conditioneval
###################################################
#line 238 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Use the table and addmargins functions to create a table displaying the count
# for each code of the QHEI status variable
cat("\nA table displaying the number of values for each level of the QHEI status variable follows:\n")
addmargins(table(IN_streams$QHEI_status))



###################################################
### chunk number 17: Conditioneval
###################################################
#line 254 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
#
# Conduct an analysis of stream condition variables
#

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

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



###################################################
### chunk number 18: Conditioneval
###################################################
#line 274 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Create the data.cat data frame, which specifies the variables to use in the
# analysis
data.cat <- data.frame(siteID=IN_streams$siteID,
                       IBI_Status=IN_streams$IBI_status,
                       QHEI_Status=IN_streams$QHEI_status)



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



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



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



###################################################
### chunk number 22: Conditionevalpop
###################################################
#line 318 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
#
# Conduct an analysis of stream 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("1"=4514.450, "2"=1443.260, "3"=740.146, "4"=660.294)



###################################################
### chunk number 23: Conditionevalpop
###################################################
#line 338 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Calculate estimates for the stream condition variables
Condition_Estimates_popsize <- cat.analysis(sites, subpop, design, data.cat,
   popsize=list(Upper_Wabash=sum(framesize),
                Strahler_Order=as.list(framesize)))



###################################################
### chunk number 24: Conditionevalpop
###################################################
#line 347 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Print the stream condition estimates for all sites combined
print(Condition_Estimates_popsize[c(1:3, 16:19),])



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



###################################################
### chunk number 26: Quanteval
###################################################
#line 372 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_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(IN_streams$IBI_score)



###################################################
### chunk number 27: Quanteval
###################################################
#line 383 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Use the summary function to summarize the data structure of the QHEI score
# variable
cat("\nSummarize the data structure of the QHEI score variable:\n")
summary(IN_streams$QHEI_score)



###################################################
### chunk number 28: Quanteval
###################################################
#line 399 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
#
# Conduct an analysis of stream 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=IN_streams$siteID,
                        IBI_Score=IN_streams$IBI_score,
                        QHEI_Score=IN_streams$QHEI_score)



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



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



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



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



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



###################################################
### chunk number 34: Quanteval
###################################################
#line 478 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Test for statistical difference between CDFs for Strahler order categories
CDF_Tests <- cont.cdftest(sites, subpop[,c(1,3)], design, data.cont,
   popsize=list(Strahler_Order=as.list(framesize)))



###################################################
### chunk number 35: Quanteval
###################################################
#line 486 "d:/Rcompile/CRANpkg/local/2.12/spsurvey/inst/doc/Linear_Analysis.Rnw"
# Print results of the statistical tests for difference between CDFs from
# Strahler order categories for IBI score
print(CDF_Tests[1:15,])



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



