# ----------------- # 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)