Goby = read.csv(file.choose(), header = TRUE) Goby_community = Goby[2:32] #columns 2-32 include all 31 species collected #13 quadrats included (2 dropped due to low abundance) #use file named Goby nMDS data.csv library(vegan) set.seed(2) example_NMDS=metaMDS(Goby_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",4),rep("Treatment2",5),rep("Treatment3",4)) ordi(example_NMDS,type="n") ordihull(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",4),rep("tan2",5),rep("gold",4)), air=0.01,cex=.7) colors=c(rep("red",4),rep("tan2",5),rep("gold",4)) ordiplot(example_NMDS,type="n") for(i in unique(treat)) { ordihull(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",4), rep("tan2",5),rep("gold",4)),air=0.01,cex=1.25) with(Goby, levels(Habitat)) scl <- 3 ## scaling = 3 colvec <- c("red3", "tan4", "darkgoldenrod1") plot(example_NMDS) ordiplot(example_NMDS,type="n") treat=c(rep("Treatment1",4),rep("Treatment2",5),rep("Treatment3",4)) 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",4),rep("tan2",5),rep("gold",4)) ordiplot(example_NMDS,type="n", xlim = c(-1.5, 1.5), ylim = c(-1.5,1.5)) 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(Goby, points(example_NMDS, display = "sites", col = colvec[Habitat], scaling = scl, pch = 21, bg = colvec[Habitat])) head(with(Goby, colvec[Habitat])) with(Goby, 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(Goby_community, index = "shannon") #Shannon's Diversity Index H' S = specnumber(Goby_community) #Species Richness S J = H/log(S) #Species Evenness J H S J ######################################### #PERMANOVA adonis(Goby_community ~ Habitat, data= Goby, permutations=999)