聚类分析

介绍

• 聚类分析是一种数据规约技术，旨在借楼一个数据集中观测值的子集。他可以把大量的观测值归约未若干类。聚类分析被广泛应用于生物和行为科学，市场以及医学研究中。例如一名心理学家可能基于抑郁症病人的症状和人口统计学数据对病人进行聚类，试图得出抑郁症的亚型，以期通过亚型来找到更加有针对性和有效的治疗方法，同时更好地理解这个疾病。医学研究人员通过对DNA微阵列数据进行聚类分析来获得基因表达模式，从而帮助他们理解人类的正常发育以及导致许多疾病的根本原因。

• 最常用的两种聚类方法是层次聚类（hierarchical agglomerative clustring)和划分聚类(partitioning clustering)。在层次聚类中，每一个观测值自成一类，这些类每次两两合并，直到所有的类被聚成一类为止。在划分聚类中，首先指定类个数K，然后观测值被随机分成K类，再重新形成聚合的类。这两种方法都对应许多可供选择的聚类算法。对于层次聚类来说，最常用的算法是单联动和全联动以及质心。对于划分聚类来说，最常用的是K均值（K-means）和围绕中心点的划分（PAM）。每一类聚类方法都有它的优点和缺点。

步骤

• 选择合适的变量
• 缩放数据

如果我们在分析中选择的变量变化范围很大，那么该变量对结果的影响也是最大的。这往往是不可取的。最常用的将数据缩放的方法是将每个变量标准化为均值为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)})

• 寻找异常点。许多聚类方法对异常值是十分敏感的，他能扭曲我们得到的聚类方案。可以通过outliers包中的函数来筛选异常单变量利群点。mvoutlier包中心包含了能识别多元变量的离群点的函数。一个替代的方案是使用对异常值稳健的聚类方法，围绕中心店的划分。
• 计算距离
• • 选择聚类算法：层次聚类对与小样本更实用，划分的方法能处理更大的数据量。
• 获得一种或者多种聚类方法
• 确定类的数目
• 获得最终的聚类方案
• 结果可视化
• 解读类
• 验证结果

计算距离

> setwd("E:\\Rwork")
> 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
> 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

层次聚类分析

如前所述，在层次聚类中，起初每一个实例或者观测值属于一类。聚类就是每一次把两类聚成新的一类，直到所有的类聚成单个类为止。算法如下：
(1) 定义每个观测值（行或单元） 为一类；

(2) 计算每类和其他各类的距离；

(3) 把距离最短的两类合并成一类，这样类的个数就减少一个；

(4) 重复步骤（2）和步骤，直到包含所有观测值的类合并成单个的类为止；

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

• 结果如下：

> setwd("E:\\Rwork")
> 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)

• 如果最终目的是这些食品分配的类较少，需要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 = "number of cluster", ylab = "number of criteria",
        main = "number of cluster chosen by 26 cruteria")
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 
******************************************************************* 
* 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 

• 显示性状之间的关系

clusters <- cutree(fit.average, k=5)
table(clusters)
plot(fit.average, hang = -1, cex = .8,
     main = "average linkage clustering \ n5 cluster solution")

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