R语言如何对模型进行比较 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相似法的分类结果