SigClust计算聚类的显著性

评估聚类效果的指标非常多，甚至有工具包clusterCrit提供了多种内部和外部评价指标。

外部指标：

Czekanowski_Dice, Folkes_Mallows, Hubert, Jaccard, Kulczynski, McNemar, Phi, Precision, Rand, Recall, Rogers_Tanimoto, Russel_Rao, Sokal_Sneath1, Sokal_Sneath2，

内部指标：

Ball_Hall, Banfeld_Raftery, C_index, Calinski_Harabasz, Davies_Bouldin, Det_Ratio, Dunn, Gamma, G_plus, GDI11, GDI12, GDI13, GDI21, GDI22, GDI23, GDI31, GDI32, GDI33, GDI41, GDI42, GDI43, GDI51, GDI52, GDI53, Ksq_DetW, Log_Det_Ratio, Log_SS_Ratio, McClain_Rao, PBM, Point_Biserial, Ray_Turi, Ratkowsky_Lance, Scott_Symons, SD_Scat, SD_Dis, S_Dbw, Silhouette, Tau, Trace_W, Trace_WiB, Wemmert_Gancarski, Xie_Beni

可以基于这些指标确定最优的聚类数量，大部分的聚类分析都假设至少有两个聚类，从k=2开始计算各种指标，那么如何说明k=2比k=1，或者说数据在有聚类的情况下能被更好的解释？SigClust提供了一种思路。

原理和实现

假设从N(0,1)分布中有n次观测，并且假设我们将最小的n/2次观测作为一类，其余观测作为另一类(如下图所示)。使用t-检验测试这两个极端分类时，发现它们之间存在显著差异(p = 3.78e-13)，但是类别之间显著的差异并不能说明聚类是真实存在的。这两个聚类来自于单个高斯分布，因此不认为这样的聚类是显著的。

事实上，对于任何数据集(即使像图1中一样来自单个高斯分布)，都可以应用一种聚类算法得到这样极端的分组。这种类别之间差异并不能衡量聚类算法的显著性。因此作者提出了SigClust方法，其零假设是数据来自于单一的高斯分布，而备择假设是数据来自于非高斯分布(例如高斯混合):

• H0：数据来自于单个d维高斯分布

• Ha：数据来自于d维非高斯分布

若拒绝零假设，说明聚类显著(有效)。

简单而言，该方法先对原始数据做一个k-means的二分类，计算聚类内平方和与总平方和之比(CI/聚类指数)作为统计量；随后使用正态分布拟合原始数据，这里的关键是估计协方差矩阵(详见Liu 2008原文2.2)，还有不同的方法可以估计协方差矩阵(可以参考sigclust文档)；随后从该正态分布中重复采样，并对每个样本计算CI值获取Null分布。

(Liu2008)

用该方法计算上面极端聚类的例子，表明这样的聚类不显著。

图中

• 竖线表示真实的CI值

• 点表示从拟合正态分布中抽样计算的CI值

• 实线和虚线分别对应于对模拟CI进行的非参数密度估计和高斯密度拟合

  非参数密度估计是指在不对数据分布做任何假设的情况下，通过计算样本直方图、核密度估计等方法来估计数据的概率密度函数。高斯密度拟合是指假设数据分布服从高斯分布，然后通过最大似然估计等方法来拟合数据的概率密度函数。

• 在软件报告的图中还有一个p-vNorm值，是通过高斯密度拟合计算出的p值。
  

具体操作

sigclust(x, nsim, nrep=1, labflag=0, label=0, icovest=1)

• x: 一个矩阵或数据框，其中每行对应一个样本，每列对应一个变量。可做标准化，但不能包含缺失值。

• nsim: 模拟高斯样本数量。

• nrep: kmeans的重复数，这里默认的是1，可以提高速度。若labflag=0，该设置有效。

• labflag: 是(1)否(0)使用用户定义的标签。如果labflag = 0，则SigClust使用由2-kmeans聚类生成的标签。

• label: 如果labflag = 1，由用户指示给定聚类标签(1或2)。

• icovest: 协方差估计类型：1.使用软阈值方法(默认；Huang 2015); 2.使用样本协方差估计(建议在诊断失败时使用); 3.使用原始背景噪声阈值(“硬阈值”)。

示例1，单高斯分布，数据可能被分为高低组，但是无法通过二分类被更好的解释。

# 设置均值为5
mu <- 5
# 样本数为30
n <- 30
# 变量数为500
p <- 500
# 生成2n * p的正态分布数据
dat <- matrix(rnorm(p*2*n), 2*n, p)
# 模拟次数为1000
nsim <- 1000
# Kmeans重复次数为1
nrep <- 1
# 协方差矩阵估计方法为3
icovest <- 3
# 计算p值
pvalue <- sigclust(dat, nsim = nsim, nrep = nrep, labflag = 0, icovest = icovest)
# 绘制p值图
plot(pvalue)

示例2，双高斯分布，由于设置两组数据来自不同的高斯分布，因此聚类结果显著。数据可以通过聚类被更好地解释。

# 设置均值为5
mu <- 5
# 样本数为30
n <- 30
# 变量数为500
p <- 500
# 生成2n * p的正态分布数据
dat <- matrix(rnorm(p*2*n), 2*n, p)
# 将前n行第一列的数据加上mu
dat[1:n,1] <- dat[1:n,1] + mu
# 将后n行第一列的数据减去mu
dat[(n+1):(2n),1] <- dat[(n+1):(2n),1] - mu
# 模拟次数为1000
nsim <- 1000
# Kmeans重复次数为1
nrep <- 1
# 协方差矩阵估计方法为3
icovest <- 3
# 计算p值
pvalue <- sigclust(dat, nsim = nsim, nrep = nrep, labflag = 0, icovest = icovest)


# 绘制p值图
plot(pvalue)

以上两个例子使用的均是默认kmeans(k=2)的聚类方法。实际使用中，也可以使用用户自定义的聚类方法，具体实现就是将获取的标签作为函数的输入即可。

示例3，用户自定义聚类标签

pvalue <- sigclust(dat, nsim = nsim, nrep = nrep, labflag = 1, label=group, icovest = icovest)

注意几点

• nrep在使用默认kmeans时才有用，这个选项应该是和kmeans中的nstart对应，指的是kmeans要运行几次，为了运行速度软件默认是1。但是为了获得稳定的kmeans结果，一般nstart设置30以上。如果自己提供分组标签，该函数可以拓展到任何的聚类算法，nrep此时无效。

• icovest是估计样本协方差矩阵(特征值)的方法，共有三个选项，1.软阈值处理(推荐用于高维度且达到高斯性假设的情况，2015年文章); 2.样本特征值(推荐用于低维度且不满足高斯性假设的情况); 3.硬阈值处理(对应2008文章)。作者指出第二个选项一般而言过于保守，使用下来的感受也是如此，例如上面双高斯分布的例子都无法达到显著。

• 使用plot(res,arg='qq'); plot(res, arg='background'); plot(res, arg='diag')可以查看拟合情况以及用于诊断。

• 工具主要适用于高维小样本的数据，对于低维大样本的数据，作者推荐用icovest=2，但是结果会过于保守。在假设符合的情况下，应该也可以使用icovest=1或icovest=3进行分析(存疑)。
  

参考：https://cran.r-project.org/web/packages/sigclust/sigclust.pdf

拓展

1. sigclustTest来自于CancerSubtypes包，同样是调用sigclust，但它检测两两聚类的显著性。比如group中有三个聚类，它会分别计算1-2，1-3， 2-3聚类的显著性。唯一注意的是它输入的Data，列为sample，行为variable。

源码如下

sigclustTest=function(Data,group, nsim=1000, nrep=1, icovest=1)
{
  groupN=sort(unique(group))
  len=length(groupN)
  pvalue=matrix(data = NA, nrow =len , ncol = len)
  name=paste("Subtype",groupN)
  rownames(pvalue)=name
  colnames(pvalue)=name
  
  group_temp=sort(group,index.return = TRUE)
  data=t(Data[,group_temp$ix])
  group_temp=group_temp$x
  
  for(i in 1: (len-1))
  {
    for(j in (i+1): len)
    {
      index=which(group_temp==groupN[i] | group_temp==groupN[j])
      data1=data[index,]
      label=group_temp[index]
      label[which(label == min(label))]=1
      label[which(label == max(label))]=2 
      
      res<-sigclust(x=data1, nsim=nsim, nrep=nrep, labflag=1, label=label, icovest=icovest)
      
      plot(res,arg="pvalue",sub=paste("subtype",groupN[i],"and","subtype",groupN[j]))
      pvalue[i,j]=as.numeric(res@pval) 
      pvalue[j,i]=pvalue[i,j]
    }
  }
  diag(pvalue)=1
  pvalue
  }

https://rdrr.io/bioc/CancerSubtypes/man/sigclustTest.html

2. sigclust2工具包，将sigclust整合到了分层聚类中。算的内容都大同小异，可视化换到了ggplot。

https://github.com/pkimes/sigclust2

library(sigclust2)
set.seed(1508)
n1 <- 60; n2 <- 40; n3 <- 50; n <- n1 + n2 + n3
p <- 100
data <- matrix(rnorm(n*p), nrow=n, ncol=p)
data[, 1] <- data[, 1] + c(rep(2, n1), rep(-2, n2), rep(0, n3))
data[, 2] <- data[, 2] + c(rep(0, n1+n2), rep(sqrt(3)*3, n3))
data_pc <- prcomp(data)
par(mfrow=c(1, 2))
plot(data_pc$x[, 2], data_pc$x[, 1], xlab="PC2", ylab="PC1")
plot(data_pc$x[, 3], data_pc$x[, 1], xlab="PC3", ylab="PC1")


shc_result <- shc(data, metric="euclidean", linkage="ward.D2")

3. 发现某个补充材料有对sigclust简明扼要的介绍

https://aparajita-k.github.io/pdf/Supplementary/TCYB2990112%20Appendix.pdf

参考文献

• Kimes PK, Liu Y, Hayes DN, and Marron JS. (2017). "Statistical significance for hierarchical clustering." Biometrics.

• Huang H, Liu Y, Yuan M, and Marron JS. (2015). "Statistical significance of clustering using soft thresholding." Journal of Computational and Graphical Statistics.

• Liu Y, Hayes DN, Nobel A, and Marron JS. (2008). "Statistical significance of clustering for high-dimension, low–sample size data." Journal of the American Statistical Association.

• Xu, T., Le, T. D., Liu, L., Su, N., Wang, R., Sun, B., ... & Li, J. (2017). CancerSubtypes: an R/Bioconductor package for molecular cancer subtype identification, validation and visualization. Bioinformatics

2023/03/27