本章代码
聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集。它可以把大量的观测 值归约为若干个类。这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间 相似度高。这不是一个精确的定义，从而导致了各种聚类方法的出现。
聚类分析被广泛用于生物和行为科学、市场以及医学研究中。例如，一名心理学研究员可能 基于抑郁症病人的症状和人口统计学数据对病人进行聚类，试图得出抑郁症的亚型，以期通过亚型来找到更加有针对性和有效的治疗方法，同时更好地了解这种疾病。营销研究人员根据消费者 的人口统计特征与购买行为的相似性制定客户细分战略，并基于此对其中的一个或多个子组制定 相应的营销战略。医学研究人员通过对DNA微阵列数据进行聚类分析来获得基因表达模式，从而 帮助他们理解人类的正常发育以及导致许多疾病的根本原因。
最常用的两种聚类方法是层次聚类（hierarchical agglomerative clustering）和划分聚类 （partitioning clustering）。在层次聚类中，每一个观测值自成一类，这些类每次两两合并，直到所有的类被聚成一类为止。在划分聚类中，首先指定类的个数K，然后观测值被随机分成K类，再重新形成聚合的类。
这两种方法都对应许多可供选择的聚类算法。对于层次聚类来说，最常用的算法是单联动 （single linkage）、全联动（complete linkage ）、平均联动（average linkage）、质心（centroid）和 Ward方法。对于划分聚类来说，最常用的算法是K均值（K-means）和围绕中心点的划分（PAM）。 每个聚类方法都有它的优点和缺点，我们将在本章讨论。
这一章的例子围绕食物和酒（也是我的爱好）。我们对 flexclust 包中的营养数据集 nutrient 作层次聚类，以期回答以下问题。
- 基于五种营养标准的27类鱼、禽、肉的相同点和不同点是什么？
- 是否有一种方法能把这些食物分成若干个有意义的类？
我们再用划分聚类来分析178种意大利葡萄酒样品的13种化学成分。数据在 rattle 包的 wine 数据集中。这里要解决的问题如下：
- q 这些意大利葡萄酒样品能继续分成更细的组吗？
- 如果能，有多少子组？它们的特征是什么？
事实上，样品中共有三个品种的酒（记为类）。这可以帮助我们评估聚类分析能否辨别这一 结构。
尽管聚类方法种类各异，但是它们通常遵循相似的步骤。我们在16.1节描述了这些步骤。16.2 节主要探讨层次聚类，16.3节则探讨划分聚类。最后，16.4节提出了一些相关建议。为了保证本 章的代码能正常运行，你必须确保安装了以下软件包： cluster 、 NbClust 、 flexclust 、 fMultivar 、 ggplot2 和 rattle 。第17章也将用到 rattle 包。
1 聚类分析的一般步骤
像因子分析（第14章）一样，有效的聚类分析是一个多步骤的过程，这其中每一次决策都可能影响聚类结果的质量和有效性。本节介绍一个全面的聚类分析中的11个典型步骤。
(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 包中包含了能识别多元变量的离群点的函数。一个替代的方法是使用对异常值稳健的聚类方法，围绕中心点的划分（16.3.2节）可以很好地解释这种方法。
(4) 计算距离
尽管不同的聚类算法差异很大，但是它们通常需要计算被聚类的实体之间的距离。两个观测值之间最常用的距离量度是欧几里得距离，其他可选的量度包括曼哈顿距离、兰氏距离、非对称二元距离、最大距离和闵可夫斯基距离（可使用 ?dist 查看详细信息）。在这一 章中，计算距离时默认使用欧几里得距离。计算欧几里得距离的方法见16.1.1节。
(5) 选择聚类算法
接下来选择对数据聚类的方法，层次聚类对于小样本来说很实用（如150 个观测值或更少），而且这种情况下嵌套聚类更实用。划分的方法能处理更大的数据量，但是需要事先确定聚类的个数。一旦选定了层次方法或划分方法，就必须选择一个特定的聚类算法。这里再次强调每个算法都有优点和缺点。16.2节和16.3节中介绍了最常用的几种方法。你可以尝试多种算法来看看相应结果的稳健性。
(6) 获得一种或多种聚类方法
这一步可以使用步骤(5)选择的方法。
(7) 确定类的数目
为了得到最终的聚类方案，你必须确定类的数目。对此研究者们也提出了很多相应的解决方法。常用方法是尝试不同的类数（比如2～K）并比较解的质量。在 NbClust包中的NbClust() 函数提供了30个不同的指标来帮助你进行选择（也表明这个问题有多么难解）。本章将多次使用这个包。
(8) 获得最终的聚类解决方案
一旦类的个数确定下来，就可以提取出子群，形成最终的聚类方案了。
(9) 结果可视化
可视化可以帮助你判定聚类方案的意义和用处。层次聚类的结果通常表示 为一个树状图。划分的结果通常利用可视化双变量聚类图来表示。
(10) 解读类
一旦聚类方案确定，你必须解释（或许命名）这个类。一个类中的观测值有何相似之处？不同的类之间的观测值有何不同？这一步通常通过获得类中每个变量的汇总统计来完成。对于连续数据，每一类中变量的均值和中位数会被计算出来。对于混合数据（数据中包含 分类变量），结果中将返回各类的众数或类别分布。
(11) 验证结果
验证聚类方案相当于问：“这种划分并不是因为数据集或聚类方法的某种特性，而是确实给出了一个某种程度上有实际意义的结果吗？”如果采用不同的聚类方法或不同的样本，是否会产生相同的类？ fpc 、clv 和 clValid 包包含了评估聚类解的稳定性的函数。 因为观测值之间距离的计算是聚类分析的一部分，所以我们将在下一节详细讨论。
2 计算距离
聚类分析的第一步都是度量样本单元间的距离、相异性或相似性。两个观测值之间的欧几里得距离定义为：
这里i和j代表第i和第j个观测值，p是变量的个数。
查看在 flexclust 包中的营养数据集，它包括对27种肉、鱼和禽的营养物质的测量。最初的几个观测值由下面的代码给出：
data(nutrient, package="flexclust")
head(nutrient, 4)
energy protein fat calcium iron
BEEF BRAISED 340 20 28 9 2.6
HAMBURGER 245 21 17 9 2.7
BEEF ROAST 420 15 39 7 2.0
BEEF STEAK 375 19 32 9 2.6
前两个观测值（ BEEF BRAISED 和 HAMBURGER ）之间的欧几里得距离为：
R软件中自带的 dist() 函数能用来计算矩阵或数据框中所有行（观测值）之间的距离。格式是 dist(x, method=)，这里的 x 表示输入数据，并且默认为欧几里得距离。函数默认返回一个下三角矩阵，但是 as.matrix()函数可使用标准括号符号得到距离。
对于营养数据集的数据框来说，前四行的距离为：
d <- dist(nutrient)
as.matrix(d)[1:4,1:4]
BEEF BRAISED HAMBURGER BEEF ROAST BEEF STEAK
BEEF BRAISED 0.00000 95.6400 80.93429 35.24202
HAMBURGER 95.64000 0.0000 176.49218 130.87784
BEEF ROAST 80.93429 176.4922 0.00000 45.76418
BEEF STEAK 35.24202 130.8778 45.76418 0.00000
观测值之间的距离越大，异质性越大。观测值和它自己之间的距离是0。
混合数据类型的聚类分析
欧几里得距离通常作为连续型数据的距离度量。但是如果存在其他类型的数据，则需要相异的替代措施，你可以使用 cluster 包中的 daisy() 函数来获得包含任意二元（binary）、名义（nominal）、有序（ordinal）、连续（continuous）属性组合的相异矩阵。 cluster 包中的其他函数可以使用这些异质性来进行聚类分析。例如 agnes() 函数提供了层次聚类， pam() 函数提供了围绕中心点的划分的方法。
需要注意的是，在营养数据集中，距离很大程度上由能量（ energy ）这个变量控制，这是因为该变量变化范围更大。缩放数据有利于均衡各变量的影响。在下一节中，你可以对该数据集应用层次聚类分析。
3 层次聚类分析
如前所述，在层次聚类中，起初每一个实例或观测值属于一类。聚类就是每一次把两类聚成新的一类，直到所有的类聚成单个类为止，算法如下：
(1) 定义每个观测值（行或单元）为一类；
(2) 计算每类和其他各类的距离；
(3) 把距离最短的两类合并成一类，这样类的个数就减少一个；
(4) 重复步骤(2)和步骤(3)，直到包含所有观测值的类合并成单个的类为止。
在层次聚类算法中，主要的区别是它们对类的定义不同（步骤(2)）。表16-1给出了五种最常见聚类方法的定义和其中两类之间距离的定义。
单联动聚类方法倾向于发现细长的、雪茄型的类。它也通常展示一种链式的现象，即不相似的观测值分到一类中，因为它们和它们的中间值很相像。全联动聚类倾向于发现大致相等的直径紧凑类。它对异常值很敏感。平均联动提供了以上两种方法的折中。相对来说，它不像链式，而 且对异常值没有那么敏感。它倾向于把方差小的类聚合。
Ward法倾向于把有少量观测值的类聚合到一起，并且倾向于产生与观测值个数大致相等的类。它对异常值也是敏感的。质心法是一种很受欢迎的方法，因为其中类距离的定义比较简单、 易于理解。相比其他方法，它对异常值不是很敏感。但是它可能不如平均联动法或Ward方法表现 得好。
层次聚类方法可以用 hclust()函数来实现，格式是hclust(d, method=)，其中 d是通过 dist() 函数产生的距离矩阵，并且方法包括"single" 、 "complete" 、 "average" 、 "centroid"和 "ward"。 在本节中，你可以使用平均联动聚类方法处理16.1.1节介绍的营养数据。我们的目的是基于27种食物的营养信息辨别其相似性、相异性并分组。下面的代码清单提供了实施聚类的代码。
营养数据的平均联动聚类
data(nutrient, package="flexclust")
row.names(nutrient) <- tolower(row.names(nutrient))
nutrient.scaled <- scale(nutrient)
d <- dist(nutrient.scaled) fit.average <- hclust(d, method="average") plot(fit.average, hang=-1, cex=.8, main="Average Linkage Clustering")
首先载入数据，同时将行名改为小写（因为我讨厌大写的标签）。由于变量值的变化范围很大，我们将其标准化为均值为0、方差为1。27种食物之间的距离采用欧几里得距离，应用的方法是平均联动。最后，结果用树状图来展示（参见图16-1）。 plot() 函数中的 hang 命令展示观测值的标签（让它们在挂在0下面）。
树状图应该从下往上读，它展示了这些条目如何被结合成类。每个观测值起初自成一类，然 后相距最近的两类（beef braised和smoked ham）合并。其次，pork roast和pork simmered合并， chicken canned和tuna canned合并。再次，beef braised/smoked ham这一类和pork roast/pork simmered 这一类合并（这个类目前包含四种食品）。合并继续进行下去，直到所有的观测值合并成一类。高度刻度代表了该高度类之间合并的判定值。对于平均联动来说，标准是一类中的点和其他类中 的点的距离平均值。
如果你的目的是理解基于食物营养成分的相似性和相异性，图16-1就足够了。它提供了27种 食物之间的相似性/异质性的层次分析视图。tuna canned和chicken canned是相似的，但是都和clams canned有很大的不同。但是，如果最终目标是这些食品分配到的类（希望有意义的）较少，我们 需要额外的分析来选择聚类的适当个数。
NbClust包提供了众多的指数来确定在一个聚类分析里类的最佳数目。不能保证这些指标得出的结果都一致。事实上，它们可能不一样。但是结果可用来作为选择聚类个数K值的一个参考。NbClust()函数的输入包括需要做聚类的矩阵或是数据框，使用的距离测度和聚类方法，并考虑 最小和最大聚类的个数来进行聚类。它返回每一个聚类指数，同时输出建议聚类的最佳数目。
下面的代码清单使用该方法处理营养数据的平均联动聚类。
nc <- NbClust(nutrient.scaled, distance="euclidean",
min.nc=2, max.nc=15, method="average")
table(nc$Best.n[1,])
0 1 2 3 4 5 9 10 13 14 15
2 1 4 4 2 4 1 1 2 1 4
barplot(table(nc$Best.n[1,]),
xlab="Numer of Clusters", ylab="Number of Criteria",
main="Number of Clusters Chosen by 26 Criteria")
这里，四个评判准则赞同聚类个数为2，四个判定准则赞同聚类个数为3，等等。结果在图16-2中。
你可以试着用“投票”个数最多的聚类个数（2、3、5和15）并选择其中一个使得解释最有意义。下面的代码清单展示了五类聚类的方案。
clusters <- cutree(fit.average, k=5)
table(clusters)
aggregate(nutrient, by=list(cluster=clusters), median)
aggregate(as.data.frame(nutrient.scaled), by=list(cluster=clusters),median)
pdf(file = "test.pdf")
plot(fit.average, hang=-1, cex=.8,
main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=5)
dev.off()
cutree() 函数用来把树状图分成五类➊。第一类有7个观测值，第二类有16个观测值，等等。aggregate() 函数用来获取每类的中位数➋，结果有原始度量和标准度量两种形式。最后，树状图被重新绘制， rect.hclust() 函数用来叠加五类的解决方案➌。结果展示在图16-3中。
4 划分聚类分析
在划分方法中，观测值被分为K组并根据给定的规则改组成最有粘性的类。这一节讨论两种。方法：K均值和基于中心点的划分（PAM）。
4.1 K 均值聚类
最常见的划分方法是K均值聚类分析。从概念上讲，K均值算法如下：
(1) 选择K个中心点（随机选择K行）；
(2) 把每个数据点分配到离它最近的中心点；
(3) 重新计算每类中的点到该类中心点距离的平均值（也就说，得到长度为p的均值向量，这里的p是变量的个数）；
(4) 分配每个数据到它最近的中心点；
(5) 重复步骤(3)和步骤(4)直到所有的观测值不再被分配或是达到最大的迭代次数（R把10次 作为默认迭代次数）。
这种方法的实施细节可以变化。R软件使用Hartigan & Wong（1979）提出的有效算法，这种 算法是把观测值分成k组并使得观测值到其指定的聚类中心的平方的总和为最小。也就是说，在 步骤(2)和步骤(4)中，每个观测值被分配到使下式得到最小值的那一类中：
表示第i个观测值中第j个变量的值。
表示第
个类中第
个变量的均值，其中
是变量的个数。
K均值聚类能处理比层次聚类更大的数据集。另外，观测值不会永远被分到一类中。当我们 提高整体解决方案时，聚类方案也会改动。但是均值的使用意味着所有的变量必须是连续的，并 且这个方法很有可能被异常值影响。它在非凸聚类（例如U型）情况下也会变得很差。
在R中K均值的函数格式是 kmeans(x,centers) ，这里 x 表示数值数据集（矩阵或数据框）， centers 是要提取的聚类数目。函数返回类的成员、类中心、平方和（类内平方和、类间平方和、 总平方和）和类大小。
由于K均值聚类在开始要随机选择k个中心点，在每次调用函数时可能获得不同的方案。使用 set.seed() 函数可以保证结果是可复制的。此外，聚类方法对初始中心值的选择也很敏感。 kmeans() 函数有一个 nstart 选项尝试多种初始配置并输出最好的一个。例如，加上 nstart=25 会生成25个初始配置。通常推荐使用这种方法。
不像层次聚类方法，K均值聚类要求你事先确定要提取的聚类个数。同样，NbClust包可以用来作为参考。另外，在K均值聚类中，类中总的平方值对聚类数量的曲线可能是有帮助的。可根据图中的弯曲（类似于14.2.1节描述的卵石试验弯曲）选择适当的类的数量。
图像可以用下面的代码产生：
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")
}
data 参数是用来分析的数值数据， nc 是要考虑的最大聚类个数，而 seed 是一个随机数种子。
让我们用K均值聚类来处理包含178种意大利葡萄酒中13种化学成分的数据集。该数据最初来自于UCI机器学习库（http://www.ics.uci.edu/~mlearn/MLRepository.html），但是可以通过 rattle包获得。在这个数据集里，观测值代表三种葡萄酒的品种，由第一个变量（类型）表示。你可以放弃这一变量，进行聚类分析，看看是否可以恢复已知的结构。
> data(wine, package="rattle")
> head(wine)
Type Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color Hue Dilution Proline
1 1 14.23 1.71 2.43 15.6 127 2.80 3.06 0.28 2.29 5.64 1.04 3.92 1065
2 1 13.20 1.78 2.14 11.2 100 2.65 2.76 0.26 1.28 4.38 1.05 3.40 1050
3 1 13.16 2.36 2.67 18.6 101 2.80 3.24 0.30 2.81 5.68 1.03 3.17 1185
4 1 14.37 1.95 2.50 16.8 113 3.85 3.49 0.24 2.18 7.80 0.86 3.45 1480
5 1 13.24 2.59 2.87 21.0 118 2.80 2.69 0.39 1.82 4.32 1.04 2.93 735
6 1 14.20 1.76 2.45 15.2 112 3.27 3.39 0.34 1.97 6.75 1.05 2.85 1450
标准化数据：
> df <- scale(wine[-1])
> df[1:4,1:4]
Alcohol Malic Ash Alcalinity
[1,] 1.5143408 -0.56066822 0.2313998 -1.1663032
[2,] 0.2455968 -0.49800856 -0.8256672 -2.4838405
[3,] 0.1963252 0.02117152 1.1062139 -0.2679823
[4,] 1.6867914 -0.34583508 0.4865539 -0.8069748
决定聚类的个数：
> wssplot(df)
Hit <Return> to see next plot: library(NbClust)
> set.seed(1234)
> devAskNewPage(ask=TRUE)
> nc <- NbClust(df, min.nc=2, max.nc=15, method="kmeans")
Hit <Return> to see next plot: table(nc$Best.n[1,])
*** : The Hubert index is a graphical method of determining the number of clusters.
In the plot of Hubert index, we seek a significant knee that corresponds to a
significant increase of the value of the measure i.e the significant peak in Hubert
index second differences plot.
Hit <Return> to see next plot: barplot(table(nc$Best.n[1,]),xlab="Number of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria")
*** : The D index is a graphical method of determining the number of clusters.
In the plot of D index, we seek a significant knee (the significant peak in Dindex
second differences plot) that corresponds to a significant increase of the value of
the measure.
*******************************************************************
* Among all indices:
* 2 proposed 2 as the best number of clusters
* 19 proposed 3 as the best number of clusters
* 1 proposed 14 as the best number of clusters
* 1 proposed 15 as the best number of clusters
***** Conclusion *****
* According to the majority rule, the best number of clusters is 3
*******************************************************************
进行K均值聚类分析：
> set.seed(1234)
> fit.km <- kmeans(df, 3, nstart=25)
> fit.km$size
[1] 62 65 51
> fit.km$centers
Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color
1 0.8328826 -0.3029551 0.3636801 -0.6084749 0.57596208 0.88274724 0.97506900 -0.56050853 0.57865427 0.1705823
2 -0.9234669 -0.3929331 -0.4931257 0.1701220 -0.49032869 -0.07576891 0.02075402 -0.03343924 0.05810161 -0.8993770
3 0.1644436 0.8690954 0.1863726 0.5228924 -0.07526047 -0.97657548 -1.21182921 0.72402116 -0.77751312 0.9388902
Hue Dilution Proline
1 0.4726504 0.7770551 1.1220202
2 0.4605046 0.2700025 -0.7517257
3 -1.1615122 -1.2887761 -0.4059428
> aggregate(wine[-1], by=list(cluster=fit.km$cluster), mean)
cluster Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color Hue
1 1 13.67677 1.997903 2.466290 17.46290 107.96774 2.847581 3.0032258 0.2920968 1.922097 5.453548 1.0654839
2 2 12.25092 1.897385 2.231231 20.06308 92.73846 2.247692 2.0500000 0.3576923 1.624154 2.973077 1.0627077
3 3 13.13412 3.307255 2.417647 21.24118 98.66667 1.683922 0.8188235 0.4519608 1.145882 7.234706 0.6919608
Dilution Proline
1 3.163387 1100.2258
2 2.803385 510.1692
3 1.696667 619.0588
因为变量值变化很大，所以在聚类前要将其标准化➊。下一步，使用wssplot()和 Nbclust()函数确定聚类的个数➋。图16-4表示从一类到三类变化时，组内的平方总和有一个明显的下降趋势。三类之后，下降的速度减弱，暗示着聚成三类可能对数据来说是一个很好的拟合。在图16-5中， NbClust 包中的24种指标中有14种建议使用类别数为三的聚类方案。需要注意的是并非30个指标都可以计算每个数据集。
使用 kmeans()函数得到的最终聚类中，聚类中心也被输出了➌。因为输出的聚类中心是基于标准化的数据，所以可以使用 aggregate() 函数和类的成员来得到原始矩阵中每一类的变量均值。
K均值可以很好地揭示类型变量中真正的数据结构吗？交叉列表类型（葡萄酒品种）和类成员由下面的代码表示：
> ct.km <- table(wine$Type, fit.km$cluster)
> ct.km
1 2 3
1 59 0 0
2 3 65 3
3 0 0 48
你可以用 flexclust 包中的兰德指数（Rand index）来量化类型变量和类之间的协议：
> library(flexclust)
> randIndex(ct.km)
[1] 0.897
调整的兰德指数为两种划分提供了一种衡量两个分区之间的协定，即调整后机会的量度。它的变化范围是从–1（不同意）到1 （完全同意）。葡萄酒品种类型和类的解决方案之间的协定是0.9。结果不坏，那我们应该喝点酒？
4.2 围绕中心点的划分
因为K均值聚类方法是基于均值的，所以它对异常值是敏感的。一个更稳健的方法是围绕中心点的划分（PAM）。与其用质心（变量均值向量）表示类，不如用一个最有代表性的观测值来 表示（称为中心点）。K均值聚类一般使用欧几里得距离，而PAM可以使用任意的距离来计算。 因此，PAM可以容纳混合数据类型，并且不仅限于连续变量。
PAM算法如下：
(1) 随机选择K个观测值（每个都称为中心点）；
(2) 计算观测值到各个中心的距离/相异性；
(3) 把每个观测值分配到最近的中心点；
(4) 计算每个中心点到每个观测值的距离的总和（总成本）；
(5) 选择一个该类中不是中心的点，并和中心点互换； (6) 重新把每个点分配到距它最近的中心点；
(7) 再次计算总成本；
(8) 如果总成本比步骤(4)计算的总成本少，把新的点作为中心点；
(9) 重复步骤(5)～(8)直到中心点不再改变。
在PAM方法中应用基础数学的一个例子可以在这里找到：http://en.wikipedia.org/wiki/k- medoids（我不经常引用维基百科，但这确实是一个很好的例子）。
你可以使用 cluster 包中的 pam() 函数使用基于中心点的划分方法。格式是pam(x, k, metric="euclidean", stand=FALSE)，这里的 x 表示数据矩阵或数据框，k 表示聚类的个数， metric 表示使用的相似性/相异性的度量，而 stand 是一个逻辑值，表示是否有变量应该在计算 该指标之前被标准化。图16-6中列出了使用PAM方法处理葡萄酒的数据。
pdf(file = "test.pdf")
set.seed(1234)
fit.pam <- pam(wine[-1], k=3, stand=TRUE)
fit.pam$medoids
clusplot(fit.pam, main="Bivariate Cluster Plot")
dev.off()
注意，这里得到的中心点是葡萄酒数据集中实际的观测值。在这种情况下，分别选择36、107和175个观测值来代表三类。通过从13个测定变量上得到的前两个主成分绘制每个观测的坐标来创建二元图（参见第14章）。每个类用包含其所有点的最小面积的椭圆表示。还需要注意的是，PAM在如下例子中的表现不如K均值：
> ct.pam <- table(wine$Type, fit.pam$clustering)
> 1 2 3
> 1 59 0 0
> 2 16 53 2
> 3 0 1 47
> randIndex(ct.pam)
> [1] 0.699
调整的兰德指数从（K均值的）0.9下降到了0.7
5 避免不存在的类
在结束讨论前，还要提出一点注意事项。聚类分析是一种旨在识别数据集子组的方法，并且在此方面十分擅长。事实上，它甚至能发现不存在的类。
请看以下的代码:
library(fMultivar)
set.seed(1234)
df <- rnorm2d(1000, rho=.5)
df <- as.data.frame(df)
plot(df, main="Bivariate Normal Distribution with rho=0.5")
随后，使用 wssplot()和Nbclust()函数来确定当前聚类的个数：
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")
}
pdf(file = "test1.pdf")
wssplot(df)
wssplot() 函数建议聚类的个数是3，然而 NbClust 函数返回的准则多数支持2类或3类。如
果利用PAM法进行双聚类分析：
pdf(file = "test.pdf")
library(ggplot2)
library(cluster)
fit <- pam(df, k=2)
df$clustering <- factor(fit$clustering)
p <- ggplot(data=df, aes(x=V1, y=V2, color=clustering, shape=clustering))
p <- p + geom_point() + ggtitle("Clustering of Bivariate Normal Data")
print(p)
dev.off()
你将得到如图16-10所示的双聚类图像（ ggplot() 函数是综合图像 ggplot2 包的一部分，第19章包含 ggplot2 的详细内容）。
很明显划分是人为的。实际上在这里没有真实的类。那么你怎样避免这种错误呢？虽然并非万无一失，但是我发现 NbClust 包中的立方聚类规则（Cubic Cluster Criteria，CCC）往往可以帮助我们揭示不存在的结构。代码是：
plot(nc$All.index[,4], type="o", ylab="CCC",
xlab="Number of clusters", col="blue")
结果展示在图16-11中。当CCC的值为负并且对于两类或是更多的类递减时，就是典型的单峰分布。
聚类分析（或你对它的解读）找到错误聚类的能力使得聚类分析的验证步骤很重要。如果你试图找出在某种意义上“真实的”类（而不是一个方便的划分），就要确保结果是稳健的并且是可重复的。你可以尝试不同的聚类方法，并用新的样本复制结果。如果同一类持续复原，你就可以对得出的结果更加自信。
6 小结
在这一章里，我们回顾了把观测值凝聚成子组的常见聚类方法。首先，我们回顾了常见聚类分析的一般步骤。接下来，描述了层次聚类和划分聚类的常见方法。
最后，如果寻求的不只是更方便的划分方法，我强调了需要验证所产生的类。聚类分析是一个宽泛的话题，而R有一些最全面的方案来实施现有的方法。想要了解更多，可以参考CRAN任务视图的聚类分析和有限混合模型部分（http://cran.r-project.org/web/views/Cluster.html）。
除此之外，Tan、Steinbach & Kumar（2006）写了一本关于数据挖掘技术的好书。它有一章清晰地讲解了聚类分析，可以自由下载（www-users.cs.umn.edu/~kumar/dmbook/ch8.pdf）。最后，Everitt、Landau、Leese & Stahl（2011）就这个问题写了一本实用性的课本，评价很高。聚类分析方法用于发现数据集中有凝聚力的观测值分组。在下一章中，我们将考虑已经定义好分组的情况，你的目标是用一个准确的方法把观测值分入这些组。
