# Potrebne knihovny - pokud je dosud nemate, je treba je nejdrive nainstalovat, napr. funkci "install.packages"
library (agricolae)	#	pro Scheffeho test
library (car)	#	pro Levenuv test

#	Nacteme data
tabulka <- read.csv2 (file = "../data/seno.csv")
# Podivame se na strukturu
summary (tabulka)
str (tabulka)
# Zkontrolujeme spravne prirazeni typu "factor" ke kategorialnim promennym
# (Mozna) se nam bude hodit i prevedeni na jeden kombinovany factor
tabulka$group <- factor (paste (tabulka$puda, tabulka$hnojeni, sep = "."))


# Podivame se na cetnosti, prumery a mediany
mean (tabulka$vynos) 
tapply (tabulka$vynos, tabulka$puda, mean)
tapply (tabulka$vynos, tabulka$hnojeni, mean)
tapply (tabulka$vynos, list (tabulka$puda, tabulka$hnojeni), mean)
median (tabulka$vynos) 
tapply (tabulka$vynos, tabulka$puda, median)
tapply (tabulka$vynos, tabulka$hnojeni, median)
tapply (tabulka$vynos, list (tabulka$puda, tabulka$hnojeni), median)

#	Grafy 
plot.design (vynos ~ puda + hnojeni, data = tabulka, fun = mean)
plot.design (vynos ~ puda + hnojeni, data = tabulka, fun = median)

# Boxplot a overeni homogenity rozptylu
boxplot (vynos ~ puda, data = tabulka, xlab = "typ pudy", ylab = "vynos")
points (tabulka$puda, tabulka$vynos)
boxplot (vynos ~ hnojeni, data = tabulka, xlab = "zpusob hnojeni", ylab = "vynos")
points (tabulka$hnojeni, tabulka$vynos)
boxplot (vynos ~ puda + hnojeni, data = tabulka, xlab = "typ pudy a zpusob hnojeni", ylab = "vynos")
boxplot (vynos ~ group, data = tabulka, xlab = "typ pudy a zpusob hnojeni", ylab = "vynos")
points (tabulka$group, tabulka$vynos)

# Bartlettuv test
bartlett.test (vynos ~ group, data = tabulka)

# Levenuv test
leveneTest (vynos ~ group, data = tabulka)

# Rozdeleni datove tabulky podle jednotlivych hdodnot faktoru, napr. takto: 
skupiny <- lapply (levels (tabulka$group), function (L) {
	return (subset (tabulka, group == L))
	})

# Overeni normality - QQ-ploty
par (mfrow = c (2, 3))
sapply (skupiny, function (x) {
	qqnorm (x$vynos, main = x$group[1])
	qqline (x$vynos)
	})
par (mfrow = c (1, 1))

# Overeni normality - Shapiruv-Wilkeuv a Kolmogorovuv-Smirnovuv test
lapply (skupiny, function (x) {
	shapiro.test (x$vynos)
	})
lapply (skupiny, function (x) {
	X <- x$vynos
	U <- (X - mean (X)) / sd (X)
	ks.test (U, "pnorm")
	})


# Dvojne trideni bez interakci


# ANOVA pomoci funkce "aov"
aov.model <- aov (vynos ~ puda + hnojeni, data = tabulka)
summary (aov.model)
# Odhady efektu, kontrastu a strednich hodnot
aov.model$coefficients
model.tables (aov.model, type = "effects")
model.tables (aov.model, type = "means")

# Mnohonasobne porovnavani

# Tukeyuv test
TukeyTest <- TukeyHSD (aov.model)
TukeyTest
plot (TukeyHSD (aov.model, "puda"), las = 1, col = 2)
plot (TukeyHSD (aov.model, "hnojeni"), las = 1, col = 3)

# Scheffeho testScheffeTest <- scheffe.test (aov.model, "puda")
ScheffeTest.puda <- scheffe.test (aov.model, "puda")
ScheffeTest.puda$groups
ScheffeTest.hnojeni <- scheffe.test (aov.model, "hnojeni")
ScheffeTest.hnojeni$groups


# Dvojne trideni s interakcemi


# ANOVA pomoci funkce "aov"
aov.model <- aov (vynos ~ puda + hnojeni + puda:hnojeni, data = tabulka)
aov.model
summary (aov.model)
# anebo zkracene 
aov.model <- aov (vynos ~ puda * hnojeni, data = tabulka)
aov.model
summary (aov.model)
# Odhady efektu, kontrastu a strednich hodnot
aov.model$coefficients
model.tables (aov.model, type = "effects")
model.tables (aov.model, type = "means")

# Mnohonasobne porovnavani

# Tukeyuv test
TukeyTest <- TukeyHSD (aov.model)
TukeyTest
par (mar = c (5.1, 6.1, 4.1, 2.1))	#	upraveni vnejsich okraju grafu pro zvetseni mista pro popisy na osach
plot (TukeyHSD (aov.model, "puda"), las = 1, col = 2)
plot (TukeyHSD (aov.model, "hnojeni"), las = 1, col = 3)
plot (TukeyHSD (aov.model, "puda:hnojeni"), las = 1, col = 4)

# Scheffeho testScheffeTest <- scheffe.test (aov.model, "puda")
ScheffeTest.puda <- scheffe.test (aov.model, "puda")
ScheffeTest.puda$groups
ScheffeTest.hnojeni <- scheffe.test (aov.model, "hnojeni")
ScheffeTest.hnojeni$groups



# Jednoduche trideni pro spojene faktory


# ANOVA pomoci funkce "aov"
aov.model <- aov (vynos ~ group, data = tabulka)
aov.model
summary (aov.model)
# Odhady efektu, kontrastu a strednich hodnot
aov.model$coefficients
model.tables (aov.model, type = "effects")
model.tables (aov.model, type = "means")

# Mnohonasobne porovnavani

# Tukeyuv test
TukeyTest <- TukeyHSD (aov.model)
TukeyTest
par (mar = c (5.1, 6.1, 4.1, 2.1))	#	upraveni vnejsich okraju grafu pro zvetseni mista pro popisy na osach
plot (TukeyHSD (aov.model, "group"), las = 1)

# Scheffeho testScheffeTest <- scheffe.test (aov.model, "puda")
ScheffeTest.puda <- scheffe.test (aov.model, "group")
ScheffeTest.puda$groups