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

# nacteni balicku factoextra
install.packages("devtools",'vegan')
devtools::install_github("kassambara/factoextra")
install.packages("factoextra")
install.packages("ggplot2")
library("ggplot2")
library("factoextra")
library("vegan")

install.packages(c('vegan', 'gclus', 'cluster', 'vegan', 'ade4','rafalib'))

library(ade4)
library(vegan)
library(rafalib)

# nastaveni cesty
setwd("C:/Users/ucitel/Desktop/3. cviceni")
getwd()
dir()

lokality<-read.table('lokality.csv', sep=";", header=TRUE, dec=".")

head(lokality,10)
dim(lokality)
is(lokality)
names(lokality)

row.names(lokality)<-lokality[,1]				# ulozeni nazvu lokalit do nazvu radku
lokality<-lokality[,-1]
head(lokality,10)

# PCA pro charakteristiky prostredi ----------------------------------------------------
# Priklad 1: --------------------------------------------------------------
# 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 shepparduv 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. Jaky má vztah k hodnote vlastniho cisla? = odmocnina z vlastniho cisla
lokality.pca$rotation	#vlastni vektory ve sloupcich - normalizovany na jednotkovou delku
sum(lokality.pca$rotation[,1]^2)
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 
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

summary(lokality.pca)		#% rozptylu popsaneho novymi osami - 2. a 3. radek

eig.val <- get_eigenvalue(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

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

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

var <- get_pca_var(lokality.pca) #extrahujeme si vysledky do listu
var

var$coord			#koordinaty promennych na novych osach = korelace promenne s novou osou
var$cor			#korelace puvodnich promennych s novymi osami
var$cos2			#korelace na druhou - vyjadruje % rozptylu promenne popsane novou osou = communalities
var$contrib			#prispevek promennych k definici novych os 

# 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)) 
fviz_pca_var(lokality.pca,axes=c(2,3))

# 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

# communalities
var$cos2 #korelace na druhou - vyjadruje rozptyl puvodni promenne popsane novou osou
sum(var$cos2[1,])
# soucet v radku musi byt roven jedne - nove osy jsou na sobe nezavisle 
# napr. vzdalenost od pramene je nejlepe popsana prvni osou, zatimco pH je nejlepe popsano az 3. faktorovou osou

# naopak % rozptylu faktorovych os, ktere je vysvetleno puvodnimi promennymi = soucet ve sloupci dava 100 %
var.cos2<-var$cos2
comp.cos2 <- apply(var.cos2, 2, sum) # rozptyl popsany osou = vlastni cislo
contrib <- function(var.cos2, comp.cos2){(var.cos2/comp.cos2)*100}
var.contrib <- t(apply(var.cos2,1, contrib, comp.cos2))
var.contrib
sum(var.contrib[,1])
#soucet ve sloupci dava 100 % - rika nam, ktere puvodni promenne definuji smer os 
# 3. osa je definovana z 82 % na zaklade pH

#biplot vzdalenosti

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")) 
fviz_pca_ind(lokality.pca,axes=c(1,3),geom=c("point", "text")) 

#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 24 ma extremne vysokou kladnou hodnotu na 1. faktorove ose. 
# Tedy bude nabyvat vysokych hodnot u promennych, ktere jsou kladne korelovane s touto 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 shepparduv diagram.

lokality.stan<-decostand(lokality, "standardize")
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) #na 1. ose

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

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

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

## 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")
distance.result
pairs(distance.result,data=distance.result, diag.panel = panel.hist, main="Shepparduv diagram")

cor(distance.result) # [1,5] - jednotkova korelace = pokud vebereme pocet os rovny poctu promennych v puvodnim souboru, popiseme 100 % variability dat

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

# Priklad 2 ---------------------------------------------------------------
# Srovnejte vystup z PCA s vystupem ze shlukove analyzy.
# Samostatny ukol: Objekty shlukujte na zaklade Euklidovske vzdalenosti pomoci Wardovy metody do 4 shluku. 
# Prislusnost lokalit do shluku ulozte do vektoru s nazvem "gr".

lokality.stan<-decostand(lokality,"standardize")
dist.lokality.stan<-dist(lokality.stan,method='euclidean')
as.matrix(dist.lokality.stan)
dist.lokality.ward<-hclust(dist.lokality.stan, method="ward.D2") # vypocet shlukove analyzy
                                                    
plot(dist.lokality.ward)    # data rozdelime do 4 shluku
gr<-cutree(dist.lokality.ward, k=4) # rozdeleni souboru
myplclust(dist.lokality.ward,labels=row.names(lokality),lab.col=gr)

site.sc<-lokality.pca$x	    # extrahujeme si score pro lokality

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

# zaver: pomoci PCA lze nalezt shluky v datech. Vyhoda: pri vyhodnoceni korelaci s novymi osami zname i promenne, podle kterych se objekty shlukuji.

# ------------------------------------------------------------------------------------
# Samostatny ukol 1 ------------------------------------------------------------------
# Provedte analogicky PCA (od promenne "Sales in thousands") na datovem souboru car_sales.csv

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

# dale si vlozime nazev auta z promenne Model jako nazev radku -> abychom v biplotu meli nazvy aut
rownames(cars)<-cars$Model
cars[cars$Model=="Neon",] # mame tam duplicitu, dve auta maji stejny nazev "Neon", je potreba jej prejmenovat

# prejmenovani se musi udelat tak, ze nejdrive je potreba pridat novy nazev mezi faktory promenne (nazvy tech aut)
levels(cars$Model) <- c(levels(cars$Model),"Neon_2")
cars["78",2]<-"Neon_2" # pak uz jen jednoduse nazev prejmenujeme

rownames(cars)<-cars$Model # vlozime nazvy aut z promenne Model jako nazvy radku
head(cars)

# PCA od "Sales in thousands"
head(cars)
cars<-cars[,3:ncol(cars)]

# -> 0) V datech jsou chybejici hodnoty - pred analyzou je potreba tyto odstranit. 
summary(cars)
which(apply(is.na(cars),1,sum)>0)

cars<-cars[which(apply(is.na(cars),1,sum)==0),] # vybere pouze kompletne vyplnena data
dim(cars)

# -> 1) Spocitejte PCA. Jakou matici pouzijete na vstupu? 
summary(cars)
cars.pca<-prcomp(cars, scale=TRUE)#pouzijeme matici korelaci
cars.pca	 

cars.pca$sdev 		#Standard deviations: smerodatna odchylka pocitana na novych osach = odmocnina z vlastniho cisla
cars.pca$rotation		#vlastni vektory ve sloupcich - normalizovany na jednotkovou delku
sum(cars.pca$rotation[,1]^2)
cars.pca$center		#prumery puvodnich promennych
cars.pca$scale		#smerodatna odchylka puvodnich promennych
cars.pca$x			#pozice objektu na novych osach

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

# vlastni cisla a % rozptylu popsaneho novymi osami

eig.val <- get_eigenvalue(cars.pca) 
head(eig.val) #interpretujeme prvni 2 hlavni komeponenty, ktere popisuji 80 % rozptylu dat

# scree plot
screeplot(cars.pca,type="lines",main="Scree plot") #na ose y mame rozptyl na novych osach = vlastni cisla

eig[eig > mean(eig)]  # prvni 2 hlavni komeponenty popisuji vice rozptylu nez 1 puvodni promenna


# -> 3) Vyhodnotte vztah promennych k novym faktrovym osam - communalities a % rozptylu nove osy popsane puvodnimi promennymi. 

Otazka c. 1: Kterou osou je nejlepe popsano 'Horsepower' a 'Wheelbase'?
Otazka c. 2: Ktera promenna nejvice prispiva k definici 1. a 2. faktorove osy? 

var <- get_pca_var(cars.pca) #extrahujeme si vysledky do listu
var

var$cos2
# napr. 'Horsepower', 'engine size' jsou nejlepe popsany prvni osou, zatimco 'Wheelbase' 2. faktorovou osou
var$contrib
# 1. osa je z nejvetsi casti popsana promennou 'Curb.weight' a 2. osa promennou 'X4.year.resale.value'

# -> 4) Vykreslete biplot korelaci a vzdalenosti.
# Vyhodnotte velikost motoru ('engine size') a efektivitu paliva ('Fuel efficiency') modelu S-class na zaklade pozic modelu na 
# faktorovych osach a korelaci puvodnich promennych s novymi faktorovymi osami.

#biplot - promenne i objekty zaroven
biplot(cars.pca, cex = 0.8, col = c("black", "red") )

# interpretace: promenne, ktere se spolu vyskytuji spolecne v jednom kvadrantu - kladna korelace
# napr. cim je auto delsi, tim bude tezsi (kladna korelace)
# promenne, ktere se vyskytuji na opacnych stranach od pocatku - zaporna korelace (napr. ucinnost paliva vs. velikost motoru)

# promenne korelovane s prvnimi osami predstavuji nejvetsi zdroj variability (napr. Fuel efficiency, fuel capacity)
# u promennych, ktere nejsou korelovany s faktorovymi osami nebo jsou korelovany s poslednimi osami, lze zvazit jejich vyrazeni z analyzy

# model S-Class bude mit veliky motor (model ma zapornou hodnotu na 1. faktorove ose a promenna 'engine size' je 
# zaporne korelovana s touto osou - hodnota v tomto paramteru bude tedy opacna = vysoka) a nizkou efektivitu paliva
#(promenna 'Fuel efficiency' je kladne korelovana s osou - hodnota v tomto parametru bude stejneho smeru jako je hodnota na ose - nizka)


# -> 5) Vykreslete shepparduv diagram. 
# Interpretujte, zda Vami zvoleny pocet os v novem prostoru dobre odpovida vzdalenostem v puvodnim prostoru.

summary(cars)
cars.stan<-decostand(cars, "standardize")
dist.cars<-dist(cars.stan,method='euclidean') # vypocet asociacni matice na zaklade Euklidovske vzdalenosti
dist.cars<-as.vector(dist.cars) #puvodni vzdalenosti v prostoru

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

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

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

dim(cars)
dist.cars.pca11<-dist(cars.pca$x[,1:11],method='euclidean')
dist.cars.pca11<-as.vector(dist.cars.pca11) #na vsech faktorovych osach

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

#vidime, ze vzdalenost objektu pri interpretaci prvnich 2 os dobre koresponduje vzdalenostem v puvodnim prostoru na vsech promennych


