#read data in
aymf <- read.csv("~/Dropbox (Yale_FES)/Angus_Herps/Drafts/PeerJ submit/Code and Data Files/YMF_Herp_Data_Update.csv", header = T)
View(aymf)

dim(aymf)
#[1] 14 17
library(YaleToolkit)
whatis(aymf)

row.names(aymf) <- aymf$Stand_Name
row.names(aymf)

#check for correlation to chose variables
variables=as.matrix(aymf[, c(4:16)])
corr <- cor(variables)
round(corr, 2)

library(MASS)
library(arm)
library(lm.beta)
library(car)
library(rsq)

###red backed salamander model
RBS <- glm.nb(Red_Backed_Salamander~Year_Since_Cut+Volume_m3_per_ha+number_sapling_per_ha+Overstory_basal_area, data = aymf)
summary(RBS)
lm.beta(RBS)
vif(RBS)
rsq(RBS)

Newt <- glm.nb(Eastern_Newt~Year_Since_Cut+Volume_m3_per_ha+number_sapling_per_ha+Overstory_basal_area, data = aymf)
summary(Newt)
lm.beta(Newt)
vif(Newt)
rsq(Newt)


###### Graphically exploring some abiotic relatoinships
  #regression of pieces of CWD against stand age
m1 <- lm(pieces_of_wood ~ Year_Since_Cut, data = aymf)
summary(m1)

  #regresion of decay class against stand age
m2 <- lm(Average_Wood_Decay_Class_per_Stand ~ Year_Since_Cut, data = aymf)
summary(m2) #


###### Red-Backed Salamander Plots (syntax for plot annotated on first plot only)

rbs_p_1 <- ggplot(aymf, aes(x=Year_Since_Cut, y = Red_Backed_Salamander)) + #set up the plot variables
  geom_point()+ #add points
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) # remove gridlines
rbs_p_1
rbs_p_2 <- rbs_p_1 + theme_bw()+ # create a white background
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + #play with x-axis text style
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) + # play wiht y-axis text style
  labs(x="Year Since Cut", y="# of Red-Backed Salamander")+ #axis labels
  geom_hline(yintercept=c(8.4), linetype="dotted")+ # add horizontal line
  ylim(0,40)+ #set y-axis limits
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black") #insert poisson regression curve
rbs_p_2


  #regression of CWD against Years Since Cut
year_wood <- lm(pieces_of_wood ~ Year_Since_Cut, data = aymf)
resid.year_wood <- resid(year_wood) #extracting residuals from this regression


  #plotting red-backed salamander abundance against the residuals of CWD:years
rbs_p_1b <- ggplot(aymf, aes(x=resid.year_wood, y = Red_Backed_Salamander)) + geom_point()+
theme(panel.grid.major=element_blank(),
      panel.grid.minor=element_blank())
rbs_p_2b <- rbs_p_1b + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="residuals(Pieces of Course Woody Debris ~ Year Since Cut)", y="# of Red-Backed Salamander")+
  ylim(0,40)+
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
rbs_p_2b


  #plotting CWD against Years Since Cut
rbs_p_1c <- ggplot(aymf, aes(x=Year_Since_Cut, y = pieces_of_wood)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
rbs_p_2c <- rbs_p_1c + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Years Since Cut", y="Pieces of Course Woody Debris")+
  geom_smooth(method = "lm", colour = "black")
rbs_p_2c


  #plotting RBS against CWD
rbs_p_3 <- ggplot(aymf, aes(x=pieces_of_wood, y = Red_Backed_Salamander)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
rbs_p_4 <- rbs_p_3 + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Pieces of Course Woody Debris", y="# of Red-Backed Salamander")+
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
rbs_p_4



  #plotting RBS against basal area
rbs_p_5 <- ggplot(aymf, aes(x=Overstory_basal_area, y = Red_Backed_Salamander)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
rbs_p_6 <- rbs_p_5 + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Overstory Basal Area (m2 / ha)", y="# of Red-Backed Salamander")+
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
rbs_p_6



### Newt Plots

  #plotting Newts against saplings
en_p_1 <- ggplot(aymf, aes(x=number_sapling_per_ha, y = Eastern_Newt)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
en_p_2 <- en_p_1 + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Saplings per Hectare", y="# of Eastern Newts")+ ylim(0,15)+
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
en_p_2


  #plotting newts against basal area
en_p_3 <- ggplot(aymf, aes(x=Overstory_basal_area, y = Eastern_Newt)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
en_p_4 <- en_p_3 + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Overstory Basal Area (m2 / ha)", y="# of Eastern Newts")+ylim(0,15)+
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
en_p_4


  #plotting newts against stand age
en_p_5 <- ggplot(aymf, aes(x=Year_Since_Cut, y = Eastern_Newt)) + geom_point()+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())
en_p_6 <- en_p_5 + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.x = element_text(face="bold", size=15),
        axis.text.x  = element_text(size=12)) + 
  theme(axis.title.y = element_text(face="bold", size=17, vjust=0.1),
        axis.text.y  = element_text(size=10)) +
  labs(x="Year Since Cut", y="# of Eastern Newts")+ 
  geom_hline(yintercept=c(8.2), linetype="dotted")+ylim(0,40) +
  geom_smooth(method = "glm", method.args=list(family="poisson"), colour = "black")
en_p_6