## Code for all fish collected

Fish = read.csv(file.choose(), header = TRUE)
###use file named All fish nMDS data.csv

Fish_community = Fish[2:74]
#columns 2-74 include all species collected
#all 15 quadrats included


library(vegan)

set.seed(2)
example_NMDS=metaMDS(Fish_community, 
                     k=2) 

##########FOR nMDS###############################

plot(example_NMDS)
ordiplot(example_NMDS,type="n")
##orditorp(example_NMDS,display="species",col="black",air=.01, cex=.75)
orditorp(example_NMDS,display="sites", pch = 2, col = "blue", cex=.7)


treat=c(rep("Treatment1",5),rep("Treatment2",5),rep("Treatment3",5))
ordi(example_NMDS,type="n")
ordiellipse(example_NMDS,groups=treat,draw="polygon",col="grey90",label=F)
###orditorp(example_NMDS,display="species",col="black",air=0.01, cex=.6)
orditorp(example_NMDS,display="sites",col=c(rep("red",5),rep("tan2",5),rep("gold",5)),
         air=0.01,cex=.7)

colors=c(rep("red",5),rep("tan2",5),rep("gold",5))
ordiplot(example_NMDS,type="n")
for(i in unique(treat)) {
  ordiellipse(example_NMDS$point[grep(i,treat),],draw="polygon",
           groups=treat[treat==i],col=colors[grep(i,treat)],label=F) } 
###orditorp(example_NMDS,display="species",col="black",air=0.01, cex=.75)
orditorp(example_NMDS,display="sites",col=c(rep("red",5),
                                            rep("tan2",5),rep("gold",5)),air=0.01,cex=1.25)



with(Fish, levels(Habitat))
scl <- 3 ## scaling = 3
colvec <- c("red3", "tan4", "darkgoldenrod1")
plot(example_NMDS)
ordiplot(example_NMDS,type="n")


treat=c(rep("Treatment1",5),rep("Treatment2",5),rep("Treatment3",5))
ordiplot(example_NMDS,type="n")
ordiellipse(example_NMDS,groups=treat,draw="polygon", conf= .95,col="white",label=F)
###orditorp(example_NMDS,display="species",col="black",air=0.01, cex=.75)

colors=c(rep("red",5),rep("tan2",5),rep("gold",5))
## TO edit plot dimensions##
ordiplot(example_NMDS,type="n", xlim = c(-1.25, 1.25), ylim = c(-1.25,1.25))
for(i in unique(treat)) {
  ordiellipse(example_NMDS$point[grep(i,treat),],draw="polygon", conf = .95,
           groups=treat[treat==i],col=colors[grep(i,treat)],label=F) } 

with(Fish, points(example_NMDS, display = "sites", col = colvec[Habitat],
                  scaling = scl, pch = 21, bg = colvec[Habitat]))
head(with(Fish, colvec[Habitat]))

with(Fish, legend("topright", legend = levels(Habitat), bty = "n",
                  col = colvec, pch = 21, pt.bg = colvec))
##orditorp(example_NMDS,display="species",col="black",air=0.01, cex=.75)



################# STRESS ######################################

example_NMDS$stress

################# COMMUNITY METRICS ##################

H = diversity(Fish_community, index = "shannon") #Shannon's Diversity Index H'
S = specnumber(Fish_community) #Species Richness S
J = H/log(S) #Species Evenness J

H
S
J


#########################################

#PERMANOVA TEST

adonis(Fish_community ~ Habitat, data= Fish, permutations=999)