数量生态学：R语言的应用 第四章 聚类分析3—非层次聚类
在聚类分析中层次聚类被经常使用，层次聚类通过某种相似性测度计算节点之间的相似性，并按相似度由高到低排序，逐步重新连接个节点。而另一类非层次聚类是对一组对象进行简单分组的方法，尽量使组内对象之间比组间对象之间的相似度更高。但是需要我们自己来决定分组的组数k。
通常是设定不同的初始结构，然后通过大量的迭代以找到最佳的解决方法。本次介绍两种非层次聚类算法：k-均值划分和围绕中心点划分（PAM）。这些都是欧氏距离空间属性的算法。
1. 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-均值划分是一种线性模型的方法，因此不适合含很多零值的原始数据。遇到这样的数据，有两种方法进行处理。
为了与之前的Ward聚类结果比较，设定k=4组进行k-均值划分，并将其结果与ward结果进行比较。
#预转化后物种数据k-均值划分
spe.kmeans <- kmeans(spe.norm, centers = 4, nstart = 100)
spe.kmeans
# 比较当前分4组的分类结果与之前Ward聚类的结果。
table(spe.kmeans$cluster, spech.ward.g)
#这两个聚类结果是否非常相似？哪个（或哪些）对象有差别？
spe.kmeans$cluster
spech.ward.g
kmeans（）函数每次分析只能产生一个简单的预先设定组数的分组结果。如果需要尝试不同组数的结果，需要重新运行命令。但到底多少组是最好的分类方案呢？
最好方案的标准有很多，cclust程序包clustIndex（）函数内有一部分标准。推荐使用最大化Calinski-Harabasz指数（比较分类结果组间与组内平方和的F统计量），尽管在分类组所含对象数量差异比较大情况下Calinski-Harabasz指数的值可能比较低。“ssi”（简单结构指数）也是衡量分类结果的另一个不错的指标。
以下使用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)
到底多少组数是最佳方案？如果倾向于较大的组数，哪个是最佳方案呢？
函数casecadeKM()也提供数据结果。其中“result”给出每种组数k条件下的TESS统计量和准则(calinski或ssi)的值，其中partition包含每个聚类簇所含对象的信息。如果对象的地理信息可用，就可以绘制对象空间分布地图，同时以不同的颜色带标识对象的归属。
summary(spe.KM.cascade)
spe.KM.cascade$results
确定样方聚类簇后，可以查看聚类簇的内容。最简单的方法是基于分组图形定义样方的亚组合计算基本统计量。以下对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组中各个物种的平均多度矩阵
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)
startobjects <- spe.norm[c(2, 17, 21, 23), ]
spe.kmeans3 <- kmeans(spe.norm, centers = startobjects)
#将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(1, 1))
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组分布图")
2. 围绕中心点划分（PAM）
围绕中心点划分：“从所有的数据观测点寻找k个代表性的对象或形心点，这些代表性的对象应该反映数据的主体结构。k个形心点选定后，将每个观测点分配给某个形心点构建k个聚类簇，不断寻找最佳的k个代表性对象，使对象之间的相异性总和最小”。k-均值法是最小化组内欧氏距离平方和。因此，k-均值法是传统的最小二乘法，但PAM不是。
以下代码两的最后部分利用双轮廓图比较k-均值法和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)
# 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)
#k=4组下与Ward聚类和k-均值划分进行比较
dev.new(
title = "轮廓宽度值 - k-均值 and PAM",
width = 12,
height = 8,
noRStudioGD = TRUE
)
par(mfrow = c(1, 2))
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)
)
基于此图，可以分辨PAM和k-均值法哪个有更好的平均轮廓值。我们可以将当前结果与Ward聚类比较轮廓宽度图比较。除了聚类簇成员不同外，在k-均值和最优化的Ward聚类之间还有哪些不同?我们可以通过检查两个聚类比较的列联表进行分析
#比较k-均值划分和最优化后的Ward聚类
table(spe.kmeans.g, spech.ward.gk)
从前面的例子，可以看出即使是两种目标相同且属于同一大类（k-均值法和PAM同属于非层次聚类）的聚类方法也可能产生完全不同的结果。我们需要根据哪种方法的结果能够提供更多有用的信息，或能更好地被环境变量解释来选择合适的方法。
好了，第四章聚类分析的非层次聚类就这么多内容，主要就是k-均值法和PAM这两种方法，如果还不能好好理解，可以多回顾几次，下次内容是聚类分析中的使用环境变量来解释来比较，敬请期待。
谢谢你的阅读。
《数量生态学：R语言的应用》第二版第三章-关联测度与矩阵------Q模式
