R语言实战(16)——聚类分析

聚类分析是一种数据归约技术,旨在揭露一个数据集中观测值的子集。它可以把大量的观测值归约为若干个类。
为若干个观测值组成的群组,群组内观测值的相似度比群间相似度高。

最常用的两种聚类方法是层次聚类(hierarchical agglomerative clustering)划分聚类(partitioning clustering)

  • 在层次聚类中,每一个观测值自成一类,这些类每次两两合并,直到所有的类被聚成一类为止。最常用的算法是单联动(single linkage)、全联动(complete linkage )、平均联动(average linkage)、质心(centroid)和Ward方法。
  • 在划分聚类中,首先指定类的个数K,然后观测值被随机分成K类,再重新形成聚合的类。最常用的算法是K均值(K-means)和围绕中心点的划分(PAM)。

尽管聚类方法种类各异,但是它们通常遵循相似的步骤。

16.1 聚类分析的一般步骤

有效的聚类分析是一个多步骤的过程,这其中每一次决策都可 能影响聚类结果的质量和有效性。

1. 选择合适的变量——最重要

选择可能对识别和理解数据中不同观测值分组有重要影响的变量。

2. 缩放数据

如果在分析中选择的变量变化范围很大,那么该变量对结果的影响也是最大的。因此往往在分析之前要缩放数据。

  1. 最常用的方法是将每个变量标准化为均值为0和标准差为1的变量。
  2. 其他的替代方法包括每个变量被其最大值相除或该变量减去它的平均值并除以变量的平均绝对偏差。

可用下面的代码来解释:

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)})
PS:可以使用scale()函数来将变量标准化到均值为0和标准差为1的变量,和第一个代码片段(df1)等价。

3. 寻找异常点

许多聚类方法对于异常值是十分敏感的,它能扭曲得到的聚类方案。

  1. 可以通过outliers包中的函数来筛选(和删除)异常单变量离群点。
  2. mvoutlier包中包含了能识别多元变量的离群点的函数。

4. 计算距离

尽管不同的聚类算法差异很大,但是它们通常需要计算被聚类的实体之间的距离。两个观测值之间最常用的距离量度是欧几里得距离,其他可选的量度包括曼哈顿距离、兰氏距离、非对称二元距离、最大距离和闵可夫斯基距离。

5. 选择聚类算法

  1. 层次聚类适用于小样本(如150个观测值或更少),而且这种情况下嵌套聚类更实用。
  2. 划分聚类能处理更大的数据量,但是需要事先确定聚类的个数。

一旦选定了层次方法或划分方法,就必须选择一个特定的聚类算法。

6. 获得一种或多种聚类方法

可以使用步骤(5)选择的方法。

7. 确定类的数目

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

8. 获得最终的聚类解决方案

一旦类的个数确定下来,就可以提取出子群,形成最终的聚类方案。

9. 结果可视化

可视化可以判定聚类方案的意义和用处。

  1. 层次聚类的结果通常表示为一个树状图。
  2. 划分结果通常利用可视化双变量聚类图来表示。

10. 解读类

一旦聚类方案确定,必须解释(或许命名)这个类。通常通过获得类中每个变量的汇总统计来完成。

  1. 对于连续数据,每一类中变量的均值和中位数会被计算出来。
  2. 对于混合数据(数据中包含分类变量),结果中将返回各类的众数或类别分布。

11. 验证结果

如果采用不同的聚类方法或不同的样本,是否会产生相同的类?
fpc、clv和clValid包包含了评估聚类解的稳定性的函数。

16.2 计算距离

聚类分析的第一步都是度量样本单元间的距离、相异性或相似性。
两个观测值之间的欧几里得距离定义为:

  1. i和j代表第i和第j个观测值
  2. p是变量的个数。

例:flexclust包中的营养数据集,包括对27种肉、鱼和禽的营养物质的测量

给出前4个观测值:

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)之间的欧几里得距离为:
d = √(340-245)^2 + (20-21)^2 + (28-17)^2 + (9-9)^2 + (26-27)^2 = 95.64
dist()函数能用来计算矩阵或数据框中所有行(观测值)之间的距离。格式是:

dist(x, method=)

x表示输入数据,并且默认为欧几里得距离。
函数默认返回一个下三角矩阵,但是as.matrix()函数可使用标准括号符号得到距离。

计算前4个观测值之间的距离: 

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

观测值之间的距离越大,异质性越大。 dist()函数计算出的前两个观测值之间的距离与手算一致。

需要注意的是,在营养数据集中,距离很大程度上由能量(energy)这个变量控制,这是因为该变量变化范围更大。缩放数据有利于均衡各变量的影响。

16.3 层次聚类分析

层次聚类算法如下:

  1. 定义每个观测值(行或单元)为一类;
  2. 计算每类和其他各类的距离;
  3. 把距离最短的两类合并成一类,这样类的个数就减少一个;
  4. 重复步骤(2)和步骤(3),直到包含所有观测值的类合并成单个的类为止。

在层次聚类算法中,主要的区别是它们对类的定义不同(步骤(2))。
表16-1 常见层次聚类方法

聚类方法 两类之间的距离定义
单联动 一个类中的点和另一个类中的点的最小距离
全联动 一个类中的点和另一个类中的点的最大距离
平均联动 一个类中的点和另一个类中的点的平均距离(也称作UPGMA,即非加权对组平均)
质心 两类中质心(变量均值向量)之间的距离。对单个的观测值来说,质心就是变量的值
Ward法 两个类之间所有变量的方差分析的平方和
  • 单联动聚类方法倾向于发现细长的、雪茄型的类。它也通常展示一种链式的现象,即不相似的观测值分到一类中,因为它们和它们的中间值很相像。
  • 全联动聚类倾向于发现大致相等的直径紧凑类。它对异常值很敏感。
  • 平均联动提供了以上两种方法的折中。相对来说,它不像链式,而且对异常值没有那么敏感。它倾向于把方差小的类聚合。
  • Ward法倾向于把有少量观测值的类聚合到一起,并且倾向于产生与观测值个数大致相等的类。它对异常值也是敏感的。
  • 质心法是一种很受欢迎的方法,因为其中类距离的定义比较简单、易于理解。相比其他方法,它对异常值不是很敏感。但是它可能不如平均联动法或Ward方法表现得好。

层次聚类方法可以用hclust()函数来实现,格式是:

hclust(d, method=)

其中d是通过dist()函数产生的距离矩阵,并且方法包括 “single”、”complete”、”average”、 “centroid”和”ward”。

基于27种食物的营养信息辨别其相似性、相异性并分组:

#营养数据的平均联动聚类data(nutrient, package="flexclust"#载入数据    row.names(nutrient) <- tolower(row.names(nutrient)) #将行名改为小写    nutrient.scaled <- scale(nutrient) #将其标准化为均值为0、方差为1    d <- dist(nutrient.scaled) #采用欧几里得距离    fit.average <- hclust(d, method="average"#应用平均联动方法    plot(fit.average, hang=-1, cex=.8,    #hang命令展示观测值的标签(让它们在挂在0下面)     main="Average Linkage Clustering"#结果用树状图来展示,见图16-1

图16-1 营养数据的平均联动聚类

树状图应该从下往上读,它展示了这些条目如何被结合成类。

  1. 每个观测值起初自成一类,然后相距最近的两类(beef braised和smoked ham)合并。
  2. pork roast和pork simmered合并,chicken canned和tuna canned合并。
  3. beef braised/smoked ham这一类和pork roast/pork simmered这一类合并(这个类目前包含四种食品)。
  4. 合并继续进行下去,直到所有的观测值合并成一类。
  5. 高度刻度代表了该高度类之间合并的判定值。对于平均联动来说,标准是一类中的点和其他类中的点的距离平均值。
NbClust包提供了众多的指数来确定在一个聚类分析里类的最佳数目。结果可用来作为选择聚类个数K值的一个参考。
输入:需要做聚类的矩阵或是数据框,使用的距离测度和聚类方法,最小和最大聚类的个数。
返回:每一个聚类指数,建议聚类的最佳数目。
选择聚类的个数:
library(NbClust)    devAskNewPage(ask=TRUE)    nc <- NbClust(nutrient.scaled, distance="euclidean"              min.nc=2, max.nc=15, method="average")    按<Return>键来看下一个图: #见图16-2    *** : 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.      按<Return>键来看下一个图: #见图16-3    *** : 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:                                                    4 proposed 2 as the best number of clusters     4 proposed 3 as the best number of clusters     2 proposed 4 as the best number of clusters     4 proposed 5 as the best number of clusters     1 proposed 9 as the best number of clusters     1 proposed 10 as the best number of clusters     2 proposed 13 as the best number of clusters     1 proposed 14 as the best number of clusters     4 proposed 15 as the best number of clusters                    ***** Conclusion *****                                 * According to the majority rule, the best number of clusters is  2 
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")按<Return>键来看下一个图: #见图16-4

图16-2 Hubert指数确定聚类个数的图形方法

在Hubert指数图中,寻找一个显著的拐点来对应测量值的显著增加,即Hubert指数第二差异图中的显著峰值。

图16-3 D指数确定聚类个数的图形方法

在D指数图中,寻找一个显著的拐点(D指数第二差分图中的显著峰值),对应于测量值的显著增加。

图16-4 使用Nbclust包中提供的26个评判准则得到的推荐聚类个数

四个评判准则赞同聚类个数为2,四个判定准则赞同聚类个数为3,等等。
可以试着用“投票”个数最多的聚类个数(2、3、5和15)并选择其中一个使得解释最有意义。

获取最终的聚类方案(五类聚类的方案):

clusters <- cutree(fit.average, k=5#分配情况,把树状图分成五类    table(clusters)    clusters     1  2  3  4  5      7 16  1  2  1 
aggregate(nutrient, by=list(cluster=clusters), median)  #描述聚类,获取每类的中位数,结果以原始度量形式输出        cluster energy protein fat calcium iron    1       1  340.0      19  29       9 2.50    2       2  170.0      20   8      13 1.45    3       3  160.0      26   5      14 5.90    4       4   57.5       9   1      78 5.70    5       5  180.0      22   9     367 2.50
aggregate(as.data.frame(nutrient.scaled),  #描述聚类,获取每类的中位数,结果以标准度量形式输出                    by=list(cluster=clusters), median)        cluster     energy    protein        fat    calcium        iron    1       1  1.3101024  0.0000000  1.3785620 -0.4480464  0.08110456    2       2 -0.3696099  0.2352002 -0.4869384 -0.3967868 -0.63743114    3       3 -0.4684165  1.6464016 -0.7534384 -0.3839719  2.40779157    4       4 -1.4811842 -2.3520023 -1.1087718  0.4361807  2.27092763    5       5 -0.2708033  0.7056007 -0.3981050  4.1396825  0.08110456
plot(fit.average, hang=-1, cex=.8,  #结果绘图             main="Average Linkage Clusteringn5 Cluster Solution")    rect.hclust(fit.average, k=5#叠加五类的解决方案,见图16-5

图16-5 使用五类解决方案得到营养数据的平均联动聚类

  • sardines canned形成自己的类,因为钙比其他食物组要高得多。
  • beef heart也是单独成类,是因为富含蛋白质和铁。
  • clams类是低蛋白和高铁的。
  • 从beef roast到pork simmered的类中,所有项目都是高能量和高脂肪的。
  • 最大的类(从mackerel canned到bluefish baked)含有相对较低 的铁。

分层算法中,一旦一个观测值被分配给一个类,它就不能在后面的过程中被重新分配。另外,层次聚类难以应用到有数百甚至数千观测值的大样本中。

16.4 划分聚类分析

在划分方法中,观测值被分为K组并根据给定的规则改组成最有粘性的类。
常用算法有K均值和基于中心点的划分(PAM)。

1. K均值聚类

最常见的划分方法是K均值聚类分析。
从概念上讲,K均值算法如下:

  1. 选择K个中心点(随机选择K行);
  2. 把每个数据点分配到离它最近的中心点;
  3. 重新计算每类中的点到该类中心点距离的平均值(也就说,得到长度为p的均值向量,这里的p是变量的个数);
  4. 分配每个数据到它最近的中心点;
  5. 重复步骤(3)和步骤(4)直到所有的观测值不再被分配或是达到最大的迭代次数(R把10次作为默认迭代次数)。

在R中K均值的函数格式是:

kmeans(x, centers)
  1. x表示数值数据集(矩阵或数据框)
  2. centers是要提取的聚类数目

函数返回:类的成员、类中心、平方和(类内平方和、类间平方和、总平方和)和类大小。

  • 由于K均值聚类在开始要随机选择k个中心点,在每次调用函数时可能获得不同的方案。使用set.seed()函数可以保证结果是可复制的

  • 此外,聚类方法对初始中心值的选择也很敏感。kmeans()函数有一个nstart选项尝试多种初始配置并输出最好的一个。例如,加上nstart=25 会生成25个初始配置。通常推荐使用这种方法。

  • 另外,在K均值聚类中,类中总的平方值对聚类数量的曲线可能是有帮助的。可根据图中的弯曲选择适当的类的数量。

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是一个随机数种子。

例:包含178种意大利葡萄酒中13种化学成分的数据集
观测值代表三种葡萄酒的品种,由第一个变量(类型)表示。可以放弃这一变量,进行聚类分析,看看是否可以恢复已知的结构。
葡萄酒数据的K均值聚类:
data(wine, package="rattle")    head(wine)        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      
df <- scale(wine[-1]) #标准化数据  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
wssplot(df) #使用wssplot()决定聚类的个数    按<Return>键来看下一个图: #见图16-6  
library(NbClust) #使用Nbclust()决定聚类的个数    set.seed(1234)    devAskNewPage(ask=TRUE)
nc <- NbClust(df, min.nc=2, max.nc=15, method="kmeans")    按<Return>键来看下一个图: #见图16-7    *** : 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.      按<Return>键来看下一个图: #见图16-8    *** : 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 
table(nc$Best.n[1,])     0  1  2  3 14 15      2  1  2 19  1  1
barplot(table(nc$Best.n[1,]),                  xlab="Number of Clusters", ylab="Number of Criteria",        main="Number of Clusters Chosen by 26 Criteria")    按<Return>键来看下一个图: #见图16-9  

图16-6 画出组内的平方和和提取的聚类个数的对比

从一类到三类下降得很快(之后下降得很慢),建议选用聚类个数为三的解决方案。

图16-7 Hubert指数确定聚类个数

图16-8 D指数确定聚类个数

图16-9 使用NbClust包中的26个指标推荐的聚类个数

使用kmeans()函数得到最终聚类:

set.seed(1234)    fit.km <- kmeans(df, 3, nstart=25#进行K均值聚类分析    fit.km$size    [162 65 51    fit.km$centers              Alcohol      Malic        Ash Alcalinity   Magnesium    1  0.8328826 -0.3029551  0.3636801 -0.6084749  0.57596208    2 -0.9234669 -0.3929331 -0.4931257  0.1701220 -0.49032869    3  0.1644436  0.8690954  0.1863726  0.5228924 -0.07526047                Phenols  Flavanoids Nonflavanoids Proanthocyanins    1  0.88274724  0.97506900   -0.56050853      0.57865427    2 -0.07576891  0.02075402   -0.03343924      0.05810161    3 -0.97657548 -1.21182921    0.72402116     -0.77751312                  Color        Hue   Dilution    Proline    1  0.1705823  0.4726504  0.7770551  1.1220202    2 -0.8993770  0.4605046  0.2700025 -0.7517257    3  0.9388902 -1.1615122 -1.2887761 -0.4059428
aggregate(wine[-1], by=list(cluster=fit.km$cluster), mean) #得到原始矩阵中每一类的变量均值        cluster  Alcohol    Malic      Ash Alcalinity Magnesium    1       1 13.67677 1.997903 2.466290   17.46290 107.96774    2       2 12.25092 1.897385 2.231231   20.06308  92.73846    3       3 13.13412 3.307255 2.417647   21.24118  98.66667          Phenols Flavanoids Nonflavanoids Proanthocyanins    Color    1 2.847581  3.0032258     0.2920968        1.922097 5.453548    2 2.247692  2.0500000     0.3576923        1.624154 2.973077    3 1.683922  0.8188235     0.4519608        1.145882 7.234706                    Hue Dilution   Proline    1 1.0654839 3.163387 1100.2258    2 1.0627077 2.803385  510.1692    3 0.6919608 1.696667  619.0588

交叉列表类型(葡萄酒品种)和类成员由下面的代码表示:

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)              ARI     0.897495

调整的兰德指数为两种划分提供了一种衡量两个分区之间的协定,即调整后机会的量度。它的变化范围是从–1(不同意)到1 (完全同意)。葡萄酒品种类型和类的解决方案之间的协定是0.9。

2. 围绕中心点的划分(PAM)

因为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)
  1. x表示数据矩阵或数据框,
  2. k表示聚类的个数,
  3. metric表示使用的相似性/相异性的度量,
  4. stand是一个逻辑值,表示是否有变量应该在计算该指标之前被标准化。

对葡萄酒数据使用基于质心的划分方法:

library(cluster)    set.seed(1234)    fit.pam <- pam(wine[-1], k=3, stand=TRUE) #聚类数据的的标准化    fit.pam$medoids #输出中心点              Alcohol Malic  Ash Alcalinity Magnesium Phenols    [1,]   13.48  1.81 2.41       20.5       100    2.70[2,]   12.25  1.73 2.12       19.0        80    1.65[3,]   13.40  3.91 2.48       23.0       102    1.80     Flavanoids Nonflavanoids Proanthocyanins Color  Hue    [1,]       2.98          0.26            1.86   5.1 1.04[2,]       2.03          0.37            1.63   3.4 1.00[3,]       0.75          0.43            1.41   7.3 0.70     Dilution Proline    [1,]     3.47     920[2,]     3.17     510[3,]     1.56     750clusplot(fit.pam, main="Bivariate Cluster Plot")   #画出聚类的方案,见图16-10

图16-10 基于意大利葡萄酒数据使用PAM算法得到的三组聚类图

注意,这里得到的中心点是葡萄酒数据集中实际的观测值。
在这种情况下,分别选择36、107和175个观测值来代表三类。通过从13个测定变量上得到的前两个主成分绘制每个观测的坐标来创建二元图。每个类用包含其所有点的最小面积的椭圆表示。

还需要注意的是,PAM在如下例子中的表现不如K均值:

ct.pam <- table(wine$Type, fit.pam$clustering)    ct.pam               1  2  3      1 59  0  0      2 16 53  2      3  0  1 47    randIndex(ct.pam)                ARI     0.6994957 

调整的兰德指数从(K均值的)0.9下降到了0.7。

16.5 避免不存在的类

聚类分析是一种旨在识别数据集子组的方法,它甚至能发现不存在的类。

library(fMultivar)    set.seed(1234)    df <- rnorm2d(1000, rho=.5) #从相关系数为0.5的二元正态分布中抽取1000个观测值    df <- as.data.frame(df)    plot(df, main="Bivariate Normal Distribution with rho=0.5"

图16-11 二元正态数据(样本量n=1000),该数据集中无类

使用wssplot()和Nbclust()函数来确定当前聚类的个数:

wssplot(df) #使用wssplot()函数确定聚类个数,见图16-12    library(NbClust)    nc <- NbClust(df, min.nc=2, max.nc=15, method="kmeans")    按<Return>键来看下一个图: #见图16-13       *** : 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.      按<Return>键来看下一个图: #见图16-14    *** : 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:                                                    * 8 proposed 2 as the best number of clusters     * 4 proposed 3 as the best number of clusters     * 1 proposed 4 as the best number of clusters     * 1 proposed 5 as the best number of clusters     * 4 proposed 9 as the best number of clusters     * 1 proposed 10 as the best number of clusters     * 1 proposed 13 as the best number of clusters     * 3 proposed 15 as the best number of clusters                      ***** Conclusion *****                                 * According to the majority rule, the best number of clusters is  2  barplot(table(nc$Best.n[1,]),                  xlab="Number of Clusters", ylab="Number of Criteria"        main="Number of Clusters Chosen by 26 Criteria"#见图16-15

图16-12 二元数据的组内平方和和K均值聚类个数的对比

图16-13 Hubert指数确定聚类个数

图16-14 D指数确定聚类个数

图16-15 使用Nbclust包中的判别准则推荐的二元数据的聚类数据,推荐的类数为2或3

wssplot()函数建议聚类的个数是3,而NbClust函数返回的准则多数支持2类或3类。

利用PAM法进行双聚类分析: 

library(ggplot2)    library(cluster)   fit <- pam(df, k=2)    df$clustering <- factor(fit$clustering)    ggplot(data=df, aes(x=V1, y=V2, color=clustering, shape=clustering)) +        geom_point() + ggtitle("Clustering of Bivariate Normal Data"

图16-16 对于二元数据的PAM聚类分析,提取两类。注意类中的数据任意分割

很明显划分是人为的。实际上在这里没有真实的类。

NbClust包中的立方聚类规则(Cubic Cluster Criteria,CCC)往往可以助我们揭示不存在的结构。

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

图16-17 二元正态数据的CCC图,它正确地表明没有类存在

当CCC的值为负并且对于两类或是更多的类递减时,就是典型的单峰分布。

 

聚类分析找到错误聚类的能力使得聚类分析的验证步骤很重要。
如果试图找出在某种意义上“真实的”类(而不是一个方便的划分),就要确保结果是稳健的并且是可重复的。可以尝试不同的聚类方法,并用新的样本复制结果。如果同一类持续复原,就可以对得出的结果更加自信。

16.6 小结

在本章的学习中,我们学习了常见聚类分析的一般步骤,描述了层次聚类和划分聚类的常见方法,以及验证不存在的类。至此,我们总共16章的《R语言实战》专题系列学习分享就结束了,不知道大家对于R语言以及它的数据处理和绘图功能是否已经有了基本的了解呢?关于R语言还有更多高级的方法技能可以通过阅读原版书籍学习到,本公众号也会继续推出各种实用的R包学习分享,欢迎大家继续关注,一起学习,共同进步~

生物信息学

R语言实战(15)——时间序列

2019-12-17 21:52:28

生物信息学

R语言实战系列——汇总贴

2019-12-17 21:57:59

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索