聚类

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

归约

先谈一下归约（Reduction），貌似看英文比中文更容易理解这个词

举个栗子 👇

一个物理学家和一个数学家正坐在教师休息室里。突然间，休息室里的咖啡机着火了。物理学家就拿了一个垃圾桶，把里面的垃圾清空，跑到水池前，给垃圾桶灌满水，随后扑灭了火。由于这个咖啡机着过一次火了，大家都同意把垃圾桶装满水放在这个咖啡机旁边。
第二天，同样的两个人又坐在同样的休息室里，咖啡机又一次着火了。这一次，数学家站了起来，把装满水的垃圾桶拿了起来，把里面的水倒掉，又放了一些垃圾在里面，交给了物理学家。这样就把问题归约到了一个之前已经解决过的问题上。

归约也就是说，我们现在遇到了个问题，可以把它转化到一个某个已解决的问题上，而不是一定要直接解决这个问题。用数学公式来表示，即对于问题 聚类分析 - 图1 和问题 聚类分析 - 图2 ，如果存在一个可以计算的函数 聚类分析 - 图3 ，使得对于任意问题 聚类分析 - 图4 的实例都有

聚类分析 - 图5%20%3D%20B(f(x))%0A#card=math&code=A%28x%29%20%3D%20B%28f%28x%29%29%0A&id=rQ58O)

那么我们就说，问题 聚类分析 - 图6 可以被归约到问题 聚类分析 - 图7 上。

聚类方法

回到聚类的定义，由于定义的不精确，导致了各种聚类方法的出现。

聚类分析 - 图8

两种定义的原理

  • 层次聚类（hierarchical agglomerative clustering）：每一个观测值自成一类，这些类两两合并，直到所有的类被聚成一类为止。

  • 划分聚类（partitioning clustering）：首先指定类的个数K，然后观测值被随机分成K类，根据某种启发式算法迭代重置数据点所属的类别，最后达到“类内的数据点足够近，类间的数据点足够远”的效果。

聚类分析的一般步骤

  1. 选择合适的变量
    高级的聚类方法也不能弥补聚类变量选不好的问题

  2. 缩放数据
    如果在分析中选择的变量变化范围很大，那么该变量对结果的影响也是最大的。缩放数据的方法： 

  1. # 将每个变量标准化为均值为0和标准差为1的变量
  2. df1 = apply(mydata, 2, function(x){(x-mean(x))/sd(x)})
  3. df.scale = sacle(mydata)
  4. # 每个变量被其最大值相除
  5. df2 = apply(mydata, 2, function(x){x/max(x)})
  6. # 每个变量减去其平均值并除以变量的平均绝对偏差
  7. df3 = apply(mydata, 2, function(x){(x-mean(x))/mad(x)})

  1. 寻找异常点
    不同的聚类方法对异常值的敏感度不同，如果聚类方法对异常值敏感，那么则可能扭曲聚类方案。可以用 outliers包 来筛选异常单变量离群点；mvoutlier包 中包含了能识别多元变量的离群点的函数；也可以使用对异常值稳健的聚类方法。

  2. 计算距离
    各种聚类算法都需要计算被聚类的实体之间的距离，用 dist()函数 来做，最常用的是欧几里得距离，其他包括曼哈顿距离、兰氏距离、非对称二元距离、最大距离、闵可夫斯基距离。


欧几里得距离定义： 聚类分析 - 图9%5E2%7D%0A#card=math&code=d%7Bij%7D%20%3D%20%5Csqrt%7B%5Csum%7Bp%3D1%7D%5Ep%20%28x%7Bip%7D%20-%20x%7Bjp%7D%29%5E2%7D%0A&id=kciKm) 

  1. dist(x, method = ) # x表示输入数据，method为距离选择，默认为欧几里得距离

聚类分析 - 图10和聚类分析 - 图11代表第聚类分析 - 图12个和第聚类分析 - 图13个观测值，聚类分析 - 图14是变量的个数

  1. dist(x, method = ) # x表示输入数据，method为距离选择，默认为欧几里得距离

  1. 选择聚类方法
    层次聚类：小样本（< 150 个观测值）
    划分聚类：大样本，但需事先确定聚类的个数
    每种聚类方法都有优缺点，可以尝试多种算法来看看结果的稳健性。

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

  3. 确定类的数目
    常用的方法是尝试不同的类数并比较解的质量， NbClust()函数 提供了30个不同的指标进行选择 

  1. NbClust(x.scaled,         # 缩放后的输入矩阵
  2.               distance = ,    # 距离选择
  3.                 min.ce = ,       # 最小聚类个数
  4.               max.ce = ,       # 最大聚类个数
  5.               method = )     # 聚类方法

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

  2. 结果可视化

  3. 解读类

  4. 验证结果

层次聚类

聚类算法

  1. 定义每一个观测值（行）为一类
  2. 计算每类和其他各类的距离
  3. 把距离最短的两类合并成一类，用 hclust()函数 
  1. hclust(d, method = ) # d为dist()函数产生的距离矩阵，method包括'single', 'complete', 'average', 'centroid', ward'
  1. 重复步骤2和3，直到所有观测值的类合并成一个

在各种层次聚类算法中，主要的区别是它们对类的定义不同

层次聚类方法

聚类方法 两类之间的距离定义 特点
单联动 一个类中的点和另一个类中的点的最小距离 发现细长、雪茄型的类。通常展示一种链式现象，即不相似的观测值分到一类，因为它们和它们的中间值很像
全联动 一个类中的点和另一个类中的点的最大距离 发现大致相等的直径紧凑类，对异常值敏感
平均联动 一个类中的点和另一个类中的点的平均距离（UPGMA，非加权对组平均） 倾向于把方差小的类聚合，不像链式，且对异常值不很敏感
质心 两类中质心（变量均值向量）之间的距离 对异常值不很敏感，但可能表现不如平均联动或Ward
Ward法 两个类之间所有变量的方差分析的平方和 倾向于把少量观测值的类聚合到一起，并且产生与观测值个数大致相等的类。对异常值敏感

层次聚类代码

  1. library(flexclust)
  2. library(NbClust)
  4. data(nutrient, package = 'flexclust') # 输入数据
  5. row.names(nutrient) = tolower(row.names(nutrient)) # 把行名改成小写，纯粹个人爱好 hhh
  6. nutrient.scaled = scale(nutrient) # 缩放数据
  7. d = dist(nutrient.scaled) # 计算所有观测值（行）之间的距离，默认为欧几里得距离
  8. fit.average = hclust(d, method = 'average') # 层次聚类，这里选择平均联动聚类
  9. plot(fit.average, hang = -1, cex = .8, main = 'Average Linkage Clustering')
  10. nc = NbClust(nutrient.scaled, # 确定聚类分析里类的最佳数目
  11.              distance = 'euclidean', 
  12.              min.nc = 2,  # 最小聚类个数
  13.              max.nc = 15, # 最大聚类个数
  14.              method = 'average')
  15. table(nc$Best.nc[1,]) # 第一行为聚类数，第二行为支持该种聚类数的评判准则数
  16. barplot(table(nc$Best.nc[1,]),
  17.         xlab = 'Number of Clusters',
  18.         ylab = 'Number of Criteria',
  19.         main = 'Number of Clusters Chosen by 26 Criteria')
  20. clusters = cutree(fit.average, k = 5) # 将hc的结果cut成自己想要的类数
  21. aggregate(nutrient, by=list(clusters), median)          
  22. aggregate(as.data.frame(nutrient.scaled), by=list(clusters), median)
  23. plot(fit.average, hang = -1, cex = .8,
  24.      main = 'Average Linkage Clustering\n5 Cluster Solution')
  25. rect.hclust(fit.average, k=5)

划分聚类

在划分方法中，观测值被分为K组并根据给定的规则改组成最有粘性的类

（一）K均值聚类

聚类算法

  1. 选择K个中心点（随机选择K行）
  • K值的选取可以参考 NbClust() 函数
  1. nc = NbClust(df, min.nc = 2, max.nc = 15, method = 'kmeans')
  2. table(nc$Best.nc[1,])
  • 也可以用自定义函数 wssplot()
  1. wssplot <- function(data, nc=15, seed=1234){
  2.   wss <- (nrow(data)-1)*sum(apply(data,2,var))
  3.   for (i in 2:nc){
  4.     set.seed(seed)
  5.     wss[i] <- sum(kmeans(data, centers=i)$withinss)}
  6.   plot(1:nc, wss, type="b", xlab="Number of Clusters",
  7.        ylab="Within groups sum of squares")}
  8. wssplot(df)
  1. 把每个数据点分配到离它最近的中心点
    把观测值分成K组并使得观测值到其指定的聚类中心的平方和为最小，即每个观测值被分配到使下式得到最小值的那一类中 聚类分析 - 图15%20%3D%20%5Csum%7Bi%3D1%7D%5En%20%5Csum%7Bj%3D0%7D%5Ep%20(x%7Bij%7D-%5Coverline%7Bx%7D%7Bkj%7D)%5E2%0A#card=math&code=ss%28k%29%20%3D%20%5Csum%7Bi%3D1%7D%5En%20%5Csum%7Bj%3D0%7D%5Ep%20%28x%7Bij%7D-%5Coverline%7Bx%7D%7Bkj%7D%29%5E2%0A&id=brqFj)

聚类分析 - 图16 表示第聚类分析 - 图17个观测值中第聚类分析 - 图18个变量的值，聚类分析 - 图19表示第聚类分析 - 图20个类中的第聚类分析 - 图21个变量的值，其中 聚类分析 - 图22 是变量的个数

  1. 重新计算每类中的点到该类中心点距离的平均值

  2. 分配每个数据到它最近的中心点

  3. 重复步骤3和4，直到所有的观测值不再被分配或是达到最大的迭代次数

一些注意事项

  • 均值的使用意味着所有的变量必须是连续的，并且会被异常值影响；

  • 在非凸聚类的情况下可能变得很差；

  • K均值聚类在开始是随机选择k个中心点，在每次调用函数时可能获得不同的方案；

  • 为使结果可重复，要使用 set.seed() 函数进行设置。

  • 当实际分组情况已知时，可以通过 flexclust 包中的兰德指数（Rank index）来量化类型变量和类之间的协议（这一步我看来没什么必要，因为分组已知的话，还聚类干什么呢？）

    兰德指数（RI）：RI取值范围为[0,1]，值越大意味着聚类结果与真实情况越吻合

代码

  1. library(rattle)
  2. library(NbClust)
  3. library(flexclust)
  5. data(wine, package = 'rattle') # 导入数据
  6. df = scale(wine[,-1]) # 缩放数据，第1列是实际分类type
  7. set.seed(1234) # 设seed，为使结果可重复
  8. nc = NbClust(df, min.nc = 2, max.nc = 15, method = 'kmeans') # 确定聚类的个数
  9. table(nc$Best.nc[1,])                                   
  10. barplot(table(nc$Best.nc[1,]),
  11.         xlab = 'Number of Clusters',
  12.         ylab = 'Number of Criteria',
  13.         main = 'Number of Clusters Chosen by 26 Criteria')
  14. set.seed(1234)
  15. fit.km = kmeans(df, 3, nstart = 25) # 这里选择分为3类
  16. fit.km$size
  17. fit.km$centers
  18. clusters = fit.km[["cluster"]] # 聚类算法分出来的类
  19. ct.km = table(wine$Type, clusters) # 前代表行，后代表列     
  20. ct.km
  21. randIndex(ct.km) # 用兰德指数看一下聚类结果和实际分类之间的吻合程度

（二）围绕中心点的划分

中心点：最具代表性的观测值来表示类

与K均值不同之处：

  • K均值聚类基于均值，对异常值敏感；PAM基于观测值，对异常值不很敏感
  • K均值聚类一般使用欧几里得距离；PAM可以使用任意距离
  • K均值聚类仅限于连续变量；PAM可以容纳混合数据类型

聚类算法

  1. 随机选择K个观测值（每个都成为中心点）
  2. 计算观测值到各个中心的距离/相异性
  3. 把每个观测值分配到最近的中心点
  4. 计算每个中心点到每个观测值距离的总和（总成本）
  5. 选择一个改类中不是中心的点，并和中心点互换
  6. 重新把每个点分配到据它最近的中心点
  7. 再次计算总成本
  8. 如果总成本比步骤4计算的总成本少，把新的点作为中心点
  9. 重复步骤5~8，直到中心点不再改变

怎么说呢，划分聚类这两种方法我看得时候是有点抽象的，但是通过代码也可以理解其中的思想，就可意会不可言传吧，也不知道怎么解释

关于PAM函数

cluster包里有个 pam() 函数，可以用来做PAM

  1. pam(x, k, metric = 'euclidean', stand = FALSE) # x表示数据矩阵或向量，k表示聚类个数，metric表示使用的相似性/相异性的度量（欧几里得距离或曼哈顿距离，stand表示是否需要标化）

代码

  1. library(rattle)
  2. library(cluster)
  4. data(wine, package = 'rattle') # 导入数据
  5. set.seed(1234)
  6. fit.pam = pam(wine[-1], k = 3, stand = TRUE) # 这里就不用scale标化了
  7. fit.pam$medoids 
  8. fit.pam$id.med # 3个中心点
  9. clusplot(fit.pam, main = 'Bivariate Cluster Plot')
  10. ct.pam = table(wine$Type, fit.pam$clustering)
  11. randIndex(ct.pam) # pam表现不如kmeans

避免不存在的类

举个例子 👇

从相关系数为0.5的二元正态分布中抽取1000个观测值

  1. library(fMultivar)
  2. set.seed(1234)
  3. df = rnorm2d(1000, rho=.5)
  4. df = as.data.frame(df)
  5. plot(df, main = 'Bivariate Normal Distribution with rho = 0.5')

聚类分析 - 图23

可见实际上观测中并没有类

接着分别用 wssplot() 和 NbClust() 来确定聚类个数

  1. wssplot <- function(data, nc=15, seed=1234){
  2.   wss <- (nrow(data)-1)*sum(apply(data,2,var))
  3.   for (i in 2:nc){
  4.     set.seed(seed)
  5.     wss[i] <- sum(kmeans(data, centers=i)$withinss)}
  6.   plot(1:nc, wss, type="b", xlab="Number of Clusters",
  7.        ylab="Within groups sum of squares")}
  8. wssplot(df)

聚类分析 - 图24

wssplot建议聚类的个数为3

wssplot的图看下降突然减慢的那个点

  1. library(NbClust)
  2. nc = NbClust(df, min.nc = 2, max.nc = 15, method = 'kmeans')
  3. barplot(table(nc$Best.nc[1,]),
  4.         xlab = 'Number of Clusters',
  5.         ylab = 'Number of Criteria',
  6.         main = 'Number of Clusters Chosen by 26 criteria')

聚类分析 - 图25

NbClust返回的准则多数支持2类

NbClust结果进一步做pam

  1. library(ggplot2)
  2. library(cluster)
  3. fit = pam(df, k = 2)
  4. df$clustering = factor(fit$clustering)
  5. ggplot(df, aes(V1, V2, color = clustering, shape = clustering))+
  6.   geom_point()+
  7.   ggtitle('Clustering of Bivariate Normal Data')

聚类分析 - 图26

可以看到这样的两类是人为划分的，实际上并没有这样的类

可以用 NbClust包 立方聚类规则（Cubic Cluster Criteria，CCC）帮助揭示不存在的类

  1. plot(nc$All.index[,4],
  2.      type = 'o',
  3.      ylab = 'CCC',
  4.      xlab = 'Number of clusters',
  5.      col = 'blue')

聚类分析 - 图27

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