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

#############################################################################
# 1. cviceni
#############################################################################

# Nejprve musime nacist data. Jednim ze zpusobu je vyuzit 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 "anova-means-skull.txt" otevrete napr. v poznamkovem bloku.

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

# Pomoci funkce summary se podivame, jak nase data vypadaji. Jsou v nich NA hodnoty?
summary(...)

# Odstraneni NA hodnot provedeme pomoci funkce na.omit(). Vysledek ulozime do promenne skulls
skulls <- na.omit(...)

# Aby dalsi prikazy v R spravne fungovaly, musi byt nezavisla promenna typu faktor. Pomoci $ vybereme sloupec pop
# promenne skulls:
is.factor(skulls$...)

# Pro vytvoreni promenne typu faktor muzeme pouzit funkci as.factor()
skulls$pop <- as.factor(...)

# Znovu se podivame na nase data. Vsimnete si rozdilu ve sloupecku "pop":
summary(...)

# Celkovy prumer pro vsechny populace naraz. Vyuzijeme funkci mean
...(skulls$upface.H)

# Pro vypocet prumeru jednotlivych populaci muzeme vyuzit funkce tapply() a mean(). Funkce tapply
# aplikuje danou funkci (v nasem pripade mean) na numerickou promennou na prvnim miste podle faktoru
# na druhem miste. Tj. jinak receno: "Pouzij mean() na skulls$upface.H podle skulls$pop":
(group.means <- tapply(skulls$upface.H, skulls$pop, mean))
## Kulate zavorky kolem zpusobuji okamzity vypis na obrazovku

# Smerodatne odchylky spocteme analogicky pomoci funkce sd(). Zkuste doplnit sami
(group.sd <- tapply(..., ..., ...))

# Pokud bychom chteli pocitat po jedne, museli bychom vyuzit funkci which:
# mean(skulls$upface.H[which(skulls$pop == "ban")])
# a obdobne pro zbyle populace

# Pro spocteni rozsahu jednotlivych populaci je uzitecna funkce table: 
(sizes <- ...(skulls$pop))

# Zakladni charakteristiky pro cely datovy soubor
mean(skulls$upface.H)
sd(skulls$upface.H)
length(skulls$upface.H)

# Nastaveni popisku pro graf:
ylabel <- "Vyska horniho obliceje [mm]"

## Pro vykreslovani grafu naraz si vytvorime vektor obsahujici jen jmena populaci. Vyuzijeme funkci levels:
populations <- ...(skulls$pop)

# Nastvaeni parametru grafu (chceme jedno okno 1 x 1)
par(mfrow=c(1,1))

# Pro vykresleni boxplotu pouzijme funkci boxplot(). V ni zadame formuli typu "zavisla promenna" ~ "nezavisla promenna"
# Dalsi argumenty jsou popisky os, to jestli se ma sirka jednotlivych krabic ridit proporcne poctem pozorovani
# a take, zda maji mit krabice zarezy. 
boxplot(skulls$... ~ skulls$..., ylab=ylabel, xlab="Populace", varwidth=T, notch=T)

# Body pridame pomoci funkce points(). Prvnim argumentem (x-ova souradnice) je jen ciselny vektor udavajici cislo populace, druhym
# argumentem (y-ova souradnice) jsou prumery jednotlivych populaci. Nasleduje typ bodu (16 je plne male kolecko) a barva.
points(1:length(levels(skulls$pop)), group.means, pch=..., col="...")
## boxplot populace automaticky abecedne seradi

# Jake jsou predpoklady pro analyzu rozptylu?
# 1 ...
# 2 ...
# 3 ...

# Zacneme overenim predpokladu normality

# Chceme obrazek 2x3
par(mfrow=c(..., ...))

# Postupne budeme pridavat QQ ploty pro jednotlive populace. Pomoci funkce which() a hranatych zavorek 
# vybereme vzdy prislusnou populaci (vybirame z "ban" "cin" "mal" "nem" "per").
# Zbyle QQ ploty zkuste vykreslit sami
qqnorm(skulls$upface.H[which(skulls$pop=="...")], main="", 
       sub=paste("Populace ", "...", sep=""), xlab="Teoreticke kvantily",
       ylab="Vyberove kvantily")
qqline(skulls$upface.H[which(skulls$pop=="...")])


# Pro zjednoduseni prace muzeme vyuzit cyklus. Promenna j nabyva hodnot od 1 do 5 a my v kazdem kroku cyklu 
# vybereme prislnou populaci z vektoru populations. Vsimnete si, ze napr. populations[1] je "ban"
for (j in 1:length(populations))
{
  qqnorm(skulls$upface.H[which(skulls$pop==populations[j])], main="", 
         sub=paste("Populace ",populations[j], sep=""), xlab="Teoreticke kvantily",
         ylab="Vyberove kvantily")
  qqline(skulls$upface.H[which(skulls$pop==populations[j])])
  
}

# Overeni normality pomoci Shapiro-Wilkova testu
shapiro.test(skulls$upface.H[which(skulls$pop=="...")])
shapiro.test(skulls$upface.H[which(skulls$pop=="...")])
shapiro.test(skulls$upface.H[which(skulls$pop=="...")])
shapiro.test(skulls$upface.H[which(skulls$pop=="...")])
shapiro.test(skulls$upface.H[which(skulls$pop=="...")])
# Jak testy dopadly?

# Muzeme rovnez vyuzit Lilieforsovu variantu Kolmogorovova – Smirnovova testu
library(nortest)
lillie.test(skulls$upface.H[which(skulls$pop=="...")])
lillie.test(skulls$upface.H[which(skulls$pop=="...")])
lillie.test(skulls$upface.H[which(skulls$pop=="...")])
lillie.test(skulls$upface.H[which(skulls$pop=="...")])
lillie.test(skulls$upface.H[which(skulls$pop=="...")])

ad.test(skulls$upface.H[which(skulls$pop=="ban")])
ad.test(skulls$upface.H[which(skulls$pop=="cin")])
ad.test(skulls$upface.H[which(skulls$pop=="mal")])
ad.test(skulls$upface.H[which(skulls$pop=="nem")])
ad.test(skulls$upface.H[which(skulls$pop=="per")])

# Normalita je hranicni, u nekterych skupin je zamitnuta.
# Anova vsak neni prilis citliva na poruseni normality, proto budeme postupovat dal.
# Jinak bychom museli pouzit Kruskal-Wallisuv test (neparametricka obdoba anova).

# Overeni predpokladu homogenity rozptylu (nutne pro ANOVA)

# Vyuzijeme Bartlettuv test. Argumentem je formule typu "zavisla promenna" ~ "nezavisla promenna" 
# stejne jako u vykreslovani boxplotu
bartlett.test(..., data=skulls)

# Muzeme pouzit i Levenuv test
library("car")
leveneTest(..., data=skulls)
## nezamitame homogenitu rozptylu

# A nyni provedem samotnou analyzu rozptylu - ziskame anova tabulku. Opet provedeme zadani pomoci formule.
m.skull.1 <- aov(..., data=skulls)
anova(...)
## nezamitame hypotezu o shode strednich hodnot

# Vyzkousejte si i jine zpusoby ziskani anova tabulky v R.

# Tukeyho metoda mnohonasobneho porovnavani - neni treba delat, protoze jsme hypozetu o shode strednich hodnot nezamitli
TukeyHSD(aov(upface.H ~ pop, data=skulls))