生信代码:层次聚类和K均值聚类

2021-01-25 10:17:38 浏览数 (1)

1. 层次聚类

层次聚类 (hierarchical clustering)是一种对高维数据进行可视化的常见方法。

层次聚类常用方法是聚合法 (agglomerative approach),它是一种自下而上的方法,把数据当做一些独立的点,计算数据点之间的距离,然后按照一定的合并策略,先找出数据集中最近的两点,把它们合并到一起看作一个新的点,重复这个过程,得到一棵数据树——树状图 (dendrogram),展示数据聚类结果。

➢距离度量

1. 欧氏距离(euclidean distance):两点之间的直线距离。

2.相似距离(correlation similarity):两点之间的相关性,相似度。

3.曼哈顿距离 (Manhattan distance):两点在标准坐标系上的轴距离之差的绝对值的和。

i和j代表第i和第j个观测值,p是维度。

➢层次聚类的合并策略

・Average Linkage聚类法:计算两个簇中的每个数据点与其他簇的所有数据点的距离。将所有距离的均值作为两个簇数据点间的距离。

・Complete Linkage聚类法:两个簇中距离最远的两个点间的距离作为这两个簇的距离。

例:

生成一个数据集,text( )给每个点加标签

代码语言:javascript复制
set.seed(1234) 
par(mar = c(0, 0, 0, 0))
x <- rnorm(12, mean = rep(1:3, each = 4), sd = 0.2)
y <- rnorm(12, mean = rep(c(1, 2, 1), each = 4), sd = 0.2)
plot(x, y, col = "blue", pch = 19, cex = 2)
text(x   0.05, y   0.05, labels = as.character(1:12))

dist( )计算数据框中不同⾏所表示的观测值之间的距离,返回距离矩阵 (distance matrix),默认计算欧⽒距离。

hclust( )进行聚类分析

代码语言:javascript复制
distxy <- dist(dataFrame)
hClustering <- hclust(distxy)
par(mar = c(6, 4, 2, 0))
plot(hClustering)

聚类分析返回树状图,展示了这些点是怎样聚成簇的,但没有说明一共有多少个簇。

在y=2.0的位置截断,会碰到两个分支,表明大致有两个簇;在y=1.0的位置截断,碰到三个分支,这说明有三个簇。目前没有规则确定要从哪儿截断,一旦在某个位置截断,就可以从层次聚类中得到各个簇的情况,必须截断在合适的位置。

myplclust函数
代码语言:javascript复制
myplclust <- function(hclust, lab = hclust$labels, lab.col = rep(1, length(hclust$labels)), 
                      hang = 0.1, ...) {
  y <- rep(hclust$height, 2)
  x <- as.numeric(hclust$merge)
  y <- y[which(x < 0)]
  x <- x[which(x < 0)]
  x <- abs(x)
  y <- y[order(x)]
  x <- x[order(x)]
  plot(hclust, labels = FALSE, hang = hang, ...)
  text(x = x, y = y[hclust$order] - (max(hclust$height) * hang), labels = lab[hclust$order],
       col = lab.col[hclust$order], srt = 90, adj = c(1, 0.5), xpd = NA, ...)
}

dataFrame <- data.frame(x = x, y = y)
distxy <- dist(dataFrame)
hClustering <- hclust(distxy)
myplclust(hClustering, lab = rep(1:3, each = 4), lab.col = rep(1:3, each = 4))

myplclust( )输出一个聚类树状图,每个簇里边的所有点都会由它们所在簇的标签来标记,并且会由不同的颜色来表现。注意,需要在实际标注不同颜色的"1" "2" "3" 之前指明一共有多少类。

heatmap函数
代码语言:javascript复制
dataFrame <- data.frame(x = x, y = y) 
set.seed(143)
dataMatrix <- as.matrix(dataFrame)[sample(1:12), ] 
heatmap(dataMatrix)

heatmap( )对行进行聚类分析,将列看作为观测值,生成热图,根据层次聚类算法对表格中的行和列进行重排。行的左侧有一个聚类树状图,说明可能存在三个簇。

2. K均值聚类

K均值聚类 (K-means clustering)是一种迭代求解的聚类分析算法,可以用于整理高维数据,了解数据的规律,寻找最佳的数据模式,但前提需要确定簇的数量(肉眼判断,交叉验证,信息理论等方法),因此需要进行多次尝试计算,选择最佳的结果。

➢基本方法

确定将数据分为K组,随机选取K个几何中心(centroid),计算每个数据点到这些几何中心的距离,把所有点分配给距离它最近的中心,然后重新计算每一簇的几何中心,再重新分配所有点,反复操作直到K均值聚类算法得到一个对于几何中心位置的最终估计并说明每个观测值分配到哪一个几何中心。

如果运行了3次K均值算法,每次得到的模式都不同,那就表示这个算法或许不能对这个数据产生稳定的判断,因此K均值用在这一类的数据集上可能是有问题的。

以上文使用的数据集为例,选取3个随机的点作为几何中心

读取数据点分配给最近的几何中心,重新计算几何中心,如通过计算这个簇的平均值,重新读取数据点分配给最近的几何中心。

重复以上计算几何中心及分配数据点的过程,直到得到一个对于几何中心位置的最终估计。

kmeans( )执行K均值算法

代码语言:javascript复制
dataFrame <- data.frame(x,y)
kmeansObj <- kmeans(dataFrame,centers=3) 
#对 dataFrame 调用 kmeans() 并说明有3个几何中心,返回一个列表
names(kmeansObj)
[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"   
[7] "size"         "iter"         "ifault"

kmeansObj$cluster #指明数据点属于哪个簇
[1] 3 3 3 3 2 2 2 2 1 1 1 1

kmeansObj$centers #几何中心在空间中的位置
          x         y
1 2.8534966 0.9831222
2 1.9906904 2.0078229
3 0.8904553 1.0068707

绘制k均值聚类结果

代码语言:javascript复制
par(mar=rep(0.2,4)) 
plot(x,y,col=kmeansObj$cluster,pch=19,cex=2) 
#绘制数据,col 参数等于簇的数量,根据数据点所在的簇上色
points(kmeansObj$centers,col=1:3,pch=3,cex=3,lwd=3)
#添加几何中心

set.seed(1234)
dataMatrix <- as.matrix(dataFrame)[sample(1:12),] 
kmeansObj <- kmeans(dataMatrix,centers=3) 
par(mfrow=c(1,2), mar = c(2, 2, 2, 2)) 
image(t(dataMatrix)[,nrow(dataMatrix):1],yaxt="n") 
image(t(dataMatrix)[,order(kmeansObj$cluster)],yaxt="n")

图左为原始数据的热图,图右为按照簇进行归类重新排列过的数据。

0 人点赞