head <- 'TO-DO' # nacteni dat head.txt, uprava struktury, odstraneni NA

## vytvoreni jednoduche funkce pro vypocet zakl. charakteristik

charakteristiky <- function(x) {
  # funkce pocitajici pocet pozorovani, prumer a smerodatnou odchylku
  # argument: x ... vektor
  # vraci: vektor (pocet pozorovani, prumer, smerodatna odchylka)
  n <- length(x)
  m <- mean(x)
  s <- sd(x)
  return(c(n, m, s))
}

tapply(head$body.H, head$sex, charakteristiky) # charakteristiky pro body.H delene podle pohlavi
'TO-DO' # charakteristiky pro ostatni promenne

## vykresleni boxplotu pro jednotlive promenne podle pohlavi
par(mfrow = c(2, 3)) # nastaveni vice grafu
boxplot(head$body.H ~ head$sex,
  varwidth = T, notch = T, xlab = "Vyska [mm]",
  ylab = "Pohlavi", horizontal = T
)

'TO-DO' # doplnte dalsi boxploty pro zbyvajici promenne

## chceme testovat, zda se zeny a muzi lisi v danych rozmerech, tj. testujeme stredni hodnotu rozmeru pro muze a zeny
## to provadime pomoci "klasickeho" t-testu o rozdilu strednich hodnot
## ten ovsem plati za jistych predpokladu - musime tedy overit normalitu a shodu rozptylu

## overeni predpokladu - normalita
## zvlast pro muze a pro zeny
nortest::lillie.test(head$body.H[head$sex == "f"])
nortest::ad.test(head$body.H[head$sex == "m"])

par(mfrow = c(1, 2)) # overeni graficky pomoci qqplotu
qqnorm(head$body.H[head$sex == "f"], main = "Telesna vyska - zeny")
qqline(head$body.H[head$sex == "f"])
qqnorm(head$body.H[head$sex == "m"], main = "Telesna vyska - muzi")
qqline(head$body.H[head$sex == "m"])

## overeni predpokladu - rovnost rozptylu
var.test(head$body.H ~ head$sex)

## t-test o rozdilu strednich hodnot pro muze a zeny
t.test(head$body.H ~ head$sex, var.equal = T)


'TO-DO' # zopakujte overeni predpokladu a t-testy pro vsechny dalsi promenne


## sestaveni modelu
## provedeme pomoci fce glm s nastavenim family = binomial(logit)
m.head <- glm(sex ~ body.H + head.L + head.W + bigo.W + bizyg.W,
              family = binomial(logit), data = head)
summary(m.head) # vysledky modelu

## test vyznamnosti modelu - porovnani s modelem konstanty
m0 <- glm(sex ~ 1, family = binomial(logit), data = head)
anova(m0, m.head, test = "Chisq")

## sestaveni noveho modelu s vynechanim nevyznamnych promennych
m.head2 <- 'TO-DO'
  
anova(m.head2, m.head, test = "Chisq") # test dostatecnosti mensiho modelu

## vypocet intervalu spolehlivosti
coef(m.head)
lower <- coef(m.head) - qnorm(0.975) * summary(m.head)$coefficients[, 2]
upper <- coef(m.head) + qnorm(0.975) * summary(m.head)$coefficients[, 2]
cbind(lower, upper)

## pomer sanci a intervaly spolehlivosti pro pomer sanci
exp(coef(m.head))
exp(cbind(lower, upper))

## vyber modelu pomoci STEPWISE procedury
step(glm(sex ~ body.H + head.L + head.W + bigo.W + bizyg.W,
         family = binomial(logit), data = head),
     direction = "both"
)

'TO-DO' # analogicky vytvorte dalsi stepwise modely

## koeficient determinace
## ruzne typy - potrebujeme knihovnu rsq
library(rsq)
rsq(m.head, type = "n") # nagelkerke
rsq(m.head, type = "kl") # mcfadden
rsq(m.head, type = "lr") # cox and snell


## sestaveni klasifikacni tabulky pro zarazovani pozorovani mezi muze/zeny
fitted <- predict(m.head, newdata = head, type = "response")
fitted.cat <- ifelse(fitted < 0.5, "f", "m")
tab <- table(fitted.cat, head$sex)
tab

sum(diag(tab)) / sum(tab) # relativni cetnost spravne zarazenych pozorovani

## hodnoceni modelu pomoci ROC krivky a hodnoty AUC 
## funkce z knihovny ROCR
library(ROCR)
preds <- prediction(fitted, as.numeric(head$sex))
roc <- performance(preds, "tpr", "fpr")
auc <- performance(preds, "auc")
auc.value <- round(as.numeric(auc@y.values), 4)
plot(roc, main = "ROC Curve", lwd = 2, col = "blue")
grid()
lines(c(0, 1), c(0, 1))
text(0.5, 0.1, paste("AUC = ", auc.value), col = "darkred")



'TO-DO' # samostatne vypracujte nereseny priklad