聚类分析笔记

前言

以前看过StatQuest的聚类分析笔记,可以参考这篇笔记《StatQuest学习笔记17——聚类》

聚类的概率与类型

聚类分析法的概念

聚类分析法(cluster analysis)是研究“物以类聚”的一种现代统计分析方法,聚类的目的是把分类对象按照一定的规则分成若干类,这些类不是事先给定的,而是根据数据的特征确定的,在同一类中这些对象在某种意义上趋向于彼此相似,而在不同类中对象趋于不相似。

聚类分析法的类型

在实际问题中,经常要对一些东西进行分类,例如在古生物研究中,通过挖掘出来的一些骨髓的形状和大小对它们进行科学的分类,在这里,骨髓的形状和大小是我们用来分类的依据,称为指标或变量,用$X_{1},X_{2},\cdots,X_{p}$表示,其中,p是变量的个数,需要进行分类的骨髓为样品,用1,2,3,…n表示,其中n是样本的个数,聚类分析的数据结果见下表:

mark

在聚类分析中,基本的思想是认为所研究的样本或变量存在着程度不同的相似性(亲疏关系)。于是,根据一批样品中的多个观测指标,具体找出一些能够度量样本之间相似程度的统计量,以这些统计量为划分类型的依据,把一些相似程度较大的样本聚合为一类,把另外一些彼此之间相似程度较大的样品又聚合为另一类,关系密切的聚合到一个小的分类单位,关系疏远的聚合到一个大的分类单位,直到把所有的样本都又聚合完毕,把不同类型一一划分出来,形成一个由小到大的分类系统,最后把整个分类系统画成一张聚类图,用它把所有样本间的亲疏关系表示出来。

通常根据分类对象的不同,可以将聚类分为两类:一类是对样品进行分类处理,叫Q型;一类是对变量进行分类处理,叫R型。Q型聚类又叫样品分类,就是对观测对象进行聚类,是根据被观测的对象各种特征进行分类。在生物信息学中,经常使用的是R型,例如有对照组1,2,3三个样本,实验组1,2,3这三个样本,把这6个样本把基因表达谱的数据进行聚类。

聚类统计量

聚类分手析的基本原则是将有较大相似性的对象归为同一类,而将差异较大的个体归入不同的类。为了将样本聚类,就需要研究样品之间的关系,一种方法是将每一个样本看做p维空间的一个点,并在空间定义距离,距离较近的点归为一类,距离较远的点属于不同的类。对变量通常计算它们的相似系数,性质越接近的变量,它们的相似系数越接近于1或-1,彼此无关的变量的相似系数越接近于0。可进行聚类的统计量有距离和相似系数。距离有欧氏距离,马氏距离,兰氏距离,相似系数有夹角余弦,相关系数。

对样品进行聚类时,我们将把样品间的“靠近”程度用某种距离来刻画对指标的聚类,往往用某种相似系数来刻画。当选用n个样本,p个指标时,就可以得到一个$n\times p$的数据矩阵,即$X=(x_{ij})n \times p$,该距离的元素$x_{ij}$表示第i个样本的第j个变量值。

距离和相似系数

对样品或变量进行分类时,我们常用距离和相似系数对样品或变量之间的相似性进行度量。距离常用来度量样本之间的相似性,而相似系数常用来度量变量间的相似性。

距离多用于样品的分类,信$d_{ij}$表示样本$x_{i}$和$x_{j}$的距离,一般要求$d_{ij}$满足以下四个条件,如下所示:

mark

在聚类分析中并不严格要求定义的距离都满足这四条,一般来说前三条大部分是能满足的,一些不能满足(4),但是在广义的角度上也称其为距离。

设$x_{ij}(i=1,2,\cdots,n;j=1,2,\cdots,p)$为第i个样本的第j个指标的观测数据。即若每个样本有p个变量,则每个样本都可以看成p维空间中的一个点,n个样品就是p维空间中的n个点,定义$d_{ij}$为样本xj和xj的距离,于是得到一个$n\times n$的距离矩阵,如下所示:

mark

样品聚类都是基于此距离矩阵进行的,现在我们看一个简单的案例:

距离:案例A:样本聚类

以下列举5个观察值,2个变量数据的平面散点图,如下所示:

1
2
3
4
x1 <- c(5,7,3,6,6)
x2 <- c(7,1,2,5,6)
plot(x1,x2)
text(x1,x2,labels=c(1:5),adj=-.5)

由于只有2个变量,所以从散点图上就可以直观地将这5个样本分为几类,但当变量较多时,这种方法显著是不行的。为了计算平面上各点之间的距离$d_{ij}$,在聚类分析中对连续型变量常用的距离有以下3种,如下所示:

明氏距离(Minkowski):

$d_{i j}(q)=\left[\sum_{k=1}^{p}\left(x_{i k}-x_{j k}\right)^{q}\right]^{\frac{1}{q}}$

q=1时,$d_{i j}(1)=\sum_{k=1}^{p}\left|x_{i k}-x_{j k}\right|$,称为绝对值距离(manhattan);

q=2时,$d_{i j}(2)=\left[\sum_{k=1}^{p}\left(x_{i k}-x_{j k}\right)^{2}\right]^{\frac{1}{2}}$,称为欧氏距离(euclidean);

当$q=\infty$时,$d_{i j}(\infty) \max _{1 \leqslant k \leqslant D}\left|x_{i k}-x_{j k}\right|$,称为切比雪夫距离(maximun)。

马氏距离(Mahalanobi):

$d_{i j}(M)=\left(x_{i}-x_{j}\right)^{\prime} \Sigma^{-1}\left(x_{i}-x_{j}\right)$

其中:$x_{i}$为样本$i$的p什指标组成的行徽号$\sum$为协方差阵。

优点:马氏距离既排除了各指标间的相关性干扰,又消除了各指标的量纲。

缺点:样本协方差中阵在聚类过程中不变,这点不合理。

兰氏距离(Canberra)

$d_{i j}(L W)=\frac{1}{p} \sum_{k=1}^{p} \frac{\left|x_{i k}-x_{j k}\right|}{x_{i k}+x_{j k}} \quad\left(x_{i j}>0\right)$

下面我们使用dist()函数计算一下案例A的欧氏距离与马氏距离,用法如下所示:

1
2
3
4
5
6
dist(x,mdehtod="euclidean",diag=FALSE,upper=FALSE,p=2)
# x是数据矩阵
# method是计算方法,包括euclidean,maximum,manhattan,canberra,binary,minkowski
# diag为是否包含对角线元素
# upper为是否需要上三角
# p为Minkowski距离的幂次

案例A的计算结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
> X<-as.matrix(data.frame(x1,x2))
> dist(X) # 默认情况下是欧氏距离
1 2 3 4
2 6.324555
3 5.385165 4.123106
4 2.236068 4.123106 4.242641
5 1.414214 5.099020 5.000000 1.000000
> dist(X,diag=T,upper=T) # 默认情况下是欧氏距离
1 2 3 4 5
1 0.000000 6.324555 5.385165 2.236068 1.414214
2 6.324555 0.000000 4.123106 4.123106 5.099020
3 5.385165 4.123106 0.000000 4.242641 5.000000
4 2.236068 4.123106 4.242641 0.000000 1.000000
5 1.414214 5.099020 5.000000 1.000000 0.000000
> dist(X,method="manhattan",diag=T,upper=T) #马氏距离
1 2 3 4 5
1 0 8 7 3 2
2 8 0 5 5 6
3 7 5 0 6 7
4 3 5 6 0 1
5 2 6 7 1 0
> dist(X,method="minkowski",diag=T,upper=T,p=2)
1 2 3 4 5
1 0.000000 6.324555 5.385165 2.236068 1.414214
2 6.324555 0.000000 4.123106 4.123106 5.099020
3 5.385165 4.123106 0.000000 4.242641 5.000000
4 2.236068 4.123106 4.242641 0.000000 1.000000
5 1.414214 5.099020 5.000000 1.000000 0.000000

相似系数

对两个变量之间的相似程度可以用相似系数来刻画,用$C_{ij}$表示第i个变量与第j个变量之间的相似系数,$C_{ij}$的绝对值越接近1,表示指标i与指标j的关系越密切,$C_{ij}$的绝对值越接近于0,表示指标i与指标j的关系越疏远,常用的相似系数有夹角余弦和相关系数,如下所示:

夹角余弦:
相关系数:

相似系数之间的转换,一般来说,距离越小,两样本之间关系越密切,而相似系数越大,两变量之间的关系越密切,为了聚类分析方便起见,可以用下面的通用公式得到变量间的距离,即$d_{ij}^2=1-C_{ij}^2$。

系统聚类法

系数聚类法的基本思想

确定了距离和相似系数后就要进行分类。分类有许多种方法,最常用的一种方法是在样品距离的基础上定义类与类之间的距离。首选将n个样品分成n类,每个样品自成成一类,然后每将将具有最小距离的两类合并,合并后重新计算类与类之间的距离,这个过程一直持续到将所有的样品归为一类为止,并把这个过程画成一张聚类图,由聚类图可以方便地进行分类。因为聚类图很像一张系统图,所以这种方法就叫层次聚类法(hierachical clustering method)。层次聚类方法是目前应用最多的一种方法。从上面的分析可以看出,虽然已经给出了计算样品之间的距离,但在实际计算过程中还要定义类与类之间的距离。如何定义类与类之间的距离,也有很多方法,不同的方法就产生了不同的层次聚类,常用的方法有六种:

  1. 最短距离法:类与类之间的距离等于两类最近样品之间的距离;
  2. 最长距离法:类与类之间的距离等于两类最远样品之间的距离;
  3. 类平均法:类与类之间的距离等于各类元素两两之间平方距离的平均;
  4. 重心法:类与类之间的距离定义为对应这两类重心之间的距离。对样品分类来说,每一类的类重心就是该类样品的均值。
  5. 中间距离法:最长距离夸大了类间距离,最短距离低估了类间距离,介于两者之间的距离法即为中间距离法,中间距离法则是采用了介于最远和最近之间的距离。
  6. 离差平方和法(Ward法):基于方差分析的思想,如果类分得正确,同样样品之间的离差平方和应当较小,类与类之间的离差平方和应该较大。

层次聚类法的计算公式

最短距离法

该法使用$D_{k}(p, q)=\min \left\{d_{i j} | i \in G_{p}, j \in G_{q}\right\}$来刻画$G_p$与类$G_{1}$中最临近的两个样本的距离。若类$G_{p}$与类$G_{q}$合并为$G_{r}$,则$G_{r}$与$G_{s}$的距离为:

最长距离法

该法用$D_{k}(p, q)=\max \left\{d_{i j} | i \in G_{p}, j \in G_{q}\right\}$来刻画$G_p$与类$G_{1}$中最临近的两个样本的距离。

若类$G_{p}$与类$G_{q}$合并为$G_{r}$,则$G_{r}$与$G_{s}$的距离为:

重心法

在样本空间中,一个类用它的重心(即该类样本的均值)作代表较为合理,类与类之间的距离就用重心之间的距离来表示。

设样本之间的距离用欧氏距离,若若类$G_{p}$与类$G_{q}$合并为$G_{r}$,它们各有$n_{p}$、$n_{q}$、$n_{r}\left(n_{r}=n_{p}+n_{q}\right)$个样本,它们的重心用$\overline{x}_{P}$、$\overline{x}_{q}$和$\overline{x}_{r}$表示,显然,$\overline{x}_{r}=\frac{1}{n_{r}}\left(n_{p} \overline{x}_{P}+n_{q} \overline{x}_{q}\right)$,某一类$G_{k}$重心为$\overline{x}_{k}$,它与新类$G_{r}$的距离是:

其递推公式为:

中间距离法

该方法是对最短距离法和最长距离法的一个折中,即令类间距离的递推公式为:

类平均法

将两类之间的距离平方案底为这两类元素两两之间的平均来计算距离,即:

其递推公式为:

离差平方和法(Ward法)

这种方法是Ward提出来的,因此称为Ward法,该方法的基本思想来自于方差分析,如果分类正确,同样样品的离差平方和应当较小,类与类的离差平方和应当较大。具体做法是先将n个样品各自归为一类,然后每次缩小一类,每缩小一类,离差平方和就要增大,选择使方差增加最小的两类合并,直接所有的样本都归为一类为止。计算过程如下所示:

mark

mark

这六种层次聚类法的并类原则和过程完全相同,不同之处在于类与类之间的距离定义。当采用欧氏距离时,Lance和Williams于19678年将这些方法统一成如下的递推公式,即:

mark

mark

层次聚类法的基本步骤

  1. 计算n个样本两两的距离$\{d_{ij}\}$,记作D;
  2. 构造n个类,每个类只包含一个样品;
  3. 合并距离最近的两类为一个新类;
  4. 计算新类与当前各类的距离,若类个数为1,进行第5步,否则回到第3步;
  5. 画聚类图;
  6. 决定类的个数和类。

层次聚类的函数

R中进行层次聚类的函数是hclust(),用法如下所示:

1
2
3
hclust(d,mehtod="complete",...)
# d为相似矩阵
# method是层次聚类方法,包括ward,single,complete,average,mcquitty,median,centroid

现在使用案例A中的数据进行层次聚类。开始有5类,即每个样品自成一类,即$\{G_{1},G_{2},G_{3},G_{4},G_{5}\}$,这5类之间的距离就等于5个样品之间的距离,距离矩阵记为$D_{0}$,其中最小的元素为$D_{0}(4,5)=1$,因此将$G_{4}$和$G_{5}$聚合并成一个新类,即$G_{6}=\{G_{4},G_{5}\}$,然后再过计算$G_{6}$与$G_{1},G_{2},G_{3}#之间的距离。

应用公式$G_{1}(6,i)=min{D_{0}(4,i),D_{0}(5,i)}$,求其最近相邻的距离,如下所示:

mark

同理可知,其它类与类之间的距离如下所示:

mark

现在绘制层次聚类图,如下所示;

1
2
3
4
5
6
7
8
> hc <-hclust(dist(X),"single")
> cbind(hc$merge,hc$height)
[,1] [,2] [,3]
[1,] -4 -5 1.000000
[2,] -1 1 1.414214
[3,] -2 2 4.123106
[4,] -3 3 4.123106
> plot(hc)

mark

ward法(采用欧氏距离),如下所示:

1
2
3
4
5
6
7
8
9
hc2 <-hclust(dist(X),"ward")
cbind(hc2$merge,hc$height)
> cbind(hc2$merge,hc$height)
[,1] [,2] [,3]
[1,] -4 -5 1.000000
[2,] -1 1 1.414214
[3,] -2 -3 4.123106
[4,] 2 3 4.123106
plot(hc2)

mark

案例分析

为了研究我们31个省、市、自治区2001年城镇居民生活消费的分布规律,根据调查资料作区域消费类型划分,其中变量如下所示:

X1:人均食品支出;

X2:人均衣着商品支出

X3:人均家庭设备用品及服务支出

X4:人均医疗保健支出

X5:人均交通和通信支出

X6:人均娱乐教育文化服务支出

X7:人均居住支出

X8:人均杂项商品中和服务支出

表格如下所示(完整的数据在github上):

mark

现在使用六种聚类方法进行聚类,如下所示:

第1种:欧氏距离,euclidean

1
2
3
4
5
6
7
x_raw<-read.table("https://raw.githubusercontent.com/20170505a/raw_data/master/wbh_cluster_072.txt",header=T,sep=",")# 完整数据
rownames(x_raw)<-as.vector(x_raw$X)
x_raw<-x_raw[,-1]
x_dist1<-dist(x_raw,method="euclidean")
hc1 <-hclust(x_dist1,"single")
fviz_dend(hc1) #factoextra包

mark

还可以对分成几类进行颜色区分,如下所示:

1
2
3
4
5
fviz_dend(hc1, k = 4,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
rect = TRUE)

mark

第二种:欧氏距离,complete

1
2
hc2 <-hclust(x_dist1,"complete")
fviz_dend(hc2)

mark

第三种:欧氏距离,euclidean,median,如下所示:

1
2
hc3 <-hclust(x_dist1,"median")
fviz_dend(hc3)

mark

第四种:欧氏距离,euclidean,average,如下所示:

1
2
hc4 <-hclust(x_dist1,"average")
fviz_dend(hc4)

mark

第五种:欧氏距离,centroid,如下所示:

1
2
hc5 <-hclust(x_dist1,"centroid")
fviz_dend(hc5)

mark

第六种:欧氏距离,ward,如下所示:

1
2
hc6 <-hclust(x_dist1,"ward")
fviz_dend(hc6)

mark

划分类别,如下所示::

1
2
3
4
5
fviz_dend(hc6, k = 3,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800"),
color_labels_by_k = TRUE,
rect = TRUE)

mark

kmeans聚类法

kmeans聚类法的概念

层级聚类法需要计算出不同样本或变量的距离,还要在聚类的每一步计算出类间距离,相应的计算量比较大,尤其是当样本量比较大的时候,非常消耗计算机资源。而kmeans是一类快速聚类法。kmeans法(K均值法)的具体算法至少包括以下三个步骤:

  1. 将所有的样本分成k个初始类;
  2. 通过欧氏距离将某个样本划入离中心最近的类中,并对获得的样本与失去样品的类重新计算中心坐标
  3. 重复步骤2,直到所有的样品都不能再分配为止。

kmeans法和层级聚类法一样,都是以距离的远近亲疏为标准进行聚类的。但两者的不同之处在于:层级聚类对没的类数产生一系列的聚类结果,而K均值法只能产生指定类数的聚类结果。具体类数的确定,离不开实践经验的积累。有时候也要借助层级聚类法,以一部分样本为对象进行聚类,其结果为Kmeans法提供参考。

kemans聚类的原理与计算

kmeans算法以k为参数,把n个对象分为k个聚类,以使聚类内具有较高的相似度,而聚类间的相似度较低。相似度的计算是根据一个聚类中对象的均值来进行的。kmeans算法的处理流程如下所示:

  1. 随机地选择k个对象,每个对象初始地代表了一个簇的平均值或中心,对剩余的每个对象,根据其与各个聚类中心的距离将其赋给最近的簇;
  2. 重新过计算每个簇的平均值作为聚类中心进行聚类。

这个过程不断重复,直到准则函数收敛,通常,采用平方误差准则,其定义如下所示:

mark

其中,E为数据中所有对象与相应聚类中心的均方差之和,p为代表对象空间中的一个点,$m_{i}$为类$C_{i}$的均值(p和$m_{i}$都是多维的)。该式所示聚类的标准旨在使所有获得的聚类有以下特点:

各本身尽可能地紧凑,而各类之间尽可能地分开,kmean迭代图例如下所示:

mark

根据聚类中的均值进行聚类划分的kmeans算法如下:

  1. 从n个数据对象中取任意k个对象作为初始簇中心;
  2. 循环下述流程3到4,直到每个聚类不再发生变化为止;
  3. 根据每个簇中对象的均值(中心对象),计算每个对象与这些中心对象的距离,并根据最小距离重新对相应对象进行划分;
  4. 重新计算每个(有变化)簇的均值。

kmeans()

R中的kmeans()函数用于进行kmeans聚类,用法为:

1
2
3
kmeans(x,centers,..)
其中,x为数据矩阵或数据框
centers为聚类数或初始聚类中心

案例分析

kmeans算法的R语言实现及模拟。本例模拟正态随机变量$x~N(\mu,\delta^2)$。

  1. 用R模拟1000个均值为0,标准差为0.3的正态分布随机数据,再把这些随机数转换为10个变量、100个对象的矩阵;
  2. 用同样的方法模拟1000个均值为1,标准差为0.3的正态分布随机数,再转化为10个变量、100个对象的矩阵;
  3. 把这两个矩阵合并为10个变量、200个样本的数据矩阵;
  4. 利用kmeans聚类法将其聚成两类,观察其效果,如下所示:

先看一下层级聚类:

1
2

mark

再看一下kmeans聚类:

1
2
3
library(factoextra)
cl<-kmeans(x_norm,2)
fviz_cluster(cl,data=x_norm)

mark

我们通过fviz_nbclust函数查看一下聚成几类比较好,如下所示:

1
2
x_norm2<-scale(x_norm)
fviz_nbclust(x_norm2, kmeans, method = "wss") + geom_vline(xintercept = 4, linetype = 2)

mark

从图上可以看出,聚成2类比较好,因为1到2变化幅度最大,2之后的变化幅度就非常小了。

可视化二:可以将不同的聚类呈现出来:

使用silhouette函数来计算一下要绘图的轮廓。

1
2
3
x.res <- kmeans(x_norm, 2)
sil <- silhouette(km.res$cluster, dist(x_norm))
fviz_silhouette(sil)

mark

其它R语言聚类函数

eclust()是另外的一个聚类函数,此函数简化了聚类分析的工作流程,可以用于计算层次聚类和分区聚类,eclust()自动计算最佳聚类簇数。 自动提供Silhouette plot,可以结合ggplot2绘制图形,使用eclust()的K均值聚类如下所示:

1
2
x_eclust <- eclust(x_norm,"hclust")
fviz_dend(x_eclust,rect=TRUE)

mark

现在计算一下前面的那个案例,也就是各个省市的聚类,如下所示:

1
2
x_eclust2 <-eclust(x_raw,"hclust",k=4)
fviz_dend(x_eclust2,rect=TRUE)

mark

各省市的Kmean聚类图,如下所示:

1
2

代码汇总

参考资料

  1. 多元统计分析及R语言建模.王斌会

  2. R语言聚类分析-kmeans聚类分析实战

  3. R语言学习笔记之聚类分析