library(ISLR2)
?Hitters

dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters <- na.omit(Hitters)
dim(Hitters)

## best subset selection
library(leaps)
?regsubsets

regfit.full <- regsubsets(Salary ~ ., Hitters, nvmax=20)
summary(regfit.full)
# npr. CRBI se nalazi u prvih 6 modela, ali ne i u sedmom. kod FSS to nije slucaj.


# ?plot.regsubsets
# plot(regfit.full, scale="bic")

## forward stepwise selection
regfit.fwd <- regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "forward")
summary(regfit.fwd)

regfit.bwd <- regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "backward")
summary(regfit.bwd)

## za p=1,..,6 sve metode daju iste modele, za p=7 to nije slucaj

coef(regfit.full, 7)
coef(regfit.fwd, 7)
coef(regfit.bwd, 7)


## unakrsna validacija za bss (k=10)

predict.regsubsets <- function(object, newdata, id, ...) {
  form <- as.formula(object$call[[2]])
  mat <- model.matrix(form, newdata)
  coefi <- coef(object, id = id)
  xvars <- names(coefi)
  mat[, xvars] %*% coefi
}
## ne postoji predict funkcija za regsubsets(), buduci da su objekti koje koristimo klase regsubsets, potrebno
## je samo napisati predict

n <- nrow(Hitters)
k <- 10

set.seed(1)
folds <- sample(rep(1:k, length = n)) 
## replace=FALSE -> npr. folds[1]=10 znaci da prvi element skupa za ucenje stavljamo u 10-ti blok
cv.errors <- matrix(NA, k, 19)

for (j in 1:k) {
  best.fit <- regsubsets(Salary ~ .,
                         data = Hitters[folds != j, ],
                         nvmax = 19)
  ### kada izbacujemo j-ti blok, bitno je da odabir najboljeg modela za svaki broj kovarijata (sto je dio treninga)
  ### radimo na podacima bez j-tog bloka! Na taj nacin j-ti blok zaista igra ulogu nezavisnog testnog skupa.
  for (i in 1:19) {
    pred <- predict(best.fit, Hitters[folds == j, ], id = i)
    cv.errors[j, i] <-
      mean((Hitters$Salary[folds == j] - pred)^2) 
    ## sadrzi testnu gresku na j-tom bloku za najbolji model s i varijabli (na skupu bez j-tog bloka)
    ## kako bi dobili CV procjenu testne greske za najbolji model s i varijabli, uzimamo prosjek testnih
    ## gresaka na svih k blokova (vidi dolje)
  }
}

mean.cv.errors <- apply(cv.errors, 2, mean) 

plot(mean.cv.errors, type = "b")  

p0=which.min(mean.cv.errors)
p0 ## 10

reg.best <- regsubsets(Salary ~ ., data = Hitters,
                       nvmax = 19)
coef(reg.best, p0)

## slicno bi bilo za forward i backward ss 


#### ridge i lasso

x <- model.matrix(Salary ~ ., Hitters)[, -1] ## generira matricu podataka X
y <- Hitters$Salary

dim(x)


library(glmnet)
grid <- 10^seq(10, -2, length = 100)
ridge.mod <- glmnet(x, y, alpha = 0, lambda = grid)
## alpha=0 odgovara ridge regresiji, a alpha=1 je lasso (vidi help)
## ako ne zadamo vrijednosti za lambda, odredit ce se automatski (tako inace radimo)
## po deafultu se kovarijate standardiziraju (inace je kazna "nepostena")

dim(coef(ridge.mod))
## za svaki od 100 vrijednosti od lambda imamo 19+1=20 koeficijenata

plot(ridge.mod, xvar="lambda", label=TRUE) ## veci lambda znaci vece smanjenje
plot(ridge.mod, xvar="norm", label=TRUE) ## x-os je L1 norma koeficijenata (veca norma manji lambda)


# predict radi slicno kao i prije, npr. mozemo dobiti koeficijente za novu vrijednost za lambda, npr. 50,
predict(ridge.mod, s = 50, type = "coefficients")[,1]
?predict.glmnet

## testna greska
set.seed(1)
train <- sample(1:nrow(x), nrow(x) / 2)
test <- (-train)
y.test <- y[test]

# npr. za lambda=4
ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid, thresh = 1e-12) 
## zadnji parametar?
ridge.pred <- predict(ridge.mod, s = 4, newx = x[test, ])
mean((ridge.pred - y.test)^2)

# opcenito, koristimo unakrsnu validaciju
set.seed(1)
cv.out <- cv.glmnet(x[train, ], y[train], alpha = 0)
plot(cv.out) ## na gornjoj x-osi su stupnjevi slobode, a dvije vertikalne crte oznacavaju
## lambda koji minimizira CV gresku, odnosno lambda dobiven uz one SE rule (veci lambda manja kompleksnost!)

?cv.glmnet ## vidi Value
?plot.cv.glmnet

bestlam <- cv.out$lambda.min
bestlam
## ili
# bestlam <- cv.out$lambda.1se
# bestlam

ridge.pred <- predict(ridge.mod, s = bestlam, newx = x[test, ])
mean((ridge.pred - y.test)^2)  # [1] 139856.6

# odabrani model je dakle
out <- glmnet(x, y, alpha = 0)
predict(out, type = "coefficients", s = bestlam)[,1]


## lasso
lasso.mod <- glmnet(x[train, ], y[train], alpha = 1, lambda = grid)
plot(lasso.mod, xvar="lambda", label=TRUE) ## za razliku od ridge regresije, lasso radi odabir kovarijata

ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = grid, thresh = 1e-12) 
plot(ridge.mod, xvar="lambda", label=TRUE) ## lambde nisu usporedive, ali stupnjevi slobode jesu

plot(lasso.mod, xvar="norm", label=TRUE)
plot(ridge.mod, xvar="norm", label=TRUE)

set.seed(1)
cv.out <- cv.glmnet(x[train, ], y[train], alpha = 1)
plot(cv.out)

bestlam <- cv.out$lambda.min
lasso.pred <- predict(lasso.mod, s = bestlam, newx = x[test, ])
mean((lasso.pred - y.test)^2) # [1] 143673.6
## slicno kao kod ridge regresije, s tim da koristimo manje stupnjeva slobode (11 umjesto 19) te usput radimo
## i odabir kovarijata

out <- glmnet(x, y, alpha = 1, lambda = grid)
lasso.coef <- predict(out, type = "coefficients", s = bestlam)[,1]
lasso.coef
lasso.coef[lasso.coef != 0] ## tocno 11 nenul koeficijenata, a to je nepristran procjenitelj
## za stupnjeve slobode (vidi predavanja)

## obje metode daju puno bolje rezultate nego lin. regresija sa metodom najmanjih kvadrata (lambda=0)

## DZ - na trening skupu provesti best ss metodu, procijeniti testnu gresku odabranog modela na testnom skupu
## te usporediti sa ridge regresijom i lasso metodom

## PCR i PLS
library(pls)
set.seed(2)
pcr.fit <- pcr(Salary ~ ., data = Hitters, scale = TRUE, validation="CV")
summary(pcr.fit)
validationplot(pcr.fit, val.type = "MSEP")

set.seed(1)
pcr.fit <- pcr(Salary ~ ., data = Hitters, subset = train,
               scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")
which.min(pcr.fit$validation$PRESS) ## ali i manje se cini dovoljno dobro

pcr.pred <- predict(pcr.fit, x[test, ], ncomp = 5)
mean((pcr.pred - y.test)^2) ## [1] 142811.8 -> slicno kao ridge i lasso

## odabrani model je
pcr.fit <- pcr(y ~ x, scale = TRUE, ncomp = 5)
summary(pcr.fit)

## PLS vidjeti sami