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

#############################################################################
# 6. 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 "cneck.txt" otevrete napr. v poznamkovem bloku.

# Doplnete (v uvozovkach!) cestu k datovemu souboru a hodnotu T/F, podle toho, 
# jestli soubor ma hlavicku
neck <- read.table("...", header = ...)

# Podivejte se na zakladni informace o nasem datovem souboru
summary(...)

# Zacneme vykreslenim bodovych diagramu pro muze a zeny do jednoho obrazku.
# Nejprve vykreslime prazdne okno grafu pomoci funkce plot a nastavenim argumentu
# type = "n". Na ose x (nezavisla promenna) mame obvod boku, na ose y telesnou vahu.
# Doplnte nazvy promennych a typ vykreslovaneho grafu.
plot(neck$..., neck$..., xlab = "Obvod boku (mm)", ylab = "Telesna hmotnost (kg)",
     type = "...")

# Pridame body vyjadrujici jednotliva pozorovani. Nejprve vykreslime zeny cervenou barvou
# a pote muze modrou barvou. Vyber (zen/muzu) provedeme v R pomoci hranatych zavorek.
# Nyni tedy vybirame z vektoru neck$hip.C nebo neck$body.W jen ta pozorovani, pro ktera 
# plati ze promenna sex je rovna hodnote "f".
points(neck$hip.C[neck$sex == "..."], neck$body.W[neck$sex == "..."], pch = 1, col = "red")

# A obdobne pro muze
points(neck$hip.C[...], neck$body.W[...], pch = 2, col = "blue")

# Do obrazku pridame i legendu. 
legend("topleft", c("zena", "muz"), pch = c(1,2), col = c("red", "blue"))

# Vsimnete si dobre problemovych pozorovani.

# Nyni mame vypocitat pocty pozorovani, prumery a smerodatne odchylky.

# Pocty muzu a zen jsme videli uz v tabulce summary(). Muzeme vsak vyuzit i funkci 
# table(), kterou pouzijeme na promennou pohlavi.
table(neck$...)

# Pocet radku cele tabulky udava celkovy pocet pozorovani
nrow(...)

# Nejprve se budeme venovat zavisle promenne - body.W
# Vypocteme prumer (funkce mean()) a smerodatnou odchylku jen pro 
# zeny. Opet musime pouzit zapis pres hranate zavorky. 
mean(neck$body.W[neck$sex == "..."])
sd(neck$body.W[neck$sex == "..."])

# Pro muze zkuste napsat prikazy sami
mean(...)
sd(...)

# Jak bude vypadat vypocet prumeru a smerodatne odchylky promenne body.W 
# pro cely datovy soubor?
mean(...)
sd(...)

# Zcela analogicky budeme postupovat pro promennou hip.C.
# Pro zeny:
mean(neck$hip.C[neck$sex == "..."])
sd(neck$hip.C[neck$sex == "..."])

# Pro muze
mean(...)
sd(...)

# Pro cely datovy soubor
mean()
sd()

# 3 
# Mame overit predpoklad rovnobeznosti primek. To se v R provede nasledujicim zpusobem -
# vytvorime modely kdy promennou body.W vysvetlime pomoci promennych sex a hip.C. 
# V jednom modelu vsak uvazujeme interakce mezi promennymi, zatimco v druhem modelu nikoli.
# Tzn. model bez interakce (zadavame v R do formule pomoci +) uvazuje obe primky rovnobezne,
# ale model s interakci (zadavame v R do formule pomoci *) pripousti primky ruznobezne.
# Nasledne tyto modely porovneme - pokud nam vyjde, ze se nelisi (p-hodnota vyssi nez 0,05), 
# tak je zbytecne uvazovat slozitejsi model s interakci a vezmeme model bez interakce.
# Pokud vyjde, ze slozitejsi model vysvetluje data lepe (p-hodnota nizsi nez 0,05), vezmeme ten.

# Napiste formule obou modelu
model.interakce <- lm(... ~ ... * ..., data = ...)
model.bez.int <- lm(... ~ ... + ..., data = ...)

# A oba modely porovname. 
anova(model.bez.int, model.interakce)

# 4 
# Nyni overime zbyvajici predpoklady analyzy kovariance. Ktere to jsou?
# ...
# ...
# ...
# ...

# Zacneme diagnosticky grafy. Ty uz znate z minulych cviceni. Ve strucnosti:
# 1. graf - posuzuje vhodnost linearniho modelu 
# 2. graf - Q-Q plot posuzuje normalitu 
# 3. graf - posuzuje homoskedasticitu rozptylu 
# 4. graf - detekce vlivnych pozorovani
par(mfrow = c(2,2))
plot(model.bez.int)

# Zamerte se zejmena na graf 4. Ktera pozorovani nam model nejvice "kazi"?

# Znovu vykreslime graf z ukolu 1

par(mfrow = c(1,1))
plot(neck$..., neck$..., xlab = "Obvod boku (mm)", ylab = "Telesna hmotnost (kg)", type = "n")

# Nyni odstranime podezrela pozorovani 12, 23, 36. Z minula jiz vite, jak se to dela.
neck1 <- neck[-c(...),]
## Mame vytvorenou tabulku neck1. Nyni uz budeme pracovat jen s ni.

# Vykrelime zeny, muze a problemova pozorovani
points(neck1$hip.C[neck1$sex == "..."], neck1$body.W[neck1$sex == "..."], pch = 1, col = "red")
points(neck1$hip.C[...], neck1$body.W[...], pch = 2, col = "blue")

# Problemova pozorovani vykreslujeme z puvodni tabulky. Jsou to prave ony tri radky: 12, 23, 36.
points(neck$hip.C[c(...)], neck$body.W[c(...)], pch = 3, col = "green")

legend("topleft", c("zena", "muz","odlehle"), pch = c(1,2,3), 
       col = c("red", "blue","green"))

# 6
# Opet overime predpoklad rovnobeznosti primek
model.interakce <- lm(... ~ ... * ..., data = neck1)
model.bez.int <- lm(... ~ ... + ..., data = neck1)
anova(model.bez.int, model.interakce)

# Vytvorime diagnosticke grafy? Je viditelne zlepseni?
par(mfrow = c(2,2))
plot(model.bez.int)

# A overime i zbyvajici predpoklady. Ujistete se, ze vite, co kterym testem overujeme.
shapiro.test(model.bez.int$residuals)
t.test(model.bez.int$residuals)

library(car)
durbinWatsonTest(model.bez.int)

# 7
# Nyni mame overit shodnost regresnich primek. To se v R provede jednoduse
anova(model.bez.int)

# Vidime, ze obe promenne (sex i hip.C) jsou vyznamne, nebudeme tedy nas model dale 
# zjednodusovat. Na radku u faktorove promenne sex provadime test o nevyznamnosti faktoru.
# Kolik vysla testova statistika? Kolik prislusna p-hodnota? 
# ...
# Zaroven na radku u promenne hip.C provadime test o nulovosti regrese (viz prednaska).
# Kolik vysla testova statistika? Kolik prislusna p-hodnota? 
# ...



# 8
# Vykresleni pasu spolehlivosti

# Nejprve nastavime sekvenci bodu na ose x (promenna obvod boku)
xx <- seq(min(neck1$...), max(neck1$...), length = 300)

# Vypocteme odhad regresni primky i hranic pasu spolehlivosti. K tomu vyuzijeme funkci predict(),
# kterou aplikujme na nas model. Nejprve pro zeny:
interval.spol.zeny <- predict(...,newdata = data.frame(hip.C = xx, sex = "..."), 
                              interval = "confidence")
# A analogicky pro muze
interval.spol.muzi <- predict(...)

# Nastavime okno grafu na 1x1
par(mfrow = c(1,1))

# Vykreslime prazdny graf
plot(neck1$..., neck1$..., xlab = "Obvod boku (mm)", ylab = "Telesna hmotnost (kg)", 
     type = "...")

# Pridame body. Nejprve zeny
points(neck1$hip.C[neck1$sex == "..."], neck1$body.W[neck1$sex == "..."], pch = 1, col = "red")

# Analogicky pridame i muze
points(..., pch = 2, col = "blue")

# Vykreslime predikci stredni hodnoty telesne hmotnosti pro zeny (skryva se v prvnim sloupci
# promenne interval.spol.zeny)
lines(xx, interval.spol.zeny[,...])

# Pasy spolehlivosti jsou ve sloupcich 2 a 3. Rovnez zmente typ cary na 2.
lines(xx, interval.spol.zeny[,...], lty = ...)
lines(xx, interval.spol.zeny[,...], lty = ...)

# Stejnym zpusobem postupujeme pro muze. Zkuste sami
lines(...)
lines(...)
lines(...)

legend("topleft", c("zena", "muz"), pch = c(1,2), col = c("red", "blue"))

# Pokud bychom chteli znat sklon primky (beta) nasich rovnobezek, muzeme ho videt 
# na radku u hip.C ve sloupci Estimate, pouzijeme-li funkci summary na nas model.
summary(model.bez.int)