# Grebena regresija (Ridge regression)

# Kako je kod multikolinearnosti matrica X^T * X blizu singularne, dobijamo
# nestabilne ocene koeficijenata. Jedan nacin da se resi taj problem je da se
# umesto X^T * X gleda matrica (X^T * X + lambda*I), koja, za neko lambda, biti
# regularna i imati stabilan inverz. Tako dobijamo ocenu koeficijenata ridge
# regresijom beta = (X^T * X + lambda*I)^{-1} * X^T * Y.

# U kontekstu masinskog ucenja, ovo se naziva L2 regularizacijom i trazenje
# koeficijenta je ekvivalentno minimizovanju SSE/n + L * ||beta||_2, odnosno
# minimizujemo srednjekvadratnu gresku (kao u obicnoj MNK), uz uvodjenje
# regularizacionog izraza koji predstavlja L2 norma koeficijenata beta. Ono sto
# se time dobija je da kaznjavamo velike vrednosti koeficijenata, cime
# povecavamo stabilnost ocena.

# Preporuceno je za ridge regresiju u R koristiti glmnet paket
# (moze i ridge.lm iz MASS)
library(glmnet)

# Prvo pravimo matricu prediktora iz baze
X <- as.matrix(longley[, -7])
# Model pravimo funkcijom glmnet, koja podrazumevano pravi LASSO model, a da bi
# napravila Ridge model, potrebno je proslediti parametar alpha = 0.
# Njoj se prosledjuju *matrica* prediktora i vektor zavisne promenljive
ridge <- glmnet(X, longley$Employed, alpha = 0)
# Ovaj funkcija pravi vise modela za razne vrednosti regularizacionog parametra
# lambda, koji odredjuje koliku tezinu dajemo regularizacionom izrazu nasuprot
# srednjekvadratnoj gresci.

# Od ovih raznih modela, treba odrediti onaj koji je najbolji. Jedan nacin da se
# to odradi je koriscenje unakrsne validacije. Implementiracemo rucno jedan od
# najjednostavnijih primera.

# Implementiracemo funkciju koja unakrsnom validacijom odredjuje ocene
# srednjekvadratne greske za svako pojedinacno lambda, pa vraca i model za koji
# je vrednost ocene greske najmanja. Najjednostavnija varijanta unakrsne
# validacije je ona sa izbacivanjem samo jednog elementa skupa podataka
# (leave-one-out cross validation). Detaljan opis se moze naci u skripti za
# Masinsko ucenje od profesora Mladena Nikolica (ml.matf.bg.ac.rs).
cv_longley <- function(lambdas) {
  cv_results <- list(lambda = lambdas, mse = numeric(length(lambdas)))
  X <- as.matrix(longley[, -7])
  
  for (k in seq_along(lambdas)) {
    
    lambda <- lambdas[k]
    preds <- numeric(nrow(longley))
    
    for (i in 1:nrow(longley)) {
      m <- glmnet(X[-i, ], longley$Employed[-i], lambda = lambda, alpha = 0)
      preds[i] <- predict(m, as.matrix(longley[i, -7]))
    }
    cv_results$mse[k] <- mean( (longley[ ,"Employed"] - preds)^2 )
  }
  best <- cv_results$lambda[which.min(cv_results$mse)]
  cv_results$best_model <- lm.ridge(Employed ~ ., longley[-i,], lambda = best)
  cv_results
}

results <- cv_longley(seq(0, 0.1, 0.001))
# Nije lose nacrtati grafik zavisnosti greske od lambda
plot(results$lambda, results$mse, type = "l")
# vidi se ocigledan minimum i lambda za koje se dostize je...
results$best_model$lambda # = 0.001

ridge_best2 <- results$best_model
coef(ridge_best2)

# Najbolji model, tj. najbolje lambda za Ridge model mozemo odrediti koriscenjem
# unakrsne validacije koja je implementirana u funkciji cv.glmnet. Glmnet
# implementira kompleksniju varijantu unakrsne validacije od leave-one-out.
# Takodje joj prosledjujemo matricu prediktora i zavisnu promenljivu
cv_model <- cv.glmnet(X, longley$Employed, alpha = 0, lambda = seq(0, 0.1, 0.001))
# ignorisite grouped=FALSE warning

# Crtanjem grafika unakrsne validacije, na kom se vidi srednjekvadratna greska
# za razlicite vrednosti lambda mozemo oceniti optimalnu vrednost parametra
# lambda. Dve vertikalne sive linije oznacavaju log(lambda) za koju je MSE
# najmanja, i najvece log(lambda) takvo da MSE nije veca od 1 standardne
# devijacije od minimalne MSE.
plot(cv_model)
# Ove vrednosti lambda mozemo da dobijemo, redom, pomocu
cv_model$lambda.min
cv_model$lambda.1se
# Preporuceno je da se koristi lambda.1se
# Koeficijenttakvog modela su
coef(cv_model)
# Vidimo da su izabrani isti prediktori koje smo i mi izabrali ispitivanjem
# multikolinearnosti prosli cas.


## LASSO

# Za razliku od ridge regresije, koja koristi L2 regularizaciju uz minimizaciju
# srednjekvadratne greske, LASSO regresija pretpostavlja L1 regularizaciju, cime
# se koeficijenti traze minimizovanjem SSE/n + L * ||beta||_1. Glavna razlika
# izmedju Ridge i LASSO je ta sto LASSO moze da se koristi za odabir prediktora,
# buduci da regularizacioni izraz, koji je jednak sumi apsolutnih vrednosti
# koeficijenata moze da svede koeficijente na tacno 0, dok se kod L2
# regularizacije desava da se koeficijenti mogu pribliziti nuli, mada je ne
# dostizu nikad.

# u R-u je za LASSO preporuceno koriscenje glmnet paketa
library(glmnet)

# Opet cemo gledati longley podatke i videti kako se LASSO moze koristiti i za
# selekciju prediktora.

# Prvo pravimo matricu standardizovanih prediktora iz baze
Xs <- as.matrix(scale(longley[, -7]))
# Model pravimo funkcijom glmnet, koja podrazumevano pravi LASSO model
# njoj se prosledjuju *matrica* prediktora i vektor zavisne promenljive
lasso <- glmnet(Xs, longley$Employed)
# Ovaj funkcija pravi vise modela za razne vrednosti regularizacionog parametra
# lambda, koji odredjuje koliku tezinu dajemo regularizacionom izrazu nasuprot
# srednjekvadratnoj gresci.
# Mozemo videti tok kretanja koeficijenata modela u zavisnosti od vrednosti
# parametra lambda
plot(lasso, xvar="lambda", label=TRUE)
# Svaka linija predstavlja vrednost jednog koeficijenta. Vidimo da za veliko
# lambda svi koeficijenti su nule.

# Najbolji model, tj. najbolje lambda za LASSO model mozemo odrediti koriscenjem
# unakrsne validacije, koja je implementirana u funkciji cv.glmnet
# Takodje joj prosledjujemo matricu prediktora i zavisnu promenljivu
cv_model <- cv.glmnet(Xs, longley$Employed) # ignorisite grouped=FALSE warning

# Crtanjem grafika unakrsne validacije, na kom se vidi srednjekvadratna greska
# za razlicite vrednosti lambda mozemo oceniti optimalnu vrednost parametra
# lambda. Dve vertikalne sive linije oznacavaju log(lambda) za koju je MSE
# najmanja, i najvece log(lambda) takvo da MSE nije veca od 1 standardne
# devijacije od minimalne MSE.
plot(cv_model)
# Ove vrednosti lambda mozemo da dobijemo, redom, pomocu
cv_model$lambda.min
cv_model$lambda.1se
# Preporuceno je da se koristi lambda.1se
# Koeficijenttakvog modela su
coef(cv_model)
# Vidimo da su izabrani isti prediktori koje smo i mi izabrali ispitivanjem
# multikolinearnosti prosli cas.