#  ------------------------------------------------------------------------
#  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. cvičení/")
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)
summary(USArrests)

################################
######## maticovy graf #########
################################

# histogram na diagonale
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")
plot(USArrests) # nebo jednoduse pomoci funkce plot

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

# -> 1)

summary(USArrests)						    # promenne maji odlisny rozsah -> data je potreba standardizovat/normalizovat
USArrests.norm<-apply(USArrests, 2, function(x) (x-min(x))/(max(x)-min(x)))
head(USArrests.norm,10)
summary(USArrests.norm)					 	# min = 0, max =1
dim(USArrests.norm)					     	# zadne chybejici hodnoty = dimenze datoveho souboru zustava stejna

dist.USArrests.stan<-dist(USArrests.norm,method='euclidean') # vypocet asociacni matice na normalizovanych 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)					      # 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 do dat vektor prirazeni objektu do shluku 
head(USArrests)
table(USArrests$ward.4)                       # pocty v jednotlivych shlucich
USArrests[order(USArrests$ward.4),] 
# POZOR - v dendrogramu interpretujeme vzdalenosti, pri kterych dochazi ke shlukovani, na poradi shluku ale nezalezi 
# = tj. ruzne funkce mohou shluky seradit odlisne!

# Funkce fviz_dend pouzila pro 3. skupinu zelenou barvu (3. barvu palety), ale funkce cutree skupinu techto statu bere jako 4. shluk!
# Proto pro dalsi vykresleni napevno nadefinujeme sloupec "col" do dat, aby dalsi graficke vystupy byly srovnatelne s dendrogramem!

USArrests$col<-NA
USArrests$col[which(USArrests$ward.4==1)]<-brewer.pal(4,"Set1")[1]
USArrests$col[which(USArrests$ward.4==2)]<-brewer.pal(4,"Set1")[2]
USArrests$col[which(USArrests$ward.4==3)]<-brewer.pal(4,"Set1")[4]
USArrests$col[which(USArrests$ward.4==4)]<-brewer.pal(4,"Set1")[3]
head(USArrests)

# Interpretace shluku = tj. v kterych parametrech se skupiny statu lisi?

# maticovy graf
windows()
pairs(USArrests[,1:4], panel = panel.smooth, main="Simple Scatterplot Matrix",col=USArrests$col) # maticovy graf

# jednorozmerna statistika: boxplot + popis rozlozeni + testovani hypotez
boxplot(USArrests$Murder~USArrests$ward.4) 
aggregate(USArrests$Murder ~ USArrests$ward.4, FUN = function(x) c(Median=median(x),Q25=quantile(x,0.25),Q75=quantile(x,0.75)))
kruskal.test(Murder~ward.4,data=USArrests)

# profil shluku
windows(width=600,height=400)
par(mfrow=c(2,2))
boxplot(USArrests$Murder~USArrests$ward.4,main="Murder",xlab="cluster") 
boxplot(USArrests$Assault~USArrests$ward.4,main="Assault",xlab="cluster") 
boxplot(USArrests$UrbanPop~USArrests$ward.4,main="UrbanPop",xlab="cluster") 
boxplot(USArrests$Rape~USArrests$ward.4,main="Rape",xlab="cluster") 

# Samostatny ukol:
# Provedte shlukovou analyzu na datovem souboru mtcars:
# 1) Zvolte vhodnou miru asociace mezi objekty, spocitejte asociacni matici. Data pred timto standardizujte na z-skore.
# 2) Vyzkousejte algritmus nejvzdalenejsiho souseda a Wardovu metodu. Vykreslete dendrogramy. 
#    Ktera metoda dava lepsi 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)

################################################
#### Vyber vhodneho shlukovaciho algoritmu #####
################################################

# Pomoci kofeneticke matice
# 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 (srovnatelne s Wardovou metodou)

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

################################################
########## Vyber vhodneho poctu 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.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 Wardovu metodu.

####################
# 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 2 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 = 2,						# staty budou rezdeleny do 3 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 Wardovu metodu, pomoci Mantelova testu.
# Porovnejte rozdeleni dat podle Mantelova testu a indexu siluety.