博文
数量生态学笔记||聚类物种关联
|
组内数据简单统计
对聚类结果获得的样方组进行简单统计（如平均多度），寻找每个样方组数量最多，频度最高或最具代表性的物种。
# 加载所需的程序包
library(ade4)
library(vegan) # 应该先加载ade4后加载vegan以避免冲突
library(gclus)
library(cluster)
library(RColorBrewer)
library(labdsv)
library(mvpart)
library(MVPARTwrap) # MVPARTwrap这个程序包必须从本地zip文件安装
# 导入CSV格式的数据
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)
# 删除无物种数据的样方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")
# 预转化后物种数据k-均值划分
# ****************************
spe.kmeans <- kmeans(spe.norm, centers=4, nstart=100)# k-均值样方聚类簇平均多度
# ***********************
groups <- as.factor(spe.kmeans$cluster)
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[,i], spe.kmeans$cluster, mean)
}
# 每组的物种平均多度
group1 <- round(sort(spe.means[,1], decreasing=TRUE), 2)
group2 <- round(sort(spe.means[,2], decreasing=TRUE), 2)
group3 <- round(sort(spe.means[,3], decreasing=TRUE), 2)
group4 <- round(sort(spe.means[,4], decreasing=TRUE), 2)
> # 显示多度大于平均值的物种
> group1.domin <- which(group1>mean(group1))
> group1
TRU VAI LOC OMB CHE CHA BLA VAN GOU BRO TAN BAR PER GAR HOT TOX SPI BOU PSO ROT CAR BCO PCH GRE BBO
4.00 4.00 3.58 1.00 1.00 0.92 0.75 0.58 0.50 0.42 0.33 0.25 0.25 0.08 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
ABL ANG
0.00 0.00
> group1.domin
TRU VAI LOC OMB CHE CHA BLA
1 2 3 4 5 6 7Kendall 共性系数（W）
这个方法的目标是找到最全面的物种集合，即利用最少组数并尽可能让显著性正相关的物种分在同一组。此处给一个简单的案例：首先从鱼类数据中提取出多度最大的一部分物种，使用K-均值法将这些物种分为几组，然后运行全局检验（Kendall,glpba()）识别是否所有物种组显著关联。如果显著关联，将对每个组物种进行后验概率检验（Kendall，post()）验证同一组内部时候具有共性。
# Kendall一致性系数（W）
# ********************
# 提取多度最大的物种
sp.sum <- apply(spe, 2, sum)
spe.sorted <- spe[,order(sp.sum, decreasing=TRUE)]
spe.small <- spe.sorted[,1:20]
# 物种数据转化和矩阵转置
spe.small.hel <- decostand(spe.small, "hellinger")
spe.small.std <- decostand(spe.small.hel, "standardize")
spe.small.t <- t(spe.small.std)
# 物种k-均值划分
spe.t.kmeans.casc <- cascadeKM(spe.small.t, inf.gr=2, sup.gr=8,
iter=100, criterion="calinski")
plot(spe.t.kmeans.casc, sortg=TRUE)
#结果显示分2组是不错的选择，如果分更多的组，可能有些组只有一个物种。
> # 在object $partition中显示分为2组
> clusters <- spe.t.kmeans.casc$partition[,1]
> clusters
LOC VAI GAR TRU ABL CHE GOU TAN VAN BAR BRO GRE PER BOU BBO PSO SPI ANG TOX BCO
1 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
> #现在有两个物种组，对这两组运行全局的Kendall W检验
> spe.kendall.global <- kendall.global(spe.small.hel, clusters)
> spe.kendall.global
$`Concordance_analysis`
Group.1 Group.2
W 8.363838e-01 4.545404e-01
F 1.022373e+01 1.333306e+01
Prob.F 3.486929e-13 5.222120e-43
Corrected prob.F 3.486929e-13 1.044424e-42
Chi2 7.025624e+01 2.163612e+02
Prob.perm 1.000000e-03 1.000000e-03
Corrected prob.perm 2.000000e-03 2.000000e-03
$Correction.type
[1] "holm"
attr(,"class")
[1] "kendall.global"
> spe.kendall.post <- kendall.post(spe.small.hel, clusters, nperm=9999)
> spe.kendall.post
$`A_posteriori_tests_Group`
$`A_posteriori_tests_Group`[[1]]
LOC VAI TRU
Spearman.mean 0.7414122 0.8232272 0.6976134
W.per.species 0.8276082 0.8821515 0.7984089
Prob 0.0001000 0.0001000 0.0001000
Corrected prob 0.0020000 0.0020000 0.0020000
$`A_posteriori_tests_Group`[[2]]
GAR ABL CHE GOU TAN VAN BAR BRO GRE PER
Spearman.mean 0.3776566 0.4562921 -0.09825354 0.4707461 0.4287597 0.2341254 0.4527895 0.3087706 0.5156226 0.4342510
W.per.species 0.4142650 0.4882749 -0.03365039 0.5018787 0.4623621 0.2791768 0.4849784 0.3494312 0.5441154 0.4675303
Prob 0.0012000 0.0001000 0.75610000 0.0001000 0.0005000 0.0389000 0.0001000 0.0091000 0.0001000 0.0004000
Corrected prob 0.0048000 0.0020000 0.75610000 0.0020000 0.0025000 0.0778000 0.0020000 0.0273000 0.0020000 0.0024000
BOU BBO PSO SPI ANG TOX BCO
Spearman.mean 0.5823022 0.5171591 0.5636835 0.5353123 0.5848444 0.4729478 0.5419302
W.per.species 0.6068727 0.5455615 0.5893492 0.5626469 0.6092653 0.5039508 0.5688755
Prob 0.0001000 0.0001000 0.0001000 0.0001000 0.0001000 0.0001000 0.0001000
Corrected prob 0.0020000 0.0020000 0.0020000 0.0020000 0.0020000 0.0020000 0.0020000
$Correction.type
[1] "holm"
attr(,"class")
[1] "kendall.post"基于有无数据的物种集合
p值可以当做距离：共发生度越高的两个物种，其p值越小。
# 基于有-无数据的物种集合
# ***********************
source("test.a.R") # 函数必须在当前工作目录下
spe.pa <- decostand(spe,"pa") # 将数据转化为有-无数据
res <- test.a(spe.pa, nperm=99999)
Computation time = 246.690000 sec
summary(res)
Length Class Mode
a 729 -none- numeric
p.a 729 -none- numeric
p.a.dist 351 dist numeric
#输出的结果是res$p.a.dsit包含一个p值的矩阵。下一步计算向量化的p值矩阵
#的Holm校正数（见7.2.6节）。
res.p.vec <- as.vector(res$p.a.dist)
adjust.res <- p.adjust(res.p.vec,method="holm")
#在Holm校正后的p值中，可以发现0.05或比0.05小但很接近的值。例如，数
#值0.04878的未校正p值应该是0.00018。Holm校正后的p值矩阵内大约有83
#个值小于0.05（可能每次运行这个数值有轻微变化）。
#因此0.00018或更小概率才能称为显著。接下来将相异矩阵内将所有大于#0.00018的值都换成1。
res.pa.dist <- res$p.a.dist
res.pa.dist[res.pa.dist>0.00018] <- 1
#由p值构成的相异矩阵可以用热图表示（见第3章）
source("coldiss.R") # 函数必须在当前工作目录下
coldiss(res.pa.dist, nc=10, byrank=TRUE, diag=TRUE)
> # 运行模糊非层次聚类（见4.12节）。尝试使用多k值获得最佳的分组。
> res.pa.fuz <- fanny(res.pa.dist, k=5, memb.exp=1.5)
> summary(res.pa.fuz)
Fuzzy Clustering object of class 'fanny' :
m.ship.expon. 1.5
objective 3.084468
tolerance 1e-15
iterations 73
converged 1
maxit 500
n 27
Membership coefficients (in %, rounded):
[,1] [,2] [,3] [,4] [,5]
CHA 99 0 0 0 0
TRU 22 26 13 15 23
VAI 22 26 13 15 23
LOC 22 26 13 15 23
OMB 99 0 0 0 0
BLA 22 26 13 15 23
HOT 0 0 100 0 0
TOX 0 0 100 0 0
VAN 22 26 13 15 23
CHE 22 26 13 15 23
BAR 0 0 100 0 0
SPI 0 0 100 0 0
GOU 9 70 6 6 9
BRO 9 70 6 6 9
PER 0 95 4 0 0
BOU 0 0 100 0 0
PSO 0 0 100 0 0
ROT 0 0 0 100 0
CAR 0 0 100 0 0
TAN 9 63 11 7 10
BCO 0 0 2 98 0
PCH 1 2 1 95 1
GRE 3 4 16 64 13
GAR 0 1 0 0 98
BBO 0 0 0 100 0
ABL 0 0 2 0 97
ANG 0 0 100 0 0
Fuzzyness coefficients:
dunn_coeff normalized
0.7335123 0.6668904
Closest hard clustering:
CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE BAR SPI GOU BRO PER BOU PSO ROT CAR TAN BCO PCH GRE GAR BBO ABL ANG
1 2 2 2 1 2 3 3 2 2 3 3 2 2 2 3 3 4 3 2 4 4 4 5 4 5 3
Silhouette plot information:
cluster neighbor sil_width
CHA 1 2 0.9999500
OMB 1 2 0.9999500
BRO 2 1 0.1111100
GOU 2 1 0.1110933
TRU 2 1 0.0000000
VAI 2 1 0.0000000
LOC 2 1 0.0000000
BLA 2 1 0.0000000
VAN 2 1 0.0000000
CHE 2 1 0.0000000
TAN 2 3 -0.1562075
PER 2 3 -0.6249531
SPI 3 5 0.9999800
TOX 3 2 0.9999794
CAR 3 5 0.9999771
BOU 3 4 0.9999750
ANG 3 4 0.9999750
BAR 3 4 0.9999732
HOT 3 4 0.9999250
PSO 3 5 0.8241758
ROT 4 1 0.9999375
BBO 4 5 0.9999050
BCO 4 3 0.9998201
PCH 4 1 0.7499250
GRE 4 3 0.3333267
GAR 5 3 0.9999543
ABL 5 3 0.9998400
Average silhouette width per cluster:
[1] 0.99995000 -0.05589573 0.97799507 0.81658285 0.99989716
Average silhouette width of total data set:
[1] 0.5684301
Available components:
[1] "membership" "coeff" "memb.exp" "clustering" "k.crisp" "objective" "convergence" "diss"
[9] "call" "silinfo"> plot(silhouette(res.pa.fuz), main="Silhouette plot - Fuzzy clustering", cex.names=0.8, col=res.pa.fuz$silinfo$widths+1)
> res.pa.fuz$silinfo # 轮廓信息
$`widths`
cluster neighbor sil_width
CHA 1 2 0.9999500
OMB 1 2 0.9999500
TRU 2 1 0.0000000
VAI 2 1 0.0000000
LOC 2 1 0.0000000
BLA 2 1 0.0000000
VAN 2 1 0.0000000
CHE 2 1 0.0000000
GOU 2 1 0.0000000
GAR 2 3 -0.2222000
PSO 3 5 0.9999656
ANG 3 5 0.9999625
BOU 3 5 0.9999594
HOT 3 5 0.9999219
CAR 3 4 0.9999138
SPI 3 4 0.9999138
BAR 3 1 0.8749662
TOX 3 4 0.8124859
ABL 3 5 0.3749281
PER 4 3 0.9999400
BRO 4 1 0.4999900
TAN 4 3 0.3571576
ROT 5 1 0.9999300
BCO 5 3 0.9998875
BBO 5 3 0.9998575
PCH 5 1 0.7499225
GRE 5 3 0.2500275
$clus.avg.widths
[1] 0.9999500 -0.0277750 0.8957797 0.6190292 0.7999250
$avg.width
[1] 0.5813493
#此轮廓信息包含一个聚类簇中最可能的成员（"cluster"）、最近邻体聚类
#簇（"neighbor"）和轮廓宽度值("sil_width")。
InVal 物种指示值
使用InVal方法寻找的指示种与分组类型有充分必要的关系，即如果找到该物种，就可找到该类型；如果找到该类型，就能找到该物种。物种指示值高表示该物种在该样方内平均度大于其他样方（特异性），并且该物种在该组内绝大部分样方都存在。
IndVal 值根据以下公式进行计算( Dufrêne and Legendre，1997) :
Aij = Nindividualsij /Nindividualsi
Bij = Nsitesij /Nsitesj
IndValij = Aij × Bij × 100
式中:Aij反映了物种i 对样地组j 的专一性，是物种i 在样地组j 的平均多度( 个体数) 与该物种在所有样地组的平均多度之和的比值;
Nindividualsij表示物种i 在样地组j 的平均个体数，
而Nindividualsi是物种i 在所有样地组的平均个体数之和。
Bij反映了物种i 对样地组j 的忠实性，是样地组j 中存在物种i 的样地数与该样地组的所有样地数的比值;
Nsitesij表示样地组j 中存在物种i的样地数，
Nsitesj表示样地组j 中的所有样地数。
# Dufrene & Legendre物种指示值
# ****************************
# 依据das（离源头距离）环境变量将样方分为4组
das.D1 <- dist(data.frame(das=env[,1], row.names=rownames(env)))
dasD1.kmeans <- kmeans(das.D1, centers=4, nstart=100)
dasD1.kmeans$cluster
# 样方组的指示种
(iva <- indval(spe, dasD1.kmeans$cluster) )
#输出结果包括以下表格：
#-relfrg=物种在每个组的相对频度
#-relabu=物种在组间的相对多度
#-indval=每个物种的指示值
#下面两项内容将从indval表格内提取出来（含有最高指示值的组）和置换检
#验的结果。
# 显著指示物种的表格
gr <- iva$maxcls[iva$pval<=0.05]
iv <- iva$indcls[iva$pval<=0.05]
pv <- iva$pval[iva$pval<=0.05]
fr <- apply(spe>0, 2, sum)[iva$pval<=0.05]
fidg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fidg <- fidg[order(fidg$group, -fidg$indval),]
fidg
group indval pvalue freq
CHA 1 0.8000000 0.001 8
BLA 1 0.6300000 0.002 8
OMB 1 0.6125000 0.003 8
VAI 1 0.5075922 0.001 20
LOC 1 0.4269231 0.009 24
TRU 2 0.4560811 0.018 17
PCH 4 0.8571429 0.001 7
ANG 4 0.7346939 0.002 11
BCO 4 0.7034483 0.002 9
BBO 4 0.6936416 0.001 10
GRE 4 0.6857143 0.001 12
CAR 4 0.6818182 0.001 12
PSO 4 0.6623377 0.001 13
BOU 4 0.6551724 0.003 11
ROT 4 0.5919283 0.002 11
BRO 4 0.5454545 0.003 18
TAN 4 0.5351171 0.007 17
SPI 4 0.5238095 0.008 12
ABL 4 0.5067568 0.012 14
BAR 4 0.4951456 0.007 14
GOU 4 0.4754717 0.010 20
GAR 4 0.4698206 0.011 18
PER 4 0.4581498 0.024 15
CHE 4 0.3836164 0.033 25https://blog.sciencenet.cn/blog-1835014-1127449.html
上一篇：ggplot2调整X轴标签顺序
下一篇：数量生态学||多元回归树
