####################################################################
##### Laboratorio 4- MODELLI PER DATI ORDINALI- Aprile 2025######
###############################################################################

###Librerie###

library(plyr)
library(MASS)
##############################################
###Analisi della qualità dei vini#############
vino <- read.csv("wine.csv", row.names="X")
str(vino)

summary(vino)

##variabile risposta
# Analisi esplorativa ----------------------------------------------

barplot(table(vino$quality),  xlab="", main="quality", col="darkblue")

# ricodifica-trasformo in variabile qualitativa
vino$quality <- as.factor(vino$quality) 

#accorpamento di livelli 
vino$quality <- revalue(vino$quality, c("3"="Pessimo", "4"="Pessimo", "5"="Cattivo", "6"="Medio",
                                        "7"="Buono", "8"="Ottimo", "9"="Ottimo"))
str(vino$quality)

# ordinamento dei valori (fattore ordinato)
vino$quality <- ordered(vino$quality, levels=c("Pessimo", "Cattivo", "Medio", "Buono", "Ottimo"))
str(vino$quality)
summary(vino$quality)

par(mfrow=c(1,1))
barplot(table(vino$quality), xlab="", main="quality", col="darkblue")
#

###boxplot di variabili per qualità

par(mfrow=c(3,4))
for(i in c(1:11))  boxplot(vino[,i]~vino$quality, outline=F, names=c("P","C","M","B","O"), main=colnames(vino[i]))


#boxplot di variabili per colore

par(mfrow=c(3,4))
for(i in c(1:11))  boxplot(vino[,i]~vino$type, outline=F, names=c("Bianco","Rosso"), main=colnames(vino[i]))


############################################
# Regressione modello con logit cumulati ###
############################################

?polr  # package MASS

# primo modello
mod <- polr(quality~., data=vino)
summary(mod)


#stepwise AIC

mod2<- step(mod)
summary(mod2)
#
?step
#modifico il criterio di informazione e uso BIC
modstepbic<- step(mod, k=log(nrow(vino)))
summary(modstepbic)
#questo modello ha eliminato anche la variabile 'chlorides'


# classificazione
summary(mod2)
p.mod2<- predict(mod2)

#matrice confusione
t1<-table(p.mod2, vino$quality) 
t1

#tasso di corretta classificazione complessivo (da confrontare con una classificazione casuale, ma qui devo considerare 5/25, essendo le categorie 5. Non è una tabella 2x2. Ottengo circa 0.5 che è da confrontare con 0.2)
a <-sum(diag(t1))/sum(t1) 
a

##############################################################################
####################################
####Analisi del rating di film (dataset di recensioni)
###############################

dati <- read.csv("movies.csv", stringsAsFactors=TRUE)
str(dati)
###Analisi esplorativa del dataset
#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 alcune 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 
dati$release_date<-as.numeric(dati$release_date)



library(MASS)

# preprocessing, creazione di variabile ordinale 

dati$vote_classes <- ordered(dati$vote_classes, levels=c("0 a 5", "5 a 6", "6 a 7", "7 a 8", "8 a 10"))

str(dati$vote_classes)

barplot(table(dati$vote_classes), space=rep(0,5), xlab="", main="Rating", col="blue")


#per comodità di interpretazione trasformo le data di uscita del film in anni (con i decimali; attualmente è misurata in giorni).
#per normalizzare su 100 anni basta contare gli anni dal minimo 
#e dividere per il numero di giorni che ci sono ogni anno=365.
#
#In totale ci sono circa 100 anni (dal 1916 al 2016)

dati$release_date_years<- (dati$release_date-min(dati$release_date))/365

summary(dati$release_date_years)



# regressione con logit cumulati (variabile vote_classes)

######stepwise backward
ordreg <- polr(vote_classes ~ .- release_date-vote_average, data=dati)
summary(ordreg)

ordreg1<- step(ordreg)
summary(ordreg1)