###HIERARCHICKA KLAS
library(vegan)

?dist
?vegdist

head(forest)

spe <- read.delim("forest_JMK_species.txt", header=T, 
                  row.names = 1)
spe[1:5,1:5]

spe.pa <- decostand(spe, method="pa")
spe.pa[1:5,1:5]

dist.jac <- vegdist(spe.pa, method = "jaccard")
class(dist.jac)
dist.jac

dist.jac.m <- as.matrix(dist.jac)
dist.jac.m[1:5,1:5]

dist.bray <- vegdist(spe, method = "bray")

plot(dist.jac, dist.bray)

dist.geo <- vegdist(forest[, c("X", "Y")], method = "euclidean")
dist.geo.m <- as.matrix(dist.geo)
dist.geo.m[1:5,1:5]

plot(dist.geo, dist.jac)
plot(dist.geo, dist.bray)

dist.mantel <- mantel(dist.geo, dist.bray, method="pearson", 
                      permutations=99)
dist.mantel

?hclust

clust.single <- hclust(dist.jac, method = "single")
plot(clust.single, hang=-1)

clust.comp <- hclust(dist.jac, method = "complete")
plot(clust.comp, hang=-1)

clust.avg <- hclust(dist.jac, method = "average")
plot(clust.avg, hang=-1)

clust.avg.k5 <- cutree(clust.avg, k=5)

plot(forest[, c("X", "Y")], pch=21, 
     bg=my.col[clust.avg.k5],
     main= "Jaccard + UPGMA", cex=1.2)

boxplot(ext$TWI ~ clust.avg.k5, 
        col=my.col[unique(clust.avg.k5)])

library(vioplot)
vioplot(ext$Altitude ~ clust.avg.k5, 
        col=my.col[unique(clust.avg.k5)])


clust.avg.ab <- hclust(dist.bray, method = "average")
clust.avg.ab.k5 <- cutree(clust.avg.ab, k=5)

plot(forest[, c("X", "Y")], pch=21, 
     bg=my.col[clust.avg.ab.k5],
     main= "Bray + UPGMA", cex=1.2)
table(clust.avg.ab.k5)

library(ade4)
is.euclid(dist.jac)
is.euclid(sqrt(dist.jac))

dist.jac.c <- cailliez(dist.jac, cor.zero = F)
head(dist.jac.c)
head(dist.jac.c - dist.jac)

clust.ward <- hclust(dist.jac.c, method = "ward.D2")
plot(clust.ward)

rect.hclust(clust.ward, k=5, border="blue")
clust.ward.k5 <- cutree(clust.ward, k=5)

plot(forest[, c("X", "Y")], pch=21, 
     bg=my.col[clust.ward.k5],
     main= "Ward + UPGMA", cex=1.2)
table(clust.ward.k5)

plot(env$Altitude, axes=F, col=terrain.colors(100))
points(forest[, c("X", "Y")], pch=21, 
       bg=my.col[clust.ward.k5])
title("Lesy v JMK")

legend("topleft", legend = 1:5, pch=21, 
       pt.bg=my.col[unique(clust.ward.k5)], bty="n",
       pt.cex = 1.5)

boxplot(forest$Diversity ~ clust.ward.k5)
boxplot(rowSums(spe.pa) ~ clust.ward.k5)
