R语言聚类分析

自己整理编写的R语言常用数据分析模型的模板,原文件为Rmd格式,直接复制粘贴过来,作为个人学习笔记保存和分享。部分参考薛毅的《统计建模与R软件》和《R语言实战》

聚类分析是一类将数据所研究对象进行分类的统计方法,这一类方法的共同特点是:事先不知道类别的个数和结构,据以进行分析的数据是对象之间的相似性或相异性的数据。将这些相似(相异)性数据看成是对象之间的“距离”远近的一种度量,将距离近的变量归为一类,不同类之间的对象距离较远。这就是聚类分析方法的共同思路。

setwd("C://Users//DELL//Documents//WD")

需要安装的包

install.packages("NbClust")
install.packages("flexclust")

1. 数据的预处理

1.1 R里的距离计算公式

在度量空间中,可自行定义距离公式。常用的定量变量间的距离公式有:绝对值距离、Euclide距离、Minkowski距离、Chebyshev距离、Mahalanobis距离、Lance和Williams距离。除此之外还有对定型变量的距离公式。

在R软件中,dist()函数给出了各种距离的计算结果,其使用格式为:

  dist(x, method = " ",diag = FALSE, upper = FALSE, p = 2)

其中x是样本构成的数据矩阵(样本按行输入)或数据框,method表示计算距离选用的方法,缺省值为Euclide距离,其他可选用的方法有:

  • “euclidean” – Euclide距离
  • “maximum” – Chebyshev距离
  • “manhattan” – 绝对值距离
  • “canberra” – Lance距离
  • “minkowski” – Minkowski距离,其中参数p是Minkowski距离的阶数,即公式中的q
  • “binary” – 定型变量的距离
    (以上距离的数学公式可参见薛毅书P467)

diag是逻辑变量,当diag = TRUE时,给出对角线上的距离;upper是逻辑变量,当upper = TRUE时,给出上三角矩阵的值(缺省值仅给出下三角矩阵的值)

1.2 数据中心化与标准化变换

在作聚类分析过程中,大多数数据往往是不能直接参与运算的,需要先将数据作中心化或标准化处理。

* 1.2.1 中心化变换*

Xij* = Xij - Xj均, i=1,2,…,p 为中心化变换,
变换后数据的均值为0,方差阵不变

* 1.2.2 标准化变换*

Xij* = (Xij - Xj均)/Sj, i=1,2,…,p 为中心化变换,
变换后数据,每个变量的样本均值均为0,标准差为1,而且标准化后的数据与变量的量纲无关

在R软件中,用scale()函数作数据的中心化或标准化,使用格式为:

 scale(x, center = TRUE, scale = TRUE)

其中x是样本构成的数据矩阵,center是逻辑变量,TRUE(缺省值)表示对数据作中心化变换,FALSE表示不作变换;scale是逻辑变量,TRUE(缺省值)表示表示对数据作标准化变换,FALSE表示不作变换。
1.2.1中公式的计算函数为:X* = scale(x,scale = FALSE);
1.2.2中公式的计算函数为:X* = scale(x)

* 1.2.3 极差标准化变换*

Xij* = (Xij - Xj均)/Rj, i=1,2,…,n,j=1,2,…,p 为极差标准化变换,
变换后的数据,每个变量的样本均值都为0,极差为1,且|Xij*|<1,在以后的分析计算中可以减少误差的产生,同时变换后的数据也是无量纲的量

在R软件中,可用sweep()函数作极差标准化变换,变换过程如下:

 center <- sweep( x, 2, apply(x, 2, mean) )
 R <- apply(x, 2, max) - apply(x, 2, min)
 x_star <- sweep(center, 2, R, "/")

其中x是样本构成的数据矩阵,第一行是将数据中心化,第二行是计算极差Rj,j=1,2,…,p,第三行是将中心化后的数据除以极差,得到数据的极差标准化数据。

其中用到的sweep()函数是对数组或矩阵运算,格式为:

 sweep(x, MARGIN, STATS, FUN = "-" , ...)

其中x是数组或矩阵,MARGIN是运算的区域,对于矩阵来讲,1表示行,2表示列。STATS是统计量,如apply(x,2,mean)表示各列的均值。FUN表示函数的运算,缺省值为减法运算。

从sweep()函数的用法可知,将上述命令第三行改为x_star <- sweep(center, 2, sd(x), “/”),得到的就是(普通)标准化变换后的结果

* 1.2.4 极差正则化变换*

Xij* = [Xij - MIN 1<=k<=n (Xkj)]/Rj, i=1,2,…,n,j=1,2,…,p 为极差正则化变换
变换后数据 0 <= Xij* <= 1,极差为1,也是无量纲的量

利用sweep()函数,可以得到数据的极差正则化变换,其变换过程如下:

 center <- sweep( x, 2, apply(x, 2, min) )
 R <- apply(x, 2, max) - apply(x, 2, min)
 x_star <- sweep(center, 2, R, "/")

1.3 相似系数

聚类分析不仅用来对样本进行分类,也可以对变量来进行分类。在对变量进行分类时,常用相似系数来度量变量之间的相似程度。

* 1.3.1 夹角余弦*

计算公式与计算向量余弦的算法一样,变量Xi的n次观测值作为一个向量,Xj的n个观测值做另一个向量,这里不作展示

在R软件中,可用scale()函数完成两向量余弦的计算,公式如下:

 y <- scale(x, center = F, scale = T)/sqrt(nrow(x)-1)
 C <-t(y) %*% y

其中x是样本构成的数据矩阵,C是计算得出的相似系数构成的矩阵

* 1.3.2 相关系数*

相关系数就是对数据作标准化处理后的夹角余弦,也就是变量Xi和变量Xj的相关系数rij,这里记作cij,公式这里不作展示

在R软件中,cij的计算更加方便,即样本的相关矩阵:C <- cor(x)

变量之间常借助于相似系数来定义距离,如令 dij^2 = 1 - cij^2,有时也用相似系数来度量变量间的相似程度。 (上述两公式见薛毅书P473)


2.层次聚类法(也叫系统聚类法)

既可处理分类变量,也可处理连续变量,但不能同时处理两种变量类型,不需要指定类别数。聚类结果间存在着嵌套,或者说层次的关系。系统聚类方法是聚类分析诸方法中用得最多的一种,其基本思想是:开始将n个样本各自作为一类,并规定样本之间的距离和类和类之间的距离,然后将距离最近的两类合并为一个新类,计算新类与其他类的距离;重复进行两个类的合并,每次减少一类,直至所有的样本合并为一类。可用于定义“距离”的统计量包括了欧氏距离(euclidean)、马氏距离(manhattan)、两项距离(binary)、明氏距离(minkowski)。还包括相关系数和夹角余弦。在计算类间距离时则有六种不同的方法,分别是最短距离法、最长距离法、类平均法、重心法、中间距离法、离差平方和法。

在R软件中,hclust()函数提供了系统聚类的计算,plot()函数可画出系统聚类的树形图(或称为谱系图)

hclust()函数使用格式为:

 hclust(d, method = "complete", members = NULL)

其中d是由“dist”构成的结构。method是系统聚类方法(缺省值是最长距离法),可选用的其他参数有:

  • “single” – 最短距离法
  • “complete” – 最长距离法
  • “median” – 中间距离法
  • “mcquitty” – Mcquitty相似法
  • “average” – 类平均法
  • “centroid” – 重心法
  • “ward“ – 离差平方和法

members的缺省值为NULL,或与d有相同变量长度的向量,具体使用方法见帮助

plot()函数画出谱系图的格式为:

 plot(x, labels = NULL, hang = 0.1, axes = TRUE, frame.plot = FALSE, ann = TRUE, 
 main = "Cluster Dendrogram", sub = NULL, xlab = NULL, ylab = "Height",...)

其中x是由hclust()函数生成的对象,hang是表明谱系图中各类所在的位置,当hang取负值时,谱系图中的类从底部画起。其他参数见plot帮助文件

下面用实例来说明系统聚类方法

2.1 例1. 城镇居民消费水平(对变量进行分类) (cluster1)

城镇居民消费水平通常用8项指标来描述:人均粮食之处(元/人)x1,人均副食支出(元/人)x2,人均烟、酒、茶支出(元/人)x3,人均其他副食支出(元/人)x4,人均衣着商品支出(元/人)x5,人均日用品支出(元/人)x6,人均燃料支出(元/人)x7,人均非商品支出(元/人)x8。这8项指标间存在着一定的相关性。为了研究城镇居民的消费结构,需将相关性强的指标归并到一起,这实际上是对指标聚类。

* 2.1.1 绘制谱系图*

#用sas中相同的数据源,输入数据,列名称是x1-x8,行名称是1-30的序号
x <- read.csv("cluster1.csv")
#相关系数矩阵
r <- cor(x)
#生成距离结构
d<-as.dist(1-r)  #method缺省值默认是Euclide距离
#用四种不同的距离方法生成系统聚类
hc1<-hclust(d, "single"); hc2<-hclust(d, "complete") #最短距离法和最长距离法
hc3<-hclust(d, "median"); hc4<-hclust(d, "mcquitty") #中间距离法和Mcquitty相似法
#绘制出以上四种方法的树形结构图,并以2x2的形式画在一张图上.图像距离边界的距离
opar <- par(mfrow = c(2, 2))
#hang是表明谱系图中各类所在的位置 当hang取负值时,谱系图中的类从底部画起  生成谱系图
plot(hc1,hang=-1); plot(hc2,hang=-1)
plot(hc3,hang=-1); plot(hc4,hang=-1)
par(opar)

R中与绘制谱系图有关的函数还有as.dendrogram(object, hang = -1, …)

其中object是是由hclust得到的对象,此时plot()函数用法为:

 plot(x, type = c("rectangle","triangle"), center = FALSE, 
 edge.root = is.leaf(x)||!is.null(attr(x,"edgetest"), 
 nodePar = NULL, edgePar = list(),
 leaflab = c("perpendicular","textlike","none"),
 dLeaf = NULL, xlab = "", ylab = "", xaxt = "n", yaxt = "s",
 horiz = FALSE, frame.plot = FALSE,...)

其中x是由dendrogram得到的对象,type表示画谱系图的类型,“rectangle”表示矩形(缺省值),“triangle”为三角形;horiz是逻辑变量,当horiz = TRUE时,表示谱系图水平放置。其他参数见帮助文件

* 2.1.2 设置谱系图形状*

#以hc1为例绘制不同参数下的谱系图
dend1<-as.dendrogram(hc1)
#生成两行一列  图像距离边界的距离
opar <- par(mfrow = c(2, 2),mar = c(4,3,1,2))
#通过设置参数分别画出四种不同的图
plot(dend1)  #原始图像
plot(dend1, nodePar=list(pch = c(1,NA), cex=0.8, lab.cex = 0.8),
     type = "t", center=TRUE)  #斜线
plot(dend1, edgePar=list(col = 1:2, lty = 2:3), dLeaf=1, edge.root = TRUE)  #虚线
plot(dend1, nodePar=list(pch = 2:1,cex=.4*2:1, col = 2:3), horiz=TRUE)      #纵向
par(opar)

2.2 例2. 女中学生体型指标分类(对变量进行分类)(一个小函数让图形好看一点)

对305名女中学生测量8个体型指标,给出相应的相关矩阵。将相关系数看成相似系数,定义距离为:dij = 1 - rij,用最长距离法作系统聚类分析

# 输入相关矩阵. 
x<- c(1.000, 0.846, 0.805, 0.859, 0.473, 0.398, 0.301, 0.382,
      0.846, 1.000, 0.881, 0.826, 0.376, 0.326, 0.277, 0.277, 
      0.805, 0.881, 1.000, 0.801, 0.380, 0.319, 0.237, 0.345, 
      0.859, 0.826, 0.801, 1.000, 0.436, 0.329, 0.327, 0.365, 
      0.473, 0.376, 0.380, 0.436, 1.000, 0.762, 0.730, 0.629, 
      0.398, 0.326, 0.319, 0.329, 0.762, 1.000, 0.583, 0.577, 
      0.301, 0.277, 0.237, 0.327, 0.730, 0.583, 1.000, 0.539, 
      0.382, 0.415, 0.345, 0.365, 0.629, 0.577, 0.539, 1.000)
names<-c("身高 x1", "手臂长 x2", "上肢长 x3", "下肢长 x4", "体重 x5", 
         "颈围 x6", "胸围 x7", "胸宽 x8")
r<-matrix(x, nrow=8, dimnames=list(names, names))
# 作系统聚类分析, 
# as.dist()的作用是将普通矩阵转化为聚类分析用的距离结构. 
d<-as.dist(1-r); hc<-hclust(d); dend<-as.dendrogram(hc)
# 写一段小程序, 其目的是在绘图命令中调用它, 使谱系图更好看.
nP<-list(col=3:2, cex=c(2.0, 0.75), pch= 21:22, 
         bg= c("light blue", "pink"),
         lab.cex = 1.0, lab.col = "tomato")
addE <- function(n){
   if(!is.leaf(n)) {
       attr(n, "edgePar") <- list(p.col="plum")
       attr(n, "edgetext") <- paste(attr(n,"members"),"members")
   }
   n
}
# 画出谱系图.
op<-par(mfrow=c(1,1), mar=c(4,3,0.5,0))
de <- dendrapply(dend, addE); plot(de, nodePar= nP)
par(op)

从图中容易看出,变量x2(手臂长)与x3(上肢长)最先合并成一类,解析来是变量x1(身高)与x4(下肢长)合并成一类.再合并就是将新得到的两类合并成一类(可称为“长”类).后面要合并的是x5(体重)与x3(颈围).再合并就是将x7(胸围)加到新类中.再加就是x8(胸宽).最后合并为一类.


2.3 类的个数的确定

在聚类问题中如何确定合适的分类个数,是一个非常困难的问题,至今没有特别令人满意的答案。

* 2.3.1 rect.hclust()函数确定分类个数*

在R软件中,与确定类的个数有关的函数是rect.hclust()函数,它的本质是由给定类的个数或给定阈值(类与类间的距离)来确定分类个数

 rect.hclust(tree, k = NULL, which = NULL, x = NULL, h = NULL, border = 2, cluster = NULL)

其中tree是由hclust生成的结构,k是类的个数,h是谱系图中的阈值,要求分成的各类的距离大于h.border是数或向量,表示矩形框的颜色

对2.2中例2的八个体型指标的聚类分析中,将变量分为三类,即k=3,其程序和计算结果如下:

plclust(hc,hang = -1)   #与plot类似用法
re <- rect.hclust(hc,k=3)   #加入k=3的类

根据谱系图确定分类个数的准则:
A.各类重心的距离必须很大
B.确定的类中,各类所包含的元素都不要太多
C.类的个数必须符合实用目的
D.若采用集中不同的聚类方法处理,则在各自的聚类途中因发现相同的类。

得到身高(x1),手臂长(x2),上肢长(x3),下肢长(x4)分为第一类;胸宽(x8)为第二类;体重(x5),颈围(x6),胸围(x7)分为第三类。

上边程序中,plclust()函数是另一种绘谱系图的函数,与plot()函数所画图形略有差别,使用格式如下:

 plclust(tree, hang = 0.1, unit = FALSE, level = FALSE, hmin = 0,
         square = TRUE, labels = NULL, plot. = TRUE,
         axes = TRUE, frame.plot = FALSE, ann = TRUE,
         main = "", sub = NULL, xlab = NULL, ylab = "Height")

其中tree是由hclust()函数生成的对象,其他参数设置与plot()函数中的相同

* 2.3.2 例3.NbClust包确定在一个聚类分析里类的最佳数目*(用flexclust包里自带数据作为例子,clust1是对变量聚类,数据是相关系数矩阵,无法使用)

par(ask=TRUE)
opar <- par(no.readonly=FALSE)
#输入flexclust包自带数据nutrient,是27个不同种类的肉的5个成分含量
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, hang=-1, cex=.8, main="Average Linkage Clustering")
#用NbClust()函数判断生成类的个数
library(NbClust)
nc <- NbClust(nutrient.scaled, distance="euclidean", 
              min.nc=2, max.nc=15, method="average")
#设置图形
par(opar)
table(nc$Best.n[1,])
#标准图
barplot(table(nc$Best.n[1,]), xlab="Numer of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria") 

函数给出最佳类别数量为2,选择2类

#使用前先清除历史图片,或者将第二个代码块运行两次
#每个种类的个数
clusters <- cutree(fit.average, k=2) 
table(clusters)
#每个种类的中位数
aggregate(nutrient, by=list(cluster=clusters), median) 
#标准化后的中位数
aggregate(as.data.frame(nutrient.scaled), by=list(cluster=clusters),median)
#画图,k选择2
plot(fit.average, hang=-1, cex=.8,main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=2)

根据实际情况,选择5类

clusters <- cutree(fit.average, k=5) 
table(clusters)
aggregate(nutrient, by=list(cluster=clusters), median) 
aggregate(as.data.frame(nutrient.scaled), by=list(cluster=clusters),median)
plot(fit.average, hang=-1, cex=.8,  main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=5)   #k选择5

2.4 例4. 一个具体的实例来总结前面介绍的聚类分析的方法(cluster2、对样本分类)

列出1999年全国31个省市自治区,的城镇居民家庭平均每人全年消费性支出的八个主要指标(变量)数据,这八个变量是x1食品,x2衣着,x3家庭设备用品及服务,x4医疗保健,x5交通与通讯,x6娱乐教育文化服务,x7居住,x8杂项商品和服务

例. 某研究者收集了24种菌株,其中17~22号为已知的标准菌株,它们分别取自牛、羊、犬、猪、鼠、绵羊,其他的为xx菌株。测得各菌株的16种脂肪酸百分含量,试作样品聚类分析,以便了解哪些为止菌株与已知的标准菌株在全部指标上最为接近。

#用sas中相同的数据源,输入数据,列名称是x1-x8,行名称是1-30的序号
y <- read.csv("cluster2.csv")
#生成距离结构
bacterium <- dist(scale(y))  #method缺省值默认是Euclide距离,scale函数先将y作标准化变换
#用四种不同的距离方法生成系统聚类
hc1<-hclust(bacterium, "complete")  #最长距离法
hc2<-hclust(bacterium, "average")   #类平均法
hc3<-hclust(bacterium, "centroid")  #重心法
hc4<-hclust(bacterium, "ward.D2")      #离差平方和法
#以上聚类模型画图

* 2.4.1 rect.hclust()函数,根据谱系图,结合实际情况选定类别个数*

k <- 3   #根据谱系图,结合实际情况,自行选定类别个数
#生成两行一列  图像距离边界的距离
opar<-par(mfrow=c(2,1), mar=c(2,4,2,2))
#hang是表明谱系图中各类所在的位置 当hang取负值时,谱系图中的类从底部画起  生成谱系图
plot(hc1,hang=-1);re1<-rect.hclust(hc1,k,border="red")  #将分类结果分成5类 用红色矩形笔迹标记
plot(hc2,hang=-1);re2<-rect.hclust(hc2,k,border="red")
#在活动设备中返回所有图形参数和他们的值
par(opar)
opar<-par(mfrow=c(2,1), mar=c(2,4,2,2))
plot(hc3,hang=-1);re3<-rect.hclust(hc3,k,border="red")
plot(hc4,hang=-1);re4<-rect.hclust(hc4,k,border="red")
par(opar)

由结果的谱系图看,k=3时的最长距离法分类较为合理,分类结果是:
第一类:3、4、5、7、8、12、14、15、16、17、18、22、23、24、26、27、28、29、30; 第二类:1、2、6、10、11、25; 第三类:9、13、19、20、21

* 2.4.2 NbClust包确定类别个数*

#library(NbClust)
#默认index参数为all,会报错,出现极小值需要设置tolerance。但是这个函数过于复杂,不作手动修改
#nc <- NbClust(scale(y),min.nc=2,max.nc=8,method = "average")
#index=“hubert”,通过图形判断类别个数
nc <- NbClust(scale(y),min.nc=2,max.nc=8,method = "average",index = "hubert")

没有查找到Hubert统计值的具体含义,猜想左边图k=3时极小值,右图k=3到k=3+1有明显下降趋势,所以设置k=3.


3.动态聚类法(这里介绍的是用K-均值方法的kmeans()函数)例5. K-均值法作cluster2的聚类分析

 系统聚类法以此形成类之后就不能改变,这就要求一次分类分得比较准确,对分类的方法要求较高,相应的计算量也会增大。尤其是对于大样本问题,需要很长的计算时间,还会占据足够大的计算机内存。基于这种情况,产生了动态聚类,即动态聚类法。

 动态聚类法又称逐步聚类法,其基本思想是,开始先粗略的分一下类,然后按照某种最优原则修改不合理的分类,直至类分得比较合理为止,形成最终的分类结果。这种方法具有计算量较小,占计算机内存较少和方法简单的优点,适用于大样本的Q型聚类分析。

 这里介绍用于动态聚类的R函数kmeans()函数,它采用的是K-均值方法。针对连续变量,也可处理有序分类变量,运算很快,但需要指定类别数。K-均值聚类法不会自动对数据进行标准化处理,需要先自己手动进行标准化分析

 kmeans(x, centers, iter.max = 10, nstart = 1, algorithm = c("Hartigan-Wong","Lloyd","Forgy","MacQueen"))

其中x是由数据构成的矩阵或数据框,centers是聚类的个数或者是初始类的中心,iter.max为最大迭代次数(缺省值为10),nstart随机集合的个数(当centers为聚类的个数时),algorithm为动态聚类的算法(缺省值为Hartigan-Wong方法)

3.1 令类的个数k=5作聚类分析

与例4.一样,为消除数据数量级的影响,先对数据作标准化处理,然后再用kmeans()函数作动态聚类,为与前面的方法作比较,类的个数选为5

y <- read.csv("cluster2.csv")
#算法选择Hartigan-Wong方法,即缺省状态
km <- kmeans(scale(y),5,nstart = 20);km
#便于看出聚类后的分类情况,sort()函数对分类情况排序
sort(km$cluster)
library(fpc)
plotcluster(scale(y), km$cluster)    #生成聚类图

“cluster”是一个整数向量,用于表示记录所属的聚类
“centers”是一个矩阵,表示每聚类中各个变量的中心点
“totss”表示所生成聚类的总体距离平方和
“withinss”表示各个聚类组内的距离平方和
“tot.withinss”表示聚类组内的距离平方和总量
“betweenss”表示聚类组间的聚类平方和总量
“size”表示每个聚类组中成员的数量

(between_SS / total_SS = 80.9 %),组间的距离平方和占了整体距离平方和的的80.9 %,也就是说各个聚类间的距离做到了最大

结合实际情况17-22为标准类,发现k=5并不符合实际情况。

3.2 类别个数k的选择

在SAS中PROC CLUSTER会在输出聚类结果的同时还是输出一大坨统计量或者伪统计量,主要是这四个:
RSQ、SPRSQ、PSF以及PST2,分别指的是R^2、半偏R^2(伪R^2)、伪F和伪t^2 。

先简单介绍下这四个量,如果把样本分为k个类,R2统计量的思想其实和回归分析中的可决系数系数,即类间偏差平方和的总和与总离差平方总和的比值,如果值越大就说明k个类越能够区分开,即聚类效果越好。但需要注意的是,跟回归中类似,R2的值总是随着k的减少而减小,也就是说单纯的根据R2来判断是没有意义,因此我们更关注的是他的变化情况,如果从k到k-1有非常显著的下降,那么可以认为分为k类是比较合适。基于此有人就又提出了半偏R2统计量,其实就是相邻两个R2的差值。F的值越大分类效果越显著,但唯一的缺点是他只是看起来像F统计量,但实际上却不服从F分布。当然伪t2同样如此。
数学公式参考网页 https://www.52ml.net/8663.html/comment-page-1#comment-185

下面提供上述统计量的计算和图像

* 3.2.1 先看within-cluster sum of squares组内平方和的趋势图*

y <- read.csv("cluster2.csv")
mydata <- y
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
  for (i in 2:13) wss[i] <- sum(kmeans(mydata,centers=i)$withinss)
###这里的wss(within-cluster sum of squares)是组内平方和
plot(1:13, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")

* 3.2.2 R方和伪R方*

#先对k=2的时候计算,再构造循环,计算k=1-7的情况,并画图
#分类个数i为2时,R方的公式为:
R2 <- kmeans(mydata,centers=2)$betweenss / kmeans(mydata,centers=2)$tot.withinss
   for (i in 2:7)R2[i] <- kmeans(mydata,centers=i)$betweenss / kmeans(mydata,centers=i)$tot.withinss
fakeR2 <- R2[2]-R2[1]
   for (i in 2:7)fakeR2[i] <- R2[i]-R2[i-1]
plot(R2,type="b",xlab="Number of Clusters", ylab="R Square")
plot(fakeR2,type="b",xlab="Number of Clusters", ylab="Fake R Square")

* 3.2.3 伪F和伪t^2*

#先对k=2的时候计算,再构造循环,计算k=1-7的情况,并画图
fakeF <- 1 -  kmeans(mydata,centers=2)$betweenss*(nrow(mydata)-2) /
             (kmeans(mydata,centers=2)$tot.withinss-kmeans(mydata,centers=2)$betweenss)*(2-1)
   for (i in 2:7)fakeF[i] <- 1 -  kmeans(mydata,centers=i)$betweenss*(nrow(mydata)-i) /
                              (kmeans(mydata,centers=i)$tot.withinss-kmeans(mydata,centers=i)$betweenss)*(i-1)
#伪t方没有找到计算公式,略
plot(fakeF,type="b",xlab="Number of Clusters", ylab="fake F")

* 3.2.4 结合以上统计量一起画图比较*

#四行一列的图,调整图形左右上下的间距大小
opar<-par(mfrow=c(4,1), mar=c(4,5,3,2))
plot(1:13, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")
plot(R2,type="b",xlab="Number of Clusters", ylab="R Square")
plot(fakeR2,type="b",xlab="Number of Clusters", ylab="Fake R Square")
plot(fakeF,type="b",xlab="Number of Clusters", ylab="fake F")
par(opar)

图一wws看出,分类个数为2-4个。从3类变成2类时,伪R方急剧减小,聚类结果的解释度迅速下降,说明2类不合适,3类较合适。此时wss趋于平稳,R方也在3类变成2类时明显波动,伪F值也是波动明显,最后选定分类个数k=3

* 3.2.5 NbClust包来确定类别数量*

library(NbClust)
#默认index参数为all,会报错,出现极小值需要设置tolerance。但是这个函数过于复杂,不作手动修改
#nc <- NbClust(scale(y),min.nc=2,max.nc=15,method = "kmeans")
#index=“hubert”,通过图形判断类别个数
nc <- NbClust(scale(y),min.nc=2,max.nc=15,method = "kmeans",index = "hubert")

没有查找到Hubert统计值的具体含义,猜想左边图k=3时极小值,右图k=3到k=3+1有明显上升趋势,所以设置k=3.


3.3 选择3.2里的类别个数,类的个数k=3作聚类分析**

y <- read.csv("cluster2.csv")
#算法选择Hartigan-Wong方法,即缺省状态
km <- kmeans(scale(y),centers = 3, nstart = 20);km
#便于看出聚类后的分类情况,sort()函数对分类情况排序
sort(km$cluster)
library(fpc)
plotcluster(scale(y), km$cluster)   #生成聚类图

24个菌株分为三类,第一类:1、2、3、4、5、7、9、10、11、16、18、22
第二类:6、8、12、13、14、15、17、19、21、23、24
第三类:20
和张佳玉sas的结果完全一致。


4. 结果导出(对应sas三个过程)

* 4.1 例1.(2.1) cluster1数据的对变量分类* (k=3)

x <- read.csv("cluster1.csv")
r <- cor(x)
d<-as.dist(1-r) 
hc1<-hclust(d, "single"); hc2<-hclust(d, "complete")
hc3<-hclust(d, "median"); hc4<-hclust(d, "mcquitty")
n1 <- cutree(hc1,k=3);n2 <- cutree(hc2,k=3);n3 <- cutree(hc3,k=3);n4 <- cutree(hc4,k=3)
n1=as.data.frame(n1);n2=as.data.frame(n2);n3=as.data.frame(n3);n4=as.data.frame(n4)
n <- rbind(t(n1),t(n2),t(n3),t(n4))
x[31,] <- n[1,];x[32,] <- n[2,];x[33,] <- n[3,];x[34,] <- n[4,]
write.csv(x,"cluster1_result.csv")

* 4.2 例4.(2.4) cluster2数据对样本分类* (k=3)

y <- read.csv("cluster2.csv")
bacterium <- dist(scale(y))
hc1<-hclust(bacterium, "complete");hc2<-hclust(bacterium, "average")
hc3<-hclust(bacterium, "centroid");hc4<-hclust(bacterium, "ward.D2")
n1 <- cutree(hc1,k=3);n2 <- cutree(hc2,k=3);n3 <- cutree(hc3,k=3);n4 <- cutree(hc4,k=3)
n1=as.data.frame(n1);n2=as.data.frame(n2);n3=as.data.frame(n3);n4=as.data.frame(n4)
n <- cbind(n1,n2,n3,n4)
y$n1 <- n$n1;y$n2 <- n$n2;y$n3 <- n$n3;y$n4 <- n$n4
write.csv(y,"cluster2_result1.csv")

* 4.3 例5.(3.3) K-均值法作cluster2的聚类分析* (k=3)

y <- read.csv("cluster2.csv")
km <- kmeans(scale(y),centers = 3, nstart = 20)
y$cluster <- km$cluster
write.csv(y,"cluster2_result2.csv")

n1,n2,n3,n4(31,32,33,34)分别为最短距离法,最长距离法,中间距离法,Mcquitty相似法的分类结果

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值