Capitolo 6 Clustering Gerarchico (Agglomerativo) (HC)

6.1 Introduzione

Come già accennato nel precedente paragrafo, il clustering gerarchico (aggregativo, nel nostro caso) produce raggruppamenti di item (variabili o oggetti) seguendo un processo gerarchico: parte da tutti gli item separati in gruppi individuali, e quindi procede iterativamente con l’aggregazione di coppie di gruppi “che si somigliano di più”, arrivando al termine a produrre un unico gruppo contenente tutti gli item.

Il punto di partenza di ogni metodo di Clustering Gerarchico (HC), quindi, è il calcolo della dissimilarità (“dissomiglianza”; o “distanza”, se vale la disuguaglianza triangolare) di ciascun item rispetto a tutti gli altri. Quale definizione di distanza usare (euclidea, Manhattan, Canberra, ecc.), spesso dipende dall’applicazione specifica o è una scelta soggettiva. Alcune misure di distanza sono appropriate solo per certi tipi di dato, mentre altre sono state introdotte per raggruppare variabili (colonne di dati) invece che osservazioni (righe di dati).

Uno strumento generale in R per il calcolo delle dissimilarità è la funzione dist(). Per avere aiuti sulla funzione potete digitare:

?dist

nella console di R.
La funzione dist() produce un oggetto di classe dist, e sono proprio oggetti di classe dist quelli che devono essere forniti alla funzione per il clustering gerarchico.
Per ottenere le dissimilarità tra un insieme di item è anche possibile usare la funzione daisy() del package cluster. Questa funzione include una metrica (Gower) che permette di calcolare “dissimilarità” anche su misture di variabili quantitative e qualitative.

In ogni modo, dopo aver ottenuto la matrice di dissimilarità, per applicare un clustering gerarchico si deve specificare come si vuole siano calcolate le distanze tra gruppi di osservazioni (cluster) durante le iterazioni dell’algoritmo di clustering stesso. Anche in questo caso ci sono diverse scelte da fare. Tra queste, c’è la scelta del cosiddetto metodo di linkage (legame):

  • single linkage ==> (legame singolo) la distanza tra due gruppi è definita come il valore più piccolo delle distanze tra gli item dei due gruppi;
  • complete linkage ==> (legame completo) la distanza tra due gruppi è definita come il valore più grande di distanza tra gli item dei due gruppi;
  • average linkage ==> (legame medio) questo è un compromesso tra i due precedenti approcci, ottenuto come la media delle corrispondenti distanze.

Un altro approccio popolare per il calcolo delle distanze (o dissimilarità) tra gruppi è il Metodo di Ward, che unisce gruppi che ottengono la crescita minima di una data misura di eterogeneità (“varianza entro gruppi”).
Il criterio della varianza minima di Ward minimizza quindi la varianza entro-gruppi totale.
Per implementare questo metodo, ad ogni step si trova la coppia di cluster che porta alla crescita minima della varianza entro-gruppi totale dopo l’unione. Questo incremento è funzine della distanza quadratica ponderata tra centri dei cluster.
Nello step iniziale, tutti i cluster sono “singleton” (cluster contenenti un solo punto). Per applicare un algoritmo ricorsivo con questa funzione obiettivo, la distanza inziale tra singoli oggetti deve essere (proporzionale a) la distanza euclidea al quadrato.
Il metodo della minima varianza di Ward può essere definito ed implementato ricorsivamente tramite un algoritmo di Lance–Williams, dove, per cluster disgiunti \(C_i\), \(C_j\) e \(C_k\) con dimensioni \(n_i\), \(n_j\) e \(n_k\) rispettivamente, la distanza dal cluster “riunito” \(C_i \cup C_j\) ed un terzo cluster \(C_k\) è calcolata come:

\(d(C_i \cup C_j, C_k) = \frac{n_i+n_k}{n_i+n_j+n_k}\;d(C_i,C_k) + \frac{n_j+n_k}{n_i+n_j+n_k}\;d(C_j,C_k) - \frac{n_k}{n_i+n_j+n_k}\;d(C_i,C_j)\).

Tutti gli approcci mostrati qui sopra producono potenzialmente cluster diversi, e quale usare nella pratica dipende essenzialmente da una scelta soggettiva.

6.2 Esempio: Dati Utility

Come illustrazione, consideriamo il dataset utilities, che contiene informazioni su 22 aziende americane fornitrici di energia elettrica, con le seguenti variabili:

  • coverage ==> tasso di copertura (reddito/debito)
  • return ==> tasso di remunerazione del capitale
  • cost ==> costo per kW di capacità
  • load ==> fattore di carico annuale
  • peak ==> crescita del picco di dimanda kWh dal 1974 al 1975
  • sales ==> vendite (kWh uso per anno)
  • nuclear ==> percentuale di nucleare
  • fuel ==> costo totale del fonti energetiche (cent per kWh)
  • company ==> nome completo dell’azienda
  • comp_short ==> nome breve dell’azienda

Carichiamo quindi i dati che useremo per gli esempi, e rappresentiamoli graficamente come una matrice di scatterplot (si noti che rimuovo le ultime 2 colonne del data frame perché non numeriche).

## Classes 'tbl_df', 'tbl' and 'data.frame':    22 obs. of  10 variables:
##  $ coverage  : num  1.06 0.89 1.43 1.02 1.49 1.32 1.22 1.1 1.34 1.12 ...
##  $ return    : num  9.2 10.3 15.4 11.2 8.8 13.5 12.2 9.2 13 12.4 ...
##  $ cost      : int  151 202 113 168 192 111 175 245 168 197 ...
##  $ load      : num  54.4 57.9 53 56 51.2 60 67.6 57 60.4 53 ...
##  $ peak      : num  1.6 2.2 3.4 0.3 1 -2.2 2.2 3.3 7.2 2.7 ...
##  $ sales     : int  9077 5088 9212 6423 3300 11127 7642 13082 8406 6455 ...
##  $ nuclear   : num  0 25.3 0 34.3 15.6 22.5 0 0 0 39.2 ...
##  $ fuel      : num  0.628 1.555 1.058 0.7 2.044 ...
##  $ company   : chr  "Arizona Public Service" "Boston Edison Co." "Central Louis"..
##  $ comp_short: chr  "Arizona" "Boston" "Central" "Common" ...
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

Ci sono molte funzioni per eseguire l’analisi dei cluster in R (per un elenco completo, si veda la relativa task view del CRAN all’URL http://cran.r-project.org/web/views/Cluster.html). Qui descriveremo le funzioni basilari disponibili nella versione base di R, in particolare la funzione hclust().

Per evitare effetti di scala (cioè, dare più importanza alle variabili con valori assoluti o variabilità più elevati), standardizziamo prima le variabili, e quindi proviamo ad applicare il clustering gerarchico con i quattro metodi indicati sopra:

Come si vede, la funzione hclust() richiede in input una matrice di dissimilarità o distanza.

Una volta applicati i metodi di clustering come sopra, l’approccio grafico standard per rappresentare la soluzione di un algoritmo HC agglomerativo è quello del dendrogramma, che rappresenta la sequenza di fusioni successive dei diversi gruppi, nonché le distanze su cui sono state prodotte le aggregazioni Un dendrogramma in R è disponibile con il metodo plot per gli oggetti di classe (tipo) hclust:

Come si vede, il metodo del legame singolo tende a creare grandi gruppi aggiungendo in successione item a gruppi già creati (il cosiddetto effetto di “concatenazione”). Gli altri tre metodi ritornano invece risultati simili tra loro, con cluster più strutturati.

Una volta scelto il numero di gruppi (più avanti vedremo anche questo aspetto), si può usare la funzione cutree() per ottenere a quale gruppo (cluster) appartiene ciascun item:

##  [1] 1 1 1 1 1 1 2 3 1 1 3 2 1 1 2 3 2 1 1 1 2 1

Come si vede, il parametro k definisce il numero di gruppi che si vogliono ottenere.
L’appartenenza ai cluster, che otteniamo dalla chiamata della funzione cutree(), è utile per produrre una profilazione dei cluster, cioè per fornire una descrizione più dettagliata delle caratteristiche principali dei cluster identificati. Col prossimo grafico proviamo a visualizzare la profilezione:

## 
##  1  2  3 
## 14  5  3
## utilities$member: 1
##     coverage         return           cost            load            peak       
##  Min.   :0.890   Min.   : 8.80   Min.   : 96.0   Min.   :49.80   Min.   :-2.200  
##  1st Qu.:1.062   1st Qu.:10.53   1st Qu.:121.8   1st Qu.:53.17   1st Qu.: 1.100  
##  Median :1.135   Median :11.90   Median :159.5   Median :54.35   Median : 2.450  
##  Mean   :1.171   Mean   :11.71   Mean   :155.2   Mean   :55.31   Mean   : 2.429  
##  3rd Qu.:1.290   3rd Qu.:12.68   3rd Qu.:187.5   3rd Qu.:57.60   3rd Qu.: 3.475  
##  Max.   :1.490   Max.   :15.40   Max.   :202.0   Max.   :60.40   Max.   : 7.200  
##      sales          nuclear           fuel       
##  Min.   : 3300   Min.   : 0.00   Min.   :0.5270  
##  1st Qu.: 6636   1st Qu.: 0.00   1st Qu.:0.6300  
##  Median : 8742   Median :19.05   Median :0.7820  
##  Mean   : 8355   Mean   :18.20   Mean   :0.9699  
##  3rd Qu.: 9988   3rd Qu.:32.38   3rd Qu.:1.2078  
##  Max.   :13507   Max.   :50.20   Max.   :2.0440  
## ------------------------------------------------------------- 
## utilities$member: 2
##     coverage         return           cost            load            peak      
##  Min.   :0.760   Min.   : 6.40   Min.   :136.0   Min.   :61.00   Min.   :-0.10  
##  1st Qu.:0.960   1st Qu.: 7.60   1st Qu.:164.0   1st Qu.:61.90   1st Qu.: 2.20  
##  Median :1.040   Median : 8.60   Median :175.0   Median :62.00   Median : 3.50  
##  Mean   :1.022   Mean   : 9.14   Mean   :171.4   Mean   :62.94   Mean   : 3.66  
##  3rd Qu.:1.130   3rd Qu.:10.90   3rd Qu.:178.0   3rd Qu.:62.20   3rd Qu.: 3.70  
##  Max.   :1.220   Max.   :12.20   Max.   :204.0   Max.   :67.60   Max.   : 9.00  
##      sales         nuclear          fuel      
##  Min.   :5714   Min.   :0.00   Min.   :1.400  
##  1st Qu.:6154   1st Qu.:0.00   1st Qu.:1.652  
##  Median :6468   Median :0.00   Median :1.897  
##  Mean   :6526   Mean   :1.84   Mean   :1.797  
##  3rd Qu.:6650   3rd Qu.:0.90   3rd Qu.:1.920  
##  Max.   :7642   Max.   :8.30   Max.   :2.116  
## ------------------------------------------------------------- 
## utilities$member: 3
##     coverage         return           cost            load            peak      
##  Min.   :0.750   Min.   :7.500   Min.   :173.0   Min.   :51.50   Min.   :3.300  
##  1st Qu.:0.925   1st Qu.:8.350   1st Qu.:209.0   1st Qu.:53.75   1st Qu.:4.900  
##  Median :1.100   Median :9.200   Median :245.0   Median :56.00   Median :6.500  
##  Mean   :1.003   Mean   :8.867   Mean   :223.3   Mean   :54.83   Mean   :6.333  
##  3rd Qu.:1.130   3rd Qu.:9.550   3rd Qu.:248.5   3rd Qu.:56.50   3rd Qu.:7.850  
##  Max.   :1.160   Max.   :9.900   Max.   :252.0   Max.   :57.00   Max.   :9.200  
##      sales          nuclear       fuel       
##  Min.   :13082   Min.   :0   Min.   :0.3090  
##  1st Qu.:14536   1st Qu.:0   1st Qu.:0.4645  
##  Median :15991   Median :0   Median :0.6200  
##  Mean   :15505   Mean   :0   Mean   :0.5657  
##  3rd Qu.:16716   3rd Qu.:0   3rd Qu.:0.6940  
##  Max.   :17441   Max.   :0   Max.   :0.7680

Nel grafico qui sopra, la torta in basso a destra rappresenta la legenda per interpretare le “fette” delle altre 3 torte, che rappresentano i 3 gruppi scelti tremite la funzione cutree(). Per esempio, la “fetta” blu rappresenta il valore della variabile sales nei 3 gruppi.
La funzione cutree() può essere usata per generare anche più di un solo vettore alla volta di appartenenza ai gruppi, coe vediamo qui sotto:

##      clus4
## clus2 1 2 3 4
##     1 7 7 0 0
##     2 0 0 5 3

Un ulteriore approccio grafico per ricavare l’appartenenza ai cluster è dato dalla funzione rect.hclust(); questa aggiunge al dendrogramma dei rettangoli che mostrano i cluster identificati. Il grafico che segue identifica 5 gruppi con legame completo:

E’ possibile usare la funzione identify() per selezionare direttamente sul dendrogramma il numero di cluster da usare. Potete usare un codice come questo:

e poi fare clic sulle diramazioni del dendrogramma per ottenere i raggruppamenti degli item in una lista (in questo caso, nell’oggetto r).
Per uscire dalla modalità interattiva bisogna premere il tasto Esc della tastiera.

Infine, il package clusterSim fornisce molte funzioni per confrontare la bontà delle diverse soluzioni di clustering. Tra queste funzioni troviamo:

  • index.DB() ==> Calcola l’indice di Davies-Bouldin
  • index.G1() ==> Calcola la statistica pseudo-F di Calinski-Harabasz
  • index.G2() ==> Calcola l’indice G2 di qualità interna dei cluster
  • index.G3() ==> Calcola l’indice G3 di qualità interna dei cluster
  • index.Gap() ==> Calcola l’indice gap di Tibshirani, Walther e Hastie
  • index.H() ==> Calcola l’indice di Hartigan
  • index.KL() ==> Calcola l’indice di Krzanowski-Lai
  • index.S() ==> Calcola l’indice Silhouette di qualità interna dei cluster di Rousseeuw

Per avere maggiori dettagli si vedano le pagine di help di queste funzioni.
Come esempio, consideriamo la statistica pseudo-F di Calinski-Harabasz. Questa statistica dovrebbe ottenere il valore massimo in corispondenza della configurazione di gruppi migliore.
Proviamo quindi a rappresentare il valore della statistica al variare del numero di cluster individuati:

Seguendo i valori della statistica pseudo-F di Calinski-Harabasz, possiamo concludere che secondo questo indice la soluzione con 4 cluster è la migliore.

Un’ulteriore rappresentazione grafica che permette di valutare la qualità della rappresentazione dei cluster individuati è il “Silhouette plot”; in questo caso vogliamo rappresentare in forma grafica i valori di Silhouette per la soluzione a 4 gruppi individuata sopra.

##                  util_cl_m
##  [1,] "Arizona"  "1"      
##  [2,] "Boston"   "2"      
##  [3,] "Central"  "1"      
##  [4,] "Common"   "2"      
##  [5,] "Consolid" "2"      
##  [6,] "Florida"  "1"      
##  [7,] "Hawaiian" "3"      
##  [8,] "Idaho"    "4"      
##  [9,] "Kentucky" "1"      
## [10,] "Madison"  "2"      
## [11,] "Nevada"   "4"      
## [12,] "NewEngla" "3"      
## [13,] "Northern" "2"      
## [14,] "Oklahoma" "1"      
## [15,] "Pacific"  "3"      
## [16,] "Puget"    "4"      
## [17,] "SanDiego" "3"      
## [18,] "Southern" "1"      
## [19,] "Texas"    "1"      
## [20,] "Wisconsi" "2"      
## [21,] "United"   "3"      
## [22,] "Virginia" "2"

Il grafico qui sopra mostra i valori di silhouette per i singoli item. Ogni barra orizzontale rappresenta il valore di silhouette (\(s_i\)) di un item all’interno del suo cluster assegnato. \(s_i\) è sempre compreso tra -1 e 1. Un valore di \(s_i\) prossimo a 1 indica che l’item è assegnato ad un cluster che “lo rappresenta” molto bene. Un valore di \(s_i\) prossimo a -1 indica che l’item \(i\)-mo è più simile ad item che appartengono ad altri cluster vicini. Un \(s_i\) prossimo a zero indica che l’item \(i\)-mo si trova in una posizione sostanzialmente intermedia (cioè, sul bordo) tra due o più cluster. Nel nostro caso, sembra che la società alla riga 2 (“Boston”), cassificata come appartenente al secondo cluster, non sia “ben rappresentata” in quel gruppo. Potrebbe essere oggetto di analisi successive.

In ogni modo, dopo aver identificato una “buona” soluzione, possiamo salvare la corrispondente classificazione in cluster ed usarla a fini grafico/analitici:

E quindi produrre le statistiche descrittive e lo star plot per descrivere i cluster:

## 
## 1 2 3 4 
## 7 7 5 3
## utilities$member: 1
##     coverage         return           cost            load            peak       
##  Min.   :1.050   Min.   : 9.20   Min.   : 96.0   Min.   :49.80   Min.   :-2.200  
##  1st Qu.:1.075   1st Qu.:11.85   1st Qu.:107.5   1st Qu.:53.50   1st Qu.:-0.350  
##  Median :1.160   Median :12.60   Median :113.0   Median :54.40   Median : 1.600  
##  Mean   :1.207   Mean   :12.49   Mean   :127.6   Mean   :55.47   Mean   : 1.714  
##  3rd Qu.:1.330   3rd Qu.:13.25   3rd Qu.:150.5   3rd Qu.:58.35   3rd Qu.: 3.050  
##  Max.   :1.430   Max.   :15.40   Max.   :168.0   Max.   :60.40   Max.   : 7.200  
##      sales          nuclear            fuel       
##  Min.   : 8406   Min.   : 0.000   Min.   :0.5880  
##  1st Qu.: 9144   1st Qu.: 0.000   1st Qu.:0.6320  
##  Median : 9673   Median : 0.000   Median :0.8620  
##  Mean   :10163   Mean   : 3.214   Mean   :0.8744  
##  3rd Qu.:10634   3rd Qu.: 0.000   3rd Qu.:1.0830  
##  Max.   :13507   Max.   :22.500   Max.   :1.2410  
## ------------------------------------------------------------- 
## utilities$member: 2
##     coverage         return           cost            load            peak      
##  Min.   :0.890   Min.   : 8.80   Min.   :148.0   Min.   :51.20   Min.   :0.300  
##  1st Qu.:1.045   1st Qu.: 9.80   1st Qu.:171.0   1st Qu.:53.35   1st Qu.:1.600  
##  Median :1.120   Median :11.20   Median :192.0   Median :54.30   Median :2.700  
##  Mean   :1.134   Mean   :10.93   Mean   :182.9   Mean   :55.14   Mean   :3.143  
##  3rd Qu.:1.175   3rd Qu.:12.10   3rd Qu.:198.0   3rd Qu.:56.95   3rd Qu.:4.700  
##  Max.   :1.490   Max.   :12.70   Max.   :202.0   Max.   :59.90   Max.   :6.400  
##      sales          nuclear           fuel       
##  Min.   : 3300   Min.   :15.60   Min.   :0.5270  
##  1st Qu.: 5756   1st Qu.:25.95   1st Qu.:0.6615  
##  Median : 6455   Median :34.30   Median :0.7020  
##  Mean   : 6546   Mean   :33.19   Mean   :1.0653  
##  3rd Qu.: 7233   3rd Qu.:40.15   3rd Qu.:1.4305  
##  Max.   :10093   Max.   :50.20   Max.   :2.0440  
## ------------------------------------------------------------- 
## utilities$member: 3
##     coverage         return           cost            load            peak      
##  Min.   :0.760   Min.   : 6.40   Min.   :136.0   Min.   :61.00   Min.   :-0.10  
##  1st Qu.:0.960   1st Qu.: 7.60   1st Qu.:164.0   1st Qu.:61.90   1st Qu.: 2.20  
##  Median :1.040   Median : 8.60   Median :175.0   Median :62.00   Median : 3.50  
##  Mean   :1.022   Mean   : 9.14   Mean   :171.4   Mean   :62.94   Mean   : 3.66  
##  3rd Qu.:1.130   3rd Qu.:10.90   3rd Qu.:178.0   3rd Qu.:62.20   3rd Qu.: 3.70  
##  Max.   :1.220   Max.   :12.20   Max.   :204.0   Max.   :67.60   Max.   : 9.00  
##      sales         nuclear          fuel      
##  Min.   :5714   Min.   :0.00   Min.   :1.400  
##  1st Qu.:6154   1st Qu.:0.00   1st Qu.:1.652  
##  Median :6468   Median :0.00   Median :1.897  
##  Mean   :6526   Mean   :1.84   Mean   :1.797  
##  3rd Qu.:6650   3rd Qu.:0.90   3rd Qu.:1.920  
##  Max.   :7642   Max.   :8.30   Max.   :2.116  
## ------------------------------------------------------------- 
## utilities$member: 4
##     coverage         return           cost            load            peak      
##  Min.   :0.750   Min.   :7.500   Min.   :173.0   Min.   :51.50   Min.   :3.300  
##  1st Qu.:0.925   1st Qu.:8.350   1st Qu.:209.0   1st Qu.:53.75   1st Qu.:4.900  
##  Median :1.100   Median :9.200   Median :245.0   Median :56.00   Median :6.500  
##  Mean   :1.003   Mean   :8.867   Mean   :223.3   Mean   :54.83   Mean   :6.333  
##  3rd Qu.:1.130   3rd Qu.:9.550   3rd Qu.:248.5   3rd Qu.:56.50   3rd Qu.:7.850  
##  Max.   :1.160   Max.   :9.900   Max.   :252.0   Max.   :57.00   Max.   :9.200  
##      sales          nuclear       fuel       
##  Min.   :13082   Min.   :0   Min.   :0.3090  
##  1st Qu.:14536   1st Qu.:0   1st Qu.:0.4645  
##  Median :15991   Median :0   Median :0.6200  
##  Mean   :15505   Mean   :0   Mean   :0.5657  
##  3rd Qu.:16716   3rd Qu.:0   3rd Qu.:0.6940  
##  Max.   :17441   Max.   :0   Max.   :0.7680

Dal grafico sopra possiamo identificare il cluster 2 come quello caratterizato da elevata quantità di energia proveniente dal nucleare (nuclear), il cluster 2 come quello con grandi costi per combustibili (fuel), il cluster 4 come quello caratterizzato da elevati livelli di valori in vendita (sales) e di picco di domanda (peak), il cluster 1 come quello caratterizzato da livelli bassi di peak, nuclear, fuel (fonti energetiche alternative?).