聚类分析方法学习
本文章仅供参考，如果有什么不对，请多多指教！
分类
1.层次聚类分析（150个观测值或者更小）2.划分聚类分析（大样本分析）：K-Mean聚类分析和基于中心点的划分
注意
1. 层次聚类分析
1. 选择合适的变量；
2.选择变量前先去掉部分不可取变量，缩放数据，可用标准化，最大最小标准化等方法；
3.寻找异常点，通过outliers包中的函数进行筛选和删除异常单变量离群，或者用mvoutliers包围绕着中心点进行划分；
4.计算距离，一般是欧几里得距离，其中包括曼哈顿距离、最大距离等，可以调用dist（）函数进行计算；
5.选择聚类算法，样本少的采用层次聚类分析，样本多的采用划分聚类分析；
6.获得一种或者多种的聚类分析方法；
7.确定类的数目，一般调用NbClust（）函数去确定；
8.获得最终的聚类解决方案；
9.结果可视化，一般观察树状图；
10.解读类；
11.验证结果，利用fpc、clv和clValid包包含评估聚类解的稳定性的函数。
2. 划分聚类分析：K-均值聚类和基于中心点的划分（PAM）
2.1 K-均值聚类分析
1.选择k个中心点；
2.把每个数据点分配到离它最近的中心点；
3.重新计算每类中的点到该类中心点距离的平均值；
4.分配每个数据到它最近的中心点；
5.重复3和4的步骤直到所有观测值不再被分配或是达到最大的迭代次数（默认：10次）；
其中停止的时候使得观测值到其指定的聚类中心的平方总和最小，即
最小。
6.注意：平方总和最小的条件限制了所有变量必须是连续的，而且容易被异常值所影响；同时对初始中心值的选择十分敏感，可用nstart()函数尝试多种初始配置并输出最好的那个。其中k的取值跟层次分析一样重要，可以用NbClust()函数进行参考，也可以借助类中总的平方值对聚类数量的曲线，也就是代码中的wssplot()函数就可以绘画出SSE函数。
2.2 基于中心点的划分
K-均值聚类分析是利用欧几里得距离进行聚类分析，而Pam聚类分析是采用任意距离，可以容纳混合数据类型，不仅限于连续变量。PAM算法如下：
1.随机选择K个观测值（每个都称为中心点）；
2.计算观测值到各个中心的距离/相异性；
3.把每个观测值分配到最近的中心点；
4.计算每个中心点到每个观测值的距离的总和（总成本）；
5.选择一个该类中不是中心的点，并和中心点互换；
6.重新把每个点分配到距离它最近的中心点；
7.再次计算总成本；
8.如果总成本比步骤（4）计算的总成本少，把新的点作为中心点；
9.重复步骤（5）-（8）直到中心点不再改变.
R语言代码
1. 计算距离
data(nutrient,package="flexclust")
head(nutrient,4)
d<-dist(nutrient)
as.matrix(d)[1:4,1:4]#强制转化为矩阵
注：可用cluster包中的daisy()函数获得包含任意二元、名义、有序属性组合的相异矩阵。
2. 层次聚类分析
data(nutrient,package="flexclust")
row.names(nutrient)<-tolower(row.names(nutrient))#tolower()将字体全部改成小写，toupper()将字体全部改成大写
nutrient.scaled<-scale(nutrient)#对数据进行标准化
d<-dist(nutrient.scaled)
fit.average<-hclust(d,method="average")#对距离进行层次聚类分析
dev.new()
plot(fit.average,hang=-1,cex=.8,main="Average Linkage Clustering")
注：层次聚类分析就是将距离最短的两类合并成一类，从而把所有类合并成单个类为止，其中距离定义有：1.单联动（细长、雪茄型数据）；2.全联动（大致相等的直径紧凑型）；3.平均联动（偏向于方差较小的聚类）；4.质心（变量均值向量之间的距离）；5.Ward法（所有变量的方差分析的平方和）。层次聚类方法用hclust()函数，其中method=single,complete,average,controid,ward。
-
选择聚类的个数
library(NbClust)
dev.new()
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="Name of Clusters",ylab="Name of Criteria",main="Number of Cluster Chosen by 26 Criteria")
注：根据最后一个直方图，值最大的就用K的数目，观察上图，可以发现K=2、3、5、15是最好的选择，最后需要测试一下哪个K值最好。
-
最后的方案
clusters<-cutree(fit.average,k=5)#对进行聚类分析的变量进行绘画树状图，并分为5类
table(clusters)
aggregate(nutrient,by=list(cluster=clusters),median)#对组合后的数据进行求中位数
aggregate(as.data.frame(nutrient.scaled),by=list(cluster=clusters),median)
dev.new()
plot(fit.average,hang=-1,cex=.8,main="Average Linkage Clustering\n5 Cluster Solution")#对进行聚类分析的变量进行重新绘制树状图
rect.hclust(fit.average,k=5)#对进行聚类分析的变量叠加 5类的解决方案
注：生物科学中经常使用层次聚类分析方法，不过这种算法是贪婪的，一旦一个观测值被分配给一个类，它就不能在后面进行重新分配，如果在大样本中，用划分聚类分析比较合适。
3. 划分聚类分析
3.1 K-均值聚类分析
-
生成计算SSE函数
wssplot <- function(data, nc=15, seed=1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))#此处使用apply函数
for (i in 2:nc){
set.seed(seed)#保证结果是可复制的
wss[i] <-sum(kmeans(data,centers=i)$withiness)#withiness是组内平方和
}
plot(1:nc,wss,type="b",xlab="Number of Clusters", ylab="Within groups sum of squares")}
-
读取葡萄酒数据
library(rattle)
data(wine,package="rattle")
head(wine)
-
确定K的个数
library(NbClust)
df <- scale(wine[-1])
wssplot(df)
set.seed(1234)
dev.new()
devAskNewPage(ask=TRUE)
nc <- NbClust(df,min.nc=2,max.nc=15,method="kmeans")
table(nc$Best.n[1,])
barplot(nc$Best.n[1,],xlab="Number of Clusters",ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria")
-
从上述的图表，可以看到聚类数目K=3，对df进行K均值聚类
set.seed(1234)
fit.km <- kmeans(df,3,nstart=25)
fit.km$size
fit.km$centers
aggregate(wine[-1],by=list(cluster=fit.km$cluster),mean)#总结该变量的平均值
ct.km <- table(wine$Type,fit.km$cluster)
ct.km
library(flexclust)
randIndex(ct.km)#兰德指数：为两种划分提供了一种衡量两种分区之间的协定：K-均值聚类分析和PAM聚类分析
3.2 PAM聚类分析
library(cluster)
set.seed(1234)
fit.pam<-pam(wine[-1],k=3,stand=TRUE)#stand表示该变量是否标准化
fit.pam$medoids#输出中心点
clusplot(fit.pam,main="Bivariate Cluster Plot")
ct.pam <- table(wine$Type,fit.pam$clustering)#此处的cluster和clustering一样
ct.pam
randIndex(ct.pam)#由于pam中的兰德指数比K-均值的兰德指数低，所以建议选择K-均值聚类分析
-
避免不存在的类
library(fMultivar)
library(NbClust)
library(cluster)
library(ggplot2)
set.seed(1234)
df <- rnorm2d(1000,rho=.5)#rnorm()生成标准正态分布的随机数，而rnorm2d()生成相关系数为0.5的正态分布
df <- as.data.frame(df)
plot(df,main="Bivariate Normal Distribution with rho=0.5")
wssplot(df)
nc <- NbClust(df,min.nc=2,max.nc=15,method="kmeans")
dev.new()
barplot(table(nc$Best.n[1,]),xlab="Number of Clusters",ylab="Number of Criteria",main="Number of Chosen by 26 Criterias")#建议k=3
fit <- pam(df,k=2)
df$clustering <- factor(fit$clustering)
ggplot(data=df,aes(x=V1,y=V2,color=clustering,shape=clustering))+geom_point()+ggtitle("Clustering of Bivariate Normal Data")
plot(nc$All.index[,4],type="o",ylab="CCC",xlab="Number of clusters",col="blue")
注：最后一个plot中的CCC是调用了NbCluster包中的立方聚类规则，揭示不存在的类，通过分析，该图中的CCC值都是负值并且对于两类或是更多的类递减时，就是典型的单峰分布，就是偏左正态分布或者偏右正态分布。
