使用SigClust对Kmeans二聚类结果进行显著性检验

动机

1. 假设从N(0,1)分布中有n次观测（n个数据点），将最小的n/2个数据点作为一类，其余数据点作为另一类
2. 使用双样本t-检验比较该二分类的均值时，发现它们之间存在显著差异，但是类别之间显著的差异并不能说明聚类是真实存在的
3. 实际上这两个聚类来自于单个高斯分布总体N(0,1)，因此不能认为这样的聚类是显著的

SigClust原理

1. SigClust方法，其零假设是数据来自于单一的高斯分布，而备择假设是数据来自于非高斯分布
2. 若拒绝零假设，说明聚类显著（有效）
3. 简单而言，该方法先对原始数据做一个k-means的二分类，计算聚类指数CI=聚类内平方和与总平方和之比作为统计量
4. 随后使用正态分布拟合原始数据，这里的关键是估计协方差矩阵

R代码

#加载gdata库
install.packages("gdata")
library(gdata)

#读取数据
dat <- read.xls("Cluster_data.xlsx") 

#特征归一化
for (col in 1:4) {
  dat[,col] = scale(dat[,col])
}

#向量标准化
for (row in 1:44) {
  dat[row,] = dat[row,]/norm(dat[row,],type=c("2"))
}

plot(dat)

#Kmeans聚类
km <- kmeans(dat,2,nstart=1000) #聚类两类，重复次数1000次
km_cluster <- unlist(km["cluster"]) #提取聚类标号

plot(km_cluster)

#存聚类标签
install.packages("openxlsx")
library(openxlsx)

km_cluster <- data.frame(km_cluster) #创建data frame

openxlsx::write.xlsx(km_cluster , "Cluster_data_label.xlsx" , sheetName = "Sheet1" , append = TRUE )

#SigClust检验

nsim <- 100000 #模拟次数100000

icovest <- 2 #协方差估计方法2，适用低纬度且不满足高斯性假设

pvalue <- sigclust(dat , nsim = nsim , labflag = 1, label = km_cluster , icovest = icovest) #计算p值

plot(pvalue)