聚類分析4—環(huán)境數(shù)據(jù)來解釋 (數(shù)量生態(tài)學：R語言的應(yīng)用-第四章)
在這之前我們學習了聚類分析的基本概念呜魄、幾種計算層次聚類的方法师倔、進一步解讀和比較層次聚類結(jié)果以及非層次聚類构韵，這些聚類方法都是基于物種多度數(shù)據(jù)對樣方進行分組，當然這些聚類方法也可以用于其他類型數(shù)據(jù)趋艘，特別是環(huán)境數(shù)據(jù)疲恢，所以本次就是介紹用環(huán)境數(shù)據(jù)來進行聚類分析。
本次的內(nèi)容不多瓷胧，主要分為兩個部分：
- 用外部數(shù)據(jù)進行類型比較（方差分析途徑）
- 雙類型比較（列聯(lián)表分析）
1.1 加載所需的包和數(shù)據(jù)
#加載包和數(shù)據(jù)
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)
#加載函數(shù)
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")
#從聚類結(jié)果獲得二元差異矩陣的函數(shù)
grpdist <- function(X)
{
require(cluster)
gr <- as.data.frame(as.factor(X))
distgr <- daisy(gr,"gower")
distgr
}
#導(dǎo)入Doubs數(shù)據(jù)
load("Doubs.RData")
#剔除無物種數(shù)據(jù)的樣方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
latlong <- latlong[-8,] #經(jīng)緯度
1.2數(shù)據(jù)預(yù)處理
#先計算樣方之間的弦距離矩陣
spe.norm <- decostand(spe,"normalize")
spe.ch <- vegdist(spe.norm,"euc")
# 計算Ward最小方差聚類
spe.ch.ward <- hclust(spe.ch, method="ward.D2")
#Ward聚類同表型相關(guān)
spe.ch.ward.coph <- cophenetic(spe.ch.ward)
cor(spe.ch, spe.ch.ward.coph)
# 設(shè)定聚類組的數(shù)量
k <- 4
# 根據(jù)上面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
dev.new(
title = "輪廓-優(yōu)化分區(qū)",
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))
說明：這里數(shù)據(jù)處理方法，都是在聚類分析這一章學過的搓萧，不明白的同學杂数，可以在回頭去看看
這里用到的數(shù)據(jù)是非層次距離中得到的spech.ward.gk
2 用外部數(shù)據(jù)進行類型比較（方差分析途徑）
我們之前學習的主要是內(nèi)部的準則（例如輪廓法或其他聚類質(zhì)量指數(shù)）都是僅僅依賴物種數(shù)據(jù)，還不足以選擇最佳樣方聚類結(jié)果矛绘。選擇最終的聚類結(jié)果有時也需要基于生態(tài)學解釋耍休。生態(tài)學解釋可視為樣方聚類的外部驗證刃永。
- 利用外部獨立的解釋變量驗證聚類結(jié)果（作為響應(yīng)數(shù)據(jù)）可以用線性判別式分析（這個在后面我們會學習到）货矮。
- 可以將樣方分組當作因子對解釋變量（當作響應(yīng)變量值）進行方差分析，了解解釋變量在各組間是否有顯著差異斯够。
下面的我們將學習用樣方聚類簇為因子去對解釋變量進行方差分析囚玫。
步驟：
- 檢驗?zāi)骋画h(huán)境變量是否符合方差分析假設(shè)（殘差正態(tài)性和方差齊性）
- 用傳統(tǒng)的（參數(shù)估計）單因素方差分析或非參數(shù)的Kruskal-Wallis檢驗解釋變量在組間是否有顯著差異
盡管在方差分析中，是將物種組成數(shù)據(jù)獲得的聚類的分組結(jié)果作為解釋變量读规，但是從生態(tài)學角度去分析抓督，實際上是尋找環(huán)境因子對樣方的分組的解釋 。
基于最優(yōu)化的Ward聚類（分4組）的環(huán)境變量（為提高正態(tài)性進行某些簡單的轉(zhuǎn)化）箱線圖
#繪制基于最優(yōu)化的Ward聚類（分四組）的環(huán)境變量線箱線圖
with(env, {
dev.new(
title = "定量環(huán)境變量箱線圖",
noRStudioGD = TRUE
)
par(mfrow = c(2, 2))
boxplot(
sqrt(ele) ~ spech.ward.gk,
main = "海拔",
las = 1,
ylab = "sqrt(alt)",
col = (1:k) + 1,
varwidth = TRUE
)
boxplot(
log(slo) ~ spech.ward.gk,
main = "坡度",
las = 1,
ylab = "log(slo)",
col = (1:k) + 1,
varwidth = TRUE
)
boxplot(
oxy ~ spech.ward.gk,
main = "含氧量",
las = 1,
ylab = "oxy",
col = (1:k) + 1,
varwidth = TRUE
)
boxplot(
sqrt(amm) ~ spech.ward.gk,
main = "銨濃度",
las = 1,
ylab = "sqrt(amm)",
col = (1:k) + 1,
varwidth = TRUE
)
})
# 方差分析假設(shè)檢驗
with(env, {
# 殘差的正態(tài)性
shapiro.test(resid(aov(sqrt(ele) ~ as.factor(spech.ward.gk))))
shapiro.test(resid(aov(log(slo) ~ as.factor(spech.ward.gk))))
shapiro.test(resid(aov(oxy ~ as.factor(spech.ward.gk))))
shapiro.test(resid(aov(sqrt(amm) ~ as.factor(spech.ward.gk))))
#檢驗結(jié)果表明sqrt(alt)束亏、log(pen)铃在、oxy和sqrt(amm)的殘差是正態(tài)分布。嘗試為其他的環(huán)境變量尋找好的標準化
# 方差齊性
bartlett.test(sqrt(ele), as.factor(spech.ward.gk))
bartlett.test(log(slo), as.factor(spech.ward.gk))
bartlett.test(oxy, as.factor(spech.ward.gk))
bartlett.test(sqrt(amm), as.factor(spech.ward.gk))
#變量sqrt(alt)的方差不齊碍遍，所以參數(shù)檢驗的方差分析不適用
# 可檢驗變量的方差分析A
summary(aov(log(slo) ~ as.factor(spech.ward.gk)))
summary(aov(oxy ~ as.factor(spech.ward.gk)))
summary(aov(sqrt(amm) ~ as.factor(spech.ward.gk)))
# 海拔Kruskal-Wallis 檢驗
kruskal.test(ele ~ as.factor(spech.ward.gk))
})
注意：
- 在一系列分析的開頭使用with（）來避免在每次分析中重復(fù)輸入對象env的名稱定铜。這比使用attach（）和detach（）方便，因為當你的R控制臺中有幾個數(shù)據(jù)集怕敬，有些數(shù)據(jù)集碰巧有名字相同的變量時揣炕，attach（）可能會導(dǎo)致混淆。
- Shapiro檢驗的零假設(shè)是變量正態(tài)分布东跪；Bartlett檢驗的零假設(shè)是組間方差相等畸陡。因此鹰溜，對于這兩個檢驗，只有當p值大于顯著性水平丁恭，即p>0.05時接受零假設(shè)曹动，才能滿足方差分析的假設(shè)。
- 參數(shù)Bartlett檢驗對偏離正態(tài)分布很敏感牲览。對于非正態(tài)分布數(shù)據(jù)仁期，可以使用bartlett.perm.R的函數(shù)來計算參數(shù)、置換和自助法（即替換的置換）Bartlett檢驗竭恬。
以下可以使用作者編寫的通用函數(shù)跛蛋，執(zhí)行方差分析的多重比較和顯示帶有字母的環(huán)境變量分組后箱線圖多重比較結(jié)果。不同字母表示組間有顯著差異（按中位線遞減順序組）痊硕。
- 使用boxplert（）函數(shù)來執(zhí)行方差分析和LSD多重比較檢驗
- 使用boxplerk（）函數(shù)執(zhí)行Kruskal-Wallis 檢驗及其相應(yīng)的多重比較（兩者都使用“Holm”的p值校正）
dev.new(
title = "ANOVA and Kruskal-Wallis 檢驗",
noRStudioGD = TRUE
)
par(mfrow = c(2, 2))
#使用boxplert（）或boxplerk（）繪制事后檢測的結(jié)果
#使用boxplert（）函數(shù)畫矯正后的多重比較結(jié)果
with(env, {
boxplerk(
ele,
spech.ward.gk,
xlab = "",
ylab = "ele",
main = "海拔",
bcol = (1:k) + 1,
p.adj = "holm"
)
boxplert(
log(slo),
spech.ward.gk,
xlab = "",
ylab = "log(slo)",
main = "坡度",
bcol = (1:k) + 1,
p.adj = "holm"
)
boxplert(
oxy,
spech.ward.gk,
xlab = "",
ylab = "oxy",
main = "含氧量",
bcol = (1:k) + 1,
p.adj = "holm"
)
boxplert(
sqrt(amm),
spech.ward.gk,
xlab = "",
ylab = "sqrt(amm)",
main = "銨濃度",
bcol = (1:k) + 1,
p.adj = "holm"
)
})
基于上面這些分析和圖示赊级，能確定這組魚類群落的生態(tài)習性。
當然岔绸，我們也可以基于環(huán)境變量對樣方進行聚類（類似獲得生境類型的分組）理逊，然后通過指示種分析（以后會講）檢驗不同生境內(nèi)物種分布是否有差異。指示種分析過程中基于不同的生境類型物種需要逐個分析盒揉。因此晋被，需要考慮多個物種指示種分析時會產(chǎn)生多重檢驗的統(tǒng)計學問題。
另外刚盈，作為替代方案羡洛，之后第6章會提出基于排序的多元方法，也可以直接描述和檢驗物種-生境關(guān)系藕漱。請期待欲侮。
3 雙類型比較（列聯(lián)表分析）
要是想直接比較分別基于物種數(shù)據(jù)和環(huán)境數(shù)據(jù)的樣方聚類結(jié)果該怎么辦呢？
- 直接生成一個含兩種結(jié)果的列聯(lián)表
- 用列聯(lián)表Fishr精準檢驗比較兩種樣方聚類結(jié)果是否有顯著差異
# 基于環(huán)境變量的樣方聚類
env2 <- env[, -1]
env.de <- vegdist(scale(env2), "euc")
env.kmeans <- kmeans(env.de, centers = 4, nstart = 100)
env.kmeans.g <- env.kmeans$cluster
# 比較從物種和環(huán)境數(shù)據(jù)獲得聚類4組的結(jié)果
table(spe.kmeans.g, env.kmeans.g)
#兩種聚類結(jié)果是否相同肋联？
# 列聯(lián)表 Fisher精確檢驗
fisher.test(table(spe.kmeans.g, env.kmeans.g))
# 用卡方檢驗分析兩種聚類之間的差異
chisq.test(table(spe.kmeans.g, env.kmeans.g))
#改用置換的方法進行卡方檢驗分析
chisq.test(table(spe.kmeans.g, env.kmeans.g),
simulate.p.value = TRUE)
可以上面的列聯(lián)表分析得出威蕉，物種和環(huán)境數(shù)據(jù)獲得聚類4組的結(jié)果具有顯著差異
同時列聯(lián)表分析同樣適用于比較分別基于物種數(shù)據(jù)和分類（定性）解釋變量數(shù)據(jù)的樣方聚類結(jié)果。
用環(huán)境數(shù)據(jù)進行比較的內(nèi)容就是這些橄仍，雖然不是很多韧涨，但是要聯(lián)系之前學習的才能更好的掌握它，所以還是有難度的侮繁，主要是用外部數(shù)據(jù)進行類型比較（方差分析途徑）和雙類型比較（列聯(lián)表分析兩部分內(nèi)容虑粥，好好學習掌握他。
這些案例的源代碼我都已經(jīng)上傳到getee鼎天，可以在我的公眾號回復(fù)：數(shù)量生態(tài)學獲得倉庫鏈接
謝謝你的閱讀舀奶，請期待下一期數(shù)量生態(tài)學：R語言的應(yīng)用 第四章 聚類分析5—聚類物種集合
如有不足或錯誤之處，請批評指正斋射。
有什么不明白的也歡迎留言討論育勺。
歡迎關(guān)注微信公眾號：fafu 生信 小蘑菇
往期內(nèi)容：
《數(shù)量生態(tài)學：R語言的應(yīng)用》第三章-R模式
《數(shù)量生態(tài)學：R語言的應(yīng)用》第二版第三章-關(guān)聯(lián)測度與矩陣------Q模式
《數(shù)量生態(tài)學：R語言的應(yīng)用》第二版筆記2
《數(shù)量生態(tài)學——R語言的應(yīng)用》第二版閱讀筆記--緒論和第二章（一部分）
感謝你的閱讀！＝е痢腹躁！你的點贊關(guān)注轉(zhuǎn)發(fā)是對我最大的鼓勵。
