Analisi dei gruppi (clustering gerarchico)

L'analisi dei gruppi si applica a dati multivariati ed è un metodo statistico di tassonomia numerica che riveste un ruolo importante nella analisi esplorativa dei dati.

L'obiettivo dell'analisi dei gruppi (cluster analysis o clustering) è concettualmente semplice: verificare la possibile esistenza, in un insieme di oggetti, di sottoinsiemi di oggetti particolarmente simili tra loro (gruppi/cluster).

Nonostante alla base dell'analisi dei gruppi vi sia un'idea semplice e logica, vi sono numerosi modi per realizzarla [1], qui ci occupiamo dell'implementazione del clustering gerarchico con metodi esclusivi nelle due versioni:

→ clustering agglomerativo o bottom-up (con il metodo di Ward e con AGglomerative NESting, AGNES);

→ clustering divisivo o top-down (con DIvisive ANAlysis, DIANA).

Come dati impieghiamo i valori di BMI (indice di massa corporea) rilevati a livello europeo alcuni anni fa e pubblicati dall'Istat [2].

Per proseguire è necessario:

→ effettuare il download del file di dati bmi.csv

→ salvare il file nella cartella C:\Rdati\

Per il file di dati trovate link e modalità di download alla pagina Dati, ma potete anche semplicemente copiare i dati riportati qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarli in C:\Rdati\ in un file di testo denominato bmi.csv (assicuratevi che il file sia effettivamente salvato con l'estensione .csv).

Nazione;sottopeso;normale;sovrappeso;obeso

Austria;2.4;49.6;33.3;14.7

Belgio;2.7;48.0;35.3;14.0

Bulgaria;2.2;43.8;39.2;14.8

Cipro;3.9;47.8;33.8;14.5

Croazia;1.9;40.7;38.7;18.7

Danimarca;2.2;50.0;32.9;14.9

Estonia;2.2;43.9;33.5;20.4

Finlandia;1.2;44.1;36.4;18.3

Francia;3.2;49.6;31.9;15.3

Germania;1.8;46.1;35.2;16.9

Grecia;1.9;41.3;39.4;17.3

Irlanda;1.9;42.3;37.0;18.7

Lettonia;1.7;41.8;35.2;21.3

Lituania;1.9;42.5;38.3;17.3

Lussemburgo;2.8;49.3;32.4;15.6

Malta;2.0;37.0;35.0;26.0

Olanda;1.6;49.0;36.0;13.3

Polonia;2.4;42.9;37.5;17.2

Portogallo;1.8;44.6;36.9;16.6

Regno Unito;2.1;42.2;35.6;20.1

Repubblica Ceca;1.1;42.1;37.6;19.3

Romania;1.3;42.9;46.4;9.4

Slovacchia;2.1;43.6;38.0;16.3

Slovenia;1.6;41.8;37.4;19.2

Spagna;2.2;45.4;35.7;16.7

Svezia;1.8;48.3;35.9;14.0

Ungheria;2.9;41.9;34.0;21.2

Inoltre è necessario scaricare dal CRAN il pacchetto aggiuntivo cluster [3] e il pacchetto aggiuntivo dendextend [4]. 

Iniziamo con un esempio di custering gerarchico che impiega esclusivamente i valori di default. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING GERARCHICO con i valori di default

#  

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

windows() # apre e inizializza una nuova finestra grafica

#

mat <- dist(z) # calcola la matrice delle distanze

clus <- hclust(mat) # effettua il clustering gerarchico

plot(clus) # traccia il dendrogramma

#  

Con la prima riga di codice i dati sono importati nell'oggetto mydata.

Con la funzione scale() della seconda riga per ciascun dato viene calcolata la deviata normale standardizzata z. In pratica questa funzione prima calcola per i dati di ciascuna colonna/variabile la media e la deviazione standard, poi calcola per ciascuno dato x la corrispondente deviata normale standardizzata z come

z = (x – media) / deviazione standard

I valori di z sono salvati nel nuovo oggetto qui denominato per comodità mnemonica z. Se nella Console di R digitate z 

> z 

sono mostrati i dati standardizzati, in coda ai quali sono riportate  la media 

attr(,"scaled:center")

 sottopeso    normale sovrappeso      obeso 

  2.103704  44.537037  36.240741  17.111111 

e la deviazione standard

attr(,"scaled:scale")

 sottopeso    normale sovrappeso      obeso 

 0.6111033  3.3717318  2.8989732  3.2322097 

impiegate per effettuare la standardizzazione dei dati.

Dopo aver aperto e inizializzato una nuova finestra grafica con windows(), si passa al clustering gerarchico, realizzato in modo molto semplice, con tre sole funzioni, che impiegano per i loro argomenti i valori di default:

→ la funzione dist() che, sui valori standardizzati z, calcola la matrice delle distanze;

→ la funzione hclust() che dalla matrice delle distanze genera il clustering gerarchico; 

→ la funzione plot() che dal clustering gerarchico traccia il dendrogramma.

I principali argomenti delle tre funzioni sono evidenziati negli script successivi e al bisogno sono illustrati nella documentazione che si può richiamare con help(nomedellafunzione).

In questo secondo esempio viene illustrato il metodo agglomerativo forse più estesamente impiegato nella pratica, e cioè il clustering gerarchico con il metodo di Ward.

Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING GERARCHICO con distanza euclidea e metodo di Ward

#  

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

windows() # apre e inizializza una nuova finestra grafica

#

mat <- dist(z, method="euclidean") # calcola la matrice delle distanze 

clus <- hclust(mat, method="ward.D2") # effettua il clustering gerarchico 

plot(clus, cex=0.8, hang=-1, main="Dendrogramma - metodo di Ward", xlab="BMI nei Paesi europei", sub="", ylab="Distanza nei valori di BMI dei cluster") # traccia il dendrogramma 

#

rect.hclust(clus, k=4, border=c("goldenrod", "red", "green", "blue")) # evidenzia 4 gruppi/cluster 

#  

Qui le funzioni necessarie diventano quattro:

→ la matrice delle distanze mat è calcolata con la funzione dist() sui dati standardizzati z impiegando l'argomento method="euclidean" che prevede di misurare la distanza tra due punti con il teorema di Pitagora e può assumere in alternativa uno dei seguenti valori: "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski";

→ il clustering gerarchico clus è effettuato con la funzione hclust() sulla matrice delle distanze mat impiegando l'argomento method=ward.D2 che specifica come costruire i cluster, e può assumere in alternativa uno dei seguenti valori: "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC), "centroid" (= UPGMC), il valore di default per l'argomento method è "complete";

→ a partire dal clustering gerarchico contenuto nell'oggetto clus con la funzione plot() viene tracciato il dendrogramma. L'argomento hang=-1 fa si che rami del dendrogramma ed etichette risultino allineati in basso. Potete togliere questo argomento se preferite il posizionamento dei rami del dendrogramma e delle etichette previsto di default (vedere il dendrogramma precedente);

→ con la funzione rect.hclust() e l'argomento k=4 viene stabilito il numero dei gruppi da evidenziare nel dendrogramma, mentre i riquadri che li delimitano sono tracciati con i colori definiti nell'argomento border=c(...).

Questo è il dendrogramma ottenuto con il metodo di Ward:

Ora vediamo come il clustering agglomerativo può essere realizzato anche con il metodo di AGglomerative NESting sostituendo la funzione hclust() dello script precedente con la funzione agnes().

Copiate questo nuovo script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING GERARCHICO agglomerativo (AGglomerative NESting)

#  

library(cluster) # carica il pacchetto

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

windows() # apre e inizializza una nuova finestra grafica

#

mat <- dist(z, method="euclidean") # calcola la matrice delle distanze

clus <- agnes(mat, method="ward") # effettua il clustering gerarchico agglomerativo (AGNES)

plot(clus, which.plots=2, cex=0.8, hang=-1, main="Dendrogramma - clustering con AGglomerative NESting", xlab="BMI nei Paesi europei", sub="", ylab="Distanza nei valori di BMI dei cluster") # traccia il dendrogramma

#

rect.hclust(clus, k=4, border=c("goldenrod","blue","green","red")) # evidenzia 4 gruppi/cluster 

# 

Le differenze di codice necessarie per realizzare il clustering agglomerativo in questo script rispetto al precedente sono poche, ma ovviamente determinanti:

→ è richiesto il pacchetto cluster;

→ per il clustering gerarchico agglomerativo, viene impiegata la funzione agnes() che prevede per l'argomento method i valori "average", "single", "complete", "ward", "weighted", "gaverage";

→ nella funzione plot() l'argomento which.plot=2 serve per rappresentare solamente il dendrogramma, mentre l'argomento hang=-1 di nuovo allinea rami del dendrogramma ed etichette  in basso.

In effetti come vedete impiegando l'argomento method="ward" i cluster sono uguali a quelli ottenuti con il metodo di Ward dello script precedente.

Per realizzare il clustering divisivo (top-down) con la DIvisive ANAlysis impieghiamo quest'altro script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING GERARCHICO divisivo (DIvisive ANAlysis)

#  

library(cluster) # carica il pacchetto

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

windows() # apre e inizializza una nuova finestra grafica

#

mat <- dist(z, method="euclidean") # calcola la matrice delle distanze

clus <- diana(mat) # effettua il clustering gerarchico divisivo (DIANA)

plot(clus, which.plots=2, cex=0.8, hang=-1, main="Dendrogramma - clustering con DIvisive ANAlysis", xlab="BMI nei Paesi europei", sub="", ylab="Distanza nei valori di BMI dei cluster") # traccia il dendrogramma

#

rect.hclust(clus, k=4, border=c("goldenrod","blue","green","red")) # evidenzia 4 gruppi/cluster 

# 

Nella documentazione della funzione diana() per il clustering con la DIvisive ANAlysis si legge che questo metodo "... is probably unique in computing a divisive hierarchy, whereas most other software for hierarchical clustering is agglomerative". 

La funzione prevede un unico algoritmo di clusterizzazione, per cui si riduce all'espressione diana(mat) mentre il resto del codice, a parte naturamente titolo ed etichette, resta identico a quello dello script precedente.

Per evidenziare le differenze tra i risultati ottenuti con il clustering divisivo e quelli ottenuti con il precedente clustering agglomerativo, impieghiamo la possibilità piuttosto interessante di effettuare il confronto tra due dendrogrammi, sia graficamente sia calcolando un indice numerico di corrispondenza.

Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CLUSTERING GERARCHICO confronto tra i dendrogrammi AGNES e DIANA

#  

library(dendextend) # carica il pacchetto

library(cluster) # carica il pacchetto

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

z <- scale(mydata) # calcola la deviata normale standardizzata dei dati

windows() # apre e inizializza una nuova finestra grafica

#

mat <- dist(z, method="euclidean") # calcola la matrice delle distanze

clusAgnes <- agnes(mat, method="ward") # effettua il clustering gerarchico agglomerativo (AGNES)

clusDiana <- diana(mat) # effettua il clustering gerarchico divisivo (DIANA)

dendAgnes <- as.dendrogram(clusAgnes) # predispone il dendrogramma

dendDiana <- as.dendrogram(clusDiana) # predispone il dendrogramma

tanglegram(dendAgnes, dendDiana, main_left="AGNES", main_right="DIANA") # traccia il grafico di confronto tra i due dendrogrammi

entanglement(dendAgnes, dendDiana) # calcola la diffferenza tra i due dendrogrammi

#

Dopo avere caricato anche il pacchetto dendextend, i passaggi prevedono di:

→ calcolare come sempre la matrice delle distanze con la funzione dist(); 

→ effettuare il clustering gerarchico con i due algoritmi che vogliamo confrontare, che sono rispettivamente agnes() per il clustering agglomerativo e diana() per il clustering divisivo;

→ con la funzione as.dendrogram() estrarre e salvare i due dendrogrammi che vogliamo confrontare;

→ impiegare la funzione tanglegram() per tracciare il grafico che confronta visivamente i due dendrogrammi;

→ impiegare la funzione entanglement() per calcolare l'indice numerico che misura la differenza tra due dendrogrammi.

 

Il grafico risultante

mette a confronto i due dendrogrammi.

L'indice numerico che misura la differenza tra i due dendrogrammi, che nella documentazione del pacchetto è denominato "entanglement", ha un valore compreso tra 0 e 1, dove:

→ 0 indica che non vi è alcuna differenza tra i due dendrogrammi;

→ 1 indica che i due dendrogrammi non hanno nulla in comune e sono totalmente differenti.

Nel nostro caso abbiamo

> entanglement(dendAgnes, dendDiana) # calcola la differenza tra i due dendrogrammi

[1] 0.09419058

quindi il valore ottenuto 0.09419058 conferma che i due dendrogrammi sono molto simili, anche se determinano un raggruppamento in cluster lievemente diverso.

Trovate il seguito e le strategie alternative di clustering e di analisi dei dati multivariati nei post:

→ Analisi delle componenti principali 

→ Analisi dei gruppi (clustering non gerarchico) 

→ Analisi dei gruppi (clustering non esclusivo) 

----------

[1] Nonostante alla base dell'analisi dei gruppi vi sia un'idea semplice e logica, vi sono numerosi modi per realizzarla, infatti esistono:
→ metodi gerarchici, che danno luogo a una suddivisione ad albero (dendrogramma) in base alla distanza tra i singoli oggetti dell'insieme;
→ metodi non gerarchici, nei quali l'appartenenza di un oggetto (dell'insieme) ad uno specifico sottoinsieme/gruppo/cluster viene stabilita sulla base della sua distanza dal centro o dalla media dei dati o da un punto rappresentativo del cluster;
→ metodi bottom-up noti anche come metodi agglomerativi nei quali all'inizio del processo di classificazione ad ogni oggetto viene fatto corrispondere un cluster. In questo stadio gli oggetti sono considerati tutti dissimili tra di loro. Al passaggio successivo i due oggetti più simili sono raggruppati nello stesso cluster. Il numero dei cluster risulta quindi pari al numero di oggetti diminuito di uno. Il procedimento viene ripetuto ciclicamente, fino ad ottenere (all'ultimo passaggio) un unico cluster;
→ metodi top-down noti anche come metodi divisivi nei quali inizialmente tutti gli oggetti sono considerati come appartenenti ad un unico cluster, che viene via via suddiviso in cluster fino ad avere un numero di cluster uguale al numero degli oggetti;

→ metodi esclusivi, che prevedono che un oggetto possa appartenere esclusivamente a un cluster;
→ metodi non esclusivi (fuzzy) che prevedono che un oggetto possa appartenere, in modo quantitativamente diverso, a più di un cluster. 

[2] Vedere il post Indice di massa corporea (BMI).

[3] Trovate la documentazione nel manuale di riferimento del pacchetto: Package 'cluster'.

https://cran.r-project.org/web/packages/cluster/cluster.pdf

[4] Trovate la documentazione nel manuale di riferimento del pacchetto: Package ‘dendextend’.

https://cran.r-project.org/web/packages/dendextend/dendextend.pdf

Calcolare la media mobile e la mediana mobile

Data una sequenza di valori numerici la media mobile è la media calcolata in corrispondenza di ciascun dato su un numero fisso di dati, determinato da una finestra di ampiezza prestabilita che scorre i dati e nella quale il dato più vecchio viene ogni volta sostituito con un nuovo dato.

Così ad esempio per una finestra con l'ampiezza di 5 dati, in corrispondenza dei primi 4 dati la media mobile non può essere calcolata, in corrispondenza del dato numero 5 viene calcolata la media dei dati da 1 a 5, in corrispondenza del dato numero 6 viene calcolata la media dei dati da 2 a 6, in corrispondenza del dato numero 7 viene calcolata la media dei dati da 3 a 7 e così via. La mediana mobile a sua volta è la mediana calcolata in corrispondenza di ciascun dato nello stesso modo.

Per calcolare la media mobile e la mediana mobile facciamo ricorso al pacchetto zoo, che deve essere scaricato dal CRAN.

Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# CALCOLARE LA MEDIA MOBILE E LA MEDIANA MOBILE 

#

library(zoo) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediamob <- rollmean(valore_y, k=7, fill=NA, align='right') # calcola la media mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata$medianamob <- rollmedian(valore_y, k=7, fill=NA, align='right') # calcola la mediana mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media mobile e mediana mobile", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone il grafico della media mobile

lines(mydata$rileva_x, mydata$medianamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone il grafico della mediana mobile

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media mobile dei valori osservati", "Mediana mobile dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto zoo i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media mobile e la mediana mobile;

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna/vettore mydata$mediamob nella quale è riportata la media mobile calcolata in corrispondenza di ogni dato;

→ la colonna/vettore mydata$medianamob nella quale è riportata la mediana mobile calcolata in corrispondenza di ogni dato.

La media mobile e la mediana mobile sono calcolate rispettivamente con le funzioni rollmean() e rollmedian() che impiegano i seguenti argomenti:

→ la variabile (valore_y) sulla quale calcolare la media mobile;

→ l'ampiezza della finestra (k=7) che specifica che la media (o mediana) mobile va calcolata, in corrispondenza di ogni dato, su quel dato e sui sei dati che lo precedono;

→ fill=NA che specifica di riportare NA in corrispondenza dei valori non computabili;

→ align='right' che - in quanto con k=7 la media e la mediana possono essere calcolate solamente a partire dal settimo dato - specifica di riportare i valori non computabili NA nelle prime sei posizioni.

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto zoo [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media mobilea (mediamob) e mediana mobile (medianamob)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y mediamob medianamob

1         1        5       NA         NA

2         2        9       NA         NA

3         3       12       NA         NA

4         4        7       NA         NA

5         5        8       NA         NA

6         6       13       NA         NA

7         7       26 11.42857          9

8         8        7 11.71429          9

9         9        8 11.57143          8

10       10       19 12.57143          8

11       11       10 13.00000         10

12       12        7 12.85714         10

13       13        5 11.71429          8

14       14       36 13.14286          8

15       15       12 13.85714         10

16       16       18 15.28571         12

17       17       31 17.00000         12

18       18       17 18.00000         17

19       19       16 19.28571         17

20       20       43 24.71429         18

21       21        7 20.57143         17

22       22       12 20.57143         17

23       23        9 19.28571         16

24       24       27 18.71429         16

25       25       15 18.42857         15

26       26       18 18.71429         15

27       27       24 16.00000         15

28       28       29 19.14286         18

29       29       32 22.00000         24

30       30       23 24.00000         24

31       31       29 24.28571         24

32       32       12 23.85714         24

33       33        7 22.28571         24

34       34       48 25.71429         29

35       35       34 26.42857         29

36       36       45 28.28571         29

37       37       19 27.71429         29

38       38       12 25.28571         19

39       39       16 25.85714         19

40       40       18 27.42857         19

41       41       22 23.71429         19

42       42       31 23.28571         19

43       43       38 22.28571         19

44       44        7 20.57143         18

45       45       11 20.42857         18

46       46        4 18.71429         18

47       47       38 21.57143         22

48       48       36 23.57143         31

49       49       43 25.28571         36

50       50       15 22.00000         15

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media mobile con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana mobile con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplicativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto zoo: S3 Infrastructure for Regular and Irregular Time Series (Z's Ordered Observations).

https://cran.r-project.org/web/packages/zoo/index.html

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.