实验：使用R评估kmeans聚类的最优K
baorui
2017年8月20日
本文目的
最近这几天一直在研究如何评估Kmeans聚类算法中的最优K值。主要理论依据是《数据挖掘导论》8.5.5节中介绍的SSE和Silhouette Coefficient系数的方法评估最优K。现在记录整个实验过程，作为备忘。不过，体验过程中，由于R软件使用的还不太熟练，实现过程中有些地方可能不准确，还请大牛指点。
SSE(Sum of Squares for Error)即误差项平方和。
轮廓系数（Silhouette Coefficient），是聚类效果好坏的一种评价方式。最早由 Peter J. Rousseeuw 在 1986 提出。它结合内聚度和分离度两种因素。可以用来在相同原始数据的基础上用来评价不同算法、或者算法不同运行方式对聚类结果所产生的影响。
实验步骤概述
1.下载实验数据，点击这里。
2.取k值范围，计算出SSE，并绘制出曲线图，观察规律
3.取步骤2同样的范围，计算出Silhouette Coefficient，并绘制出曲线图，观察规律
4.根据步骤2，3观察的规律，评估出最优K值
5.验证最优聚类
1.数据简介
mdata <- read.table("http://data.galaxystatistics.com/blog_data/clustering/question_cluster.txt", header=T)
dim(mdata)## [1] 97 214head(mdata, 3)## w0 w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12 w13 w14 w15 w16 w17 w18 w19
## 1 0 0 0 0 0 0 1 0 0 0 44 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 2 0 2 0 0 0 0 0 0 0
## 5 0 58 18 0 0 3 0 0 0 0 18 2 0 0 0 1 0 0 2 1
## w20 w21 w22 w23 w24 w25 w26 w27 w28 w29 w30 w31 w32 w33 w34 w35 w36 w37
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## 3 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 2 0
## w38 w39 w40 w41 w42 w43 w44 w45 w46 w47 w48 w49 w50 w51 w52 w53 w54 w55
## 1 0 0 0 0 0 0 1 1 0 0 0 0 2 0 0 0 0 0
## 3 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 7
## 5 0 3 0 2 0 2 24 0 0 0 2 2 0 0 0 0 0 0
## w56 w57 w58 w59 w60 w61 w62 w63 w64 w65 w66 w67 w68 w69 w70 w71 w72 w73
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 3 0 4 0 0 1 0 0 0 2 0 8
## w74 w75 w76 w77 w78 w79 w80 w81 w82 w83 w84 w85 w86 w87 w88 w89 w90 w91
## 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 1 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## 5 0 0 0 0 0 3 0 0 0 0 0 0 1 0 0 4 0 1
## w92 w93 w94 w95 w96 w97 w98 w99 w100 w101 w102 w103 w104 w105 w106 w107
## 1 0 2 0 1 2 0 0 0 0 0 0 0 2 0 0 0
## 3 0 0 0 0 19 0 0 0 0 0 0 0 2 0 0 0
## 5 4 0 0 0 24 0 0 0 0 0 0 0 4 0 4 0
## w108 w109 w110 w111 w112 w113 w114 w115 w116 w117 w118 w119 w120 w121
## 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0
## 3 0 0 0 0 0 0 1 0 0 1 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 2 1
## w122 w123 w124 w125 w126 w127 w128 w129 w130 w131 w132 w133 w134 w135
## 1 0 0 0 0 2 0 0 0 0 0 0 0 0 0
## 3 0 0 2 0 0 0 0 0 0 0 0 0 0 0
## 5 1 0 0 0 0 6 0 0 0 0 2 0 0 0
## w136 w137 w138 w139 w140 w141 w142 w143 w144 w145 w146 w147 w148 w149
## 1 0 0 0 0 1 0 0 0 0 0 1 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 3 0
## 5 0 0 0 0 0 0 1 0 0 1 0 0 2 6
## w150 w151 w152 w153 w154 w155 w156 w157 w158 w159 w160 w161 w162 w163
## 1 1 1 1 1 0 1 0 0 0 1 0 0 0 0
## 3 0 0 0 0 0 0 4 0 0 0 0 0 0 0
## 5 0 21 0 0 0 0 0 4 0 0 6 0 0 6
## w164 w165 w166 w167 w168 w169 w170 w171 w172 w173 w174 w175 w176 w177
## 1 0 0 1 3 0 0 0 0 0 1 0 0 0 0
## 3 0 2 0 2 0 8 1 0 0 0 0 0 0 0
## 5 0 1 0 2 0 1 27 0 0 2 0 0 0 0
## w178 w179 w180 w181 w182 w183 w184 w185 w186 w187 w188 w189 w190 w191
## 1 0 0 0 0 0 0 0 0 0 5 0 0 0 0
## 3 0 0 0 0 17 0 0 0 0 0 0 0 1 0
## 5 0 2 0 0 0 0 0 0 2 2 0 0 0 1
## w192 w193 w194 w195 w196 w197 w198 w199 w200 w201 w202 w203 w204 w205
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 10
## 3 0 0 0 0 1 0 0 0 0 0 1 0 0 18
## 5 0 0 0 0 1 2 0 0 0 0 16 3 0 2
## w206 w207 w208 w209 w210 w211 w212 w213
## 1 0 2 0 0 0 0 0 0
## 3 1 0 0 0 19 0 0 0
## 5 0 2 2 3 4 0 1 0数据代表的是一些文本的特征矩阵，第一行是代表的词语，出于保密原因，词语被转义了。每一行代表一个文本。第二行开始的第一列是文本的ID。所以，上图中红框中的部分才是文本的特征矩阵。矩阵中每个元素的值代表所在文本当前词语的词频。
步骤2：计算SSE
主要采用R软件计算不同K值的SSE，R脚本如下：
# 开始与结果边界
begin <- 1
length <- 15
count <- 50
end <- begin + length - 1
# 结果容器
result <- c()
result[begin:end] <- 0
# 遍历计算kmeans的SSE
qc <- read.table("http://data.galaxystatistics.com/blog_data/clustering/question_cluster.txt", header=T)
for(i in begin:end) {
# 计算SSE
tmp <- c()
tmp[1:count] <- 0
for(j in 1:count) {
kcluster <- kmeans(qc, i)
tmp[j] <- kcluster$tot.withinss
}
result[i] <- mean(tmp)
}
# 绘制结果
plot(result, type="o", xlab="Number of Cluster", ylab="Sum of Squer Error");
此脚本按照K从1到15，计算不同的聚类的SSE，由于kmeans算法中的随机因数，每次结果都不一样，为了减少时间结果的偶然性，对于每个k值，都重复运行50次，求出平均的SSE，最后绘制出SSE曲线，如上所示。
步骤3：计算Silhouette Coefficient
仍然采用R脚本计算，脚本如下：
# 开始与结果边界
begin <- 2
length <- 15
count <- 50
end <- begin + length - 1
# 结果容器
result <- c()
result[begin:end] <- -1
# 遍历计算kmeans的SSE
library(cluster);
qc <- read.table("http://data.galaxystatistics.com/blog_data/clustering/question_cluster.txt", header=T)
for(i in begin:end) {
# Silhouette coefficient
tmp <- c()
tmp[1:count] <- 0
for(j in 1:count) {
kcluster <- pam(qc, i)
tmp[j] <- kcluster$silinfo$avg.width
}
result[i] <- mean(tmp)
}
# 绘制结果
plot(result, type="o", xlab="Number of Cluster", ylab="Silhouette Coefficient");
K从2到15（k=1时无法计算），重复执行50次，得到结果如上。
步骤4：估算K值
SSE图中没有什么特点，但是充Silhouette Coefficient图中可以明显看到K=8与K=9之间有一个巨大的深沟，根据Silhouette Coefficient的定义，值较大时的K较优，所以估算最优K=8为最优点。但是，书上给出的示意图是可以看到SSE中对应的地方应该改也有一个比较大的特征点，可能是kmeans的实现算法和pam具有一定的差异造成的。
步骤5：验证最优K
采用如下代码计算8个聚类的kmeans
qc <- read.table("http://data.galaxystatistics.com/blog_data/clustering/question_cluster.txt", header=T)
kc <- kmeans(qc, centers=8)
# plot(qc, col = kc$cluster)
# points(kc$centers, col = 1:2, pch = 8, cex = 2)
plot(kc$centers, col = 1:2, pch = 8, cex = 2)
使用下面的代码将实验数据从N维空间降低到2维平面，并绘制图像
qc <- read.table("http://data.galaxystatistics.com/blog_data/clustering/question_cluster.txt", header=T)
eucQcd <- dist(qc, method="euclidean")
kc <- kmeans(qc, centers=8)
mds <- cmdscale(eucQcd)
plot(mds, col=kc$cluster, main="8 clusters")
kmeans - example
# ?kmeans
# ?pam
require(graphics)
# a 2-dimensional example
x <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2), matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
colnames(x) <- c("x", "y")
(cl <- kmeans(x, 3))## K-means clustering with 3 clusters of sizes 51, 22, 27
##
## Cluster means:
## x y
## 1 0.04243926 0.08835139
## 2 1.26043172 0.98926433
## 3 0.76246993 1.09788393
##
## 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 2 2 3 2 2 3 3 3 3 2 2 2 3 2 3 2 3 2 3 3
## [71] 3 3 3 3 2 3 2 2 3 2 3 2 3 2 2 3 3 3 2 3 1 2 2 3 3 3 3 2 2 3
##
## Within cluster sum of squares by cluster:
## [1] 10.705345 1.816357 1.717204
## (between_SS / total_SS = 77.3 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"plot(x, col = cl$cluster)
points(cl$centers, col = 1:2, pch = 8, cex = 2)
cl$tot.withinss## [1] 14.23891sum(cl$withinss)## [1] 14.23891pam - example
## generate 25 objects, divided into 2 clusters.
x <- rbind(cbind(rnorm(10,0,0.5), rnorm(10,0,0.5)),cbind(rnorm(15,5,0.5), rnorm(15,5,0.5)))
pamx <- pam(x, 3)
pamx # Medoids: '7' and '25' ...## Medoids:
## ID
## [1,] 7 0.4439699 -0.08594416
## [2,] 20 5.0627011 4.96005779
## [3,] 18 6.0983947 5.82482429
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 2 2 2 2 2 3 2
## Objective function:
## build swap
## 0.4854602 0.4456745
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"summary(pamx)## Medoids:
## ID
## [1,] 7 0.4439699 -0.08594416
## [2,] 20 5.0627011 4.96005779
## [3,] 18 6.0983947 5.82482429
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 2 2 2 2 2 3 2
## Objective function:
## build swap
## 0.4854602 0.4456745
##
## Numerical information per cluster:
## size max_diss av_diss diameter separation
## [1,] 10 1.067577 0.53376320 1.436048 6.1089348
## [2,] 13 1.269843 0.43664111 1.788741 0.9142011
## [3,] 2 0.127897 0.06394851 0.127897 0.9142011
##
## Isolated clusters:
## L-clusters: character(0)
## L*-clusters: [1] 1 3
##
## Silhouette plot information:
## cluster neighbor sil_width
## 7 1 2 0.9132609
## 10 1 2 0.9129675
## 3 1 2 0.9102523
## 5 1 2 0.9053379
## 4 1 2 0.9029237
## 6 1 2 0.8968459
## 1 1 2 0.8882844
## 8 1 2 0.8837715
## 9 1 2 0.8738197
## 2 1 2 0.8579587
## 14 2 3 0.6450287
## 20 2 3 0.6408126
## 11 2 3 0.6299775
## 16 2 3 0.6146810
## 22 2 3 0.6046692
## 17 2 3 0.5932189
## 23 2 3 0.5860739
## 21 2 3 0.5850457
## 15 2 3 0.5469291
## 12 2 3 0.5418216
## 13 2 3 0.3182719
## 19 2 3 0.3043982
## 25 2 3 0.1325566
## 18 3 2 0.9120997
## 24 3 2 0.9086218
## Average silhouette width per cluster:
## [1] 0.8945422 0.5187296 0.9103608
## Average silhouette width of total data set:
## [1] 0.7003852
##
## 300 dissimilarities, summarized :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.11106 0.74407 4.01430 4.01840 7.07700 9.22750
## Metric : euclidean
## Number of objects : 25
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"plot(pamx)
pamx$silinfo$avg.width## [1] 0.7003852pamx$silinfo$clus.avg.widths## [1] 0.8945422 0.5187296 0.9103608# head(x)
# dim(x)