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"]


library(randomForest)

?randomForest  ## dosta opcija
######## bagging (=RF uz m=p)
B=2000
set.seed(131)  
### buduci da stabla ovise o bootstrap uzorcima koji su 
##  slucajni, svaki put kada prilagodimo rf dobit
##  cemo drugaciji rezultat (ako ne koristimo set.seed)
bag.boston <- randomForest(medv ~ ., data = Boston,
                           subset = train, mtry = 12, 
                           xtest=Boston[-train, -13], 
                           ytest=Boston[-train, 13],
                           importance = TRUE, ntree=B)

### mtry=m (broj kovarijata koje biramo pri svakom dijeljenju)
### ntree=B (broj bootstrap stabala)
### nodesize -- minimalna velicina listova (veci broj, manja stabla)
##  default za regresiju je 5, za klasifikaciju 1
bag.boston

oob.procjene_bag=bag.boston$mse 
### vektor duljine B=broj stabala,
##  na i-tom mjestu je OOB procjena greske
##  za procjenitelj temeljen na prvih i stabala,
##  za i=1,...,B

test_err_bag=bag.boston$test$mse
### isto kao gore samo za testni skup umjesto OOB podataka

### Koliki B je dovoljno velik?
##  Tipicno pratimo dok se OOB procjena greske
##  ne stabilizira (vidi (8.24))


### crtamo procjene u ovisnosti o broju koristenih stabala
plot(1:B, test_err_bag, type = "l", 
     col="darkgreen", ylim=c(12, 16)) 
points(1:B, oob.procjene_bag, type="l")
abline(v=1000, lty=2)
## obje greske se stabiliziraju vec oko 1000 stabala
## ili
plot(bag.boston, ylim=c(12, 16))  
### crta OOB procjene greske



###### probajmo RF model s m=2,6,10
### za sve slucajeve crtamo OOB procjene greske
set.seed(134)
fit2=randomForest(medv ~ ., data = Boston,
                 subset = train, mtry = 2, 
                 importance = FALSE, ntree=B)
fit6=randomForest(medv ~ ., data = Boston,
                  subset = train, mtry = 6, 
                  importance = FALSE, ntree=B)
fit10=randomForest(medv ~ ., data = Boston,
                  subset = train, mtry = 10, 
                  importance = FALSE, ntree=B)

plot(1:B, fit2$mse, type="l", ylim=c(12,21), col=2)
points(1:B, fit6$mse, type="l", col=6)
points(1:B, fit10$mse, type="l", col=10)
points(1:B, bag.boston$mse, type="l", col=12)

legend("top", legend = paste("m=", c(2,6,10,12)), 
       col =  c(2,6,10,12), pch=16, horiz=TRUE,
       cex=0.5)

### m=6 se cini najbolji odabir, 
##testnu gresku onda mozemo procijeniti na testnom skupu
yhat.rf <- predict(fit6, newdata = Boston[-train, ])
sqrt(mean((yhat.rf - boston.test)^2)) ## 3.402226
mean(boston.test) ## 22.24704


### znacajnost varijabli
importance(fit6, type=2)
?importance
### ili bolje graficki
varImpPlot(fit6, type=2)
?Boston
### rm=prosjecan broj soba po nekretnini u podrucju
### lstat= the percentage of homeowners 
##  in the neighborhood considered lower class

### ove dvije kovarijate se cine daleko najznacajnije

### nox=nitrogen oxides 
##  concentration (parts per 10 million).
### dis=weighted mean of distances 
##  to five Boston employment centres.

prosjek=mean(Boston[train, "medv"])
### parcijalna zavisnost medv o 4 najznacajnije kovarijate
partialPlot(fit6, Boston[train,], 
            x.var="lstat", rug=TRUE)
abline(h=prosjek, lty=2)
### na x-osi su naznaceni decili podataka
##  van tog raspona ne mozemo procjene
##  uzimati ozbiljno
partialPlot(fit6, Boston[train,], 
            x.var="rm", rug=TRUE)
abline(h=prosjek, lty=2)

partialPlot(fit6, Boston[train,], 
            x.var="nox", rug=TRUE)
abline(h=prosjek, lty=2)
### Objasnjenje? Trebalo bi dodatno istraziti
##  (npr. vidjeti postoje li neke interakcije)

partialPlot(fit6, Boston[train,], 
            x.var="dis", rug=TRUE)
abline(h=prosjek, lty=2)