# Load libraries
library(vegan)
library(clustsig)
library(car)

# load the data without adding a column, no row names
short_term<-read.csv("intertidal_short_term.csv",header=T)
# make a new object out of the first column, turn it into text
grps<-as.character(short_term[,1])
# parse the values, separating the dat and replicate number
grp_matrix<-as.data.frame(matrix(unlist(strsplit(grps,split = "-")),nrow=42,byrow = T))
# remove the first column
community<-short_term[,-c(1)]

# nmds
nmds<-metaMDS(community,k=2,autotransform = F)
# bray curtis disance
distance<-vegdist(community, method="bray")
# ordination distance
ord_dist<-vegdist(nmds$points[,1:2],method="euclidean")
# estimate variance explained
cor(distance,ord_dist)^2

# plot
# colors and symbols for 14 dates
symb_colors<-c(rainbow(7),rainbow(7))
symb<-(rep(c(22,16),7))

plot(rbind(nmds$points[,1:2],nmds$species[,1:2]),type="n",xlab="Axis 1", ylab="Axis 2")
points(nmds$points[,1:2],pch=symb[as.factor(grp_matrix$V1)],cex=1.1,col=symb_colors[as.factor(grp_matrix$V1)])
text(nmds$species,labels=make.cepnames(row.names(nmds$species)),col="red")
legend("topright",legend = levels(as.factor(grp_matrix$V1)),pch=symb,col=symb_colors, ncol=3, bty = "n")
# simprof precedure as described in the paper
simp_res<-simprof(community,method.distance="braycurtis")
# put polygons around the groups from simprof
for(g in 1:simp_res$numgroups)
{
  p<-nmds$points[unlist(simp_res$significantclusters[g]),1:2]
  ch<-chull(p[,1],p[,2])
  polygon(p[ch,1],p[ch,2])
}


# alternative without using simprof
upgma_cluster<-hclust(vegdist(temp), method="average")
for(g in 1:3)
{
  p<-nmds$points[(cutree(upgma_cluster,k=3)==g),1:2]
  ch<-chull(p[,1],p[,2])
  polygon(p[ch,1],p[ch,2])
}

# perform the ANOSIM
anosim(community,as.factor(grp_matrix$V1))

# stress value
nmds$stress
# what exactly is this?
# Extract the residuals from a monotonic regression (recall nmds works with ranks)
dhat<-isoreg(distance, ord_dist)
# plot the monotonic regression
plot(dhat)
# this is the same plot you get from this function
stressplot(nmds)
# The non-metric fit is 1-stress^2
1-nmds$stress^2
# The linear fit is a simple correlation
cor(distance,ord_dist)
# The goodness of fit values (GOF) are also directly relatated to stress
# sum(GOF^2) = Stress^2
sum(goodness(nmds)^2)
nmds$stress^2

# % variation
cor(ord_dist<-dist(nmds$points),vegdist(community))^2
# axis 1
cor(ord_dist<-dist(nmds$points[,1]),vegdist(community))^2
# axis 2
cor(ord_dist<-dist(nmds$points[,2]),vegdist(community))^2