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

#############################################################################
# 10. cviceni
#############################################################################
# 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 = T, na.strings = "0")

# Z cele obrovske tabulky cranio chceme pouze muze (cranio$Sex == "M"), 
# populace ZALAVAR, BERG, ESKIMO a NORSE 
# (Population %in% c("ZALAVAR", "BERG", "ESKIMO", "NORSE")) 
# a promenne BPL, NPH, OBH. 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 == "..." & cranio$Population %in% c("...", "...", "...", "..."),
               c("Population", "BPL", "NPH", "OBH")]

# 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 <- factor(data$...)

# Nyni je vypis podle nasich predstav
summary(...)

# (1)

# Pocet pozorovani pro kazdou populaci ziskame pomoci funkce table()
# aplikovane na data$Population.
table(...)


# Prumery jednotlivych promennych pro dane populace ziskame s vyuzitim
# funkce colMeans() - musime ji aplikovat jen na ty data, ktera odpovidaji nasi
# populaci a zaroven neobsahuji prvni sloupec, ktery neni numericky.
...(data[... == "ZALAVAR",-1])

# Pro ostatni populace zkuste sami
colMeans(data[data$Population == "...",-1])
colMeans(data[data$Population == "...",-1])
colMeans(data[data$Population == "...",-1])

# Pro vyberovou varinacni matici vse funguje analogicky. Jen pouzijeme funkci
# var() namisto colMeans()
var(data[data$Population == "...",-1])
var(data[data$Population == "...",-1])
...
...

# (2)

# Nejprve nastavime obrazek na 1x3
par(mfrow = c(1,3))

# Chceme diagramy promennych podle populaci. Tedy vlevo budou promenne, 
# za vlnkou populace:
boxplot(data$BPL ~ data$..., ylab = "delka oblicejove casti lebky (mm)")

# Dalsi boxploty vytvorte analogicky:
boxplot(data$NPH ~ data$..., ylab = "vyska horni casti oblicejoveho skeletu (mm)")
boxplot(... ~ ..., ylab = "vyska ocnice leve strany (mm)")

# (3)

# Pro overeni vicerozmerne normality pouzijeme Mardiuv test. Stejne jako v minulem cviceni, 
# nacteme knihovnu MVN a vyuzijeme funkci mvn(). Testujeme pro mnohorozmernou normalitu 
# pro jednotlive populace, proto subset nastavime na "Population", mvnTest bude roven "mardia" 
# a pro vykresleni QQ-plotu dame multivariatePlot = "qq"
library(...)

# Nastavime obrazek na 2x2
par(mfrow = c(2,2))
mvn(data, subset = "...", mvnTest = "...", multivariatePlot = "...")$multivariateNormality

# U kterych populaci tento test mnohorozmernou normalitu zamitl?

# I presto, ze testy u nekterych populaci normalitu zamitly, budeme pokracovat dal.

# Obdobne muzeme proves i Henzeuv-Zirkleruv test. Jen nastavime mvnTest na "hz". Tentokrat
# nebudeme vykreslovat QQ-ploty.
par(mfrow = c(2,2))
mvn(data, subset = "...", mvnTest = "...")$multivariateNormality

# Nebo muzeme pouzit i funkci HZ.test(), ktera vsak testuje populace jednotlive
library(mvnTest)
HZ.test(data[data$Population == "...",-1])
HZ.test(data[data$Population == "...",-1])
HZ.test(...)
HZ.test(...)

# (4)

# Pro overeni shodnosti variancnich matic (coz je dalsi predpoklad), 
# pouzijeme Boxuv M-test z knihovny biotools. Seskupujeme opet podle populaci.
library(...)
boxM(data[,-1], grouping = data$...)

# Interpretujte vysledky testu:
# ...

# (5)

# Sestavime model mnohorozmerne analyzy rozptylu a otestujeme v R pomoci 
# ruznych testu hypotezu o tom, ze vektory strednich hodnot (promennych BPL, NPH a OBH)
# jsou ve vsech populacich stejne. Vyuzijeme funkci manova(). Do formule zadame jako 
# zavislou promennou nase data bez prvniho sloupce (protoze ten neni numericky) a  
# zaroven musime nastavit, abychom tato data uvazovali jako matici, proto pouzijme
# funkci as.matrix(). Jako nezavislou promennou ve formuli (tj. tu vpravo) nastavime
# populaci.
model <- manova(...(data[,...]) ~ data$...)

# Pomoci funkce summary() provedeme postupne vsechny pozadovane testy.
summary(..., test = "Wilks")
summary(...,test = "Pillai")
summary(...,test = "Hotelling-Lawley")
summary(...,test = "Roy")

# Ktere testy zamitaji hypotezu o shode vektoru strednich hodnot?

# (6)

# Zamitli jsme shodu vektoru strednich hodnot. Nyni chceme zjistit ve kterych 
# konkretnich promennych se populace lisi.

# Vypocteme matici meziskupinove variability
SSB <- summary(model, test = "Wilks")$SS[[1]]

# Vypocteme matici rezidualni variability
SSE <- summary(model, test = "Wilks")$SS[[2]]

# A jejich souctem ziskame matici celkove variability
SST <- ... + ...

# Vypocteme testove statistiky pro simultanni testy (vzorecek viz prednaska)
n <- nrow(data) # pocet pozorovani
k <- ... # pocet promennych
r <- ... # pocet skupin
const <- (n- (k+r)/2 -1)
K <- -const * log(diag(SSE)/diag(SST))
K

# Melo by vam vyjit:
# BPL NPH OBH
# 77.25455 29.20098 78.23122 

# Vypocteme prislusny kvantil pro stanoveni kritickeho oboru
kvantil <- qchisq(0.95, df = k*(r-1))
kvantil

# Ktere promenne zpusobuji rozdily mezi populacemi (skupinami)?
# ...

# (7)

# Nyni chceme zjistit, ktere dvojice populaci se lisi. K tomu 
# vyuzijeme Hotellingovy testy pro dvojice populaci, avsak musime upravit hladinu 
# vyznamnosti. 

# Puvodni alpha vydelime kombinacnim cislem "pocet populaci nad 2". Kombinacni 
# cislo v R vytvorime pomoci choose()
alpha.korig <- 0.05 / ...(r,2)
alpha.korig

# Nejprve nacteme knihovnu ICSNP, nasledne pro porovnani nasich populaci
# (vyber provedeme pomoci hranatych zavorek) pouzijeme funkci HotellingsT2().
# Mame 4 pop., celkem tedy musime provest 6 porovnani.
library(...)
HotellingsT2(data[data$... == "ZALAVAR",-1], data[data$... == "BERG",-1])
HotellingsT2(data[data$Population == "ZALAVAR",-1], data[data$Population == "ESKIMO",-1])
HotellingsT2(data[data$Population == "ZALAVAR",-1], data[data$Population == "NORSE",-1])
HotellingsT2(data[...], data[...]) # BERG vs ESKIMO
HotellingsT2(data[...], data[...]) # BERG vs NORSE
HotellingsT2(data[...], data[...]) # ESKIMO vs NORSE

# Ktere dvojice populaci se lisi (p-hodnotu porovnavejte s alpha.korig)?
# ...

# (8)

# Nyni chceme zjistit, ktere promenne zpusobuji rozdily mezi jednotlivymi
# dvojicemi skupin. Opet musime korigovat hladinu vyznamnosti
# (Techto testu je v nasem pripade 18, tedy budeme alpha delit 18).

alpha.korig2 <- 0.05/(k*r*(r-1)/2)
alpha.korig2

# Jelikoz budeme overovat pomoci dvouvyberovych t-testu, pro vsech techto 18 testu 
# i predpoklad rovnosti variancnich matic pomoci funkce var.test()

# Postupujeme promennou po promenne, proto musime vzdy vybrat jen jeden sloupec 
# (zde testujeme na hl. vyzn 0.05)
var.test(data[data$Population == "ZALAVAR",2], data[data$Population == "BERG",2]) # BPL
var.test(data[data$Population == "ZALAVAR",3], data[data$Population == "BERG",3]) # NPH
var.test(data[data$Population == "ZALAVAR",4], data[data$Population == "BERG",4]) # OBH
var.test(data[data$Population == "ZALAVAR",2], data[data$Population == "ESKIMO",2])
var.test(data[data$Population == "ZALAVAR",3], data[data$Population == "ESKIMO",3]) 
var.test(data[data$Population == "ZALAVAR",4], data[data$Population == "ESKIMO",4]) 
var.test(data[data$Population == "ZALAVAR",2], data[data$Population == "NORSE",2])
var.test(data[data$Population == "ZALAVAR",3], data[data$Population == "NORSE",3]) 
var.test(data[data$Population == "ZALAVAR",4], data[data$Population == "NORSE",4]) 
var.test(data[data$Population == "BERG",2], data[data$Population == "ESKIMO",2])
var.test(data[data$Population == "BERG",3], data[data$Population == "ESKIMO",3]) 
var.test(data[data$Population == "BERG",4], data[data$Population == "ESKIMO",4]) 
var.test(data[data$Population == "BERG",2], data[data$Population == "NORSE",2])
var.test(data[data$Population == "BERG",3], data[data$Population == "NORSE",3]) 
var.test(data[data$Population == "BERG",4], data[data$Population == "NORSE",4]) 
var.test(data[data$Population == "ESKIMO",2], data[data$Population == "NORSE",2])
var.test(data[data$Population == "ESKIMO",3], data[data$Population == "NORSE",3]) 
var.test(data[data$Population == "ESKIMO",4], data[data$Population == "NORSE",4]) 

# Zadny test nazamitl rovnost rozptylu. V t-testech tedy vyuzijeme variantu
# var.equal = T. Nyni porovnavame s hodnotou alpha.korig2 = 0.002777778.
# Zkuste si sami overit, ktere populace se lisi a ve kterech promenych.

# ZALAVAR vs BERG
t.test(data[data$Population == "ZALAVAR",2], data[data$Population == "BERG",2], var.equal = T) # BPL # zam.
t.test(data[data$Population == "ZALAVAR",3], data[data$Population == "BERG",3], var.equal = T) # NPH
t.test(data[data$Population == "ZALAVAR",4], data[data$Population == "BERG",4], var.equal = T) # OBH

# ZALAVAR vs ESKIMO
t.test(data[data$Population == "ZALAVAR",2], data[data$Population == "ESKIMO",2], var.equal = T) # zam.
t.test(data[data$Population == "ZALAVAR",3], data[data$Population == "ESKIMO",3], var.equal = T)  # zam.
t.test(data[data$Population == "ZALAVAR",4], data[data$Population == "ESKIMO",4], var.equal = T) # zam.

# ZALAVAR vs NORSE
t.test(data[data$Population == "ZALAVAR",2], data[data$Population == "NORSE",2], var.equal = T)
t.test(data[data$Population == "ZALAVAR",3], data[data$Population == "NORSE",3], var.equal = T) 
t.test(data[data$Population == "ZALAVAR",4], data[data$Population == "NORSE",4], var.equal = T)

# BERG vs ESKIMO
t.test(data[data$Population == "BERG",2], data[data$Population == "ESKIMO",2], var.equal = T) # zam.
t.test(data[data$Population == "BERG",3], data[data$Population == "ESKIMO",3], var.equal = T) # zam.
t.test(data[data$Population == "BERG",4], data[data$Population == "ESKIMO",4], var.equal = T) # zam.

# BERG vs NORSE
t.test(data[data$Population == "BERG",2], data[data$Population == "NORSE",2], var.equal = T) # zam.
t.test(data[data$Population == "BERG",3], data[data$Population == "NORSE",3], var.equal = T) 
t.test(data[data$Population == "BERG",4], data[data$Population == "NORSE",4], var.equal = T) 

# ESKIMO vs NORSE - posledni zkuste sami
t.test(...) # zam.
t.test(...) # zam.
t.test(...) # zam.