############################################################### #Laboratorio 8- MAGGIO 2025#################################### #Market Bastket Analysis (Regole associative) ############################################################# ################################### ##Regole associative############### ######Market Basket Analysis ###### library(arules) ########Toy example contenuto in 'arules' data("Groceries") Groceries inspect(head(Groceries, 10)) groc.rules <- apriori(Groceries, parameter=list(supp=0.001, conf=0.5)) inspect(head(groc.rules, n=3, by="lift")) library(arulesViz) plot(groc.rules) plot(groc.rules, method="two-key plot") #eventualemnte interattivo per ispezionare punti o fare zoom plot(groc.rules, method="two-key plot", engine="interactive") ##order indica il numero di items presenti in una regola subrules<- head(groc.rules, n=10, by="lift") inspect(subrules) plot(subrules, method="graph") #i punti sono le regole, come da legenda: dimensione=support, colore= lift plot(subrules, method="graph", engine="interactive") #analogo a sopra, ma si possono spostare punti e prodotti nel grafico per una visualizzazione #personalizzata (dimensione=support, colore= lift) #altri "engine", da provare plot(subrules, method="graph", engine="visNetwork") plot(subrules, method="graph", engine="igraph") ########################################################### lastfm <- read.csv("lastfm.csv", stringsAsFactors=TRUE) ### osserviamo il data set lastfm[1:50,] dim(lastfm) str(lastfm) lastfm$user <- factor(lastfm$user) nlevels(lastfm$user) ###15000 utenti nlevels(lastfm$artist) ####1004 artisti # creo una lista con 'split' # creo lista, un elemento per ogni utente, in cui gli elementi sono l'insieme degli artisti scelti da quell'utente playlist <- split(x=lastfm[,"artist"],f=lastfm$user) playlist[1:2] #### cosa ascoltano i primi 2 utenti (1 e 3) length(playlist) #cerco eventuali duplicati, ad esempio in (sono gli unici!) playlist$`6980` #qui si vede che m.i.a. compare due volte, una in posizione 7 e l'altra in posizione 12 playlist$`9753` #qui è james brown ascoltato due volte, uno in posizione 8 e l'altro in posizione 17 ## rimuovo duplicati perché un artista potrebbe essere nominato più di una volta, ed è importante rimuovere i duplicati prima di creare la matrice di incidenza playlist <- lapply(playlist,unique) # trasformo la lista in oggetto di classe 'transactions' (pacchetto arules) playlist <- as(playlist,"transactions") ###lista di 'transazioni' playlist # supporto: numero di volte che un artista compare nella playlist dei 15000 utenti itemFrequency(playlist)[1:10] #li ordino in maniera decrescente per vedere i 10 piu' ascoltati sort(itemFrequency(playlist),decreasing=T)[1:10] #grafico itemFrequencyPlot(playlist) #visualizziamo le frequenze solo di quelli con supporto almeno 0.08 itemFrequencyPlot(playlist,support=.08) # costruiamo le regole con algoritmo apriori musicrules <- apriori(playlist,parameter=list(support=.01,confidence=.5)) inspect(musicrules) #support= P(A e B), confidence= P(B|A), coverage= P(A), lift = P(B|A)/P(B) ###T.I., pseudonimo di Clifford Joseph Harris Jr. (Atlanta, 25 settembre 1980), è un rapper e attore statunitense, leader del gruppo rap chiamato P$C # filtro con il lift > 5 e ordino secondo la fiducia inspect(subset(musicrules, subset=lift > 5)) inspect(sort(subset(musicrules, subset=lift > 5), by="confidence")) #library(arulesViz) plot(musicrules) plot(musicrules, method = "graph") plot(musicrules, method = "graph", engine="interactive") plot(musicrules, method = "graph", engine="igraph") plot(musicrules, method = "graph", engine="visNetwork") #utile per vedere ogni regola singolarmente ####################################################################### ####Modelli aggiuntivi per segmentazione################################################# ####################################################################### # Model-based clustering ###################################################################### ###################################################################### ## Modelli a misture finite ###################################################################### ## Misture di distribuzioni normali library(flexclust) #per i dati library(mclust) data("vacmot") #dataset da "flexclust" vacmet <- vacmotdesc[, c("Obligation", "NEP", "Vacation.Behaviour")] vacmet <- na.omit(vacmet) ###rimuoviamo dati mancanti pairs(vacmet, pch = 19, col = rgb(0, 0, 0, 0.2)) ####visualizziamo i dati vacmet.m18 <- Mclust(vacmet, G = 1:8) ####stimiamo tutti i 14 modelli (default). Il miglior modello rispetto al BIC viene selezionato. G = numero di segmenti da estrarre (da 1 a 8) ##Miglior modello VEE: ellissoidale, uguale forma e orientamento vacmet.m18.equal <- Mclust(vacmet, G = 1:8, modelNames = c("EEI", "EII", "EEE")) #####in alternativa, stimiamo solo modelli in cui le matrici di covarianza hanno uguale forma, dimensione e orientamento tra i segmenti. ###EEI: Diagonal, equal volume and shape ### EII: Spherical, equal volume ## EEE: Ellipsoidal, equal volume, shape, and orientation vacmet.m18 vacmet.m18.equal plot(vacmet.m18, what = "classification") plot(vacmet.m18.equal, what = "classification") plot(vacmet.m18, what = "uncertainty") plot(vacmet.m18, what = "BIC") plot(vacmet.m18.equal, what = "uncertainty") plot(vacmet.m18.equal, what = "BIC") ###In entrambi i modelli mistura selezionati, le matrici di covarianza hanno lo stesso orientamento e la stessa forma. Questo implica che la struttura di correlazione tra le variabili è la stessa nei vari segmenti. ### ######## Misture di regressioni############################# install.packages("flexmix") library("flexmix") vacmet[, c("Obligation", "NEP")] <- scale(vacmet[, c("Obligation", "NEP")]) #####standardizziamo le variabili Obligation e NEP (migliore visualizzazione dei risultati) ###stimiamo prima di tutto un modello di regressione lineare vacmet.lm <- lm(Vacation.Behaviour ~ Obligation + NEP, data = vacmet) summary(vacmet.lm) ####proviamo a stimare una mistura di modelli di regressione, per capire se la relazione può essere diversa tra i segmenti set.seed(1234) vacmet.m15 <- stepFlexmix(Vacation.Behaviour ~ ., data = vacmet, k = 1:5, nrep = 10, verbose = FALSE, control = list(iter.max = 1000)) ####selezioniamo il modello migliore secondo il BIC vacmet.m2 <- getModel(vacmet.m15) vacmet.m2 summary(refit(vacmet.m2)) ####la funzione refit ci serve per ottenere gli standard error delle stime. Refit usa la soluzione ottenuta con l'algoritmo EM e usa un ottimizzatore per ottenere l'informazione relativa allo standard error. par(mfrow = c(1, 2)) plot(Vacation.Behaviour ~ Obligation, data = vacmet, pch = 20, col = clusters(vacmet.m2)) abline(parameters(vacmet.m2)[1:2, 1], col = 1, lwd = 2) abline(parameters(vacmet.m2)[1:2, 2], col = 2, lwd = 2) plot(Vacation.Behaviour ~ NEP, data = vacmet, pch = 20, col = clusters(vacmet.m2)) abline(parameters(vacmet.m2)[c(1, 3), 1], col = 1, lwd = 2) abline(parameters(vacmet.m2)[c(1, 3), 2], col = 2, lwd = 2)