聚类分析的簇权重

#聚类分析

#聚类分析是一种数据归约技术，旨在揭露一个数据集中观测值的子集
#它可以把大量的观测值规约为若干个类
#这里的类被定义为若干个观测值组成的群组，群组内观测值的相似度比群间相似度高

#最常用的聚类方法：
#方法一，层次聚类
#在层次聚类中，每一个观测值自成一类，这些类每次两两合并，直到所有类被聚成一类为止
#常用算法包括：单联动，全联动，平均联动，质心，和Ward方法

#方法二，划分聚类
#在划分聚类中，首先指定类的个数k,然后观测值被随机分成k类，再重新形成聚合的类
#常用算法包括：k-means(k均值)，围绕中心点的划分（PAM）

#聚类分析的一般步骤
#第一步：选择合适的变量
#第二步：缩放数据
#df1 <- apply(mydata,2,function(x){(x-mean(x)/sd(x))})  #标准化,u=0,sd=1
#df2 <- apply(mydata,2,fuction(x){x/max(x)})#除以最大值
#df3 <- apply(mydata,2,fuction(x){(x-mean(x)/mad(x))}) #减去均值并除以变量的平均绝对偏差

#第三步：寻找异常点
#第四步：计算距离
#第五步：选择聚类算法
#第六步：获得一种或多种聚类方法
#第七步：确定类的数目
#第八步：获得最终的聚类解决方案
#第九步：解读类
#第十步：验证结果

#批量安装包,感觉像是并行操作
install.packages(c("cluster","Nbclust","flexclust","fMultivar","ggplot2","rattle"))

#计算距离

library(lattice)
library(modeltools)
library(stats4)
library(flexclust)  #library支持加载函数，加载数据无效object 'nutrient' not found

data(nutrient,package="flexclust") #加载进入全局环境变量
head(nutrient,4) #查看前四行数据
dim(nutrient) #返回 27,5，即27行，5列

#dist(x,method=)用来计算矩阵或数据框中所有行之间的距离
#x         输入数据
#method    默认为欧几里得距离

#补充说明:
#欧几里得距离通常做为连续型数据的距离度量
#但是如果存在其他类型的数据，则需要相异的替代措施，你可以使用cluster包中
#daisy()函数来获得包含任意二元(binary),名义(nominal)，有序(ordinal),
#连续(continuous)属性组合的相异矩阵

d <- dist(nutrient)
class(d) #返回dist对象
m <- as.matrix(d)#强制转换成matirx对象
dim(m) #返回27行，27列
m[1:4,1:4] #查看距离计算后结果矩阵的前四行
#结果解读：
#以BEEF BRAISED，HAMBURGER数据为例，计算欧几里得距离
#原始数据：
#energy protein fat calcium iron
#BEEF BRAISED    340      20  28       9  2.6
#HAMBURGER       245      21  17       9  2.7
#d方=(340-245)^2+(20-21)^2+(28-17)^2+(9-9)^2+(2.6-2.7)^2=95.64^2
#d=95.64

#观测值之间的距离越大，异质性越大
#观测值和它自己之间的距离为0

#层次聚类
#优点：
#当需要嵌套聚类和有意义的层次结构时，层次聚类或许特别有用
#在生物科学中这种情况很常见
#在某种意义上分层算法是贪婪的，一旦一个观测值被分配给一个类，它就不能在后面的
#过程中被重新分配
#层次聚类难以应用到数百甚至数千观测值的大样本中

#算法
#1定义每个观测值（行或单元）为一类
#2计算每类和其他各类的距离
#3把距离最短的两类合并成一类，这样类的个数就减少一个
#4重复步骤2，步骤3，直到包含所有观测值的类合并成单个的类为止

#方法
#single单联动：一个类中的点和另一个类中的点的最小距离
#倾向于发现细长，雪茄型的类。它通常展示一种链式的现象，即不相似的观测值分到一类中
#因为它们和它们的中间值很想像

#complete全联动：一个类中的点和另一个类中的点的最大距离
#倾向于发现大致相等的直径紧凑类
#它对异常值很敏感

#average平均联动：一个类中的点和另一个类中的点的平均距离
#平均联动提供了以上两种方法的折中
#它不像链式，而且对异常值没有那么敏感
#它倾向于把方差小的类聚合

#centroid质心：两类中质心(变量均值向量)之间的距离。对单个的观测值来说，质心就是变量的值
#很受欢迎
#但是它不能如平均联动法，或ward法表现的好

#Ward法：两个类之间所有变量的方差分析的平方和
#倾向于把少量观测值的类聚合到一起，并且倾向于产生与观测值个数大致相等的类

#范例一  营养数据的平均联动聚类
data(nutrient,package="flexclust")
row.names(nutrient) <- tolower(row.names(nutrient)) #将大写行名替换成小写行名
nutrient.scaled <- scale(nutrient) #将每一列值，标准化
class(nutrient.scaled) #返回矩阵对象，作为新的数据源，而不在使用原始数据源
dim(nutrient.scaled ) #返回 27行，5列
d <- dist(nutrient.scaled) #基于标准化后的数据源，计算欧几里得距离
fit.average <- hclust(d,method="average") #基于平均联动average,生成层次聚类模型
fit.average

#树状图应该从下往上读，它展示了这些条目如何被结合成类
#每一个观测值起初自成1类，然后基于距离，将相距最近的两类（beef braised和smoked ham）合并
#然后合并继续，直到所有的观测值都合并成一类
#高度刻度代表了该高度类之间合并的判断值

plot(fit.average,hang=-1,cex=0.8,
     main="Average Linkage Clustering")#基于模型画图，hang=-1,表示将标签悬挂在y=-1处理
#结果解读：
#tuna canned 和 chiicken canned是相似的
#但它们与clams canned有很大的不同
#辅助于理解基于食物营养成分的相似性和相异性

#范例二 选择聚类的个数
library(NbClust)
#Prompt before New Page
#options("device.ask.default")
devAskNewPage(ask=TRUE) #弹出窗口;Hit <Return> to see next plot
devAskNewPage(ask=FALSE)

#NbClust(),确定在一个聚类分析里面类的最佳数目
#输入包括需要做聚类的矩阵或是数据框
#distance,使用的距离测度
#method,使用的聚类方法
#min.nc,最小的聚类个数
#max.nc,最大的聚类个数
nc <- NbClust(nutrient.scaled,distance="euclidean",
              min.nc=2,max.nc=15,method="average")
class(nc) #返回list

#在min.nc和max.nc之间使用多个集群获得的数据集的每个分区的索引值。
nc$All.index
dim(nc$All.index)  #返回14*26

#每个指标提出的最优聚类数及其对应的指标值。
nc$Best.nc 
dim(nc$Best.nc) #返回 2*26
class(nc$Best.nc) #返回matirx
nc$Best.nc[1,] #第一行，对应Number_clusters信息，观测值对应组别信息
nc$Best.nc[2,] #第二行，对应Value_Index信息

#基于组别信息，统计频数值
table(nc$Best.nc[1,]) 
#0  1  2  3  4  5  9 10 13 14 15 --第一行，聚类个数
#2  1  4  4  2  4  1  1  2  1  4 --第二行，判断原则赞同聚类个数
#第三列：四个评判准则赞同聚类个数为2
#第四列：四个评判准则赞同聚类个数为3
#第六列：四个评判准则赞同聚类个数为5
#第十一列：四个评判准则赞同聚类个数为15

#一个观测值，一行观测数据，即对应一条评判规则
#第二行求和：2+1+4+4+2+4+1+1+2+1+4=26，共计26行观测值，对应26条评判规则
barplot(table(nc$Best.nc[1,]),
        xlab="Number of Clusters",
        ylab="Number of Criteria",
        main="Number of Clusters Chosen by 26 Criteria")

#范例三  获取最终的聚类方案，在2,3,5,15组中进行取舍判断,假设取5

#cutree(),将树状分成五类
clusters <- cutree(fit.average,k=5)
class(clusters) #返回向量
length(clusters)#返回27
rownames(clusters) #返回NULL
colnames(clusters) #返回NULL
is.vector(clusters) #返回TRUE
names(clusters)  #返回向量各对象的名称，对应27个观测值的名字

clusters[1] 
#beef braised 
#1 

clusters["beef braised"]
#beef braised 
#1

table(clusters)
#结果解读：
#1  2  3  4  5 
#7 16  1  2  1

#第一类，有7个观测值
#第二类，有16个观测值
#第三类，有1个观测值
#第四类，有2个观测值
#第五类，有1个观测值

#基于cluster对象中包含的27个观测值，分别对应的组别信息，进行分组
#利用aggregate,基于原始数据集群的分组情况，计算每组，每一列的中位数
aggregate(nutrient,by=list(cluster=clusters),median)

#利用aggregate,基于标准化数据集群的分组情况，计算每组，每一列的中位数
aggregate(as.data.frame(nutrient.scaled),by=list(Cluster=clusters),median)

#模型，原树状图
plot(fit.average,hang=-1,cex=0.8,
     main="Average Linkage Clustering\n5 Cluster Solution")
#在原树状图上，叠加五类的解决方案，即画框图
#各个类下观察值的个数，与table(clusters)的频数统计结果一致
rect.hclust(fit.average,k=5)

#划分聚类分析
#

#k均值聚类
#算法：
#1 选择k个中心点（可随机选择k行）
#2 把每个数据点分配到离它最近的中心点
#3 重新计算每类中的点到该类中心点距离的平均值
#4 分配每个数据到它最近的中心点
#5 重复步骤3，步骤4直到所有的观测值不再被分配或是达到最大的迭代次数，R把10次作为默认迭代次数

#优点：
#k均值聚类能处理比层次聚类更大的数据集
#另外，观测值不会永远被分到一类中
#当我们提高整体解决方案时，聚类方案也会改动
#但是均值的使用意味着所有的变量必须是连续的，并且这个方法很有可能被异常值影响

#范例一 k均值聚类
data(wine,package="rattle")#从数据包加载进全局环境变量
dim(wine)  #返回178行,14列(1列对应观测值名，13列对应化学成分)

head(wine) #默认加载前六行数据

dim(wine[,-1]) #删除酒名列，返回178,13
dim(wine[-1])  #删除酒名列，返回178,13
df <- scale(wine[,-1])#将原数据集，进行标准化操作
dim(df) #返回178,13

#apply(df,2,var) 针对每一列求方差，返回向量，每一个值代表当前列的方差
#sum(apply(df,2,var)) 所有列方差求和
wss <- (nrow(df)-1)*sum(apply(df,2,var)) #返回总平方和，等价于k$totss

#由于k均值聚类每次都是随机选择k个点，所以每次调用函数可能获得不同的方案
#使用set.seed函数可以保证结果是可复制的
set.seed(1234)
#kmeans(x,centers,nstart)
#x         表示数值型数据集
#centers   表示要提取的聚类数目
#nstart    表示尝试多种初始配置并输出最好的一个

k <- kmeans(df,centers=2)
#返回向量列表，利用类别编号对观测值进行标识
#A vector of integers (from 1:k) indicating the cluster to which each point is allocated.
k$cluster

#返回类中心
#A matrix of cluster centres.
k$centers

#返回总平方和，总平方和=所有组内平方和+所有组间平方和
#The total sum of squares.
k$totss        #返回 2301 （884.3435 + 765.0965 + 651.56=2301）

#返回各组内平方和
#Vector of within-cluster sum of squares, one component per cluster.
k$withinss     #返回  884.3435 765.0965

#返回各组内平方和的总和
#Total within-cluster sum of squares, i.e. sum(withinss).
k$tot.withinss #返回  1649.44(884.3435+765.0965=1649.44)

#返回组间平方和
#The between-cluster sum of squares, i.e. totss-tot.withinss.
k$betweenss    #返回  651.56

sum(kmeans(df,centers=2)$withinss) #返回1649.44，等价于k$tot.withinss 

#第一种方法，盲选，合适的聚类个数

#肘部法则
#如果问题中没有指定的值，可以通过肘部法则这一技术来估计聚类数量。
#肘部法则会把不同值的成本函数值画出来。
#随着值的增大，平均畸变程度会减小；每个类包含的样本数会减少，于是样本离其重心会更近。
#但是，随着值继续增大，平均畸变程度的改善效果会不断减低。值增大过程中，畸变程度的改善效果下降幅度最大的位置对应的值就是肘部。

wssplot <- function(data,nc=15,seed=1234){
    wss <- (nrow(data)-1)*sum(apply(data,2,var)) #得到一个数值
    for(i in 2:nc){
        #为了保证结果的可复制性，必须设定随机种子
        set.seed(seed)
        #得到以i为聚类个数(等价于组别个数)，求得所有组内平方和，存储到矩阵中
        wss[i] <- sum(kmeans(data,centers=i)$withinss) 
    }
    #基于聚类个数，以及当前聚类下素有组内平方和画图
    plot(1:nc,wss,type="b",xlab="Number of Clusters",
         ylab="Within groups sum of squares")
}

wssplot(df)
#结果解读：
#从一类到三类变化时，组内的平方总和有一个明显的下降趋势
#三类以后，下降的速度减弱
#暗示聚成三类可能对数据来说是一个很好的拟合

#第二种方法，利用NbClust来选择合适的聚类个数
library(NbClust)
set.seed(1234)
devAskNewPage(ask=TRUE)
nc <- NbClust(df,min.nc = 2,max.nc=15,method="kmeans")
#                   ***** Conclusion *****                            
#* According to the majority rule, the best number of clusters is  3
#NbClust，同样建议选3个聚类

table(nc$Best.nc[1,])
# 0  1  2  3 10 12 14 15 
# 2  1  4 15  1  1  1  1 
# 分3组时，有15个判断原则投票

barplot(table(nc$Best.n[1,]),
        xlab="Number of Clusters",
        ylab="Number of Clusters Chosen by 26 Criteria")

#将确定后的centers=3带入kmeans函数，生成模型
#为了保证结果的可复制性，设定随机种子1234
set.seed(1234)
fit.km <- kmeans(df,3,nstart=25)
fit.km$size #返回每一个组的观测值数量，62 65 51
fit.km$centers #返回每一个组的centers,基于标准化的数据运算
#基于原始数据，计算每一组下每一列的mean值
aggregate(wine[-1],by=list(clusters=fit.km$cluster),mean)

#范例  
class(wine) #返回 data.frame
wine[1]     #返回wine数据源的，第一列，即type列,但包含了标题行
class(wine[1]) #返回 "data.frame"
dim(wine[1])  #返回178行*1列

wine$Type  #返回wine数据源的，第一列，即type列
class(wine$Type) #返回"factor"
dim(wine$Type)   #Null
is.vector(wine$Type) #返回FALSE
length(wine$Type) #返回178

fit.km$cluster            #返回每个观测值对应组号的向量列表
class(fit.km$cluster)     #返回 "integer"
is.vector(fit.km$cluster) #返回TRUE
length(fit.km$cluster)    #返回178

#table(var1,var2,...),要求var1...n必须是因子变量
#基于Type列，fit.km$cluster值生成二维表格
#wine$Type，对应行维度
#fit.km$cluster，对应列维度
ct.km <- table(wine$Type,fit.km$cluster) 
ct.km

#范例  兰德指数(Rand index)来量化类型变量和类之间的协议
#调整的兰德指数为两种划分提供了一种衡量两个分区之间的协定，即调整后机会的量度
#它的变化范围从-1(不同意)到1(完全同意)
library(flexclust)
randIndex(ct.km)
#     ARI 
#0.897495 
#结果解读;
#葡萄酒品种类型（wine$Type）和聚类的解决方案之间的协定是0.9
#结果不坏

#围绕中心点的划分

#因为k均值聚类方法是基于均值的，所以它对异常值是敏感的
#一个更稳健的方法是围绕中心点的划分。
#k均值聚类一般使用欧几里得距离
#而PAM可以使用任意的距离来计算，它可以容纳混合数据类型，并且不仅限于连续变量

#PAM算法如下：
#1 随机选择k个观测值（每个都称为中心点）
#2 计算观测值到各个中心的距离/相异性
#3 把每个观测值分配到最近的中心点
#4 计算每个中心到每个观测值的距离的总和（总成本）
#5 选择一个该类中不是中心的点，并和中心点互换
#6 重新把每一个点分配到距它最近的中心点
#7 再次计算总成本
#8 如果总成本比步骤（4）计算的总成本少，把新的点作为中心点
#9 重复步骤5~8直到中心点不再改变

#范例一  对葡萄酒数据使用基于质心的划分方法

library(cluster)
set.seed(1234) #设定随机种子，保证结果的复制
#pam(x,k,metirc="euclidean",stand=FALSE)
#x 代表数据矩阵或数据框
#k 表示聚类的个数
#metric 表示相异性/相似性的度量
#stand  表示一个逻辑值，表示是否有变量应该在计算该指标之前被标准化

#wine[-1] , 表示wine数据去除第一列type值后，作为输入数据传入
#stand=TRUE 表示基于原始值进行的中心点方法
fit.pam <- pam(wine[-1],k=3,stand=TRUE)
fit.pam$medoids #输出中心点
fit.pam$clustering #输出各个观测值，被组标签标记后的向量
fit.pam$clustering[36] #第36个观测值对应的组号，1
fit.pam$clustering[107] #第36个观测值对应的组号，2
fit.pam$clustering[175] #第36个观测值对应的组号，3

table(fit.pam$clustering)
# 1  2  3 
#75 54 49 
#分成3类，分别包含75个观测值，54个观测值，49个观测值
#75+54+49=178，总观测值

clusplot(fit.pam,main="Bivariate Cluster Plot")
#结果解读

#范例  兰德指数

#table(var1,var2,...),要求var1...n必须是因子变量
#基于Type列，fit.pam$cluster值生成二维表格
#wine$Type，对应行维度
#fit.pam$cluster，对应列维度
ct.pam <- table(wine$Type,fit.pam$clustering)
ct.pam

#范例  兰德指数(Rand index)来量化类型变量和类之间的协议
#调整的兰德指数为两种划分提供了一种衡量两个分区之间的协定，即调整后机会的量度
#它的变化范围从-1(不同意)到1(完全同意)
randIndex(ct.pam)
#     ARI 
#0.6994957  
#结果解读;
#葡萄酒品种类型（wine$Type）和聚类的解决方案之间的协定是0.6994957 
#结果不好，从0.9减低到了0.6994957 

par(mfrow=c(1,1))

#避免不存在的类

#范例

#定义随机的数据源
library(timeDate)
library(timeSeries)
library(fBasics)
library(fMultivar)

set.seed(1234)
#双变量正态随机偏差
#从相关系数为0.5的二元正态分布中抽取1000个观测值
df <- rnorm2d(1000,rho=0.5) 
df <- as.data.frame(df)
dim(df) #返回 1000行 2列
#完全随机，没有规律，所以不存在任何类
plot(df,main="Bivariate Normal Distribution with rho=0.5")

#方法一  肘部法则，盲选聚类个数
wssplot(df) 
#建议选三个，后续变化比较平缓

#方法二，NbClust选取聚类个数
library(NbClust)
nc <- NbClust(df,min.nc=2,max.nc=15,method="kmeans")
#***** Conclusion *****                            
    
#    * According to the majority rule, the best number of clusters is  2 
#NbClust推荐选2个聚类变量

dev.new()  #新打开画布窗口，画条形图
barplot(table(nc$Best.nc[1,]),
        xlab="Number of Clusters",
        ylab="Number of Criteria",
        main="Number of Clusters Chosen by 26 Criteria")

#范例 使用PAM法进行双聚类分析

library(ggplot2)
library(cluster)
#基于确认的2聚类，进行中心点的划分，生成中心点划分模型
fit <- pam(df,k=2)
#将模型中各观测值对应的组号，向量因子化，追加到数据源的新列
df$clustering <- factor(fit$clustering)
dim(df) #返回1000行，3列

#ggplot()+geom_point()+ggtile(),
#注意+必须位于第一个函数右括号旁，当时另起一行，将失败
#将没有类的随机数据，任意划分成两类
ggplot(data=df,aes(x=V1,y=V2,color=clustering,
                   shape=clustering))+
geom_point()+
    ggtitle("Clustering of Bivariate Normal Data")

#返回14*26
dim(nc$All.index) #返回14*26
plot(nc$All.index[,4],type="o",
     ylab="CCC",xlab="Number of clusters",
     col="blue")    
#nc$All.index[,4]第四列值对应CCC列
#立方聚类规则(Cubic Cluster Criteria CCC)往往可以帮助我们揭示不存在的结构
#当CCC值为负，并且对于两类或是更多的类递减时，就是典型的单峰分布
#如果你视图找出的某种意义上的“真实的”类，而不是一个方便的划分，就要确保结果
#是文件的并且是可重复的。
#如果同一类持续复原，你就可以对得出的结果更加自信