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

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

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

# Podivejte se na zakladni udaje o datech
summary(...)

# Pokud chceme v R pouzit logistickou regresi, musi byt vysvetlovana promenna
# typu faktor. To zjistime pomoci funkce is.factor(), kterou aplikujeme na 
# promennou sex
is.factor(...)

# 1

# Abychom vypocetli zakladni charakteristiky, vyuzijeme fukci charakteristiky(), 
# kterou mate uvedou i v resenem prikladu. Oznacte celou funkci mysi a stisknete
# Ctrl + Enter. Vpravo v Global Enviroment by se vam mela objevit pod Functions.
charakteristiky <- function(x){
  # funkce pocitajici pocet pozorovani, prumer a smerodatnou odchylku
  # argument: x ... vektor
  # vraci: vektor (pocet pozorovani, prumer, smerodatna odchylka)
  x <- na.omit(x) #odstrani chybejici hodnoty
  n <- length(x)
  m <- mean(x)
  s <- sd(x)
  return(c(n, m, s))
}

# Budeme tedy aplikovat tuto funkci postupne na nase data vzdy zvlast pro zeny a
# zvlast pro muze. To udelame pomoci odkazovani se na cast vektoru skrze hranate zavorky.

# Zacneme s promennou mc2.L:
charakteristiky(data$...[data$sex == "f"])
charakteristiky(data$...[data$sex == "m"])

# Pro promennou mc2.CW zkuste spocitat tyto charakteristiky sami:
charakteristiky(...)
charakteristiky(...)

# Obdobne pro promenne mc2.MW a mc2.BW:
charakteristiky(...)
charakteristiky(...)

charakteristiky(...)
charakteristiky(...)

# 2

# Nyni budeme vykreslovat krabicove diagramy pomoci funkce boxplot(). 

# Nejprve nastavime grafy na okno 2x2
par(mfrow = c(2,2))

# Vykreslujeme boxplot pro kazdou vysvetluji promennou v zavislosti na pohlavi, 
# tedy pro mc2.L:
boxplot(data$mc2.L ~ data$sex, varwidth = T, notch = T, xlab = "Delka k.z. [mm]",
        ylab = "Pohlavi", horizontal = T)
## Pouzijte help(boxplot) a zjistete, co znamenaji argumenty
# varwidth, notch a horizontal.

# Pro promenne mc2.CW, mc2.MW a mc2.BW zkuste sestrojit formuli sami:
boxplot(data$... ~ data$..., varwidth = T, notch = T, xlab = "Sirka hlavicky k.z. [mm]",
        ylab = "Pohlavi", horizontal = T)
boxplot(data$... ~ data$..., varwidth = T, notch = T, xlab = "Sirka stredu tela k.z. [mm]",
        ylab = "Pohlavi", horizontal = T)
boxplot(data$... ~ data$..., varwidth = T, notch = T, xlab = "Sirka baze k.z. [mm]",
        ylab = "Pohlavi", horizontal = T)

# V boxplotech vsak vidime nekolik odlehlych pozorovani, a tak je nyni odstranime.
# Pouzijeme k tomu funkce:
# which() - specifikuje pro ktere prvky zkoumaneho vektoru plati zadana podminka, vyhodi poradi techto prvku
# which.max() - specifikuje maximalni prvek vektoru
# which.min() - specifikuje minimalni prvek vektoru

# Budeme postupovat od posledniho boxplotu (vpravo dole). Pro kontrolu je spravne poradi uvedeno za #.
# Najdeme poradi maximalniho prvku sirky baze (mc2.BW):
which.max(data$..) # 43

# Najdeme poradi maximalniho prvku sirky hlavicky (mc2.CW):
... 
# 43

# Najdeme poradi minimalniho prvku delky (mc2.L):
...
# 13

# Nyni zjistime, ktera jsou tri odlehla pozorovani u muzske populace v boxplotu 
# o delce zaprstni kosti. Jsou to podle grafu ta, ktera maji hodnotu delky vetsi nez 72.7:
which(data$... > ...) 
# 33 60 85

# A odhalime take posledni odlehle pozorovani. Jde o maximum delky k.z. pro zeny:
which.max(data$...[data$sex == "..."]) 
# 19

# Odlehla pozorovani odstranime z naseho datoveho souboru (z minulych cviceni vite, jak se provede),
# tj. odstranime prislusne radky v nasi tabulce a vytvorime tabulku novou, se kterou budeme nadale pracovat.
data1 <- data[-c(13, 19, 33, 43, 60, 85),]

# Pro kontrolu znovu vykreslime krabicove diagramy. Nyni by uz nikde nemela byt odlehla pozorovani.
par(mfrow = c(2,2))
boxplot(...)
boxplot(...)
boxplot(...)
boxplot(...)

# 3

# Pomoci dvouvyberovych t-testu budeme chtit zkoumat vztah mezi vysvetlovanou promennou
# a jednotlivymi vysvetlujicimi promennymi. Predpoklady na pouziti t-testu jsou vsak
# normalita a homogenita rozptylu obou testovanych skupin.

# Zacneme s prvni vysvetlujici (nezavislou) promennou "mc2.L" - delkou zaprstni kosti

# Pomoci Shapiro-Wilkova testu a Q-Q plotu nejprve overime normalitu
par(mfrow = c(1,2))
shapiro.test(data1$...[data1$sex == "f"])
qqnorm(data1$...[data1$sex == "f"], main = "Delka k.z. [mm] - zeny")
qqline(data1$...[data1$sex == "f"])

# A stejne pro muze
shapiro.test(data1$...[data1$sex == "..."])
qqnorm(data1$...[data1$sex == "..."], main = "Delka k.z. [mm] - muzi")
qqline(data1$...[data1$sex == "..."])

# Pomoci funkce var.test() overime shodnost rozptylu
var.test(data1$mc2.L ~ data1$sex)

# Oba predpoklady jsou splneny.

# A konecne po splneni vsech predpokladu pristoupime k t-testu
t.test(data1$mc2.L ~ data1$sex)
## Zamitame hypotezu o tom, ze prumerna delka k.z. pro muze a zeny je stejna. 

# Sami zkuste postup pro dalsi promenne:
par(mfrow = c(1,2))
shapiro.test(...)
qqnorm(..., main = "Sirka hlavicky k.z. [mm] - zeny")
qqline(...)

shapiro.test(...)
qqnorm(..., main = "Sirka hlavicky k.z. [mm] - muzi")
qqline(...)

var.test(...)
t.test(...)
# Okomentujte splneni predpokladu a vysledek t-testu:
# ...

par(mfrow = c(1,2))

shapiro.test(...)
qqnorm(..., main = "Sirka stredu tela k.z. [mm] - zeny")
qqline(...)

shapiro.test(...)
qqnorm(..., main = "Sirka stredu tela k.z. [mm] - muzi")
qqline(...)

var.test(...)
t.test(...)
# Okomentujte splneni predpokladu a vysledek t-testu:
# ...

par(mfrow = c(1,2))

shapiro.test(...)
qqnorm(..., main = "Sirka baze k.z. [mm] - zeny")
qqline(...)

shapiro.test(...)
qqnorm(..., main = "Sirka baze k.z. [mm] - muzi")
qqline(...)

var.test(...)
t.test(...)
# Okomentujte splneni predpokladu a vysledek t-testu:
# ...

# 4

# Dvouvyberove t-testy ukazuji rozdily u vsech promennych. Sestavime proto plny model.
# Logisticka regrese se v R tvori pomoci funkce glm(). Formule modelu se zadava podobne, jak jsme
# zvykli z drivejska. Dalsimi argumenty jsou "family" - typ logitove funkce (my budeme pouzivat vzdy 
# binomial(logit) a "data"). Zadejte tedy vsechny vysvetlujici promenne vpravo, vysvetlovana promenna
# sex bude ve formuli vlevo.
m.data1 <- glm(... ~ ... + ... + ... + ..., family = binomial(logit), data = data1)

# A ihned si muzeme vypsat informace o modelu a zjistit tak dilcimi testy vyznamosti, 
# ktere promenne jsou vyznamne.
summary(m.data1)

# 5

# Pro test vyznamnosti modelu jako celku, musime nejprve sestrojit model konstanty a nasledne
# tento model porovnat s nasim modelem.

# Model konstanty - do vysvetlujicich promennych zadejte 1
m0 <- glm(sex ~ ..., family = binomial(logit), data = data1)

# Test vyznamnosti modelu jako celku (u log. regrese musime zadat test = "Chisq")
anova(m0, m.data1, test = "Chisq")
# Hodnota testove statistiky (rozdil devianci) je 71,277, p-hodnota je 1.22e-14. 
# Zamitame tedy hypotezu o tom, ze by byl dostacuji model konstanty.

# 6

# Nyni budeme postupne vynechavat nevyznamne promenne. Muzeme vynechat bud 
# promennou mc2.CW, nebo mc2.BW. Vyzkousime postupne obe varianty.
# Nejdriv vynechame mc2.BW
m.data2 <- glm(sex ~ ... + ... + ..., family = binomial(logit), data = data1)
summary(m.data2)

# Porovname modely m.data1 a m.data2 mezi sebou
anova(..., ..., test = "Chisq")
# Ukazuje se, ze modely nejsou rozdilne, tudiz nam staci jednodussi model.

# Sestrojime model pouze s vysvetlujicimi promennymi mc2.L a mc2.MW
m.data3 <- glm(sex ~ ... + ...,family = binomial(logit), data = data1)
summary(m.data3)

# A porovname modely 2 a 3.
anova(..., ..., test = "Chisq")
# Ovsem tento model uz je nedostacujici.

# Vratime se na zacatek k modelu 1 a zkusime nyni vynechat promennou mc2.CW:
# Vytvorime tedy model m.data4 s regresory mc2.L, mc2.MW a mc2.BW
m.data4 <- glm(sex ~ ...,family = binomial(logit), data = data1)
summary(m.data4)

# A porovname jej s modelem 1
anova(..., m.data4, test = "Chisq")
# A ukazuje se, ze bychom mohli zvolit take tento model

# 7

# Stejne jako v pripade mnohorozmerne linearni regrese muzeme vyuzit stepwise procedury, 
# ktera nam pomuze s vyberem modelu. Syntaxe je uplne stejna (viz cviceni o mnohorozmerne 
# linearni regresi), akorat nyni pouzivame funkci glm() namisto lm().

# V metode backward zadavame do formule plny model
step(glm(sex ~ ... + ... + ... + ..., family = binomial(logit), data = data1),
     direction = "backward")

# V metode forward zadavame do formule model konstanty, do argumentu scope plny model
step(glm(sex ~ ..., family = binomial(logit), data = data1),
     scope = ~ ... + ... + ... + ... , direction = "forward")

# Metodu both muzeme zadat dvema zpusoby (od plneho modelu ke konstante, nebo naopak)
# Od plneho modelu ke konstante:
step(glm(sex ~ ... + ... + ... + ..., family = binomial(logit), data = data1),
     direction = "both")

# Od konstanty k plnemu modelu:
step(glm(sex ~ ..., family = binomial(logit), data = data1),
     scope = ~ ... + ... + ... + ..., direction = "both")

# Vsechny druhy stepwise procedury vybraly maximalni model (m.data1)

# 8

# Vypocteme Nagelkerkeho koeficient determinace pro vsechny modely.
# Nactete nejdrive knihovnu rsq (pripadne ji nainstalujte)
library(...)

# Nagelkerkeho koeficient determinace se vypocte nasledovne: Do funkce rsq
# se zada pozadovany model, argument type nastavime na "n". Vypoctete pro vsechny tri modely.
rsq(..., type = "...") 
## m1: nagelkerke 0.7193647

rsq(..., type = "...") 
## m2: nagelkerke 0.7042164

rsq(..., type = "...") 
## m4: nagelkerke 0.7025608

# Podle Nagelkerkeho koeficientu determinace vychazi o malinko lepe plny model.

# 9

# Sestavime klasifikacni tabulku pro vsechny nase modely
# m.data1:

# Nejprve pomoci funkce predict() predikujeme odezvu. Aplikujeme ji na nas model, ovsem 
# u log. regrese musime zadat type = "response".
# Zacneme s m.data1:
fitted1 <- predict(..., newdata = data1, type = "...")

# Tam kde je hodnota predikce nizsi nez 0,5 (predikujeme pravdepodobnost), bereme
# nase pozorovani jako "neuspech" (zaradime jako "f"), v opacnem pripade (vyssi nez 0,5)
# zaradime jako "m". To provedeme pomoci funkce ifelse. Do ni vlozime podminku a napiseme
# vysledek v pripade jeji splneni a nesplneni.
fitted1.cat <- ifelse(fitted1 < 0.5, "f", "m")

# Z predikovanych hodnot vytvorime kontigencni tabulku - funkce table(). Jako druhy argument 
# musime zadat skutecne hodnoty pohlavi z nasich dat
tab1 <- table(fitted1.cat, data1$sex)

# Prejmenujeme jeste radky (predikovane hodnoty) a sloupce (skutecne hodnoty) 
colnames(tab1) <- c("skut. f", "skut. m")
row.names(tab1) <- c("pred. f", "pred. m")

# Tabulku vypiseme
tab1

# A vypocteme relativni cetnosti spravne urcenych jedincu
sum(diag(tab1))/sum(tab1) 
# 0.8478261

# Analogicky postupujeme pro m.data2:
fitted2 <- predict(..., newdata = data1, type = "...")
fitted2.cat <- ifelse(... < 0.5, "f", "m")
tab2 <- table(... ,data1$sex)
colnames(tab2) <- c("skut. f", "skut. m")
row.names(tab2) <- c("pred. f", "pred. m")
tab2

sum(diag(tab2))/sum(tab2) 
# 0.8695652

# Analogicky postupujeme pro m.data4:
fitted4 <- predict(..., newdata = data1, type = "...")
fitted4.cat <- ifelse(... < 0.5, "f", "m")
tab4 <- table(...,data1$sex)
colnames(tab4) <- c("skut. f", "skut. m")
row.names(tab4) <- c("pred. f", "pred. m")
tab4

sum(diag(tab4))/sum(tab4)
# 0.826087

# Nejlepe vychazi model m.data2.

# 10

# Nyni vykrelime ROC krivky pro vsechny tri modely. K tomu budeme potrebovat knihovnu ROCR.
library(...)

# Nastavime obrazek na 1x3
par(mfrow = c(1,3))

# Nejprve pro m.data1.

# Toto vykresleni je ponekud komplikovanjesi. Nejprve vyuzijeme funkci prediction() na 
# nami predikovane hodnoty a skutecnou odezvu, kterou bereme jako numericky vektor.
preds1 <- prediction(fitted1, as.numeric(data1$sex))

# Pomoci funkce performance() dostaneme z promenne preds1 pozadovane udaje.
roc1 <- performance(preds1, "tpr", "fpr")
auc1 <- performance(preds1, "auc")

# Hodnotu obsahu plochy pod krivkou zaokrouhlime na 4 des. mista
auc.value1 <- round(as.numeric(auc1@y.values),4)

# A vykreslime graf
plot(roc1, main = "ROC Curve - m1", lwd = 2, col = "blue")

# Funkce grid() prida teckovanou mrizku do grafu
grid()

# Pridame graf nejhorsiho modelu:
lines(c(0,1),c(0,1))

# A doplnime do grafu text s hodnotou AUC:
text(0.5, 0.1, paste("AUC = ",auc.value1), col = "darkred")

# Stejnym zpusobem postupujeme pro dalsi dva modely:

preds2 <- prediction(fitted2, as.numeric(data1$sex))
roc2 <- performance(preds2, "tpr", "fpr")
auc2 <- performance(preds2, "auc")
auc.value2 <- round(as.numeric(auc2@y.values),4)
plot(roc2, main = "ROC Curve - m2", lwd = 2, col = "blue")
grid()
lines(c(0,1),c(0,1))
text(0.5, 0.1, paste("AUC = ",auc.value2), col = "darkred")

preds4 <- prediction(fitted4, as.numeric(data1$sex))
roc4 <- performance(preds4, "tpr", "fpr")
auc4 <- performance(preds4, "auc")
auc.value4 <- round(as.numeric(auc4@y.values),4)
plot(roc4, main = "ROC Curve - m4", lwd = 2, col = "blue")
grid()
lines(c(0,1),c(0,1))
text(0.5, 0.1, paste("AUC = ",auc.value4), col = "darkred")

# AUC je vzdy v intervalu (0.9 - 1], tedy predikcni schopnost vsech tri modelu je 
# vyborna (viz tabulka z prednasky).

# 11

# Volbu modelu nemusi byt vzdy jasna zalezitost a nyni se to ukazuje. Kterykoliv model, 
# ktery byste zvolili, je v podstate spravny. Zalezi na vas a na potrebach dane situace.

# Osobne bych vybral model m.data2. Hodnoty Nagelkerkoho koeficientu jsou srovnatelne 
# pro vsechny modely, stejne tak plocha pod ROC krivkou (AUC). Je vsak jednodussi nez plny model
# m.data1 a vychazi u nej nejlepsi relativni cetnost urceni jedincu.

# Vypis informaci o modelu provedeme pomoci tabulky summary:
summary(m.data2)

# Odhady regresnich koeficentu modelu
coef(m.data2)

# A jejich 95% intervaly spolehlivosti
lower <- coef(m.data2) - qnorm(0.975) * summary(m.data2)$coefficients[,2]
upper <- coef(m.data2) + qnorm(0.975) * summary(m.data2)$coefficients[,2]
cbind(lower,upper)

# V logisticke regresi vsak tyto koeficienty nijak neinterpretujme. Muzeme vsak 
# interpretovat pomery sanci, ktere vypocteme nasledovne:

# Pripomenme jeste, ze sance je definovana jako pravdepodobnost uspechu (v nasem 
# pripade jde o "m") deleno pravdepodobnosti neuspechu (v nasem pripade jde o "f").

# Pomery sanci 
exp(coef(m.data2))

# Ovsem exp(beta.0) = Intercept, v nasem pripade 9.662301e-23, jako jediny neudava pomer sanci, 
# ale primo sanci na uspech, pri nulovych hodnotach ostatnich regresoru (promennych).
# Tzn. pokud bychom uvazali hypotetickou situaci jedince se vsemi nulovymi rozmery zaprstnich
# kosti, tak sance ze pujde (podle naseho modelu) o muze je 9.662301e-23.

# Interpretace zbylych clenu vektoru exp(coef(m.data2)) - tj. pomeru sanci:
# Znamena to tedy, ze pokud se zvetsi hodnota delka zaprstni kosti (promenna mc2.L) o 1 mm
# a hodnoty ostatnich regresoru (promennych) zustanou stejne, tak sance na to, ze dane 
# pozorovani zaradime jako muze vzroste 1.45-krat.

# Obdobne: Pokud se zvetsi hodnota sirka hlavicky zaprstni kosti (promenna mc2.CW) o 1 mm
# a hodnoty ostatnich regresoru (promennych) zustanou stejne, tak sance na to, ze dane 
# pozorovani zaradime jako muze vzroste 2.1-krat.

# Podobne muzeme interpretovat i obraceny pomer sanci:
exp(coef(m.data2))^(-1)
# Pokud se ZMENSI hodnota sirka hlavicky zaprstni kosti (promenna mc2.CW) o 1 mm
# a hodnoty ostatnich regresoru (promennych) zustanou stejne, tak sance na to, ze dane 
# pozorovani zaradime jako muze "vzroste" 0.476-krat. (Jde o jedno a to same jako v predchozim 
# pripade - vime, ze nasobeni cislem mensim nez 1 dane cislo zmensuje, tj. stejne se zmensuje
# i sance v tomto pripade).

# Intervaly spolehlivosti pro pomery sanci
exp(cbind(lower,upper))