library(ISLR2)
?Boston
head(Boston)
### Y=medv
dim(Boston)
# [1] 506  13
### dakle, n=506, p=12


### podijelit cemo uzorak na dva jednaka dijela
##  jedan za trening, drugi za test
set.seed(1221)
train <- sample(1:nrow(Boston), nrow(Boston) / 2)
boston.test=Boston[-train, "medv"]


######### gradient boosting
library(gbm)
?gbm  
### distribution (funkcija gubitka),  
### n.trees (=M), def=100
### interaction.depth (=J-1), def=1 ("decision stumps")
### shrinkage (=lambda), def=0.1
### bag.fraction (stochastic GB), def=0.5

### buduci da koristimo bag.fraction<1, rezultati 
### ovise o generiranim poduzorcima u svakom koraku
set.seed(1)
M=4000
### prvo cemo prilagoditi aditivan model (J=2), tj.
##  rezultirajuca funkcija je oblika 
##  f_hat(x)=h_1(x_1)+...+h_p(x_p)
boost.boston <- gbm(medv ~ ., data = Boston[train, ],
                    distribution = "gaussian", n.trees = M,
                    interaction.depth = 1, cv.folds = 5)

### crtamo CV procjenu greske s obzirom na broj
##  koristenih stabala
plot(boost.boston$cv.error, type="l", ylim=c(15,25))
M_min=which.min(boost.boston$cv.error)
points(M_min, boost.boston$cv.error[M_min], pch=16)

### smanjimo "shrinkage" parametar na 0.05
set.seed(4)
boost.boston_2 <- gbm(medv ~ ., data = Boston[train, ],
                    distribution = "gaussian", n.trees = M,
                    interaction.depth = 1, cv.folds = 5,
                    shrinkage = 0.05)

points(boost.boston_2$cv.error, type="l", col="darkgreen")
M_min=which.min(boost.boston_2$cv.error)
points(M_min, boost.boston_2$cv.error[M_min], 
       pch=16, col="darkgreen")
### probajmo i 0.01
set.seed(7)
boost.boston_3 <- gbm(medv ~ ., data = Boston[train, ],
                      distribution = "gaussian", 
                      n.trees = M,
                      interaction.depth = 1, cv.folds = 5,
                      shrinkage = 0.01)

points(boost.boston_3$cv.error, type="l", col="darkblue")
M_min=which.min(boost.boston_3$cv.error)
points(M_min, boost.boston_3$cv.error[M_min], 
       pch=16, col="darkblue")

### kako smanjujemo shrinkage algoritam "sporije uci"
##  te za optimalni M tipicno dobivamo
##  manju testnu gresku

### ipak, povecava se i optimalni M sto nekada moze postati
##  prakticno nemoguce za prilagoditi

### tipicno, izabiremo shrinkage sto manji mozemo
##  s obzirom na racunalne mogucnosti

### mi cemo ostati s 0.05


### probajmo mijenjati J, tj. dopustati interakcije izmedu
##  varijabli
### sada je f_hat oblika
##  f_hat(x)=sum_j h_j(x_j) + sum_{j,k} h_jk(x_j,x_k)
##  +...+ sve do funkcija koje ovise o najvise J-1 kovarijata

set.seed(10)
boost.boston_4 <- gbm(medv ~ ., data = Boston[train, ],
                      distribution = "gaussian", n.trees = M,
                      interaction.depth = 4, cv.folds = 5,
                      shrinkage = 0.05)

plot(boost.boston_4$cv.error, type="l", 
     col="darkorange", ylim=c(13,25))
M_min=which.min(boost.boston_4$cv.error)
points(M_min, boost.boston_4$cv.error[M_min], 
       pch=16, col="orange")

### dodajmo i CV gresku za aditivan GB model
points(boost.boston_2$cv.error, type="l", col="darkgreen")
M_min=which.min(boost.boston_2$cv.error)
points(M_min, boost.boston_2$cv.error[M_min], 
       pch=16, col="darkgreen")
### buduci da koristimo kompleksnije bazne procjenitelje
##  algoritam je brze konvergirao, tj. ranije se
##  poceo previse prilagodavati podacima

## usporedimo i testne greske
M_adit=which.min(boost.boston_2$cv.error)  ## 1519
yhat.boost <- predict(boost.boston_2,
                      newdata = Boston[-train, ], 
                      n.trees = M_adit)
sqrt(mean((yhat.boost - boston.test)^2)) ## [1] 3.839342


M_4=which.min(boost.boston_4$cv.error)  ## 301
yhat.boost <- predict(boost.boston_4,
                      newdata = Boston[-train, ], 
                      n.trees = M_4)
sqrt(mean((yhat.boost - boston.test)^2)) ## 3.470567

### dakle, dodavanjem interakcija se dobije
## bolji model

### sa slucajnim sumama smo postigli 3.402226

### Boosting cesto daje nesto bolje rezultate nego RF, 
##  ali RF ima puno manje hiperparametara za odabrati


########## intepretacija -- slicno kao kod slucajnih suma

par(mfrow=c(1,1))
## aditivan GB model
summary(boost.boston_2, n.trees = M_adit)  
### crta variable importance plot 
##  slicni rezultati kao kod RF

## GB model s interakcijama
summary(boost.boston_4, n.trees = M_4)  
### slicno, ali povecao se utjecaj varijabli
##  lstat, dis i age 
### to moze biti posljedica
##  znacajne interakcije/ija izmedu tih varijabli

##### crtamo grafove parcijalne zavisnosti
### usporedivat cemo rezultate oba GB modela

### u aditivnom modelu, buduci da je
##  f_hat(x)=h_1(x_1)+...+h_p(x_p),
##  graf parcijalne ovisnosti o podskupu varijabli A
##  je f_A(x_A) = sum_{j iz A} h_j(x_j) + const.
##  (vidi predavanja)


### za orijentaciju
prosjek=mean(Boston[train, "medv"])
prosjek ## [1] 22.81858

?plot.gbm

### rm
plot(boost.boston_2, i = "rm", n.trees = M_adit)
### nema opcije za da se naznace decili podataka, ali 
##  raspon na x-osi je vec odreden s podacima:
quantile(Boston[train, "rm"],probs = (0:10)/10)

plot(boost.boston_4, i = "rm", n.trees = M_4)

## lstat
plot(boost.boston_2, i = "lstat", n.trees = M_adit)
plot(boost.boston_4, i = "lstat", n.trees = M_4)

## dis
plot(boost.boston_2, i = "dis", n.trees = M_adit)
### pogledat cemo interakcije kasnije
plot(boost.boston_4, i = "dis", n.trees = M_4)

## crim
plot(boost.boston_2, i = "crim", n.trees = M_adit)
plot(boost.boston_4, i = "crim", n.trees = M_4)

## nox
plot(boost.boston_2, i = "nox", n.trees = M_adit)
plot(boost.boston_4, i = "nox", n.trees = M_4)

## age
plot(boost.boston_2, i = "age", n.trees = M_adit)
## kasnije cemo pogledati interakcije
plot(boost.boston_4, i = "age", n.trees = M_4)
### parcijalna ovisnost medv o age se ne cini monotona:
##  za velike vrijednosti od age ponovno pocinje rasti

### sve se cini u skladu s ocekivanjima

### pogledajmo i potencijalne interakcije

library(viridis)
## ovako to izgleda u aditivnom modelu 
## dakle, bez interakcija
## funkcija parcijalne ovisnosti je do na 
## konstatnu oblika
## h_1(x_lstat) + h_2(x_dis) gdje su to,
## opet do na konstatnu,
## zapravo funkcije parcijalne ovisnosti 
## za 1d slucaj (crtali smo ih gore)
plot(boost.boston_2, i = c("dis","lstat"), 
     level.plot = TRUE, n.trees = M_adit)
plot(boost.boston_2, i = c("dis","lstat"), 
     level.plot = FALSE, n.trees = M_adit)
### u aditivnom modelu, efekt kovarijate lstat na medv 
##  s obzirom na vrijednost od dis je isti do 
##  na konstantu 

### pogledajmo i model s interakcijama
plot(boost.boston_4, i = c( "dis","lstat"), 
     level.plot = TRUE, n.trees = M_4)
plot(boost.boston_4, i = c("dis","lstat"), 
     level.plot = FALSE, n.trees = M_4)
### efekt dis na medv je drugaciji s obzirom
##  na to je li lstat malen i li velik:
##  za male vrijednsti lstat, utjecaj od dis
##  puno brze raste kako dis ide prema 0 (logicno?)

### lstat i age 
plot(boost.boston_4, i = c("age", "lstat"), 
     level.plot = TRUE, n.trees = M_4)
plot(boost.boston_4, i = c("age", "lstat"), 
     level.plot = FALSE, n.trees = M_4)

### za sve osim jako malih vrijednosti lstat
##  parcijalni utjecaj na medv opada kako
##  se povecava age
### za jako mali lstat, imamo porast utjecaja
## za jako velike vrijednosti od age

### age i dis
plot(boost.boston_4, i = c("age", "dis"), 
     level.plot = TRUE, n.trees = M_4)
plot(boost.boston_4, i = c("age", "dis"), 
     level.plot = FALSE, n.trees = M_4)
### Osobno ne vidim ovdje neku posebno znacajnu
##  interakciju buduci da je utjecaj age na medv 
##  do na konstantu isti za sve vrijednosti dis.