# -----------------
# Danelle Haake
# Volunteer Data Analysis
# First written: August 10,2020
# Last revised: April 7, 2021
# -----------------

#### LOCATING THE DATASET ####
setwd("C://Users/strea/Documents/My Publications/Post Dissertation 3-Cit Sci/stats")
#import file "Chloride.csv"
#import file "Invertebrates.csv"
#import file "Watersheds.csv"

#### BASIC SUMMARY STATS AND DATA MANAGEMENT ####

attach(Chloride)

length(ChlorideMgL)   #2234 Cl measurements
ClMean = tapply(ChlorideMgL, SiteNo, mean, na.rm=T)
ClMax = tapply(ChlorideMgL, SiteNo, max, na.rm=T)
ClMed = tapply(ChlorideMgL, SiteNo, median, na.rm=T)
ClIQR = tapply(ChlorideMgL, SiteNo, IQR, na.rm=T)
ClMin = tapply(ChlorideMgL, SiteNo, min, na.rm=T)
ClN = tapply(ChlorideMgL, SiteNo, length)
Cl25th = tapply(ChlorideMgL, SiteNo, quantile, prob=0.25, na.rm=T)
Cl75th = tapply(ChlorideMgL, SiteNo, quantile, prob=0.75, na.rm=T)

DataSummary1 = data.frame(ClN, ClMax, Cl75th, ClMed, Cl25th, ClMin, ClMean, ClIQR)
write.csv(DataSummary1, file="ChlorideSummary.csv") 
summary(DataSummary1)

# colors below to match WQ Rating in next figure 
boxplot(Chloride$ChlorideMgL~Chloride$SiteNo, log="y", ylim=c(10,10000),
        ylab="Chloride (mg/l)", xlab="Site Number",
        cex.lab=1.4, col=c("tomato","tomato","bisque","lightblue1","lightblue1","lightblue1","bisque",
              "tomato","tomato","lightblue1","bisque","bisque","tomato",
              "lightblue1","bisque","bisque","tomato","lightblue1","tomato","bisque",
              "lightblue1","tomato","bisque","bisque","tomato","tomato","tomato",
              "tomato","lightblue1","tomato","skyblue2","tomato","bisque"))
abline(h=230, col="orange", lwd=2.5)
abline(h=860, col="red", lty=2, lwd=2.5)
    # save as Cl_Boxplot.jpeg


attach(Invertebrates)

length(Invertebrates$EVENTDATE)   #218 monitoring events
df1 = data.frame(CADDIS, HELLGRA, MAYFLY, GILLSN, RIFLBEE, STONE, WATRPEN, OTHRBEE, CLAM, CRANEFL,
                 CRAYFSH, DRAGON, DAMSEL, SCUD, SOWBUG, FISHFLY, ALDRFLY, WATRSNI, AQUAWRM, 
                 BLACKFL, LEECH, MIDGE, POUCHSN, OTHERSN)
Invertebrates$TOTALINV = rowSums(df1, na.rm=T)
sum(Invertebrates$TOTALINV)   # 28732 invertebrates collected
df2 = data.frame(CADDIS, MAYFLY, STONE)
Invertebrates$TOTALEPT = rowSums(df2, na.rm=T)
Invertebrates$CPUE = Invertebrates$TOTALINV / TIMEPICK
InvertMean = tapply(Invertebrates$TOTALINV, SITENO, mean, na.rm=T)
InvertMax = tapply(Invertebrates$TOTALINV, SITENO, max, na.rm=T)
CPUEMax = tapply(Invertebrates$CPUE, SITENO, max, na.rm=T)
CPUEMean = tapply(Invertebrates$CPUE, SITENO, mean, na.rm=T)
WQMean = tapply(WQRATING, SITENO, mean, na.rm=T)
WQMax = tapply(WQRATING, SITENO, max, na.rm=T)
WQN = tapply(Invertebrates$TOTALINV, SITENO, length)

DataSummary2 = data.frame(WQN, WQMax, WQMean, InvertMax, InvertMean, CPUEMax, CPUEMean)
DataSummary2[9,6] = NA # change Site 17 CPUE Max from "-Inf" to "NA"
DataSummary2[9,7] = NA # change Site 17 CPUE Max from "NaN" to "NA"
write.csv(DataSummary2, file="InvertebrateSummary.csv") 
summary(DataSummary2)

# above 23 considered excellent, 18-23 considered good, 12-17 considered fair, and ratings below 12 considered poor
par(mfrow=c(2,1), oma = c(3,2.5,0,0), mar = c(3,4,1,1))
boxplot(WQRATING~SITENO, data=Invertebrates, ylab="Water Quality Rating", xlab="Site Number",
        cex.lab=1.4, col=c("tomato","tomato","bisque","lightblue1","lightblue1","lightblue1","bisque",
              "tomato","tomato","lightblue1","bisque","bisque","tomato",
              "lightblue1","bisque","bisque","tomato","lightblue1","tomato","bisque",
              "lightblue1","tomato","bisque","bisque","tomato","tomato","tomato",
              "tomato","lightblue1","tomato","skyblue2","tomato","bisque"))
    # save as WQrating_Boxplot.jpeg
abline(h=23)
abline(h=18)
abline(h=12)

boxplot(TOTALINV~SITENO, data=Invertebrates, cex.axis=0.75, ylab="Total Invertebrates", xlab="Site Number",
        col=c("tomato","skyblue2","bisque","bisque","bisque","bisque","tomato",
              "bisque","red","skyblue2","skyblue2","skyblue2","tomato",
              "bisque","lightblue1","bisque","tomato","lightblue1","tomato","bisque",
              "skyblue2","skyblue2","lightblue1","tomato","tomato","tomato","lightblue1",
              "lightblue1","lightblue1","bisque","skyblue2","skyblue2","skyblue2"))

boxplot(TOTALEPT~SITENO, data=Invertebrates, cex.axis=0.75, ylab="EPT", xlab="Site Number",
        col=c("tomato","skyblue2","bisque","bisque","bisque","bisque","tomato",
              "bisque","red","skyblue2","skyblue2","skyblue2","tomato",
              "bisque","lightblue1","bisque","tomato","lightblue1","tomato","bisque",
              "skyblue2","skyblue2","lightblue1","tomato","tomato","tomato","lightblue1",
              "lightblue1","lightblue1","tomato","skyblue2","skyblue2","skyblue2"))

boxplot(CPUE~SITENO, data=Invertebrates, cex.axis=0.75, ylab="CPUE", xlab="Site Number",
        col=c("tomato","skyblue2","bisque","bisque","bisque","bisque","tomato",
              "bisque","red","skyblue2","skyblue2","skyblue2","tomato",
              "bisque","lightblue1","bisque","tomato","lightblue1","tomato","bisque",
              "skyblue2","skyblue2","lightblue1","tomato","tomato","tomato","lightblue1",
              "lightblue1","lightblue1","tomato","skyblue2","skyblue2","skyblue2"))


attach(Watersheds)

summary(Watersheds)


#### GENERATING FULL SUMMARY DATA TABLE ####
Metrics = cbind(Watersheds, DataSummary1, DataSummary2)
write.csv(Metrics, file="Metrics.csv") 


#### CORRELATIONS ####
library(corrplot)

CorrelationData = Metrics[, c(3:9, 11:17, 19:24)]
colnames(CorrelationData) <- c("Area", "% Road", "% Impervious", "% Dev-All", "% Dev-LMH", "% Dev-MH", "% Green",
                               "Cl Max", "75th Pct Cl", "Cl Med", "25th Pct Cl", "Cl Min", "Cl Mean", "Cl IQR", 
                               "WQ Rate Max", "WQ Rate Mean", "Invert Max", "Invert Mean", "CPUE Max", "CPUE Mean")
plot(CorrelationData)
cor1 = cor(CorrelationData, method="pearson")

plot(CorrelationData[,1:7])
corrplot.mixed(cor1[1:7, 1:7], lower.col="black", tl.col="black")
    # save as LULC_Correlations.jpeg

plot(PctRoad~PctDev_MH, data=Metrics)
lm1 = lm(PctRoad~PctDev_MH, data=Metrics)
abline(lm1)
anova(lm1)  #p=3.065e-6
summary(lm1)

plot(CorrelationData[,c(8:10, 13:14)])
corrplot.mixed(cor1[c(8:10, 13:14), c(8:10, 13:14)], lower.col="black", tl.col="black")
    # save as Cl_Correlations.jpeg


plot(ClMax~ClMed, data=Metrics)
lm2 = lm(ClMax~ClMed, data=Metrics)
abline(lm2)
anova(lm2)  #p=0.001966
summary(lm2)

#### CHLORIDE AS A FUNCTION OF LAND USE ####

## Test for Normality ##

hist(Metrics$ClMax)
hist(Metrics$ClMed)
shapiro.test(Metrics$ClMax)  #p=0.0001046 NOT NORMAL
shapiro.test(Metrics$ClMed)  #p=5.114e-5 NOT NORMAL
hist(log10(Metrics$ClMax))
hist(log10(Metrics$ClMed))
shapiro.test(log10(Metrics$ClMax))  #p=0.4762 NORMAL
shapiro.test(log10(Metrics$ClMed))  #p=0.08913 NORMAL

par(mfrow=c(2,2))
plot(log10(ClMax) ~ PctRoad * PctDev_MH, data=Metrics)
plot(log10(ClMed) ~ PctRoad * PctDev_MH, data=Metrics)

## MULTIPLE REGRESSION ##

lm1 = lm(log10(ClMax) ~ PctRoad * PctDev_MH + WatershedAreaKm2, data=Metrics)
summary(lm1)
plot(lm1)

lm1b = lm(log10(ClMax) ~ PctRoad, data=Metrics)
summary(lm1b)
plot(lm1b)

lm2 = lm(log10(ClMed) ~ PctRoad * PctDev_MH + WatershedAreaKm2, data=Metrics)
summary(lm2)
plot(lm2)

lm2a = lm(log10(ClMed) ~ PctRoad, data=Metrics)
lm2b = lm(log10(ClMed) ~ PctDev_MH, data=Metrics)

par(mfrow=c(1,1), mar=c(5,5,2,2))
#plot(log10(ClMax) ~ PctRoad, data=Metrics, ylab="log (Cl Max)", xlab="% Road", pch=16)
#abline(lm1b)
plot((ClMax) ~ PctRoad, data=Metrics, ylab="Cl Max", xlab="% Road", cex.lab=2, cex.axis=1.5, pch=16)
curve((10^2.28986)*(10^(0.12038*x)), add=T, lwd=3) 
abline(h=860, col="red", lty=2, lwd=3)

par(mfrow=c(2,1), oma = c(3,2.5,0,0), mar = c(4,4,1,1))
plot((ClMed) ~ PctRoad, data=Metrics, log="y", ylab="Cl Med", xlab="% Road", pch=16)
abline(lm2a)
plot((ClMed) ~ PctDev_MH, data=Metrics, log="y", ylab="Cl Med", xlab="% MH", pch=16)
abline(lm2b)

## Understanding the Interaction ##
library(interactions)
library(ggplot2)
interact_plot(lm2, pred=PctRoad, modx=PctDev_MH, modx.values = "plus-minus", interval = T, int.width=0.9, x.label="% Road", y.label="log (Cl Med)") + 
  theme_classic() +
  theme(axis.text=element_text(size=14), axis.title=element_text(size=18))
interact_plot(lm2, pred=PctRoad, modx=PctDev_MH, modx.values = "plus-minus", interval = T, int.width=0.9, x.label="% Road", y.label="log (Cl Med)", plot.points=T)
interact_plot(lm2, pred=PctRoad, modx=PctDev_MH, plot.points=T)

interact_plot(lm2, pred=PctDev_MH, modx=PctRoad, modx.values = "plus-minus", interval = T, int.width=0.9, x.label="% MH", y.label="log (Cl Med)")
interact_plot(lm2, pred=PctDev_MH, modx=PctRoad, plot.points=T)

#### INVERTEBRATES AS A FUNCTION OF CHLORIDE ####

library(quantreg)


## GENERATING FULL DATA TABLE BY INVERTEBRATE SAMPLE ##

yyy = data.frame(Invertebrates$SITENO, Invertebrates$TOTALINV, Invertebrates$CPUE, Invertebrates$WQRATING)
xx = data.frame(Metrics$SiteNo, Metrics$ClMax, Metrics$Cl75th, Metrics$ClMed)
SiteMetrics <- merge(yyy, xx, by.x=1, by.y=1, all.x=T)
xxx = data.frame(SiteMetrics$Metrics.ClMax, SiteMetrics$Metrics.Cl75th, SiteMetrics$Metrics.ClMed)


## QUANTILE REGRESSION ##

par(mfrow=c(2,2))
PDFPath = "C://Users/strea/Documents/My Publications/Post Dissertation 3-Cit Sci/Inverts_x_Chloride_Both.pdf"
pdf(file=PDFPath)
for(i in 1:ncol(yyy)){
  for(j in 1:ncol(xxx)){
    rq95 = rq(yyy[,i]~xxx[,j], tau=0.95)
    rq90 = rq(yyy[,i]~xxx[,j], tau=0.90)
    plot(yyy[,i]~xxx[,j], pch=19, ylab=names(yyy[i]), xlab=names(xxx[j]), na.rm=T)
    abline(lm(yyy[,i]~xxx[,j]), lwd=1)
    abline(rq95, col="blue", lwd=3, lty=3)
    abline(rq90, col="blue", lwd=2, lty=1)
    legend("topright", legend=c("95th Pct", "90th Pct", "Mean", "Chronic", "Acute"),
           col=c("blue", "blue", "black", "orange", "red"), lty=c(3, 1, 1, 1, 3), lwd=c(3, 2, 1, 1.5, 1.5))
    abline(v=230, col="orange", lwd=1.5)
    abline(v=860, col="red", lty=3, lwd=1.5)
  }
}
dev.off()


par(mfrow = c(3,2), oma = c(3,3,0,0), mar = c(3,3,1,1))
#Tot Invert
rq95 = rq(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMax, tau=0.95)
rq90 = rq(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMax, tau=0.90)
rlm = lm(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMax)
summary(rlm)  # Adj R2 = 0.1352; p = 1.3e-8 
plot(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMax, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
#legend("topright", legend=c("95th Pct", "90th Pct", "Mean", "Chronic", "Acute"),
#       col=c("blue", "blue", "black", "orange", "red"), lty=c(3, 1, 1, 1, 3), lwd=c(3, 2, 1, 1.5, 1.5))
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)

rq95 = rq(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMed, tau=0.95)
rq90 = rq(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMed, tau=0.90)
rlm = lm(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMed)
summary(rlm)  # Adj R2 = 0.0566; p = 2.3e-4 
plot(yyy$Invertebrates.TOTALINV~xxx$SiteMetrics.Metrics.ClMed, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)


#CPUE
rq95 = rq(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMax, tau=0.95)
rq90 = rq(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMax, tau=0.90)
rlm = lm(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMax)
summary(rlm)  # Adj R2 = 0.09588; p = 2.0e-6 
plot(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMax, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)

rq95 = rq(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMed, tau=0.95)
rq90 = rq(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMed, tau=0.90)
rlm = lm(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMed)
summary(rlm)  # Adj R2 = 0.1006; p = 1.1e-6 
plot(yyy$Invertebrates.CPUE~xxx$SiteMetrics.Metrics.ClMed, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)


#WQ Rating
rq95 = rq(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMax, tau=0.95)
rq90 = rq(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMax, tau=0.90)
rlm = lm(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMax)
summary(rlm)  # Adj R2 = 0.07044; p = 4.3e-5 
plot(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMax, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
rect(-150, 23, 8320, 35.1, col="grey88", border=NA)
points(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMax, pch=19)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)

rq95 = rq(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMed, tau=0.95)
rq90 = rq(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMed, tau=0.90)
rlm = lm(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMed)
summary(rlm)  # Adj R2 = 0.05949; p = 1.6e-4 
plot(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMed, pch=19, cex.axis=1.5)
#abline(rlm, lwd=1)
rect(20, 23, 910, 35.1, col="grey88", border=NA)
points(yyy$Invertebrates.WQRATING~xxx$SiteMetrics.Metrics.ClMed, pch=19)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)
abline(v=230, col="orange", lwd=2)
abline(v=860, col="red", lty=3, lwd=2.5)

mtext(text=" Maximum Chloride (mg/l)        Median Chloride (mg/l)", cex=1.5, side=1, line=0, outer=TRUE)
mtext(text="  WQ Rating          CPUE (animals/min)        # of Animals", cex=1.5, side=2, line=0, outer=TRUE)


#### INVERTEBRATES AS A FUNCTION OF URBAN DEVELOPMENT ####

## GENERATING FULL DATA TABLE BY INVERTEBRATE SAMPLE ##

xxb = data.frame(Metrics$SiteNo, Metrics$PctRoad, Metrics$PctDev_MH)
SiteMetricsB <- merge(yyy, xxb, by.x=1, by.y=1, all.x=T)
xxxb = data.frame(SiteMetricsB$Metrics.PctRoad, SiteMetricsB$Metrics.PctDev_MH)

## QUANTILE REGRESSION ##

par(mfrow = c(3,2), oma = c(3,3,0,0), mar = c(3,3,1,1))
#Tot Invert
rq95 = rq(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.95)
rq90 = rq(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.90)
rlm = lm(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctRoad)
summary(rlm)  # Adj R2 = 0.1221; p = 7.0e-8 
plot(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctRoad, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)

rq95 = rq(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.95)
rq90 = rq(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.90)
rlm = lm(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctDev_MH)
summary(rlm)  # Adj R2 = 0.04403; p = 1.1e-3 
plot(yyy$Invertebrates.TOTALINV~xxxb$SiteMetricsB.Metrics.PctDev_MH, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)


#CPUE
rq95 = rq(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.95)
rq90 = rq(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.90)
rlm = lm(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctRoad)
summary(rlm)  # Adj R2 = 0.05543; p = 2.8e-4 
plot(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctRoad, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)

rq95 = rq(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.95)
rq90 = rq(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.90)
rlm = lm(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctDev_MH)
summary(rlm)  # Adj R2 = 0.02165; p = 1.7e-2 
plot(yyy$Invertebrates.CPUE~xxxb$SiteMetricsB.Metrics.PctDev_MH, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)


#WQ Rating
rq95 = rq(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.95)
rq90 = rq(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctRoad, tau=0.90)
rlm = lm(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctRoad)
summary(rlm)  # Adj R2 = 0.06862; p = 5.36e-5 
plot(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctRoad, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
rect(0.89, 23, 10.45, 35.2, col="grey88", border=NA)
points(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctRoad, pch=19)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)

rq95 = rq(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.95)
rq90 = rq(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctDev_MH, tau=0.90)
rlm = lm(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctDev_MH)
summary(rlm)  # Adj R2 = 0.03344; p = 3.91e-3 
plot(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctDev_MH, pch=19, cex.axis=1.2)
#abline(rlm, lwd=1)
rect(-1, 23, 39.27, 35.2, col="grey88", border=NA)
points(yyy$Invertebrates.WQRATING~xxxb$SiteMetricsB.Metrics.PctDev_MH, pch=19)
abline(rq95, col="blue", lwd=3.5, lty=3)
abline(rq90, col="blue", lwd=2.5, lty=1)

mtext(text="    % Road                                % Med/High Density Dev", cex=1.5, side=1, line=0, outer=TRUE)
mtext(text="   WQ Rating               CPUE (animals/min)               # of Animals", cex=1.5, side=2, line=0, outer=TRUE)
