# Tento pomocny skript nemusite pouzit, jeho pouziti je dobrovolne.
# Obsahuje reseni v R s vynechanymi castmi oznacenymi "...".
# Na mista "..." musite doplnit spravnou cast kodu, aby vam skript fungoval.

# Pro spousteni primo ze skriptu, nastavte kurzor na pozadovany radek 
# a stisknete Ctrl + Enter.

#############################################################################
# 8. cviceni
#############################################################################

# Nejprve musime nacist data. Jako obyvkle vyuzijeme funkci read.table(). 
# Jejimi argumenty v nasem pripade jsou cesta k datovemu soboru 
# a argument header udavajici, zda nas datovy soubor obsahuje hlavicku. 
# To zjistite, pokud si soubor "cneck.txt" otevrete napr. 
# v poznamkovem bloku.

# Doplnete (v uvozovkach!) cestu k datovemu souboru a hodnotu T/F, podle toho, 
# jestli soubor ma hlavicku
cneck <- read.table("...", header = ...)

# Podivejte se na zakladni udaje o datech
summary(...)

# Budeme pracovat jen se spojitymi promennymi, proto pomoci hranatych zavorek vybereme
# sloupecky 3 az 8.
data <- cneck[,...]

# 1

# Vykresleni bodoveho diagramu pro data provedeme pomoci funkce plot
plot(...)

# Vidime, ze mezi nekterymi promennymi je silna zavislost - napr. antb.C a neck.C
# "Skloneni bodu" z leveho dolniho rohu do praveho horniho rohu indikuje silnou
# pozitivni korelaci, naopak "Skloneni bodu" z leveho horniho rohu do praveho dolniho
# rohu by indikovalo silnou negativni korelaci.

# 2

# K vykresleni krabicovych diagramu pouzijeme funkci boxplot()
boxplot(...)
## Slouzi pouze pro ilustraci, protoze nemame standardizovana data. Tj. kazda promenna
## je zatim v jinych jednotkach.

# 3

# Vypocitame korelacni matici - vyuzijeme funkci cor()
R <- cor(...)

# Provedeme Bartlettuv test o uplne nezavislosti promennych.

# Nejprve jej provedeme maunualne. Pomoci funkce nrow() ziskame pocet
# pozorovani (= pocet radku naseho souboru), funkce ncol() zjisti pocet sloupecku, 
# tj. pocet promennych.
n <- nrow(...) 
k <- ncol(...) 

## Testovou statistiku vypocteme podle vzorecku z prednasky
test.stat <- -n*log(det(R))*(1- (2*k+11)/(6*n))
## 597.0546

# Prislusny kvantil vypocteme pomoci funkce qchisq(). Prvnim argumentem je koeficient
# spolehlivosti 0.95, druhym stupne volnosti
kvantil <- qchisq(..., df = k*(k-1)/2)
## 24.99579
## kriticky obor je W = < 24.99579; Inf )
## Sami si interpretujte vysledky testu:
## ...

# Tento test muzeme take provest v R pomoci funkce cortest.bartlett() z knihovny psych.
library(psych)
cortest.bartlett(R,n = n,diag = T)

# 4

# PCA s v R provadi pomoci funkce prcomp(), do ktere vlozime nas datovy soubor.
# Parametr center = T udava, zda chceme pracovat s centrovanymi daty (vetsinou ano, protoze
# nas datovy soubor neni vycentrovany) a parametr scale. = T zase udava, zda chceme pracovat 
# s korelacni matici.
data.PCA <- prcomp(data, center = ..., scale. = ...)

# 5

# Podil variability a kumulativni podil variability pro jednotlive komponenty
# ziskame nejjednoduseji pomoci funkce summary()
summary(...)

# 6

# Pokud bychom chteli vysvetlit alespon 80% celkove variability, ihned vidime, ze
# staci pracovat pouze s ... hlavnimi komponentami.

# Vykreslime si do grafu rozptyly jednotlivych komponent (tzv. sutovy/sutinovy graf).
# Pouzijeme funkci plot na list obsahujici hlavni komponenty a zadame typ grafu jako "l".
plot(..., type = "...")

# Abychom dobre videli, ktere rozptyly jsou nad 1 (kvuli Kaiserova kriteria), vykreslime
# jeste vodorovnou primku v 1. Vyuzijeme funkci abline, argument h nastavime na 1, a typ 
# primky (lty) jako 2.
abline(h = ..., lty = ...)

# Kolik komponent bychom vybrali pomoci Kaiserova kriteria?
# Kolik komponent byste vybrali subjektivnim posouzenim sutoveho grafu?
# ...

# 7

# Budeme tedy pracovat s prvnimi dvema komponentami. Pro vypocet korelace komponent
# a puvodnich promennych nejprve vezmeme pozorovani v souradnich hlavnich komponent.
# To najdeme pod $x v nasem listu.

data.in.pc <- data.PCA$... 
## pozorovani v souradnicich hlavnich komponent

# Dale spocteme korelaci techto pozorovani pro prvni dve komponenty s datovym souborem.
# Opet vyuzijeme funkci cor().
cor(..., ...[,1:2])

# 8

# Pro lepsi nahled na celou situaci si vykrelime i pozorovani a promenne v rovine
# prvnich dvou hlavnich komponent, tzv. biplot. Argumentem stejnojmene funkce je 
# list obsahujici hlavni komponenty a parametr scale, kterym graf skalujeme. Nastavime 
# jej na nulu, pouzivame-li korelacni matici pro vypocet komponent. Rovnez nastavime
# popisky pozorovani pomoci xlabs. Chceme aby se pozorovani rozlisila podle pohlavi, 
# vyuzijeme tedy faktoru sex v puvodni datove tabulce cneck.
biplot(..., scale = ..., xlabs = cneck$...)

# Z korelaci i z biplotu vidime, ze prvni komponenta ma vysoke zaporne korelace s vahou, 
# obvodem pasu, obvodem predlokti ci obvodem krku (sipky jsou "priblizne rovnobezne" 
# s osou PC1 - cim "vice jsou sipky rovnobeznejsi", tim vetsi ma dana promenna korelaci
# s PC1). Na druhou hlavni komponentu ma zase nejvetsi vliv telesna vyska a obvod boku.

# 9

# K vypoctu reprodukovane korelacni matice budeme potrebovat nami vybrany pocet hl. komponent
# a prislusne rozptyly (tj. vlastni cisla) k nim. Tyto rozptyly muzeme ziskat dvema zpusoby.
# Budto vezmeme smerodatne odchylky z naseho listu data.PCA a umocnime je na druhou, anebo muzeme
# vyuzit funkci eigen() na korelacni matici. Z jejiho vystupu vezmeme cast pod $values.
# Sami si zkuste, ze obe varianty vam daji stejny vysledek.

vl.cis <- (data.PCA$...^2)[1:2]
vl.cis <- eigen(R)$...[1:2]

# Hlavni komponenty jsou uvedeny jako sloupecky data.PCA$rotation, vezmeme proto jen prvni dva.
# Vyber provedeme pomoci hranatych zavorek. Vlastni cisla musime naskaladat do diagonalni matice, 
# a tak vyuzijme funkci diag(). Vse spolu vynasobime maticovym nasobenim (funkce t() vytvari transpozici).
R.reproduced <- ...[,1:2] %*% diag(...) %*% t(...[,1:2])

# Rezidualni korelacni matici ziskame jako rozdil korelacni matice a reprodukovane
# korelacni matice
R.residual <- ...

# Chceme, aby prvky matice R.residual byly co nejmensi. Pokud si matici vypiseme, 
# tak vidime, ze maximalni prvek v absolutni hodnote je 0.21