# 29 Mar 2022 - Analysis Rio Hondo springvents Inverts and Water
library(lme4)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(graphics)

# Set the working directory.
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/HondoManuscript_AquaticConservation_Code-Data/Analysis1and5/Analyses1and5_27Feb22/Analyses_RioHondo_Springvents")

#setwd("C:/Users/mbutler/Desktop/Rio Hondo/Hondo Spring-system Through Time (Analysis 1, 5)")
#setwd("C:/A_USFWS/Bitter Lake NWR/Bitter Lake Invert Draft Protocol Code/")

# Analysis 1 Noels amphipods in Rio Hondo Springvents
# Analysis 1 dataset.
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis1_InvertDensity_RioHondo_ThroughTime_Streamvents_19Nov21.csv",header=T,sep=",")

# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
data$Trap_size_sq_m[data$Trap_size_sq_m==0.006]=(10.16*3.81)

(1/((5.08*3.81)/10000))/(1/0.003)
(1/((10.16*3.81)/10000))/(1/0.006)
# 1.550003

names(data)

#Note: Subsample Location A=Bank/Vent, B=Quarter/midrun, C=Mid/endrun, Relabled as ABC so output graphs pertaining to subsample would be in order

# Noels Amphipods
# glmer without addressing scale of variables, you get "fail to converge error"

# Linear, quadratic and cubic models for Amphipods
summary(glmer(Noels_Alive~I(Year-2014)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Noels_Alive~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Noels_Alive~I((Year-2014)/10)+I(((Year-2014)/10)^2)+I(((Year-2014)/10)^3)+Season+Subsample+(1|Sample),data=data,family=poisson(link="log"),offset=log(Trap_size_sq_m)))

# Make list of candidate models
Cand.models.amp=list()
Cand.models.amp[["Year"]]=glmer(Noels_Alive~I((Year-2014)/10)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.amp[["Year^2"]]=glmer(Noels_Alive~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.amp[["Year^3"]]=glmer(Noels_Alive~I((Year-2014)/10)+I(((Year-2014)/10)^2)+I(((Year-2014)/10)^3)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))

# Make AICc table.
AIC.Table.amp=model.sel(Cand.models.amp,rank="AICc")
AIC.Table.amp

# Save the cubic model.
model1=glmer(Noels_Alive~I((Year-2014)/10)+I(((Year-2014)/10)^2)+I(((Year-2014)/10)^3)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
summary(model1)

# Make data for graph of prediction.
GRAPH.DATA=data.frame(Subsample=rep("A",16),Year=rep(seq(2014,2021,1),2),Season=c(rep("Winter",8),rep("Summer",8)),Trap_size_sq_m=rep(1,16),Sample="")
GRAPH.DATA$Pred=predict(model1,newdata=GRAPH.DATA,type="response",allow.new.levels=TRUE)

# Bootstrap the predictions in order to get CIs. NO R IS NOT BROKEN, this takes a long, long while depending on computer speed
model1.boot=bootMer(model1,function(x) predict(x,newdata=GRAPH.DATA,type="response",allow.new.levels=TRUE),nsim=1000)

# Get 95% CIs.
model1.boot.CI=data.frame(t(apply(model1.boot$t,2,quantile,c(0.025,0.975))))
names(model1.boot.CI)=c("LCL","UCL")

# Make and save a high resolution figure.
jpeg("Noel_Trend_RHvent.jpg",width=6.5,height=4.5,units="in",res=675)
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(5,5,2,2))
plot(2014:2021,GRAPH.DATA$Pred[9:16],ylim=c(0,0.6),type="l",lty=2,ylab=expression(~italic("Gammarus desperatus")~plain(~"no. cm"^-2~"")),xlab="Year",las=1)
points(2014:2021,GRAPH.DATA$Pred[9:16],pch=19)
lines(2014:2021,GRAPH.DATA$Pred[1:8],lty=3)
points(2014:2021,GRAPH.DATA$Pred[1:8],pch=19)
points(2014:2021,model1.boot.CI$LCL[1:8],pch=151)
points(2014:2021,model1.boot.CI$UCL[1:8],pch=151)
points(2014:2021,model1.boot.CI$LCL[9:16],pch=151)
points(2014:2021,model1.boot.CI$UCL[9:16],pch=151)
for (i in 1:8) {
  lines(rep(GRAPH.DATA$Year[i],2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
  lines(rep(GRAPH.DATA$Year[i+8],2),c(model1.boot.CI$LCL[i+8],model1.boot.CI$UCL[i+8]))
}
legend("topright",legend=c("Summer","Winter"),lty=c(2,3))
dev.off()

# Multi comparisons Noel's amphipods
emmeans(model1,pairwise~Season,offset=1,type="response",adjust="tukey",level=0.95)
emmeans(model1,pairwise~Subsample,offset=1,type="response",adjust="tukey",level=0.95)

# plot of Multiple comparison test: Season 
TUKEYAmphipods=emmeans(model1,list(pairwise~Season),offset=1,type="response",adjust="tukey",level=0.95)
plot(TUKEYAmphipods, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: Noel's Amphipod", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))


# plot of Multiple comparison test: Subsample 
# Get the data from the emmeans.
TUKEYAmphipods1=emmeans(model1,list(pairwise~Subsample),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYAmphipods1.plot=plot(TUKEYAmphipods1,comparison=TRUE,plotIt=FALSE,adjust="tukey")

# Save the plot
jpeg("Noel_ABC_RHvent.jpg",width=6.5,height=4.5,units="in",res=675)
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(5,6,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
#plot(1:3,TUKEYAmphipods1.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~"Gammarus"~desperatus~"(no. cm"^-2~")\n"),ylim=c(0,0.2),xlim=c(0.5,3.5),xlab="Springrun Location")
#box()
plot(1:3,TUKEYAmphipods1.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~italic("Gammarus desperatus")~plain(~"no. cm"^-2~"")),ylim=c(0,0.2),xlim=c(0.5,3.5),xlab="Springrun Location")
box()
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(0,0.2),col="black",las=1,yaxp=c(0,0.2,10))
axis(1,at=1:3,labels=c("Vent","Mid-run","End-run"),col="black",las=1,cex=2)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$asymp.LCL[i],TUKEYAmphipods1.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:3,TUKEYAmphipods1.plot$data$the.emmean,pch=19)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$the.emmean[i],TUKEYAmphipods1.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$the.emmean[i],TUKEYAmphipods1.plot$data$rcmpl[i]),lty=1)
}
dev.off()

# Make and save a high resolution figure.
jpeg("Figure3_NoelsAmphipods_RHvent.jpg",width=6.5,height=9,units="in",res=675)
layout(rbind(1,2),heights=c(2,2),widths=c(2,2)); layout.show(2)
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(5,6,1.6,2))
plot(2014:2021,GRAPH.DATA$Pred[9:16],ylim=c(0,0.6),type="l",lty=2,ylab=expression(~bolditalic("Gammarus desperatus")~bold(~"(no. cm"^-2~")")),xlab="Year",las=1)
points(2014:2021,GRAPH.DATA$Pred[9:16],pch=19)
lines(2014:2021,GRAPH.DATA$Pred[1:8],lty=3)
points(2014:2021,GRAPH.DATA$Pred[1:8],pch=19)
points(2014:2021,model1.boot.CI$LCL[1:8],pch=151)
points(2014:2021,model1.boot.CI$UCL[1:8],pch=151)
points(2014:2021,model1.boot.CI$LCL[9:16],pch=151)
points(2014:2021,model1.boot.CI$UCL[9:16],pch=151)
for (i in 1:8) {
  lines(rep(GRAPH.DATA$Year[i],2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
  lines(rep(GRAPH.DATA$Year[i+8],2),c(model1.boot.CI$LCL[i+8],model1.boot.CI$UCL[i+8]))
}
legend("topright",legend=c("Summer","Winter"),lty=c(2,3))
plot(1:3,TUKEYAmphipods1.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~bolditalic("Gammarus desperatus")~bold(~"(no. cm"^-2~")")),ylim=c(0,0.2),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=0.58,las=1,side=3,line=+21.0, cex=1.2,font=1)
mtext("(b)",col="black",at=0.58,las=1,side=3,line=-1.5, cex=1.2,font=1)
box()
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(0,0.2),col="black",las=1,yaxp=c(0,0.2,5))
axis(1,at=1:3,labels=c("Vent","Mid-run","End-run"),col="black",las=1,cex=2)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$asymp.LCL[i],TUKEYAmphipods1.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:3,TUKEYAmphipods1.plot$data$the.emmean,pch=19)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$the.emmean[i],TUKEYAmphipods1.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:3){
  lines(rep(i,2),c(TUKEYAmphipods1.plot$data$the.emmean[i],TUKEYAmphipods1.plot$data$rcmpl[i]),lty=1)
}
dev.off()

##Start Springsnails here - Analysis 1 springsnails in Rio Hondo springvents

# Linear, Quadratic and Cubic models springsnails (cubic doesn't converge, so scaled for AICc candidate set below)
summary(glmer(Springsnails~I(Year-2014)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+I(((Year-2014)/10)^3)+Season+Subsample+(1|Sample),data=data,family=poisson(link="log"),offset=log(Trap_size_sq_m)))

#AICc table for springsnails
Cand.models.ss=list()
Cand.models.ss[["Year"]]=glmer(Springsnails~I((Year-2014)/10)+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["Year+Season"]]=glmer(Springsnails~I((Year-2014)/10)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["Year^2"]]=glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["Year^2+Season"]]=glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["log(Year)"]]=glmer(Springsnails~I(log(Year-2013))+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["log(Year)+Season"]]=glmer(Springsnails~I(log(Year-2013))+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
Cand.models.ss[["Year^3"]]=glmer(Springsnails~scale(Year)+I(scale(Year)^2)+I(scale(Year)^3)+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
#Cand.models.ss[["Year^3+Season"]]=glmer(Springsnails~scale(Year)+I(scale(Year)^2)+I(scale(Year)^3)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))

# Make AICc table.
AIC.Table.ss=model.sel(Cand.models.ss,rank="AICc")
AIC.Table.ss

# Save the best model for table values (Appendix 2a Springsnails).
modelbest=glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
summary(modelbest)

# Save the quadratic model with season for plotting.
model2=glmer(Springsnails~I((Year-2014)/10)+I(((Year-2014)/10)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log"))
summary(model2)

# Make data for graph of prediction - quadratic with season and subsample.
GRAPH.DATA.ss=data.frame(Subsample=rep("A",16),Year=rep(seq(2014,2021,1),2),Season=c(rep("Winter",8),rep("Summer",8)),Trap_size_sq_m=rep(1,16),Sample="")
GRAPH.DATA.ss$Pred=predict(model2,newdata=GRAPH.DATA.ss,type="response",allow.new.levels=TRUE)

# Bootstrap the predictions in order to get CIs, this takes a long time depending on computer speed. 
model2.boot=bootMer(model2,function(x) predict(x,newdata=GRAPH.DATA.ss,type="response",allow.new.levels=TRUE),nsim=1000)

# Get 95% CIs.
model2.boot.CI=data.frame(t(apply(model2.boot$t[-63,],2,quantile,c(0.025,0.975))))
names(model2.boot.CI)=c("LCL","UCL")

# Make and save a high resolution figure.
jpeg("Snail_Trend_RHvent.jpg",width=6.5,height=4.5,units="in",res=600)
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(4.5,5,2,2))
plot(2014:2021,GRAPH.DATA.ss$Pred[9:16],ylim=c(0,0.2),type="l",lty=2,ylab=expression(~"Springsnails (no. cm"^-2~")"),xlab="Year",las=1)
points(2014:2021,GRAPH.DATA.ss$Pred[9:16],pch=19)
lines(2014:2021,GRAPH.DATA.ss$Pred[1:8],lty=3)
points(2014:2021,GRAPH.DATA.ss$Pred[1:8],pch=19)
points(2014:2021,model2.boot.CI$LCL[1:8],pch=151)
points(2014:2021,model2.boot.CI$UCL[1:8],pch=151)
points(2014:2021,model2.boot.CI$LCL[9:16],pch=151)
points(2014:2021,model2.boot.CI$UCL[9:16],pch=151)
for (i in 1:8) {
  lines(rep(GRAPH.DATA.ss$Year[i],2),c(model2.boot.CI$LCL[i],model2.boot.CI$UCL[i]))
  lines(rep(GRAPH.DATA.ss$Year[i+8],2),c(model2.boot.CI$LCL[i+8],model2.boot.CI$UCL[i+8]))
}
legend("topright",legend=c("Summer","Winter"),lty=c(2,3))
dev.off()

# Multi comparisons Springsnails
emmeans(model2,pairwise~Season,offset=1,type="response",adjust="tukey",level=0.95)
emmeans(model2,pairwise~Subsample,offset=1,type="response",adjust="tukey",level=0.95)

# plot of Multiple comparison test: Season 
TUKEYSpringsnails=emmeans(model2,offset=1,list(pairwise~Season),type="response",adjust="tukey",level=0.95)
plot(TUKEYSpringsnails, comparison = TRUE, adjust = "tukey", xlab="Estimated Marginal Means: Springsnails", horizontal = FALSE, colors = "black")+theme_bw()+theme(tex=element_text(face="bold",size = 15,color="black"))

# plot of Multiple comparison test: Subsample 
# Get the data from the emmeans.
TUKEYSpringsnails1=emmeans(model2,offset=1,list(pairwise~Subsample),type="response",adjust="tukey",level=0.95)

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

# Save the plot
jpeg("SS_ABC_RHvent.jpg",width=6.5,height=4.5,units="in",res=600)
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(5,6,2,2))
par(mex=1.6)
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:3,TUKEYSpringsnails1.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~"Springsnails (no. cm"^-2~")"),ylim=c(0,0.0016),xlim=c(0.5,3.5),xlab="Springrun Location")
mtext("a",col="black",at=1,las=1,side=1,line=0,cex=1.2,font=1)
box()
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(0,0.0016),las=1,col="black",yaxp=c(0,0.0016,4))
axis(1,at=1:3,labels=c("Vent","Mid-run","End-run"),col="black",las=1,cex=2)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$asymp.LCL[i],TUKEYSpringsnails1.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:3,TUKEYSpringsnails1.plot$data$the.emmean,pch=19)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$the.emmean[i],TUKEYSpringsnails1.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$the.emmean[i],TUKEYSpringsnails1.plot$data$rcmpl[i]),lty=1)
}
dev.off()

jpeg("Figure4_Springsnails_RHvents.jpg",width=6.5,height=9,units="in",res=675)
layout(rbind(1,2),heights=c(2,2),widths=c(2,2)); layout.show(2)
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(5,6,0.1,2))
par(mex=1.6)
plot(2014:2021,GRAPH.DATA.ss$Pred[9:16],ylim=c(0,0.2),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2014:2021,GRAPH.DATA.ss$Pred[9:16],pch=19)
lines(2014:2021,GRAPH.DATA.ss$Pred[1:8],lty=3)
points(2014:2021,GRAPH.DATA.ss$Pred[1:8],pch=19)
points(2014:2021,model2.boot.CI$LCL[1:8],pch=151)
points(2014:2021,model2.boot.CI$UCL[1:8],pch=151)
points(2014:2021,model2.boot.CI$LCL[9:16],pch=151)
points(2014:2021,model2.boot.CI$UCL[9:16],pch=151)
for (i in 1:8) {
  lines(rep(GRAPH.DATA.ss$Year[i],2),c(model2.boot.CI$LCL[i],model2.boot.CI$UCL[i]))
  lines(rep(GRAPH.DATA.ss$Year[i+8],2),c(model2.boot.CI$LCL[i+8],model2.boot.CI$UCL[i+8]))
}
legend("topright",legend=c("Summer","Winter"),lty=c(2,3))
plot(1:3,TUKEYSpringsnails1.plot$data$the.emmean,axes=FALSE,pch="",ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,0.0016),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=0.6,las=1,side=3,line=+13.2, cex=1.2,font=1)
mtext("(b)",col="black",at=0.6,las=1,side=3,line=-0.9, cex=1.2,font=1)
box()
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(0,0.0016),las=1,col="black",yaxp=c(0,0.0016,4))
axis(1,at=1:3,labels=c("Vent","Mid-run","End-run"),col="black",las=1,cex=2)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$asymp.LCL[i],TUKEYSpringsnails1.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:3,TUKEYSpringsnails1.plot$data$the.emmean,pch=19)
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$the.emmean[i],TUKEYSpringsnails1.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:3){
  lines(rep(i,2),c(TUKEYSpringsnails1.plot$data$the.emmean[i],TUKEYSpringsnails1.plot$data$rcmpl[i]),lty=1)
}
dev.off()

#Analysis 5 Water parameters in springvents starts here
# Designate data set.
water.data=read.csv ("Analysis5_WaterParameters_ThroughTime_RioHondo_Springvents_4Jun21.csv", header=TRUE)
names(water.data)

# A Comparison of water parameter across springvents using using MANOVA.
# MANOVA function is below.
modelManova=manova(cbind(Water_Temp_C,DO_mg_L,SAL_ppt,pH,Water_depth_cm)~I(Year)+I((Year)^2)+Season,data=water.data)
summary(modelManova,'Wilks')

modelManova2=manova(cbind(Water_Temp_C,DO_mg_L,SAL_ppt,pH,Water_depth_cm)~I(Year-2014)+Season,data=water.data)
summary(modelManova2,'Wilks')

#view the results.
summary(modelManova,'Wilks')

# view the results with response variable by response variable with ANOVA because all were significant.
# summary.aov(modelManova)

#Investigate change in WaterTemp parameter Hondo Springruns using lm
ModelSeasonHondoH20Temp=lm(Water_Temp_C~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelSeasonHondoH20Temp)

ModelSeasonHondoH20Temp2=lm(Water_Temp_C~I(Year)+Season,data=water.data)
summary(ModelSeasonHondoH20Temp2)

#AICc comparison of linear and singular
Cand.models.Temp=list()
Cand.models.Temp[["Year+Season"]]=lm(Water_Temp_C~I(Year)+Season,data=water.data)
Cand.models.Temp[["Year^2+Season"]]=lm(Water_Temp_C~I(Year)+I((Year)^2)+Season,data=water.data)
# Make AICc table.
AIC.Table.Temp=model.sel(Cand.models.Temp,rank="AICc")
AIC.Table.Temp

# Make data.frame with predictions and CIs.
graph.Water_Temp_C=data.frame(predict(Cand.models.Temp[["Year^2+Season"]],newdata=data.frame(Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8))),interval="confidence"),Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8)))

#Investigate change of in Dissolved Oxygen in Rio Hondo Springruns using lm
ModelSeasonHondoDO=lm(DO_mg_L~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelSeasonHondoDO)

ModelSeasonHondoDO2=lm(DO_mg_L~I(Year)+Season,data=water.data)
summary(ModelSeasonHondoDO2)

#AICc comparison of linear and singular
Cand.models.DO=list()
Cand.models.DO[["Year+Season"]]=lm(DO_mg_L~I(Year)+Season,data=water.data)
Cand.models.DO[["Year^2+Season"]]=lm(DO_mg_L~I(Year)+I((Year)^2)+Season,data=water.data)
# Make AICc table.
AIC.Table.DO=model.sel(Cand.models.DO,rank="AICc")
AIC.Table.DO

# Make data.frame with predictions and CIs.
graph.DO=data.frame(predict(Cand.models.DO[["Year^2+Season"]],newdata=data.frame(Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8))),interval="confidence"),Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8)))

#Investigate change in pH parameter in Hondo Springvents using lm
ModelSeasonHondopH=lm(pH~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelSeasonHondopH)

ModelSeasonHondopH2=lm(pH~I(Year)+Season,data=water.data)
summary(ModelSeasonHondopH2)

#AICc comparison of linear and singular
Cand.models.pH=list()
Cand.models.pH[["Year+Season"]]=lm(pH~I(Year)+Season,data=water.data)
Cand.models.pH[["Year^2+Season"]]=lm(pH~I(Year)+I((Year)^2)+Season,data=water.data)
# Make AICc table.
AIC.Table.pH=model.sel(Cand.models.pH,rank="AICc")
AIC.Table.pH

# Make data.frame with predictions and CIs.
graph.pH=data.frame(predict(Cand.models.pH[["Year^2+Season"]],newdata=data.frame(Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8))),interval="confidence"),Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8)))

#Investigate change in Salinity parameter Hondo Springvents using lm
ModelSeasonHondoSalt=lm(SAL_ppt~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelSeasonHondoSalt)

ModelSeasonHondoSalt2=lm(SAL_ppt~I(Year)+Season,data=water.data)
summary(ModelSeasonHondoSalt2)

#AICc comparison of linear and singular
Cand.models.salt=list()
Cand.models.salt[["Year+Season"]]=lm(SAL_ppt~I(Year)+Season,data=water.data)
Cand.models.salt[["Year^2+Season"]]=lm(SAL_ppt~I(Year)+I((Year)^2)+Season,data=water.data)
# Make AICc table.
AIC.Table.salt=model.sel(Cand.models.salt,rank="AICc")
AIC.Table.salt

# Make data.frame with predictions and CIs.
graph.salt=data.frame(predict(Cand.models.salt[["Year^2+Season"]],newdata=data.frame(Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8))),interval="confidence"),Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8)))

#Investigate change in Water Depth parameter Hondo Springvents using lm
ModelH20cm=lm(Water_depth_cm~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelH20cm)

ModelH20cm2=lm(Water_depth_cm~I(Year)+Season,data=water.data)
summary(ModelH20cm2)

# One water depth of zero!
water.data$Water_depth_cm[water.data$Water_depth_cm==0]=0.0001

#AICc comparison of linear and singular
Cand.models.depth=list()
Cand.models.depth[["Year+Season"]]=lm(I(log(Water_depth_cm))~I(Year)+Season,data=water.data)
Cand.models.depth[["Year^2+Season"]]=lm(I(log(Water_depth_cm))~I(Year)+I((Year)^2)+Season,data=water.data)
# Make AICc table.
AIC.Table.depth=model.sel(Cand.models.depth,rank="AICc")
AIC.Table.depth


# Make data.frame with predictions and CIs.
graph.depth=data.frame(predict(Cand.models.depth[["Year+Season"]],newdata=data.frame(Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8))),interval="confidence"),Year=c(2014:2021,2014:2021),Season=c(rep("Winter",8),rep("Summer",8)))

jpeg("Appendix3b_WaterParamaters_RHvent.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(5,5,2,2))
plot(2014:2021,graph.Water_Temp_C$fit[9:16],ylim=c(15,20),type="l",lty=2,ylab="Temperature (°C)",xlab="Year",las=1)
points(2014:2021,graph.Water_Temp_C$fit[9:16],pch=19)
lines((2014:2021)+0.1,graph.Water_Temp_C$fit[1:8],lty=3)
points((2014:2021)+0.1,graph.Water_Temp_C$fit[1:8],pch=19)
points((2014:2021)+0.1,graph.Water_Temp_C$lwr[1:8],pch=151)
points((2014:2021)+0.1,graph.Water_Temp_C$upr[1:8],pch=151)
points(2014:2021,graph.Water_Temp_C$lwr[9:16],pch=151)
points(2014:2021,graph.Water_Temp_C$upr[9:16],pch=151)
for (i in 1:8) {
  lines(rep(graph.Water_Temp_C$Year[i]+0.1,2),c(graph.Water_Temp_C$lwr[i],graph.Water_Temp_C$upr[i]))
  lines(rep(graph.Water_Temp_C$Year[i+8],2),c(graph.Water_Temp_C$lwr[i+8],graph.Water_Temp_C$upr[i+8]))
}
# legend("bottomright",legend=c("Summer","Winter"),lty=c(2,3))

plot(2014:2021,graph.DO$fit[9:16],ylim=c(0,12),type="l",lty=2,ylab=expression(~bold("Dissolved oxygen (mg L"^-1~")")),xlab="Year",las=1)
points(2014:2021,graph.DO$fit[9:16],pch=19)
lines((2014:2021)+0.1,graph.DO$fit[1:8],lty=3)
points((2014:2021)+0.1,graph.DO$fit[1:8],pch=19)
points((2014:2021)+0.1,graph.DO$lwr[1:8],pch=151)
points((2014:2021)+0.1,graph.DO$upr[1:8],pch=151)
points(2014:2021,graph.DO$lwr[9:16],pch=151)
points(2014:2021,graph.DO$upr[9:16],pch=151)
for (i in 1:8) {
  lines(rep(graph.DO$Year[i]+0.1,2),c(graph.DO$lwr[i],graph.DO$upr[i]))
  lines(rep(graph.DO$Year[i+8],2),c(graph.DO$lwr[i+8],graph.DO$upr[i+8]))
}
# legend("topright",legend=c("Summer","Winter"),lty=c(2,3))

plot(2014:2021,graph.pH$fit[9:16],ylim=c(6.6,7.4),type="l",lty=2,ylab="pH",xlab="Year",las=1)
points(2014:2021,graph.pH$fit[9:16],pch=19)
lines((2014:2021)+0.1,graph.pH$fit[1:8],lty=3)
points((2014:2021)+0.1,graph.pH$fit[1:8],pch=19)
points((2014:2021)+0.1,graph.pH$lwr[1:8],pch=151)
points((2014:2021)+0.1,graph.pH$upr[1:8],pch=151)
points(2014:2021,graph.pH$lwr[9:16],pch=151)
points(2014:2021,graph.pH$upr[9:16],pch=151)
for (i in 1:8) {
  lines(rep(graph.pH$Year[i]+0.1,2),c(graph.pH$lwr[i],graph.pH$upr[i]))
  lines(rep(graph.pH$Year[i+8],2),c(graph.pH$lwr[i+8],graph.pH$upr[i+8]))
}
# legend("bottomright",legend=c("Summer","Winter"),lty=c(2,3))

plot(2014:2021,graph.salt$fit[9:16],ylim=c(2,8),type="l",lty=2,ylab="Salinity (ppt)",xlab="Year",las=1)
points(2014:2021,graph.salt$fit[9:16],pch=19)
lines((2014:2021)+0.1,graph.salt$fit[1:8],lty=3)
points((2014:2021)+0.1,graph.salt$fit[1:8],pch=19)
points((2014:2021)+0.1,graph.salt$lwr[1:8],pch=151)
points((2014:2021)+0.1,graph.salt$upr[1:8],pch=151)
points(2014:2021,graph.salt$lwr[9:16],pch=151)
points(2014:2021,graph.salt$upr[9:16],pch=151)
for (i in 1:8) {
  lines(rep(graph.salt$Year[i]+0.1,2),c(graph.salt$lwr[i],graph.salt$upr[i]))
  lines(rep(graph.salt$Year[i+8],2),c(graph.salt$lwr[i+8],graph.salt$upr[i+8]))
}
# legend("bottomright",legend=c("Summer","Winter"),lty=c(2,3))

plot(2014:2021,exp(graph.depth$fit[9:16]),ylim=c(0,20),type="l",lty=2,ylab="Water depth (cm)",xlab="Year",las=1)
points(2014:2021,exp(graph.depth$fit[9:16]),pch=19)
lines((2014:2021)+0.1,exp(graph.depth$fit[1:8]),lty=3)
points((2014:2021)+0.1,exp(graph.depth$fit[1:8]),pch=19)
points((2014:2021)+0.1,exp(graph.depth$lwr[1:8]),pch=151)
points((2014:2021)+0.1,exp(graph.depth$upr[1:8]),pch=151)
points(2014:2021,exp(graph.depth$lwr[9:16]),pch=151)
points(2014:2021,exp(graph.depth$upr[9:16]),pch=151)
for (i in 1:8) {
  lines(rep(graph.depth$Year[i]+0.1,2),c(exp(graph.depth$lwr[i]),exp(graph.depth$upr[i])))
  lines(rep(graph.depth$Year[i+8],2),c(exp(graph.depth$lwr[i+8]),exp(graph.depth$upr[i+8])))
}
legend("topleft",legend=c("Summer","Winter"),lty=c(2,3))
dev.off()

