#  ------------------------------------------------------------------------
# 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()
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)
head(USArrests)

# maticovy graf

panel.hist <- function(x, ...) ## funkce pro vykreslení maticového 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[,1])						# 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 nejblizsiho souseda
dist.USA.single <- hclust(dist.USArrests.stan, method="single") 
plot(dist.USA.single, hang=-1)	    # 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
	main = "Cluster Dendogram - single linkage")
	  
# shlukova analyza pomoci algoritmu complete linkage = metoda nejvzdalenejsiho souseda 
dist.USA.com <- hclust(dist.USArrests.stan, method="complete") 
plot(dist.USA.com, hang=-1)		
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,
	main = "Cluster Dendogram - complete linkage")

# 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, hang=-1)			
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,
	main = "Cluster Dendogram - Ward´s method")

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


# pokud by nás zajímalo, které skupiny se od sebe liší, použijeme následující příkazy.
install.packages("FSA")           # nacteni balicku pro neparametrický párový test
library(FSA)
dunnTest(Murder~ward.4,data=USArrests, method="bh")
dunnTest(Assault~ward.4,data=USArrests, method="bh")
dunnTest(UrbanPop~ward.4,data=USArrests, method="bh")
dunnTest(Rape~ward.4,data=USArrests, method="bh")



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

# pro připomenutí: dist.USA.single <- hclust(dist.USArrests.stan, method="single")
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 (complete linkage) podava nejlepsi vysledky

# 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 
silueta
summary(silueta)$avg.width
fviz_silhouette(silueta) # vyobrazeni 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 hodnotami "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)						

	
					
####################
# 2) Mantel test

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

###tvorba funkce pro výpočet asociační matice založené na Gowerově indexu

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)