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

# instalace balicku

install.packages("cluster") # pro vypocet siluety
install.packages("magrittr") # pro funkci %>%
install.packages("ggplot2") # pro vykresleni hezcich dendrogramu
install.packages("factoextra") # pro vykresleni hezcich dendrogramu
install.packages("RColorBrewer") # balicek nabizi velke mnozstvi palet barev, na http://colorbrewer2.org/ moznost nahledu
install.packages("vegan") # pro vypocet asociacni matice

library(cluster)
library(magrittr)
library(ggplot2)
library(factoextra)
library(RColorBrewer)
library(vegan)

# nastaveni cesty

getwd()
setwd("C:/Users/rcx2ucitel/Desktop/2. cviceni")
dir()

# nacteni dat
data("USArrests")
dim(USArrests) # pocty zatcenych v kazdem ze statu USA
names(USArrests)
# "Murder"  = pocet zatcenych za vrazdu (na 100 000 obyvatel)
# "Assault" = pocet zatcenych za napadeni (na 100 000 obyvatel)
# "UrbanPop" = procento obyvatel zijicich ve meste
# "Rape" = pocet zatcenych za znasilneni (na 100 000 obyvatel)

# maticovy graf

panel.hist <- function(x, ...) ## funkce pro vykresleni histogramu
{
    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(USArrests, panel = panel.smooth, diag.panel = panel.hist, main="Simple Scatterplot Matrix")

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

# -> 1)

########## standardizace
summary(USArrests)						# promenne maji odlisny rozsah, data je potreba standardizovat
USArrests.stan <- USArrests%>%
  na.omit() %>%          					# odstrani chybejici hodnoty
  scale()								# standardizuje promenne
summary(USArrests.stan)					 	# prumer promennych je roven 0
var(USArrests.stan[,2])					  	# rozptyl promennych je roven 1
dim(USArrests.stan)						# zadne chubejici hodnoty = dimenze datoveho souboru zustava stejna

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

# -> 2)

?hclust
# shlukova analyza pomoci algoritmu single linkage = metoda nejblisiho souseda
dist.USA.single <- hclust(dist.USArrests.stan, method="single") 
plot(dist.USA.single)					# vykresli dendrogram
fviz_dend(dist.USA.single,				# v dendrogramu je vykreslen vystup funkce hclust
	k = 4,						# pomoci k definujeme pocet skupin, nebo pomoci h vzdalenost, kde maji byt shluky vytvoreny
	cex = 0.8,						# velikost popisku
	k_colors = brewer.pal(4,"Set1"),		# barvy pro jednotlive skupiny
	color_labels_by_k = TRUE, 			# obarveni popisku
	rect = TRUE,					# ohraniceni shluku
	horiz=T)						# vykresli dendrogram horizontalne

# shlukova analyza pomoci algoritmu complete linkage = metoda nejvzdalenejsiho souseda 
dist.USA.com <- hclust(dist.USArrests.stan, method="complete") 
plot(dist.USA.com)					# vykresli dendrogram
fviz_dend(dist.USA.com,					# v dendrogramu je vykreslen vystup funkce hclust
	k = 4,						# pomoci k definujeme pocet skupin, nebo pomoci h vzdalenost, kde maji byt shluky vytvoreny
	cex = 0.8,						# velikost popisku
	k_colors = brewer.pal(4,"Set1"),		# barvy pro jednotlive skupiny
	color_labels_by_k = TRUE, 			# obarveni popisku
	rect = TRUE,					# ohraniceni shluku
	horiz=T)						# vykresli dendrogram horizontalne

# shlukova analyza pomoci Wardovy metody -> porovnava soucty ctvercu prirustku vzdalenosti od shlukoveho prumeru
dist.USA.ward <- hclust(dist.USArrests.stan, method="ward.D2") 
plot(dist.USA.ward)					# vykresli dendrogram
fviz_dend(dist.USA.ward,				# v dendrogramu je vykreslen vystup funkce hclust
	k = 4,						# pomoci k definujeme pocet skupin, nebo pomoci h vzdalenost, kde maji byt shluky vytvoreny
	cex = 0.8,						# velikost popisku
	k_colors = brewer.pal(4,"Set1"),		# barvy pro jednotlive skupiny
	color_labels_by_k = TRUE, 			# obarveni popisku
	rect = TRUE,					# ohraniceni shluku
	horiz=T)						# vykresli dendrogram horizontalne

# -> 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.

USArrests$ward.4<-cutree(dist.USA.ward, k=4) #pridame vektor prirazeni objektu do shluku 
head(USArrests)
table(USArrests$ward.4)

windows()
pairs(USArrests[,1:4], panel = panel.smooth, main="Simple Scatterplot Matrix",col=brewer.pal(4,"Set1")[USArrests$ward.4])

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

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

summary(ryby)
ryby.mis <- ryby%>% na.omit()       			# odstrani chybejici hodnot

# -> 1) jedna se data abundanci, musime tedy  pouzit nektery koeficient podobnosti - > Sorensenův asymetrický koeficient / Jaccarduv koefieficient 
dist.ryby.b<-vegdist(ryby.mis, "jac", binary=T)		# vypocet asociacni matice na ne/pritomnosti druhu (binarnich prom.) = 1 - Jaccard 

# -> 2) Vykreslime histogram asociacni matice
hist(as.vector(dist.ryby.b))  # histogram neni bimodalni - > nepredpokladame v datech dobre oddelitelne shluky

# -> 3)

# shlukova analyza pomoci algoritmu single linkage = metoda nejblizsiho souseda
dist.ryby.single <- hclust(dist.ryby.b, method="single") 
windows()
plot(dist.ryby.single)				# vykresli dendrogram
fviz_dend(dist.ryby.single,			# v dendrogramu je vykreslen vystup funkce hclust
	main="Algoritmus nejblizsiho souseda",	# specifikace algoritmu
	k = 3,						# pomoci k definujeme pocet skupin, nebo pomoci h vzdalenost, kde maji byt shluky vytvoreny
	cex = 0.8,						# velikost popisku
	k_colors = brewer.pal(3,"Set1")[1:3],		# barvy pro jednotlive skupiny
	color_labels_by_k = TRUE, 			# obarveni popisku
	rect = TRUE,					# ohraniceni shluku
	horiz=T)					# vykresli dendrogram horizontalne

# shlukova analyza pomoci algoritmu complete linkage = metoda nejvzdalenejsiho souseda 
dist.ryby.com <- hclust(dist.ryby.b, method="complete") 
plot(dist.ryby.com)					# vykresli dendrogram
fviz_dend(dist.ryby.com,				# v dendrogramu je vykreslen vystup funkce hclust
	main="Algoritmus nejvzdalenejsiho souseda",# specifikace algoritmu
	k = 3,						# pomoci k definujeme pocet skupin, nebo pomoci h vzdalenost, kde maji byt shluky vytvoreny
	cex = 0.8,						# velikost popisku
	k_colors = brewer.pal(3,"Set1"),		# barvy pro jednotlive skupiny
	color_labels_by_k = TRUE, 			# obarveni popisku
	rect = TRUE,					# ohraniceni shluku
	horiz=T)						# vykresli dendrogram horizontalne

#  ------------------------------------------------------------------------
#  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.USA.single.cop<-cophenetic(dist.USA.single) # vypocet kofeneticke matice
coe.single<-cor (dist.USArrests.stan, dist.USA.single.cop, method="spearman") # korelacni koeficient
dist.USA.complete.cop<-cophenetic(dist.USA.com) 
coe.complete<-cor (dist.USArrests.stan, dist.USA.complete.cop, method="spearman") 
dist.USA.ward.cop<-cophenetic(dist.USA.ward)
coe.ward<-cor (dist.USArrests.stan, dist.USA.ward.cop, method="spearman")
coe.single;coe.complete;coe.ward # algoritmus nejvzadelenejsiho souseda podava nejlepsi vysledky

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

dist.ryby.com.cop<-cophenetic(dist.ryby.com) 
dist.ryby.ward.cop<-cophenetic(dist.ryby.ward) 

cor (dist.ryby.b, dist.ryby.com.cop, method="spearman") # algoritmus nejvzadelenejsiho souseda podava nejlepsi vysledky
cor (dist.ryby.b, dist.ryby.ward.cop, method="spearman") 

# 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   

silueta<-silhouette(cutree(dist.USA.com, k=2), dist.USArrests.stan) # vypocet siluety pro 4 shluky 
summary(silueta)$avg.width
fviz_silhouette(silueta) # obrazeni siluety pro 4 shluky

# Priklad 3
# Zjistete optimalni pocet shluku pro datovy soubor USArrests pro metodu nejvzdalenejsiho souseda pro vsechny mozne pocty shluku 

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

##naplnit vektor hodontami "siluety"
for (k in 2:(nrow(USArrests)-1)){
  sil <-silhouette(cutree(dist.USA.com, k=k), dist.USArrests.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(USArrests), 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 
fviz_dend(dist.USA.com,					
	k = 2,						# staty budou rezdeleny do 2 skupin
	cex = 0.8,						
	k_colors = brewer.pal(4,"Set1")[1:2],	
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)						

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

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

##naplnit vektor hodontami "siluety"
for (k in 2:(nrow(ryby.mis)-1)){
  sil <-silhouette(cutree(hclust(dist.ryby.b, method="complete"), 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.mis), asw, type="h", main="Silhouette-optimal number of cluster",
     xlab="k (number of groups)", 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)

silueta<-silhouette(cutree(hclust(dist.ryby.b, method="complete"), k=13), dist.ryby.b) # vypocet siluety pro 4 shluky 
summary(silueta)$avg.width
fviz_silhouette(silueta) # obrazeni siluety pro 4 shluky

# Finalni vykresleni dendrogramu rozdeleneho do optimalniho poctu shluku 
fviz_dend(dist.ryby.com,				# v dendrogramu je vykreslen vystup funkce hclust
	main="Algoritmus nejvzdalenejsiho souseda",# specifikace algoritmu
	k = 13,						# 13 skupin 
	cex = 0.8,						
	k_colors = brewer.pal(13,"Set1"),		
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)	
					
####################
# 1) Mantel test

# Priklad 4
# Zjistete optimalni pocet shluku pro datovy soubor USArrests 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(USArrests), r=0)
kt

for (i in 2: (nrow(USArrests)-1)){ 
  gr<-cutree(dist.USA.com,i)
  distgr<-grpdist(gr)
  mt<-cor(dist.USArrests.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 jsou 3 shluky

# 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)

# Finalni vykresleni dendrogramu rozdeleneho do optimalniho poctu shluku
fviz_dend(dist.USA.com,					
	k = 3,						# staty budou rezdeleny do 3 skupin
	cex = 0.8,						
	k_colors = brewer.pal(4,"Set1")[1:3],	
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)						

# 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.  

kt<-data.frame(k=1:nrow(ryby.mis), r=0)
for (i in 2: (nrow(ryby.mis)-1)){ 
  gr<-cutree(hclust(dist.ryby.b, method="complete"),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 4 

# 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)

# Finalni vykresleni dendrogramu rozdeleneho do optimalniho poctu shluku 
fviz_dend(dist.ryby.com,				# v dendrogramu je vykreslen vystup funkce hclust
	main="Algoritmus nejvzdalenejsiho souseda",# specifikace algoritmu
	k = 4,						# 4 shluky
	cex = 0.8,						
	k_colors = brewer.pal(4,"Set1"),		
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)