Il primo passo in una cluster analysis è scegliere in che modo misurare la distanza tra le osservazioni.
Quando calcoliamo la distanza tra le osservazioni, è preferibile che le variabili siano scalate in modo che siano più facilmente confrontabili.
Esistono diversi metodi per misurare la distanza tra le osservazioni.
I metodi classici sono la distanza Euclidea e la distanza di Manhattan. Esistono altre misure di dissimilarità che si basano sulla correlazione tra le variabili. La correlazione di Pearson è la più utilizzata tra le tecniche parametriche mentre la correlazione di Kendall e di Spearman tra le non parametriche (si utilizzano con variabili qualitative). Nella maggior parte dei software la distanza di default è la distanza Euclidea.
Se standardizzate le variabili, le distanze Euclidea,Manahattan e Correlation based diventano simili
set.seed(123)
summary(USArrests)
## Murder Assault UrbanPop Rape
## Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
## 1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
## Median : 7.250 Median :159.0 Median :66.00 Median :20.10
## Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
## 3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
## Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00
# prendiamo solo 15 casi per semplicitÃ
ss = sample(1:50, 15)
df = USArrests[ss,]
df.scaled = scale(df)
diversi pacchetti e funzioni calcolano le distanze tra le osservazioni
1. dist() del pacchetto stats accetta solo numeric come input
2. get_dist() del pacchetto factorextra anche, ma supporta anche le distanze correlation based quindi Pearson, Kendall e Sperman
3. daisy() del pacchetto cluster consente di calcolare il Gower coefficient, che è il più usato per misurare la distanza quando i dati sono di diversa classe
dist.eucl = dist(df.scaled, method = "euclidean")
# per visualizzare meglio i dati
round(as.matrix(dist.eucl), 1) # queste sono le distanze tra le righe
## Iowa Rhode Island Maryland Tennessee Utah Arizona Mississippi
## Iowa 0.0 2.8 4.1 3.3 2.4 4.3 3.9
## Rhode Island 2.8 0.0 3.6 3.7 2.0 3.3 4.4
## Maryland 4.1 3.6 0.0 1.5 3.0 1.3 2.4
## Tennessee 3.3 3.7 1.5 0.0 2.8 2.4 2.0
## Utah 2.4 2.0 3.0 2.8 0.0 2.5 4.2
## Arizona 4.3 3.3 1.3 2.4 2.5 0.0 3.6
## Mississippi 3.9 4.4 2.4 2.0 4.2 3.6 0.0
## Wisconsin 0.7 2.2 4.0 3.4 2.0 4.1 4.1
## Virginia 2.1 2.6 2.1 1.3 1.8 2.5 2.5
## Maine 0.7 3.0 4.2 3.5 2.9 4.6 3.8
## Texas 3.7 2.9 1.6 1.6 2.2 1.6 3.1
## Louisiana 3.9 3.5 1.2 1.2 3.1 2.2 1.8
## Montana 1.2 3.0 3.1 2.2 2.3 3.6 2.8
## Michigan 4.4 3.9 1.2 1.7 2.9 1.2 3.3
## Arkansas 2.4 3.3 2.2 1.4 2.7 3.0 1.8
## Wisconsin Virginia Maine Texas Louisiana Montana Michigan
## Iowa 0.7 2.1 0.7 3.7 3.9 1.2 4.4
## Rhode Island 2.2 2.6 3.0 2.9 3.5 3.0 3.9
## Maryland 4.0 2.1 4.2 1.6 1.2 3.1 1.2
## Tennessee 3.4 1.3 3.5 1.6 1.2 2.2 1.7
## Utah 2.0 1.8 2.9 2.2 3.1 2.3 2.9
## Arizona 4.1 2.5 4.6 1.6 2.2 3.6 1.2
## Mississippi 4.1 2.5 3.8 3.1 1.8 2.8 3.3
## Wisconsin 0.0 2.1 1.3 3.4 3.8 1.5 4.3
## Virginia 2.1 0.0 2.4 1.8 1.8 1.2 2.4
## Maine 1.3 2.4 0.0 4.0 4.0 1.3 4.7
## Texas 3.4 1.8 4.0 0.0 1.4 2.9 1.4
## Louisiana 3.8 1.8 4.0 1.4 0.0 2.8 1.8
## Montana 1.5 1.2 1.3 2.9 2.8 0.0 3.5
## Michigan 4.3 2.4 4.7 1.4 1.8 3.5 0.0
## Arkansas 2.6 1.1 2.3 2.6 2.0 1.2 2.8
## Arkansas
## Iowa 2.4
## Rhode Island 3.3
## Maryland 2.2
## Tennessee 1.4
## Utah 2.7
## Arizona 3.0
## Mississippi 1.8
## Wisconsin 2.6
## Virginia 1.1
## Maine 2.3
## Texas 2.6
## Louisiana 2.0
## Montana 1.2
## Michigan 2.8
## Arkansas 0.0
round(as.matrix(dist.eucl)[1:3,1:3], 1)
## Iowa Rhode Island Maryland
## Iowa 0.0 2.8 4.1
## Rhode Island 2.8 0.0 3.6
## Maryland 4.1 3.6 0.0
le colonne sono le variabili, per fare confronti pairwaise tra variabili bisogna trasporre i dati.
La distanza viene costruita sottraendo 1 al coefficiente di correlazione
library(factoextra)
dist.cor = get_dist(df.scaled, method = "pearson")
round(as.matrix(dist.cor)[1:3,1:3],1)
## Iowa Rhode Island Maryland
## Iowa 0.0 0.4 1.9
## Rhode Island 0.4 0.0 1.5
## Maryland 1.9 1.5 0.0
la funzione daisy fornisce una soluzione se ci sono variabili non numeriche, accetta anche variabili qualitative, e ordinali.
library(cluster)
data("flower")
str(flower)
## 'data.frame': 18 obs. of 8 variables:
## $ V1: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 2 2 ...
## $ V2: Factor w/ 2 levels "0","1": 2 1 2 1 2 2 1 1 2 2 ...
## $ V3: Factor w/ 2 levels "0","1": 2 1 1 2 1 1 1 2 1 1 ...
## $ V4: Factor w/ 5 levels "1","2","3","4",..: 4 2 3 4 5 4 4 2 3 5 ...
## $ V5: Ord.factor w/ 3 levels "1"<"2"<"3": 3 1 3 2 2 3 3 2 1 2 ...
## $ V6: Ord.factor w/ 18 levels "1"<"2"<"3"<"4"<..: 15 3 1 16 2 12 13 7 4 14 ...
## $ V7: num 25 150 150 125 20 50 40 100 25 100 ...
## $ V8: num 15 50 50 50 15 40 20 15 15 60 ...
dd = daisy(flower)
round(as.matrix(dd)[1:3,1:3],1)
## 1 2 3
## 1 0.0 0.9 0.5
## 2 0.9 0.0 0.5
## 3 0.5 0.5 0.0
fviz_dist(dist.eucl)
attenzione, queste distanze sono tre le osservazioni del dataframe ovvero sono relative alle righe non alle colonne. Il raggruppamento avviene tra le osservazioni e non tra le variabili
K-MEANS, PAM, CLARA sono algoritmi per classificare le osservazioni di un data set in gruppi; si basano sulla similarità .
Gli algoritmi richiedono che sia individuato un numero di cluster a priori da generare
CLARA è per dataset molto grandi
è l’algoritmo più usato per partizionare il dataset in gruppi. crea i cluster in modo che ogni gruppo sia composto da oggetti simili l’obiettivo è minimizzare la varianza intra-cluster
OGNI OSSERVAZIONE è ASSEGNATA A UN CLUSTER IN MODO CHE LA SOMMA DEI QUADRATI DELLE DISTANZE TRA LE OSSERVAZIONI E IL CENTRO DEL CLUSTER SIA MINIMA
Prima cosa, indicare il numero di cluster da generare
data("USArrests")
df = scale(USArrests)
head(df ,3)
## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
# kmeans(x deve essere numeric, centers, k o centri, iter max, nstart)
COME scelgo il numero di cluster??
questi sono tre modi per per individuare il numero di cluster.
library(ggplot2)
fviz_nbclust(df, kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)
fviz_nbclust(df, kmeans, method = "silhouette") +
geom_vline(xintercept = 4, linetype = 2)
fviz_nbclust(df, kmeans, method = "gap_stat") +
geom_vline(xintercept = 4, linetype = 2)
Il primo è un grafico che rappresenta la varianza within cluster. diminuisce con l’aumentare del numero di cluster, ma a un certo punto non cambia di molto da k = 4 in poi in questo caso
set.seed(123)
km.res = kmeans(df, 4, nstart = 25) # 4 è il numero di cluster
# nstart = 25 dice a R di provare 25 diversi punti di partenza e selezionare
# quello con la within cluster variance inferiore. più è grande più è stabile il risuktato
print(km.res)
## K-means clustering with 4 clusters of sizes 13, 16, 13, 8
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 2 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 3 0.6950701 1.0394414 0.7226370 1.27693964
## 4 1.4118898 0.8743346 -0.8145211 0.01927104
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 4 3 3 4 3
## Colorado Connecticut Delaware Florida Georgia
## 3 2 2 3 4
## Hawaii Idaho Illinois Indiana Iowa
## 2 1 3 2 1
## Kansas Kentucky Louisiana Maine Maryland
## 2 1 4 1 3
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 3 1 4 3
## Montana Nebraska Nevada New Hampshire New Jersey
## 1 1 3 1 2
## New Mexico New York North Carolina North Dakota Ohio
## 3 3 4 1 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 4
## South Dakota Tennessee Texas Utah Vermont
## 1 4 3 2 1
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 1 1 2
##
## Within cluster sum of squares by cluster:
## [1] 11.952463 16.212213 19.922437 8.316061
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"
aggregate(USArrests, by=list(cluster = km.res$cluster), mean)
## cluster Murder Assault UrbanPop Rape
## 1 1 3.60000 78.53846 52.07692 12.17692
## 2 2 5.65625 138.87500 73.87500 18.78125
## 3 3 10.81538 257.38462 76.00000 33.19231
## 4 4 13.93750 243.62500 53.75000 21.41250
# per aggiungere il cluster al data frame
dd = cbind(USArrests, cluster = km.res$cluster)
head(dd)
## Murder Assault UrbanPop Rape cluster
## Alabama 13.2 236 58 21.2 4
## Alaska 10.0 263 48 44.5 3
## Arizona 8.1 294 80 31.0 3
## Arkansas 8.8 190 50 19.5 4
## California 9.0 276 91 40.6 3
## Colorado 7.9 204 78 38.7 3
è una buona idea fare un grafico dei cluster per osservare meglio e per vedere se la scelta dei cluster è coerente
sarebbe interessante creare uno scatter plot e colorare ogni dato con il cluster
ma le variabili sono più di due, una possibile soluzione è ridurre il numero di dimensioni del dataset con un algoritmo di riduzione delle dimensioni
l’analisi delle componenti principali (PCA) riduce le dimensioni del dataset a due singole dimensioni
fviz_cluster(km.res, data = df,
palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "euclid", # Concentration ellipse
star.plot = TRUE, # Add segments from centroids to items
repel = TRUE, # Avoid label overplotting (slow)
ggtheme = theme_minimal()
)
k-means è veloce e semplice, ha delle debolezze però
assume una conoscenza dei dati a priori, perchè bisogna scegliere il numero di cluster
la soluzione finale è sensibile al centro del cluster scelto all’inizio
è sensibile agli outliers
se i dati sono ordinati diversamente potrebbero emergere soluzioni differenti
soluzioni: variare il numero di cluster e il centro iniziale e vedere la total within sum of square per scegliere la soluzione migliore
è correlato al k-means; il medoids è un data point che è più al centro del cluster di tutti gli altri, ce n’è uno per ogni cluster ed è l’esempio rappresentativo dei membri di quel cluster, questa cosa può essere utile in alcune situazioni
è più robusto del k-means perchè utilizza come centro del cluster non la media ma il medoids e quindi è meno sensibile agli outliers.
l’algoritmo necessita di specificare il numero di cluster che si può stabilire tramite i metodi specificati precedentemente
l’algoritmo k-medoids più usato è il PAM (partitioning around medoids)
L’algoritmo funziona così, in modo semplice:
seleziona k oggetti per diventare i medoids
calcola la matrice di dissimilarità se non viene fornita
assegna ogni oggetto a un cluster
per ogni cluter cerca se alcuni oggetti del cluster riducono il coefficiente totale di dissimilarità , se sì selezionano quell oggetto e lo rendono il medodide del cluster
se un medoide è cambiato torna a (3). altrimenti l’algoritmo si arresta
l’algoritmo PAM lavora con la matrice di dissimilarità , e per calcolare tale matrice si possono utilizzare due metriche: 1. la distanza euclidea, ovvero la radice quadrata della somma dei quadrati degli scarti 2. la distanza di Manhattan, che è la somma delle distanze assolute
Se i dati contengono OUTLIERS, la distanza di Manhattan è più robuta mentre la dista euclidea è un po’ più sensibile ai dati inusuali
data("USArrests")
df = scale(USArrests)
library("fpc","cluster")
library(factoextra)
fviz_nbclust(df, pam, method = "silhouette") +
theme_classic()
pam.res = pam(df, 2)
print(pam.res)
## Medoids:
## ID Murder Assault UrbanPop Rape
## New Mexico 31 0.8292944 1.3708088 0.3081225 1.1603196
## Nebraska 27 -0.8008247 -0.8250772 -0.2445636 -0.5052109
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 1 1 2 1
## Colorado Connecticut Delaware Florida Georgia
## 1 2 2 1 1
## Hawaii Idaho Illinois Indiana Iowa
## 2 2 1 2 2
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 1 2 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 2 1 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 2 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 1 2 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 1 2 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 2 2 2
## Objective function:
## build swap
## 1.441358 1.368969
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"
dd = cbind(USArrests, cluster = pam.res$cluster)
head(dd)
## Murder Assault UrbanPop Rape cluster
## Alabama 13.2 236 58 21.2 1
## Alaska 10.0 263 48 44.5 1
## Arizona 8.1 294 80 31.0 1
## Arkansas 8.8 190 50 19.5 2
## California 9.0 276 91 40.6 1
## Colorado 7.9 204 78 38.7 1
# visualizzare, anche qui riduce le dimensioni con la PCA. sceglie le due dimensioni
fviz_cluster(pam.res,
palette = c("#00AFBB", "#FC4E07"), # color palette
ellipse.type = "t", # Concentration ellipse
repel = TRUE, # Avoid label overplotting (slow)
ggtheme = theme_classic()
)
il k-medoids è un alternativa robusta al k-means, ogni cluster è rappresentato da un oggetto selezionato dentro il cluster che è il più rappresentativo del cluster. l’algoritmo necessita del numero di cluster per data set molto grandi il pam richiede troppa memoria e tempo di computazione.
set.seed(1234)
df = rbind(cbind(rnorm(200,0,8), rnorm(200,0,8)),
cbind(rnorm(300, 50, 8), rnorm(300,50,8)))
colnames(df) = c("x","y")
rownames(df) = paste0("S", 1:nrow(df))
head(df, 5)
## x y
## S1 -9.656526 3.881815
## S2 2.219434 5.574150
## S3 8.675529 1.484111
## S4 -18.765582 5.605868
## S5 3.432998 2.493448
# clara() è la funzione, ammette NAs
# numero ottimale di cluster
#### COMPUTE CLARA IN R ####
library(cluster)
library(factoextra)
fviz_nbclust(df, clara, method = "silhouette") +
theme_classic()
clara.res = clara(df, 2, samples = 50, pamLike = TRUE)
print(clara.res)
## Call: clara(x = df, k = 2, samples = 50, pamLike = TRUE)
## Medoids:
## x y
## S121 -1.531137 1.145057
## S455 48.357304 50.233499
## Objective function: 9.87862
## Clustering vector: Named int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "names")= chr [1:500] "S1" "S2" "S3" "S4" "S5" "S6" "S7" ...
## Cluster sizes: 200 300
## Best sample:
## [1] S37 S49 S54 S63 S68 S71 S76 S80 S82 S101 S103 S108 S109 S118
## [15] S121 S128 S132 S138 S144 S162 S203 S210 S216 S231 S234 S249 S260 S261
## [29] S286 S299 S304 S305 S312 S315 S322 S350 S403 S450 S454 S455 S456 S465
## [43] S488 S497
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"
dd = cbind(df, cluster = clara.res$cluster)
fviz_cluster(clara.res,
palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
geom = "point", pointsize = 1,
ggtheme = theme_classic())
```