Capitolo 7 Clustering Non Gerarchico (K-Means) (NHC)
7.1 Introduzione
Dato un numero preassegnato \(K\) di gruppi, i metodi di clustering non gerarchico (NHC) cercano la paritizione dei dati (“righe” o item) in \(K\) cluster in maniera tale che gli item entro ciascun cluster o gruppo siano quanto possibile simili tra loro, mentre gli item appartenenti a cluster diversi siano quanto possibile diversi.
Un possibile approccio per ottenere questo risultato potrebbe passare attraverso l’elencazione di tutti i possibili raggruppamenti in \(K\) gruppi costruibili con gli item e quindi scegliere come migliore soluzione il raggruppamento che ottimizza un qualche criterio predefinito. Sfortunatamente, un tale approccio diventerebbe rapidamente inapplicabile, specialmente per grandi dataset, poiché richiederebbe una quantità enorme di tempo macchina e di spazio di memoria. Di conseguenza tutte le tecniche di clustering disponibili sono iterative e operano solo su un numero molto ristretto di enumerazioni.
Tra i molti algoritmi di clustering non gerarchici sviluppati finora, il più popolare
è il K-means. Nella sua implementazione fondamentale, l’algoritmo
K-means inizia assegnando gli item a uno dei \(K\) cluster predeterminati
e quindi calcolando i \(K\) centroidi di gruppo, oppure pre-specificando
i \(K\) centroidi di gruppo. I centroidi pre-specificati possono essere item selezionati casualmente
oppure possono essere ottenuti “taglliando” un dendrogramma ad una altezza appropriata.
In seguito, tramite una procedura iterativa, l’algoritmo cerca di minimizzare la somma dei
quadrati entro gruppi (WGSS) su tutte le variabili, riassegnando (“spostando”) gli item sui diversi
cluster. La procedura si ferma quando non si ottengono miglioramenti di WGSS con ulteriori riassegnamenti.
La soluzione così ottenuta può non essere unica: l’algoritmo potrebbe trovare solo un minimo locale di WGSS. Si suggerisce pertanto di lanciare più volte l’algoritmo con assegnazioni casuali iniziali differenti degli item ai \(K\) cluster (oppure selezionando casualmente i \(K\) centroidi iniziali), così da trovare un probabile minimo globale di WGSS e, quindi, la migliore soluzione di clustering basata su \(K\) cluster.
7.2 Esempio: Dati Utility
Come esempio, consideriamo ancora il dataset utilities:
library(dplyr)
library(ggplot2)
library(GGally)
## Dataset dal package
library(pdataita)
data(utilities)
set.seed(10) # Usato per riproducibilità
utilities_tmp <- utilities[, 1:8]
# Matrice di scatterplot
ggpairs(utilities_tmp)
Se calcoliamo prima le varianze delle variabili troviamo i seguenti:
## coverage return cost load peak sales
## 3.404437e-02 5.035758e+00 1.696727e+03 1.990184e+01 9.723485e+00 1.260239e+07
## nuclear fuel
## 2.819686e+02 3.092451e-01Le varianze sono molto diverse, e l’uso del K-means sui dati grezzi potrebbe fornire
risutati poco utili perché “dominati” dalle caratteristiche più variabili.
Come nel caso del clustering gerarchico, quindi, dovremo normalizzare i dati in qualche modo;
in questo caso scegliamo di normalizzare ogni variabile dividendone i valori per il rispettivo range. Dopo la normalizzazione
le varianze diventano:
# calcola il range di ogni variabile
rge <- sapply(utilities_tmp, function(x) diff(range(x)))
# Divide ogni valore di ogni variabile per il suo range
utilities_s <- sweep(x = utilities_tmp, MARGIN = 2, STATS = rge, FUN = "/")
# Calcola la varianza di tutte le variabili "normalizzate"
sapply(utilities_s, var)## coverage return cost load peak sales nuclear
## 0.06217015 0.06216985 0.06972088 0.06281353 0.07481906 0.06302205 0.11189051
## fuel
## 0.09470796Molto più confrontabili.
Possiamo quindi procedere con il clustering dei dati normalizzati. Dapprima tracciamo i valori di WGSS per soluzioni di
raggruppamento con un numero di cluster che varia tra 1 e 8 per vedere se possiamo ottenere indicazioni sul
numero di gruppi (che usualmente non è noto a priori). Il grafico può essere ottenuto come segue:
k_max <- 8
wss <- sapply(1:k_max,
function(k,data) kmeans(data, centers = k)$tot.withinss,
data=utilities_s)
ggp <- ggplot(data = data.frame(x=1:k_max, y=wss), mapping = aes(x=x,y=y)) +
geom_point(colour = "red") +
geom_line(colour = "blue") +
xlab("Numero di gruppi") +
ylab("Somma dei quadrati entro gruppi (WGSS)")
print(ggp)
Mano a mano che il numero di gruppi aumenta, la somma dei quadrati necessariamente
si riduce; se si presenta un “gomito” nel grafico, questo potrebbe essere
indicativo della soluzione più utile che l’analista dovrebbe osservare in dettaglio. Nel nostro caso
un possibile “gomito” (anche se non particolarmente evidente) potrebbe essere nella soluzione con
quattro gruppi.
Analizzeremo quindi questa soluzione. Le medie di gruppo per questa soluzione sono
calcolabili tramite:
## coverage return cost load peak sales
## 1 1.07000 11.725714 0.8477473 38.774639 0.1993734 4.723675e+00
## 2 12.20270 13931.503704 12.8846154 43561.694757 5.0000000 1.550467e+04
## 3 254.47876 69.642540 127.5714286 156.441894 23.4586466 3.607876e+01
## 4 26.74811 2.304929 20.2189744 5.622003 5.8708772 9.567733e-01
## nuclear fuel
## 1 0.1055037 11.34142
## 2 0.0000000 4426.72514
## 3 9.9886170 24.29237
## 4 13.5733865 0.77160Volendo, possiamo provare l’indice di Calinski-Harabasz anche con il clustering non gerarchico, per trovare il numero “ottimale” di cluster:
require(clusterSim)
minC <- 2
maxC <- 10
res <- sapply(minC:maxC,
function(nc, data) index.G1(data, kmeans(data,centers = nc)$cluster),
data = utilities_s)
ggp <- ggplot(data=data.frame(x=2:(length(res)+1), y= res), mapping = aes(x=x,y=y)) +
geom_point() +
geom_line() +
xlab("Numero di cluster") +
ylab("Statistica pseudo-F di Calinski-Harabasz")
print(ggp)
Questo indice fornisce anche in questo caso una chiara indicazione che il numero di cluster da scegliere debba essere uguale a 4.
Ora, con l’obiettivo di spiegare o comprendere le caratteristiche e i significati dei gruppi, possiamo applicare i metodi PCA per ridurre la dimensionalità dei dati:
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 1.4740918 1.3785018 1.1504236 0.9983701 0.80561801
## Proportion of Variance 0.2716183 0.2375334 0.1654343 0.1245929 0.08112755
## Cumulative Proportion 0.2716183 0.5091517 0.6745860 0.7991789 0.88030644
## Comp.6 Comp.7 Comp.8
## Standard deviation 0.75608141 0.46529886 0.4115657
## Proportion of Variance 0.07145739 0.02706288 0.0211733
## Cumulative Proportion 0.95176383 0.97882670 1.0000000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## coverage 0.446 0.232 0.555 0.401 0.206 0.481
## return 0.571 0.101 0.332 -0.336 0.133 -0.150 -0.629
## cost -0.349 -0.161 0.467 0.409 0.269 -0.538 -0.118 -0.303
## load -0.289 0.409 -0.143 0.334 -0.680 -0.299 0.248
## peak -0.355 -0.283 0.281 0.391 -0.163 0.719 0.122
## sales -0.603 -0.332 0.191 -0.132 -0.150 0.661 -0.103
## nuclear 0.168 0.738 -0.333 -0.250 0.489
## fuel -0.336 0.540 -0.134 0.293 0.252 0.489 -0.433
Applicando il criterio della “proporzione di varianza spiegata”, il numero di componenti da tenere dovrebbe essere almeno pari a 4. Anche il criterio di Kaiser indicherebbe un numero di componenti pari a 4 (si veda il grafico qui sopra). In ogni modo, per questo esempio, lavoreremo con le sole prime due componenti, e quindi analizzeremo la correlazione tra gli score e le variabili originarie.
La prima componente principale sembra contrapporre misure di “uscita” (fuel, load, cost) con misure di valore più finanziario (coverage and return), mentre la seconda componente principale sembra contrapporre le vendite (sales) con i costi (fuel).
Proviamo quindi a tracciare le informazioni sulla suddivisione in gruppi sul grafico degli score della semplice analisi PCA a 2 dimensioni.
Cluster <- as.character(km_4$cluster)
utilities_scores <- cbind(data.frame(utilities_pca$scores[, c("Comp.1", "Comp.2")]),
company=utilities$comp_short,
Cluster=Cluster)
ggp <- ggplot(data= utilities_scores,
mapping = aes(x = Comp.1, y=Comp.2, label=company, colour=Cluster)) +
geom_point() +
xlab("Prima dimensione PCA") +
ylab("Seconda dimensione PCA") +
geom_text(hjust=0.5, vjust=-0.5, size=3)
print(ggp)
Questo grafico mostra come i 4 gruppi si dispongano all’interno delle prime due componenti della PCA. Il grafico mostra che il gruppo 2 e 3 sono ben separati.
Gli altri due gruppi sono in qualche modo sovrapposti; probabilmente potrebbero risultare meglio separati usando la terza componente principale.
Senza grafici 3D, il grafico che segue può aiutare nel trovare alcune differenze tra i rimanenti due gruppi:
Il prossimo grafico aiuterà ad “identificare” i gruppi come è stato fatto per il clustering gerarchico:
utilities$member <- Cluster
util.summ <- summarise(group_by(utilities, member),
coverage = mean(coverage),
return = mean(return),
cost = mean(cost),
load = mean(load),
peak = mean(peak),
sales = mean(sales),
nuclear = mean(nuclear),
fuel = mean(fuel))
palette(rainbow(8))
to.draw <- apply(util.summ[, -1], 2, function(x) x/max(x))
stars(to.draw, draw.segments = TRUE, scale = FALSE, key.loc = c(4.6,2.0), nrow=3,
ncol=2,
labels = c("CLUSTER 1", "CLUSTER 2","CLUSTER 3", "CLUSTER 4"),
main = "Dati utility (profili cluster)", cex = .75, ylim=c(0,8),
flip.labels = TRUE)
Si vede come il cluster 2 sia caratterizzato dall’uso elevato di energia fossile,
mentre il cluster 4 dall’uso importante di energia nucleare; il custer 3 è caratterizzato
da valori elevati di peak e di sales, mentre il cluster 1 tende ad avere valori bassi di
peak, sales, nuclear e fuel.