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

#############################################################################
# 11. cviceni / Linearni diskriminacni analyza
#############################################################################
# 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 "logistic-metacarpals.txt" otevrete napr. 
# v poznamkovem bloku.

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

# Pokud chceme v R pouzit diskriminacni analyzu, musi byt promenna, kterou 
# chceme rozlisovat typu faktor (v tomto pripade jde o promennou pohlavi). 
# Zjistime to pomoci funkce is.factor(), kterou aplikujeme na promennou "sex".
is.factor(data$...)

# (1)

# Pocty pozorovani zjistime pomoci funkce table
...(data$sex)

# Prumery jednotlivych promennych ziskame s vyuzitim funkce colMeans() - 
# - musime ji aplikovat jen na ta data, ktera odpovidajivybranemu pohlavi
# a zaroven neobsahuji prvni dva sloupece, ktere obsahuji promemne, 
# jez prumerovat nechceme.
...(data[data$sex == "f",3:6])

# Pro vyberovou varinacni matici vse funguje analogicky. Jen pouzijeme funkci
# cov() namisto colMeans(). Muzete pouzit i funkci var(), dostanete stejny vysledek.
...(data[data$sex == "f", 3:6])

# Pro muze zkuste sami
colMeans(...)
cov(...)

# (2)

# Orientacni overeni linearity provedeme tak, ze si vykreslime bodove grafy
# zvlast pro zeny a zvlast pro muze.
plot(data[data$... == "f", 3:6], main = "Zeny")
plot(data[data$... == "m", 3:6], main = "Muzi")

# Pomoci funkce max(), kterou aplikujeme napr. na promennou mc2.BW u zen, 
# nejprve zjistime nejvyssi hodnotu u mc2.BW u zen.
...(data$mc2.BW[data$sex == "f"])

# Dale pomoci hranatych zavorek najdeme toto pozorovani v celem datovem souboru.
# Musime pouzit logickou spojku "a zaroven"
data[(data$mc2.BW == 16.8 & data$sex == "f"),]

# Vidime, ze jde o pozorovani 43. Vytvorime tedy novy datovy soubor bez tohoto 
# pozorovani a nadale budeme pracovat s novym souborem.
data2 <- data[-...,]

# (3)

# Rovnez odstranime sloupec "id"
data2 <- data2[,-1]

# (4)

# K overeni normality vyuzijeme Mardiuv test i H-Z test a k overeni shodnosti 
# variancnich matic Boxuv M test.

# Z minulych cviceni uz znate syntaxi i vyznam jednotlivych parametru.
# Nacteme knihovnu MVN.
library(...)

# Nastavime obrazek na 1x2.
par(mfrow = c(1,2))

# A provedeme Mardiuv i H-Z test.
mvn(data2, subset = "sex", mvnTest = "...", multivariatePlot = "qq")$multivariateNormality
mvn(data2, subset = "sex", mvnTest = "...")$multivariateNormality
## Interpretujte vysledky obou testu:
## ...

# Pro overeni variancnich matic vsak musime z data2 vybrat jen numericke sloupce.
library(biotools)
boxM(data2[,2:5], grouping = ...)
## Interpretujte sami vysledky testu:
## ...

# (5)

# K otestovani hypotezy o shodnosti strednich hodnot vektoru pro obe pohlavi
# pouzijeme Hotellinguv T2 test z knihovny ICSNP. Zademe pozorovani prislusne
# zenam a muzum - musime vsak vybrat jen numericke sloupce.
library(...)
HotellingsT2(..., ...)
## Interpretujte vysledky testu
## ...

# (6)

# K sestaveni linearni diskriminacni funkce se v R pouziva funkce lda() 
# z knihovny MASS.
# Do ni zadame formuli podobnym zpusobem jako drive napriklad do funkce lm().
library(...)
model.lda <- lda(... ~ ... + ... + ... + ..., data = data2)

# A hned si nas model prohledneme
model.lda
# Jake jsou hodnoty koeficientu?
# ...

# (7)

# Nyni se podivame, jak nas model urcuje pohlavi na zaklade vsech promennych 
# a nasich dat. Vyuzijeme funkci predict().

fitted <- ...(model.lda)

# To jak nas model rozhodl najdeme pod $class. Vypiseme si tedy vsechny udaje
# do tabulky, kde radky budou predikce a sloupecky skutecne hodnoty
tab <- table(fitted$class, data2$sex)
colnames(tab) <- c("skut. f", "skut. m")
row.names(tab) <- c("pred. f", "pred. m")
tab

# Vypocitame podil spravne a mylne zarazenych pozorovani.

# Spravne klasifikovana pozorovani jsou na hlavni diagonale nasi tabulky.
# Secteme je a podelime souctem vsech pozorovani v tabulce.
sum(diag(...)) / sum(...)

# (8)

# Nespravne klasifikovana pozorovani jsou mimo hl. diagonalu
(...[1,2] + ...[2,1]) / sum(tab)

# Tento podil nespravne zarazenych pozorovani muzeme porovnat s nahodnou
# klasifikaci, kde bereme v potaz jen apriorni pravdepodobnosti.

p <- model.lda$prior
2*p[1]*p[2]

# Linearni diskriminacni analyza snizila podil spatneho zarazeni z ... na ...

# Pozn. Pokud bychom uvazovali jen na zaklade apriornich pravdepodobnosti, 
# tj. na zaklade relativnich cetnosti skupin, tak procento SPRAVNE zarazenych 
# pozorovani by bylo:
1 - 2*p[1]*p[2]

# (9)

# Nyni provedeme zarazeni neznameho (noveho) pozorovani. Tj. sitauce je takova, ze
# namerime hodnoty nasich promennych a chteli bychom rozhodnout, zda jde o zenu, ci muze.
# Rozhodujeme na zaklade naseho modelu. Vyuzijeme funkci predict(), s argumentem
# newlist. Do nej vlozime list s hodnotami naseho noveho pozorovani.

...(model.lda, ... = ...(mc2.L = 64, mc2.CW = 14, mc2.MW = 8, mc2.BW = 10))
## Jake jsou posteriorni pravdepodobnosti pro zarazeni do skupiny "f" a "m"?
## ...
## Jak nas model toto nove pozorovani zaradil? 
## ...

# (10)

# Pokud chceme vybrat promenne do naseho modelu automaticky (procedura
# vybere jen ty dulezite prom. pro LDA), muzeme pouzit funkci greedy.wilks()
# z knihovny "klar". Zpusob zadani je stejny jako do funkce lda().
library(...)
greedy.wilks(... ~ ... + ... + ... + ..., data = data2)
## Vidime vsak, ze dopredna krokova metoda vybrala vsechny promenne.

# Priklad 2

# Nacteme datovy soubor. Kvuli odlisnemu formatu vstupnich dat tentokrat vyuzijeme
# funkci read.csv(). Argumentem je datovy soubor Howell.csv, pritomnst hlavicky (T/F)
# a take na.strings, ktery rika, jake hodnoty v datech maji byt brane jako NA - 
# - nastavime na "0".
cranio <- read.csv("...", header = ..., na.strings = "0")

# Z cele obrovske tabulky cranio chceme pouze muze (cranio$Sex == "M"), 
# populace BERG, BURIAT a PERU 
# (Population %in% c("BERG", "BURIAT", "PERU")) 
# a promenne XFB, NPH, NLH, OBH, OBB, MAB, EKB. Vyber provedeme pomoci hranatych zavorek.
# Umistime do nich tedy podminku ve forme logickeho vyrazu a vektor oznacujici, 
# ktere sloupecky tabulky chceme vybrat.
data <- cranio[cranio$Sex == "M" & cranio$Population %in% c("...", "...", "..."),
               c("Population", "XFB", "...", "...", "...", "...", "...", "EKB")]

# Vypiseme si prehled o datech
summary(...)

# Jeste bychom chteli odstranit prebytecne populace z vypisu 
# (ty, co maji cetnost nula). To provedeme tak, ze ze sloupecku populace vytvorime
# novy faktor a jim nahradime puvodni faktor:
data$Population <- ...(data$Population)

# Nyni je vypis podle nasich predstav
summary(data)

# Nasledujici postupy jsou temer analogicke s predchozim pripadem. 
# Zkuste proto postupovat samostatneji.

# (1)

# Pocty pozorovani
...(data$Population)

...(data[data$Population == "BERG",2:8])
...(data[data$Population == "BURIAT",2:8])
...

# Variancni matice
cov(...)
cov(...)
cov(...)

# (2)

plot(data[data$Population == "BERG", 2:8], main = "BERG")
plot(..., main = "BURIAT")
plot(..., main = "PERU")

# (3)

library(...)
par(mfrow = c(1,3))
mvn(data, subset = "Population", mvnTest = "...", multivariatePlot = "qq")$multivariateNormality
mvn(data, subset = "Population", mvnTest = "...")$multivariateNormality
## Interpretujte vysledky obou testu:
## ...

library(biotools)
boxM(data[,...], grouping = ...)
## Interpretujte sami vysledky testu:
## ...

# (4)

# K otestovani hypotezy o shodnosti strednich hodnot vektoru pro vsechny populace
# pouzijeme Wilksovu statistiku. Vyuzijeme funkci manova() jako na minulem cviceni.
model <- ...(as.matrix(data[,-1]) ~ data$Population)
summary(model, test = "Wilks")
## Interpretujte vysledky testu.
## ...

# (5)

# Sestaveni linearni diskriminacni funkce 
library(MASS)
model.lda <- lda(... ~ ... + ... + ... + ... + ... + ... + ..., data = data)
model.lda
## Ve vystupu vidime obe linearni diskriminacni funkce pro vsechny promenne.

# (6)

fitted <- predict(...)

tab <- table(fitted$class, data$Population)
colnames(tab) <- c("skut. BERG", "skut. BURIAT", "skut. PERU")
row.names(tab) <- c("pred. BERG", "pred. BURIAT", "pred. PERU")
tab

# Vypocitame podil spravne zarazenych pozorovani.

sum(...) / sum(...)

# (7)

# Predikce noveho pozorovani
predict(..., newdata = ...(XFB = 124, NPH = 70, NLH = 53, OBH = 34, OBB = 41, MAB = 64, EKB =  100))
## Jake jsou posteriorni pravdepodobnosti pro zarazeni do skupin "BERG", "BURIAT" a "PERU"?
## ...
## Jak nas model toto nove pozorovani zaradil? 
## ...

# (8)

# Dopredna krokova metoda
library(klaR)
greedy.wilks(... ~ ... + ... + ... + ... + ... + ... + ..., data = data)
## Kterou promennou dopredna krokova metoda vynechala?
## ...