# Nacteme datovy soubor pomoci read.csv(),
# na.strings rika, jake hodnoty v datech maji byt brany jako NA - 
# - nastavime na "0".
cranio <- read.csv("Data/Howell.csv", header = T, na.strings = "0")

# vybereme pouze muze danych populaci a prislusne promenne
# pozor na kodovani, v tomto souboru jsou muzi kodovani "M"
data <- cranio[cranio$Sex == "M" & cranio$Population %in% c("AINU", "N JAPAN", "PERU"),
               c("Population", "XFB", "ZYB", "ZMB")]
'TO-DO' # uprava data$Population na faktor
'TO-DO' # prohlednuti dat, prip. odstraneni NA



# (1)

# Pocet pozorovani pro kazdou populaci ziskame pomoci funkce table()
# aplikovane na data$Population nebo zjistime v summary
table(data$Population)


# vyberove prumery a vyberove variancni matice
colMeans(data[data$Population == "AINU", -1])
var(data[data$Population == "AINU", -1])

'TO-DO' # ostatni populace

# (2) krabicove diagramy podle populaci

par(mfrow = c(1, 3))
boxplot(data$XFB ~ data$Population,
        ylab = "max. transversalni sirka cela (mm)", xlab = "")
"TO-DO" # ostatni promenne

# (3) Overeni vicerozmerne normality pro jednotlive populace

library(MVN)
par(mfrow = c(1, 3))
mvn(data, subset = "Population", mvnTest = "mardia", multivariatePlot = "qq")$multivariateNormality
'TO-DO' # Henzeuv-Zirkleruv test

# (4) Overeni shodnosti variancnich matic

biotools::boxM(data[,-1], grouping = data$Population)

# (5) model mnohorozmerne analyzy rozptylu 
# pozor, vstup musi byt matice - pouzijeme as.matrix

model <- manova(as.matrix(data[,-1]) ~ data$Population)

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

# (6) ktere promenne se lisi
summary.aov(model)

# nebo postup podle vzorce a rozhodnuti pomoci kritickeho oboru
SSB <- summary(model, test = "Wilks")$SS[[1]] # meziskupinova variabilita
SSE <- summary(model, test = "Wilks")$SS[[2]] # rezidualni variabilita
SST <- SSB + SSE # celkova variabilita

# testove statistiky pro simultanni testy (vzorec viz prednaska)
n <- "TO-DO" # pocet pozorovani
k <- "TO-DO" # pocet promennych
r <- "TO-DO" # pocet skupin
const <- (n - (k + r) / 2 - 1)
K <- -const * log(diag(SSE) / diag(SST))
K # hodnoty test. statistik pro jednotlive promenne

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

# (7) Hotellingovy testy pro dvojice populaci
# porovnavame populace mezi sebou po dvojicich

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

library(ICSNP)
HotellingsT2(data[data$Population == "AINU",-1], 
             data[data$Population == "N JAPAN",-1])
'TO-DO' # ostatni dvojice

# (8) 

# Nyni chceme zjistit, ktere promenne zpusobuji rozdily mezi jednotlivymi
# dvojicemi populaci. 

# Korigovana hladina vyznamnosti - delime poctem testu
alpha.korig2 <- 0.05 / (k * r * (r - 1) / 2)
alpha.korig2

# overeni predpokladu rovnosti rozptylu pro dvojice populaci


# Postupujeme promennou po promenne, proto musime vzdy vybrat jen jeden sloupec 
# (! zde testujeme na hl. vyzn 0.05, korigovanou pouzijeme az v t.testu)
var.test(data[data$Population == "AINU", 2], data[data$Population == "PERU", 2])
var.test(data[data$Population == "AINU", 3], data[data$Population == "PERU", 3])
var.test(data[data$Population == "AINU", 4], data[data$Population == "PERU", 4])

'TO-DO' # analogicky pro dalsi dvojice populaci a promenne

# t-testy pro porovnani populaci a promennych
# AINU vs PERU
t.test(data[data$Population == "AINU", 2], data[data$Population == "PERU", 2], var.equal = T)
t.test(data[data$Population == "AINU", 3], data[data$Population == "PERU", 3], var.equal = T)
t.test(data[data$Population == "AINU", 4], data[data$Population == "PERU", 4], var.equal = T)
# AINU vs N JAPAN
'TO-DO'

# N JAPAN vs PERU
'TO-DO'



'TO-DO' # nereseny priklad samostatne