Loading...
墨滴

fafu生信小蘑菇

2021/05/09  阅读:48  主题:橙心

数量生态学:R语言的应用 第四章 聚类分析3—非层次聚类

数量生态学:R语言的应用 第四章 聚类分析3—非层次聚类

在聚类分析中层次聚类被经常使用,层次聚类通过某种相似性测度计算节点之间的相似性,并按相似度由高到低排序,逐步重新连接个节点。而另一类非层次聚类是对一组对象进行简单分组的方法,尽量使组内对象之间比组间对象之间的相似度更高。但是需要我们自己来决定分组的组数k。

通常是设定不同的初始结构,然后通过大量的迭代以找到最佳的解决方法。本次介绍两种非层次聚类算法:k-均值划分和围绕中心点划分(PAM)。这些都是欧氏距离空间属性的算法。

注意:不同纲量的变量在进行非层次聚类之前应该进行标准化,因为方差的纲量等于变量纲量的平方,如果不标准化,不同量纲的方差相加无意义。

1. k-均值划分

k-均值划分使用数据局部结构构建聚类簇:通过确认数据高密度区构建分类组。

其步骤是:预估将数据分为K组,则随机选取K个对象作为初始的聚类中心,然后计算每个对象与各个种子聚类中心之间的距离,把每个对象分配给距离它最近的聚类中心 ,聚类中心以及分配给它们的对象就代表一个聚类 ,每分配一个样本,聚类的聚类中心会根据聚类中现有的对象被重新计算,这个过程将不断重复直到满足某个终止条件 ,终止条件可以是没有(或最小数目)对象被重新分配给不同的聚类,没有(或最小数目)聚类中心再发生变化,误差平方和局部最小。

k-均值划分算法以欧式距离作为相似度测度,采用误差平方和准则函数作为聚类准则函数。

加载所需的包和数据

#加载包和数据
library(ade4)
library(adespatial)
library(vegan)
library(gclus)
library(cluster)
library(pvclust)
library(RColorBrewer)
library(labdsv)
library(rioja)
library(indicspecies)
library(mvpart)
library(MVPARTwrap)
library(dendextend)
library(vegclust)
library(colorspace)
library(agricolae)
library(picante)

#加载函数
source("drawmap.R")
source("drawmap3.R")
source("hcoplot.R")
source("test.a.R")
source("coldiss.R")
source("bartlett.perm.R")
source("boxplerk.R")
source("boxplert.R")

#从聚类结果获得二元差异矩阵的函数
grpdist <- function(X)
{
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr,"gower")
  distgr
}

#导入Doubs数据
load("Doubs.RData")
#剔除无物种数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
latlong <- latlong[-8,] #经纬度

#先计算样方之间的弦距离矩阵
spe.norm <- decostand(spe,"normalize")

1.1 随机开始的K-均值划分

在已经想好设定的组数后,使用stats包的 kmeans()函数运行k-均值划分。kmeans()函数会以随机设定的初始结构为基础,不断的自动重复迭代到设定次数(参数nstart)为止,并获得当前次数下最佳的分类方案(即最小SSE值)。

k-均值划分是一种线性模型的方法,因此不适合含很多零值的原始数据。遇到这样的数据,有两种方法进行处理。

  • 使用非欧氏相异矩阵,如百分数差异(又名Bray-Curtis);这样的非欧氏相异矩阵应进行平方根变换并提交主坐标分析(PCoA)以获得欧氏属性的主坐标用于k-均值划分。
  • 对数据进行预转化。为了与前面聚类分析的对象保持一致,可以使用之前已经范数标准化后(弦转化)的物种数据作为案例。为了与之前的Ward聚类结果比较,设定k=4组进行k-均值划分,并将结果与Ward层次聚类4组样方的结果进行比较。

为了与之前的Ward聚类结果比较,设定k=4组进行k-均值划分,并将其结果与ward结果进行比较。

#预转化后物种数据k-均值划分
spe.kmeans <- kmeans(spe.norm, centers = 4, nstart = 100)
spe.kmeans
预转化后物种数据k-均值划分
预转化后物种数据k-均值划分

注意:即使给定的nstart相同,每次运行上述命令,所产生的结果也不一定完全相同,因为每次运算设定的初始结构是随机的。

# 比较当前分4组的分类结果与之前Ward聚类的结果。
table(spe.kmeans$cluster, spech.ward.g)

#这两个聚类结果是否非常相似?哪个(或哪些)对象有差别?
spe.kmeans$cluster
spech.ward.g
比较当前分4组的分类结果与之前Ward聚类的结果
比较当前分4组的分类结果与之前Ward聚类的结果
比较当前分4组的分类结果与之前Ward聚类的结果
比较当前分4组的分类结果与之前Ward聚类的结果

对于那些无法从原始数据转化后进行欧氏距离计算的距离指数(例如百分数差异指数,又名Bray-Curtis相异指数),必须先将原始的距离矩阵通过主坐标分析(PCoA)获得一个n行的矩形数据表格,然后再进行k-均值划分。

对于百分数差异相异指数,必须先取平方根进行PCoA分析,或使用校对负特征根的PCoA函数,才能获得完全的欧氏属性数据表格。

kmeans()函数每次分析只能产生一个简单的预先设定组数的分组结果。如果需要尝试不同组数的结果,需要重新运行命令。但到底多少组是最好的分类方案呢?

最好方案的标准有很多,cclust程序包clustIndex()函数内有一部分标准。推荐使用最大化Calinski-Harabasz指数(比较分类结果组间与组内平方和的F统计量),尽管在分类组所含对象数量差异比较大情况下Calinski-Harabasz指数的值可能比较低。“ssi”(简单结构指数)也是衡量分类结果的另一个不错的指标。

我们不必手动多次运行kmeans()函数来获得不同组数的分类结果,可以使用vegan包里的函数cascadeKM()。cascadeKM()函数可以一次运行生成组数k从少(参数Inf.gr)到多(参数sup.gr)的聚类结果。

以下使用cascadeKM()函数,组数设定从2组到10组,并使用ssi评估聚类的质量,同时绘制聚类的结果。

#k-均值划分,2组到10组
spe.KM.cascade <-
  cascadeKM(
    spe.norm,
    inf.gr = 2,
    sup.gr = 10,
    iter = 100,
    criterion = "ssi"
  )
summary(spe.KM.cascade)
spe.KM.cascade$results
spe.KM.cascade$partition
dev.new(
  title = "CascadeKM",
  width = 10,
  height = 6,
  noRStudioGD = TRUE
)
#在plot函数中,sortg=TRUE代表在每个聚类簇内按照对象之间的紧密程度重新排列对象。
plot(spe.KM.cascade, sortg = TRUE)
显示不同组数条件下每个对象归属的k-均值划分
显示不同组数条件下每个对象归属的k-均值划分

该图显示每个对象在每种分类组数下的归属(图上每行代表一种组数)。图内的表格有不同的颜色,每行两种颜色,代表分两组k=2,三种颜色代表k=3,依此类推。右图代表不同k值条件下的中止标准的统计量。

到底多少组数是最佳方案?如果倾向于较大的组数,哪个是最佳方案呢?

函数casecadeKM()也提供数据结果。其中“result”给出每种组数k条件下的TESS统计量和准则(calinski或ssi)的值,其中partition包含每个聚类簇所含对象的信息。如果对象的地理信息可用,就可以绘制对象空间分布地图,同时以不同的颜色带标识对象的归属。

summary(spe.KM.cascade)
spe.KM.cascade$results
image-20210508172756428
image-20210508172756428

这里,最小SSE值用于确定在给定组数k下的最佳方案,而calinski和ssi指标用于确定最佳k值,两个指标解决不同的问题。

记住不同的组数k={2,3,…,10},每次都是独立运行。因此上图从下到上结构互相独立,与层次聚类树的嵌套结构不同。

确定样方聚类簇后,可以查看聚类簇的内容。最简单的方法是基于分组图形定义样方的亚组合计算基本统计量。以下对4组k-均值划分结果进行分析。

# 按照k-均值划分结果重新排列样方
spe.kmeans.g <- spe.kmeans$cluster
spe[order(spe.kmeans.g), ]

# 使用函数vegemite()重新排列样方-物种矩阵
ord.KM <- vegemite(spe, spe.kmeans.g)
spe[ord.KM$sites, ord.KM$species]

1.2 使用k-均值划分优化一个独立获得的分类

k-均值划分的出发点可以以k x p的平均值矩阵(k为聚类获得的组数,p为变量数量,矩阵里面的值为该变量在该组中的平均值),或是以k个对象(每个对象视为每个组的典型代表)作为列表,这个典型的代表即下一节要讨论的分组方法所称的“形心(medeids)”。

让我们应用这个想法来验证由鱼类数据Ward聚类获得的4个组,然后被k-均值修改分组。

以Ward 聚类获得4组中各个物种的平均多度矩阵作为初始结构进行k-均值划分

# 计算Ward聚类获得4组中各个物种的平均多度矩阵
groups <- as.factor(spech.ward.g)
spe.means <- matrix(0, ncol(spe), length(levels(groups)))
row.names(spe.means) <- colnames(spe)
for (i in 1:ncol(spe)) {
  spe.means[i, ] <- tapply(spe.norm[, i], spech.ward.g, mean)
}
# 平均物种多度矩阵作为初始点
startpoints <- t(spe.means)
# 基于初始点的k-均值划分
spe.kmeans2 <- kmeans(spe.norm, centers = startpoints)
image-20210509110909213
image-20210509110909213

与前面的方法有点不同,这里需要返回看一下聚类(层次结构聚类)的结果,确认每组最典型的对象(轮廓宽度值最大),将这些对象作为k-均值划分的初始结构(参数centers)

startobjects <- spe.norm[c(2172123), ]
spe.kmeans3 <- kmeans(spe.norm, centers = startobjects)
image-20210509112421431
image-20210509112421431
#将k-均值划分结果与Ward聚类分析结果进行比对
table(spe.kmeans2$cluster, spech.ward.g)

# 两个k-均值划分优化后的4组进行比较
table(spe.kmeans2$cluster, spe.kmeans3$cluster)

聚类分析结果进行比对
聚类分析结果进行比对
# 最终分组的轮廓宽度值图
spech.ward.gk <- spe.kmeans2$cluster
dev.new(
  title = "轮廓-优化分区",
  noRStudioGD = TRUE
)
par(mfrow = c(11))
k <- 4
sil <- silhouette(spech.ward.gk, spe.ch)
rownames(sil) <- row.names(spe)
plot(sil,
     main = "轮廓宽度值 - Ward & k均值",
     cex.names = 0.8,
     col = 2:(k + 1))
最终分组的轮廓宽度值图
最终分组的轮廓宽度值图

从代码中可以看出,可以通过列联表法进行原始分组和优化分组之间或两个优化分组之间的比较。两种优化(基于物种的平均值和典型的初始对象)为我们选择的4个“典型”对象产生相同的结果。注意这个选择至关重要。与原来的Ward聚类4组结果比较显示,只有一个对象第15样方已从组2移动到组1。这样移动可以改进轮廓图并略微改变了沿着河流分布的地图。请注意,样方19的成员归属仍不清楚:它的轮廓宽度值仍为负值。

# 沿DOUBS河流优化后Ward聚类4组分布图
dev.new(title = "河流优化后Ward聚类4组分布图",
        width = 9,
        noRStudioGD = TRUE)
drawmap(xy = spa,
       clusters = spech.ward.gk,
       main = "沿DOUBS河流优化后Ward聚类4组分布图")
沿DOUBS河流优化后Ward聚类4组分布图
沿DOUBS河流优化后Ward聚类4组分布图

2. 围绕中心点划分(PAM)

围绕中心点划分:“从所有的数据观测点寻找k个代表性的对象或形心点,这些代表性的对象应该反映数据的主体结构。k个形心点选定后,将每个观测点分配给某个形心点构建k个聚类簇,不断寻找最佳的k个代表性对象,使对象之间的相异性总和最小”。k-均值法是最小化组内欧氏距离平方和。因此,k-均值法是传统的最小二乘法,但PAM不是

在R里,pam()函数(cluster程序包)的输入可以是原始数据,也可以是相异矩阵(与kmean()相比,pam()的优势是可以输入更多类型的关联测度),并且允许通过轮廓宽度值确定最佳的分组数量。

以下代码两的最后部分利用双轮廓图比较k-均值法和PAM法的分组结果。

基于弦距离的围绕中心点划分(PAM)

#聚类簇数量的选择
#循环计算2至28组分类数下平均轮廓宽度
asw <- numeric(nrow(spe))
for (k in 2:(nrow(spe) - 1))
  asw[k] <- pam(spe.ch, k, diss = TRUE)$silinfo$avg.width
k.best <- which.max(asw)
dev.new(title = "PAM", noRStudioGD = TRUE)
plot(
  1:nrow(spe),
  asw,
  type = "h",
  main = "聚类簇数量的选择",
  xlab = "k (群集数)",
  ylab = "平均轮廓宽度"
)
axis(
  1,
  k.best,
  paste("optimum", k.best, sep = "\n"),
  col = "red",
  font = 2,
  col.axis = "red"
)
points(k.best,
       max(asw),
       pch = 16,
       col = "red",
       cex = 1.5)
聚类簇数量的选择
聚类簇数量的选择

当k=2时,PAM具有最佳的方案(asw=0.3841),这并不是我们期望的结果。如果选择常用的4组分法,从轮廓宽度角度分析结果并不好(asw=0.2736)。尽管如此,我们还是需要分析PAM分4组的情况。

# PAM分4组情况
spe.ch.pam <- pam(spe.ch, k = 4, diss = TRUE)
summary(spe.ch.pam)
spe.ch.pam.g <- spe.ch.pam$clustering
spe.ch.pam$silinfo$widths

# 将当前的分类结果与之前Ward聚类和k-均值划分进行比较
table(spe.ch.pam.g, spech.ward.g)
table(spe.ch.pam.g, spe.kmeans.g)
当前的分类结果与之前Ward聚类和k-均值划分进行比较
当前的分类结果与之前Ward聚类和k-均值划分进行比较

PAM结果与Ward聚类和k-均值划分结果显著不同

#k=4组下与Ward聚类和k-均值划分进行比较
dev.new(
  title = "轮廓宽度值 - k-均值 and PAM",
  width = 12,
  height = 8,
  noRStudioGD = TRUE
)
par(mfrow = c(12))
k <- 4
sil <- silhouette(spe.kmeans.g, spe.ch)
rownames(sil) <- row.names(spe)
plot(sil,
     main = "轮廓宽度值 - k-均值",
     cex.names = 0.8,
     col = 2:(k + 1))
plot(
  silhouette(spe.ch.pam),
  main = "轮廓宽度值 - PAM",
  cex.names = 0.8,
  col = 2:(k + 1)
)
轮廓宽度值 - k-均值 and PAM
轮廓宽度值 - k-均值 and PAM

基于此图,可以分辨PAM和k-均值法哪个有更好的平均轮廓值。我们可以将当前结果与Ward聚类比较轮廓宽度图比较。除了聚类簇成员不同外,在k-均值和最优化的Ward聚类之间还有哪些不同?我们可以通过检查两个聚类比较的列联表进行分析

#比较k-均值划分和最优化后的Ward聚类
table(spe.kmeans.g, spech.ward.gk)

比较k-均值划分和最优化后的Ward聚类
比较k-均值划分和最优化后的Ward聚类

PAM法表现更稳定,一方面因为它是最小化相异系数总和,代替了最小化欧氏距离的平方和;另一方面在给定分组数k的条件下,PAM相异系数总和更容易收敛。但对于特定的研究目的,并不能保证PAM一定是最合适的聚类方法。

从前面的例子,可以看出即使是两种目标相同且属于同一大类(k-均值法和PAM同属于非层次聚类)的聚类方法也可能产生完全不同的结果。我们需要根据哪种方法的结果能够提供更多有用的信息,或能更好地被环境变量解释来选择合适的方法

好了,第四章聚类分析的非层次聚类就这么多内容,主要就是k-均值法和PAM这两种方法,如果还不能好好理解,可以多回顾几次,下次内容是聚类分析中的使用环境变量来解释来比较,敬请期待。

谢谢你的阅读。

如有不足或错误之处,请批评指正。 有什么不明白的也欢迎留言讨论。

欢迎关注微信公众号:fafu 生信 小蘑菇

fafu 生信 小蘑菇
fafu 生信 小蘑菇

往期内容:

《数量生态学:R语言的应用》第三章-R模式

《数量生态学:R语言的应用》第二版第三章-关联测度与矩阵------Q模式

《数量生态学:R语言的应用》第二版笔记2

《数量生态学——R语言的应用》第二版阅读笔记--绪论和第二章(一部分)

R语言 pheatmap 包绘制热图(基础部分)

R语言pheatmap包绘制热图进阶教程

使用PicGo和gitee搭建图床

组间分析—T检验、R语言绘图

Rmarkdown的xaringan包来制作PPT

htlm文件部署到个人网站

感谢你的阅读!!!你的点赞关注转发是对我最大的鼓励。

fafu生信小蘑菇

2021/05/09  阅读:48  主题:橙心

作者介绍

fafu生信小蘑菇