setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/HondoManuscript_AquaticConservation_Code-Data/Analysis1and5/Analyes1and5_27Feb22/Analyses_RioHondo_Midstream")
# Analysis 1 Springsnails in Midstream
# Analysis 1 dataset.
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis1_InvertDensity_RioHondo_ThroughTime_Midstream19Nov21.csv",header=T,sep=",")
# Set the working directory.
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/HondoManuscript_AquaticConservation_Code-Data/Analysis1and5/Analyses1and5_27Feb22/Analyses_RioHondo_Midstream")
# Analysis 1 Springsnails in Midstream
# Analysis 1 dataset.
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis1_InvertDensity_RioHondo_ThroughTime_Midstream19Nov21.csv",header=T,sep=",")
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
names(data)
# Springsnails in Rio Hondo Midstream.
# Quadratic springsnaiis
summary(glmer(Springsnails~I(Year-2015)+I((Year-2015)^2)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
summary(glmer(Springsnails~I(Year-2015)+Season+Subsample+offset(log(Trap_size_sq_m))+(1|Sample),data=data,family=poisson(link="log")))
# Make list of candidate models.
Cand.models.SS=list()
Cand.models.SS[["Year"]]=glmer(Springsnails~I(Year-2015)+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-2015)+I((Year-2015)^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-2014))+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 quadratic model.
model1=glmer(Springsnails~I(Year-2015)+I((Year-2015)^2)+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.SS=data.frame(Subsample=rep("A",14),Year=rep(seq(2015,2021,1),2),Season=c(rep("Winter",7),rep("Summer",7)),Trap_size_sq_m=rep(1,14),Sample="")
GRAPH.DATA.SS$Pred=predict(model1,newdata=GRAPH.DATA.SS,type="response",allow.new.levels=TRUE)
# Bootstrap the predictions in order to get CIs.This will take a very long time.
model1.boot=bootMer(model1,function(x) predict(x,newdata=GRAPH.DATA.SS,type="response",allow.new.levels=TRUE),nsim=1000)
# Get 95% CIs.
model1.boot.CI=data.frame(t(apply(model1.boot$t[!is.na(model1.boot$t[,1]),],2,quantile,c(0.025,0.975))))
names(model1.boot.CI)=c("LCL","UCL")
# Make and save a high resolution figure.
jpeg("SS_Trend_RHmidstream.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,5,2,2))
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",legend=c("Summer","Winter"),lty=c(2,3))
dev.off()
# Multi comparisons springsnails
summary(model1)
emmeans(model1,pairwise~Season,offset=1,type="response",adjust="tukey",level=0.95)
# plot of Multiple comparison test: Season
TUKEYSpringsnails=emmeans(model1,list(pairwise~Season),offset=1,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
emmeans(model1,pairwise~Subsample,offset=1,type="response",adjust="tukey",level=0.95)
TUKEYSpringsnails1=emmeans(model1,list(pairwise~Subsample),offset=1,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_RHmidstream.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(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,0.02),xlim=c(0.5,3.5),xlab="Channel Location")
box()
# axis make new axis with customizations.
# yaxp=c(0,1.6,8) sets the ticks.
axis(2,ylim=c(0,0.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=3,las=1,side=3,line=+21.0, cex=1.2,font=1)
mtext("(b)",col="black",at=3,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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=4,las=1,side=3,line=+15.0, cex=1.2,font=1)
mtext("(b)",col="black",at=4,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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=3.5,las=1,side=3,line=+15.0, cex=1.2,font=1)
mtext("(b)",col="black",at=3.5,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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=3.3,las=1,side=3,line=+15.8, cex=1.2,font=1)
mtext("(b)",col="black",at=3.3,las=1,side=3,line=-1.4, 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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
# Make and save a high resolution figure.
jpeg("Figure5_Springsnails_RHmidstream.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.35)
plot(2015:2021,GRAPH.DATA.SS$Pred[8:14],ylim=c(0,3),type="l",lty=2,ylab=expression(~bold("Springsnails (no. cm"^-2~")")),xlab="Year",las=1)
points(2015:2021,GRAPH.DATA.SS$Pred[8:14],pch=19)
lines((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],lty=3)
points((2015:2021)+0.1,GRAPH.DATA.SS$Pred[1:7],pch=19)
points((2015:2021)+0.1,model1.boot.CI$LCL[1:7],pch=151)
points((2015:2021)+0.1,model1.boot.CI$UCL[1:7],pch=151)
points(2015:2021,model1.boot.CI$LCL[8:14],pch=151)
points(2015:2021,model1.boot.CI$UCL[8:14],pch=151)
for (i in 1:7) {
lines(rep(GRAPH.DATA.SS$Year[i]+0.1,2),c(model1.boot.CI$LCL[i],model1.boot.CI$UCL[i]))
lines(rep(GRAPH.DATA.SS$Year[i+7],2),c(model1.boot.CI$LCL[i+7],model1.boot.CI$UCL[i+7]))
}
legend("topleft",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.02),xlim=c(0.5,3.5),xlab="Sample Location")
mtext("(a)",col="black",at=3.4,las=1,side=3,line=+15.6, cex=1.2,font=1)
mtext("(b)",col="black",at=3.4,las=1,side=3,line=-1.35, 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.02),col="black",las=1,yaxp=c(0,0.02,5))
axis(1,at=1:3,labels=c("Bank","Quarter","Center"),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()
#Analsysis 5 Water Parameters for Rio Hondo midstream
# Designate data set for water parameters.
water.data=read.csv ("Analysis5_WaterParameters_ThroughTime_RioHondo_Midstream_20Aug21.csv", header=TRUE)
names(water.data)
# A Comparison of water parameter across springruns 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')
#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+Season"]],newdata=data.frame(Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7))),interval="confidence"),Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7)))
#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(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7))),interval="confidence"),Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7)))
#Investigate change in pH parameter in Hondo Springruns 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+Season"]],newdata=data.frame(Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7))),interval="confidence"),Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7)))
#Investigate change in Salinity parameter Hondo Springruns 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(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7))),interval="confidence"),Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7)))
#Investigate change in Water Depth parameter Hondo Springruns using lm
ModelH20cm=lm(I(log(Water_depth_cm))~I(Year)+I((Year)^2)+Season,data=water.data)
summary(ModelH20cm)
ModelH20cm2=lm(I(log(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(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7))),interval="confidence"),Year=c(2015:2021,2015:2021),Season=c(rep("Winter",7),rep("Summer",7)))
jpeg("Appendix3a_Water_RHmidstream.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(2015:2021,graph.Water_Temp_C$fit[8:14],ylim=c(10,25),type="l",lty=2,ylab="Temperature (°C)",xlab="Year",las=1)
points(2015:2021,graph.Water_Temp_C$fit[8:14],pch=19)
lines((2015:2021)+0.1,graph.Water_Temp_C$fit[1:7],lty=3)
points((2015:2021)+0.1,graph.Water_Temp_C$fit[1:7],pch=19)
points((2015:2021)+0.1,graph.Water_Temp_C$lwr[1:7],pch=151)
points((2015:2021)+0.1,graph.Water_Temp_C$upr[1:7],pch=151)
points(2015:2021,graph.Water_Temp_C$lwr[8:14],pch=151)
points(2015:2021,graph.Water_Temp_C$upr[8:14],pch=151)
for (i in 1:7) {
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+7],2),c(graph.Water_Temp_C$lwr[i+7],graph.Water_Temp_C$upr[i+7]))
}
plot(2015:2021,graph.DO$fit[8:14],ylim=c(0,20),type="l",lty=2,ylab=expression(~bold("Dissolved oxygen (mg L"^-1~")")),xlab="Year",las=1)
points(2015:2021,graph.DO$fit[8:14],pch=19)
lines((2015:2021)+0.1,graph.DO$fit[1:7],lty=3)
points((2015:2021)+0.1,graph.DO$fit[1:7],pch=19)
points((2015:2021)+0.1,graph.DO$lwr[1:7],pch=151)
points((2015:2021)+0.1,graph.DO$upr[1:7],pch=151)
points(2015:2021,graph.DO$lwr[8:14],pch=151)
points(2015:2021,graph.DO$upr[8:14],pch=151)
for (i in 1:7) {
lines(rep(graph.DO$Year[i]+0.1,2),c(graph.DO$lwr[i],graph.DO$upr[i]))
lines(rep(graph.DO$Year[i+7],2),c(graph.DO$lwr[i+7],graph.DO$upr[i+7]))
}
plot(2015:2021,graph.pH$fit[8:14],ylim=c(7.2,7.6),type="l",lty=2,ylab="pH",xlab="Year",las=1)
points(2015:2021,graph.pH$fit[8:14],pch=19)
lines((2015:2021)+0.1,graph.pH$fit[1:7],lty=3)
points((2015:2021)+0.1,graph.pH$fit[1:7],pch=19)
points((2015:2021)+0.1,graph.pH$lwr[1:7],pch=151)
points((2015:2021)+0.1,graph.pH$upr[1:7],pch=151)
points(2015:2021,graph.pH$lwr[8:14],pch=151)
points(2015:2021,graph.pH$upr[8:14],pch=151)
for (i in 1:7) {
lines(rep(graph.pH$Year[i]+0.1,2),c(graph.pH$lwr[i],graph.pH$upr[i]))
lines(rep(graph.pH$Year[i+7],2),c(graph.pH$lwr[i+7],graph.pH$upr[i+7]))
}
plot(2015:2021,graph.salt$fit[8:14],ylim=c(3,9),type="l",lty=2,ylab="Salinity (ppt)",xlab="Year",las=1)
points(2015:2021,graph.salt$fit[8:14],pch=19)
lines((2015:2021)+0.1,graph.salt$fit[1:7],lty=3)
points((2015:2021)+0.1,graph.salt$fit[1:7],pch=19)
points((2015:2021)+0.1,graph.salt$lwr[1:7],pch=151)
points((2015:2021)+0.1,graph.salt$upr[1:7],pch=151)
points(2015:2021,graph.salt$lwr[8:14],pch=151)
points(2015:2021,graph.salt$upr[8:14],pch=151)
for (i in 1:7) {
lines(rep(graph.salt$Year[i]+0.1,2),c(graph.salt$lwr[i],graph.salt$upr[i]))
lines(rep(graph.salt$Year[i+7],2),c(graph.salt$lwr[i+7],graph.salt$upr[i+7]))
}
plot(2015:2021,exp(graph.depth$fit[8:14]),ylim=c(0,25),type="l",lty=2,ylab="Water depth (cm)",xlab="Year",las=1)
points(2015:2021,exp(graph.depth$fit[8:14]),pch=19)
lines((2015:2021)+0.1,exp(graph.depth$fit[1:7]),lty=3)
points((2015:2021)+0.1,exp(graph.depth$fit[1:7]),pch=19)
points((2015:2021)+0.1,exp(graph.depth$lwr[1:7]),pch=151)
points((2015:2021)+0.1,exp(graph.depth$upr[1:7]),pch=151)
points(2015:2021,exp(graph.depth$lwr[8:14]),pch=151)
points(2015:2021,exp(graph.depth$upr[8:14]),pch=151)
for (i in 1:7) {
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+7],2),c(exp(graph.depth$lwr[i+7]),exp(graph.depth$upr[i+7])))
}
legend("topleft",legend=c("Summer","Winter"),lty=c(2,3))
dev.off()
