物种空间聚类
1 聚类过程
1) 数据准备:包括特征标准化和降维。
2) 特征选择:从最初的特征中选择最有效的特征,并将其存储于向量中。
3) 特征提取:通过对所选择的特征进行转换形成新的突出特征。
4) 聚类(或分组):首先选择合适特征类型的某种距离函数(或构造新的距离函数)进行接近程度的度量; 而后执行聚类或分组。
5) 聚类结果评估:是指对聚类结果进行评估.评估主要3 种，外部有效性评估、内部有效性评估和相关性测试评估.。
2 聚类
2.1 加载函数包
library(ade4)
library(vegan) #应该先加载ade4后加载vegan以避免冲突
library(gclus)
library(cluster)
library(RColorBrewer)
library(labdsv)
library(mvpart)
library(MVPARTwrap) # MVPARTwrap这个程序包必须从本地zip文件安装setwd('C:\\Users\\dell\\Desktop\\谈德林_物种空间聚类_20190624\\原始资料及数据')
# 导入CSV格式的数据
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)以上数据读入后，分别用str，查看数据后，需要对异常情况的数据进行处理。
# 剔除无物种数据的样方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]2.2 单连接聚合聚类
物种多度数据：先计算样方之间的弦距离矩阵，然后进行单连连接聚合聚类。
spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")
spe.ch.single <- hclust(spe.ch, method="single")对单接聚合聚类的结果进行可视化，使用默认参数选项绘制聚类树，如下所示：
plot(spe.ch.single)
通过上述图形可以看到，左侧的系数值显示，聚成2类的效果为0.8左右.说明9号物种和其它物种的差异较为明显，其它物种可能临近亲疏效果更好。
2.3 完全连接聚合聚类
spe.ch.complete <- hclust(spe.ch, method="complete")
plot(spe.ch.complete)
完全连接对噪声和离群点不太敏感，但是它可能使大的簇破裂，并且偏好于球形。实现了组平均。组平均：定义簇邻近度为取自不同簇的所有点的平均逐对邻近度。
2.4 计算UPGMA聚合聚类
这个UPGMA聚合聚类树是介于单连接聚类和完全连接聚类之间的一种聚类。
spe.ch.UPGMA <- hclust(spe.ch, method="average")
plot(spe.ch.UPGMA)
用来重建系统发生树时 ,其假定的前提条件是:在进化过程中 ,每一世系发生趋异的次数相同 ,即核苷酸或氨基酸的替换速率是均等且恒定的。通过 UPGMA 法所产生的系统发生树可以说是物种树的简单体现 ,在每一次趋异发生后 ,从共祖节点到 2 个 OTU 间的支的长度一样。因此 ,这种方法较多地用于物种树的重建。 UPGMA 法在算法上较简单。聚类时 ,首先将距离最小的 2 个 OTU 聚在一起 ,并形成一个新的OTU ,其分支点位于 2 个 OTU 间距离的 1/ 2 处;然后计算新的 OTU 与其它 OTU 间的平均距离 ,再找出其中的最小 2 个 OTU 进行聚类;如此反复 ,直到所有的 OTU 都聚到一起 ,最终得到一个完整的系统发生树。
2.5 形心聚类
spe.ch.centroid <- hclust(spe.ch, method="centroid")
plot(spe.ch.centroid)
该聚类通过一种翻转产生，Legendre 和Legendre解释了聚类树翻转如何产生，并建议用多分法（polychotomies），代替二分法（dichotomies）解读这种图。
2.6 Ward最小方差聚类
最小化所有聚类内的平方差总和，使得下述物种间方差最小化的优化方向。
spe.ch.ward <- hclust(spe.ch, method="ward")## The "ward" method has been renamed to "ward.D"; note new "ward.D2"plot(spe.ch.ward)
使用距离平方造成此聚类树上半部分过于膨胀。为了使聚类树比例看起来更协调而不影响结构，可以使用当前融合水平的平方根重新绘图。如下所示：
spe.ch.ward$height <- sqrt(spe.ch.ward$height)
plot(spe.ch.ward)
3 聚类效果比较
3.1 同表型相关
采用同表型相关，可以比较聚类结果的优劣，选择合适的聚类结果。
spe.ch.single.coph <- cophenetic(spe.ch.single)
cor(spe.ch, spe.ch.single.coph)## [1] 0.599193单连接聚类同表型相关，得到系数值为0.599193。
3.2 完全连接聚类同表型相关
spe.ch.comp.coph <- cophenetic(spe.ch.complete)
cor(spe.ch, spe.ch.comp.coph)## [1] 0.7655628完全连接聚类同表型相关，得到系数值为0.7655628。
3.3 平均聚类同表型相关
spe.ch.UPGMA.coph <- cophenetic(spe.ch.UPGMA)
cor(spe.ch, spe.ch.UPGMA.coph)## [1] 0.8608326平均聚类同表型相关，得到系数值为0.8608326。
3.4 Ward聚类同表型相关
spe.ch.ward.coph <- cophenetic(spe.ch.ward)
cor(spe.ch, spe.ch.ward.coph)## [1] 0.7985079Ward聚类同表型相关，得到系数值为0.7985079。
综上，平均聚类同表型相关结果较好。
3.5 同表型相关spearman秩相关
哪个聚类树保持与原始的弦距离矩阵最接近，则说明聚类效果较好。
cor(spe.ch, spe.ch.ward.coph, method="spearman")## [1] 0.76611714 解读
4.1 Gower距离
gow.dist.single <- sum((spe.ch-spe.ch.single.coph)^2)
gow.dist.single## [1] 95.41391单一聚类和原始的弦距离矩阵的距离为95.41391
gow.dist.comp <- sum((spe.ch-spe.ch.comp.coph)^2)
gow.dist.comp## [1] 40.48897完全聚类和原始的弦距离矩阵的距离为40.48897
gow.dist.UPGMA <- sum((spe.ch-spe.ch.UPGMA.coph)^2)
gow.dist.UPGMA## [1] 11.6746UPGMA聚类和原始的弦距离矩阵的距离为11.6746
gow.dist.ward <- sum((spe.ch-spe.ch.ward.coph)^2)
gow.dist.ward## [1] 532.0055聚类和原始的弦距离矩阵的距离为532.0055。
4.2 融合水平值图
绘制单连接聚类融合水平值图，随着节点高度增加，聚类数量减少。
4.3 单一连接水平值图
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)
4.4 完全连接水平值图
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)
绘制完全连接聚类融合水平值图，随着节点高度增加，聚类数量减少。当节点高度为1.4时，聚类数量为2.
4.5 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)
绘制 UPGMA聚类融合水平值图，随着节点高度增加，聚类数量减少。当节点高度为1.2时，聚类数量为2.
4.6 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)
绘制 Ward聚类融合水平值图，随着节点高度增加，聚类数量减少。当节点高度为2.0和2.5后，聚类数量为2.
4.7 选择聚类簇的数量
k <- 4
cutg <- cutree(spe.ch.ward, k=k)
#设定组数的最终聚类树
# 函数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)
通过上述图形，得知最佳聚类数量为4.
4.8 聚类树可视化
使用自编函数hcoplot（）可以快速获得最终聚类树。能绘制带不同颜色的最终聚类树。
4.9 自编函数绘制聚类树
source("C:\\Users\\dell\\Desktop\\谈德林_物种空间聚类_20190624\\原始资料及数据\\hcoplot.R") # hcoplot.R脚本必须在当前工作目前下
hcoplot(spe.ch.ward, spe.ch, k=4)
上述的参数hang=-1是设定聚类树的分支终端从0的位置开始，并且聚类簇的标识将在这个设定值之下.得到的聚类树。
4.10 转化数据k-均值划分
注意：即使给定的nstart相同，每次运行上述命令，所产生的结果也不一定 完全相同，因为每次运算设定的初始结构是随机的。
以PAM分4组情况（partitioning around medoids，PAM）进行划分。
spe.kmeans <- kmeans(spe.norm, centers=4, nstart=100)
spe.kmeans ## K-means clustering with 4 clusters of sizes 8, 3, 6, 12
##
## Cluster means:
## CHA TRU VAI LOC OMB BLA
## 1 0.00000000 0.006691097 0.02506109 0.06987391 0.006691097 0.006691097
## 2 0.00000000 0.000000000 0.00000000 0.00000000 0.000000000 0.000000000
## 3 0.06167791 0.122088022 0.26993915 0.35942538 0.032664966 0.135403325
## 4 0.10380209 0.542300691 0.50086515 0.43325916 0.114024105 0.075651573
## HOT TOX VAN CHE BAR SPI
## 1 0.10687104 0.09377516 0.14194394 0.2011411 0.24327992 0.1326062
## 2 0.05205792 0.00000000 0.07647191 0.3166705 0.00000000 0.0000000
## 3 0.06212775 0.21568957 0.25887226 0.2722562 0.15647062 0.1574388
## 4 0.00000000 0.00000000 0.06983991 0.1237394 0.02385019 0.0000000
## GOU BRO PER BOU PSO ROT
## 1 0.28386032 0.20630360 0.16920496 0.2214275 0.19066542 0.13171275
## 2 0.20500174 0.07647191 0.00000000 0.0000000 0.05205792 0.07647191
## 3 0.16822286 0.12276089 0.17261621 0.0793181 0.06190283 0.04516042
## 4 0.05670453 0.04722294 0.02949244 0.0000000 0.00000000 0.00000000
## CAR TAN BCO PCH GRE GAR
## 1 0.16019126 0.26230024 0.19561641 0.1331835 0.26713081 0.32103755
## 2 0.00000000 0.00000000 0.00000000 0.0000000 0.18058775 0.31667052
## 3 0.06190283 0.14539027 0.01473139 0.0000000 0.03192175 0.32201597
## 4 0.00000000 0.03833408 0.00000000 0.0000000 0.00000000 0.01049901
## BBO ABL ANG
## 1 0.22883055 0.3326939 0.18873077
## 2 0.05205792 0.7618709 0.00000000
## 3 0.01473139 0.1095241 0.04739636
## 4 0.00000000 0.0000000 0.00000000
##
## Clustering vector:
## 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 4 4 4 4 3 4 4 3 4 4 4 4 4 4 3 3 3 3 1 1 1 2 2 2 1
## 27 28 29 30
## 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 0.4696535 0.3560423 1.7361453 2.5101386
## (between_SS / total_SS = 66.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"spe.ch.pam <- pam(spe.ch, k=4, diss=TRUE)#default for dist or dissimilarity
summary(spe.ch.pam)## Medoids:
## ID
## [1,] "11" "12"
## [2,] "6" "6"
## [3,] "16" "17"
## [4,] "25" "26"
## Clustering vector:
## 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
## 1 1 1 2 2 2 1 2 2 1 1 1 1 1 3 3 3 3 4 4 4 4 4 4 4
## 27 28 29 30
## 4 4 4 4
## Objective function:
## build swap
## 0.4426460 0.4390712
##
## Numerical information per cluster:
## size max_diss av_diss diameter separation
## [1,] 9 0.8769750 0.4350972 1.1168811 0.3429386
## [2,] 5 0.8814766 0.4962134 1.0854578 0.3429386
## [3,] 4 0.5988903 0.3678822 0.8037373 0.4977122
## [4,] 11 0.8841837 0.4422357 1.1276080 0.4977122
##
## Isolated clusters:
## L-clusters: character(0)
## L*-clusters: character(0)
##
## Silhouette plot information:
## cluster neighbor sil_width
## 13 1 2 0.42988827
## 12 1 2 0.40576079
## 2 1 2 0.38341138
## 14 1 2 0.35676210
## 11 1 2 0.34188972
## 3 1 2 0.27044354
## 1 1 2 0.24590101
## 7 1 2 0.24135661
## 15 1 3 0.05061851
## 6 2 1 0.09343621
## 10 2 1 0.06179173
## 9 2 3 0.04916887
## 5 2 3 -0.03817150
## 4 2 1 -0.14708925
## 17 3 2 0.45357978
## 18 3 4 0.38401410
## 16 3 2 0.23661137
## 19 3 4 0.15364654
## 26 4 3 0.48031045
## 27 4 3 0.46380689
## 28 4 3 0.45849776
## 22 4 3 0.41331416
## 21 4 3 0.39375977
## 29 4 3 0.36426967
## 30 4 3 0.36420744
## 24 4 3 0.31070288
## 25 4 3 0.26599024
## 20 4 3 0.23649622
## 23 4 3 0.21029600
## Average silhouette width per cluster:
## [1] 0.302892439 0.003827211 0.306962948 0.360150134
## Average silhouette width of total data set:
## [1] 0.2736094
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call"4.11 k-均值划分与4个环境变量分析
attach(env)
# 定量环境变量箱线图
# 海拔、坡度、氧含量和铵浓度
par(mfrow=c(2,2))
boxplot(sqrt(alt) ~ spe.kmeans$cluster, main="海拔", las=1, ylab="sqrt(alt)", col=2:5, varwidth=TRUE)
boxplot(log(pen) ~ spe.kmeans$cluster, main="坡度", las=1, ylab="log(pen)", col=2:5, varwidth=TRUE)
boxplot(oxy ~ spe.kmeans$cluster, main="氧含量", las=1, ylab="oxy", col=2:5, varwidth=TRUE)
boxplot(sqrt(amm) ~ spe.kmeans$cluster, main="铵浓度", las=1, ylab="sqrt(amm)", col=2:5, varwidth=TRUE)
通过不同的k均值划分后，得到海拔、坡度、氧含量和铵浓度和环境变量之间的关系，如上所示。聚类类型为3的，海拔、坡度。含氧量都较高与其他聚类的类型。
5 检验
5.1 方差分析检验假设
假设检验（Hypothesis Testing）是数理统计学中根据一定假设条件由样本推断总体的一种方法。
5.2 检验残差的正态性
shapiro.test(resid(lm(sqrt(alt) ~ as.factor(spe.kmeans$cluster))))##
## Shapiro-Wilk normality test
##
## data: resid(lm(sqrt(alt) ~ as.factor(spe.kmeans$cluster)))
## W = 0.96238, p-value = 0.3756shapiro.test(resid(lm(log(pen) ~ as.factor(spe.kmeans$cluster))))##
## Shapiro-Wilk normality test
##
## data: resid(lm(log(pen) ~ as.factor(spe.kmeans$cluster)))
## W = 0.95182, p-value = 0.204shapiro.test(resid(lm(oxy ~ as.factor(spe.kmeans$cluster))))##
## Shapiro-Wilk normality test
##
## data: resid(lm(oxy ~ as.factor(spe.kmeans$cluster)))
## W = 0.93941, p-value = 0.09674shapiro.test(resid(lm(sqrt(amm) ~ as.factor(spe.kmeans$cluster))))##
## Shapiro-Wilk normality test
##
## data: resid(lm(sqrt(amm) ~ as.factor(spe.kmeans$cluster)))
## W = 0.95398, p-value = 0.2319检验结果表明sqrt(alt)、log(pen)、oxy和sqrt(amm)的残差是正态分布。也尝试为其他的环境变量寻找好的标准化转化。
5.3 检验方差齐性
bartlett.test(sqrt(alt), as.factor(spe.kmeans$cluster))##
## Bartlett test of homogeneity of variances
##
## data: sqrt(alt) and as.factor(spe.kmeans$cluster)
## Bartlett's K-squared = 16.953, df = 3, p-value = 0.0007227bartlett.test(log(pen), as.factor(spe.kmeans$cluster))##
## Bartlett test of homogeneity of variances
##
## data: log(pen) and as.factor(spe.kmeans$cluster)
## Bartlett's K-squared = 2.9832, df = 3, p-value = 0.3942bartlett.test(oxy, as.factor(spe.kmeans$cluster))##
## Bartlett test of homogeneity of variances
##
## data: oxy and as.factor(spe.kmeans$cluster)
## Bartlett's K-squared = 1.8258, df = 3, p-value = 0.6093bartlett.test(sqrt(amm), as.factor(spe.kmeans$cluster))##
## Bartlett test of homogeneity of variances
##
## data: sqrt(amm) and as.factor(spe.kmeans$cluster)
## Bartlett's K-squared = 5.7203, df = 3, p-value = 0.126变量sqrt(alt)的方差不齐，所以参数检验的方差分析不适用。
5.4 变量的方差分析
summary(aov(log(pen) ~ as.factor(spe.kmeans$cluster)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(spe.kmeans$cluster) 3 18.05 6.018 7.283 0.00114 **
## Residuals 25 20.66 0.826
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(oxy ~ as.factor(spe.kmeans$cluster)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(spe.kmeans$cluster) 3 103.54 34.51 26.27 6.75e-08 ***
## Residuals 25 32.84 1.31
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aov(sqrt(amm) ~ as.factor(spe.kmeans$cluster)))## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(spe.kmeans$cluster) 3 2.6126 0.8709 49.55 1.15e-10 ***
## Residuals 25 0.4394 0.0176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1从上述结果看出，坡度、含氧量和铵浓度在不同聚类簇之间是无显著性差异较为明显。
5.5 变量alt的Kruskal-Wallis检验
Kruskal-Wallis检验用来检测总体函数分布的一致性原假设和其替代假设，关于至少两个样本之间存在差异的假设。
kruskal.test(alt ~ as.factor(spe.kmeans$cluster))##
## Kruskal-Wallis rank sum test
##
## data: alt by as.factor(spe.kmeans$cluster)
## Kruskal-Wallis chi-squared = 21.526, df = 3, p-value = 8.186e-05从上述看出，p值远小于0.05，说明海拔在不同聚类簇之间是无显著性差异较为明显。
detach(env)6 本章汇总
| 参数 | 类别 | 功能 |
|---|---|---|
| ade4 | 包 | 可绘制最小生成树进行对应分析 |
| vegan | 包 | 提供非度量多维尺度分析的包装 |
| gclus | 包 | 针对聚类的散点图和平行坐标图 |
| cluster | 包 | 实现聚类 |
| FD | 包 | 计算功能多样性 |
| RColorBrewer | 包 | 颜色调节 |