本周開始我們的《數量生態(tài)學筆記》的第四章：聚類分析。聚類分析又稱群分析唁影，它是研究（樣品或指標）分類問題的一種統(tǒng)計分析方法耕陷，同時也是數據挖掘的一個重要算法。在生態(tài)學研究當中夭咬，聚類的目的是識別環(huán)境中不連續(xù)對象的子集啃炸。實際上，聚類分析是所研究對象集合的分組卓舵。
需要注意大部分聚類方法都是基于關聯矩陣進行計算南用，這也說明選擇恰當的關聯系數非常重要。
如圖所示掏湾，我們需要識別不同類型的聚類方法及其應用的條件裹虫。
# 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\\數量生態(tài)學\\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="單連接聚合聚類")
#基于第一次聚類的結果拇涤，如何描述這個數據集捣作？是簡單的單一梯度還是區(qū)
#分明顯的樣方組？能否辨認樣方的連接鏈鹅士？樣方1券躁、5和9為什么最后連
#接？
完全連接聚合聚類
允許一個對象或簇與另一組聚合的依據是最遠距離對。
# 計算完全連接聚合聚類
# ********************
spe.ch.complete <- hclust(spe.ch, method="complete")
plot(spe.ch.complete,main="聚類樹",ylab="高度",xlab="完全連接聚合聚類")
#當前所給的樣方是沿著河流分布（樣方的編號按照流向編排）也拜，這個聚類分
#析結果是否將位置相近的樣方排在同一個組呢以舒？
#兩種完全有效的聚類分析方法分析同一數據，為什么產生如此不同的聚類
#結果呢慢哈？
單一連接是一個對象很容易聚合到一個分組蔓钟，因為單次連接足以導致融合。因此單連接聚類也被稱為最親密朋友法岸军。雖然產生的分類組不清晰奋刽，但容易識別梯度瓦侮。相反艰赞，完全連接聚類產生的分類之間差異比較明顯。完全連接聚類更傾向與產生很多小的分離的組肚吏，更適合于尋找和識別數據的間斷分布方妖。
平均聚合聚類
平均聚合聚類是一類基于對象平均相異性或聚類簇中心的聚類方法。此類聚類有4種罚攀，不同的方法區(qū)別在于組的位置計算方式和計算融合時是否包含對象數量作為權重党觅。
| 權重 | 算術平均 | 形心聚類 |
|---|---|---|
| 等權重 | UPGMA（average） | UPGMC（centeoid） |
| 不等權重 | WPGMA(mcquitty) | WPGMC(median) |
最著名的當屬UPGMA方法，一個對象加入一個組的依據是這個對象與該組的每個成員之間的平均距離斋泄。
# 計算UPGMA聚合聚類
# ***********************
spe.ch.UPGMA <- hclust(spe.ch, method="average")
plot(spe.ch.UPGMA,main="聚類樹",ylab="高度",xlab="UPGMA聚合聚類")
#這個UPGMA聚合聚類樹看起來介于單連接聚類和完全連接聚類之間杯瞻。這種
#情況經常發(fā)生。
需要注意的是 UPGMC和WPGMC有時會導致樹的翻轉炫掐，分類結果難以解讀魁莉。
# 計算魚類數據的形心聚類
# ***********************
spe.ch.centroid <- hclust(spe.ch, method="centroid")
plot(spe.ch.centroid,main="聚類樹",ylab="高度",xlab="UPGMC聚合聚類")
#這種聚類樹對生態(tài)學家來說簡直是噩夢。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聚類")
#使用距離平方造成此聚類樹上半部分過于膨脹。為了使聚類樹比例看起來
#更協(xié)調而不影響結構祷嘶，可以使用當前融合水平的平方根重新繪圖（圖4.5）
spe.ch.ward$height <- sqrt(spe.ch.ward$height)
plot(spe.ch.ward,main="聚類樹",ylab="高度",xlab="Ward聚類")
解讀和比較聚類分析結果
需要牢記的是聚類分析是一種探索分析屎媳，而非統(tǒng)計檢驗。影響聚類結果的因素包括聚類方法本省和用于聚類分析的關聯系數论巍。
同表型相關
對已經完成層次聚類任意兩個對象烛谊，在聚類樹上從一個對象向上走，到達與另一個對象交回節(jié)點向下走环壤，勢必會到達第二個對象晒来。交匯節(jié)點所在的層次水平即是兩個對象的同表型距離 。
# 單連接聚類同表型相關
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 (節(jié)點高度）)", 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 (節(jié)點高度）", 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 (節(jié)點高度）", 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 (節(jié)點高度）", 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個組里。
確定分組數
輪廓寬度圖
輪廓寬度是描述一個對象與所屬聚類歸屬程度的測度驶乾，是一個對象同組內其他對象的平均距離與該對象同最鄰近聚類簇內所有對象的平均距離比較邑飒。
# 依據輪廓寬度圖選擇最優(yōu)化的聚類簇數量（Rousseeuw質量指數）
# ********************************************************
# 繪制所有分類水平（除了k=1組的情況）輪廓寬度值（Ward 聚類）
# 首先產生一個長度等于樣方數量的空向量asw
asw <- numeric(nrow(spe))
# 其次循環(huán)獲得輪廓寬度值并依次填入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="輪廓寬度-最優(yōu)聚類簇數（Ward聚類）",
xlab="k (組數）", ylab="平均輪廓寬度")
axis(1, k.best, paste("最優(yōu)",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
# 屏幕上將輸出如下內容：
# 輪廓寬度最優(yōu)的聚類簇數量k=2，此時平均輪廓寬度為0.365819
#輪廓寬度法經常選擇2組作為最優(yōu)的分類數量级乐。然而疙咸，從生態(tài)學角度分析
#分4組似乎更合理。
距離矩陣與代表分組的二元矩陣的比較
# 依據Mantel統(tǒng)計（Pearson）相關選擇最優(yōu)分組數量
# **********************************************
# 編寫計算代表分類水平的二元距離矩陣的函數
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-最優(yōu)聚類簇數（Ward聚類）",
xlab="k (組數）", ylab="Pearson 相關")
axis(1, k.best, paste("最優(yōu)", 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="樣方")
參考：
常用聚類算法有哪些？六大類聚類算法詳細介紹
無監(jiān)督學習-- 聚類（Clustering）
百科||聚類算法
聚類分析
