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

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

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

# nastaveni cesty

getwd()
setwd("c:/Users/brozova/Desktop/vyuka/vicerozmerky_podzim_2016_bimat/2. cviceni/")
dir()

lokality<-read.table('lokality.csv', sep=";", header=TRUE, dec=".")
row.names(lokality)<-lokality[,1]
lokality<-lokality[,-1]
head(lokality,10)

ryby<-read.table('ryby.csv', sep=";", header=TRUE, dec=".")
row.names(ryby)<-ryby[,1]
ryby<-ryby[,-1]
head(ryby,10)

# maticovy graf

## put histograms on the diagonal
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
pairs(lokality[1:5],data=lokality, diag.panel = panel.hist, main="Simple Scatterplot Matrix")

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

# 1) Spocitejte asociacni matici zalozenou na Euklidovske vzdalenosti
# 2) Vyzkousejte ruzne algoritmy shlukovani
# 3) Rozdelte data do optimalniho poctu shluku

# -> 1)

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

# -> 2)

?hclust
dist.lokality.single<-hclust(dist.lokality.stan, method="single") # provede shlukovou analyzu pomoci algoritmu single linkage = metoda nejblisiho souseda 
plot(dist.lokality.single) # vykresli dendrogram

dist.lokality.com<-hclust(dist.lokality.stan, method="complete") # provede shlukovou analyzu pomoci algoritmu complete linkage = metoda nejvzdalenejsiho souseda 
plot(dist.lokality.com) # vykresli dendrogram

dist.lokality.ward<-hclust(dist.lokality.stan, method="ward.D2") # provede shlukovou analyzu pomoci Wardovy metody -> porovnava soucty ctvercu prirustku vzdalenosti od shlukoveho prumeru
plot(dist.lokality.ward) # vykresli dendrogram

# -> 3)

?cutree # funkce, ktera nam vytvori sloupecek prirazeni objektu do shluku
# k = an integer scalar or vector with the desired number of groups
# h = numeric scalar or vector with heights where the tree should be cut.

install.packages("RColorBrewer") # balicek nabizi velke mnozstvi palet barev, na http://colorbrewer2.org/ moznost nahledu
library("RColorBrewer")

barvy<-brewer.pal(6,"Spectral")
display.brewer.pal(6,"Spectral")

# Algoritmus nejblizsiho souseda
plot(dist.lokality.single) # vykresli dendrogram
lokality$single.4<-cutree(dist.lokality.single, k=4) #pridame vektor prirazeni objektu do shluku 
myplclust(dist.lokality.single,labels=row.names(lokality),lab.col=barvy[lokality$single.4]) #vykresli dendrogram s barevnym odlisenim shluku

# Algoritmus nejvzdalenejsiho souseda
plot(dist.lokality.com) # vykresli dendrogram
lokality$com.5<-cutree(dist.lokality.com, k=5) #pridame vektor prirazeni objektu do shluku 
lokality$com.5<-cutree(dist.lokality.com, h=5) #na vzdalenosti 5 provede rez - zde ekvivaletni prikazu vyse 
abline(h=5,col="red") # jak bude provedena klasifikace do shluku
myplclust(dist.lokality.com,labels=row.names(lokality),lab.col=barvy[lokality$com.5]) #vykresli dendrogram s barevnym odlisenim shluku

# Wardova metoda
plot(dist.lokality.ward) # vykresli dendrogram
lokality$ward.5<-cutree(dist.lokality.ward, k=5) #pridame vektor prirazeni objektu do shluku 
myplclust(dist.lokality.ward,labels=row.names(lokality),lab.col=barvy[lokality$ward.5]) #vykresli dendrogram s barevnym odlisenim shluku

## Vytvoreni hezciho grafu
?reorder.hclust
dist.lokality.single.reor<-reorder.hclust(dist.lokality.single,dist.lokality.stan) # seradi lokality ve shlucich podle vzdalenosti

plot(dist.lokality.single.reor) # vykresli serazeny dendrogram
plot(dist.lokality.single.reor, hang=-1) # uprava dendrogramu, funguje i u funkce plot
rect.hclust(dist.lokality.single.reor, k=4) # rozdeli dendrogram do k shluku, funguje i s parametrem h

# Samostatny ukol 1
# 1) priradte ceske ekvivalenty k nazvum algoritmu, ktere nabizi funkce hclust

?hclust 
# ward.D/ward.D2 = 
# single = 
# complete = 
# average(= UPGMA) = 
# mcquitty(= WPGMA) = 
# median(= WPGMC) = 
# centroid(= UPGMC) = 

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



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

# Algoritmus ------------------------------------------------------
#1) Pomoci kofenetickeho indexu - korelace 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.stan, dist.lokality.single.cop, method="spearman")  #vypocet korelacniho koeficientu
plot (dist.lokality.stan, dist.lokality.single.cop)  #z grafu lze videt, ze hodnoty kofeneticke matice silne koreluji s puvodnimi vzdalenostmi

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

#vypocet korelacniho koeficientu pro vsechny algoritmy
cor (dist.lokality.stan, dist.lokality.single.cop, method="spearman")  
cor (dist.lokality.stan, dist.lokality.com.cop, method="spearman")  
cor (dist.lokality.stan, dist.lokality.ward.cop, method="spearman")

# v porovnani korelacnich koeficientu nejlepe vychazi metoda nejvzdalenejsiho souseda (s nevyssi hodnotou korelacniho koeficientu) 

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


# 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.stan) # ukazka vypoctu siluety pro 4 shluky 
summary(silueta)$avg.width


# Priklad 3
# Zjistete optimalni pocet shluku pro datovy soubor 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.stan)
  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 clusters)", ylab="Average 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 rozdeleneho do optimalniho poctu shluku 
plot(dist.lokality.com, hang=-1) # uprava dendrogramu
rect.hclust(dist.lokality.com, k=2) # rozdeli dentdrogram do k shluku

# Samostatny ukol 3 
# Zjiste optimalni pocet shluku pro datovy soubor ryby pro algoritmus nejvzdalenejsiho souseda pomoci metody siluety.


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

# Priklad 4
# Zjistete optimalni pocet shluku pro datovy soubor lokality pro metodu nejvzdalenejsiho souseda pro vsechny mozne pocty shluku 

###tvorba funkce 

grpdist<- function(X)
{
  require(cluster)
  gr<-as.data.frame(as.factor(X))
  distgr<-daisy(gr, "gower") # 0 pokud objekty nejsou spolu ve shluku, 1 pokud jsou spolu ve shluku - R prevadi na vzdalenosti (expressed as a dissimilarity)
  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.stan,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 clusters)", ylab="Spearman 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)

##### Porovnani vystupu ze siluety a Mantel testu
par(mfrow=c(1,2))
## 2 shluky -> silueta
plot(hclust(dist.lokality.stan, method="complete"), hang=-1)
rect.hclust(hclust(dist.lokality.stan, method="complete"), k=2)

## 5 shluku -> Mantel test
plot(hclust(dist.lokality.stan, method="complete"), hang=-1)
rect.hclust(hclust(dist.lokality.stan, method="complete"), k=5)

# Jako lepsi vysledek se zde jevi rozdeleni do 5 shluku

# Samostatny ukol 4 
# Zjiste optimalni pocet shluku pro datovy soubor ryby, pro algoritmus nejvzdalenejsiho souseda pomoci Mantelova testu.
# Porovnejte rozdeleni dat podle Mantelova testu a indexu siluety.  

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

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

?kmeans # partition the points into k groups such that the sum of squares from points to the assigned cluster centres is minimized.
lokality.kmeans<-kmeans(dist.lokality.stan, centers=5, nstart=100) # vypocet shlukove analyzy pomoci algoritmu K-means
table(lokality.kmeans$cluster, cutree(dist.lokality.com, k=5)) # Srovnani vysledku alforitmu K-means a nejvzdalenejsiho souseda

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