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

#############################################################################
# 5. cviceni
#############################################################################

# Nejprve musime nacist data. Stejne jako na minulem cviceni 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 "mlrm-fat.txt" otevrete napr. v poznamkovem bloku.

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

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

# Stejne jako na minulem cviceni odstranime 36. pozorovani, protoze jsme zjistili,
# ze je odlehle. To se provede odstranenim 36. radku. Vysledna tabulka by mela mit
# 50 pozorovani.
mfat1 <- mfat[...,]

# Nyni vykrelime teckove diagramy pro dvojice promennych, ktere nas zajimaji. 
# Tyto diagramy mohou castou napovedet o pripadne zavislosti danych promennych.
# Musime vsak vybrat jen ty sloupecky tabulky mfat1, ktere nas zajimaji. Za carku 
# v hranatych zavorkach proto vlozime vektor, ktery tyto sloupecky oznacuje 
# (zajimaji nas sloupce 1, 2, 4, 5, 6, 7). Pokud zadame vektor pomoci negativnich 
# cisel, rikame tim, ktere sloupecky nechceme. Proto zadani
# plot(mfat1[,c(1:2,4:7)]) a plot(mfat1[,-3]) je jedno a to same.
plot(mfat1[,c(...)])

# Pro dvojice vybranych promennych spocitame korelacni koeficienty. Pocitame 
# pearsonuv korelacni koef., proto doplnete prislusnou metodu do uvozovek.
cor(mfat1[,-3], method = "...")

# Nyni budeme chtit otestovat nulovost techto koeficientu. Na to potrebujeme funkci
# rcorr() z balicku Hmisc. 

# Mozna budete muset balicek Hmisc nejprve nainstalovat. To se provede kliknutim 
# na zalozku "Packages" vpravo, nasledne na "Install". Do policka napisete Hmisc
# a potvrdite.

# Nacteme balicek Hmisc
library(...)

# Funkce rcorr se musi pouzit pouzit na datovy typ "matice", proto z nasi tabulky
# vytovrime matici pomoci funkce as.matrix(). Argument type je opet "pearson".
rcorr(...(mfat1[,...]), type = "...")
## Interpretujte vysledky testu (nemusite pro vsechny, ujistete se jen, ze vite, 
## co znamena p-hodonta mensi nez 0.05 a vetsi nez 0.05)
# ...


# Nyni budeme pocitat parcialni korelacni koeficienty. To lze provest dvema zpusoby.
# Bud manualne po jednom, nebo naraz pomoci funkce pcor(). Nejprve manualne:

# Budeme pocitat parcialni korelacni koeficient mezi telesnou vahou (body.W) a telesnou
# vyskou (body.H). Znamena to, ze nas zajima vztah techto dvou promenneych pri vylouceni 
# vlivu vsech ostatnich promennych (viz prednaska).

# Nejprve vytvorime model linearni zavislosti telesne vahy na ostatnich prom. (m1)
# a model linearni zavislosti telesne vysky na ostatnich prom. (m2). Doplnete sami 
# vsechny 4 ostatni promenne a datovy soubor.
m1 <- lm(body.W ~ ... + ... + ... + ..., data = ...)
m2 <- lm(body.H ~ ... + ... + ... + ..., data = ...)

# Dale vezmeme rezidua techto modelu. Z drivejska uz vite, jak jsou v R ulozena.
y.res <- m1$...
z.res <- m2$...

# Nyni dopocteme parcialni korelacni koeficent
r_yz.x <- cor(y.res,z.res)

# Dale podle vzorecku dopocteme prislusnou statistiku
n <- nrow(mfat1)
k <- 4
t.obs <- r_yz.x*sqrt(n-k-2)/sqrt(1-r_yz.x^2)
## Hodnota prislusne statistiky pro testovani nulovosti tohoto koeficientu

qt(0.975,n-k-2) 
##Kriticky obor je (-Inf,  -2.015368) u (2.015368, Inf) 

# A prislusna p-hodnota
2*(1-pt(abs(t.obs), n-k-2)) 
# Jak dopadl test?
# ...

# Analogicky muzeme postupovat pro dalsi promenne. Nyni budeme pocitat parcialni korelacni koeficient 
# mezi telesnou vahou (body.W) a tloustkou kozni rasy ve vysi 10. žebra (rib.F).
# Postup je stejny, uz bez komentare:
m3 <- lm(body.W ~ body.H + abdo.F + hip.F + quad.H, data = mfat1)
y.res <- m3$residuals

m4 <- lm(rib.F ~ body.H + abdo.F + hip.F + quad.H, data = mfat1)
z.res <- m4$residuals

r_yz.x <- cor(y.res,z.res)
r_yz.x

n <- nrow(mfat1)
k <- 4
t.obs <- r_yz.x*sqrt(n-k-2)/sqrt(1-r_yz.x^2)
t.obs

qt(0.975,n-k-2) ## 2.015368
2*(1-pt(abs(t.obs), n-k-2)) # p-hodnota
# Interpretujete vysledky testu
# ...

# Dalsi tri par. kor. koefiecienty bychom spocetli analogicky.

# Nyni provedeme vypocet pro vsechny naraz.

# Nacteme (pripadne i nainstalujeme) knihovnu ppcor
library(...)
# Opet doplnime pearson do pouzite metody
pcor(mfat1[,-3], method = "...")
## Najdete si vysledne p-hodnoty (muzete zkontrolovat spravnost s manualnim vypoctem)
## a interpretujte vysledky testu.
# ...
# Interpretujte take jednotlive par. kor. koeficienty.
# ...

# Nakonec spocteme koeficient vicenasobne korelace telesne vahy 
# s ostatnimi promennymi (tj. body.H, rib.F, abdo.F, hip.F, quad.H)

# Nejprve spocteme (Pearsonovu) korelaci mezi body.W a zbylymi promennymi. 
# Pozor na zadani vektoru urcujici sloupce tabulky mfat1. Zajimaji nas pouze 
# 2, 4, 5, 6, 7.
cor.yx <- cor(mfat1$..., mfat1[,c(...)], method = "...")
cor.yx

# Nyni spocteme korelacni matici zbylych promennych (tzn, zadavame to stejne
# jako na radku 334, jen bez prvniho argumentu)
cor.xx <- cor(mfat1[,c(...)], method = "...")
cor.xx

# Podle vzorecku dopocteme koeficient vicemasobne korelace. Funkce solve()
# pocita v R inverzni matici.
r.yx <- sqrt(cor.yx %*% solve(cor.xx) %*% t(cor.yx))
r.yx
## Interpretujte koeficient
## ...

# Jiny zpusob vypoctu:

# Sestavime linearni model zavislosti body.W na zbylych promennych
model <- lm(body.W ~ ... + ... + ... + ... + ..., data = ...)

# Ziskame odhad hodnot telesne vahy (body.W) z naseho modelu (najdeme pod fitted.values)
y.odhad <- model$...

# Vypocitame korelacni koeficient mezi telesnou vahou a nasim odhadem
r.yx <- cor(mfat1$..., y.odhad, method = "...")
## Dostavame stejny vysledek

# Vsimneme si jeste souvislosti mezi koeficientem vicenasobne korelace 
# a indexem determinace naseho modelu. Index determinace najdeme
# v summary(model) pod nazvem r.squared
ID <- summary(model)$...
sqrt(ID)

# Nyni podle vzorecku vypocteme hodnotu statistiky pro testovani nulovosti
# koeficientu vicenasobne korelace.
n <- nrow(mfat1)
k <- 5
F.obs <- (n-k-1) / k * r.yx^2 / (1-r.yx^2)

qf(0.95, df1 = k, df2 = n-k-1) 
# Kriticky obor je W = (2.42704, Inf). Statistika F.obs patri do W.
# Interpretujte vysledky testu
# ...

# Vypocet p-hodnoty. Pouzije funkci pro hodnotu distribucni funkce (p)
# F-rozdeleni (f) => tedy pf.
(1-...(F.obs, df1 = k, df2 = n-k-1))

# Hodnotu statistiky i p-hodnotu najdeme take ve vypisu summary
summary(...)