mercoledì 30 giugno 2010

Cluster Analysis in R #1: Hierarchical Clustering

Dietro suggerimento e con il prezioso aiuto di Mirko Modenese, statistico presso l’European Academy of Research (www.eurac.edu) e docente di R presso l'Università Ca'Foscari di Venezia (Facoltà di Economia), che qui ringrazio per la sua disponibilità, ho deciso di sviluppare in due posts il tema della Cluster Analysis in R.
Questo è l'indice degli argomenti svolti:


Introduzione alla Cluster Analysis \
Hierarchical clustering |
Esempio 1. |
Esempio 1.R | CA in R #1
Esempio 2.R |
Varianti |
Il taglio del dendrogramma /

Partitional clustering \
Definizione di distanza |
Distanza di Minkowski |
Distanza euclidea |
Algoritmo K-means |
Esempio 1 | CA in R #2
Esempio 2 |
Esempio 1.R |
Esempio 2.R |
Varianti /





Introduzione alla Cluster Analysis

La ricerca e la semplificazione di popolazioni di oggetti (unità statistiche) in base alle loro caratteristiche (variabili) è da sempre oggetto di studio dell’uomo: dalla classificazione delle specie animali, alla catalogazione della flora, del patrimonio bibliotecario, dei gruppi etnici, del patrimonio genetico, etc. Il processo di identificazione di gruppi (cluster) può avvenire sostanzialmente secondo due differenti approcci: supervisionato e non supervisionato. Con il primo, andiamo a catalogare gli oggetti in base a dei criteri prestabiliti; con il secondo la classificazione è meramente esplorativa e si basa sulle caratteristiche rilevate per singola unità statistica (questo è oggetto dei presenti articoli).
Lo scopo della cluster analysis è quello di raggruppare le unità sperimentali in classi secondo criteri di (dis)similarità (similarità o dissimilarità sono concetti complementari, entrambi applicabili nell’approccio alla cluster analysis), cioè determinare un certo numero di classi in modo tale che le osservazioni siano il più possibile omogenee all'interno delle classi (coesione interna) ed il più possibile disomogenee tra le diverse classi (divisione esterna).
Gli approcci classici all’unsupervisioned clustering sono di due tipi: il clustering partitivo e il clustering gerarchico. Con la prima tecnica, si definisce l’appartenenza ad un gruppo utilizzando la dissimilarità da un punto rappresentativo del cluster (centroide); questo è l’argomento del prossimo post.
Con il clustering gerarchico viene invece creata una gerarchia di partizioni caratterizzata da un numero crescente o decrescente di gruppi, visualizzabile mediante un diagramma ad albero (dendrogramma).
I principali algoritmi del hierarchical clustering sono due:
  1. possiamo seguire un approccio bottom-up: si creano tanti clusters quanti sono gli oggetti (singleton), e per passaggi successivi si uniscono gli oggetti (leggasi clusters) tra loro più vicini, fino ad ottenere un gruppo contenente tutti gli oggetti (hierarchical agglomerative clustering)

  2. oppure utilizziamo un approccio top-down: da un unico cluster contenente tutti gli oggetti, si splitta l'insieme fino ad ottenere clusters contenenti ciascuno un solo oggetto (hierarchical divisive clustering).

Il metodo più usato è quello bottom-up, ed è quello che qui verrà spiegato.

Iniziamo subito con gli esempi numerici.




Hierarchical clustering

Esempio 1.
Supponiamo di avere le distanze tra 7 città, organizzate in una matrice triangolare come la seguente (matrice delle distanze):



Nota matematica: se abbiamo m elementi, calcolando la matrice delle distanze (che è una matrice simmetrica triangolare con diagonale principale –banalmente- nulla), otteniamo distanze (nell’esempio 7*6/2= 21 distanze).

Individuiamo la distanza minima tra tutte le distanze presenti (in questo caso è quella tra Torino e Milano):



A questo punto creiamo un cluster contenente le due città tra loro più vicine -quindi avremo il cluster (To,Mi)-, e lo sistemiamo nella matrice delle distanze come segue:



Come legare i cluster che si formano man mano per fusione (con numerosità di gruppo pari a p) con i precedenti gruppi (con numerosità di gruppo p-1)? Si utilizza la formula di Lance-Williams, così definita:



che calcola la dissimilarità tra i gruppi ed . Alfa, beta e gamma sono i coefficienti che identificano la scelta del legame (singolo, completo, centroide, medio, Ward, etc...), ognuno dei quali conferisce delle caratteristiche differenti ai gruppi così formati (maggior coesione interna, maggior divisione tra gruppi, etc…) – v. Tab sottostante.



Decidiamo di utilizzare un single linkage: scegliamo come distanza tra due clusters il minimo tra le distanze calcolate tra tutti gli elementi dei clusters:



Possiamo così completare lo schema:



Ripetiamo la procedura, andando ad individuare la distanza minore:



Il secondo cluster è Napoli-Roma, e come sopra uniamo le righe e le colonne relative, e calcoliamo le distanze:





Cerchiamo la distanza minore:



Uniamo i valori, e calcoliamo le distanze:





Cerchiamo la distanza minore:



Uniamo i valori e calcoliamo le distanze:





Cerchiamo la distanza minore:



Uniamo i valori e calcoliamo le distanze:



Il raggruppamento in un cluster unico è:
(((Rm,Na),Ba),Pa),((To,Mi),Ve)

Da qui possiamo disegnare il dendrogramma.
Il dendrogramma non è altro che il diagramma costruito sulle matrici precedentemente calcolate sulla base della formula di L&W, che definiscono il legame tra cluster (legame singolo, legame completo, legame medio, etc...). Sull’asse delle ordinate si ha l’altezza di fusione tra due elementi (dal basso verso l’alto) e sull'asse delle ascisse (sempre dal basso verso l’alto) si denota graficamente dove vengono raggruppati i singleton. Per capire quale sia la relazione tra due elementi, tracciate un percorso da uno dei due all'altro, seguendo i rami del diagramma ad albero e scegliendo la strada più breve. La distanza dall'origine alla linea verticale più esterna toccata dal percorso rappresenta il grado di somiglianza tra i due elementi.




Esempio 1.R

Per disegnare il dendrogramma in R, possiamo usare la funzione hclust, del pacchetto stats, la quale riceve in input la matrice delle distanze (ovvero l’output della funzione dist).
Per creare la matrice delle distanze, poichè come dati abbiamo direttamente le distanze, procediamo come segue:


mydata <- c(0,363,486,418,534,177,397,
363,0,787,446,864,223,602,
486,787,0,884,125,658,247,
418,446,884,0,902,307,814,
534,864,125,902,0,711,367,
177,223,658,307,711,0,534,
397,602,247,814,367,534,0)

m1 <- matrix(mydata, 7,7)
rownames(m1) <- c("Roma", "Bari", "Milano", "Palermo", "Torino", "Napoli", "Venezia")
colnames(m1) <- c("Roma", "Bari", "Milano", "Palermo", "Torino", "Napoli", "Venezia")

m1
Roma Bari Milano Palermo Torino Napoli Venezia
Roma 0 363 486 418 534 177 397
Bari 363 0 787 446 864 223 602
Milano 486 787 0 884 125 658 247
Palermo 418 446 884 0 902 307 814
Torino 534 864 125 902 0 711 367
Napoli 177 223 658 307 711 0 534
Venezia 397 602 247 814 367 534 0

m1.dist <- as.dist(m1)

m1.dist
Roma Bari Milano Palermo Torino Napoli
Bari 363
Milano 486 787
Palermo 418 446 884
Torino 534 864 125 902
Napoli 177 223 658 307 711
Venezia 397 602 247 814 367 534


Adesso possiamo usare la funzione hclust, specificando il metodo single-linkage, e disegnare il dendrogramma:


hc <- hclust(m1.dist, "single")
plot(hc)




Il risultato è analogo a quello ottenuto by-hand. In ordinate vi è il livello di fusione dei cluster. Ad esempio abbiamo calcolato la distanza tra Milano e Torino, che è pari a 125: la riga che unisce queste due città cade proprio sull'ordinata 125, e in questo caso corrisponde alla distanza tra i due punti. Analogamente la distanza tra il cluster (Rm,Na) e il cluster (Ba) era 223, e la riga che unisce i due cluster cade proprio su 223 (per distanza qui intendo il livello a cui avviene la fusione tra clusters).
Possiamo richiamare i livelli di fusione così calcolati in questo modo (sono disposti nell'ordine delle righe orizzontali):


hc$height
[1] 125 177 223 247 307 397


Possiamo altresì visualizzare quali oggetti vengono inglobati e a quale livello (altezza) di fusione:


cbind(hc$merge,hc$height)
[,1] [,2] [,3]
[1,] -3 -5 125
[2,] -1 -6 177
[3,] -2 2 223
[4,] -7 1 247
[5,] -4 3 307
[6,] 4 5 397


Dove le prime due colonne rappresentano la fusione tra singleton (con il simbolo "-" ) e cluster (senza simbolo), mentre la terza colonna indica il livello di fusione. Appunto al primo raggruppamento si fondono in un unico gruppo Milano (-3) e Torino (-5). Al terzo step di raggruppamento vengono inglobati in un unico gruppo il cluster "2" (Rm-Na) e il singleton "-2" (Bari), e così via…


hc$labels
[1] "Roma" "Bari" "Milano" "Palermo" "Torino" "Napoli" "Venezia"




Una rappresentazione alternativa che si rileverà utile per il "taglio del dendrogramma" è la seguente:


par(mfrow=c(1,2))
plot(hc)
plot(hc$height)
for(i in 1:length(hc$height)) lines(c(i,i),c(0,hc$height[i]),col=4)


Sulla sinistra il dendrogramma, sulla destra le altezze di fusione dal primo cluster formato da due unità statistiche all’ultimo cluster che riassume tutta la popolazione oggetto di studio.






Esempio 2.R

Vediamo adesso un esempio bidimensionale.
Innanzitutto generiamo alcuni punti, casuali, con la funzione rnorm:


x <- abs(round(c(rnorm(10, mean=1, sd=0.5)), 1))*10
y <- abs(round(c(rnorm(10, mean=1, sd=0.5)), 1))*10

x
[1] 19 10 4 13 2 9 1 2 16 5
y
[1] 2 12 17 12 9 11 17 14 3 3


Adesso sistemiamo i dati in una matrice, con la funzione cbind (c sta per coloumn, mentre bind sta per unire: la funzione unisce due vettori in colonna):


mydata <- cbind(x,y)
rownames(mydata) <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "l")
mydata
x y
a 19 2
b 10 12
c 4 17
d 13 12
e 2 9
f 9 11
g 1 17
h 2 14
i 16 3
l 5 3


Possiamo plottare i 10 punti con la funzione plot:


plot(mydata, pch=15)
vec.label <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "l")
text(mydata+c(0.5,0.5), vec.label)




A questo punto generiamo la matrice triangolare delle distanze, con la funzione dist:


distance <- dist(mydata, "eucl")
distance
a b c d e f g h i
b 13.453624
c 21.213203 7.810250
d 11.661904 3.000000 10.295630
e 18.384776 8.544004 8.246211 11.401754
f 13.453624 1.414214 7.810250 4.123106 7.280110
g 23.430749 10.295630 3.000000 13.000000 8.062258 10.000000
h 20.808652 8.246211 3.605551 11.180340 5.000000 7.615773 3.162278
i 3.162278 10.816654 18.439089 9.486833 15.231546 10.630146 20.518285 17.804494
l 14.035669 10.295630 14.035669 12.041595 6.708204 8.944272 14.560220 11.401754 11.000000


Quest'ultimo parametro viene passato nella funzione hclust, che ci permette di generare il dendrogramma:


hc <- hclust(distance, "single")
plot(hc)


Abbiamo così ottenuto il dendrogramma relativo ai 10 dati.






Varianti

Negli esempi precedenti, per decidere la distanza tra due clusters (o un cluster e un punto), abbiamo scelto la minima distanza tra un elemento e . Questo metodo si chiama "single linkage" (o metodo del legame singolo), e può essere così descritto:

Sia C(xy) il cluster contenente gli elementi x e y, e sia z un terzo punto. La distanza tra C(xy) e z è:


Una più formale scrittura della precedente eq.ne è la seguente.
Considerando la già citata formula di Lance e Williams ed in riferimento alla relativa tabella dei coefficienti, sostituendo , , , , otteniamo l’equazione relativa al single linkage:



Sempre in riferimento alla tabella dei coefficienti, possiamo ricavare le eq.ni relative agli altri legami noti:
- completo,
- medio,
- centroide,
- ward.
La scelta del legame è guidata dai principi base del clustering: internal cohesion e external separation (v. Gordon,A.D. "Classification", 2nd edition – Chapman & Hall).




Il taglio del dendrogramma, un esempio pratico di clustering gerarchico

Consideriamo il dataset iris:

data(iris)
X <- iris


Tale dataset si riferisce a tre specie di fiori di iris su un campione di 150 unità statistiche (50 per specie). Non considerando la variabile specie, vogliamo testare il clustering agglomerativo in maniera non supervisionata e quindi meramente esplorativo. Effettuiamo un’analisi grafica dello scatterplot, cercando di capire se le variabili descrivono o meno dei gruppi ben definiti:

plot(X[,1:4])



Dalle relazioni tra variabili si intravedono due cluster "ben separati" e a tratti coesi internamente. Cerchiamo di capire se l’analisi dei gruppi può dirci qualcosa in più.

Effettuiamo l’analisi gerarchica mediante funzione hclust su matrice di dissimilarità euclidea con legame di ward:

hc <- hclust(dist(X[,1:4]),method="ward")

ora, l’oggetto hc contiene la nostra analisi dei gruppi gerarchica.

Quanti gruppi ho identificato?
La scelta non può essere arbitraria e deve rispettare il principio con il quale si effettua il clustering, cioè la semplificazione del campione in un numero di gruppi non troppo basso (un gruppo unico equivale al campione e non ha senso) o elevato - la numerosità campionaria n è altrettanto assurda (considererei i singleton come singoli gruppi). Un numero g di gruppi ideale (dipendentemente dal caso in esame) può variare tra .
Procediamo con il taglio del dendrogramma. Tale operazione si basa su un concetto molto semplice: confrontando le altezze di fusione del dendrogramma, si valutano come possibili livelli di taglio - al di sotto dei quali si identificano i componenti dei cluster (e quindi la scelta del numero di gruppi) - i "salti" maggiori tra una fusione e l’altra. Un’alta differenza tra un livello di merging e l’altro denota una divisione netta tra i gruppi (viceversa, se bassa, si rileva un’alta coerenza interna tra gli elementi che compongono il gruppo). Il taglio è visibile dal seguente grafico:




par(mfrow=c(1,2))
plot(hc)
abline(h=hc$height[147],col=2,lwd=2,lty=3)
plot(hc$height,ylab="height",xlab="merging steps",xlim=c(120,150))
for(i in 1:length(hc$height)) lines(c(i,i),c(0,hc$height[i]),col=4)
abline(h=hc$height[147],col=2,lwd=2,lty=3)


Il "salto" più consistente avviene tra il terzultimo ed il penultimo step agglomerativo (cioè tra la divisione in 3 e 2 cluster); ciò significa che la numerosità di gruppi che hanno maggior coesione interna e divisione esterna (tra gruppi) è g = 3. Tali raggruppamenti sono visibili anche dal dendrogramma, al quale abbiamo aggiunto un’ideale linea rossa: il taglio, appunto.

E ora? Abbiamo tre gruppi, andiamo a vedere come si visualizzano nello scatterplot:

plot(X[,c(1:4)], col=cutree(hc, 3),pch=cutree(hc, 3))

il comando cutree è proprio il taglio del dendrogramma per g = 3 gruppi sull’oggetto hc. Dando una sguardo al dataset originale verifichiamo che i cerchi neri sono la specie "setosa", i triangoli rossi la specie "versicolor" e le croci verdi la specie "virginica".
L’algoritmo hclust ha funzionato perfettamente, o quasi? Ci sono degli outlier? A voi la soluzione.