library(lme4)
library(emmeans)
library(lsmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)

# Set the working directory.
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/HondoManuscript_AquaticConservation_Code-Data/Analysis7")
#setwd("C:/Users/mbutler/Desktop/Rio Hondo/Hondo Spring-system Compared to Reference systems (Analysis 3, 7)")
#setwd("C:/A_USFWS/Bitter Lake NWR/Bitter Lake Invert Draft Protocol Code/")


# Analysis 7 Water Parameters for Summer 2020-2021 accross springsystems.
water.data=read.csv ("Analysis7_WaterParameters_AcrossSpringSystems_Summers2020-2021_19Nov21.csv", header=TRUE)
names(water.data)

# A Comparison of water parameter across systems using using MANOVA.
# MANOVA function is below.
Y=cbind(water.data$Water_Temp_C,water.data$SAL_ppt,water.data$pH,water.data$Water_depth_cm,water.data$DO_mg_L)
modelManova=manova(Y~Site, data=water.data)
# view the results.
summary(modelManova, 'Wilks')

# Water Temperature by Springsystem
ModelTemp1=lm(Water_Temp_C~Site,data=water.data)
summary(ModelTemp1)

# Multiple comparison test for Water Temp
emmeans(ModelTemp1,list(pairwise~Site),adjust="tukey",level=0.95)
TUKEYTemp=emmeans(ModelTemp1,list(pairwise~Site),adjust="tukey",level=0.95)
plot(TUKEYTemp, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: Celsius", ylab="Spring-system", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# Get the info for plotting.
TUKEYTemp.plot=plot(TUKEYTemp,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# Dissolved Oxygen by Springsystem
ModelDO1=lm(DO_mg_L~Site,data=water.data)
summary(ModelDO1)

# Multiple comparison test for Dissolved Oxygen
emmeans(ModelDO1,list(pairwise~Site),adjust="tukey",level=0.95)
TUKEYDO=emmeans(ModelDO1,list(pairwise~Site),adjust="tukey",level=0.95)
plot(TUKEYDO, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: DO", ylab="Spring-system", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# Get the info for plotting.
TUKEYDO.plot=plot(TUKEYDO,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# pH by Springsystem
ModelpH1=lm(pH~Site,data=water.data)
summary(ModelpH1)

# Multiple comparison test for Dissolved Oxygen
emmeans(ModelpH1,list(pairwise~Site),adjust="tukey",level=0.95)
TUKEYpH=emmeans(ModelpH1,list(pairwise~Site),adjust="tukey",level=0.95)
plot(TUKEYpH, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: pH", ylab="Spring-system", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# Get the info for plotting.
TUKEYpH.plot=plot(TUKEYpH,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# Salinity by Springsystem
ModelSalt1=lm(SAL_ppt~Site,data=water.data)
summary(ModelSalt1)

# Multiple comparison test for Dissolved Oxygen
emmeans(ModelSalt1,list(pairwise~Site),adjust="tukey",level=0.95)
TUKEYSalt=emmeans(ModelSalt1,list(pairwise~Site),adjust="tukey",level=0.95)
plot(TUKEYSalt, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: Salt", ylab="Spring-system", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# Get the info for plotting.
TUKEYSalt.plot=plot(TUKEYSalt,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# Water depth by Springsystem
ModelDepth1=lm(Water_depth_cm~Site,data=water.data)
summary(ModelDepth1)

# Multiple comparison test for Dissolved Oxygen
emmeans(ModelDepth1,list(pairwise~Site),adjust="tukey",level=0.95)
TUKEYDepth=emmeans(ModelDepth1,list(pairwise~Site),adjust="tukey",level=0.95)
plot(TUKEYDepth, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: Depth", ylab="Spring-system", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# Get the info for plotting.
TUKEYDepth.plot=plot(TUKEYDepth,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# Save the plot - this if manuscript Figure 8
jpeg("Figure8_Water_systems_summer.jpg",width=6.5,height=9,units="in",res=675)
layout(cbind(c(1,2,3),c(4,5,6)),heights=c(2,2,2),widths=c(2,2,2)); layout.show(6)
par(cex.lab=1.4)
par(font.lab=2)
par(cex.main=1.6)
par(font.main=2)
par(cex.axis=1.2)
par(lwd=2)
par(mar=c(7,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYTemp.plot$data$the.emmean,axes=FALSE,pch="",ylab="Temperature (°C)",ylim=c(15,35),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(14.3,35.6,35.6,14.3),col=adjustcolor("gray",alpha.f=0.31),border=NA)
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(15,35),col="black",las=1,yaxp=c(15,35,4))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSprings","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYTemp.plot$data$lower.CL[i],TUKEYTemp.plot$data$upper.CL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYTemp.plot$data$the.emmean,pch=19)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYTemp.plot$data$the.emmean[i],TUKEYTemp.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
  lines(rep(i,2),c(TUKEYTemp.plot$data$the.emmean[i],TUKEYTemp.plot$data$rcmpl[i]),lty=1)
}

plot(1:5,TUKEYDO.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~bold("Dissolved oxygen (mg L"^-1~")")),ylim=c(0,8),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-0.3,10.3,10.3,-0.3),col=adjustcolor("gray",alpha.f=0.31),border=NA)
axis(2,ylim=c(0,8),col="black",las=1,yaxp=c(0,8,4))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSprings","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDO.plot$data$lower.CL[i],TUKEYDO.plot$data$upper.CL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYDO.plot$data$the.emmean,pch=19)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDO.plot$data$the.emmean[i],TUKEYDO.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDO.plot$data$the.emmean[i],TUKEYDO.plot$data$rcmpl[i]),lty=1)
}

plot(1:5,TUKEYpH.plot$data$the.emmean,axes=FALSE,pch="",ylab="pH",ylim=c(6.5,8.5),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(6.2,8.7,8.7,6.2),col=adjustcolor("gray",alpha.f=0.31),border=NA)
axis(2,ylim=c(6.5,8.5),col="black",las=1,yaxp=c(6.5,8.5,4))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSprings","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYpH.plot$data$lower.CL[i],TUKEYpH.plot$data$upper.CL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYpH.plot$data$the.emmean,pch=19)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYpH.plot$data$the.emmean[i],TUKEYpH.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
  lines(rep(i,2),c(TUKEYpH.plot$data$the.emmean[i],TUKEYpH.plot$data$rcmpl[i]),lty=1)
}

plot(1:5,TUKEYSalt.plot$data$the.emmean,axes=FALSE,pch="",ylab="Salinity (ppt)",ylim=c(0,10),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-0.3,10.4,10.4,-0.3),col=adjustcolor("gray",alpha.f=0.31),border=NA)
axis(2,ylim=c(0,10),col="black",las=1,yaxp=c(0,10,5))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSprings","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYSalt.plot$data$lower.CL[i],TUKEYSalt.plot$data$upper.CL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSalt.plot$data$the.emmean,pch=19)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYSalt.plot$data$the.emmean[i],TUKEYSalt.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
  lines(rep(i,2),c(TUKEYSalt.plot$data$the.emmean[i],TUKEYSalt.plot$data$rcmpl[i]),lty=1)
}

plot(1:5,TUKEYDepth.plot$data$the.emmean,axes=FALSE,pch="",ylab="Water depth (cm)",ylim=c(0,35),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-1,37,37,-1),col=adjustcolor("gray",alpha.f=0.31),border=NA)
axis(2,ylim=c(0,35),col="black",las=1,yaxp=c(0,35,7))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSprings","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDepth.plot$data$lower.CL[i],TUKEYDepth.plot$data$upper.CL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYDepth.plot$data$the.emmean,pch=19)
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDepth.plot$data$the.emmean[i],TUKEYDepth.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
  lines(rep(i,2),c(TUKEYDepth.plot$data$the.emmean[i],TUKEYDepth.plot$data$rcmpl[i]),lty=1)
}
dev.off()

