本周开始我们的《数量生态学笔记》的第四章：聚类分析。聚类分析又称群分析，它是研究（样品或指标）分类问题的一种统计分析方法，同时也是数据挖掘的一个重要算法。在生态学研究当中，聚类的目的是识别环境中不连续对象的子集。实际上，聚类分析是所研究对象集合的分组。

需要注意大部分聚类方法都是基于关联矩阵进行计算，这也说明选择恰当的关联系数非常重要。

如图所示，我们需要识别不同类型的聚类方法及其应用的条件。

# Load required libraries 
library(ade4)
library (vegan)  # should be loaded after ade4 to avoid some conflicts!
library (gclus)
library(cluster)
library(RColorBrewer)
library(labdsv)
library(mvpart) #  在cran 没了，在想办法安装
library(MVPARTwrap) # Packages MVPARTwrap and rdaTest must be installed from # zip files  
#  在cran 没了，在想办法安装
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")

# Import the data from CSV files
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)
# Remove empty site 8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]

基于连接的层次聚类

单连接聚合聚类

单连接聚合聚类也称作最近临体聚类，该方法聚合对象的依据是最短的成对距离。每个对象或簇首次连接的列表成为主链接，也成为最小拓展树。

# 物种多度数据：先计算样方之间的弦距离矩阵，然后进行单连接聚合聚类
# **************************************************************
spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")
spe.ch.single <- hclust(spe.ch, method="single")
# 使用默认参数选项绘制聚类树
plot(spe.ch.single,main="聚类树",ylab="高度",xlab="单连接聚合聚类")
#基于第一次聚类的结果，如何描述这个数据集？是简单的单一梯度还是区
#分明显的样方组？能否辨认样方的连接链？样方1、5和9为什么最后连
#接？

完全连接聚合聚类

允许一个对象或簇与另一组聚合的依据是最远距离对。

# 计算完全连接聚合聚类
# ********************
spe.ch.complete <- hclust(spe.ch, method="complete")
plot(spe.ch.complete,main="聚类树",ylab="高度",xlab="完全连接聚合聚类")
#当前所给的样方是沿着河流分布（样方的编号按照流向编排），这个聚类分
#析结果是否将位置相近的样方排在同一个组呢？
#两种完全有效的聚类分析方法分析同一数据，为什么产生如此不同的聚类
#结果呢？

单一连接是一个对象很容易聚合到一个分组，因为单次连接足以导致融合。因此单连接聚类也被称为最亲密朋友法。虽然产生的分类组不清晰，但容易识别梯度。相反，完全连接聚类产生的分类之间差异比较明显。完全连接聚类更倾向与产生很多小的分离的组，更适合于寻找和识别数据的间断分布。

平均聚合聚类

平均聚合聚类是一类基于对象平均相异性或聚类簇中心的聚类方法。此类聚类有4种，不同的方法区别在于组的位置计算方式和计算融合时是否包含对象数量作为权重。

最著名的当属UPGMA方法，一个对象加入一个组的依据是这个对象与该组的每个成员之间的平均距离。

# 计算UPGMA聚合聚类
# ***********************
spe.ch.UPGMA <- hclust(spe.ch, method="average")
plot(spe.ch.UPGMA,main="聚类树",ylab="高度",xlab="UPGMA聚合聚类")
#这个UPGMA聚合聚类树看起来介于单连接聚类和完全连接聚类之间。这种
#情况经常发生。

需要注意的是 UPGMC和WPGMC有时会导致树的翻转，分类结果难以解读。

# 计算鱼类数据的形心聚类
# ***********************
spe.ch.centroid <- hclust(spe.ch, method="centroid")
plot(spe.ch.centroid,main="聚类树",ylab="高度",xlab="UPGMC聚合聚类")
#这种聚类树对生态学家来说简直是噩梦。Legendre 和Legendre（1998，
#第341页）解释了聚类树层级顺序倒置如何产生，并建议用多分法
#（polychotomies）代替二分法（dichotomies）解读这种图。

Ward 最小方差聚类

这是一种基于最小二乘法线性模型准则的聚类方法，分组依据是是组内平方和（即方差分析的方差）最小化。

# 计算Ward最小方差聚类
# ***********************
spe.ch.ward <- hclust(spe.ch, method="ward.D")
plot(spe.ch.ward,main="聚类树",ylab="高度",xlab="Ward聚类")
#使用距离平方造成此聚类树上半部分过于膨胀。为了使聚类树比例看起来
#更协调而不影响结构，可以使用当前融合水平的平方根重新绘图（图4.5）
spe.ch.ward$height <- sqrt(spe.ch.ward$height)
plot(spe.ch.ward,main="聚类树",ylab="高度",xlab="Ward聚类")

解读和比较聚类分析结果

需要牢记的是聚类分析是一种探索分析，而非统计检验。影响聚类结果的因素包括聚类方法本省和用于聚类分析的关联系数。

同表型相关

对已经完成层次聚类任意两个对象，在聚类树上从一个对象向上走，到达与另一个对象交回节点向下走，势必会到达第二个对象。交汇节点所在的层次水平即是两个对象的同表型距离 。

# 单连接聚类同表型相关
 spe.ch.single.coph <- cophenetic(spe.ch.single)
 cor(spe.ch, spe.ch.single.coph)
[1] 0.599193
 # 完全连接聚类同表型相关
 spe.ch.comp.coph <- cophenetic(spe.ch.complete)
 cor(spe.ch, spe.ch.comp.coph)
[1] 0.7655628
 # 平均聚类同表型相关
 spe.ch.UPGMA.coph <- cophenetic(spe.ch.UPGMA)
 cor(spe.ch, spe.ch.UPGMA.coph)
[1] 0.8608326
 # Ward聚类同表型相关
 spe.ch.ward.coph <- cophenetic(spe.ch.ward)
 cor(spe.ch, spe.ch.ward.coph)
[1] 0.7985079
 #哪个聚类树保持与原始的弦距离矩阵最接近的关系？
 #同表型相关也可以用spearman秩相关或Kendall秩相关表示
 cor(spe.ch, spe.ch.ward.coph, method="spearman")
[1] 0.7661171

为了描述一个距离矩阵与通过不同聚类方法得到的同表型矩阵之间的相关性，可以绘制原始距离对阵同表型距离的Shepard
图。

# Shepard图
# ***********
par(mfrow=c(2,2))
plot(spe.ch, spe.ch.single.coph, xlab="弦距离", 
  ylab="同表型距离", asp=1, xlim=c(0,sqrt(2)), ylim=c(0,sqrt(2)), 
  main=c("单连接",paste("同表型相关",
  round(cor(spe.ch, spe.ch.single.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.single.coph), col="red")
plot(spe.ch, spe.ch.comp.coph, xlab="弦距离", 
    ylab="同表型距离", asp=1, xlim=c(0,sqrt(2)), ylim=c(0,sqrt(2)),
  main=c("完全连接", paste("同表型相关",
  round(cor(spe.ch, spe.ch.comp.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.comp.coph), col="red")
plot(spe.ch, spe.ch.UPGMA.coph, xlab="弦距离", 
    ylab="同表型距离", asp=1, xlim=c(0,sqrt(2)), ylim=c(0,sqrt(2)), 
  main=c("UPGMA", paste("同表型相关",
  round(cor(spe.ch, spe.ch.UPGMA.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.UPGMA.coph), col="red")
plot(spe.ch, spe.ch.ward.coph, xlab="弦距离", 
    ylab="同表型距离", asp=1, xlim=c(0,sqrt(2)), 
  ylim=c(0,max(spe.ch.ward$height)),
  main=c("Ward聚类", paste("同表型相关",
  round(cor(spe.ch, spe.ch.ward.coph),3))))
abline(0,1)
lines(lowess(spe.ch, spe.ch.ward.coph), col="red")

Gower 距离

原始距离与同表型距离之间的差值平方和

 # Gower（1983）距离
gow.dist.single <- sum((spe.ch-spe.ch.single.coph)^2)
gow.dist.comp <- sum((spe.ch-spe.ch.comp.coph)^2)
gow.dist.UPGMA <- sum((spe.ch-spe.ch.UPGMA.coph)^2)
gow.dist.ward <- sum((spe.ch-spe.ch.ward.coph)^2)
gow.dist.single
[1] 95.41391
 gow.dist.comp
[1] 40.48897
 gow.dist.UPGMA
[1] 11.6746
 gow.dist.ward
[1] 532.0055

寻找可解读的聚类簇

为了解读和比较聚类的结果，通常需要寻找可解读的聚类簇，这意味着需要决定聚类树裁剪到哪个水平。

融合水平值图

聚类树的融合水平值是聚类树中两个分支融合处的相异性的数值。

# 融合水平值图
# *************
par(mfrow=c(2,2))
# 绘制单连接聚类融合水平值图
summary(spe.ch.single)  # 总结聚类分析的结果
plot(spe.ch.single$height, nrow(spe):2, type="S", main="融合水平值-弦距离-单连接", ylab="k （聚类簇数量）", xlab="h (节点高度）)", col="grey")
text(spe.ch.single$height, nrow(spe):2, nrow(spe):2, col="red", cex=0.8)
# 绘制完全连接聚类融合水平值图
plot(spe.ch.complete$height, nrow(spe):2, type="S", 
    main="融合水平值-弦距离-完全连接", 
    ylab="k （聚类簇数量）", xlab="h (节点高度）", col="grey")
text(spe.ch.complete$height, nrow(spe):2, nrow(spe):2, col="red", cex=0.8)
#在这里建议的组数可能不同。再次裁剪聚类树。这些解决方案是否有意义？
# 绘制UPGMA聚类融合水平值图
plot(spe.ch.UPGMA$height, nrow(spe):2, type="S", 
    main="融合水平值-弦距离-UPGMA", 
    ylab="k （聚类簇数量）", xlab="h (节点高度）", col="grey")
text(spe.ch.UPGMA$height, nrow(spe):2, nrow(spe):2, col="red", cex=0.8)
#绘制Ward聚类融合水平值图
plot(spe.ch.ward$height, nrow(spe):2, type="S", 
    main="融合水平值-弦距离-Ward", 
    ylab="k （聚类簇数量）", xlab="h (节点高度）", col="grey")
text(spe.ch.ward$height, nrow(spe):2, nrow(spe):2, col="red", cex=0.8)
#上面4个图看起来差异很大。记住，这些解决方案中任何一个都不是绝对
#正确，每个方案都可以为数据分组提供独特见解。

使用cutree()函数设定分类组的数量，并使用列联表比较分类差异。

 # 裁剪聚类树以获得k个分类组并使用列联表比较组之间的差异
 # *******************************************************
 # 设定聚类组的数量
 k <- 4  # 根据上面4个融合水平值图，可以观察到分4组水平在所有图里有小
 #的跳跃
 # 裁剪聚类树
 spebc.single.g <- cutree(spe.ch.single, k)
 spebc.complete.g <- cutree(spe.ch.complete, k)
 spebc.UPGMA.g <- cutree(spe.ch.UPGMA, k)
 spebc.ward.g <- cutree(spe.ch.ward, k)
 # 通过列联表比较分类结果
 # 单连接vs完全连接
 table(spebc.single.g, spebc.complete.g)
              spebc.complete.g
spebc.single.g 1 2 3 4
             1 1 0 0 0
             2 7 8 8 3
             3 0 1 0 0
             4 0 1 0 0
 # 单连接vs UPGMA
 table(spebc.single.g, spebc.UPGMA.g)
              spebc.UPGMA.g
spebc.single.g  1  2  3  4
             1  1  0  0  0
             2  0 15  0 11
             3  0  1  0  0
             4  0  0  1  0
 #单连接vs Ward
 table(spebc.single.g, spebc.ward.g)
              spebc.ward.g
spebc.single.g  1  2  3  4
             1  1  0  0  0
             2 10  5  8  3
             3  0  1  0  0
             4  0  1  0  0
 # 完全连接vs UPGMA
 table(spebc.complete.g, spebc.UPGMA.g)
                spebc.UPGMA.g
spebc.complete.g 1 2 3 4
               1 1 7 0 0
               2 0 9 1 0
               3 0 0 0 8
               4 0 0 0 3
 # 完全连接vs Ward
 table(spebc.complete.g, spebc.ward.g)
                spebc.ward.g
spebc.complete.g 1 2 3 4
               1 8 0 0 0
               2 3 7 0 0
               3 0 0 8 0
               4 0 0 0 3
 # UPGMA vs Ward
 table(spebc.UPGMA.g, spebc.ward.g)
             spebc.ward.g
spebc.UPGMA.g  1  2  3  4
            1  1  0  0  0
            2 10  6  0  0
            3  0  1  0  0
            4  0  0  8  3
 #如果两个聚类的结果完全一样，那么这个列联表每行和每列只有一个非零
 #数字，其他应该为0。此处并没有出现这种情况。如何解读这些列联表呢？
 #例如，单连接聚类第2组含有26个样方，这些样方在Ward聚类中被分散
 #到4个组里。

确定分组数

轮廓宽度图

轮廓宽度是描述一个对象与所属聚类归属程度的测度，是一个对象同组内其他对象的平均距离与该对象同最邻近聚类簇内所有对象的平均距离比较。

# 依据轮廓宽度图选择最优化的聚类簇数量（Rousseeuw质量指数）
 # ********************************************************
 # 绘制所有分类水平（除了k=1组的情况）轮廓宽度值（Ward 聚类）
 # 首先产生一个长度等于样方数量的空向量asw
 asw <- numeric(nrow(spe))
 # 其次循环获得轮廓宽度值并依次填入asw向量
 for (k in 2:(nrow(spe)-1)) {
        sil <- silhouette(cutree(spe.ch.ward, k=k), spe.ch)
        asw[k] <- summary(sil)$avg.width
        }
 # 选择最佳（最大）轮廓宽度值
 k.best <- which.max(asw)
 # 利用cluster程序包内函数plot.silhouette（）绘制轮廓宽度值k
 plot(1:nrow(spe), asw, type="h",
   main="轮廓宽度-最优聚类簇数（Ward聚类）",
        xlab="k (组数）", ylab="平均轮廓宽度")
 axis(1, k.best, paste("最优",k.best,sep="\n"), col="red", font=2,
   col.axis="red")
 points(k.best, max(asw), pch=16, col="red", cex=1.5)
 cat("", "Silhouette-optimal number of clusters k =", k.best, "\n",
        "with an average silhouette width of", max(asw), "\n")
 Silhouette-optimal number of clusters k = 2
 with an average silhouette width of 0.3658319
 # 屏幕上将输出如下内容：
 # 轮廓宽度最优的聚类簇数量k=2，此时平均轮廓宽度为0.365819
 #轮廓宽度法经常选择2组作为最优的分类数量。然而，从生态学角度分析
 #分4组似乎更合理。

距离矩阵与代表分组的二元矩阵的比较

# 依据Mantel统计（Pearson）相关选择最优分组数量
# **********************************************
# 编写计算代表分类水平的二元距离矩阵的函数
grpdist <- function(X)
  {
  require(cluster)
  gr <- as.data.frame(as.factor(X))
  distgr <- daisy(gr, "gower")
  distgr
  }
# 基于Ward 聚类结果运行该函数
kt <- data.frame(k=1:nrow(spe), r=0)
for (i in 2:(nrow(spe)-1)) {
    gr <- cutree(spe.ch.ward, i)
    distgr <- grpdist(gr)
    mt <- cor(spe.ch, distgr, method="pearson")
    kt[i,2] <- mt
}
head(kt)
   k         r
1   1 0.0000000
2   2 0.6554297
3   3 0.7080851
4   4 0.7154912
5   5 0.6429201
6   6 0.6741193
7   7 0.6865197
8   8 0.6879179
9   9 0.6925904
10 10 0.6498757
k.best <- which.max(kt$r)
# 通过cluster程序包内plot.silhouette函数绘制分析图
plot(kt$k, kt$r, type="h", main="Mantel-最优聚类簇数（Ward聚类）", 
    xlab="k (组数）", ylab="Pearson 相关")
axis(1, k.best, paste("最优", k.best, sep="\n"), col="red", font=2,
    col.axis="red")
points(k.best, max(kt$r), pch=16, col="red", cex=1.5)

最终分组的轮廓图

# 最终分组的轮廓图
# ****************
# 选择聚类簇的数量
k <- 4
# 轮廓图
cutg <- cutree(spe.ch.ward, k=k)
sil <- silhouette(cutg, spe.ch)
silo <- sortSilhouette(sil)
rownames(silo) <- row.names(spe)[attr(silo,"iOrd")]
plot(silo, main="轮廓宽度图-弦距离（Ward聚类）", 
    cex.names=0.8, col=cutg+1, nmax.lab=100,border="white"
)
#组1和组3最连贯，同时组2可能含有被错分的对象。

利用绘图工具修饰的最终聚类树

# 设定组数的最终聚类树
# *********************
# 函数reorder.hclust（）的作用是重新排列从函数hclust（）获得的聚类树，使
# 聚类树内对象的排列顺序与原始相异矩阵内对象的排列顺序尽可能一致。重排# 不影响聚类树的结构。
spe.chwo <- reorder.hclust(spe.ch.ward, spe.ch)
# 绘制重排后带组标识的聚类树
plot(spe.chwo, hang=-1, xlab="4 groups", sub="", ylab="Height", 
    main="Chord - Ward (reordered)", labels=cutree(spe.chwo, k=k))
rect.hclust(spe.chwo, k=k)
# 绘制带不同颜色的最终聚类树 
# 使用我们自编函数hcoplot（）可以快速获得最终聚类树:
source("hcoplot.R")           # hcoplot.R脚本必须在当前工作目前下
hcoplot(spe.ch.ward, spe.ch, k=4)

聚类结果空间分布

# 4个Ward聚类簇在Doub河的分布情况
# ********************************
# 绘制Doubs河流地图（也见第2章）
plot(spa, asp=1, type="n", main="4个Ward聚类组", 
    xlab="x坐标 (km)", ylab="y坐标 (km)")
lines(spa, col="light blue")
text(70, 10, "上游", cex=1.2, col="red")
text(15, 115, "下游", cex=1.2, col="red")
# 添加4组分类信息
grw <- spebc.ward.g
k <- length(levels(factor(grw)))
for (i in 1:k) {
   points(spa[grw==i,1], spa[grw==i,2], pch=i+20, cex=3, col=i+1, bg=i+1)
   }
text(spa, row.names(spa), cex=0.8, col="white", font=2)
legend("bottomright", paste("组",1:k), pch=(1:k)+20, col=2:(k+1), 
  pt.bg=2:(k+1), pt.cex=1.5, bty="n")
#请比较当前生成的地图与第2章生成的4种鱼类物种分布地图。

热图和群落图的群落表

# 热图（heat map）
# ***************
# 用聚类结果重排距离矩阵的热图
dend <- as.dendrogram(spe.chwo)
heatmap(as.matrix(spe.ch), Rowv=dend, symm=TRUE, margin=c(3,3))
#观察如何设定使最热的色彩（计算机输出的是黑色或红色）代表最近的相似
#性，例如对角线代表对象自身的相似性，所以颜色最深。

# 重排群落表格
# 物种按照在样方得分加权平均进行排列
or <- vegemite(spe, spe.chwo)
# 基于聚类树的双排列群落表格的热图
heatmap(t(spe[rev(or$species)]), Rowv=NA, Colv=dend,
    col=c("white", brewer.pal(5,"Greens")), scale="none", margin=c(4,4), 
    ylab="物种（样方的加权平均）", xlab="样方")

参考：
常用聚类算法有哪些？六大类聚类算法详细介绍
无监督学习— 聚类（Clustering）
百科||聚类算法
聚类分析