###############################################
#Laboratorio 6- 
#Gradient Boosting per regressione- Aprile 2025
###############################################
#
# dataset
dati <- read.csv("movies.csv", stringsAsFactors=TRUE)
str(dati)
#cancello le colonne di indicatori
dati<-dati[,-c(1,2)]
#trasformo la variabile release_date in un formato "data"
dati$release_date <- as.Date(dati$release_date, "%d/%m/%Y")

str(dati)

# variabile risposta
summary(dati$vote_average)
boxplot(dati$vote_average, col="orange", ylim=c(0,10), main="Movies", ylab="Rating")
#
hist(dati$vote_average, col="orange", main="Movies", xlab="Rating")
#
#variabili esplicative
summary(dati)

#da notare che la release_date ha come minimo un giorno del 1916 e come massimo un giorno del 2016: circa 100 anni

#
#grafici delle variabili quantitative
summary(dati[,c(1,2,5,6)])
par(mfrow=c(2,2))
for(i in c(1,2,5,6)){
  hist(dati[,i], col="orange", main=paste(colnames(dati)[i]), xlab="")
}


#trasformo le esplicative quantitative in scala logaritmica
dati$budget <- log(dati$budget)
dati$popularity <- log(dati$popularity)
dati$revenue <- log(dati$revenue)

summary(dati[,c(1,2,5,6)])

par(mfrow=c(2,2))
for(i in c(1,2,5,6)){
  hist(dati[,i], col="orange", main=paste(colnames(dati)[i]), xlab="")
}
#torno al panel originale
par(mfrow=c(1,1))

#trasformo la variabile release_date in numerica (per poter usare gbm)
dati$release_date<-as.numeric(dati$release_date) 
####numero di giorni dal 1 gennaio 1970

# train and test
set.seed(1)
train = sample (1:nrow(dati), 0.7*nrow(dati))
#### insieme da cui estrarre, da cui prendo il 70%
dati.train=dati[train ,] 
dati.test=dati[-train ,] 
## 969 osservazioni, 25 variabili 


# Gradient Boosting #######
#install.packages("gbm")
library (gbm)

# primo boosting per regressione
boost.movies=gbm(vote_average ~ .-vote_classes, data=dati.train, 
                 distribution="gaussian", n.trees=5000, interaction.depth=1)
#
#per il grafico. panel unico (se non l'avevo fatto prima)
par(mfrow=c(1,1))
#
#grafico con gli errori di previsione sull'insieme di stima
plot(boost.movies$train.error, type="l", ylab="training error")
#si nota che e', ovviamente, sempre decrescente al crescere delle iterazioni
#
#
#grafico di importanza delle variabili
summary(boost.movies) 
#il grafico non è bello: modifico i parametri grafici.
#
#Allargo lo spazio a sinistra del grafico per farci stare i nomi delle variabili
#
#creo vettore con parametri di default
mai.old<-par()$mai
mai.old
#creo vettore con nuovi parametri uguali ai vecchi
mai.new<-mai.old
#sostituisco il parametro relativo allo spazio a sinistra
mai.new[2] <- 2.1 
mai.new
#modifico i parametri grafici
par(mai=mai.new)
summary(boost.movies, las=1) 
#las=1 mette i nomi orizzontali sulle y
summary(boost.movies, las=1, cBar=10) 

#las=1 mette i nomi orizzontali sulle y
#cBar definisce quante variabili disegnare
#torno ai parametri grafici originali
par(mai=mai.old)



# previsione sul test set
#previsione del modello per ogni iterazione del gbm
#yhat.boost=predict(boost.movies, newdata=dati.test)

yhat.boost=predict(boost.movies, newdata=dati.test, n.trees=1:5000)

#### calcolo l'errore di previsione per ciascun passo del gbm
#uso apply che e' un modo veloce di fare ciclo for
#il primo argomento e' la matrice da usare (le previsioni), il secondo e' 
#la dimensione su cui fare l'operazione ("2" ovvero "per colonna") , il 
#terzo e' la funzione da calcolare che e' l'errore quadratico medio


err = apply(yhat.boost, 2, function(pred) mean((dati.test$vote_average - pred)^2))
#
plot(err, type="l")

# confronto tra errori (train e test)
plot(boost.movies$train.error, type="l")
lines(err, type="l", col=2)
#disegno il minimo dell'errore di previsione nel test set
best=which.min(err)
abline(v=best, lty=2, col=4)
#
min(err) #minimo errore di previsione


# secondo boosting - alberi con più profondità
boost.movies=gbm(vote_average ~ .-vote_classes, data=dati.train, 
                 distribution="gaussian", n.trees=5000, interaction.depth=4)
#
#grafico degli errori
plot(boost.movies$train.error, type="l")
#
#Modifica parametri grafici per importanza variabili
#
#uso il vettore di parametri grafici 'mai' ottenuto prima
#modifico i parametri grafici
par(mai=mai.new)
#summary(boost.movies, las=1) 
summary(boost.movies, las=1, cBar=10) 

#torno ai parametri grafici originali
par(mai=mai.old)
#
# previsione su test set
yhat.boost=predict(boost.movies ,newdata=dati.test,n.trees=1:5000)
err = apply(yhat.boost,2,function(pred) mean((dati.test$vote_average-pred)^2))
plot(err, type="l")

# confronto tra errori
plot(boost.movies$train.error, type="l")
lines(err, type="l", col=2)
best=which.min(err)
abline(v=best, lty=2, col=4)
min(err)


# terzo boosting - learning rate più basso

boost.movies=gbm(vote_average ~ .-vote_classes, data=dati.train, 
                 distribution="gaussian", n.trees=5000, interaction.depth=1, shrinkage=0.01)
plot(boost.movies$train.error, type="l")

par(mai=mai.new)
#summary(boost.movies, las=1) #las=1 mette i nomi orizzontali sulle y
summary(boost.movies, las=1, cBar=10) #cBar definisce quante variabili disegnare
#torno ai parametri grafici originali
par(mai=mai.old)

# previsione su test set
yhat.boost=predict(boost.movies ,newdata=dati.test,n.trees=1:5000)
err = apply(yhat.boost,2,function(pred) mean((dati.test$vote_average-pred)^2))
plot(err, type="l")

# confronto tra errori
plot(boost.movies$train.error, type="l")
lines(err, type="l", col=2)
best=which.min(err)
abline(v=best, lty=2, col=4)
min(err)


# quarto boosting - learning rate e alberi più profondi
boost.movies=gbm(vote_average ~ .-vote_classes ,data=dati.train, 
                 distribution="gaussian",n.trees=5000, interaction.depth=4, shrinkage=0.01)

plot(boost.movies$train.error, type="l")
#

par(mai=mai.new)
#summary(boost.movies, las=1) 
summary(boost.movies, las=1, cBar=10) 


#torno ai parametri grafici originali
par(mai=mai.old)

# previsione su test set
yhat.boost=predict(boost.movies ,newdata=dati.test,n.trees=1:5000)
err = apply(yhat.boost, 2, function(pred) mean((dati.test$vote_average-pred)^2))
plot(err, type="l")


#####volendo posso usare anche un'altra misura di errore (MAE)
mae = apply(yhat.boost, 2, function(pred) mean(abs(dati.test$vote_average-pred)))



# confronto tra errori
plot(boost.movies$train.error, type="l")
lines(err, type="l", col=2)
best=which.min(err)
abline(v=best, lty=2, col=4)

err.boost= min(err)
mae.boost= min(mae)


boost.movies
# effetti marginali
plot(boost.movies, i.var=1, n.trees = best)
plot(boost.movies, i.var=2, n.trees = best)
plot(boost.movies, i.var=5, n.trees = best)
plot(boost.movies, i.var=c(1,5), n.trees = best) #bivariata
#
plot(boost.movies, i.var=3, n.trees = best) # variabile indicatrice
plot(boost.movies, i.var=6, n.trees = best)

plot(boost.movies, i=23, n.trees = best)# variabile qualitativa
plot(boost.movies, i=17, n.trees = best) #variabile che non ha effetto

#
#
#

######Altri modelli da provare
####Modello lineare############

m1 <- lm(vote_average~.-vote_classes, data=dati.train)

summary(m1)
#####################
# Stepwise Regression
m2 <- step(m1, direction="both")
summary(m2)


#Previsione
p.lm <- predict(m2, newdata=dati.test)
err.lm <-mean ((p.lm-dati.test$vote_average)^2)
err.lm
err.boost
1-err.boost/err.lm #tasso di miglioramento del boosting rispetto a lm
########

mae.lm <-mean (abs(p.lm-dati.test$vote_average))
mae.lm
mae.boost
1-mae.boost/mae.lm

AIC(m2)
#######################
#######GAM#############
dati.train[,c(3,7, 10:24)]= lapply(dati.train[,c(3,7, 10:24)],factor)
dati.test[,c(3,7, 10:24)]= lapply(dati.test[,c(3,7, 10:24)],factor)

library(gam)

###################
##Stepwise GAM
#Modello lineare (df=1)
g3 <- gam(vote_average~.-vote_classes, data=dati.train)

#Effetti lineari 
par(mfrow=c(3,5))
plot(g3, se=T) 

#Selezione stepwise
#df=1 modello lineare

sc = gam.scope(dati.train[,-8], response=8, arg=c("df=2","df=3","df=4"))
g4<- step.Gam(g3, scope=sc, trace=T)
summary(g4)

AIC(g4)

par(mfrow=c(3,5))
plot(g4, se=T)

# qualche grafico
par(mfrow=c(1,1))
plot(g4, se=T, ask=T)


#Previsione
dati.test[,c(3,7, 10:24)]= lapply(dati.test[,c(3,7, 10:24)],factor)
p.gam <- predict(g4,newdata=dati.test)     
err.gam <- mean((p.gam-dati.test$vote_average)^2)
err.gam