CLUSTER ANALYSIS

Calcolare la distanza tra le osservazioni

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

EUCLIDEAN DISTANCE

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.

CORRELATION BASED DISTANCES

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

DISTANCES FOR MIXED DATA

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

VISUALIZZARE LE MATRICI DI DISTANZA

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

PARTITIONING CLUSTERING

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

K-Means Clustering

è 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

COMPUTING K-MEANS CLUSTERING

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

K-medoids

è 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:

  1. seleziona k oggetti per diventare i medoids

  2. calcola la matrice di dissimilarità se non viene fornita

  3. assegna ogni oggetto a un cluster

  4. 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

  5. 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

COMPUTE PAM

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.

Clara: Clustering Large Applications

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())

```