Данные

library(psych)
library(cluster)
library(mclust)
library(PerformanceAnalytics)

df <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", sep = ",")
colnames(df) <- c("type", "alcohol", "malic", "ash", "alcalinity", "magnesium", "phenols", "flavanoids", "nonflavanoids", 
                  "proanthocyanins", "color", "hue", "dilution", "proline")

head(df)
##   type alcohol malic  ash alcalinity magnesium phenols flavanoids
## 1    1   14.23  1.71 2.43       15.6       127    2.80       3.06
## 2    1   13.20  1.78 2.14       11.2       100    2.65       2.76
## 3    1   13.16  2.36 2.67       18.6       101    2.80       3.24
## 4    1   14.37  1.95 2.50       16.8       113    3.85       3.49
## 5    1   13.24  2.59 2.87       21.0       118    2.80       2.69
## 6    1   14.20  1.76 2.45       15.2       112    3.27       3.39
##   nonflavanoids proanthocyanins color  hue dilution proline
## 1          0.28            2.29  5.64 1.04     3.92    1065
## 2          0.26            1.28  4.38 1.05     3.40    1050
## 3          0.30            2.81  5.68 1.03     3.17    1185
## 4          0.24            2.18  7.80 0.86     3.45    1480
## 5          0.39            1.82  4.32 1.04     2.93     735
## 6          0.34            1.97  6.75 1.05     2.85    1450

Это данные химического состава вин, которые росли в одном регионе Италии, но разных сортов. Есть 3 сорта вин.

pairs.panels(df, lm = FALSE, bg=c("red","yellow","blue")[df$type], pch=21, lwd = 0)

Matrix plot с раскраской по сортам.

pairs.panels(df, lm = FALSE)

Без раскраски.

Прологорифмируем некоторые переменные

df$magnesium <- log(df$magnesium)
df$flavanoids <- log(df$flavanoids)
df$nonflavanoids <- log(df$nonflavanoids)
df$color <- log(df$color)
df$proline <- log(df$proline)
pairs.panels(df, lm = FALSE)

Кластеризация по двум признакам

Для начала расмотрим кластеризацию по двум признакам: Malic acid (яблочная кислота) и Phenols. (Эти признаки были выбраны потому, что визуально видно неоднородности в плотностях и на скатер плоте выделяются кластеры)

#### clustering on 2 parameters ####
df2 <- df_stand[,c(2,6)]
plot(df2[,1], df2[,2], xlab = "Malic", ylab = "Phenols")

Представление данных в плосткости двух признаков. Видно 2 кластера

#### clustering on 2 parameters ####
df2 <- df_stand[,c(2,6)]
plot(df2[,1], df2[,2], col = df$type ,xlab = "Malic", ylab = "Phenols")

С раскараской по сортам.

K-Means

Раздели данные на 2 кластера, которые выделяются визульно.

set.seed(5)
kmeans_fit2_2 <- kmeans(df2, 2) # k = 2
kmeans_fit2_2
## K-means clustering with 2 clusters of sizes 65, 113
## 
## Cluster means:
##        malic    phenols
## 1  0.9978360 -0.8525573
## 2 -0.5739765  0.4904091
## 
## Clustering vector:
##   [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [36] 2 2 2 2 1 2 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [71] 1 2 2 2 2 1 2 1 2 1 2 2 2 1 2 2 2 2 2 2 1 1 1 2 2 2 1 2 2 2 2 2 2 2 2
## [106] 1 2 1 2 2 2 2 1 2 2 2 2 2 1 1 2 2 1 1 1 2 2 2 2 1 2 1 1 1 2 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1
## 
## Within cluster sum of squares by cluster:
## [1] 85.38078 92.25036
##  (between_SS / total_SS =  49.8 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"
t <- table(df$type, kmeans_fit2_2$cluster)
rownames(t) <- c("1 type", "2 type", "3 type")
t
##         
##           1  2
##   1 type  5 54
##   2 type 18 53
##   3 type 42  6
plot(df2[,1], df2[,2], col = kmeans_fit2_2$cluster, xlab = "Alcohol", ylab = "Malic")

Но мы-то знаем, что у нас 3 класса на самом деле.

set.seed(5)
kmeans_fit2_3 <- kmeans(df2, 3) # k = 3
kmeans_fit2_3
## K-means clustering with 3 clusters of sizes 50, 82, 46
## 
## Cluster means:
##        malic    phenols
## 1 -0.5236095 -0.7978134
## 2 -0.4857823  0.8651787
## 3  1.4351005 -0.6750866
## 
## Clustering vector:
##   [1] 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [36] 2 2 2 2 3 2 3 2 3 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 1 2 2 1 2 1
##  [71] 1 2 1 2 2 1 1 3 1 3 2 1 1 3 1 1 1 1 1 1 1 1 1 2 2 2 1 2 2 2 1 1 2 2 1
## [106] 1 1 1 2 2 2 2 1 2 2 2 1 1 3 3 2 2 3 3 3 2 2 1 1 3 1 3 1 3 1 1 3 3 3 3
## [141] 1 1 3 3 3 3 3 3 3 3 3 1 1 3 1 3 3 3 2 2 3 3 3 3 1 3 3 3 1 3 3 1 1 3 3
## [176] 3 1 3
## 
## Within cluster sum of squares by cluster:
## [1] 24.73970 41.21231 46.08206
##  (between_SS / total_SS =  68.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"
t <- table(df$type, kmeans_fit2_3$cluster)
rownames(t) <- c("1 type", "2 type", "3 type")
t
##         
##           1  2  3
##   1 type  1 53  5
##   2 type 35 27  9
##   3 type 14  2 32
plot(df2[,1], df2[,2], col = kmeans_fit2_3$cluster, xlab = "Alcohol", ylab = "Malic")

Этот график очень похож на график с истинной раскраской по группам.

Hierarchical clustering

d2 <- dist(df2, method = "euclidean")
H_fit2 <- hclust(d2)

plot(H_fit2) 
groups <- cutree(H_fit2, k=3)
rect.hclust(H_fit2, k=3, border="red") 
rect.hclust(H_fit2, k=2, border="blue")

По иерархической кластеризации визульно можно выделить 2 или 3 кластера.

Model-Based Clustering

В качестве дополнительного метода кластеризации рассмотрим метод основенный на модели.

Каждая компонента k предполагается нормально распределенной и характеризуется параметрами:

  1. \(\mu_k\) - вектор средних

  2. \(\Sigma_k\) - ковариционная матрица

  3. Каждая точка имеет вероятность принаждлежности к каждому кластеру

Параметры модели могут быть оценены с помощью EM (Expectation-Maximization). Каждый кластер центрирован в \(\mu_k\).

Геометрические свойства кластеров (shape, volume, orientation) определяются ковариционной матрицей \(\Sigma_k\)

В пакете mclust доступны следующие идентификаторы модели: EII, VII, EEI, VEI, EVI, VVI, EEE, EEV, VEV and VVV.

Первый индетификатор относится к объему (volume), второй к форме (shape), третий к ориентации (orientation). E означает равенство (equal), V - изменчивость (variable), I - идентичность.

mod_fit2 <- Mclust(df2)
summary(mod_fit2, parameters = TRUE)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust EVI (diagonal, equal volume, varying shape) model with 3 components:
## 
##  log.likelihood   n df       BIC       ICL
##       -428.3401 178 12 -918.8617 -935.8563
## 
## Clustering table:
##   1   2   3 
## 113  17  48 
## 
## Mixing probabilities:
##         1         2         3 
## 0.6188688 0.1042878 0.2768434 
## 
## Means:
##               [,1]      [,2]       [,3]
## malic   -0.6405217 1.2358592  0.9663005
## phenols  0.3938525 0.5688192 -1.0947122
## 
## Variances:
## [,,1]
##             malic   phenols
## malic   0.1101862 0.0000000
## phenols 0.0000000 0.8151158
## [,,2]
##             malic   phenols
## malic   0.3979572 0.0000000
## phenols 0.0000000 0.2256888
## [,,3]
##             malic   phenols
## malic   0.6902664 0.0000000
## phenols 0.0000000 0.1301157

Оптимальное количество клатеров - 3. Наша модель - VVI. (Диагональные ковариционные матрицы, кластеры различны по объему и форме)

head(mod_fit2$z)
##           [,1]         [,2]         [,3]
## [1,] 0.9971941 2.805784e-03 8.438666e-08
## [2,] 0.9958820 4.115451e-03 2.529526e-06
## [3,] 0.8452349 1.547637e-01 1.426477e-06
## [4,] 0.9999537 4.632382e-05 1.041407e-21
## [5,] 0.4231833 5.768129e-01 3.819913e-06
## [6,] 0.9990321 9.679014e-04 3.825326e-13

Вероятность точки принадлежать данному кластеру

t <- table(df$type,mod_fit2$classification)
rownames(t) <- c("1 type", "2 type", "3 type")
t
##         
##           1  2  3
##   1 type 51  8  0
##   2 type 56  7  8
##   3 type  6  2 40

Видим, в первый кластер попали 1 и 2 сорт вина, а почти всеь третий сорт - в 3 кластер.

plot(mod_fit2, "BIC")

Выбор параметров модели основан на BIC. У какой модели BIC больше, ту и берем. Этот график предстваляет сравнение BIC с разным набором параметров.

plot(mod_fit2, "classification")

Представление точек с плоскости признаков с раскраской по кластерам. Эта картинка совсем не похожа на ту, которую мы получили с помощью kmeans и на истинное расположение групп.

plot(mod_fit2, "density")

Линии уровня плотности.

Кластеризация по всем признакам

K-Means

Возьмем 3 кластера.

set.seed(12)
kmeans_fit <- kmeans(df_stand, 3) # k = 3
kmeans_fit
## K-means clustering with 3 clusters of sizes 64, 65, 49
## 
## Cluster means:
##      alcohol      malic        ash alcalinity   magnesium     phenols
## 1  0.7864297 -0.3417790  0.2405124 -0.6591262  0.64806806  0.83917754
## 2 -0.9145601 -0.3437694 -0.4241500  0.2102009 -0.61840118 -0.08314348
## 3  0.1860184  0.9024258  0.2485092  0.5820616 -0.02612813 -0.98577624
##   flavanoids nonflavanoids proanthocyanins      color        hue
## 1  0.8169078   -0.59546935      0.61129895  0.2684744  0.4910412
## 2  0.2130814    0.08888258     -0.03839514 -0.9614914  0.4120431
## 3 -1.3496406    0.65985043     -0.74749896  0.9247874 -1.1879477
##     dilution    proline
## 1  0.7438950  1.0776124
## 2  0.2459502 -0.8705320
## 3 -1.2978785 -0.2527064
## 
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 1
##  [71] 2 2 2 1 1 2 2 2 1 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2
## [106] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3
## [141] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [176] 3 3 3
## 
## Within cluster sum of squares by cluster:
## [1] 384.6019 585.5024 287.9127
##  (between_SS / total_SS =  45.3 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"
t <- table(df$type, kmeans_fit$cluster)
rownames(t) <- c("1 type", "2 type", "3 type")
t
##         
##           1  2  3
##   1 type 59  0  0
##   2 type  5 65  1
##   3 type  0  0 48

Хорошие результаты кластеризации.

clusplot(df_stand, kmeans_fit$cluster, main='2D representation of the Cluster solution',color=TRUE, 
         shade=TRUE,labels=2, lines=0)

Представление кластеров в плоскости первых двух главных компонент.

Hierarchical Clustering

d <- dist(df_stand, method = "euclidean")
H_fit <- hclust(d)
plot(H_fit) 
groups <- cutree(H_fit, k=3)
rect.hclust(H_fit, k=3, border="red")

Визульно можно выделить 3 кластера

Model-Based Clustering

mod_fit <- Mclust(df_stand)
summary(mod_fit)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust EEE (ellipsoidal, equal volume, shape and orientation) model with 6 components:
## 
##  log.likelihood   n  df       BIC       ICL
##       -2231.335 178 174 -5364.301 -5368.273
## 
## Clustering table:
##  1  2  3  4  5  6 
## 60  8 22 39  9 40

Модель подобрала 6 кластеров. 2 и 5 кластер совсем маленькие.

t <- table(df$type,mod_fit$classification)
rownames(t) <- c("1 type", "2 type", "3 type")
t
##         
##           1  2  3  4  5  6
##   1 type 59  0  0  0  0  0
##   2 type  1  8 22 39  1  0
##   3 type  0  0  0  0  8 40
plot(mod_fit, "BIC")    

plot(mod_fit, "classification")

Нечитабельный рисунок с раскраской по кластерам.

Теперь ограничим модель до 3х кластеров.

mod_fit <- Mclust(df_stand, G = 3)
summary(mod_fit)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 3 components:
## 
##  log.likelihood   n  df       BIC       ICL
##       -2274.056 178 158 -5366.834 -5369.729
## 
## Clustering table:
##  1  2  3 
## 57 73 48
plot(mod_fit, "classification")