# Insurance data

library(randomForest)
library(gbm)
library(vip)
library(ggplot2)

TL <- read.csv("TL.csv", header=TRUE, sep=",", row.names=1)
dim(TL)

library(leaps)
LFACE <- log(TL$FACE)
LINCOME <- log(TL$INCOME)

TL2 <- TL[, -c(5,6)]
TL2$LINCOME <- with (TL2, LINCOME)
TL2$LFACE <- with(TL2, LFACE)

## bagging

B <- 10^3 # numero di alberi
p <- dim(TL2)[2]-1 # numero di covariate
train <- sample(1:nrow(TL2), 137) # training set
set.seed(1)
bag.tl <- randomForest(LFACE ~ ., data = TL2,
                           subset = train, 
                           ntree = B, 
                           mtry = p)
print(bag.tl)
y.test <- TL2[-train, 'LFACE']
round(importance(bag.tl), 2)

# %IncMSE (Percentuale di Incremento dell'Errore Quadratico Medio)

# - Questo metodo misura l’importanza della variabile attraverso la 
#   tecnica della permutazione (o permutation importance).
# - Funziona permutando i valori della variabile di interesse in ciascun albero e 
#   ricalcolando l’errore OOB (Out-Of-Bag error).
# - La differenza media tra l’errore OOB originale e l’errore OOB con la 
#   variabile permutata rappresenta l’importanza della variabile. 
#   Maggiore è l'aumento dell'errore, più la variabile è considerata importante.
# - %IncMSE è riportato come percentuale di incremento dell'errore quadratico medio.

# IncNodePurity (Incremento dell'Impurità del Nodo)

# - Questo metodo misura l’importanza di una variabile in base alla riduzione dell’impurità
#   dei nodi per i quali quella variabile è utilizzata come criterio di suddivisione.
# - Nella regressione, la metrica di impurità è data dalla somma dei residui quadrati 
#   (RSS, Residual Sum of Squares).
# - Per ogni nodo diviso da una variabile specifica, si calcola la riduzione dell'RSS; 
#   queste riduzioni vengono sommate lungo tutto l'albero per ottenere l'IncNodePurity della variabile.

p1 <- vip(bag.tl, num_features = p) + ggtitle("Variable Importance")
p1
yhat.bag <- predict(bag.tl, newdata = TL2[-train, ])
plot(yhat.bag, y.test)
abline(0, 1)
mean((y.test - yhat.bag )^2)  # errore sul test set
varImpPlot(bag.tl)


## rf

set.seed(1)
rf.tl <- randomForest(LFACE ~ ., data = TL2,
                       subset = train, 
                       ntree = B, 
                       mtry = sqrt(p))
print(rf.tl)
y.test <- TL2[-train, 'LFACE']
round(importance(rf.tl), 2)
p1 <- vip(rf.tl, num_features = p) + ggtitle("Variable Importance")
p1
yhat.rf <- predict(rf.tl, 
                    newdata = TL2[-train, ])
plot(yhat.rf, y.test)
abline(0, 1)
mean((y.test-  yhat.rf)^2)

# Two measures of variable importance are reported.
# The first is based upon the mean decrease of 
# accuracy in predictions on the out of bag samples
# when a given variable is permuted.
# The second is a measure of the total decrease 
# in node impurity that results from splits over 
# that variable, averaged over all trees 
# (this was plotted in Figure 8.9). 
# In the case of regression trees, the node impurity 
# is measured by the training RSS, and for classification 
# trees by the deviance. 
# Plots of these importance measures can be produced 
# using the varImpPlot() function.

round(importance(rf.tl), 2)
varImpPlot(rf.tl)


## boosting

set.seed(1)
TL2$MARSTAT <- as.factor(TL2$MARSTAT)
boost.tl <- gbm(LFACE ~ ., data = TL2[train, ],
                    distribution = "gaussian", n.trees = B,
                    interaction.depth = 1)
summary(boost.tl)

# - La funzione gbm() costruisce il modello usando alberi decisionali additivi per 
#   minimizzare una funzione di perdita, che può variare in base al tipo di problema 
#   (ad esempio, per la regressione è usata la devianza).

# - Ogni albero del modello contribuisce alla riduzione della funzione di perdita e, 
#   a ogni split dell’albero, la riduzione viene assegnata alla variabile usata per fare lo split.
# - Ogni split che riduce la perdita contribuisce quindi a determinare l’importanza della 
#   variabile: maggiore è la riduzione della funzione di perdita, 
#   maggiore è il contributo di quella variabile.

# - La riduzione della perdita ottenuta grazie a ciascuna variabile viene 
#   accumulata lungo tutti gli alberi del modello.

# - Per ogni split di ogni albero, se una variabile contribuisce alla riduzione 
#   della perdita, il valore di questa riduzione viene aggiunto al punteggio complessivo 
#   di quella variabile.

# - Dopo aver sommato le riduzioni di perdita per ciascuna variabile, 
#   questi valori vengono normalizzati rispetto alla somma totale delle riduzioni.

# - Questo significa che la somma delle influenze relative di tutte le variabili è pari al 100%.
# - La Relative Influence di ciascuna variabile è quindi espressa come una 
#   percentuale, che rappresenta l’influenza relativa della variabile rispetto alle altre.

# - Una variabile con una Relative Influence più alta ha contribuito in modo più 
#   significativo alla riduzione della perdita e, quindi, è considerata più importante 
#   per le previsioni del modello.
# - Le variabili con Relative Influence molto basse possono essere considerate poco 
#   influenti e, se necessario, eliminate dal modello per ridurre la complessità 
#   senza una grande perdita in performance.



yhat.boost <- predict(boost.tl, newdata = TL2[-train, ], n.trees = B)
mean((yhat.boost - y.test)^2)

### change lambda (parametro di shrinkage)
boost.tl <- gbm(LFACE ~ ., data = TL2[train, ],
                    distribution = "gaussian", n.trees = B,
                    interaction.depth = 1, 
                    shrinkage = 0.001, verbose = F)
yhat.boost <- predict(boost.tl, newdata = TL2[-train, ], n.trees = B)
mean(( y.test- yhat.boost)^2)
plot(boost.tl, i = "LINCOME") # partial dependence plots
plot(boost.tl, i = "AGE")
plot(boost.tl, i = "MARSTAT")


## simulazione finale all'aumentare del numero di alberi
# (potete far variare gli iperparametri a piacere, provateci!)

mse_bag <- c()
mse_rf <- c()
mse_boost <- c()
for (b in 1:B){
  # bagging
  bag <- randomForest(LFACE ~ ., data = TL2,
                                 subset = train, 
                                 ntree = b, 
                                 mtry = p)
  yhat.bag <- predict(bag, newdata = TL2[-train, ])
  mse_bag[b] <- mean((y.test-yhat.bag)^2)
  
  # random forest
  rf <- randomForest(LFACE ~ ., data = TL2,
                                subset = train, 
                                ntree = b, 
                                mtry = sqrt(p))
  yhat.rf <- predict(rf, newdata = TL2[-train, ])
  mse_rf[b] <- mean((y.test-yhat.rf)^2)
  
  # boosting
  boost <- gbm(LFACE ~ ., data = TL2[train, ],
                          distribution = "gaussian", 
                          n.trees = b,  # numero alberi
                          interaction.depth = 1,  # profondità albero
                          shrinkage = 0.1, # parametro shrinkage
                          verbose = F) 
  yhat.boost <- predict(boost, newdata = TL2[-train, ], n.trees = B)
  mse_boost[b] <- mean((y.test-yhat.boost)^2)
}

plot(mse_bag, type ="l", col = 1, ylab = "MSE", cex.lab =1.4, xlab = "Trees", ylim =c(2,4))
lines(1:B, mse_rf, col = 2)
lines(1:B, mse_boost, col =3)
legend(700,4,  lty =1, col =c(1,2,3), c("Bagging", "RF", "Boosting"))