#PERMANOVA AND SIMPER ANALYSIS FOR SPECIES CHARACTERISTICS. CATEGORICAL VARIABLES NEED TO BE ASSIGNED REPRESENTATIVE CONTINUOUS NUMBER
#PERMANOVA

library(vegan)

var <- Data[,6:30]

perm <- adonis(var ~ Species, data = Data, permutations = 999)
perm

#NMDS
nmds = metaMDS(var, distance = "bray")
nmds
plot(nmds)
data.scores = as.data.frame(scores(nmds))

data.scores$Species = Data$Species
data.scores$Location = Data$Location
head(data.scores)

library(ggplot2)
library(RColorBrewer)

nb.cols <- 18
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

grp.a <- data.scores[data.scores$Species == "sp3", ][chull(data.scores[data.scores$Species == 
                                                                         "sp3", c("NMDS1", "NMDS2")]), ]  
grp.b <- data.scores[data.scores$Species == "sp2", ][chull(data.scores[data.scores$Species == 
                                                                         "sp2", c("NMDS1", "NMDS2")]), ] 
grp.c <- data.scores[data.scores$Species == "xamachana", ][chull(data.scores[data.scores$Species == 
                                                                               "xamachana", c("NMDS1", "NMDS2")]), ]  

hull.data <- rbind(grp.a, grp.b, grp.c) 
hull.data

#Plotting NMDS

xx = ggplot(data.scores, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 4, aes( shape = Location, colour = Species))+ 
  geom_polygon(data=hull.data,aes(x=NMDS1,y=NMDS2,fill=Species,group=Species),alpha=0.30) + 
  theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"), 
        axis.text.x = element_text(colour = "black", face = "bold", size = 12), 
        legend.text = element_text(size = 12, face ="bold", colour ="black"), 
        legend.position = "right", axis.title.y = element_text(face = "bold", size = 14), 
        axis.title.x = element_text(face = "bold", size = 14, colour = "black"), 
        legend.title = element_text(size = 14, colour = "black", face = "bold"), 
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, size = 1.2),
        legend.key=element_blank()) + 
  labs(x = "NMDS1", colour = "Species", y = "NMDS2", shape = "Location")

xx

jit <- position_jitter(seed = 123)

xx + geom_jitter(aes(shape = Location, color = Species), size = 3, position = jit) +
  scale_shape_manual(values = c(0,1,2,3,4,5,6,7,8))

# Cluster analysis setting up data
mydata <- Data[,6:30]
mydata <- mydata[1:95,]
mydata <- na.omit(mydata)
mydata <- scale(mydata)

# Cluster including plot

wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
for (i in 1:10) wss[i] <- sum(kmeans(mydata,
                                     centers=i)$withinss)
plot(1:10, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares")

fit <- kmeans(mydata, 3) 
aggregate(mydata,by=list(fit$cluster),FUN=mean)
mydata <- data.frame(mydata, fit$cluster)

d <- dist(mydata, method = "euclidean") 
fit <- hclust(d, method="ward")
plot(fit, labels= Combined$SpeciesID)
rect.hclust(fit, k=2, border="red")

library(pvclust)
fit <- pvclust(mydata, method.hclust="ward",
               method.dist="euclidean")
plot(fit)
pvrect(fit, alpha=.95)

library(mclust)
fit <- Mclust(mydata)
plot(fit)
summary(fit)
fit <- kmeans(mydata, 2)


#SIMPER ANALYSIS 
library(vegan)
cont.var <- Data[,6:30]
(sim <- with(Data, simper(cont.var, Species)))
summary(sim)

#CLUSTER ANALYSIS WITHIN SPECIES
var <- xamachana[,6:30]

perm <- adonis(var ~ Location, data = xamachana, permutations = 999)
perm

mydata <- xamachana[,6:30]
mydata <- mydata[1:65,]
mydata <- na.omit(mydata)
mydata <- scale(mydata)

d <- dist(mydata, method = "euclidean") 
fit <- hclust(d, method="ward")
plot(fit, labels= xamachana$LocationID)
rect.hclust(fit, k=2, border="red")

var <- sp3[,6:30]

perm <- adonis(var ~ Location, data = sp3, permutations = 999)
perm

mydata <- sp3[,6:30]
mydata <- mydata[1:24,]
mydata <- na.omit(mydata)
mydata <- scale(mydata)


d <- dist(mydata, method = "euclidean") 
fit <- hclust(d, method="ward")
plot(fit, labels= sp3$LocationID)
rect.hclust(fit, k=4, border="red")


#SIMPER ANALYSIS FOR WITHIN SPECIES
library(vegan)
cont.var <- sp3[,6:30]
(sim <- with(sp3, simper(cont.var, Location)))
summary(sim)

cont.var <- xamachana[,6:30]
(sim <- with(xamachana, simper(cont.var, Location)))
summary(sim)