--- title: "Algoritmi non gerarchici: K-means e PAM" output: pdf_document: default html_document: df_print: paged --- Esistono molti algoritmi di clustering diversi. Questo laboratorio si concentra sul metodo delle $K$ medie o $K$-means. ## Algoritmo K-Means La funzione ``kmeans()`` implementa il metodo delle $K$ medie in R, prendendo come input i dati e il numero di clusters desiderato. L'esempio che segue utilizza i dati in ``food_spending.csv``, riguardanti il numero medio (mensile) di pasti consumati fuori casa (ad es., al ristorante) e la spesa per i pasti consumati a casa (annuale) per un insieme di 100 famiglie. ```{r} dat<-read.csv("food_spending.csv", header=TRUE) head(dat) ``` Per eseguire l'algoritmo $K$-means standardizziamo prima i dati e poniamo $K = 2$ nella funzione ``kmeans()`` (la funzione ``set.seed()`` viene utilizzata per garantire che il risultato sia riproducibile). ```{r} x.scaled<-scale(dat) set.seed(4) km.out <- kmeans(x.scaled, 2) km.out # see the output ``` L'output stampato mostra: - i vettori delle medie di gruppo (Cluster means) - il vettore del clustering, cioè un vettore di numeri interi (da $1$ a $K$) che indica il cluster di appartenenza di ciascuna unità - la devianza *within* per ciascun gruppo $$ WSS_k=\sum_{i\in C_{k}} \sum_{j=1}^{p} (x_{ij}-\bar{x}_{kj})^2 $$ ($\bar{x}_{kj}$ è la media della j-esima variabile nel gruppo $k$) e il rapporto "devianza TRA i gruppi"/"devianza TOT" ($BSS/DEV_{TOT}$), che rappresenta la percentuale di varianza spiegata tramite il clustering $$ \frac{BSS}{BSS+WSS}=\frac{BSS}{DEV_{TOT}}=\frac{\sum_{k=1}^{K} \sum_{j=1}^{p} |C_k|(\bar{x}_{kj}-\bar{x}_{j})^2}{\sum_{j=1}^{p}\sum_{i=1}^{n} (x_{ij}-\bar{x}_{j})^2} $$ dove $WSS=\sum_{k=1}^K WSS_k$. La partizione finale è contenuta nel vettore ``km.out$cluster``: ```{r} str(km.out$cluster) ``` I due cluster ottenuti contengono rispettivamente 56 e 44 osservazioni, le unità sono raggruppate in due cluster ben separati, come ci si aspettava. I due gruppi possono essere riportati con colori diversi nel grafico di dispersione: ```{r, fig.width=4.5, fig.height=4} plot(dat, col =(km.out$cluster+1), main="K-means (K=2)", xlab ="", ylab="",pch=19) ``` In questo esempio con dati simulati sapevamo che in realtà c'erano due cluster ben distinti. Tuttavia, supponendo di non conoscere il numero di gruppi, avremmo potuto eseguire l'algoritmo con $K = 3$ ```{r} set.seed(5) km.out<-kmeans(x.scaled, 3) km.out ``` ```{r, fig.width=4.5, fig.height=4} plot(dat, col =(km.out$cluster+1), main="K-means (K=3)", xlab ="", ylab="", pch=19) ``` La funzione ``kmeans()`` consente di eseguire l'assegnazione iniziale un certo numero di volte e scegliere la configurazione iniziale migliore attraverso l'argomento ``nstart'' che può essere incrementato arbitrariamente: ```{r} set.seed(6) km.out<-kmeans(dat, 3, nstart=1) km.out$tot.withinss km.out<-kmeans(dat, 3, nstart=20) km.out$tot.withinss ``` **$\color{red}{\text{Esercizio}}$** *Si considerino i dati nel file "sanita.txt" che si riferiscono ad un rilevamento ISTAT sui servizi sanitari nelle regioni italiane (i dati sono del 1998). I dati sono utilizzati in Rapallo & Rogantin 'Statistica descrittiva multivariata' (2003). Le variabili su cui effettuare l'analisi sono: il numero di medici per 10.000 abitanti; il numero di posti letto ospedalieri per 10.000 abitanti; il numero di pediatri per 10.000 abitanti di età inferiore ai 14 anni; il numero annuo di interventi di pronto soccorso per 1.000 abitanti; la percentuale di persone che si sono dichiarate soddisfatte dell'assistenza medica ospedaliera. Si importino i dati con il seguente comando* ```sanita<-read.table("sanita.txt", header =TRUE, row.names=1)``` *Si effettui una analisi di raggruppamento utilizzando il metodo delle $K$-medie.* ```{r, echo=FALSE, eval=FALSE} sanita<-read.table("sanita.txt", header =TRUE, row.names=1) str(sanita) sanita_st<-scale(sanita) #### standardizzazione delle variabili d <- dist(sanita_st, method = "euclidean") num_clus<-3 agg_kmeans<-kmeans(sanita_st,center=num_clus) agg_kmeans[1:7] #Una successiva esecuzione dell'algoritmo quindi può condurre a una diversa aggregazione. # Si può scegliere di fornire una lista di k centri iniziali (come valore di center nella chiamata della funzione kmeans). centri<-sanita_st[11:13,] agg_kmeans2<-kmeans(sanita_st,center=centri) agg_kmeans2$tot.withinss # o eseguire kmeans un numero prefissato di volte, con diverse scelte casuali di punti iniziali, tramite l'opzione nstart= nella chiamata della funzione agg_kmeans3<-kmeans(sanita_st,center=num_clus,nstart=20) agg_kmeans3$tot.withinss ``` ## Partitioning Around Medoids (PAM) L'obiettivo comune degli algoritmi $K$-means e PAM è di determinare una partizione ammissibile che sia ottimale rispetto alla omogeneità interna dei cluster ed alla eterogeneità tra i gruppi. Nell'algoritmo *Partitining Around Medoids* (PAM) si adoperano le distanze tra le unità. I $K$ centroidi diventano $K$ medoidi, che sono $K$ delle unità tra le $n$ considerate. La funzione ``pam()`` del pacchetto **cluster** implementa l'algoritmo PAM prendendo la matrice dei dati oppure la matrice di dissimilarità in input, e consentendo l'eventuale standardizzazione preliminare dei dati attraverso l'argomento ``stand = TRUE`` (il valore di default è FALSE). Se l'argomento è la matrice dei dati la distanza utilizzata di default è la distanza euclidea (la documentazione fornisce le alternative possibili): ```{r} library(cluster) pam.out<-pam(dat, 2, metric="euclidean", stand=TRUE) pam.out ``` I medoidi e il vettore del clustering si ottengono come segue: ```{r} pam.out$medoids #pam.out$clustering ``` ```{r, fig.width=4.5, fig.height=4} plot(dat, col =(pam.out$cluster+1), main="PAM (K=2)", xlab ="", ylab="", pch=19) points(pam.out$medoids, pch=as.character(pam.out$cluster[pam.out$id.med])) ``` **Nota** Per dataset di dimensioni elevate si veda la funzione ``clara()``, che risulta essere computazionalmente più efficiente. ### Scelta del numero di gruppi La scelta del numero di gruppi negli algoritmi di cluster analysis è un problema cruciale per questo tipo di analisi e non sempre risulta semplice individuare una soluzione. Negli algoritmi non gerarchici, il numero di gruppi è stabilito dal principio e pertanto occorre una strategia per la validazione della procedura di clustering adottata. Vengono descritti e implementati di seguito due degli indici *interni* più utilizzati per determinare il numero ottimale di gruppi: l'indice **silhouette** e l'**indice di Calinski e Harabasz**. Si veda il pacchetto di R **fpc** per molti altri indici comunemente utilizzati. L'indice *silhouette* (Rousseeuw, 1987) si basa sul grado di coesione dell’unità su cui viene calcolata con il guppo di appartenenza e sulla sua separazione con gli altri gruppi. Tale indice viene calcolato per ogni unità statistica attraverso la formula $$ s(i) = \frac{b(i)-a(i)}{\max\{ a(i), b(i) \}} $$ dove $a(i)$ è la distanza media dall'osservazione $i$ da tutte le unità appartenenti al medesimo gruppo; $b(i)$ è la più piccola tra le distanze medie dell'$i$-esima osservazione con le unità degli altri gruppi. Si noti che $s(i)\in [-1, 1]$: un $s(i)$ vicino a $1$ significa che l'unità $i$ è ben classificata, un $s(i)$ vicino a 0 implica che l'$i$-esima osservazione è in una situazione borderline tra due gruppi, infine un $s(i)$ vicino a $-1$ implica che l'$i$-esima osservazione è classificata nel gruppo sbagliato. La *silhouette media* (ASW) $S_K$ è data da $$ S_K =\frac{1}{n}\sum_{i=1}^{n}s(i) $$ dove $K$ è il numero di gruppi della partizione. Ottenendo la *silhouette media* al variare di $K$, la scelta di $K$ dovrà ricadere sul valore che massimizza $S_K$. L'output di ``pam`` consente di ottenere facilmente l'indice ASW: ```{r, fig.height=4, fig.width=4} plot(pam.out, which=2, main="") #silhouette plot ``` ```{r, fig.height=4, fig.width=4} plot(pam(dat, 3, metric="euclidean", stand=TRUE), which=2, main="") ``` L'indice di *Calinski e Harabasz* (Calinski and Harabasz, 1974) valuta la qualità della partizione ottenuta in termini di rapporto tra la varianza entro i gruppi e la varianza tra i gruppi. Per una data partizione è definito da $$ CH_K=\frac{BSS}{WSS}\frac{n-K}{K-1} $$ dove $n$ e $K$ sono rispettivamente il numero totale di osservazioni e il numero di cluster. Il numero di gruppi $K$ che risulta ottimale secondo tale indice `è quello che lo massimizza in quanto al numeratore compare un indicatore di separazione tra i gruppi e al denominatore un un indicatore di coesione. ```{r, fig.height=4, fig.width=4} require(fpc) require(ggplot2) Diss<-pam(dat, 2, stand = TRUE, keep.diss = TRUE)$diss minC <- 2 maxC <- 8 res <- sapply(minC:maxC, function(nc, data) cluster.stats(Diss, pam(data, nc, stand = TRUE, cluster.only=TRUE))$ch, data = dat) ggp <- ggplot(data=data.frame(x=2:(length(res)+1), y=res), mapping = aes(x=x,y=y)) + geom_point() + geom_line() + xlab("k") + ylab("Indice CH") print(ggp) ``` **$\color{red}{\text{Esercizio}}$** *La funzione ```pam``` consente di utilizzare una matrice di dissimilarità arbitraria impostando ``diss=TRUE``. Si implementi il metodo PAM per effettuare l'analisi di clustering dei dati in "claims.csv", che riguardano le caratteristiche di 1000 sinistri. Si ricordi che il dataset contiene tre variabili binarie e una variabile numerica. Si utilizzi una opportuna matrice di dissimilarità.* ```{r, echo=FALSE, eval=FALSE} d<-read.csv(file="claims.csv", header = TRUE) d$litig<-as.factor(d$litig) d$soft_injury<-as.factor(d$soft_injury) d$emergency_tr<-as.factor(d$emergency_tr) str(d) mixd<-daisy(d, metric="gower") # Gower's distance dclaims<-as.matrix(mixd) # transform to matrix pam.claims<-pam(dclaims, 2, diss=TRUE) pam.claims$silinfo$avg.width # add clustering dd<-cbind(d, clusterID=as.factor(pam.claims$clustering)) #groups<-split(dd, dd$clusterID) # data by cluster table(dd$clusterID) table(dd$litig, dd$clusterID) results <- data.frame(table(dd$litig, dd$clusterID)) results library(ggplot2) ggplot(data = results, aes(x = Var2, y = Freq, fill=Var1)) + geom_bar(stat = "identity")+ labs(x = "\n Cluster", y = "Frequency\n", title = "\n var: 'litig' \n", fill = "litig") ```