library(mlbench)
library(caret)
library(glmnet)
library(corrplot)

## Regularizacija u linearnim modelima
# Radimo sa Sonar podacima sa proslog casa, imali smo problem sa
# preprilagodjavanjem modela logisticke regresije
data(Sonar)

# Odvajamo odmah deo podataka za trening i validaciju, posto cemo kasnije da 
# koristimo unakrsnu validaciju
set.seed(216)
train_val_idx = createDataPartition(Sonar$Class, p = 0.8, list = FALSE)
train_val = Sonar[train_val_idx,]
test = Sonar[-train_val_idx,]

# Obicna logisticka regresija daje warning
log_model = glm(Class ~., data = train_val, family = 'binomial')
summary(log_model)
# Mozemo da vidimo da je p-vrednost za svaki atribut prakticno 1. Takodje
# vidimo da odredjeni atributi imaju jako visoke koeficijente (red velicine 10^3)

# Imamo previse atributa na samo 208 instanci, osim toga postoji jako korelacija
# izmedju njih, tad je cesto preprilagodjavanje
corrplot(cor(Sonar[,-61]))

# Mogli bismo da izbacujemo rucno prediktore, da izbegnemo jake korelacije i
# da smanjimo broj prediktora, ali u ovom slucaju to nikako nije prakticno jer
# imamo mnogo prediktora

# Koristicemo funkciju glmnet iz istoimenog paketa, omogucava nam da
# regularizujemo bilo koji generalizovani linearni model, u koje izmedju ostalog
# spada logisticka regresija
?glmnet
# Vidimo da je potrebno da se podaci odvoje u x i y
X_train_val = train_val[, -61]
y_train_val = as.numeric(train_val[, 61])-1
X_test = test[, -61]
y_test = as.numeric(test[, 61])-1

# Hiperparametar za regularizaciju nam je lambda, sto je lambda vece, to je
# regularizacija jaca. Radi tako sto sprecava model da ima prevelike koeficijente
# uz atribute. Vrednost lambda = 0 odgovara modelu bez regularizacije.
# Standardno je da se pri regularizaciji skaliraju podaci da bi se jednakom
# merom 'kaznila' velicina koeficijenata za svaki atribut. Glmnet to 
# podrazumevano radi. Navodimo standardize = FALSE jer su nama svi podaci vec
# iste skale.
regularized_log_model = glmnet(X_train_val, y_train_val, lambda = 10, alpha = 0, family = 'binomial', standardize = FALSE)
regularized_log_model$beta #koeficijenti su sada jako mali

# Nije ocigledno kako odabrati vrednost hiperparametra lambda, a za razliku
# od KNN-a, ne postoji ni bas jasno intuitivno znacenje tog broja koje bi nas
# navelo na izbor.
# Funkcija glmnet moze da napravi modele za vise vrednosti lambda odjednom
regularized_log_models = glmnet(X_train_val, y_train_val, lambda = c(0.01, 1, 10, 100), alpha = 0, family = 'binomial', standardize = FALSE)
plot(regularized_log_models, xvar = 'lambda') #zavisnost vrednosti koeficijenata od logaritma lambde

# Da bismo odabrali optimalnu vrednost lambde, potrebno je odvojimo trening
# trening i validacioni skup i da na validacionom uporedimo modele. Unakrsna
# validacija je naravno jos bolja varijanta, ali radi jednostavnosti koda
# radimo samo obicnu podelu.

set.seed(216)
train_idx = createDataPartition(train_val$Class, p = 0.75, list = FALSE)
train = train_val[train_idx,]
val = train_val[-train_idx,]

X_train = train[,-61]
y_train = as.numeric(train[,61])-1
X_val = val[,-61]
y_val = as.numeric(val[,61])-1

# Za potencijalne vrednosti lambda cemo uzeti stepene vrednosti
lambdas = exp(seq(-10,3, length.out = 50))

regularized_log_models = glmnet(X_train, y_train, lambda = lambdas, alpha = 0, family = 'binomial', standardize = FALSE)
plot(regularized_log_models, xvar = 'lambda')

# Mozemo individualne modele da koristimo za predikciju
# Izracunacemo tacnost klasifikacije na validacionom skupu za svaku vrednost
# lambda i prikazati ih na grafiku.
accs = c()
for(lambda in lambdas)
{
  y_pred = predict(regularized_log_models, s = lambda, newx = as.matrix(X_val), type = 'class')
  accuracy = mean(y_pred == y_val)
  accs = c(accs, accuracy)
}
plot(accs ~ log(lambdas), type = 'l')
# Vidimo da nakon nule (lambda = 1) modeli postaju skoro pa neupotrebljivi,
# regularizacija je prejaka.
# Postoji mnogo naglih skokova i padova i nije jasno sta odabrati za lambda
# Bolji bi bilo da smo radili unakrsnu validaciju za konzistentnije ocene

# Srecom, glmnet vec dolazi sa podrskom za unakrsnu validaciju
# type.measure = 'class' govori da je parametar za minimozovanje greska
# misklasifikacije, odnosno 1 - tacnost
regularized_models_cv <- cv.glmnet(as.matrix(X_train_val), as.matrix(y_train_val), nfolds = 10, lambda=lambdas, type.measure = "class", family="binomial", standardize = FALSE)
plot(regularized_models_cv) #grafik performansi u odnosu na logariam lambda
# Dobijamo slicno ponasanje kao i sto smo rucno racunali, ali je jasnija
# pravilnost zbog unakrsne validacije
best_lambda = regularized_models_cv$lambda.min # lambda za koje je greska najmanja

# Dalje mozemo da vidimo performanse modela na test skupu
# Ali prvo ga natreniramo na trening + validacioni
best_regularized_model = glmnet(X_train_val, y_train_val, lambda = best_lambda, family = 'binomial', standardize = FALSE)
y_pred_reg = predict(best_regularized_model, s = best_lambda, newx = as.matrix(X_test), type = 'class')
test_acc_reg = mean(y_pred_reg == y_test)
test_acc_reg #0.7317073
y_pred_proba = predict(log_model, newdata = test, type = 'response') #ne postoji opcija za class
y_pred = as.numeric(y_pred_proba > 0.5)
test_acc = mean(y_pred == y_test)
test_acc #0.6097561


## Agregacija

# Podrazumeva pravljenje vise modela i uzimanje predvidjanja kao prosek njihovih
# predvidjanja. U praksi se koristi najcesce nad stablima odlucivanja koja se
# rade na SS4, ali radi za razne tipove modela

# Napravicemo 100 modela, gde svaki ima samo 10 atributa
# Cesto se i uzima samo deo podataka, ali necemo da komplikujemo
n_models = 100
attributes_per_model = 10

models <- list()
for(i in 1:n_models)
{
  predictor_idx = sample(1:60, size = 10)
  X_data = X_train_val[,predictor_idx]
  models[[i]] <- glm(y_train_val ~ ., data = X_data, family = 'binomial')
}

# Predikcija na testu
# Mozemo da zamislimo proces predikcije kao da glasanja svakog modela
# Brojacemo za svaku instancu koliko modela glasaju za pozitivnu klasu
n_votes = rep(0, nrow(X_test))
for(i in 1:n_models)
{
  y_pred_proba = predict(models[[i]], newdata = X_test, type = 'response')
  y_pred = as.numeric(y_pred_proba > 0.5)
  n_votes = n_votes + y_pred
}
y_pred_aggr = as.numeric(n_votes > n_models/2) #ako vise od pola modela glasaju
test_acc_aggr = mean(y_pred_aggr == y_test)
test_acc_aggr #0.7560976