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

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

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

# Nyni ziskame zakladni prehled o datech
summary(fat)

# Vykreslime si nase data do grafu, kde na ose x bude tloustka kozni rasy
# a na ose y BMI index. 

# Pomoci notace s $ se odkazujeme na jednotlive slozky nasich dat. Doplnte 
# nazvy jednotlivych promennych ve spravnem poradi (nejprve zadavame osu x, pote y)
plot(fat$..., fat$..., xlab = "Tloustka kozni rasy (mm)", ylab = "Index BMI")

# Nyni sestrojime linearni model. To se v R-ku provede prikazem lm(). Dulezite je zadat
# spravne formuli, kde vlevo od "~" je zavisla promenna a vpravo od "~" je nezavisla 
# promenna. Nezapomente pridat i z jakeho datoveho souboru nase promenne pochazeji.

m.weight <- lm(... ~ ..., data = ...)

# Abychom vsak mohli pouzit linearni model, musi byt splneny ruzne predpoklady.
# Ktere to jsou? Sami si doplnte:
# 1) ...
# 2) ...
# 3) ...
# 4) ...

# Nastavime, aby se grafy vykresloval ve formatu 2x2
par(mfrow=c(2,2))

# A vykreslime diagnosticke grafy. Jak by kazdy graf mel vypadat v idealnim pripade?
# ...
plot(m.weight)

# Nyni budeme testovat normalitu rezidui. Promenna m.weight, kterou jsme vytvorili,
# v sobe obsahuje spoustu dalsich udaju. Na kazdy z nich se opet dostaneme pomoci $.
# Rezidua najdeme pod "residuals"
shapiro.test(m.weight$...)

# Interpretujte vysledky testu. 
# ...

# Normalitu muzeme testovat i jinymy zpusoby. Muzeme pouzit napriklad Lilieforsovu 
# variantu Kolmogorovova – Smirnovova testu. K tomu vask nejprve musime nacist balicek 
# nortest.
library(nortest)

# Vyzkousejte tento test na rezidua naseho modelu.
lillie.test(...)

# Interpretujte vysledky tohoto testu.
# ...

# I kdyz je normalita rezidui na povazenou, tak budeme pokracovat.

# Pomoci t-testu overime, zda reziua nemaji systematickou chybu, tj. zda jejich 
# stredni hodnota je rovna nule.
t.test(...)

# Interpretujte vysledky tohoto testu.
# ...

# Pro Durbin – Watsonuv test potrebujeme nacist knihovnu car.
library(car)

# Overime nezavislost rezidui.
durbinWatsonTest(m.weight)

# Interpretujte vysledky tohoto testu.
# ...

# Az na mirne poruseni normality jsou predpoklady pro linerani model splneny.

# Odhady parametru, index determinace, celkovy F-test i jednotlive t-testy - to vse
# najdeme ve vystupni tabulce ze summary.

# Aplikujte summary na nas model a najdete prislusne hodnoty
summary(...)

# Index determinace udava, jakou cast variability zavisle promenne (BMI v nasem pripade)
# lze vysvetlit zvolenou regresni funkci.
# Index determinace (= koeficient determinace): ...

# Odhady parametru jsou 
# Beta0 = 17.10928 - stredni hodnota predikce pri nulovych hodnotach vsech ostatnich regresoru
# tj. index BMI pri nulove tloustce kozni rasy
# Beta1 = ... - rozdil v predikci, kdyz se tloustka kozni rasy zvedne o 1 mm

# Celkovy F-test: testuje, ze nas model je lepsi nez nahoda, tj. testuje 
# nulovou hypotezu o tom ze vsechny koeficienty naseho modelu jsou nulove.
# F-statistika: ..., p-hodnota: ...

# Doplnte p-hodnotu pro jednotlive koeficienty Beta0 a Beta1. Co jimi testujeme?
# ...

# Nyni vypocitame intervaly spolehlivosti pro koeficienty naseho modelu
confint(...)

# Nakonec budeme vykreslovat grafy. 

# Aby zobrazeni nasi regresni primky, IS a predikcnich intervalu bylo 
# dostatecne hladke, vytvorime v R dostatecne hustou sekvenci bodu na ose x.

# Tato sekvence bude zacinat v minimalni hodnote kozni rasy uvedene v nasich datech
# a koncit bude v maximalni hodnote. Delku sekvence nastavte na 300.
xx <- seq(...(...), max(fat$hip.F), length = ...)

# Nastavime vykreslovani naseho grafu pouze na jeden obrazek, tj. 1x1
par(mfrow = c(1,1))

# Intervaly spolehlivosti vypocteme pomoci funkce predict. Ta se vzdy aplikuje na nas model
# a do argumentu newdata uvedeme ty hodnoty ve kterych chceme dopocitat nasi predikci.
# V nasem pripade chceme predikovat v hodnotach xx. Typ intervalu nastavime na "confidence".
interval.spol <- predict(..., newdata=data.frame(hip.F = ...), interval = "confidence")

# Predikci interval vypocteme analogicky s tim rozdilem, ze interval nastavime na "predict".
pred.interval <- predict(..., newdata=data.frame(... = xx), interval = "...")

# Opet vykreslime nejprve bodovy graf (=teckovy diagram). Az posleze do nej budeme
# prikreslovat dalsi primky.
plot(..., ..., xlab = "Tloustka kozni rasy (mm)", ylab = "Index BMI")

# Prikazem lines doplnime do obrazku primky. Nejprve doplnime samotnou regresni primku.
# Hodnoty y teto primky najdeme v prvnim sloupecku promenne "interval.spol". Barvu primky
# nastavte na cervenou.

lines(xx, ...[,1],col = ...)

# Analogickym zpusobem doplnime IS. Dolni hranici IS a hornici IS najdeme ve druhem, resp. 
# tretim sloupecku promenne "interval.spol". Typ primky nastavte na carkovanou, tj. =2.
lines(xx,interval.spol[,2], col = "red", lty = ...)
lines(xx,interval.spol[,3], col = "red", lty = ...)

# Predikcni intervaly vykreslime stejnym zpusobem, jako intervaly spolehlivosti. 
# Mame je akorat ulozeny v promenne "pred.interval". Jejich barvu nastavte na modrou.
lines(xx, ...[,2], col = ..., lty = 2)
lines(xx, ...[,3], col = ..., lty = 2)

# Nyni jen doplnime legendu do grafu
legend("topleft",c("model","IS","pred. int."), lty = c(1,2,2),
       col = c("red", "red", "blue"))

# Priklad 2
# Postupujeme analogicky jako v prikladu 1

foot <- read.table("Data/lrm-foot.txt", header = ...)

summary(...)

plot(foot$..., foot$..., xlab = "Delka chodidla (mm)", ylab = "Telesna vyska (mm)")

m.foot <- lm(... ~ ..., data = ...)

# Diagnosticke grafy
par(mfrow = c(2,2))
plot(m.foot)
# Jejich interpretace: ...

# Testy normality
shapiro.test(...)

library(nortest)
lillie.test(...)

# Interpretujte vysledky techto testu
# ...

# K cemu slouzi tento test?
t.test(...)



library(car)
durbinWatsonTest(m.foot)

# Jsou rezidua nezavisla?

summary(m.foot)

# Interpretace odhadu parametru, indexu determinace a jejich hodnoty:
# Beta0: ...
# Beta1: ...
# ID: ...

# Celkovy F-test: Interpretace, hodnota F-statistiky a p-hodnota: ...
# Jednotlive t-testy: Interpretace a p-hodnoty: ...

# Intervaly spolehlivosti pro koeficienty beta
confint(...)

# Vykreslovani grafu:

xx <- seq(min(...), max(...), length = ...)

par(mfrow = c(1,1))

# Intervaly spolehlivosti
interval.spol <- predict(m.foot, newdata = data.frame(foot.L = xx), interval = "confidence")

# Predikcni intervaly
pred.interval <- predict(..., newdata=data.frame(... = ...), interval = ...)

# Teckovy diagram
plot(..., ..., xlab = "Delka chodidla (mm)", ylab = "Telesna vyska (mm)")

# Vykresleni regresni primky
lines(..., ...[,1],col = "red")

# Vykresleni IS a predikcnich intervalu
lines(xx, ...[,2], col = "red", lty = 2)
lines(xx,interval.spol[,3], col = "red", lty = 2)

lines(..., ...[,2], col="blue", lty = 2)
lines(..., ...[,3], col="blue", lty = 2)

# Legenda
legend("topleft",c("model","IS","pred. int."), lty = c(1,2,2),
       col=c("red", "red", "blue"))