# 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, i jedan sa polinomom desetog
# stepena.

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 = 10, raw = TRUE), data = train_set) #za ekstremniji primer staviti degree = 10
summary(model2) #R^2 = 0.7366
yn = predict(model2, newdata = data.frame(weight = xn))
lines(xn, yn, col = "blue", lwd = 3)

# R^2 modela desetog 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 nepotrebno ici preko drugog stepena.
# 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 specificni 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 predvidjanja 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, sto bi nas navelo na pogresan zakljucak.

## 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 svaka opservacija nekad biti u tacno jednom od test skupova.
library(dplyr) #biblioteka je potrebna za funkciju bind_rows

# 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 = 10 #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(Auto[fold_idx,])
})

# Nakon sto smo izdvojili foldove, odvajamo trening i test k puta i belezimo
# predvidjanja u svakoj iteraciji. Cilj je da za kvadratni model potrosnje od
# mase automobila izvrsimo sto precizniju evaluaciju.
cross_validated_preds = lapply(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 vrsimo predvidjanje na test skupu 
  model <- lm(mpg ~ weight + I(weight^2), data = train)
  predicted_values = predict(model, newdata = test)
  return(predicted_values)
})
# Spajamo listu u jedan vektor:
cross_validated_preds = unlist(cross_validated_preds)
y_pred = cross_validated_preds

# Sada mozemo da racunamo mere kvaliteta nad svim podacima.
# Izracucemo, primera radi, R^2 i MAPE
y = Auto[shuffled_idx,]$mpg
e = y - y_pred
SSE = sum(e^2)
SSTO = sum((y - mean(y))^2)
R2 = 1 - SSE/SSTO
R2
# Model objasnjava 71% disperzije u ciljnoj promenljivoj

MAPE = mean(abs((y - y_pred)/y))
MAPE
# Model u proseku masi 13%

# Nakon sto izvrsimo unakrsnu validaciju, izvrsimo trening modela nad svim
# podacima da dobijemo finalni model
final_model = lm(mpg ~ weight + I(weight^2), data = Auto)
summary(final_model)

# Napomena: Korak pravljenja foldova smo mogli da preskocimo primenom funkcije
# createFolds() iz paketa caret, sto cemo u buducnosti i da radimo.