#  ------------------------------------------------------------------------
# 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/brozoval/Desktop/3. cviceni/")
dir() # vypise vsechny soubory v nastavenem adresari

# 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 na diagonale maticoveho grafu
{
    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

# -> 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)		
fviz_dend(dist.USA.com,			
	k = 4,					
	cex = 0.8,					
	k_colors = brewer.pal(4,"Set1"),
	color_labels_by_k = TRUE, 	
	rect = TRUE,			
	horiz=T)				

# 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)			
fviz_dend(dist.USA.ward,				
	k = 4,			
	cex = 0.8,		
	k_colors = brewer.pal(4,"Set1"),
	color_labels_by_k = TRUE, 
	rect = TRUE,
	horiz=T)	

# -> 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)
USArrests[order(USArrests$ward.4),]

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

# v kterych parametrech se skupiny statu lisi?

aggregate(USArrests$Murder ~ USArrests$ward.4, FUN = function(x) c(Median=median(x)))
aggregate(USArrests$Assault ~ USArrests$ward.4, FUN = function(x) c(Median=median(x)))
aggregate(USArrests$UrbanPop ~ USArrests$ward.4, FUN = function(x) c(Median=median(x)))
aggregate(USArrests$Rape ~ USArrests$ward.4, FUN = function(x) c(Median=median(x)))

kruskal.test(Murder~ward.4,data=USArrests)
kruskal.test(Assault~ward.4,data=USArrests)
kruskal.test(UrbanPop~ward.4,data=USArrests)
kruskal.test(Rape~ward.4,data=USArrests)

# Samostatny ukol:
# Provedte shlukovou analyzu na datovem souboru mtcars:
# 1) Zvolte vhodnou miru asociace mezi objekty, spocitejte asociacni matici.
# 2) Vyzkousejte algritmus nejvzdalenejsiho souseda a Wardovu metodu. Vykreslete dendrogramy. 
#    Ktera metoda dava nejlepsi vysledek? 
# 3) Rozdelte modely aut do 4 shluku na zaklade Wardovy metody. Uvedte pocet modelu v jednotlivych shlucich.

#nacteni dat
data("mtcars")
?mtcars # popis dat
dim(mtcars)
head(mtcars)

# -> 1) 

summary(mtcars)						    # promenne maji odlisny rozsah, data je potreba standardizovat
mtcars.stan <- mtcars%>%
  na.omit() %>%          			# odstrani chybejici hodnoty
  scale()								      # standardizuje promenne
summary(mtcars.stan)					# prumer promennych je roven 0

dist.mtcars.stan<-dist(mtcars.stan,method='euclidean')# vypocet asociacni matice na standardizovanych datech

# -> 2)

# shlukova analyza pomoci algoritmu complete linkage = metoda nejvzdalenejsiho souseda
dist.mtcars.complete <- hclust(dist.mtcars.stan, method="complete") 
plot(dist.mtcars.complete)

# shlukova analyza pomoci algoritmu single linkage = metoda nejblizsiho souseda
dist.mtcars.ward<- hclust(dist.mtcars.stan, method="ward.D2") 
plot(dist.mtcars.ward)

# -> 3)
			
fviz_dend(dist.mtcars.ward,
	main="Wardova metoda",
	k = 4,
	cex = 0.8,
	k_colors = brewer.pal(4,"BrBG"),
	color_labels_by_k = TRUE, 
	rect = TRUE,		
	horiz=T)

mtcars$ward.4<-cutree(dist.mtcars.ward, k=4)
table(mtcars$ward.4)

#  ------------------------------------------------------------------------
#  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 statu USA pomoci korelace kofeneticke matice s matici asociacni.

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 mtcars

dist.mtcars.complete.cop<-cophenetic(dist.mtcars.complete) 
dist.mtcars.ward.cop<-cophenetic(dist.mtcars.ward) 

cor (dist.mtcars.stan, dist.mtcars.complete.cop, method="spearman") 
cor (dist.mtcars.stan, dist.mtcars.ward.cop, method="spearman") # vychazi temer totozne

# 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=4), 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 mtcars pro algoritmus nejvzdalenejsiho souseda pomoci Wardovy metody.

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

##naplnit vektor hodontami "siluety"
for (k in 2:(nrow(mtcars)-1)){
  sil <-silhouette(cutree(dist.mtcars.ward, k=k), dist.mtcars.stan)
  asw[k]<-summary(sil)$avg.width
}
k.best<-which.max(asw)
k.best  # Optimalni pocet shluku pro datovy soubor je 10 shluku

# Graficke zobrazeni vysledku siluety
plot(1:nrow(mtcars), 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(dist.mtcars.ward, k=10), dist.mtcars.stan) # vypocet siluety pro 10 shluku 
summary(silueta)$avg.width
fviz_silhouette(silueta) # obrazeni siluety pro 10 shluku

# Finalni vykresleni dendrogramu rozdeleneho do optimalniho poctu shluku 
fviz_dend(dist.mtcars.ward,			
	main="Wardova metoda",
	k = 10,						# 10 skupin 
	cex = 0.8,						
	k_colors = brewer.pal(10,"BrBG"),		
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)	
					
####################
# 2) 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 mtcars, pro Wardovu metodu, pomoci Mantelova testu.
# Porovnejte rozdeleni dat podle Mantelova testu a indexu siluety.  

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

k.best<-which.max(kt$r)
k.best
# Optimalni pocet shluku pro datovy soubor je 5 

# 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.mtcars.ward,			
	main="Wardova metoda",
	k = 5,				
	cex = 0.8,						
	k_colors = brewer.pal(6,"BrBG")[1:5],		
	color_labels_by_k = TRUE, 			
	rect = TRUE,					
	horiz=T)