#  ------------------------------------------------------------------------

# Shlukova analyza --------------------------------------------------------

#  ------------------------------------------------------------------------


install.packages(c('vegan', 'gclus', 'cluster', 'vegan', 'ade4'))

library(gclus)
library(cluster)
library(vegan)
library(ade4)


# nastaveni cesty

getwd()
setwd("F:/FSTA_2015/06_cviceni")
dir()


ryby<-read.table('ryby.csv', sep=";", header=TRUE, dec=".")
lokality<-read.table('lokality.csv', sep=";", header=TRUE, dec=".")



# Priklad 1
# Provedte shlukovou analyzu na datech lokality.csv.

# 1) Spocitejte asociacni matici zalozenou na Euklidovske vzdalensoti
# 2) Vyzkousejte ruzne algorytmy
# 3) Rozdelte data do optimalniho poctu shluku

# -> 1)
fix(lokality)

lokality.stan<-decostand(lokality, "standardize")
dist.lokality<-dist(lokality.stan,method='euclidean') # vypocet asociacni matice na zaklade Euklidovske vzdalenosti
hist(as.vector(dist.lokality)) # histogram hodnot asociacnich koeficientu -> zjistime zda v datech mame dobre oddelitelne shluky

# -> 2)
?hclust
dist.lokality.single<-hclust(dist.lokality, method="single") # provede shlukovou analyzu pomoci algoritmu single linkage/metoda nejblisiho souseda 
plot(dist.lokality.single) # vykresli dendrogram, zde je otazka do kolika shluku dana data rozdelit(moznosti jsou 2 x 4 )

dist.lokality.com<-hclust(dist.lokality, method="complete") # provede shlukovou analyzu pomoci algoritmu complete linkage/metoda nejvzdalenejsiho souseda 
plot(dist.lokality.com) # vykresli dendrogram, zde je otazka do kolika shluku dana data rozdelit(moznosti jsou 2 x 4 x 5 x .... )


dist.lokality.ward<-hclust(dist.lokality, method="ward.D") # provede shlukovou analyzu pomoci Wradovy metody -> porovnava soucty ctvercu prirustku vzdalensoti od shlukoveho prumeru
plot(dist.lokality.ward) # vykresli dendrogram, zde je otazka do kolika shluku dana data rozdelit(moznosti jsou 2 x 4 x 5 x .... )


# -> 3)
?cutree # funkce ktera nam vztvori sloupecek prirazeni objektu do shluku
# Algoritmus nejblizsiho souseda
dist.lokality.single.cut<-cutree(dist.lokality.single, k=4) # ziskame vektor prirazeni objektu do shluku
plot (dist.lokality.single)
lokality$single.4<-cutree(dist.lokality.single, k=4) #vlozime si ho k datum, dale jiz muzeme postupovat pomoci kasicke statistiky a pouzivat testy
fix(lokality)

# Algoritmus nejvzdalenejsiho souseda
dist.lokality.com.cut<-cutree(dist.lokality.com, k=5) # ziskame vektor prirazeni objektu do shluku
lokality$complete.5<-cutree(dist.lokality.com, k=5) #vlozime si ho k datum, dale jiz muzeme postupovat pomoci kasicke statistiky a pouzivat testy

# Wardova metoda
dist.lokality.ward.cut<-cutree(dist.lokality.ward, k=5) # ziskame vektor prirazeni objektu do shluku
lokality$ward.5<-cutree(dist.lokality.ward, k=5) #vlozime si ho k datum, dale jiz muzeme postupovat pomoci kasicke statistiky a pouzivat testy


## Vytvoreni hezciho grafu
?reorder.hclust
dist.lokality.single.reor<-reorder.hclust(dist.lokality.single,dist.lokality) # seradi lokality

plot(dist.lokality.single.reor) # vykresli serazeny dendrogram
plot(dist.lokality.single.reor, hang=-1) # uprava dendrogramu
rect.hclust(dist.lokality.single.reor, k=4) # rozdeli dentdrogram do k shluku


# Ukol 1
# Provedte shlukovou analyzu na datovem souboru ryby:
# 1) Zvolte vhodnou miru asociace mezi objekty
# 2) Jsou v datech dobre odelitelne shluky?
# 3) Vyzkousejte metodu nejvzdalenejsiho a nejblizsiho souseda a wardovu metodu. Vykreslete dendrogramy. 
#    Ktera metoda dava nejlepsi vysledek?
# 4) Do kolika shluku byste soubor rozdelili.


# -> 1) jedna se data abundanci, musime tedy  pouzit nektery koeficient podobnosti - > Jaccarduv koefieficient


# -> 2) toto yjistime pomoci histogramu
hist(as.vector(dist.ryby.b))  # histogram neni bimodalni - > nepredpokladame v datech dobre oddelitelne shluky
                              # objekty jsou od sebe daleko, nejsou si moc pobne 

# -> 3)
plot(hclust(dist.ryby.b, method="single"))    # hodne zretezeny dendrogram -> dostali bychom nekolik shluku, ktere by obsahovali pouze po jedne lokalite
plot(hclust(dist.ryby.b, method="complete"))  # jiz mene zretezeny dendrogram (nez predchozi) -> ale opet bychom dostali nekolik shluku, ktere by obsahovali pouze po jedne lokalite
plot(hclust(dist.ryby.b, method="ward.D"))    # POZOR "ward" x "ward.D"
                                              # Nejhezci dendrogram - > zde dostaneme rozumne velke shluky, nedostaneme zde zadny shluk, ktery by obsahoval pouze jeden objekt (pokud nebudeme delit na 8 a vice shluku)

# -> 4)
# Podle Wardovy metody do 3,4,5 nebo 6 shluku
# U ostanich algoritmu je to slozite urcit neb dendrogramy jsou znacne zretezene  

#  ------------------------------------------------------------------------

# Validace shlukove analyzy -----------------------------------------------


# Algoritmus ------------------------------------------------------
#1) Pomoci kofenetickeho indexu - korelce mezi puvodni matici vzdalenosti a kofenetickou matici, ktera predstavuje matici vzdalenosti 
#   objektu v okamziku, kdy byl objekt poprve zarazen do shluku.  

?cophenetic

# Priklad 2 
# Provnejte predchozi agloritmy shlukovani na datech lokality, ktery nejlepe odpovida puvodnim datum?

dist.lokality.single.cop<-cophenetic(dist.lokality.single) # vypocet kofeneticke matice - nejblizsi soused
cor (dist.lokality, dist.lokality.single.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota 
plot (dist.lokality, dist.lokality.single.cop)  #vypocet korelacniho indexu. Neni uvedena p-hondota
#-> nelze interpretovat statistickou vyznamnost -> matice jsou na sobe zavisle

dist.lokality.single.cop<-cophenetic(dist.lokality.single) # vypocet kofeneticke matice - nejblizsi soused
dist.lokality.com.cop<-cophenetic(dist.lokality.com) # vypocet kofeneticke matice - nejvzdalenejsi soused
dist.lokality.ward.cop<-cophenetic(dist.lokality.ward) # vypocet kofeneticke matice - Wardova metoda

cor (dist.lokality, dist.lokality.single.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota
cor (dist.lokality, dist.lokality.com.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota
cor (dist.lokality, dist.lokality.ward.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota

# v porovnani korelacnich indexu nejlepe vychazi metoda nejvzdalenejsi souseda, s nevyssi hodnotou korelacniho indexu 

# Ukol 2
# Urcete nejlepsi argoritmus pro shlukovou analyzu provedenou na datovem souboru ryby

dist.ryby.single.cop<-cophenetic(hclust(dist.ryby.b, method="single")) # vypocet kofeneticke matice - nejblizsi soused
dist.ryby.com.cop<-cophenetic(hclust(dist.ryby.b, method="complete")) # vypocet kofeneticke matice - nejvzdalenejsi soused
dist.ryby.ward.cop<-cophenetic(hclust(dist.ryby.b, method="ward.D")) # vypocet kofeneticke matice - Wardova metoda

cor (dist.ryby.b, dist.lokality.single.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota
cor (dist.ryby.b, dist.lokality.com.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota
cor (dist.ryby.b, dist.lokality.ward.cop, method="spearman")  #vypocet korelacniho indexu. Neni uvedena p-hondota

# v porovnani korelacnich indexu nejlepe vychazi Wardova metoda

# Optimalni pocet shluku --------------------------------------------------

####################
# 1) Metoda siluety
# Jedna se o metodu, ktera pocita sirku siluety pro jednotlive objekty, prumernou sirku siluety pro shluk a prumernou sirku siluety
# pro cely soubor   

?silhouette
silueta<-silhouette(cutree(dist.lokality.com, k=4), dist.lokality) # ukazka vypoctu siluety pro 4 shluky 
summary(silueta)$avg.width


# Priklad 3
# Zjistet optimalni pocet shluku pro datovzsoubor lokality pro metodu nejvzdalenejsiho souseda pro vsechny mozne pocty shluku 

##vytvorit prazdny vektor
asw<-numeric(nrow(lokality)) 
asw

##naplnit vektor hodontami "siluety"
for (k in 2:(nrow(lokality)-1)){
  sil <-silhouette(cutree(dist.lokality.com, k=k), dist.lokality)
  asw[k]<-summary(sil)$avg.width
}
asw
k.best<-which.max(asw)
k.best

# Optimalni pocet shluku pro datovy soubor je 2 shluky

# Graficke zobrazeni vysledku siluety
plot(1:nrow(lokality), asw, type="h", main="Silhouette-optimal number of cluster",
     xlab="k (number of groups)", ylab="Avraege silhouette width")
axis(1, k.best, paste ("optimum", k.best, sep="\n"), col="red", font=2, col.axis="red")
points(k.best, max(asw), pch=16, col="red", cex=1.5)

# Finalni vykresleni dendrogramu roydeleneho do optimalniho poctu shluku 
plot(dist.lokality.com, hang=-1) # uprava dendrogramu
rect.hclust(dist.lokality.com, k=2) # rozdeli dentdrogram do k shluku

# Ukol 3 
# Zjiste optimalni pocet shluku pro datovy soubor ryby, pro nejlepsi algoritmus z predchoziho ukolu - pomoci metody siluety 

asw<-numeric(nrow(ryby)) ##vytvorit prazdny vektor

##naplnit vektor hodontami "siluety"
for (k in 2:(nrow(ryby)-1)){
  sil <-silhouette(cutree(hclust(dist.ryby.b, method="ward.D"), k=k), dist.ryby.b)
  asw[k]<-summary(sil)$avg.width
}

k.best<-which.max(asw)
k.best  # Optimalni pocet shluku pro datovy soubor je 13 shluku

# Graficke zobrazeni vysledku siluety
plot(1:nrow(ryby), asw, type="h", main="Silhouette-optimal number of cluster",
     xlab="k (number of groups)", ylab="Avraege silhouette width")
axis(1, k.best, paste ("optimum", k.best, sep="\n"), col="red", font=2, col.axis="red")
points(k.best, max(asw), pch=16, col="red", cex=1.5)

# Finalni vykresleni dendrogramu roydeleneho do optimalniho poctu shluku 
plot(hclust(dist.ryby.b, method="ward.D"), hang=-1) # uprava dendrogramu
rect.hclust(hclust(dist.ryby.b, method="ward.D"), k=k.best) # rozdeli dentdrogram do k shluku

####################
# 1) Mantel test

# Priklad 4
# Zjistet optimalni pocet shluku pro datovy soubor lokality pro metodu nejblizsiho souseda pro vsechny mozne pocty shluku 

###tvorba funkce 

grpdist<- function(X)
{
  require(cluster)
  gr<-as.data.frame(as.factor(X))
  distgr<-daisy(gr, "gower")
  distgr
}

kt<-data.frame(k=1:nrow(lokality), r=0)
kt

for (i in 2: (nrow(lokality)-1)){ 
  gr<-cutree(dist.lokality.com,i)
  distgr<-grpdist(gr)
  mt<-cor(dist.lokality, distgr, method="spearman")
  kt[i,2]<-mt
}

names(kt)
kt

k.best<-which.max(kt$r)
k.best


# Optimalni pocet shluku pro datovy soubor je 5 shluku

# Graficke zobrazeni vysledku Mantel testu
plot( kt$k, kt$r, type="h", main="Mantel-optimal number of cluster",
      xlab="k (number of groups)", ylab="Sepearman correlation")
axis(1, k.best, paste ("optimum", k.best, sep="\n"), col="red", font=2, col.axis="red")
points(k.best, max(kt$r[-1]), pch=16, col="red", cex=1.5)

# Ukol 4 
# Zjiste optimalni pocet shluku pro datovy soubor ryby, pro nejlepsi algoritmus z predchoziho ukolu - pomoci Mantel testu


kt<-data.frame(k=1:nrow(ryby), r=0)
for (i in 2: (nrow(ryby)-1)){ 
  gr<-cutree(hclust(dist.ryby.b, method="ward.D"),i)
  distgr<-grpdist(gr)
  mt<-cor(dist.ryby.b, distgr, method="spearman")
  kt[i,2]<-mt
}

k.best<-which.max(kt$r)
k.best

# Optimalni pocet shluku pro datovy soubor je 3 shluku

# Graficke zobrazeni vysledku Mantel testu
plot( kt$k, kt$r, type="h", main="Mantel-optimal number of cluster",
      xlab="k (number of groups)", ylab="Sepearman correlation")
axis(1, k.best, paste ("optimum", k.best, sep="\n"), col="red", font=2, col.axis="red")
points(k.best, max(kt$r[-1]), pch=16, col="red", cex=1.5)

# Finalni vykresleni dendrogramu roydeleneho do optimalniho poctu shluku 
plot(hclust(dist.ryby.b, method="ward.D"), hang=-1) # uprava dendrogramu
rect.hclust(hclust(dist.ryby.b, method="ward.D"), k=k.best) # rozdeli dentdrogram do k shluku

##### Porovnani vystupu ze siluety a Mantel testu
par(mfrow=c(1,2))
## 13 shluku -> silueta
plot(hclust(dist.ryby.b, method="ward.D"), hang=-1)
rect.hclust(hclust(dist.ryby.b, method="ward.D"), k=13)

## 3 shluky -> Mantel test
plot(hclust(dist.ryby.b, method="ward.D"), hang=-1)
rect.hclust(hclust(dist.ryby.b, method="ward.D"), k=3) 

# Jako lepsi vysledek se zde jevi rozdeleni do 3 shluku, nedostaneme tak shluky, ktere by oobsahovaly po jednom objektu
# tak jako v pripade shluku 13


# Nehiearchicke shlukovani ------------------------------------------------

# Priklad 5 
# Provedte shlukovou analyzu na zaklade K-means clustering

?kmeans
lokality.kmeans<-kmeans(dist.lokality, centers=5, nstart=100) # vypocet shlukove analyzy pomoci algoritmu K-means5
table(lokality.kmeans$cluster, cutree(dist.lokality.com, k=5)) # Srovnani vysledku alforitmu K-means a nejbliysiho souseda

lokality.km.cascade<-cascadeKM(dist.lokality, inf.gr=2, sup.gr=10, iter=100, criterion="ssi")
# Vypocet K means shlukovani pro pocet shluku 2 az 10, doplnene o hodntu ssi  - Simple structure index - ktera ukazuje 
# na nejvhodnejsi pocet shluku (podle nejvetsi hodnoty ssi)  
plot(lokality.km.cascade, sortg=TRUE) # Graficke vykreseni predchoziho prikazu