##############################################################
# R Script that imports Sed Yield MFA Results along with
# Physical and Flow Metrics to create correlation matrices
# among all metrics (Figure 2.11 in Dissertation
##############################################################

# Clear data
rm(list = ls())

# set working dir
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload")
#library(wq) # Sens Slope
library(MASS) # RLM

source("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/minor.ticks.axis.R")

# Import tables from various analyses. Should all have same # of rows
yield.metrics <- read.table("results_yield_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
flow.metrics <- read.table("results_flow_metrics.csv", stringsAsFactors = F, header = TRUE, sep = ",")
site.info <- read.table("site_summary.csv", stringsAsFactors = F, header = TRUE, sep = ",")
reg.results <- read.table("results_bed_regression.csv", stringsAsFactors = F, header = TRUE, sep = ",")
peaks <- read.table("results_peaks.csv", stringsAsFactors = F, header = TRUE, sep = ",")

# 
# yield.metrics$Qeff.spread <- round ((yield.metrics$Qeff.75 - yield.metrics$Qeff.25) 
#                                        / yield.metrics$Qeff.e,3)

## Changes name so that tables can be merged based on "site.no"
# colnames(reg.results)[1:2] <- c("site.no", "site.nm")

# Check to see if tables match up: Should all be TRUE
sum( c(
  sum(yield.metrics$site.no == flow.metrics$site.no) == nrow(yield.metrics),
  sum(yield.metrics$site.no == site.info$site.no) == nrow(yield.metrics),
  sum(yield.metrics$site.no == reg.results$site.no) == nrow(yield.metrics), 
  sum(yield.metrics$site.no == peaks$site.no) == nrow(yield.metrics) ) ) == 4

# Calc some new, normalized flow metrics
# NOTE: site.info$Q1.5 contains Q1.5 calc'd from inst. peak annual max series
# with some Q1.5's calc'd from daily avg flow data where peak flow data not available
# Qbf <- ifelse(!is.na(site.info$Qbf), round(site.info$Qbf, 3), round(site.info$Q1.5, 3) )
Q1.5 <- site.info$Q1.5
Q2 <- peaks$Q2 * 0.028316847
Q2[which(is.na(Q2))] <- flow.metrics$Q2[which(is.na(Q2))]
Qavg <- peaks$Qavg * 0.028316847
Qbf <- site.info$Qbf
Qbf.mean <- round( Qbf / flow.metrics$mean, 3)
Qeff.mean <- round(yield.metrics$Qeff.e / flow.metrics$mean, 3)
Qeff.Q1.5 <- round(yield.metrics$Qeff.e / Q1.5, 3)
Q1.5.mean <- round(Q1.5 / flow.metrics$mean, 3)
Q1.5.Qavg <- round(Q1.5 / Qavg, 3)
Qeff.Qbf <-  round(yield.metrics$Qeff.e / Qbf, 3)
da.km2 <- round(site.info$da.km2, 0)


# Calc hydraulic / sediment metrics
omega <- signif( sqrt(1.65*9.81*site.info$d50/1000), 3)
u.star <-  signif( sqrt(9.81 *site.info$Dbf * site.info$slope), 3 )
u.o <- signif( u.star / omega, 3 )
tau.star <- signif( site.info$Dbf * site.info$slope / (1.65 * site.info$d50/1000), 3)
omega.star <- signif( Qbf / site.info$Wbf * 9.81 * site.info$slope / 
              (9.81 * 1.65 * site.info$d50/1000)^1.5, 3)
w.d <- signif( site.info$Wbf / site.info$Dbf, 3)

hydr.metrics <- data.frame(site.no = site.info$site.no, site.info$site.nm, 
                           slope = site.info$slope, d50 = site.info$d50, d84 = site.info$d84,
                           Qbf = Qbf, tau.star = tau.star, w.d = w.d, da.km2 = da.km2)

# append a few combined flow-sed yield metrics to the yield.metrics table
yield.metrics0 <- data.frame(yield.metrics, Qbf = Qbf)
yield.metrics <- data.frame(yield.metrics, 
                            Qbf = Qbf, 
                            Qeff.mean = Qeff.mean, 
                            Q50.mean = round( yield.metrics$Q.50 / flow.metrics$mean, 2 ),
                            Qeff.Q1.5 = Qeff.Q1.5,
                            Qeff.Qbf = Qeff.Qbf,
                            Qeff.25.Q1.5 = round( yield.metrics$Qeff.25 / Q1.5, 2 ),
                            Qeff.75.Q1.5 = round( yield.metrics$Qeff.75 / Q1.5, 2 ), 
                            Q.25.Q1.5 = round( yield.metrics$Q.25 / Q1.5, 2 ),
                            Q.50.Q1.5 = round( yield.metrics$Q.50 / Q1.5, 2 ),
                            Q.75.Q1.5 = round( yield.metrics$Q.75 / Q1.5, 2 ),
                            Qeff.25.Qbf = round( yield.metrics$Qeff.25 / Qbf, 2 ), 
                            Qeff.75.Qbf = round( yield.metrics$Qeff.75 / Qbf, 2 ),
                            Q.25.Qbf = round( yield.metrics$Q.25 / Qbf, 2 ),
                            Q.50.Qbf = round( yield.metrics$Q.50 / Qbf, 2 ),
                            Q.75.Qbf = round( yield.metrics$Q.75 / Qbf, 2 ))
## First look for relationships between flow and hydraulic metrics AND a very select 
## amount of yield metrics: 
## beta, Qeff, Qeff.yield, Qeff.spread, Q.50, yield.spread, Qeff.RI, and Q50.RI

yield.sub <- yield.metrics[c("site.no", "beta", "Qeff.e", "Qeff.yield", "Qeff.spread",
                             "Qeff.RI", "Q.50", "yield.spread", "Q50.RI" )]
yield.sub.norm <- yield.metrics[c("site.no", "beta", "Qeff.mean", "Qeff.yield", "Qeff.spread",
                                  "Qeff.RI", "Q50.mean", "yield.spread", "Q50.RI" )]
flow.sub <- flow.metrics[c("site.no", "site.nm", "yrs","mean", "cv", "skewness", 
                           "spread", "flash.RB")]
flow.sub <- data.frame(flow.sub, Q1.5 = Q1.5, Q1.5.mean = Q1.5.mean, da.km2 = da.km2)

merge.flow <- merge(flow.sub, yield.sub, by = "site.no")
merge.flow.norm <- merge(flow.sub, yield.sub.norm, by = "site.no")

merge.hydr <- merge(hydr.metrics, yield.sub, by = "site.no")
merge.hydr.norm <- merge(hydr.metrics, yield.sub.norm, by = "site.no")

################### Qeff.RI and Skewness ###################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Figures")
pdf("RI v Skewness.pdf", width = 8.5, height = 5)
  par(mfrow = c(1,2), mar = c(5.1,0,0,0), pty = "s", oma = c(0.5,5.1,1,0.5))
  
  scatter.smooth(log10(flow.metrics$skewness),log10(yield.metrics$Qeff.RI), 
                 xlab = NA, 
                 ylab = NA,
                 cex.axis = 0.8,
                 axes = F,
                 ylim = c(log10(0.8), log10(25)))
  abline(coef(lm(log10(yield.metrics$Qeff.RI)~log10(flow.metrics$skewness))), lty = 2)
  minor.ticks.axis(2,9, mn = -3, mx = 2, cex = 0.8)
  minor.ticks.axis(1,9, mn = 0, mx = 2, cex = 0.8)
  box()
  text(x = log10(2.5), y = log10(20), labels = expression(italic(Q[eff])))
  mtext("Return Interval (yrs)", line = 2.5, outer = F, side = 2)

  scatter.smooth(log10(flow.metrics$skewness),log10(yield.metrics$Q50.RI), 
       xlab = NA, 
       ylab = NA,
       cex.axis = 0.8,
       axes = F,
       ylim = c(log10(0.8), log10(25)))
  text(x = log10(2.5), y = log10(20), labels = expression(italic(Q[s50])))
  abline(coef(lm(log10(yield.metrics$Q50.RI)~log10(flow.metrics$skewness))), lty = 2)
  minor.ticks.axis(1,9, mn = 0, mx = 2, cex = 0.8)
  box()
  
  mtext("Bed Load Sites", col = "seagreen", line = -4, outer = T, side = 3)
  mtext("Skewness", line = -7.5, outer = T, side = 1)

dev.off()

##### Calc and Plot Correlation Matrices in logspace ##############################

label.flow <- names(merge.flow[3:10])
label.hydr <- names(merge.hydr[3:10])
label.yield <- names(merge.flow[12:19])
label.yield.norm <- names(merge.flow.norm[12:19])

cor.flow.s <- matrix(data = 0, ncol = length(label.yield), nrow = length(label.flow))
rownames(cor.flow.s) <- label.flow
colnames(cor.flow.s) <- label.yield.norm

cor.hydr.s <- matrix(data = 0, ncol = length(label.yield), nrow = length(label.hydr))
rownames(cor.hydr.s) <- label.hydr
colnames(cor.hydr.s) <- label.yield.norm

cor.mix.s <- matrix(data = 0, ncol = length(label.hydr), nrow = length(label.flow))
rownames(cor.mix.s) <- label.flow
colnames(cor.mix.s) <- label.hydr

# clean up flow metrics to remove zero values or INF values
merge.flow[which(merge.flow == Inf, arr.ind = T)] <- NA
merge.flow[which(merge.flow == 0, arr.ind = T)] <- NA
merge.hydr[which(merge.hydr == Inf, arr.ind = T)] <- NA
merge.hydr[which(merge.hydr == 0, arr.ind = T)] <- NA
merge.flow.norm[which(merge.flow == Inf, arr.ind = T)] <- NA
merge.flow.norm[which(merge.flow == 0, arr.ind = T)] <- NA
merge.hydr.norm[which(merge.hydr == Inf, arr.ind = T)] <- NA
merge.hydr.norm[which(merge.hydr == 0, arr.ind = T)] <- NA

# create correlation matrices
for (i in 1:length(label.flow)){
  for (j in 1:length(label.yield)) {
    cor.flow.s[i,j] <- round( cor( log(merge.flow.norm[,2+i]), 
                                   log(merge.flow.norm[,length(label.flow)+3+j]),
                                   use = "complete.obs", method = "s"), 2)
  }
}

for (i in 1:length(label.hydr)){
  for (j in 1:length(label.yield)) {
    cor.hydr.s[i,j] <- round( cor( log(merge.hydr.norm[,2+i]), 
                                   log(merge.hydr.norm[,length(label.hydr) + 1 +j]),
                                   use = "complete.obs", method = "s"), 2)
  }
}

for (i in 1:length(label.flow)){
  for (j in 1:length(label.hydr)) {
    cor.mix.s[i,j] <- round( cor( log(merge.flow[,2+i]), log(merge.hydr[,2+j]),
                                  use = "complete.obs", method = "s"), 2)
  }
}

# create vectors for labels and plots

# Flow-Yield plotting position
plot.x.f <- matrix(data = 1:length(label.flow), ncol = length(label.yield), 
                   nrow = length(label.flow))
plot.y.f <- matrix(data = 1:length(label.yield), ncol = length(label.yield), 
                   nrow = length(label.flow), byrow = T)

# Hydr-Yield plotting position
plot.x.h <- matrix(data = 1:length(label.hydr), ncol = length(label.yield), 
                   nrow = length(label.hydr))
plot.y.h <- matrix(data = 1:length(label.yield), ncol = length(label.yield), 
                   nrow = length(label.hydr), byrow = T)

# Flow-Hydr plotting position
plot.x.m <- matrix(data = 1:length(label.flow), ncol = length(label.hydr), 
                   nrow = length(label.flow))
plot.y.m <- matrix(data = 1:length(label.hydr), ncol = length(label.hydr), 
                   nrow = length(label.flow), byrow = T)

### Plot abbreviated Pairs in log-space ###########################################
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/Figures")
pdf("Pairs_metrics.pdf", width = 11, height = 8.5)
par(cex = .5)
pairs(flow.sub[3:ncol(flow.sub)], log = "xy", cex = 0.5, cex.axis = 0.5, 
      cex.lab = 0.8, cex.sub = 0.5)
pairs(hydr.metrics[3:ncol(hydr.metrics)], log = "xy", cex = 0.5, cex.axis = 0.5, 
      cex.lab = 0.8, cex.sub = 0.5)
pairs(merge.flow[3:ncol(merge.flow)], log = "xy", cex = 0.5, cex.axis = 0.5, 
      cex.lab = 0.8, cex.sub = 0.5)
pairs(merge.hydr[3:ncol(merge.hydr)], log = "xy", cex = 0.5, 
      cex.axis = 0.5, cex.lab = 0.8, cex.sub = 0.5)
dev.off()
###################################################################################
########## PLOT CORRELATION MATRICES ############################### 
pdf("Correlation_Matrices_Norm1.5.pdf", width = 11, height = 6)
par(pty = "s", mfrow = c(1,3), mar= c(5.1, 5.1, 4.1, 2.1))

# Flow - Yield
plot(plot.x.f, plot.y.f,
     cex = abs(cor.flow.s)*7,
     pch = 21,
     bg = "gray30",
     col = "gray30",
     axes = F,
     xlim = c(.5, length(label.flow) + 0.5),
     ylim = c(.5, length(label.yield) + 0.5),
     xlab = "",
     ylab = "",
     main = "Flow : Yield")

points(plot.x.f[which(cor.flow.s > 0, arr.ind = T)], 
       plot.y.f[which(cor.flow.s > 0, arr.ind = T)],
       cex = abs(cor.flow.s[which(cor.flow.s > 0, arr.ind = T)])*7,
       pch = 21,
       bg = "white")

text(plot.x.f[which(cor.flow.s >= .40, arr.ind = T)], 
     plot.y.f[which(cor.flow.s >= .40, arr.ind = T)],
     labels = cor.flow.s[which(cor.flow.s >= .40, arr.ind = T)],
     cex = 0.5)
text(plot.x.f[which(cor.flow.s <= -.40, arr.ind = T)], 
     plot.y.f[which(cor.flow.s <= -.40, arr.ind = T)],
     labels = cor.flow.s[which(cor.flow.s <= -.40, arr.ind = T)],
     col = "white", cex = 0.5)

axis(1, labels = label.flow, at = 1:length(label.flow), cex.axis = .8, las = 2)
axis(2, labels = label.yield.norm, at = 1:length(label.yield), cex.axis = .8, las = 2)
box()

# Hydr - Yield
plot(plot.x.h, plot.y.h,
     cex = abs(cor.hydr.s)*7,
     pch = 21,
     bg = "gray30",
     col = "gray30",
     axes = F,
     xlim = c(.5, length(label.hydr) + 0.5),
     ylim = c(.5, length(label.yield) + 0.5),
     xlab = "",
     ylab = "",
     main = "Physical : Yield")

points(plot.x.h[which(cor.hydr.s > 0, arr.ind = T)], 
       plot.y.h[which(cor.hydr.s > 0, arr.ind = T)],
       cex = abs(cor.hydr.s[which(cor.hydr.s > 0, arr.ind = T)])*7,
       pch = 21,
       bg = "white")

text(plot.x.h[which(cor.hydr.s >= .40, arr.ind = T)], 
     plot.y.h[which(cor.hydr.s >= .40, arr.ind = T)],
     labels = cor.hydr.s[which(cor.hydr.s >= .40, arr.ind = T)],
     cex = 0.5)
text(plot.x.h[which(cor.hydr.s <= -.40, arr.ind = T)], 
     plot.y.h[which(cor.hydr.s <= -.40, arr.ind = T)],
     labels = cor.hydr.s[which(cor.hydr.s <= -.40, arr.ind = T)],
     col = "white", cex = 0.5)

axis(1, labels = label.hydr, at = 1:length(label.hydr), cex.axis = .8, las = 2)
axis(2, labels = label.yield.norm, at = 1:length(label.yield), cex.axis = .8, las = 2)
mtext("Spearman Log-Correlation Matrices \nBed Load Sites, Normalized",
      col = "seagreen", 
      line = 3, font = 4)
box()

# Flow - Hydr
plot(plot.x.m, plot.y.m,
     cex = abs(cor.mix.s)*5,
     pch = 21,
     bg = "gray30",
     col = "gray30",
     axes = F,
     xlim = c(.5, length(label.flow) + 0.5),
     ylim = c(.5, length(label.hydr) + 0.5),
     xlab = "",
     ylab = "",
     main = "Flow : Physical")

points(plot.x.m[which(cor.mix.s > 0, arr.ind = T)], 
       plot.y.m[which(cor.mix.s > 0, arr.ind = T)],
       cex = abs(cor.mix.s[which(cor.mix.s > 0, arr.ind = T)])*7,
       pch = 21,
       bg = "white")

text(plot.x.m[which(cor.mix.s >= .40, arr.ind = T)], 
     plot.y.m[which(cor.mix.s >= .40, arr.ind = T)],
     labels = cor.mix.s[which(cor.mix.s >= .40, arr.ind = T)],
     cex = 0.5)
text(plot.x.m[which(cor.mix.s <= -.40, arr.ind = T)], 
     plot.y.m[which(cor.mix.s <= -.40, arr.ind = T)],
     labels = cor.mix.s[which(cor.mix.s <= -.40, arr.ind = T)],
     col = "white", cex = 0.5)

axis(1, labels = label.flow, at = 1:length(label.flow), cex.axis = .8, las = 2)
axis(2, labels = label.hydr, at = 1:length(label.hydr), cex.axis = .8, las = 2)
box()

dev.off()
##########################################################################

############ CALC SLOPE OF RELATIONSHIPS - Log-Linear OLS #############

label.flow <- names(merge.flow[3:10])
label.hydr <- names(merge.hydr[3:10])
label.yield <- names(merge.flow[12:19])
label.yield.norm <- names(merge.flow.norm[12:19])

slope.flow <- matrix(data = 0, ncol = length(label.yield), nrow = length(label.flow))
rownames(slope.flow) <- label.flow
colnames(slope.flow) <- label.yield.norm
p.flow <- slope.flow
r2.flow <- slope.flow

slope.hydr <- matrix(data = 0, ncol = length(label.yield), nrow = length(label.hydr))
rownames(slope.hydr) <- label.hydr
colnames(slope.hydr) <- label.yield.norm
p.hydr <- slope.hydr
r2.hydr <- slope.hydr

slope.mix <- matrix(data = 0, ncol = length(label.hydr), nrow = length(label.flow))
rownames(slope.mix) <- label.flow
colnames(slope.mix) <- label.hydr
p.mix <- slope.mix
r2.mix <- slope.mix

# clean up flow metrics to remove zero values or INF values
merge.flow[which(merge.flow == Inf, arr.ind = T)] <- NA
merge.flow[which(merge.flow == 0, arr.ind = T)] <- NA
merge.hydr[which(merge.hydr == Inf, arr.ind = T)] <- NA
merge.hydr[which(merge.hydr == 0, arr.ind = T)] <- NA
merge.flow.norm[which(merge.flow.norm == Inf, arr.ind = T)] <- NA
merge.flow.norm[which(merge.flow.norm == 0, arr.ind = T)] <- NA
merge.hydr.norm[which(merge.hydr.norm == Inf, arr.ind = T)] <- NA
merge.hydr.norm[which(merge.hydr.norm == 0, arr.ind = T)] <- NA

# create correlation matrices
for (i in 1:length(label.flow)){
  for (j in 1:length(label.yield)) {
    slope.flow[i,j] <- round( coef( lm( log(merge.flow.norm[,length(label.flow)+3+j]) ~
                                         log(merge.flow.norm[,2+i]) ) )[2], 2)
    p.flow[i,j] <- round( summary( lm( log(merge.flow.norm[,length(label.flow)+3+j]) ~
                                          log(merge.flow.norm[,2+i]) ) )$coefficients[2,4], 3)
    r2.flow[i,j] <- round( summary( lm( log(merge.flow.norm[,length(label.flow)+3+j]) ~
                                         log(merge.flow.norm[,2+i]) ) )$r.squared, 2)
  }
}

for (i in 1:length(label.hydr)){
  for (j in 1:length(label.yield)) {
    slope.hydr[i,j] <- round( coef( lm( log(merge.hydr.norm[,length(label.hydr) + 1 +j]) ~
                                         log(merge.hydr.norm[,2+i]) ) )[2], 2)
    p.hydr[i,j] <- round( summary( lm( log(merge.hydr.norm[,length(label.hydr) + 1 +j]) ~
                                      log(merge.hydr.norm[,2+i]) ) )$coefficients[2,4], 3)
    r2.hydr[i,j] <- round( summary( lm( log(merge.hydr.norm[,length(label.hydr) + 1 +j]) ~
                                         log(merge.hydr.norm[,2+i]) ) )$r.squared, 2)
  }
}

for (i in 1:length(label.flow)){
  for (j in 1:length(label.hydr)) {
    slope.mix[i,j] <- round( coef( lm( log(merge.hydr[,2+j]) ~ log(merge.flow[,2+i]) ) )[2], 2)
    p.mix[i,j] <-  round( summary( lm( log(merge.hydr[,2+j]) ~ log(merge.flow[,2+i]) ) )$coefficients[2,4], 3)
    r2.mix[i,j] <-  round( summary( lm( log(merge.hydr[,2+j]) ~ log(merge.flow[,2+i]) ) )$r.squared, 2)
  }
}

# create vectors for labels and plots

# Flow-Yield plotting position
plot.x.f <- matrix(data = 1:length(label.flow), ncol = length(label.yield), 
                   nrow = length(label.flow))
plot.y.f <- matrix(data = 1:length(label.yield), ncol = length(label.yield), 
                   nrow = length(label.flow), byrow = T)

# Hydr-Yield plotting position
plot.x.h <- matrix(data = 1:length(label.hydr), ncol = length(label.yield), 
                   nrow = length(label.hydr))
plot.y.h <- matrix(data = 1:length(label.yield), ncol = length(label.yield), 
                   nrow = length(label.hydr), byrow = T)

# Flow-Hydr plotting position
plot.x.m <- matrix(data = 1:length(label.flow), ncol = length(label.hydr), 
                   nrow = length(label.flow))
plot.y.m <- matrix(data = 1:length(label.hydr), ncol = length(label.hydr), 
                   nrow = length(label.flow), byrow = T)

########## PLOT CORRELATION MATRICES ############################### 
setwd("C:/Users/jsholtes/Dropbox/Research/Qeff/Analysis/Bedload/Figures")
pdf("Slope_Matrices_Norm_RI_mean.pdf", width = 11, height = 6)
  par(pty = "s", mfrow = c(1,3), mar= c(5.1, 5.1, 3.1, 2.1), oma = c(1,1,1,1))
  caxis = 1
  ctxt = 0.7 
  
  # Flow - Yield

  cex.plot1 <- ifelse(abs(slope.flow)/mean(abs(slope.flow))*4 < 6,
                      abs(slope.flow)/mean(abs(slope.flow))*4,
                      8)
  cex.plot2 <- ifelse(abs(slope.flow[which(slope.flow > 0, arr.ind = T)])*4 / mean(abs(slope.flow)) < 6,
                      abs(slope.flow[which(slope.flow > 0, arr.ind = T)])*4 / mean(abs(slope.flow)),
                      8)
  plot(plot.x.f, plot.y.f,
       cex = cex.plot1,
       pch = 21,
       bg = "gray30",
       col = "black",
       axes = F,
       xlim = c(.5, length(label.flow) + 0.5),
       ylim = c(.5, length(label.yield) + 0.5),
       xlab = "",
       ylab = "",
       main = "Flow : Yield")
  
  points(plot.x.f[which(slope.flow > 0, arr.ind = T)], 
         plot.y.f[which(slope.flow > 0, arr.ind = T)],
         cex = cex.plot2,
         pch = 21,
         bg = "white")
  
  text(plot.x.f[which(p.flow < 0.05 & slope.flow > 0, arr.ind = T)], 
       plot.y.f[which(p.flow < 0.05 & slope.flow > 0, arr.ind = T)],
       labels = r2.flow[which(p.flow < 0.05 & slope.flow > 0, arr.ind = T)],
       cex = ctxt, font = 2)
  text(plot.x.f[which(p.flow < 0.05 & slope.flow < 0, arr.ind = T)], 
       plot.y.f[which(p.flow < 0.05 & slope.flow < 0, arr.ind = T)],
       labels = r2.flow[which(p.flow < 0.05 & slope.flow < 0, arr.ind = T)],
       cex = ctxt, font = 2, col = "palegreen3")
  
  axis(1, labels = label.flow, at = 1:length(label.flow), cex.axis = caxis, las = 2)
  axis(2, labels = label.yield.norm, at = 1:length(label.yield), cex.axis = caxis, las = 2)
  box()
  
  # Hydr - Yield
  cex.plot1 <- ifelse(abs(slope.hydr)/mean(abs(slope.hydr))*4 < 6,
                      abs(slope.hydr)/mean(abs(slope.hydr))*4,
                      8)
  cex.plot2 <- ifelse(abs(slope.hydr[which(slope.hydr > 0, arr.ind = T)])*4 / mean(abs(slope.hydr)) < 6,
                      abs(slope.hydr[which(slope.hydr > 0, arr.ind = T)])*4 / mean(abs(slope.hydr)),
                      8)
    
  plot(plot.x.h, plot.y.h,
       cex = cex.plot1,
       pch = 21,
       bg = "gray30",
       col = "black",
       axes = F,
       xlim = c(.5, length(label.hydr) + 0.5),
       ylim = c(.5, length(label.yield) + 0.5),
       xlab = "",
       ylab = "",
       main = "Physical : Yield")
  
  points(plot.x.f[which(slope.hydr > 0, arr.ind = T)], 
         plot.y.f[which(slope.hydr > 0, arr.ind = T)],
         cex = cex.plot2,
         pch = 21,
         bg = "white")
  
  text(plot.x.f[which(p.hydr < 0.05 & slope.hydr > 0, arr.ind = T)], 
       plot.y.f[which(p.hydr < 0.05 & slope.hydr > 0, arr.ind = T)],
       labels = r2.hydr[which(p.hydr < 0.05 & slope.hydr > 0, arr.ind = T)],
       cex = ctxt, font = 2)
  text(plot.x.f[which(p.hydr < 0.05 & slope.hydr < 0, arr.ind = T)], 
       plot.y.f[which(p.hydr < 0.05 & slope.hydr < 0, arr.ind = T)],
       labels = r2.hydr[which(p.hydr < 0.05 & slope.hydr < 0, arr.ind = T)],
       cex = ctxt, font = 2, col = "palegreen3")
  
  axis(1, labels = label.hydr, at = 1:length(label.hydr), cex.axis = caxis, las = 2)
  axis(2, labels = label.yield.norm, at = 1:length(label.yield), cex.axis = caxis, las = 2)
  mtext(bquote("Log-Linear OLS Model: Slope -"~ R^2 ~ "Plots"),
        col = "seagreen", 
        line = 4)
  mtext("Bed Load Sites, Normalized",
        col = "seagreen", 
        line = 3)
  box()
  
  # Flow - Hydr
  
  cex.plot1 <- ifelse(abs(slope.mix)/mean(abs(slope.mix))*4 < 6,
                      abs(slope.mix)/mean(abs(slope.mix))*4,
                      8)
  cex.plot2 <- ifelse(abs(slope.mix[which(slope.mix > 0, arr.ind = T)])*4 / mean(abs(slope.mix)) < 6,
                      abs(slope.mix[which(slope.mix > 0, arr.ind = T)])*4 / mean(abs(slope.mix)),
                      8)
  plot(plot.x.m, plot.y.m,
       cex = cex.plot1,
       pch = 21,
       bg = "gray30",
       col = "black",
       axes = F,
       xlim = c(.5, length(label.flow) + 0.5),
       ylim = c(.5, length(label.hydr) + 0.5),
       xlab = "",
       ylab = "",
       main = "Flow : Physical")
  
  points(plot.x.f[which(slope.mix > 0, arr.ind = T)], 
         plot.y.f[which(slope.mix > 0, arr.ind = T)],
         cex = cex.plot2,
         pch = 21,
         bg = "white")
  
  text(plot.x.f[which(p.mix < 0.05 & slope.mix > 0, arr.ind = T)], 
       plot.y.f[which(p.mix < 0.05 & slope.mix > 0, arr.ind = T)],
       labels = r2.hydr[which(p.mix < 0.05 & slope.mix > 0, arr.ind = T)],
       cex = ctxt, font = 2)
  text(plot.x.f[which(p.mix < 0.05 & slope.mix < 0, arr.ind = T)], 
       plot.y.f[which(p.mix < 0.05 & slope.mix < 0, arr.ind = T)],
       labels = r2.mix[which(p.mix < 0.05 & slope.mix < 0, arr.ind = T)],
       cex = ctxt, font = 2, col = "palegreen3")
  
  axis(1, labels = label.flow, at = 1:length(label.flow), cex.axis = caxis, las = 2)
  axis(2, labels = label.hydr, at = 1:length(label.hydr), cex.axis = caxis, las = 2)
  box()

dev.off()




