【r<-高级|聚类】层次聚类与划分聚类

聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集。它可以把大量的观测值归约为若干类。这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间相似度高。这不是一个精确的定义，从而导致了各种方法的出现。

通俗地来说，聚类分析是一种将数据集中数据进行分类的一个分析过程，分类的方法有很多，它们针对数据集中不同数据特征。所以在做聚类分析的时候，根据数据集的特征选择适当的聚类方法是非常有必要的。

最常用的两种聚类方法是层次聚类和划分聚类。在层次聚类中，每个观测值自成一类，这些类每次两两合并，直到所有的类被聚成一类为止。在划分聚类中，首先指定类的个数K，然后观测值被随机分成K类，再重新形成聚合的类。

这两种方法都对应许多可供选择的聚类算法。前者，常用的有单联动、全联动、平均联动、质心以及Ward方法。对于后者，最常用的是K均值(K-means)和围绕中心点的划分(PAM)。

这一章节以flexclust包中的营养数据集nutrient作为数据进行层次聚类示范，rattle包中的意大利葡萄酒样品数据集wine进行划分聚类分析。思考以下问题：

• 基于五种营养标准的食物相同和不同之处？
• 有什么方法可以将食物分为有意义的类？
• 葡萄酒样品还能细分吗？
• 如果能，怎么分？有什么特征？

确保安装以下包：cluster,NbClust,flexclust,fMultivar,ggplot2,rattle。

聚类分析一般步骤

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

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

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

   df1 <- apply(mydata, 2, function(x){(x-mean(x))sd(x)}) # 可以使用scale()函数
   df2 <- apply(mydata, 2, function(x){x/max(x)})
   df3 <- apply(mydata, 2, function(x){(x - mean(x))/mad(x)})
   

3. 寻找异常点。许多聚类方法对异常值十分敏感。可以通过outliers包中的函数来筛选异常单变量离群点。mvoutlier包能识别多元变量的离群点的函数。另一种方法是使用对异常值稳健的聚类方法，比如划分聚类。

4. 计算距离。虽然所使用的算法差异大，但是通常都需要计算被聚类的实体之间的距离。最常用欧几里得距离，其他可选曼哈顿距离、兰式距离、非对称二元距离、最大距离和闵可夫斯基距离（?dist查看详细信息）。

5. 选择聚类算法。层次聚类对于小样本很实用，划分的方法能处理更大的数据量。

6. 获得一种或多种聚类方法。使用步骤（5）选择的方法。

7. 确定类的数目。NbClust包中的NbClust()函数提供了30个不同的指标来帮助如何选择。

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

9. 结果可视化。

10. 解读类。

11. 验证结果。采用不同的聚类方法或补贴的样本，是否会产生相同的类？fcp,clv,clValid包包含了评估聚类解的稳定性的函数。

计算距离

两个观测值之间的欧几里得距离定义为：$d_{ij}=\sqrt{\sum_{p=1}^p(x_{ip}-x_{jp})}$

R中自带的dist()函数能够用来计算矩阵或数据框中所有行之间的距离。格式是dist(x, method=)，这里的x表示输入数据，并且默认为欧几里得距离。函数默认返回一个下三角矩阵，但是as.matirx()函数可使用标准括号符号得到距离。例如营养数据集的前四行：

# data in flexclust package
data(nutrient, package="flexclust")
head(nutrient)
d <- dist(nutrient)
d
as.matrix(d)[1:4, 1:4]

> 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

混合数据类型的聚类分析

如果存在其他类型的数据，则需要相异的替代措施，可以使用cluster包中的daisy()函数来获得包含任意二元、名义、有序、连续属性组合的相异矩阵。cluster包中的其他函数可以使用这些异质性来进行聚类分析。例如agnes()函数提供了层次聚类，pam()函数提供了围绕中心点的划分的方法。

需要注意的是，在营养数据集中，能量由于变化范围比其他变量更大，因而会影响数据的平衡。缩放数据有利于均衡各变量的影响。

层次聚类分析

算法：

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

在层次聚类算法中，主要区别在于第二步骤对类的定义不同，下表列出五种

聚类方法	两类之间的距离定义
单联动	一个类中的点和另一个类中的点的最小距离
全联动	一个类中的点和另一个类中的点的最大距离
平均联动	一个类中的点和另一个类中的点的平均距离（也称为UPGMA，非加权对组平均）
质心	两类中质心（变量均值向量）之间的距离。对于单个观测值来说，质心就是变量的值
Ward法	两个类之间所有变量的方差分析的平方和

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

​ hclust(d, method=) d为dist()产生的距离矩阵

​ method包括single, complete, average, centroid, ward

例如：

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")
png("Average Link Clustering.png")
plot(fit.average, hang = -1, cex=.8, main = "Average Linkage Clustering")
dev.off()

层次聚类

hang命令显示观测值的标签。

树状图应该从下往上读，它展示了这些条目如何被结合成类。每个观测值起初自成一类，然后相聚最近的两类合并。

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,])
png("numberofcluster.png")
barplot(table(nc$Best.n[1,]),
        xlab = "Number of Clusters", ylab = "Number of Criteria",
        main = "Number of Clusters Chosen by 26 Criteria")
dev.off()

聚类个数分布

我们测测“投票”最多的聚类个数并选择其中一个使得解释有意义。

> 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 Clustering\n5 Cluster Solution")
> rect.hclust(fit.average, k=5)  # 叠加五类的解决方案

K均值聚类

当需要嵌套聚类和有意义的层次结构时，层次聚类或许特别有用。在生物科学中这种情况很常见。在某种意义上分层算法是贪婪的，一旦一个观测值被分配给一个类，它就不能在后面的过程中被重新分配。

划分聚类分析

在划分方法中，观测值被分为K组并根据给定的规则改组成最有粘性的类。这里讨论两种方法：K均值和基于中心点的划分PAM。

K均值聚类

最常见的划分方法是K均值聚类分析。算法如下：

1. 选择K个中心点（随机选择K行）;
2. 把每个数据点分配到离它最近的中心点；
3. 重新计算每类中的点到该类中心点距离的平均值;
4. 分配每个数据到它最近的中心点;
5. 重复步骤3,4直到所有观测值不再被分配或是达到最大的迭代次数（R默认10次）。

这种方法的实施细节可以变化。R软件使用Hartigan & Wong （1979）提出的有效算法，这种算法是把观测值分成K组并使得观测值到其指定的聚类中心的平方的总和为最小。也就是说，在步骤2,4中，每个观测值被分配到使下式得到最小值的那一类中：
$$
ss(k)=\sum_{i=1}^n \sum_{j=0}^p (x_{ij}-\overline{x_{kj}})^2
$$

$x_{ij}$表示第i个观测值中第j个变量的值。$\overline{x_{kj}}$表示第k类中第j个变量的均值，其中p是变量的个数。

K均值聚类能处理比层次聚类更大的数据集，另外，观测值不会永远被分到一类中。这个方法很有可能被异常值影响。

在R中K均值的函数格式是kmeans(x, centers)，这里x表示数值数据集（矩阵或数据框），centers是要提取的聚类数目。函数返回类的成员、类中心、平方和和类的大小。

kmeans()函数有一个nstart选项尝试多种初始配置并输出最好的一个。通常推荐使用这种方法。

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

让我们用K均值聚类来处理包含178种意大利葡萄酒中13种化学成分的数据集。

data(wine, package = "rattle")
head(wine)

df <- scale(wine[-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.nc[1,]),
        xlab="Number of Clusters", ylab="Number of Criteria",
        main="Number of Clusters Chosen by 26 Criteria")

set.seed(1234)
fit.km <- kmeans(df, 3, nstart=15)
fit.km$size

因为变量值变化很大，所以在聚类之前要将其标准化。下一步，使用wssplot()和Nbclust()函数确定聚类的个数。

使用kmeans()函数得到的最终聚类中，聚类中心也被输出了。

> 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
> library(flexclust)
载入需要的程辑包：grid
载入需要的程辑包：lattice
载入需要的程辑包：modeltools
载入需要的程辑包：stats4
> randIndex(ct.km)
     ARI 
0.897495 

可以用flexclust包中的兰德指数来量化类型变量和类之间的协议。

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

围绕中心点的划分

因为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是一个逻辑值，表示是否有变量应该在计算该指标之前被标准化。

下图列出了PAM方法处理葡萄酒的数据。

PAM划分聚类

注意，这里得到的中心点是葡萄酒数据集中实际的观测值。

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

> ct.pam <- table(wine$Type, fit.pam$clustering)
> randIndex(ct.pam)
      ARI 
0.6994957 

调整的兰德指数从0.9下降到了0.7。