R语言实战 聚类分析

16．聚类分析

1.前言

聚类分析是一种数据规约技术，在与揭露一个数据集中观测值的子集，可以把大量的观测值规约为若干个类，而类即是被定义为若干个观测值组成的群组，组内相似度高于组间相似度，即是聚类.最常用的聚类时层次聚类和划分聚类。

层次聚类即是每个观测值自成一类，这些类再两两合并，直到最后只剩下一类。常用的算法有单联动算法，全联动算法，平均联动、质心和word方法。

划分聚类是制定类的K个个数，然后随机分成K类，再重新形成聚合的类。常用的算法是K均值算法和围绕中心点的划分（PAM）。

对flexclust包中的营养数据集nutrient做层次聚类。

解决的问题为：

（1）：基于五种营养标准的27类鱼，禽，肉的相同点与不同点

（2）：是否有一种方法将这些食物分为若干个有意义的类

另一个问题是用划分聚类分析葡萄酒的化学成分，用rattle包的wine数据，解决的问题有：

（1）：意大利葡萄酒样品可否分为更细的类别

（2）：若可以，有多少子组？特征分别是什么？

2.聚类分析的一般步骤

（1）：选择合适的变量

（2）：缩放数据（将每个变量标准化，变为均值为0，标准差为1的变量）

代码实现为：

df1<-apply(mydata,2,function(x){x-mean(x)}/sd(x)))

df2<-apply(mydata,2,function(x){x-mean(x)})

df3<-apply（mydata,2,function(x){(x-mean(x))/mad(x)}）

也可以使用scale()函数来讲变量标准化

（3）:寻找异常点（采用outliers包中的函数来筛选和删除单变量的离群点，也可用mvoutlier包中包含了能识别多元变量离群点的函数）

（4）：计算距离：

（5）：选择聚类算法：（层次聚类适用于小样本，划分聚类适合较大的数据量）

（6）：获得一种或多种聚类方法

（7）：确定聚类数目

（8）：得解决方案

（9）：结果可视化

（10）：解读类

（11）：验证结果

3.计算距离

采用的是欧氏距离

R语言实现

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

4.层次聚类分析

4.1算法一览：

（1）：定义每个观测值为一类

（2）：计算类间距离

（3）：合并最短距离

（4）重复2,3，直到最后只剩下一类

4.2聚类方法一览

最小距离：一个类中的点与另一个类中的点的最小距离

最大距离：一个类中的点与另一个类中的点的最大距离

平均距离：一个类中的点与另一个类中的点的平均距离

质心法：两类中质心（变量均值向量）之间的距离

Ward法：两个类间的所有变量的方差分析的平方和

较受欢迎的是质心法，层次聚类方法的代码实现是通过hclust()函数

格式为hclust(d,method=)其中d是通过dist()函数产生的距离矩阵，方法包括“single”、“complete”、“average”、“centroid”、“ward”。

采用平均联动去对营养信息聚类，R语言实现如下：

row.names(nutrient)<-tolower(row.names(nutrient))

nutrient.scaled<-scale(nutrient)

d<-dist(nutrient.scaled)

fit.average<-hclust(d,menthod="average")

plot(fit.average,hang=-1,cex=.8,main="average")

起初各成一类，后来距离最近的成为了一类，以此下去，最后将所有的观测变量聚成了一类。

倘若昂知道具体聚类分析中最佳的分类数，可用NbClust包中的众多指标来确认，具体实现如下：

nc<-NbClust(nutrient.scaled,distance="euclidean",min.nc=2,max.nc=15,method="average")

According to the majority rule, the best number ofclusters is  2

table(nc$Best.n[1,])

0  1  2 3  4  5  9 1013 14 15

2  1  4 4  2  4  1  1 2  1  4

barplot(table(nc$Best.n[1,]),xlab="number ofcluster",ylab="number of criteria",main="number of clusterchosen by 26 criteria")

clusters<-cutree(fit.average,k=5)  #将树状图分为5类

table(clusters)                                   #看出不同类的观测值数目

clusters

 1  2 3  4  5

 9 12  3 2  1

aggregate(nutrient,by=list(cluster=culsters),median) #用来获取每类的中位数

             cluster           energy   protein fat    calcium iron

1       1                    340.0    19.0 28       9 2.50

2       2                    170.0    21.0  7      13 1.35

3       3                   70.0     11.0   1     82 5.40

4       4                    137.5    16.5  7     158 1.25

5       5                    180.0    22.0  9     367 2.50

aggregate(as.data.frame(nutrient.scaled),by=list(cluster=culsters),median)#用此函数来获取每类中的中位数

 cluster     energy    protein        fat   calcium        iron

1      1  1.3101024  0.0000000 1.2897287 -0.4480464  0.08110456

2      2 -0.3696099  0.4704005 -0.5757717-0.3967868 -0.70586311

3      3 -1.3576760 -1.8816019 -1.1087718 0.4874403  2.06563171

4      4 -0.6907314 -0.5880006 -0.5757717 1.4613715 -0.77429508

5      5 -0.2708033  0.7056007-0.3981050  4.1396825  0.08110456

plot(fit.average,hang=-1,cex=.8,main="averagelinkage clustering\n5 cluster solution")

rect.hclust(fit.average,k=5) #展示了五种分类模式

5.划分聚类分析

在划分方法中，观测值被分为K组并根据给定的规则改组成最优粘性的类，即是K均值和基于中心点划分。

5.1K均值聚类

       K均值聚类算法：

（1）：选择K个中心点（随机选择K行）

（2）：把每个数据点分配到离它最近的中心点

（3）：重k新计算每类中的点到该类中心点距离的平均值

（4）：分配每个数据到他最近的中心点

（5）：重复3,4直至所有的观测值不再被分配或是达到最大的迭代次数

K均值聚类能处理比层次聚类更大的数据集，另外，观测值不会永远分到一类中，K均值函数形式在R语言中的实现：

kmeans(x,centres)

其中x表示数据集，centres是提取的聚类数目。函数返回类的成员、类中心、平方和和类的水平。

由于K均值聚类在开始随机选择k个中心店，每次调用函数可能会获得不同的方案，因此，在事先可以先行使用种子函数，set.seed()保证结果可以复制。对中心值的选择也要主要，kmeans()函数有个nstart选项尝试多种初始配置并输出最好的一个。例如nstart=25，会生成25个初始配置。且K均值聚类通常要求事先确定提取的聚类个数，可用NbClust包的函数来做参考

nc<-NbClust(x.scaled,diatance=”euclidean”,min.nc=2,max.nc=15,method=”average”)

采用k均值聚类对葡萄酒品种进行聚类

R语言实现：

devAskNewPage(ask=T)

wnc<-NbClust(df,min.nc=2,max.nc=15,method="kmeans")

According to the majority rule, the bestnumber of clusters is  3

即最终选择的聚类个数为3种。

bar.plot(wnc$Best.n[1,], xlab="number of clusters",ylab="number of criteria",main="number of cluster chosen by 26 criteria")

 

fit.km$size

[1] 62 65 51

aggregate(wine[-1],by=list(cluster=fit.km$cluster),mean)

 cluster  Alcohol    Malic     Ash Alcalinity Magnesium  PhenolsFlavanoids

1      1 13.67677 1.997903 2.466290  17.46290 107.96774 2.847581 3.0032258

2      2 12.25092 1.897385 2.231231  20.06308  92.73846 2.247692  2.0500000

3      3 13.13412 3.307255 2.417647  21.24118  98.66667 1.683922  0.8188235

 Nonflavanoids Proanthocyanins   Color       Hue Dilution   Proline

1    0.2920968        1.922097 5.4535481.0654839 3.163387 1100.2258

2    0.3576923        1.624154 2.9730771.0627077 2.803385  510.1692

3    0.4519608        1.145882 7.2347060.6919608 1.696667  619.0588

fit.km看出了具体分类情况，即每种类所包含的观测值

aggregate函数将各类的均值展现出来。

c.km<-table(wine$Type,fit.km$cluster)

c.km

    1  2  3

  159  0 0

 2  3 65  3

 3  0  0 48

randIndex(c.km)

    ARI

0.897495

    调整的改指数为两种划分提供一种衡量两个分区之间的协定，即调整后机会的亮度。变化范围由-1到1。

看类中不同交叉列表及其类成员分布。

5.2围绕中心点的划分

由于K均值聚类方法是基于均值额，所以对异常值是敏感的，而更稳健的方法是围绕中心点的划分，K均值划分通过欧式距离而通过对中心点的划分可以通过任意距离。

其算法如下：

（1）：选择K个观测值（中心点）

（2）：计算观测值到各个中心的距离/相异性

（3）：把每个观测值分配到最近的中心点

（4）：计算每个中心点到每个观测值的距离的总和

（5）：选择一个该类中不是中心的点，并和中心点互换

（6）：重新把每个点分配到距离他最近的中心点

（7）：在此计算总和

（8）：如果总和比步骤4计算的总和少，就把新的点换做中心点

（9）：重复5-8直至其不改变

对于中心点的划分方法，其R语言实现为

pam(x,k,metric=”euclidean”,stand=F)

x表示数据矩阵或者数据框，k表示聚类的个数，metric表示使用的相似性，相异性度量，stand代表的是变量是否在计算前应该被标准化。

fit.pam<-pam(wine[-1],k=3,stand=T)

fit.pam$medoids

    Alcohol Malic  Ash AlcalinityMagnesium Phenols Flavanoids Nonflavanoids Proanthocyanins

[1,]  13.48  1.81 2.41       20.5       100   2.70       2.98          0.26            1.86

[2,]  12.25  1.73 2.12       19.0        80   1.65       2.03          0.37            1.63

[3,]  13.40  3.91 2.48       23.0       102   1.80       0.75          0.43            1.41

    Color  Hue Dilution Proline

[1,]  5.1 1.04     3.47     920

[2,]  3.4 1.00     3.17     510

[3,]  7.3 0.70     1.56     750

clusplot(fit.pam,main="bivariate clusterplot")

由clusplot函数观察具体分类后，不同的类别与原观测的情况。

且这两个成分解释了55.45的可能

6.避免不存在的类

 library(fMultivar)

set.seed(1234)

df<-rnorm2d(1000,rho=.5)#从相关系数为0.5的二元正态分布抽取1000个观测值

df<-as.data.frame(df)

plot(df,main="bivariate Normaldistribution with rho=0.5")

df<0rnorm2d(1000,rho=.5)

barplot(table(nc$Best.n[1,]),xlab="numberof clusters",ylab="number of criteria",main="number ofclusters chosen by 26 criteria")

ggplot(data=df,aes(x=V1,y=V2,color=clustering,shape=clustering))+geom_point()+ggtitle("clusteringof bivariate normal data")

即是指不存在真实的类。