# Ucitavamo podatke Auto iz biblioteke ISLR
library(ISLR)
data(Auto)
View(Auto)
N = nrow(Auto)

# Zelimo da modeluje potrosnju na osnovu tezine automobila

set.seed(216) #stavljamo seed za konzistentne rezultate
# Delimo nasumicno na trening i test skup, 70% podataka ide u trening
train_indexes = sample(1:N, floor(0.7*N))
train_set = Auto[train_indexes,]
test_set = Auto[-train_indexes, ]

plot(mpg ~ weight, data = train_set) #vidi se da zavisnost nikako nije linearna

# Pravimo 2 modela, jedan sa kvadratnim clanom, a jedan sa kvadratnim i kubnim clanom, i cilj je da ih uporedimo

model1 <- lm(mpg ~ weight + I(weight^2), data = train_set)
summary(model1) #R^2 = 0.7292
xn = seq(1500, 5500, length.out = 200)
yn = predict(model1, newdata = data.frame(weight = xn))
lines(xn, yn, col = "red", lwd = 3)

model2 <- lm(mpg ~ poly(weight, degree = 3, raw = TRUE), data = train_set) #za ekstremniji primer staviti degree = 10
summary(model2) #R^2 = 0.7295
yn = predict(model2, newdata = data.frame(weight = xn))
lines(xn, yn, col = "blue", lwd = 3) #modeli su jako slicni
# R^2 modela treceg stepena mora da bude veci, jer ima sve atribute koje ima model drugog stepena
# p-vrednosti t-testova individualnih koeficijenata modela 2 su visoki
# Zbog toga i cinjenice da je R^2 drugog modela jedva veci, mozemo razumno da posumnjamo da je kubni clan nepotreban
# Ovo mozemo da potvrdimo i testom kolicnika verodostojnosti:
anova(model1, model2) #visoka p-vrednost

# Test kolicnika verodostojnosti, kao i statisticki testovi za individualne koeficijente su specificne za linearne modele
# Kada radimo sa proizvoljnim modelom, ne mozemo da upotrebimo te testove
# Zbog toga je izmedju ostalog korisna podela na trening i test
# Na test skupu mozemo 'fer' da uporedimo dva modela, posto rade sa podacima koje nisu videli

# Pravimo predikcije oba modela na trening i test skupu i racunamo srednjekvadratnu gresku
train_predictions_1 = predict(model1, newdata = train_set)
train_predictions_2 = predict(model2, newdata = train_set)
test_predictions_1 = predict(model1, newdata = test_set)
test_predictions_2 = predict(model2, newdata = test_set)

train_MSE_1 = mean((train_set$mpg - train_predictions_1)^2)
train_MSE_2 = mean((train_set$mpg - train_predictions_2)^2)
test_MSE_1 = mean((test_set$mpg - test_predictions_1)^2)
test_MSE_2 = mean((test_set$mpg - test_predictions_2)^2)

train_MSE_1
train_MSE_2
test_MSE_1
test_MSE_2
# Prvi model je losiji na trening skupu, sto smo i znali
# Ipak, bolji je na test skupu - iz cega mozemo da dodjemo do zakljucka da je zaista bolji

# Razlike u performansama dva modela su male
# Moglo je da se desi da su se nasumicno odabrani trening i test skup bas tako poklopili, da prvi model bude bolji
# Da bi smo tako nesto izbegli, mozemo da reuzorkujemo trening i test, i da napravimo m podela na trening i test
# da bismo videli koji model je zaista bolji u proseku.
# U kodu ispod 10 puta ponavljamo postupak iz proslog primera, belezimo razliku gresaka:
m = 10
diff = replicate(n = m, {
  train_indexes = sample(1:N, floor(0.7*N))
  train_set = Auto[train_indexes,]
  test_set = Auto[-train_indexes, ]
  
  model1 <- lm(mpg ~ weight + I(weight^2), data = train_set)
  model2 <- lm(mpg ~ poly(weight, degree = 3, raw = TRUE), data = train_set) 
  
  test_predictions_1 = predict(model1, newdata = test_set)
  test_predictions_2 = predict(model2, newdata = test_set)
  
  test_MSE_1 = mean((test_set$mpg - test_predictions_1)^2)
  test_MSE_2 = mean((test_set$mpg - test_predictions_2)^2)
  
  test_MSE_2 - test_MSE_1
})
barplot(diff)
# Vidimo sad da je prvi model konzistentno bolji
# Ali vidimo i da je u jednom slucaju drugi model bio bolji
# Da smo samo jednom obavili podelu na trening i test, moguce je da bismo bas takvu situaciju dobili

## Unakrsna validacija (cross validation)
# Ako vec unapred znamo da zelimo da napravimo podelu na trening i test vise puta, bolje je da to ne radimo skroz nasumicno
# Unakrsna validacija je metod koji nam omogucava da te podele uradimo dosta vise sistematski
# Osigurava da ce svaki podataka biti nekad biti u tacno jednom od test skupova

# Za naredni primer ucitavamo poznatu bazu za dijabetes
library(mlbench)
library(dplyr) #biblioteka je potrebna za funkciju bind_rows

data(PimaIndiansDiabetes)
data = PimaIndiansDiabetes
N = nrow(data)

# Algoritam unakrsne validacije deli podatke na k 'foldova' - disjunktnih skupova
# U svakoj iteraciji jedan od foldova sluzi za test skup za model treniran na uniji ostalih
k = N #broj foldova
shuffled_idx = sample(1:N) #nasumican niz od 1 do N
# Blok ispod predstavlja kod za pravljenje foldova
# Za svaku vrednost i od 1 do k izdvajamo dataframe, i te dataframe-ove pakujemo u listu
folds = lapply(1:k, function(i) {
  lower_index = floor(N*(i-1)/k) + 1 #donja granica za izmesane indekse i-tog folda
  upper_index = floor(N*(i/k)) #gornja
  fold_idx = shuffled_idx[lower_index:upper_index]
  return(data[fold_idx,])
})

# Nakon sto smo izdvojili foldove, odvajamo trening i test k puta i racunamo tacnost modela u svakoj iteraciji
cross_validated_accs = sapply(1:k, function(i) {
  # Za trening skup uzimamo uniju svih foldova osim i-tog, a za test i-ti fold
  train = bind_rows(folds[-i]) #funkcija bind_rows radi isto sto i rbind, ali radi za listu dataframe-ova
  test = folds[[i]] #obavezno duple uglaste zagrade, da bismo izdvojili element liste, jer inace izdvajamo podlistu sa jednim elementom
  # Pravimo model logisticke regresije na trening skupu i racunamo tacnost na test skupu 
  model <- glm(diabetes ~., data = train, family = 'binomial')
  predicted_proba = predict(model, newdata = test)
  y_pred = as.numeric(predicted_proba > 0.5)
  test_acc = mean(y_pred == as.numeric(test$diabetes)-1)
  return (test_acc)
})
cross_validated_accs #vidimo rezultate za svaku iteraciju unakrsne validacije, postoji znacajna varijacija
mean(cross_validated_accs) #prosek mozemo da uzmemo kao konacnu ocenu tacnosti


# Bibliotecka podrska za unakrsnu validaciju
# Kada delimo podatke na trening i na test, ili na foldove, zelimo da svaki skup predstavlja reprezentativan uzorak
# U slucaju klasifikacije, standardna praksa je da podelimo podatke tako da se odrzi odnos kategorija ciljne promenljive
library(caret)
# Funkcija createFolds vraca listu indeksa foldova
# Kao parametar predajemo kategoricku ciljnu promenljivu diabetes, i onda ce u svakom foldu biti isti odnos kategorija
fold_idx <- createFolds(data$diabetes, k = 5)
# Dalje, slicno racunamo tacnost u svakoj iteraciji algoritma
cross_validated_accs = sapply(1:k, function(i) {
  train = data[unlist(fold_idx[-i]),] #unlist listu vektora spaja u jedan
  test = data[fold_idx[[i]],]
  model <- glm(diabetes ~., data = train, family = 'binomial')
  predicted_proba = predict(model, newdata = test)
  y_pred = as.numeric(predicted_proba > 0.5)
  test_acc = mean(y_pred == as.numeric(test$diabetes)-1)
  return (test_acc)
})
cross_validated_accs
mean(cross_validated_accs)