par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6.5,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot1$data$the.emmean,axes=FALSE,pch="",main="Winter", ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,20),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,22,22,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSpirngs","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$asymp.LCL[i],TUKEYSubsample.plot1$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot1$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$rcmpl[i]),lty=1)
}
dev.off()
library(lme4)
library(lsmeans)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)
# Set the working directory.
#setwd("C:/A_USFWS/Bitter Lake NWR/Rio Hondo MS/Analysis3/")
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/Analysis3/Analsys3_27Feb22")
# Analysis 3 dataset 1 - Summer
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Summers_19Nov21.csv",header=T,sep=",")
names(data)
unique(data$Site)
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 8402.8
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 15703.2
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data,poisson(link="log")))
#no year since not part of objective. AIC is 8450.6
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log")))
#includes random effect of sample location and year. Best AIC from above three comparisons.
model2=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log"))
summary(model2)
# Multiple comparison test
emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample=emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot=plot(TUKEYSubsample,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Analysis 3 dataset 2 - Winter
# Bring data to be analyzed into RStudio.
data2=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Winters_19Nov21.csv",header=T,sep=",")
names(data2)
unique(data2$Site)
# Correct trap size.
data2$Trap_size_sq_m[data2$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 3495.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 7314.3
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
#no year since not part of objective. AIC is 3574.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data2,poisson(link="log")))
#Sample and year are random effects. AIC = 3495.7. No Subsample due to difference in collection at RH Springvents and all others; sample should capture subsample variance.
Model3=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
summary(Model3)
#Year is only random effect (AIC=7314.3), since year is not part of objective & No Subsample due to difference in collection at RH Springvents and all others.
Model4=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
summary(Model4)
#No Year as random or fixed effect. AIC = 3574.7
model5=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log"))
summary(model5)
# Multiple comparison test
emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample1=emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot1=plot(TUKEYSubsample1,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Save the plot
jpeg("Figure6_Springsnail_SpringSystems_summer_winter.jpg",width=4.5,height=6.5,units="in",res=675)
layout(rbind(1,2)); layout.show(2)
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot$data$the.emmean,axes=FALSE,pch="",main="Summer",ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,30),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,32,32,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSpirngs","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$asymp.LCL[i],TUKEYSubsample.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$rcmpl[i]),lty=1)
}
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
library(lme4)
library(lsmeans)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)
# Set the working directory.
#setwd("C:/A_USFWS/Bitter Lake NWR/Rio Hondo MS/Analysis3/")
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/Analysis3/Analsys3_27Feb22")
# Analysis 3 dataset 1 - Summer
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Summers_19Nov21.csv",header=T,sep=",")
names(data)
unique(data$Site)
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 8402.8
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 15703.2
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data,poisson(link="log")))
#no year since not part of objective. AIC is 8450.6
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log")))
#includes random effect of sample location and year. Best AIC from above three comparisons.
model2=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log"))
summary(model2)
# Multiple comparison test
emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample=emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot=plot(TUKEYSubsample,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Analysis 3 dataset 2 - Winter
# Bring data to be analyzed into RStudio.
data2=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Winters_19Nov21.csv",header=T,sep=",")
names(data2)
unique(data2$Site)
# Correct trap size.
data2$Trap_size_sq_m[data2$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 3495.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 7314.3
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
#no year since not part of objective. AIC is 3574.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data2,poisson(link="log")))
#Sample and year are random effects. AIC = 3495.7. No Subsample due to difference in collection at RH Springvents and all others; sample should capture subsample variance.
Model3=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
summary(Model3)
#Year is only random effect (AIC=7314.3), since year is not part of objective & No Subsample due to difference in collection at RH Springvents and all others.
Model4=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
summary(Model4)
#No Year as random or fixed effect. AIC = 3574.7
model5=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log"))
summary(model5)
# Multiple comparison test
emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample1=emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot1=plot(TUKEYSubsample1,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Save the plot
jpeg("Figure6_Springsnail_SpringSystems_summer_winter.jpg",width=4.5,height=6.5,units="in",res=675)
layout(rbind(1,2)); layout.show(2)
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot$data$the.emmean,axes=FALSE,pch="",main="Summer",ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,30),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,32,32,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSpirngs","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$asymp.LCL[i],TUKEYSubsample.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$rcmpl[i]),lty=1)
}
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6.5,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot1$data$the.emmean,axes=FALSE,pch="",main="Winter", ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,20),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,22,22,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
axis(1,at=1:5,labels=c("Bitter\nCreek","Rio Hondo\nMidstream","Rio Hondo\nSpring-vents","Sago\nSpirngs","Snail\nUnit"),col="black",las=2,cex=2)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$asymp.LCL[i],TUKEYSubsample.plot1$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot1$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$rcmpl[i]),lty=1)
}
dev.off()
library(lme4)
library(lsmeans)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)
# Set the working directory.
#setwd("C:/A_USFWS/Bitter Lake NWR/Rio Hondo MS/Analysis3/")
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/Analysis3/Analsys3_27Feb22")
# Analysis 3 dataset 1 - Summer
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Summers_19Nov21.csv",header=T,sep=",")
names(data)
unique(data$Site)
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 8402.8
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 15703.2
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data,poisson(link="log")))
#no year since not part of objective. AIC is 8450.6
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log")))
#includes random effect of sample location and year. Best AIC from above three comparisons.
model2=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log"))
summary(model2)
# Multiple comparison test
emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample=emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot=plot(TUKEYSubsample,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Analysis 3 dataset 2 - Winter
# Bring data to be analyzed into RStudio.
data2=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Winters_19Nov21.csv",header=T,sep=",")
names(data2)
unique(data2$Site)
# Correct trap size.
data2$Trap_size_sq_m[data2$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 3495.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 7314.3
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
#no year since not part of objective. AIC is 3574.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data2,poisson(link="log")))
#Sample and year are random effects. AIC = 3495.7. No Subsample due to difference in collection at RH Springvents and all others; sample should capture subsample variance.
Model3=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
summary(Model3)
#Year is only random effect (AIC=7314.3), since year is not part of objective & No Subsample due to difference in collection at RH Springvents and all others.
Model4=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
summary(Model4)
#No Year as random or fixed effect. AIC = 3574.7
model5=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log"))
summary(model5)
# Multiple comparison test
emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample1=emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot1=plot(TUKEYSubsample1,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Save the plot
jpeg("Figure6_Springsnail_SpringSystems_summer_winter.jpg",width=4.5,height=6.5,units="in",res=675)
layout(rbind(1,2)); layout.show(2)
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot$data$the.emmean,axes=FALSE,pch="",main="Summer",ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,30),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,32,32,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
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(TUKEYSubsample.plot$data$asymp.LCL[i],TUKEYSubsample.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$rcmpl[i]),lty=1)
}
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6.5,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot1$data$the.emmean,axes=FALSE,pch="",main="Winter", ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,20),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,22,22,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
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(TUKEYSubsample.plot1$data$asymp.LCL[i],TUKEYSubsample.plot1$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot1$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$rcmpl[i]),lty=1)
}
dev.off()
library(lme4)
library(lsmeans)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)
# Set the working directory.
#setwd("C:/A_USFWS/Bitter Lake NWR/Rio Hondo MS/Analysis3/")
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/Analysis3/Analsys3_27Feb22")
# Analysis 3 dataset 1 - Summer
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Summers_19Nov21.csv",header=T,sep=",")
names(data)
unique(data$Site)
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 8402.8
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 15703.2
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data,poisson(link="log")))
#no year since not part of objective. AIC is 8450.6
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log")))
#includes random effect of sample location and year. Best AIC from above three comparisons.
model2=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log"))
summary(model2)
# Multiple comparison test
emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample=emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot=plot(TUKEYSubsample,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Analysis 3 dataset 2 - Winter
# Bring data to be analyzed into RStudio.
data2=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Winters_19Nov21.csv",header=T,sep=",")
names(data2)
unique(data2$Site)
# Correct trap size.
data2$Trap_size_sq_m[data2$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 3495.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 7314.3
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
#no year since not part of objective. AIC is 3574.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data2,poisson(link="log")))
#Sample and year are random effects. AIC = 3495.7. No Subsample due to difference in collection at RH Springvents and all others; sample should capture subsample variance.
Model3=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
summary(Model3)
library(lme4)
library(lsmeans)
library(emmeans)
library(lmerTest)
library(MuMIn)
library(ggplot2)
library(glm2)
# Set the working directory.
#setwd("C:/A_USFWS/Bitter Lake NWR/Rio Hondo MS/Analysis3/")
setwd("C:/Users/bpjohnson/Desktop/RioHondo/VentCompare/Analysis3/Analsys3_27Feb22")
# Analysis 3 dataset 1 - Summer
# Bring data to be analyzed into RStudio.
data=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Summers_19Nov21.csv",header=T,sep=",")
names(data)
unique(data$Site)
# Correct trap size.
data$Trap_size_sq_m[data$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 8402.8
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 15703.2
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data,poisson(link="log")))
#no year since not part of objective. AIC is 8450.6
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log")))
#includes random effect of sample location and year. Best AIC from above three comparisons.
model2=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data,poisson(link="log"))
summary(model2)
# Multiple comparison test
emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample=emmeans(model2,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot=plot(TUKEYSubsample,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Analysis 3 dataset 2 - Winter
# Bring data to be analyzed into RStudio.
data2=read.csv("Analysis3_InvertDensity-BySpringsystem_HondoManuscript_Winters_19Nov21.csv",header=T,sep=",")
names(data2)
unique(data2$Site)
# Correct trap size.
data2$Trap_size_sq_m[data2$Trap_size_sq_m==0.003]=(5.08*3.81)
#Sample and Year are random effects, since not part of objectives. AIC is 3495.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
#Sample and Year are random effects, since not part of objectives. AIC is 7314.3
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
#no year since not part of objective. AIC is 3574.7
summary(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data2,poisson(link="log")))
#Sample and year are random effects. AIC = 3495.7. No Subsample due to difference in collection at RH Springvents and all others; sample should capture subsample variance.
Model3=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample)+(1|Year),data=data2,poisson(link="log")))
summary(Model3)
#Year is only random effect (AIC=7314.3), since year is not part of objective & No Subsample due to difference in collection at RH Springvents and all others.
Model4=(glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Year),data=data2,poisson(link="log")))
summary(Model4)
#No Year as random or fixed effect. AIC = 3574.7
model5=glmer(Springsnails~Site+offset(log(Trap_size_sq_m))+(1|Sample),data=data,poisson(link="log"))
summary(model5)
# Multiple comparison test
emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# plot tukeys for Site
TUKEYSubsample1=emmeans(Model3,list(pairwise~Site),offset=1,type="response",adjust="tukey",level=0.95)
# Get the info for plotting.
TUKEYSubsample.plot1=plot(TUKEYSubsample1,comparison=TRUE,plotIt=FALSE,adjust="tukey")
# Save the plot
jpeg("Figure6_Springsnail_SpringSystems_summer_winter.jpg",width=4.5,height=6.5,units="in",res=675)
layout(rbind(1,2)); layout.show(2)
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot$data$the.emmean,axes=FALSE,pch="",main="Summer",ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,30),xlim=c(0.5,5.5),xlab="")
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,32,32,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
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(TUKEYSubsample.plot$data$asymp.LCL[i],TUKEYSubsample.plot$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot$data$the.emmean[i],TUKEYSubsample.plot$data$rcmpl[i]),lty=1)
}
par(cex.lab=1.10)
par(font.lab=2)
par(cex.main=1.3)
par(font.main=2)
par(cex.axis=1.09)
par(lwd=2)
par(mar=c(6.5,5,2,2))
# axes=FALSE takes away the axis so you can customize them with code below.
plot(1:5,TUKEYSubsample.plot1$data$the.emmean,axes=FALSE,pch="",main="Winter", ylab=expression(~bold("Springsnails (no. cm"^-2~")")),ylim=c(0,20),xlim=c(0.5,5.5),xlab="")
mtext("(a)",col="black",at=0.61,las=1,side=3,line=+15.2, cex=0.83,font=1)
mtext("(b)",col="black",at=0.61,las=1,side=3,line=-1.15, cex=0.83,font=1)
box()
polygon(x=c(1.5,1.5,3.5,3.5),y=c(-3,22,22,-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(0,30),col="black",las=1,yaxp=c(0,30,6))
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(TUKEYSubsample.plot1$data$asymp.LCL[i],TUKEYSubsample.plot1$data$asymp.UCL[i]),lty=1,lwd=15,col="grey")
}
points(1:5,TUKEYSubsample.plot1$data$the.emmean,pch=19)
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$lcmpl[i]),lty=1)
}
for (i in 1:5){
lines(rep(i,2),c(TUKEYSubsample.plot1$data$the.emmean[i],TUKEYSubsample.plot1$data$rcmpl[i]),lty=1)
}
dev.off()
