library(readxl)       # read_excel
library(numDeriv)     # grad
library(shape)        # Arrows

data <- read_excel("01_nhanes-weight.xls")
data$Age <- data$Age / 12
data$Gender[data$Gender == 0] <- "F"
data$Gender[data$Gender == 1] <- "M"

# chybajuce meranie
which(!complete.cases(data))
data[!complete.cases(data), ]

# vyhodenie chybajuceho
data <- na.omit(data)

####
# Priamkova regresia ####

m1 <- lm(Weight ~ Height, data = data)
m1

plot(data$Height, data$Weight)
abline(m1$coef, col="red", lwd=2)

####
# Gradientna metoda pre LinReg ####

# Model na standardizovanych datach
Height_scale <- scale(data$Height)
y <- data$Weight
X <- cbind(1, Height_scale)
m0 <- lm(y ~ X - 1) # -1, aby sme nepridali dalsi intercept

# Priprava

# sucet stvorcov chyb
J <- function(a, b) {
  return(sum((y - X %*% c(a, b))^2))
}

f <- function(beta) {
  return(J(beta[1], beta[2]))
  # verzia s vektorovym vstupom
  # kvoli funkcii na vypocet gradientu
}

# hodnoty ucelovej fcie
z <- matrix(NA, nrow=100, ncol=100)
aa <- seq(40, 120, len=100)
bb <- seq(0, 20, len=100) 
for(i in 1:100) {
  for(j in 1:100) {
    z[i, j] <- J(aa[i], bb[j])
  }
}

# Gradientna metoda

# pociatocne riesenie
a <- 40
b <- 0
J(a, b)

# 2 obrazky pripravit
par(mfrow = c(1, 2)) 

# Data
plot(Height_scale, data$Weight, ylim=c(0, max(data$Weight)))
abline(c(a, b), lwd=2, col="red") # dokreslenie pociatocnej priamky
abline(m0$coef, lwd=2, col="blue") # dokreslenie optimalnej priamky

# Vykreslenie ucelovej funkcie
contour(aa, bb, z)
points(m0$coef[1], m0$coef[2], pch=15, col="blue", cex=1.5)

# dokreslenie pociatocneho riesenia
points(a, b, pch=16, col="red")

# dokreslenie gradientu
df <- grad(f, c(a, b)) # numericky gradient (dal by sa aj vzorcekom)
lam <- 5e-5
lam*df
Arrows(a, b, a - lam*df[1], b - lam*df[2], col="red", lwd=2)   # -lambda*gradient

# nova iteracia
a <- a - lam*df[1]
b <- b - lam*df[2]

# Opakovane spustat od Data

###
# Model na Nestandardizovanych datach

y <- data$Weight
X <- cbind(1, data$Height)   # POTOM UKAZAT ako grad. met. vyzera bez standardizacie
m0 <- lm(y ~ X - 1) # -1, aby sme nepridali dalsi intercept

J <- function(a, b) {
  return(sum((y - X %*% c(a, b))^2))
}
f <- function(beta) {
  return(J(beta[1], beta[2]))
}
z <- matrix(NA, nrow=100, ncol=100)
aa <- seq(-100, -50, len=100)
bb <- seq(0, 2, len=100) 
for(i in 1:100) {
  for(j in 1:100) {
    z[i, j] <- J(aa[i], bb[j])
  }
}

# Gradientna metoda

a <- -80
b <- 0.5
J(a, b)
par(mfrow = c(1, 2)) 

# Data
plot(data$Height, data$Weight, ylim=c(0, max(data$Weight)))
abline(c(a, b), lwd=2, col="red") # dokreslenie pociatocnej priamky
abline(m0$coef, lwd=2, col="blue") # dokreslenie optimalnej priamky

# Vykreslenie ucelovej funkcie
contour(aa, bb, z, nlevels=20)
points(m0$coef[1], m0$coef[2], pch=15, col="blue", cex=1.5)

# dokreslenie pociatocneho riesenia
points(a, b, pch=16, col="red")

# dokreslenie gradientu
df <- grad(f, c(a, b)) # numericky gradient (dal by sa aj vzorcekom)
lam <- 1e-9
lam*df
Arrows(a, b, a - lam*df[1], b - lam*df[2], col="red", lwd=2)   # -lambda*gradient

# nova iteracia
a <- a - lam*df[1]
b <- b - lam*df[2]

# Opakovane spustat od Data

dev.off() # vynulovat vykreslovanie

####
# Weight vs. vsetky (to, co v Pythone) ####

m2 <- lm(Weight ~ ., data = data[, -1])
# y ~ . ... zoberie vsetky stlpce
# data[, -1] ... vyhodili sme stlpec Age

m2

summary(m2) # ovela viac info

# Predikcie
data_new <- data.frame(Age = sample(data$Age, 5),
                       Height = sample(data$Height, 5),
                       Leg_Length = sample(data$Leg_Length, 5),
                       Arm_Length = sample(data$Arm_Length, 5),
                       Arm_circum = sample(data$Arm_circum, 5),
                       Waist = sample(data$Waist, 5))
data_new

predict(m2, newdata=data_new)

# Standardizacia
X_scaled <- data[, -c(1, 3)] # data bez Gender, Weight
X_scaled <- scale(X_scaled)
data_scaled <- as.data.frame(cbind(X_scaled, Weight=data$Weight))

m3 <- lm(Weight ~ ., data = data_scaled)

# zle predikcie
predict(m3, newdata=data_new)

# rovnako transformovat aj nove data
data_new_scaled <- t(t(data_new) - attr(X_scaled, "scaled:center"))
data_new_scaled <- t(t(data_new_scaled) / attr(X_scaled, "scaled:scale"))
# divne transponovanie je na to, aby sa priemery 
# odpocitali po stlpcoch a nie po riadkoch (podobne s.d.)

# nove data uz nebudu "standardizovane"
apply(data_new_scaled, 2, mean) # priemery
apply(data_new_scaled, 2, sd)   # std. odchylky

# spravne predikcie
predict(m3, newdata=as.data.frame(data_new_scaled))

# porovnanie s povodnymi
predict(m2, newdata=data_new) # OK

####
# Transformacie ####

# Weight vs. Waist, Waist^2

plot(data$Waist, data$Weight)

m4 <- lm(Weight ~ Waist + I(Waist^2), data = data)
m4

?I  # treat "as is" ... aby Rko nevymyslalo

# dokreslenie regresnej krivky
plot(data$Waist, data$Weight)
xx <- seq(from=min(data$Waist), to=max(data$Waist), len=100)
yy <- predict(m4, newdata=data.frame(Waist=xx))
lines(xx, yy, col="red", lwd=2)

# Dalsie mocniny

m5 <- lm(Weight ~ poly(Waist, degree=4), data = data)
m5
# poly ... polynom daneho stupna
#      ... ortogonalny polynom: vyrazne numericky stabilnejsie
# da sa aj manualne: Waist + I(Waist^2) + I(Waist^3) + ...

plot(data$Waist, data$Weight)
xx <- seq(from=min(data$Waist), to=max(data$Waist), len=100)
yy <- predict(m5, newdata=data.frame(Waist=xx))
lines(xx, yy, col="red", lwd=2)

# Mozeme aj rozne zvlastne kombinacie
m5 <- lm(Weight ~ log(Height) + Waist + I(Waist^2), data = data)
m5

# Interakcie
m6 <- lm(Weight ~ Age + Height + I(Age*Height), data = data)
m6

####
# Kategoricke premenne ####

m7 <- lm(Weight ~ Height + Gender, data = data)
m7

# vsetky:
m8 <- lm(Weight ~ ., data = data)
m8

# interakcie
m9 <- lm(Weight ~ Height + Gender + Height*Gender, data = data)
m9
# prida vsetky interakcie pre danu dvojicu (pre K kategorii: K-1 interakcii)

m9b <- lm(Weight ~ Height*Gender, data = data)
m9b
# Rko automaticky prida aj povodne premenne (ked nie je I(Height*Gender))

m9c <- lm(Weight ~ I(Height*Gender), data = data)
# I() nefunguje, lebo Gender nie su cisla a teda nevieme robit suciny

m9d <- lm(Weight ~ Height:Gender, data = data)
m9d
# iba interakcia: dvojbodka
# (divny model: chyba Height => pridalo vsetky interakcie)