Scarica BDA - Business Data Analytics (Prof. Ieva) e più Dispense in PDF di Modelli Stocastici E Analisi Dei Dati solo su Docsity!
DISPENSA DI BUSINESS DATA ANALYTICS (PARTE 1)
GENERALI
rm (list = ls ()) cat (“\014”) graphics.off () #Apertura delle libraries: library (GGally) #Per visualizzazione statistica avanzata library (rgl) library (tidyverse) #Per avere la manipolazione e visualizzazione dei dati library (car) #Per fare diagnostica dei modelli library (MASS) #Per fare modelli statistici library (mvtnorm) library (pracma) library (ggfortify) library (cluster) library (ggplot2) library(boot) library (ISLR) library (car) library (pls) library (glmet) library(titanic) library(PRROC) library (tree) library(randomForest) #Per cercare le informazioni sul funzionamento di una funzione o sul significato di un dataset: help (nome_funzione o nome_dataset) #Per caricare un dataset: dati = read.csv (“data.txt”) #assume che i campi siano separati da virgole dati = read.csv (“data.txt”, sep = “;”) o read.csv2 (“data.txt”) #campi separati da punto e virgola dati = read.csv (“data.txt”, sep = “;”) o dati = read.table (“data.txt”, header = TRUE) #campi separati da spazio #Per visualizzare il dataset: view (dati) #Per visualizzare i valori di una variabile del dataset: dati$nome_variabile #Per ottenere le prime righe del dataset: head (dati) #Per ottenere il numero di osservazioni (righe del dataset) e il numero di variabili (colonne del dataset): dim (dati) #In questo caso abbiamo 90 osservazioni (90 righe) e 4 variabili (4 colonne) #Per ottenere il nome delle variabili: names (dati)
#Per avere il numero di osservazioni e il numero e la tipologia di variabili: str (dati) #Per simulare un dataset di dati trivariati da una distribuzione gaussiana avendo il vettore delle medie e la matrice di covarianza: nobs = 200 sig = rbind (c(9, 1, 1), c(1, 4, 1), c(1, 1, 1)) mu = c(0, 2, 3) X = rmvnorm(nobs,mu,sig) sig #Per “pulire” il dataset: #Per eliminare la variabile X dal dataset dati : dati$X = NULL #Elimino la colonna X dati ( ,-X) #Per eliminare la riga Y: dati (-Y, ) #Elimino la riga Y remove (dati$Y) #Per eliminare tre righe di cui non conosco il numero, ma solo, ad esempio, il contenuto della colonna “names”: dati$names = row.names(dati) dati = dati [dati$names != “nome_dato”, ] dati = dati [dati$names != “nome_dato”, ] dati = dati [dati$names != “nome_dato”, ] #Nella rappresentazione grafica, se voglio che due grafici (due istogrammi) abbiano la stessa scala: x11() par (mfrow = c(1,2)) for (i in 1:2) { hist(dati[,i], main = names(dati)[i], xlab = names(dati)[i], prob = T, col = i, ylim = c(0,0.5)) }
PARTE 1
#VETTORI:
#Per creare vettori di numeri: v = c (NUM_1, NUM_2) v = seq (1, N, by = X) #Numeri equispaziati di X, da 1 a N v = seq (1,N, length.out = X) # Numeri da 1 a 9, equispaziati in modo che il vettore risultante abbia X elementi
#Per avere due istogrammi nel medesimo device grafico: x11 () par (mfrow = c(1, 2)) for (i in 1:2) { hist (dati [,i], main = names (dati)[i], xlab = names (dati)[i], prob = T, col = i)} #Se voglio che i due grafici abbiano la stessa scala, aggiungo che ylim = c (0, 0.1) #Per disegnare gli istogrammi della variabile “X” a seconda del genere e confrontarli: x11 () par (mfrow=c(2,1)) hist (dati$X [dati$Sesso ==’F’], xlim=range(X), col = ‘pink’) hist (dati$X [dati$Sesso ==’M’], xlim=range(X), col=’lightblue’) #Per rappresentare il BOXPLOT VERTICALE: x11 () boxplot (dati$X, main = “Boxplot”, ylab = “Livelli”, col = “gold”) abline (h = Q3+1.5IQR(dati$X), col = ‘red’, lwd=2) abline (h = Q1-1.5IQR(dati$X), col = ‘red’, lwd=2) abline (h=median(dati$X), col = ‘blue’, lty=2) #lty = 1: linea continua e lty = 2: linea tratteggiata #Per rappresentare il BOXPLOT ORIZZONTALE: x11 () boxplot (dati$X, horizontal = TRUE, main = “Boxplot”, xlab = “Livelli”) #Per rappresentare il boxplot della variabile “X” a seconda del genere e confrontarli: x11 () boxplot(X ~ Sesso, col = c(‘pink’, ‘lightblue’), main = ‘Boxplot ‘) #Per eliminare la presenza di eventuali outliers: length (dati$X) q1 = quantile (dati$X, 0.25) newdati$X = dati$X [which (dati$X > q1 – 1.5*range_interquartile)] length (newdati$X) quantile (newdati$X) mean (newdati$X) sd (newdati$X) x11 () hist (newdati$X, prob=TRUE, main=’Istogramma Density’, xlab = ‘Livelli’, ylab=’Densità’) boxplot(newdati$X, horizontal=FALSE, main=”Boxplot Density”, ylab=”Livelli”)
#ANALISI DESCRITTIVA E RAPPRESENTAZIONE GRAFICA:
#Per le VARIABILI CATEGORICHE: #Per rendere categorica (factor) una variabile carattere (chr): dati$X = factor (dati$X) femmine_X = dati$X [dati$Sesso == ‘F’] maschi_X = dati$X [dati$Sesso == ‘M] #Per costruire una tabella con le frequenze assoluta e relativa delle variabili categoriche: f_ass = table (dati$X) #Otteniamo le frequenze assolute f_rel = f_ass / length (dati$X) #Otteniamo le frequenze relative f_rel = prop.table (dati$X) barplot (f_rel, col = c ( ‘pink’, ‘lightblue’ ), ylab = ‘Frequenza relativa’ ) table (dati$X, dati$Y)
pie(f_ass, col = c (‘pink’,’lightblue’)) #Per visualizzare, ad esempio, una variabile numerica in relazione a due variabili categoriche: ggpairs(dati, aes(alpha = 0.4, color = X) #Guardiamo solo X come variabile categorica
#VARIABILI ALEATORIE E STRUMENTALI:
#Distribuzione NORMALE dati la media (mu) e la deviazione standard (sigma). #funzione di RIPARTIZIONE: t = seq (mu – 4sigma, mu + 4sigma, length.out = 2000 ) #vettore equispaziato P = pnorm (t, mean = mu, sd = sigma) t = seq (- 6 , 6 , length.out = 2000 ) P = pt (t, df = gradi) #funzione DENSITÀ f = dnorm (t, mean = mu, sd = sigma) f = dt (t, df = gradi) #funzione QUANTILE #la funzione quantile ha per dominio l’intervallo [0, 1] qualunque sia la distribuzione q = seq (1e-6,1-1e-6, length.out = 1000 ) Q = qnorm(q, mean = mu, sd = sigma) q = seq (1e-6,1-1e-6, length.out = 1000 ) Q = qt (q, df = gradi) #PLOT x11 () par (mfrow = c (1, 3)) plot (t, P, col = “forestgreen”, type = “l”, lwd = 2, main = “Ripartizione”) abline (h=c(0,1), col=”black”, lwd=1, lty=2) plot (t,f, col = “red”, type = “l”, lwd = 4, main = “Densità”) #type = ‘l’, la linea deve essere fatta di pallini e linee abline (v = mu, col = “black”, lwd = 2) plot (q,Q, col = “blue”, type = “l”, lwd = 2, main = “Quantile”) abline (v=c(0,1), col=”black”, lwd=1, lty=2)
#ESPLORAZIONE GRAFICA DI DATI MULTIVARIATI
#Per avere un’idea attraverso il PLOT TRIDIMENSIONALE: x11 () colors = as.numeric (dati$variabile_categorica) plot3d (dati$X, dati$Y, dati$Z, col = colors) #solo guardando i grafici, è già possibile farsi un’idea della struttura di correlazione di un dataset plot (dati) #la funzione ggpairs di Ggally è più flessibile. Ggpairs (dati, aes (alpha = 0.4, color = sesso)) # Nella parte superiore della tabella a incrocio, è riportato il risultato di un test di correlazione di Pearson fatto fra le variabili relative ad ogni incrocio. Il test è fatto sia per la popolazione completa, che per le sottopopolazioni relative a ciascun genere. Il valore numerico rappresenta la stima della correlazione, gli asterischi riportano invece il livello di significatività del test di correlazione [*** vuol dire MOLTO significativo]
METODO 2: QQPLOT
#consiste nel sovrapporre i quantili empirici ricavati dalle osservazioni con la retta dei quantili teorici della normale qqnorm (X) qqline (X, col = “red”) #in caso di normalità perfetta, i dati giacciono sulla retta METODO 3: SHAPIRO-WILK TEST #H0: Dati NORMALI shapiro.test (X) #MI RESTITUISCE IL P-VALUE, DA CONFRONTARE CON UNA SOGLIA DI SIGNIFICATIVITÀ SELEZIONATA A PIACERE E SOLITAMENTE ASSUNTA PARI ALL’1 O AL 5% Per p-value bassi (minori del livello di significatività), rifiutiamo l’ipotesi nulla: dobbiamo concludere che i dati non siano distribuiti normalmente. Per p-value alti, H0 non può essere rifiutata (devo accettare): dobbiamo concludere che i dati provengano da una popolazione distribuita normalmente. #SE NON C’E’ EVIDENZA STATISTICA PER RIFIUTARE H0 (P-VALUE ALTI / QUANTILI EMPIRICI DISTRIBUITI LUNGO LA RETTA / ISTOGRAMMA CONFRONTABILE CON LA “CAMPANA” DELLA DISTRIBUZIONE NORMALE), SI PUO’ LAVORARE SOTTO L’IPOTESI DI GAUSSIANITA’ DEI DATI.
#TEST DI IPOTESI
#Esistono 3 casi:
- La popolazione segue una distribuzione normale : si procede con un t-test _, utilizzando i quantili della distribuzione t di Student Se i dati NON sono normali, posso provare a fare lo shapiro.test su dati_new pari a log(dati$Colonna). Se è verificata l’ipotesi di normalità, posso fare i t-test su dati_new.
- La popolazione è_ non normale , ma le osservazioni sono numerose : si procede con uno z-test _approssimato basato sul TCL, utilizzando i quantili della distribuzione normale
- La_ popolazione è non normale, e la numerosità scarsa : non si può trarre alcuna conclusione senza ulteriori assunzioni #L’errore di primo tipo è la probabilità di rifiutare H0 quando è vera. #L’errore di secondo tipo è la probabilità di accettare H0 quando è falsa. CASO 1: T-TEST SULLA MEDIA DI UNA SINGOLA POPOLAZIONE X: #La popolazione segue una distr. NORMALE mu_0 = MEDIA_DA_CONTROLLARE n = length (dati$X) media_campionaria = mean(dati$X) dev_std_campionaria = sd (dati$X) t = (media_campionaria – mu_0)*sqrt(n)/dev_std_campionaria #Statistica test UNILATERO A SX
H0: mu >= mu0 VS H1: mu < mu
pval = pt (t, df = n-1) #Calcolo il pvalue oppure t.test (dati$X, mu = MEDIA_DA_CONTROLLARE, alternative = “less”) UNILATERO A DX
H0: mu <= mu0 VS H1: mu > mu
pval = 1 – pt (t, df = n-1) oppure
t.test (dati$X, mu = MEDIA_DA_CONTROLLARE, alternative = “greater”) BILATERO:
H0: mu = mu0 VS H1: mu != mu
pval = 2*(1-pt (abs(t), df = n-1)) oppure t.test (dati$X, mu = MEDIA_DA_CONTROLLARE, alternative = “two.sided”) o t.test(X, Y, var.equal = T)
_- p-value ALTO: l’evidenza empirica non è sufficientemente contraria all’ipotesi nulla H0, che quindi non può essere rifiutata
- p-value BASSO: l’evidenza empirica è fortemente contraria all’ipotesi H0 che quindi va rifiutata_ INTERVALLI DI CONFIDENZA alpha = 0. quantile = qt (1-alpha/2, df = n-1) c (media_campionaria – quantile/sqrt(n)dev_std_campionaria, media_campionaria + quantile/sqrt(n)dev_std_campionaria) #Se l’intervallo di confidenza non contiene mu_0, rifiuto H0; se contiene mu_0, non ho evidenza statistica per rifiutare H0: devo accettare l’ipotesi nulla. CASO 2: Z-TEST SULLA MEDIA DI UNA SINGOLA POPOLAZIONE X: #La popolazione NON normale, MA NUMEROSA SI USA LA NORMALE STD E NON LA T DI STUDENT mu_0 = MEDIA_DA_CONTROLLARE n = length(dati$X) media_campionaria = mean (dati$X) dev_std_campionaria = sd (dati$X) z = (media_campionaria – mu_0)*sqrt(n)/dev_std_campionaria #Statistica test UNILATERO A SX
H0: mu >= mu0 VS H1: mu < mu
pval = pnorm (z) UNILATERO A DX
H0: mu <= mu0 VS H1: mu > mu
pval = 1 – pnorm (z) BILATERO
H0: mu = mu0 VS H1: mu != mu
pval = 2*(1-pnorm (abs (z))) _- p-value ALTO (maggiore del livello di significatività alpha): l’evidenza empirica non è sufficientemente contraria all’ipotesi nulla H0, che quindi non può essere rifiutata, allora NON RIFIUTO H0.
- p-value BASSO (minore del livello di significatività alpha): l’evidenza empirica è fortemente contraria all’ipotesi H0 che quindi va rifiutata, allora RIFIUTO H0._ INTERVALLI DI CONFIDENZA alpha = 0. quantile = qnorm (1-alpha/2) n = length(data$X) Se l’intervallo di confidenza non contiene mu_0, rifiuto H0; se contiene mu_0, non ho evidenza statistica per rifiutare H0: devo accettare l’ipotesi nulla. #Se VOGLIO STABILIRE SE I DATI CI PERMETTONO DI AFFERMARE CON UNA CONFIDENZA DEL 95% CHE IL #TRATTAMENTO HA SUCCESSO ALMENO NEL 60% DEI CASI: p.0 = 0. alpha = 0. z.alpha = qnorm(1-alpha) p.camp = length(which(Rid_Temp==’S’))/n Z.0 = (p.camp – p.0)/ sqrt(p.0*(1 – p.0)/n) Z.0 > z.alpha #Controllo!
H0: mu >= mu0 VS H1: mu_X < mu_Y
pvalue = pt (t, df = n_X + n_Y – 2) oppure t.test (X, Y, var.equal = T, alternative = “less”)
- VARIANZE DIVERSE O NON OMOGENEE t.test (X, Y, mu = delta0, var.equal = F) t.test (X, Y, alternative = “greater”, mu = delta0, var.equal = F) t.test (X, Y, alternative = “less”, mu = delta0, var.equal = F) #In esercizi come quello sotto devo considerare distribuzioni gaussiane e varianze non note, omogenee
- VARIANZE NOTE: #le varianze possono essere note nel momento in cui il dato ci viene fornito dall’esterno, e sono VAR_X e VAR_Y. Differenza_medie = mean (X) – mean (Y) n_X = length (X) n_Y = length (Y) var_X = VAR_X var_Y = VAR_Y z = differenza_medie / sqrt (var_X/n_X + var_Y/n_Y) pvalue = 2*(1-pnorm(abs(z))) #bilatero, per gli altri guardare codici precedenti Se il p-value è ALTO (maggiore del livello di significatività alpha), l’evidenza empirica non è sufficientemente contraria all’ipotesi nulla, che non posso rifiutare: non rifiuto H0. Se il p-value è BASSO, l’evidenza empirica è fortemente contraria all’ipotesi nulla: rifiuto H0. TEST ANOVA: #Analysys Of Variance: n osservazioni divise in g gruppi con numerosità nj. H0: I gruppi hanno TUTTI stessa media VS H1: ALMENO UN gruppo ha media diversa dagli altri Effettuare uno shapiro test per verificare la normalità dei dati PER OGNI SOTTOCATEGORIA individuata tapply (dati$X, dati$VAR_CATEGORICA, shapiro.test) #Restituisce un p-value per ogni livello della variabile categorica rispetto a cui ho diviso la variabile X: per p-value ALTI si accetta l’ipotesi di normalità #La seconda Ipotesi è l’OMOSCHEDASTICITÀ, ovvero l’omogenità delle varianze. H0: Le varianze sono omogenee Se accetto H0 (pval alto) lavoriamo sotto ipotesi di omoschedasticità bartlett.test (dati$X ~ dati$VAR_CATEGORICA) #Se ho pochi dati devo anche verificare che le classi siano tra loro bilanciate : table(dati$VAR_CATEGORICA) #A questo punto è possibile iniziare col test. #Prima di tutto facciamoci un’idea se le medie sono uguali o meno graficamente: boxplot (dati$X ~ dati$VAR_CATEGORICA) #Applichiamo il test vero e proprio: model = aov (X ~ VAR_CATEGORICA, data = dati) summary(model)
#Df sono i gradi di libertà. La prima riga corrisponde a g – 1 in cui g è il numero di gruppi della variabile categorica, mentre la seconda corrisponde a n – g, in cui n è il numero totale di osservazioni. L’ultima colonna è il pvalue del test ANOVA: H0: I gruppi hanno TUTTI stessa media VS H1: ALMENO UN gruppo ha media diversa dagli altri La statistica F è il rapporto tra varianza between e within (max varianza between e minima within). La naturale prosecuzione dell’analisi consiste nel valutare quali tra le possibili coppie di gruppi presentano media significativamente diversa. Occorre dunque produrre intervalli di confidenza per le differenze delle medie dei gruppi. H0 della statistica F: tutti i coefficienti sono uguali a zero (voglio rifiutare, p-value basso) La stima della varianza comune si ricava facilmente dalla somma delle varianze nei gruppi: S = sum (model$residuals ^ 2) / model$df.residual PER FARE ESEGUIRE AL MODELLO OPPORTUNI INTERVALLI DI CONFIDENZA: source(“Confronti_multipli_da_anova.R”) alpha = 0. intervalli_anova = Confronti_multipli_da_anova(model, alpha = alpha, bonferroni = TRUE) intervalli = intervalli_anova$intervals grafico_intervalli = intervalli_anova$plot plot(grafico_intervalli) #Tale funzione fornisce gli intervalli di confidenza cercati, i p-value relativi al t-test di uguaglianza tra le medie, ed un oggetto grafico per visualizzare i risultati. Da leggere tipo boxplot: per ogni gruppo (asse x) vedo l’intervallo di confidenza (asse y). Prima vera Seconda falsa: p-value alto, i tre gruppi non sono significativamente diversi in termini di media Terza falsa: il p-value è la probabilità di osservare qualcosa di più raro di quello che ho osservato (non sono tanto raro se ho il 30% di trovare qualcosa di più raro di me) Quarta falsa: ha distribuzione F di fisher con 2 parametri Quinta falsa: non c’entra nulla la gaussianità dei residui: c’è l’assunzione di normalità, indipendenza e uguale varianza, ma il p-value non c’entra niente sul p-value della normalità dei dati.
#PRINCIPAL COMPONENT ANALYSIS
#La PCA è un metodo di riduzione dimensionale. È utile per evitare ridondanza di informazione tra le varie variabili, quindi è bene essere in grado di comprendere se può essere utile o meno effettuare una PCA prima di procedere.
#Possiamo anche visualizzare le componenti stesse che corrispondono agli autovettori della #matrice di covarianza campionaria, mentre gli autovalori associati corrispondono alle #varianze delle componenti stesse: eigen(data_scaled)$vector #PER INTERPRETARE IL SIGNIFICATO DEI LOADINGS (autovettori della matrice di covarianza campionaria): loadings = PC$loadings x11 () par (mar = c(1, 4, 0, 2), mfrow = c(n, 1)) #n = numero di componenti principali individuate for (i in 1:n) barplot (loadings [, i], ylim = c (-1, 1)) #Questo grafico ci elenca i barplot dei loadings per ciascuna componente principale: abbiamo dei #rettangolini la cui altezza coincide con il loading associato alla componente principale. Ad esempio, #lungo la prima riga le altezze sono phi11, phi12, phi13, … #In altre parole, per ogni componente plottiamo il barplot dei loadings. #Aggiunta una soglia sull’intensità dei loading da plottare, vengono plottati solo se superano un valore di soglia che stabiliamo a priori: soglia = valore_soglia x11 () par (mar = c (1, 4, 0, 2), mfrow = c (n, 1)) for (i in 1:n)
barplot(ifelse (abs (PC$loadings[,i]) < soglia, 0, PC$loadings[ ,i]), ylim = c (-1, 1)) abline (h=0) #Con lo Scree Plot visualizzo i loadings e la percentuale di varianza da loro spiegata (vedo il contributo di ogni componente alla variabilità totale): x11() plot (cumsum(PC$sd^2)/sum(PC$sd^2), type=’b’, axes=F, xlab=’numero di componenti’, ylab=’contributo alla varianza totale’, ylim=c(0,1)) abline (h=1, col=’blue’) abline (h=livello_soglia, lty=2, col=’blue’) box () axis (2, at=0:10/10, labels=0:10/10) axis (1, at=1:ncol (data_scaled), labels=1:ncol(data_scaled),las=2) dev.off() #Altri strumenti grafici per interpretare la PCA: # La media campionaria può essere vista come la “componente 0”, ovvero il punto dello #spazio che minimizza la somma dei quadrati delle distanze dei punti del campione M = colMeans(X) plot3d(X, size=4) aspect3d(“iso”) axes3d() for(i in 1:nobs) { lines3d(rbind(X[i,], M), lwd = 0.5) } points3d(t(M), col=’red’, size=12) #Per visualizzare graficamente la prima componente principale #La prima componente principale rappresenta invece la direzione della retta, passante per la media campionaria, che minimizza la somma dei quadrati delle distanze (o, equivalentemente, massimizza la varianza delle proiezioni dei punti del campione su tale retta) plot3d(X, size=4) aspect3d(“iso”) axes3d() PC1 = NULL for(I in 1:nobs){ PC1 = rbind(PC1, PC$loadings[,1]PC$scores[I,1] + M) } points3d(PC1, col=’red’, size=6) for(I in 1:nobs){ lines3d(rbind(X[I,], PC1[I,]),col=’blue’) } lines3d(rbind(M + 2PC$sdev[1] * PC$loadings[,1], M – 2*PC$sdev[1] * PC$loadings[,1]), col=’forestgreen’,lwd=0.5) #La seconda componente, ortogonale alla prima, definisce assieme ad essa un piano, passante per la media campionaria, con le stesse proprietà di approssimazione: plot3d(X, asp=1, size=4)
#Anche l’ANALISI DEGLI SCORES (con il supporto delle variabili categoriche eventualmente presenti nel dataset) è fondamentale scores = data.frame (PC$scores) names (scores) = paste0 (“Comp”, 1:n) scores = cbind (scores, data_labels) str (scores) boxplot (scores[,1:n], col = “gold”) #Poi, a seconda dei nostri scopi, possiamo stratificare l’analisi inserendo le variabili categoriche: boxplot (scores$CompN ~ scores$VAR_CATEGORICA_DI_INTERESSE) #Per capire quante e quali componenti tenere, esistono due metodologie: #METODO DEL GOMITO (con lo scree plot) x11() plot (cumsum(PC$sdev^2)/sum(PC$sdev^2), type=’b’, axes=F, xlab=’number of components’, ylab=’contribution to the total variance’, ylim=c(0,1)) abline (h=1, col=’blue’) abline (h=0.95, lty=2, col=’blue’) box () axis(2,at=0:10/10,labels=0:10/10) axis(1,at=1:ncol(dati),labels=1:ncol(dati),las=2) #Possiamo fare un’ipotesi per quanto riguarda la presenza di un gomito (un punto dopo il quale possiamo osservare un appiattimento della curva), da confermare con il controllo dei boxplot degli scores: scores = data.frame (PC$scores) x11 () boxplot (scores, col = “gold”) #In questo caso, la presenza di outliers (in ciascuno dei boxplot) non maschera una eventuale alta variabilità delle componenti successive alla seconda. Confermiamo la scelta di tenere le prime due componenti. #METODO DELLA SOGLIA x11() plot (cumsum(PC$sdev^2)/sum(PC$sdev^2), type=’b’, axes=F, xlab=’number of components’, ylab=’contribution to the total variance’, ylim=c(0,1)) abline (h=1, col=’blue’) abline (h=valore_soglia, lty=2, col=’blue’) box () axis(2,at=0:10/10,labels=0:10/10) axis(1,at=1:ncol(dati),labels=1:ncol(dati),las=2)
#CLUSTERING
#K-MEANS: ALGORITMO DI CLUSTERING
#Utilizziamo solo la metrica euclidea o euclidea standardizzata perché con metriche diverse, l’algoritmo non garantisce convergenza. L’algoritmo è finalizzato a minimizzare la variabilità NEI clusters e massimizzare la variabilità TRA I clusters #INNANZITUTTO, visualizziamo graficamente i dati: dati = read.csv(“data.txt”) x11 () plot (dati, xlab=’X.1’, ylab=’X.2’, pch = 20) #pch = 20 è il comando grafico che specifica che voglio, per indicare i punti, dei pallini, #mentre pch = 5 è il comando grafico che specifica la presenza di rombi #Se, invece, devo, ad esempio, rappresentare un dataset con tre variabili: x11 () plot (dati) #In 3D: x11 () plot3d (dati, col = “orange”, aspect = F) In R, abbiamo una funzione: kmeans, che è caratterizzata da: INPUT: dati: dati da raggruppare (primo parametro da passare) centri: numero di cluster o vettore di centroidi iniziali iter.max: numero Massimo di iterazioni da eseguire distanza: distanza da utilizzare nstart: specifica quante volte ripetiamo kmeans per diverse inizializzazioni (aiuta a trovare un minimo globale)
Qui non emerge una netta
separazione tra i gruppi.
#se il dataset ha tre variabili, la rappresentazione grafica viene così: km.output = kmeans(dati, centers = K, nstart = 20) x11 () plot (dati, col = km.out$cluster) #IN 3D plot3d (dati, size=3, col=km.output$cluster, aspect = F) points3d (km.output$centers, pch = 4, cex = 2, lwd = 4) #Per rappresentare la Total Within Sum of Squares e la Between Sum of Squares: w = b = NULL for (k in 1:10) { result.k = kmeans(dati, K, nstart = 20) w = c(w, sum(result.k$withinss)) #sum(result.k$withinss) equivale a result.k$tot.withinss b = c(b, result.k$betweenss) } x11 () par (mfrow = c(1,2)) plot (1:10, b,type=’b’, xlab=’clusters’, ylab=’betweenss’, main=’Andamento della BetweenSS’,col=’orange’) plot (1:10, w,type=’b’, xlab=’clusters’, ylab=’tot.withinss’, main=’Andamento della total WithinSS’,col=’purple’) #Per scegliere il numero totale di clusters esistono due criteri: #1) CRITERIO TOTAL WITHIN SUM OF SQUARES: Plotto il grafico della total within sum of squares (tot.withinss) con un numero diverso di clusters per confrontare le curve, cercando di individuare il valore di K che ci dia il miglior trade off tra la minimizzazione della withinss, della betweenss e evitando l’overfitting: tot.withinss = NULL #Definisco un vettore nullo che riempio a ciascuna iterazione
Si noti come i grafici siano del tutto equivalenti, poiché, come già
sottolineato, tot.withinss e betweenss hanno somma costante.
for (I in 1:10) { tot.withinss = c(tot.withinss,kmeans(dati, I, nstart=20)$tot.withinss)} x11() plot (1:10, tot.withinss, type=’b’, xlab=’K’, main=”Andamento della total within sum of squares”) Come scelgo k? REGOLA DEL GOMITO : scelgo k in modo tale da osservare una forte diminuzione prima e poi, invece, una sorta di “appiattimento” della curva. ù #Nel caso di CLUSTER SOVRAPPOSTI, solitamente, non troviamo un “gomito”, e, quindi, si preferisce guardare al grafico della tot.withinss normalizzandolo rispetto alla total sum of squares (ovvero, la variabilità totale). Si ottiene una curva del tutto equivalente ai fini interpretativi, che ha però il pregio di assumere valori tra 0 e 1: la si può interpretare come frazione della variabilità totale dovuta alla variabilità intra-cluster. Ratio = NULL for (I in 1:10) { kmeans_out = kmeans (dati, I, nstart= 20 ) ratio = c (ratio,kmeans_out$tot.withinss/kmeans_out$totss)} x11() plot (1:10,ratio,type=’b’, xlab=’K’, ylab = “tot.withinss/totss”, main=”Andamento della tot.withinss normalizzata”) #2) CRITERIO DELLA SILHOUETTE ( CLUSTER DISTINTI E SOVRAPPOSTI) #Il coefficiente di silhouette misura quanto un oggetto sia simile al suo stesso cluster rispetto agli altri clusters. #Se l’Indice di Silhouette è circa 1 significa che i cluster sono molto densi e ben separati; #Se l’Indice di Silhouette è circa 0 significa che i cluster si sovrappongono; #Se l’Indice di Silhouette è minore di 0 significa che i cluster potrebbero essere sbagliati/non corretti. Scelgo il k che mi dà un indice di silhouette (medio e locale) maggiore! X11() km.output = kmeans (dati, centers = K, nstart = 20 ) silhouette_cluster = silhouette (km.out$cluster , dist(dati) ) plot(silhouette_cluster, col = 1:K, main = ‘Silhouette plot per k=K’) #in input il risultato del calcolo dell’indice di silhouette abline (v = mean(silhouette_cluster[,K])) #Per fare il confronto con k = 2,3,4,5 cluster: x11 () par(mfrow = c(2,2)) for ( i in 2:5) { output.sim = kmeans(dati, centers = i, nstart = 20 ) silhouette.sim = silhouette(output.sim$cluster , dist(dati) )
Vediamo un chiaro guadagno nel passaggio da 1
a 2 cluster, mentre l’aggiunta di altri cluster oltre
al secondo produce guadagni sempre minori,
fino a diventare del tutto marginali.
