数量生态学：R语言的应用 第四章 聚类分析8—顺序聚类和模糊聚类
今天是数量生态学：R语言的应用 第四章 聚类分析最后一部分内容了，不知不觉已经过了一周了，内容还是很多的，不知道你掌握了多少，不会的话的好好理理，不要着急，一定要稳扎稳打，消化吸收。
今天主要是将两个部分内容：
- 顺序聚类
- 模糊聚类
1.加载数据和包及数据预处理
# 1.加载数据和包及数据预处理#####
setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
library(devtools)
library(mvpart)
library(MVPARTwrap)
library(vegan)
#导入Doubs数据
load("Doubs.RData")
#剔除无物种数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
latlong <- latlong[-8,] #经纬度
##1.2 数据预处理#####
#先计算样方之间的弦距离矩阵
spe.norm <- decostand(spe,"normalize")
spe.ch <- vegdist(spe.norm,"euc")
# 计算Ward最小方差聚类
spe.ch.ward <- hclust(spe.ch, method="ward.D2")
#Ward聚类同表型相关
spe.ch.ward.coph <- cophenetic(spe.ch.ward)
cor(spe.ch, spe.ch.ward.coph)
# 设定聚类组的数量
k <- 4
# 根据上面4个融合水平值图，可以观察到分4组水平在所有图里有小的跳跃
# 裁剪聚类树
spech.ward.g <- cutree(spe.ch.ward, k=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)
# 最终分组的轮廓宽度值图
spech.ward.gk <- spe.kmeans2$cluster
2. 顺序聚类
在某些情况下，数据本身具有空间和时间系列属性，此时对数据进行分组需要考虑数据之间的连续性。顺序聚类也可以由MRT计算，解释变量是代表采样序列变量。
对于Doubs数据的29个样方，包含数字1到29的向量（等价于变量dfs）是合适作为解释变量（约束条件）。
CONISS是用于地层研究具有邻接约束的聚类方法。这种方法的原理与Ward的最小方差聚类类似，分割的节点是空间或时间系列。CONISS可以通过rioja程序包的函数chclust（）计算。
将使用鱼类多度数据的百分数差异（又称为Bray-Curtis）矩阵。
通过比较层次分类的离差（平方和，SS）和断棍模型的离差选择聚类的数量（图1）：
library(rioja)
#基于百分数差异矩阵的CONISS聚类
spe.chcl <- chclust(vegdist(spe))
#比较层次分类的离差（平方和，SS）和断棍模型的离差
bstick(spe.chcl,10)
这个图建议分2组或4组合适（红线是从断棍模型获得，在分2组或4组的地方平方和比断棍模型高）
这里我们选择分4组：
# 选择分成4组
k <- 4
(gr4 <- cutree(spe.chcl, k = k))
# 绘制聚类树
dev.new(
title = "Constrained (sequential) clustering",
width = 12,
height = 6,
noRStudioGD = TRUE
)
plot(spe.chcl, hang = -1, main = "CONISS 聚类")
rect.hclust(spe.chcl, k = k)
如果有必要，可以按照离源头距离（dfs）这个向量去重新分配聚类树的分支
可以看到下图，你可以发现样方15-22与样方26-30中间被样本23-25强行分开，主要是因为样方23-25是污染严重的样方，鱼类组成与其他样方明显有差异（图4）。
# Doubs 河流展示聚类的结果
source("drawmap.R")
dev.new(title = "河流上分4组情况",
width = 9,
noRStudioGD = TRUE)
drawmap(xy = spa,
clusters = gr4,
main = "沿着Doubs河流的顺序聚类分4组的情况")
3. 模糊聚类
之前我们所讨论的聚类方法所产生的聚类簇是非重叠的实体，聚类方法注重数据间断分布的自然结果。然而，还有模糊聚类。在模糊聚类中，一个对象可以不同程度归属于两个组或多个组。打个比方，我们可以通过混合黄色和蓝色的涂料获得绿色，而且根据黄色和蓝色不同的比例可以调出不同程度的绿色。所以一种不同于层次和非层次法被开发出来，即“模糊聚类”。这里简单介绍一种与非层次k-均值划分类似的模糊聚类，也称为c-均值聚类。
举一个感官评价的例子，假如啤酒口味可以被品酒者判断为甜（sweet）和苦（bitter），评价同一种啤酒，30%认为是甜，70%认为是苦，这个就是啤酒模糊分类的案例。
3.1 使用clugter程序包内fanny（）函数进行c-均值模糊聚类
c-均值聚类结果：一个对象可以赋予不同的组，对象与组之间的归属程度可以通过成员值衡量。一个对象在某一组内的成员值越高，表示该对象与该组之间关系越紧密，否则关系疏松。每个对象的成员值总和为1。
运行c-均值模糊聚类的包两个：cluster包（fanny（）函数）和e1071包（cmeans（）函数）。
下面使用fanny（）函数进行c-均值聚类。
函数fanny（）输入对象可以是样方-物种矩阵或距离矩阵。对于样方-物种数据，默认的是欧氏度量。
下面直接输入前面算好的弦距离矩阵进行分析。其结果与直接输入弦转化后的物种数据“spe.ch”然后选择参数metric=
"euclidean"相同。
函数fanny（）的输出结果是两个图：聚类簇的排序图（得在之后第5章学习）和轮廓图。此处只显示轮廓图，同时将主坐标分析（PCoA）结合对象的星状图代替原始的排序图。图5中一颗小星代表一个对象，其半径与成员系数呈正比。
library(cluster)
k <- 4 # 选择聚类分组数目
spe.fuz <- fanny(spe.ch, k = k, memb.exp = 1.5)
summary(spe.fuz)
# 样方成员值
spe.fuz$membership
# 每个样方最接近的聚类簇
spe.fuz$clustering
spefuz.g <- spe.fuz$clustering
# 轮廓图
dev.new(title = "Fuzzy clustering of fish data - Silhouette plot",
noRStudioGD = TRUE)
plot(
silhouette(spe.fuz),
main = "轮廓图-模糊聚类",
cex.names = 0.8,
col = spe.fuz$silinfo$widths + 1
)
# 模糊聚类簇的主坐标 (PCoA)
# 第一步：鱼类数据弦距离矩阵的主坐标分析
dc.pcoa <- cmdscale(spe.ch)
dc.scores <- scores(dc.pcoa, choices = c(1, 2))
# 第二步：模糊聚类簇的主坐标排序
dev.new(title = "鱼类数据排序图的模糊聚类",
noRStudioGD = TRUE)
plot(dc.scores,
asp = 1,
type = "n",
main = "模糊聚类排序 (PCoA)")
abline(h = 0, lty = "dotted") #画水平和垂直线
abline(v = 0, lty = "dotted")
# 第三步：模糊聚类的体现
for (i in 1:k)
{
gg <- dc.scores[spefuz.g == i, ] #dc.scores主坐标分析数据 #spefuz.g最接近的聚类簇
hpts <- chull(gg)
hpts <- c(hpts, hpts[1])
lines(gg[hpts, ], col = i + 1)
}
stars(
spe.fuz$membership, #spe.fuz 聚类结果
location = dc.scores, #dc.scores主坐标分析数据
key.loc = c(0.6, 0.4),
key.labels = 1:k,
draw.segments = TRUE,
add = TRUE,
# scale = FALSE,
len = 0.075,
col.segments = 2:(k + 1)
)
参数memb.exp称为“模糊指数"，其值从1（接近非模糊聚类）到任何大于1的值。
图5轮廓图显示第2聚类簇分组并不理想。此外，排序图也显示（样方10、15和19）的归属不清晰。
屏幕上显示的数量结果提供每个对象的成员值。每行的和等于1。明确属于某一聚类簇的对象具有较高的成员值，例如样方2、21和23，因此它们在其他聚类簇内的成员值也相对低。相反，对于容易归类的样方，其成员值比较均匀，例如样方5、9和19。另一个结果显示每个对象最接近的聚类簇，即每个对象具有最高成员值所对应的聚类簇。
接下来可以把代表模糊隶属度的扇区添加到Doubs河的地图上（图7）：
# Doubs 河地图绘制模糊聚类的结果
dev.new(title = "河流上的四个模糊聚类",
width = 9,
noRStudioGD = TRUE)
plot(
spa,
asp = 1,
type = "n",
main = "沿着河流的模糊聚类",
xlab = "x 坐标 (km)",
ylab = "y 坐标 (km)"
)
lines(spa, col = "light blue")
text(65, 20, "上游", cex = 1.2)
text(15, 32, "下游", cex = 1.2)
#添加扇形区以表示模糊聚类成员资格
for (i in 1:k) {
stars(
spe.fuz$membership,
location = spa,
key.loc = c(150, 20),
key.labels = 1:k,
draw.segments = TRUE,
add = TRUE,
# scale = FALSE,
len = 5,
col.segments = 2:(k + 1)
)
}
3.2 使用vegclust（）函数进行噪声聚类
vegclust包能提供一个大范围的选项来执行不同模型下的群落数据非层次或层次的模糊聚类。这种聚类的算法被称为“噪声聚类”。该方法试图使模糊聚类的异常值有更明确的意义。
异常值：一旦定义了“真实聚类”形心，捕获远离“真实聚类”形心距离为δ的虚拟“噪声”聚类的对象。
δ值的选择要点：δ值太小会导致过大的异常值，即“噪声”簇中的成员数过大。
注意:如果噪声聚类簇具有比其他聚类簇更大的组内离散度，则会增加其某些合法成员被视为异常值的可能性。
为了执行噪声聚类，vegclust（）函数使用范数标准化后物种数据（设定method=“NC”）。然后将噪声聚类的结果投射到主坐标排序图中，并用扇区表示模糊成员值。
library(vegclust)
?vegclust
k <- 4
#生成分4组的噪声聚类，执行随机种子30以保持最好的解决方案
spe.nc <- vegclust(
spe.norm, #先距离处理的数据
mobileCenters = k,
m = 1.5,
dnoise = 0.75,
method = "NC",
nstart = 30
)
spe.nc
# 物种的形心
(medoids <- spe.nc$mobileCenters)
# 模糊成员矩阵
spe.nc$memb
# 模糊聚类的基数（即每个对象属于每个聚类簇的数值）
spe.nc$size
# 获取硬成员值向量，N表示无法被分组的对象
spefuz.g <- defuzzify(spe.nc$memb)$cluster
clNum <- as.numeric(as.factor(spefuz.g))
# 模糊聚类排序 (PCoA)
dev.new(title = "Noise clustering of fish data - Ordination plot",
noRStudioGD = TRUE)
plot(dc.scores, #鱼类数据弦距离矩阵的主坐标分析
main = "模糊聚类排序 (PCoA)",
asp = 1,
type = "n")
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
for (i in 1:k)
{
gg <- dc.scores[clNum == i, ]
hpts <- chull(gg)
hpts <- c(hpts, hpts[1])
lines(gg[hpts, ], col = i + 1)
}
stars(
spe.nc$memb[, 1:4],
location = dc.scores,
key.loc = c(0.6, 0.4),
key.labels = 1:k,
draw.segments = TRUE,
add = TRUE,
# scale = FALSE,
len = 0.075,
col.segments = 2:(k + 1)
)
对被包含在clNum对象中的硬分类的组进行排序
# 模拟噪声样方聚类图
dev.new(title = "Defuzzified site plot",
noRStudioGD = TRUE)
plot(
dc.pcoa,
xlab = "MDS1",
ylab = "MDS2",
pch = clNum,
col = clNum
)
legend(
"topleft",
col = 1:(k + 1),
pch = 1:(k + 1),
legend = levels(as.factor(spefuz.g)),
bty = "n"
)
将模糊聚类结果在Doubs河流地图表示
# 将模糊聚类结果在Doubs河流地图表示
dev.new(title = "Four noise clusters on river",
width = 9,
noRStudioGD = TRUE)
plot(
spa,
asp = 1,
type = "n",
main = "沿着河流的噪声聚类簇",
xlab = "x 坐标(km)",
ylab = "y 坐标 (km)"
)
lines(spa, col = "light blue")
text(65, 20, "上游", cex = 1.2)
text(15, 32, "下游", cex = 1.2)
# 给模糊聚类簇加上扇形
for (i in 1:k) {
stars(
spe.nc$memb[, 1:4],
location = spa,
key.loc = c(150, 20),
key.labels = 1:k,
draw.segments = TRUE,
add = TRUE,
# scale = FALSE,
len = 5,
col.segments = 2:(k + 1)
)
}
虽然“硬”聚类产生的分类结果往往不太符合自然界真实的情况，但在需要明确区分对象时，硬聚类的结果很有用。相比之下，模糊聚类在描述对象归类时更加谨慎，也更贴近自然界的实际情况。此外，有专门描述数据连续结构的方法，即简单排序和约束排序，之后将会学习。
聚类方法有很多，关键在于如何选择合适的方法分析数据，并获得最佳的结果。
好了，第四章的内容就全部结束了，内容很多，方法很多，还是得好好消化吸收！
下一章我们将学习非约束分析，但是呢，这一章还是得好好理理，消化后再开始，下一章的学习。
如有不足或错误之处，请批评指正。
有什么不明白的也欢迎留言讨论。
欢迎关注同名公众号
往期内容：
《数量生态学：R语言的应用》第二版第三章-关联测度与矩阵------Q模式
《数量生态学——R语言的应用》第二版阅读笔记--绪论和第二章（一部分）
感谢你的阅读！！！你的点赞关注转发是对我最大的鼓励。
