这本书关于聚类分析部分讲得相对比较简单，要做好聚类分析，还得参考其他的资料。最有收获的应该是介绍了NbClust包，用以确定聚类数目。

聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集。它可以把大量的观测值归约为若干个类。这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间相似度高。这不是一个精确的定义，从而导致了各种聚类方法的出现。最常用的两种聚类方法是层次聚类（hierarchical agglomerative clustering）和划分聚类（partitioning clustering）。在层次聚类中，每一个观测值自成一类，这些类每次两两合并，直到所有的类被聚成一类为止。在划分聚类中，首先指定类的个数K，然后观测值被随机分成K类，再重新形成聚合的类。这两种方法都对应许多可供选择的聚类算法。对于层次聚类来说，最常用的算法是单联动（single linkage）、全联动（complete linkage ）、平均联动（average linkage） 、质心（centroid）和Ward方法。对于划分聚类来说，最常用的算法是K均值（K-means）和围绕中心点的划分（PAM）。

（1)聚类分析的典型步骤

1)选择合适的变量。需要专业知识，对结果影响较大，非常重要。第一（并且可能是最重要的）步是选择你感觉可能对识别和理解数据中不同观测值分组有重要影响的变量。

2）缩放数据。如果某些变量的变化范围很大，对结果影响也很大，需要缩放。最常用的方法是将每个变量标准化为均值为0和标准差为1的变量。其他的替代方法包括每个变量被其最大值相除或该变量减去它的平均值并除以变量的平均绝对偏差。这三种方法能用下面的代码来解释：

df1 <- apply(mydata, 2, function(x){(x-mean(x))/sd(x)})

df2 <- apply(mydata, 2, function(x){x/max(x)})

df3 <- apply(mydata, 2, function(x){(x – mean(x))/mad(x)})

可以使用scale()函数来将变量标准化到均值为0和标准差为1的变量。这和第一个代码片段（df1）等价。

3）寻找异常点。很多算法对异常点敏感。通过outliers包中的函数来筛选（和删除）异常单变量离群点。 mvoutlier包中包含了能识别多元变量的离群点的函数。或者采用对异常值稳健的聚类方法。

4）计算距离。两个观测值之间最常用的距离量度是欧几里得距离，其他可选的量度包括曼哈顿距离、兰氏距离、非对称二元距离、最大距离和闵可夫斯基距离（可使用?dist查看详细信息）。

5）选择聚类算法，获得一种或多种聚类方法。可以应用多种方法来对比结果的稳健性。

6）确定类的数目。常用方法是尝试不同的类数（比如2～K）并比较解的质量。在NbClust包中的NbClust()函数提供了30个不同的指标来帮助进行选择。

7）获得最终的聚类解决方案，结果可视化，解读类，验证结果。

（2)距离计算

R软件中自带的dist()函数能用来计算矩阵或数据框中所有行（观测值）之间的距离。格式是dist(x, method=)，这里的x表示输入数据，并且默认为欧几里得距离。函数默认返回一个下三角矩阵，但是as.matrix()函数可使用标准括号符号得到距离。观测值之间的距离越大，异质性越大。观测值和它自己之间的距离是0。

d <- dist(nutrient)

as.matrix(d)[1:4,1:4]

混合数据类型的聚类分析，欧几里得距离通常作为连续型数据的距离度量。但是如果存在其他类型的数据，则需要相异的替代措施，你可以使用cluster包中的daisy()函数来获得包含任意二元（binary） 、名义（nominal） 、有序（ordinal） 、连续（continuous）属性组合的相异矩阵。 cluster包中的其他函数可以使用这些异质性来进行聚类分析。例如agnes()函数提供了层次聚类， pam()函数提供了围绕中心点的划分的方法。

（3)层次聚类分析

起初每一个实例或观测属于一类，聚类就是每一次把两类聚成新的一类，直到所有的类聚成单个类为止。层次聚类方法有：

1）单联动， 一个类中的点和另一个类中的点的最小距离；倾向于发现细长的、雪茄型的类。它也通常展示一种链式的现象，即不相似的观测值分到一类中，因为它们和它们的中间值很相像。

2） 全联动 ，一个类中的点和另一个类中的点的最大距离；倾向于发现大致相等的直径紧凑类。它对异常值很敏感。

3）平均联动 ，一个类中的点和另一个类中的点的平均距离（也称作UPGMA，即非加权对组平均）；提供了以上两种方法的折中。相对来说，它不像链式，而且对异常值没有那么敏感。它倾向于把方差小的类聚合。

4）质心， 两类中质心（变量均值向量）之间的距离。对单个的观测值来说，质心就是变量的值；质心法是一种很受欢迎的方法，因为其中类距离的定义比较简单、易于理解。相比其他方法，它对异常值不是很敏感。但是它可能不如平均联动法或Ward方法表现得好。

5） Ward法 ，两个类之间所有变量的方差分析的平方和。Ward法倾向于把有少量观测值的类聚合到一起，并且倾向于产生与观测值个数大致相等的类。它对异常值也是敏感的。

层次聚类方法可以用hclust()函数来实现，格式是hclust(d, method=)，其中d是通过dist() 函 数 产 生 的 距 离 矩 阵 ， 并 且 方 法 包 括 "single" 、 "complete" 、 "average" 、"centroid"和"ward"。

df.scaled <- scale(df)

d <- dist(df.scaled)

fit.average <- hclust(d, method="average")

plot(fit.average, hang=-1, cex=.8, main="Average Linkage Clustering")

其中，hang命令展示观测值的标签（让它们挂在0下面）。

cutree()函数用来把树状图分成五类。

clusters <- cutree(fit.average, k=5) table(clusters)

aggregate()函数用来获取每类的中位数 ，结果有原始度量和标准度量两种形式。

aggregate(df, by=list(cluster=clusters), median)

aggregate(as.data.frame(df.scaled), by=list(cluster=clusters),median)

树状图被重新绘制， rect.hclust()函数用来叠加五类的解决方案。

plot(fit.average, hang=-1, cex=.8,main="Average Linkage Clustering\n5 Cluster Solution")

rect.hclust(fit.average, k=5)

（4)确定聚类数目

NbClust包提供了众多的指数来确定在一个聚类分析里类的最佳数目。不能保证这些指标得出的结果都一致。事实上，它们可能不一样。但是结果可用来作为选择聚类个数K值的一个参考。NbClust()函数的输入包括需要做聚类的矩阵或是数据框，使用的距离测度和聚类方法，并考虑最小和最大聚类的个数来进行聚类。它返回每一个聚类指数，同时输出建议聚类的最佳数目。

library(NbClust)

devAskNewPage(ask=TRUE)

nc <- NbClust(nutrient.scaled, distance="euclidean",min.nc=2, max.nc=15, method="average")

table(nc$Best.n[1,])

barplot(table(nc$Best.n[1,]),xlab="Numer of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria")

(5)K均值聚类

选择k个中心点，把每个数据点分配到离他最近的中心点，重新计算每类的点到该类中心点距离的平均值，分配每个数据到他最近的中心点，重复迭代（R中默认迭代10次）。

在R中K均值的函数格式是kmeans(x,centers)，这里x表示数值数据集（矩阵或数据框），centers是要提取的聚类数目。函数返回类的成员、类中心、平方和（类内平方和、类间平方和、总平方和）和类大小。由于K均值聚类在开始要随机选择k个中心点，在每次调用函数时可能获得不同的方案。使用set.seed()函数可以保证结果是可复制的。此外，聚类方法对初始中心值的选择也很敏。kmeans()函数有一个nstart选项尝试多种初始配置并输出最好的一个。例如，加上nstart=25会生成25个初始配置。通常推荐使用这种方法。不像层次聚类方法， K均值聚类要求你事先确定要提取的聚类个数。同样， NbClust包可以用来作为参考。另外，在K均值聚类中，类中总的平方值对聚类数量的曲线可能是有帮助的。

自定义wssplot绘图函数，data参数是用来分析的数值数据， nc是要考虑的最大聚类个数，而seed是一个随机数种子。

wssplot <- function(data, nc=15, seed=1234){

wss <- (nrow(data)-1)*sum(apply(data,2,var))

for (i in 2:nc){

set.seed(seed)

wss[i] <- sum(kmeans(data, centers=i)$withinss) }

plot(1:nc, wss, type="b", xlab="Number of Clusters",ylab="Within groups sum of squares") }

df <- scale(df[-1])

wssplot(df)

library(NbClust)

set.seed(1234)

devAskNewPage(ask=TRUE)

nc <- NbClust(df, min.nc=2, max.nc=15, method="kmeans")

table(nc$Best.nc[1,])

barplot(table(nc$Best.n[1,]),xlab="Number of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria")

fit.km <- kmeans(df, 3, nstart=25)

fit.km$size

fit.km$centers

aggregate(df, by=list(cluster=fit.km$cluster), mean)

因为变量值变化很大，所以在聚类前要将其标准化。下一步 使用wssplot()和Nbclust()函数确定聚类的个数 。使用kmeans()函数得到的最终聚类中，聚类中心也被输出了。因为输出的聚类中心是基于标准化的数据，所以可以使用aggregate()函数和类的成员来得到原始矩阵中每一类的变量均值。

(6)围绕中心点的划分

因为K均值聚类方法是基于均值的，所以它对异常值是敏感的。一个更稳健的方法是围绕中心点的划分（PAM）。与其用质心（变量均值向量）表示类，不如用一个最有代表性的观测值来表示（称为中心点）。 K均值聚类一般使用欧几里得距离，而PAM可以使用任意的距离来计算。因此， PAM可以容纳混合数据类型，并且不仅限于连续变量。PAM算法：1) 随机选择K个观测值（每个都称为中心点）；2) 计算观测值到各个中心的距离/相异性；3) 把每个观测值分配到最近的中心点；4) 计算每个中心点到每个观测值的距离的总和（总成本）；5) 选择一个该类中不是中心的点，并和中心点互换；6) 重新把每个点分配到距它最近的中心点；7) 再次计算总成本；8) 如果总成本比步骤(4)计算的总成本少，把新的点作为中心点；9) 重复步骤5)～8）直到中心点不再改变。

cluster包中的pam()函数使用基于中心点的划分方法。格式是pam(x, k,metric="euclidean", stand=FALSE)，这里的x表示数据矩阵或数据框， k表示聚类的个数，metric表示使用的相似性/相异性的度量，而stand是一个逻辑值，表示是否有变量应该在计算该指标之前被标准化。

library(cluster)

fit.pam <- pam(df,k=7,stand = TRUE)

fit.pam$medoids

clusplot(fit.pam,main="Cluster Plot")

（7）避免不存在的类

NbClust包中的立方聚类规则（Cubic Cluster Criteria， CCC）往往可以帮助我们揭示不存在的结构。当CCC的值为负并且对于两类或是更多的类递减时，就是典型的单峰分布。

plot(nc$All.index[,4], type="o", ylab="CCC",xlab="Number of clusters", col="blue")