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

#############################################################################
# 9. 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 "lrm-foot.txt" otevrete napr. 
# v poznamkovem bloku.

# Doplnete (v uvozovkach!) cestu k datovemu souboru a hodnotu T/F, podle toho, 
# jestli soubor ma hlavicku

# 1

data.cela <- read.table("...", header = ...)
summary(data.cela)

# Chceme vzit jen zeny, proto vyuzijeme vyberu pomoci hranatych zavorek. 
# Bereme jen ty radky, kde pohlavi je rovno hodnote "f". Zaroven vybereme
# jen sloupce 2 a 3 (prvni sloupec udavajici pohlavi uz nebudeme potrebovat).
data <- data.cela[data.cela$... == "f", ...:...]

# K vypoctu prumeru pro "foot.L" a "body.H" vyuzijeme funkci colMeans() 
# aplikovanou na data
colMeans(...)

# Pro varianci matici vyuzijeme funkci var()
var(...)

# Hotellinguv test

# Nejprve musime overit predpoklad. Jaky to je?
# ...

# K tomu ucelu muzeme vyuzit napriklad Henzeuv-Zirkleruv test. Bud nacteme 
# knihovnu mvnTest a nasledne pouzijeme funkci HZ.test() na data:
library(...)
HZ.test(...)

# Anebo vyuzijme jeden z vystupu funkce mvn() z knihovny MVN. Argument
# "mvnTest" nastvaime na "hz" a multivariatePlot na "qq" pro soucasne vykresleni 
# QQ-plotu. Pozadovany vystup z funkce mvn() najdeme pod $multivariateNormality.
library(...)
mvn(data, mvnTest = "...", multivariatePlot = "...")$...

# (Pozn.: Obe funkce provadi tentyz test, ale vysledky jsou kvuli jinemu postupu 
# vypoctu malinko odlisne.)

# Jaky je vysledek HZ-testu?
# ...

# Pomoci jednovyberoveho Hotellingova testu otestujeme hypotezu, ze vektor
# strednich hodnot je roven zadanemu vektoru.
mu0 <- c(250, 1679)

# Nacteme knihovnu ICSNP a provedeme samtony test. Aplikujeme jej na data,
# vektor strednich hodnot nastavime mu nastavime na nas zadany vedktor "mu0".
library(ICSNP)
HotellingsT2(data, mu = ...)

# Jaky je vysledek? Zamitame, nebo nezamitame hypotezu o rovnosti vektoru
# strednich hodnot?
# ...

# Jelikoz jsme nulovou hypotezu ..., musime provest jednorozmerne t-testy.

# Nejprve spocteme Bonferroniho korigovanou hladinu vyznamnosti: Pro jednotlive
# testy podelime hladinu vyznamnosti poctem testovanych hypotez, tj.
alpha.korig <- 0.05 / 2

# Nyni pomoci funkce t.test() provedeme samotne testy. Argumentem je vektor
# danych pozorovani a mu, ktere udava stredni hodnotu, jez chceme testovat.
t.test(data$..., mu = mu0[...])
t.test(data$..., mu = mu0[...])
# Kvuli ktere slozky doslo k zamitnuti vicerozmerne hypotezy? 
# ...

# Priklad 2

# Nacteni dat - pracujeme se souborem skulls.txt
dataS <- read.table("...", header = ...)
summary(dataS)

# Chceme vypocitat prumery jen pro zeny. K tomu opet vyuzijeme funkci
# colMeans a hranatych zavorek. Analogicky jako v prvnim prikladu vybereme
# jen zeny a sloupce obsahujici numericke promenne (druhy az paty).
colMeans(dataS[dataS$... == "f", ...])

# Pro vypocet vyberove varinacni matice pouzijeme funkci var().
var(dataS[dataS$... == "f", ...])

# stejnym zpusobem pro muze:
colMeans(...)
var(...)

# Dvouvyberovy Hotellinguv test

# Nejprve overime jeho predpoklady. Ktere to jsou?
# ...

# Pro vicerozmernou ... muzeme opet vyuzit napriklad H-Z test. (Mohli bychom pouzit 
# i Mardiuv test - viz pracovni list)

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

# Pouziti funkce mvn() je stejne jako v pr. 1, jen navic zadame argument
# subset, ktery udava podle ktere prommene rozlisujeme datovy soubor.
# V nasem pripade to je prom. "sex".
mvn(dataS, subset = "...", mvnTest = "...", multivariatePlot = "...")$...

# Pokud bychom pouzili funkci HZ.test(), museli bychom kazdou skupinu testovat zvlast:
HZ.test(dataS[... == "f", 2:5], qqplot = T)
HZ.test(dataS[... == "m", 2:5], qqplot = T)

# Dalsim predpokladem je ...

# Njeprve si nacteme knihovnu biotools a nasledne vyuzijeme funkci boxM() na nas datovy
# soubor (jen na numericke promenne). Argument grouping nastavime na faktor udavajici pohlavi.
library(...)
boxM(...[,2:5], grouping = dataS$...)
## Lze na hladine vyznamnosti 5 % zamitnout nulovou hypotezu o shode variancnich matic?
## ...

# Pristoupime tedy k dvouvyberovemu Hotellingovu testu. Opet vyuzijeme funkci
# HotellingsT2(). Jeji argumenty budou datove tabulky - jedna pro zeny, jedna pro muze.
HotellingsT2(dataS[dataS$..., 2:5], ...[dataS$sex == "m", 2:5])

# Jaky je vysledek testu?
# ...

# Nyni pristoupime k simultanim testum a zjistime tak, ve ktre promenne se muzi
# a zeny lisi.

# Musime postupovat podle vzorecku z prednasky. Nejprve proto vypocteme
# n1 a n2 (pocty zen a muzu). Vyuzijeme funkce table() a vyberu z vektoru 
# pomoci hranatych zavorek.
n1 <- table(dataS$sex)[1] 
n2 <- table(dataS$sex)[2]

# Celkovy pocet jedincu
n <- n1 + n2

# Pocet numerickych promennych v tomto prikladu je k = 4.
k <- 4 

# Vypocteme si vektory strednich hodnot pro muze a zeny (stejne 
# jako na zacatku tohoto prikladu)
mu1 <- colMeans(...)
mu2 <- colMeans(...)

# Z vyberovych variancnich matic ziskame jejich diagonaly (funkce diag()).
# Vime, ze na diagonale se nachazi vyberove rozptyly danych promennych.
var1 <- ...(cov(dataS[dataS$sex == "f",2:5]))
var2 <- ...(cov(dataS[dataS$sex == "m",2:5]))

# Dale postupujeme podle vzorecku
var <- ( (n1-1)*var1 + (n2-1)*var2 )/(n-2)
F.stat <- n1*n2*(n-k-1) * (mu1-mu2)^2 /(var*n*k*(n-2))

# Vypocteme p-hodnotu pomoci funkce pf(). Argumentem je hodnota 
# F-statistiky a dane stupne volnosti ("k" a "n-k-1").
p.hodnota <- 1-...(F.stat, k, n-k-1)

# Prislusny kvantil vypocteme pomoci funkce qf()
kvantil <- ...(0.95, k, n-k-1)

# Vsechny hodnoty umistime do tabulky a vypiseme si ji (volime 
# zaokrouhleni na 4 des. mista).
tab <- round(rbind(F.stat,p.hodnota, kvantil), digits = ...)
rownames(tab) <- c("F","p-hodnota", "kvantil")
tab

# Ve kterych promennech se zeny a muzi lisi? Melo by vam vyjit, ze ve vsech :)