##########################################################################################################
###############################    Analýza hlavních komponent    #########################################
##########################################################################################################

# nacteni balicku 
install.packages("factoextra") # pro hezke graficke vystupy
library("factoextra")

install.packages("magrittr") # pro funkci %>%
library(magrittr)

### reseni k prikladu z prezentace
####################################
dataZK<-data.frame(matrix(c(101,105,103,98,93,16,18,42,23,6),byrow=F,ncol=2))
dataZK
var(dataZK) # kovariancni matice 
pcaZK<-prcomp(dataZK, scale=F)


(pcaZK$sdev)^2	# vlastni cisla
pcaZK$rotation	# vlastni vektory
pcaZK$x		      # pozice objektu na novych osach
summary(pcaZK)  # 1. osa popisuje 92.9% rozptylu dat, obe osy popisuji 100%, protoze mame pouze dva vstupni parametry

windows(width=5,heigh=5)
plot(dataZK[,1],dataZK[,2],xlab="Objem šedé hmoty",ylab="Objem likvoru",main="Data v původním prostoru",pch=19,cex=1.2)
windows(width=5,heigh=5)
plot(pcaZK$x[,1],pcaZK$x[,2],xlab="1. hl. komponenta",ylab="2. hl. komponenta",main="Data v prostoru nových os z PCA",pch=19,cex=1.2)
####################################

# nastaveni cesty
setwd("C:/Users/cvanova/Desktop/05_cviceni_shlukovka")
getwd()
dir()	# vypise nazvy souboru v nastavenem adresari

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

# Priklad 1: --------------------------------------------------------
# PCA pro charakteristiky lokalit(vodnich toku) 
# das - vzdalenost od pramene (km * 10)
# alt - nadmorska vyska (m)
# pen (ln(x + 1), x = sklon (per mil * 100)
# deb - minimalni prumerny prutok (m3/s *100)
# pH (* 10)
# dur - celkova tvrdost vody (mg/l of Calcium)
# pho - mnozstvi fosfatu (mg/l * 100)
# nit - mnozstvi nitratu (mg/l * 100)
# amm - mnozstvi amoniaku (mg/l * 100)
# oxy - kyslik ve vode (mg/l * 10)
# dbo - biologicky pozadavek kysliku (mg/l * 10)

# Na datovem souboru lokality provedte analyzu hlavnich komponent.
# 1) Jakou matici pouzijete na vstupu?
# 2) Kolik hlavnich komponent budete interpretovat? Vykreslete scree plot.
# 3) Vyhodnotte vztah promennych k novym faktrovym osam. Vykreslete biplot korelaci a vzdalenosti.
# 4) Vykreslete a komentujte Sheparduv diagram.

# -> 1) Jakou matici pouzijete na vstupu? 
?prcomp
lokality.pca<-prcomp(lokality, scale=TRUE)  # pouzijeme matici korelaci (pokud bychom nastavili scale=F, spocitala by se kovariancni matice)
lokality.pca	 

lokality.pca$sdev 		                      # Standard deviations: smerodatna odchylka pocitana na novych osach = odmocnina z vlastniho cisla
lokality.pca$rotation                       # vlastni vektory ve sloupcich - normalizovany na jednotkovou delku
lokality.pca$center		                      # prumery puvodnich promennych
lokality.pca$scale		                      # smerodatna odchylka puvodnich promennych
lokality.pca$x			                        # pozice objektu na novych osach

# -> 2) Kolik hlavnich komponent budete interpretovat? Vykreslete scree plot.							   

# vlastni cisla a % rozptylu popsaneho novymi osami

lokality.pca$sdev 	                      	# Standard deviations: smerodatna odchylka pocitana na novych osach = odmocnina z vlastniho cisla
eig <- (lokality.pca$sdev)^2	              # Vlastni cisla
eig 
variance <- eig*100/sum(eig)	              # % rozptylu popsaneho novymi osami
variance
cumvar <- cumsum(variance)	                # % rozptylu popsaneho novymi osami - kumulativne
cumvar 
summary(lokality.pca)		                    # % rozptylu popsaneho novymi osami - 2. a 3. radek

lokality.pca.result <- data.frame(eig = eig, variance = variance,cumvariance = cumvar)
head(lokality.pca.result)	                  # v jedne tabulce vedle sebe - vlastni cisla, % rozptylu, % rozptylu kumulativne scitano od 1. osy

eig.val <- get_eig(lokality.pca)     # nebo jednoduse pomoci funkce 'get_eigenvalue', ktera je soucasti balicku 'factoextra'
head(eig.val)

# scree plot

par(mfrow=c(1,2)) 
barplot(lokality.pca.result[, 2], names.arg=1:nrow(lokality.pca.result), ylim=c(0,60),
       main = "Scree plot",
       xlab = "Principal Components",
       ylab = "Percentage of variances",
       col ="steelblue")
# Add connected line segments to the plot
lines(x = 1:nrow(lokality.pca.result), 
      lokality.pca.result[, 2], 
      type="b", pch=19, col = "red")        # prvni tri osy popisuji temer 90 % rozptylu puvodnich dat

screeplot(lokality.pca,type="lines",main="Scree plot") # nebo pomoci funkce screeplot, zde na ose y mame rozptyl na novych osach = vlastni cisla

# Aplikace Kaiserova-Guttmanova kriteria - inerpretujeme tolik vlastnich hodnot, ktere maji vyssi hodnotu nez je prumerna hodnota vlastnich cisel v souboru
eig[eig > mean(eig)]

# pocet vlastnich hodnot > 1 indikuje pocet novych os, ktere popisuji varibilitu dat lepe nez jedna puvodni promenna: vybrali bychom 2 osy
# nebo si muzeme stanovit % z celkoveho rozptylu dat, ktere potrebujeme novymi osami popsat - napr. 80-90 %: vybrali bychom 3 osy

# stejne pomoci funkci balicku 'factoextra'
fviz_screeplot(lokality.pca, ncp=10)        # ncp - pro prvnich 10 novych os
fviz_screeplot(lokality.pca, ncp=10, choice="eigenvalue")

# -> 3) Vyhodnotte vztah promennych k novym faktrovym osam. Vykreslete biplot korelaci a vzdalenosti.
?get_pca_var
var <- get_pca_var(lokality.pca)            # extrahujeme si vysledky
var

var$coord			                              # koordinaty promennych na novych osach = korelace promenne s novou osou
var$coord^2		                              # rozptyl nove osy PCA, ktery lze popsat pomoci puvodni promenne (sdileny rozptyl)
sum((var$coord^2)[,1])                      # rozptyl na 1. ose PCA = vlastni cislo asociovane s prvni osou
sum((var$coord^2)[1,])                      # rozptyl puvodni promenne

# Jak je popsan rozptyl puvodni promenne novymi osami? 
(var$cos2	)*100                             # communalities = vyjadruje % rozptylu promenne popsane novou osou
# vzdalnost od pramene (das) je z 77% popsana  1. osou, pH je z 82% popsano 3. osou

# Ktere promenne prispeli k definici nove osy?
var$contrib			                            # prispevek promennych k definici novych os 
# % rozptylu faktorovych os, ktere je vysvetleno puvodnimi promennymi = soucet ve sloupci dava 100 %
# 1. osa je definovana parametry: nit, das, alt, pho

# biplot korelaci
a <- seq(0, 2*pi, length = 100)
plot( cos(a), sin(a), type = 'l', col="gray",xlab = "PC1",  ylab = "PC2")
abline(h = 0, v = 0, lty = 2)
arrows(0, 0, var$coord[,1], var$coord[, 2], length = 0.1, angle = 15, code = 2) #vykresli korelace promennych s novymi osami
text(var$coord, labels=rownames(var$coord), cex = 1, adj=1)				              #prida popisky promennych

fviz_pca_var(lokality.pca,axes=c(1,2))     # nebo jednim prikazem z balicku 'factoextra'
fviz_pca_var(lokality.pca,axes=c(1,3))     # 1. vs 3. osa
fviz_pca_var(lokality.pca,axes=c(2,3))     # 2. vs 3. osa

# interpretace: 
# promenne, ktere se spolu vyskytuji spolecne v jednom kvadrantu - kladna korelace
# promenne, ktere se vyskytuji na opacnych stranach od pocatku - zaporna korelace 
# napr. mnozstvi kysliku - 'oxy' bude nizsi, cim vice bude ve vode fosfatu - 'pho' a dusicnanu - 'amm', ktere souvisi s vyskytem sinic)
# cim vetsi bude nadmorska vyska - 'alt', tim kratsi bude vzdalenost od pramene - 'das'
# cim vetsi je sklon - 'pen' (lokalita se bude vyskytovat v horach), tim vyssi bude nadmorska vyska

# promenne korelovane s prvnimi osami predstavuji nejvetsi zdroj variability
# u promennych, ktere nejsou korelovany s faktorovymi osami nebo jsou korelovany s poslednimi osami, lze zvazit jejich vyrazeni z analyzy

#biplot vzdalenosti
dim(lokality.pca$x)
plot(lokality.pca$x[,1],lokality.pca$x[,2], pch = 19,
	xlab=paste("PCA 1: ",round(variance[1],1),"%",sep=""),
	ylab=paste("PCA 2: ",round(variance[2],1),"%",sep=""))
abline(h=0, v=0, lty = 2)
text(lokality.pca$x[,1],lokality.pca$x[,2], labels=rownames(lokality),cex=0.7, pos = 3)

fviz_pca_ind(lokality.pca,axes=c(1,2),geom=c("point", "text")) # nebo jednim prikazem z balicku 'factoextra'
fviz_pca_ind(lokality.pca,axes=c(2,3),geom=c("point", "text")) # 2. vs 3. osa
fviz_pca_ind(lokality.pca,axes=c(1,3),geom=c("point", "text")) # 1. vs 3. osa 

# interpretace: 
# hledame shuky/odlehle hodnoty
# cim je objekt vice vzdalen od pocatku, tim bude mit extremnejsi profil

#biplot - promenne i objekty zaroven
biplot(lokality.pca, cex = 0.8, col = c("black", "red") )
fviz_pca_biplot(lokality.pca, geom = "text") #habillage - obarveni podle shluku

# interpretujme 1. faktorovou osu. Lokalita 25 ma extremne vysokou kladnou hodnotu na 1. faktorove ose. 
# Tedy bude nabyvat vysokych hodnot u promennych, ktere jsou kladne korelovane s faktorovou osou (nit, pho, amm, das, deb, dur)
# a nizkych hodnot se zaporne korelovanymi promennymi (alt, pen, oxy). 
# Lokalita 1 ma extremne nizkou (zapornou) hodnotu na 1. faktorove ose. 
# Tedy bude nabyvat nizkych hodnot u promennych, ktere jsou kladne korelovane s faktorovou osou (nit, pho, amm, das, deb, dur)
# a vysokych hodnot se zaporne korelovanymi promennymi (alt, pen, oxy).

# -> # 4) Vykreslete Sheparduv diagram.

lokality.stan <- lokality%>%
  na.omit() %>%          					                                       # odstrani chybejici hodnoty
  scale()	       				                                                 # standardizuje data
dist.lokality<-dist(lokality.stan,method='euclidean')                    # vypocet asociacni matice na zaklade Euklidovske vzdalenosti
dist.lokality<-as.vector(dist.lokality)                                  # puvodni vzdalenosti v prostoru

dist.lokality.pca1<-dist(lokality.pca$x[,1],method='euclidean')
dist.lokality.pca1<-as.vector(dist.lokality.pca1)                        # vzdalenosti na 1. ose

dist.lokality.pca2<-dist(lokality.pca$x[,1:2],method='euclidean')
dist.lokality.pca2<-as.vector(dist.lokality.pca2)                        # vzdalenosti na 1.-2. ose

dist.lokality.pca3<-dist(lokality.pca$x[,1:3],method='euclidean')
dist.lokality.pca3<-as.vector(dist.lokality.pca3)                        # vzdalenosti na 1.-3. ose

dist.lokality.pca11<-dist(lokality.pca$x[,1:11],method='euclidean')
dist.lokality.pca11<-as.vector(dist.lokality.pca11)                      # vzdalenosti na 1.-11. ose

## 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", ...)
}

distance.result<-data.frame(cbind(dist.lokality,dist.lokality.pca1,dist.lokality.pca2,dist.lokality.pca3,dist.lokality.pca11))
colnames(distance.result)<-c("puvodni vzdalenost","PCA 1","PCA 1:2","PCA 1:3","PCA 1:11")
head(distance.result)
pairs(distance.result,data=distance.result, diag.panel = panel.hist, main="Sheparduv diagram")

cor(distance.result[c(1,5)]) 
# jednotkova korelace = pokud vybereme pocet os rovny poctu promennych v puvodnim souboru, popiseme 100 % variability dat 
# => korelace puvodnich vzdalenosti a vzdalenosti na vsech osach PCA bude rovna 1 (vzdalenosti si jsou rovny)

##########################################################################################################
############################   Kombinace shlukove analyzy a PCA     ######################################
##########################################################################################################

# Priklad 2 ---------------------------------------------------------------
# Srovnejte vystup z PCA s vystupem ze shlukove analyzy a zjistete, zda je v datech nejaka struktura.
# Objekty shlukujte na zaklade Euklidovske vzdalenosti pomoci Wardovy metody

lokality.w<-hclust(dist(lokality.stan), "ward.D2") # vypocet shlukove analyzy
                                                    
plot(lokality.w)            # data rozdelime do 4 shluku
gr<-cutree(lokality.w, k=4) # rozdeleni souboru

# Vykreslime graf z PCA s obarvenymi lokalitami dle prislusnoti do shluku
fviz_pca_ind(lokality.pca, geom=c("point", "text"), habillage=gr,addEllipses=TRUE)



# ------------------------------------------------------------------------------------
# Samostatny ukol 1 ------------------------------------------------------------------
# Provedte PCA na datovem souboru "car_sales.csv".

# Manufacturer - vyrobce
# Model - model
# Sales in thousands - pocet prodanych aut (tis.)
# 4-year resale value - hodnota auta po 4 letech
# Price in thousands - cena (tis.)
# Engine size - velikost motoru
# Horsepower - vykon motoru
# Wheelbase - delka auta (mezi koly)=rozvor
# Width - sirka auta
# Length - delka auta
# Curb weight - vaha auta
# Fuel capacity - objem palivove nadrze
# Fuel efficiency - efektivita paliva

readLines("car_sales.csv",3)
cars<-read.table('car_sales.csv', sep=";", header=TRUE, dec=".",row.names=2)
head(cars)
summary(cars)
cars <- cars%>%  na.omit()  
dim(cars)

# Ukoly:
# 1) Spocitejte PCA od promenne "Sales in thousands". Jakou matici pouzijete na vstupu?

# 2) Kolik hlavnich komponent budete interpretovat? 
#  - vypiste hodnoty vlastnich cisel
#  - vypiste kumulativni % rozptylu popsane osami 
#    (kolik os byste museli interpretovat, abyste popsali alespon 70% rozptylu puvodnich dat)
#  - dle Keiser-Gutmanova kriteria se rozhodnete, kolik os budete interpretovat
#  - Vykreslete a interpretujte scree plot (sutinovy graf) a Sheparduv diagram (do vzdalenosti na 4 hlavnich komponentach).

# 3) Vykreslete biplot korelaci (pro 1. a 2 osu PCA) a interpretujte:
#       - jaka je korelace velikosti motoru a objemu palivove nadrze? (kladna x zaporna x nelze urcit)
#       - jaka je korelace efektivity paliva a velikosti motoru? (kladna x zaporna x nelze urcit)
#       - jaka je korelace vahy a sirky automobilu? (kladna x zaporna x nelze urcit)
#       - ktera charakteristika je nejvice kladne korelovana s 1. hlavni komponentou?
#       - ktera charakteristika je nejvice zaporne korelovana s 2. hlavni komponentou?

# 4) Dle vztahu parametru k novym osam interpretujte:
#       - ktera osa nejlepe a z kolika % popisuje rozptyl "Curb.weight" = vahy automobilu?
#       - ktere 3 charakteristiky nejvice prispely k definici 2. osy?

# 5) Vykreslete biplot vzdalenosti a interpretujte:
#       - ktery model auta ma nejvyssi zapornou pozici na 1. hlavni komponente?
#       - ktery model auta ma nejvyssi kladnou pozici na 1. hlavni komponente?
#       - ktery model auta ma nejvyssi zapornou pozici na 2. hlavni komponente?
#       - ktery model auta ma nejvyssi kladnou pozici na 2. hlavni komponente?

# 6) Vykreslete biplot korelaci a vzdalenosti a interpretujte:
#       - jaka je efektivita paliva, velikost motoru a vaha modelu Viper? (vysoka x nizka)
#       - jaka je efektivita paliva, velikost motoru a vaha modelu Metro? (vysoka x nizka)
#       - vime-li, ze prodejni cena modelu Civic je po 4 letech nizka, jaka bude cena (Price.in.thousands) tohoto modelu? (vysoka x nizka)


# ------------------------------------------------------------------------------------
# Samostatny ukol 2 ------------------------------------------------------------------
# Provedte PCA na datovem souboru "normcho_237_4class.txt".

readLines("normcho_237_4class.txt",3)
geny<-read.table('normcho_237_4class.txt', sep="\t", header=TRUE, dec=".",row.names=1)
head(geny)
dim(geny)  # promenne = geny jsou popsany v radcich, tabulku musime transponovat


geny.T<-t(geny)
head(geny.T)
dim(geny.T) # 17 subjektu je popsano expresi 217 genu

# redukce datoveho souboru - cilem je snizit pocet vstupnich parametru, aby byly vystupy prehlednejsi
tmp <- cor(geny.T)                                             # korelacni matice
tmp[!lower.tri(tmp)] <- 0                                      # unikatni hodnoty pouze pod/nad diagonalou (jinde budou nuly)
geny.t <- geny.T[,!apply(tmp,2,function(x) any(abs(x) > 0.7))] # pro analyzu vybereme pouze geny s max korelaci 0.7
dim(geny.t) 
colnames(geny.t)<-paste("gen",1:40)
head(geny.t)

# Ukoly:
# 1) Spocitejte PCA. Jakou matici pouzijete na vstupu?

# 2) Kolik hlavnich komponent budete interpretovat? 
#  - vypiste hodnoty vlastnich cisel
#  - vypiste kumulativni % rozptylu popsane osami 
#    (kolik os byste museli interpretovat, abyste popsali alespon 50% rozptylu puvodnich dat)
#  - dle Keiser-Gutmanova kriteria se rozhodnete, kolik os budete interpretovat
#  - Vykreslete a interpretujte scree plot (sutinovy graf) a Sheparduv diagram (do vzdalenosti na 6 hlavnich komponentach).

# 3) Vykreslete biplot korelaci (pro 1. a 2 osu PCA) a interpretujte:
#       - jaka je korelace genu 26 a genu 32? (kladna x zaporna x nelze urcit)
#       - jaka je korelace genu 1 a genu 18? (kladna x zaporna x nelze urcit)
#       - jaka je korelace genu 11 a genu 24? (kladna x zaporna x nelze urcit)
#       - ktery gen je nejvice kladne korelovany s 1. hlavni komponentou?
#       - ktery gen je nejvice zaporne korelovany s 2. hlavni komponentou?

# 4) Dle vztahu expresi genu k novym osam interpretujte:
#       - ktera osa nejlepe a z kolika % popisuje rozptyl genu 1
#       - ktere 3 geny nejvice prispely k definici 6. osy?

# 5) Vykreslete biplot vzdalenosti a interpretujte:
#       - ktery subjekt ma nejvyssi zapornou pozici na 1. hlavni komponente?
#       - ktery subjekt ma nejvyssi kladnou pozici na 2. hlavni komponente?

# 6) Vykreslete biplot korelaci a vzdalenosti a interpretujte:
#       - jaka exprese genu 32 bude pozorovana u subjektu c1? (vysoka x nizka)
#       - jaka exprese genu 32 bude pozorovana u subjektu c5? (vysoka x nizka)
#       - jaka exprese genu 31 bude pozorovana u subjektu c2? (vysoka x nizka)
#       - vime-li, ze objekt c6 ma vysokou expresi genu 22, jakych hodnot bude dosahovat exprese genu 35 u tohoto subjektu? (vysokych x nizkych)

