數(shù)量生態(tài)學：R語言的應用 第四章 聚類分析3—非層次聚類
在聚類分析中層次聚類被經(jīng)常使用，層次聚類通過某種相似性測度計算節(jié)點之間的相似性，并按相似度由高到低排序，逐步重新連接個節(jié)點。而另一類非層次聚類是對一組對象進行簡單分組的方法，盡量使組內(nèi)對象之間比組間對象之間的相似度更高。但是需要我們自己來決定分組的組數(shù)k。
通常是設定不同的初始結構，然后通過大量的迭代以找到最佳的解決方法。本次介紹兩種非層次聚類算法：k-均值劃分和圍繞中心點劃分（PAM）。這些都是歐氏距離空間屬性的算法。
注意：不同綱量的變量在進行非層次聚類之前應該進行標準化，因為方差的綱量等于變量綱量的平方，如果不標準化，不同量綱的方差相加無意義。
1. k-均值劃分
k-均值劃分使用數(shù)據(jù)局部結構構建聚類簇：通過確認數(shù)據(jù)高密度區(qū)構建分類組。
其步驟是：預估將數(shù)據(jù)分為K組,則隨機選取K個對象作為初始的聚類中心,然后計算每個對象與各個種子聚類中心之間的距離,把每個對象分配給距離它最近的聚類中心 ,聚類中心以及分配給它們的對象就代表一個聚類 ,每分配一個樣本,聚類的聚類中心會根據(jù)聚類中現(xiàn)有的對象被重新計算,這個過程將不斷重復直到滿足某個終止條件 ,終止條件可以是沒有(或最小數(shù)目)對象被重新分配給不同的聚類,沒有(或最小數(shù)目)聚類中心再發(fā)生變化,誤差平方和局部最小。
k-均值劃分算法以歐式距離作為相似度測度，采用誤差平方和準則函數(shù)作為聚類準則函數(shù)。
加載所需的包和數(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")
#從聚類結果獲得二元差異矩陣的函數(shù)
grpdist <- function(X)
{
require(cluster)
gr <- as.data.frame(as.factor(X))
distgr <- daisy(gr,"gower")
distgr
}
#導入Doubs數(shù)據(jù)
load("Doubs.RData")
#剔除無物種數(shù)據(jù)的樣方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
latlong <- latlong[-8,] #經(jīng)緯度
#先計算樣方之間的弦距離矩陣
spe.norm <- decostand(spe,"normalize")
1.1 隨機開始的K-均值劃分
在已經(jīng)想好設定的組數(shù)后，使用stats包的 kmeans（）函數(shù)運行k-均值劃分。kmeans（）函數(shù)會以隨機設定的初始結構為基礎，不斷的自動重復迭代到設定次數(shù)（參數(shù)nstart）為止，并獲得當前次數(shù)下最佳的分類方案（即最小SSE值）。
k-均值劃分是一種線性模型的方法，因此不適合含很多零值的原始數(shù)據(jù)。遇到這樣的數(shù)據(jù)，有兩種方法進行處理。
- 使用非歐氏相異矩陣，如百分數(shù)差異（又名Bray-Curtis）；這樣的非歐氏相異矩陣應進行平方根變換并提交主坐標分析（PCoA）以獲得歐氏屬性的主坐標用于k-均值劃分。
- 對數(shù)據(jù)進行預轉化。為了與前面聚類分析的對象保持一致，可以使用之前已經(jīng)范數(shù)標準化后（弦轉化）的物種數(shù)據(jù)作為案例。為了與之前的Ward聚類結果比較，設定k=4組進行k-均值劃分，并將結果與Ward層次聚類4組樣方的結果進行比較。
為了與之前的Ward聚類結果比較，設定k=4組進行k-均值劃分，并將其結果與ward結果進行比較。
#預轉化后物種數(shù)據(jù)k-均值劃分
spe.kmeans <- kmeans(spe.norm, centers = 4, nstart = 100)
spe.kmeans
注意：即使給定的nstart相同，每次運行上述命令，所產(chǎn)生的結果也不一定完全相同，因為每次運算設定的初始結構是隨機的。
# 比較當前分4組的分類結果與之前Ward聚類的結果。
table(spe.kmeans$cluster, spech.ward.g)
#這兩個聚類結果是否非常相似？哪個（或哪些）對象有差別？
spe.kmeans$cluster
spech.ward.g
對于那些無法從原始數(shù)據(jù)轉化后進行歐氏距離計算的距離指數(shù)（例如百分數(shù)差異指數(shù)，又名Bray-Curtis相異指數(shù)），必須先將原始的距離矩陣通過主坐標分析（PCoA）獲得一個n行的矩形數(shù)據(jù)表格，然后再進行k-均值劃分。
對于百分數(shù)差異相異指數(shù)，必須先取平方根進行PCoA分析，或使用校對負特征根的PCoA函數(shù)，才能獲得完全的歐氏屬性數(shù)據(jù)表格。
kmeans（）函數(shù)每次分析只能產(chǎn)生一個簡單的預先設定組數(shù)的分組結果。如果需要嘗試不同組數(shù)的結果，需要重新運行命令。但到底多少組是最好的分類方案呢？
最好方案的標準有很多，cclust程序包clustIndex（）函數(shù)內(nèi)有一部分標準。推薦使用最大化Calinski-Harabasz指數(shù)（比較分類結果組間與組內(nèi)平方和的F統(tǒng)計量），盡管在分類組所含對象數(shù)量差異比較大情況下Calinski-Harabasz指數(shù)的值可能比較低?！皊si”（簡單結構指數(shù)）也是衡量分類結果的另一個不錯的指標。
我們不必手動多次運行kmeans（）函數(shù)來獲得不同組數(shù)的分類結果，可以使用vegan包里的函數(shù)cascadeKM（）。cascadeKM（）函數(shù)可以一次運行生成組數(shù)k從少（參數(shù)Inf.gr）到多（參數(shù)sup.gr）的聚類結果。
以下使用cascadeKM（）函數(shù)，組數(shù)設定從2組到10組，并使用ssi評估聚類的質(zhì)量，同時繪制聚類的結果。
#k-均值劃分，2組到10組
spe.KM.cascade <-
cascadeKM(
spe.norm,
inf.gr = 2,
sup.gr = 10,
iter = 100,
criterion = "ssi"
)
summary(spe.KM.cascade)
spe.KM.cascade$results
spe.KM.cascade$partition
dev.new(
title = "CascadeKM",
width = 10,
height = 6,
noRStudioGD = TRUE
)
#在plot函數(shù)中，sortg=TRUE代表在每個聚類簇內(nèi)按照對象之間的緊密程度重新排列對象。
plot(spe.KM.cascade, sortg = TRUE)
該圖顯示每個對象在每種分類組數(shù)下的歸屬（圖上每行代表一種組數(shù)）。圖內(nèi)的表格有不同的顏色，每行兩種顏色，代表分兩組k=2，三種顏色代表k=3，依此類推。右圖代表不同k值條件下的中止標準的統(tǒng)計量。
到底多少組數(shù)是最佳方案？如果傾向于較大的組數(shù)，哪個是最佳方案呢？
函數(shù)casecadeKM()也提供數(shù)據(jù)結果。其中“result”給出每種組數(shù)k條件下的TESS統(tǒng)計量和準則(calinski或ssi)的值，其中partition包含每個聚類簇所含對象的信息。如果對象的地理信息可用，就可以繪制對象空間分布地圖，同時以不同的顏色帶標識對象的歸屬。
summary(spe.KM.cascade)
spe.KM.cascade$results
這里，最小SSE值用于確定在給定組數(shù)k下的最佳方案，而calinski和ssi指標用于確定最佳k值，兩個指標解決不同的問題。
記住不同的組數(shù)k={2,3,…，10}，每次都是獨立運行。因此上圖從下到上結構互相獨立，與層次聚類樹的嵌套結構不同。
確定樣方聚類簇后，可以查看聚類簇的內(nèi)容。最簡單的方法是基于分組圖形定義樣方的亞組合計算基本統(tǒng)計量。以下對4組k-均值劃分結果進行分析。
# 按照k-均值劃分結果重新排列樣方
spe.kmeans.g <- spe.kmeans$cluster
spe[order(spe.kmeans.g), ]
# 使用函數(shù)vegemite（）重新排列樣方-物種矩陣
ord.KM <- vegemite(spe, spe.kmeans.g)
spe[ord.KM$sites, ord.KM$species]
1.2 使用k-均值劃分優(yōu)化一個獨立獲得的分類
k-均值劃分的出發(fā)點可以以k x p的平均值矩陣（k為聚類獲得的組數(shù)，p為變量數(shù)量，矩陣里面的值為該變量在該組中的平均值），或是以k個對象（每個對象視為每個組的典型代表）作為列表，這個典型的代表即下一節(jié)要討論的分組方法所稱的“形心（medeids）”。
讓我們應用這個想法來驗證由魚類數(shù)據(jù)Ward聚類獲得的4個組，然后被k-均值修改分組。
以Ward 聚類獲得4組中各個物種的平均多度矩陣作為初始結構進行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)
與前面的方法有點不同，這里需要返回看一下聚類（層次結構聚類）的結果，確認每組最典型的對象（輪廓寬度值最大），將這些對象作為k-均值劃分的初始結構（參數(shù)centers）
startobjects <- spe.norm[c(2, 17, 21, 23), ]
spe.kmeans3 <- kmeans(spe.norm, centers = startobjects)
#將k-均值劃分結果與Ward聚類分析結果進行比對
table(spe.kmeans2$cluster, spech.ward.g)
# 兩個k-均值劃分優(yōu)化后的4組進行比較
table(spe.kmeans2$cluster, spe.kmeans3$cluster)
# 最終分組的輪廓寬度值圖
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))
從代碼中可以看出，可以通過列聯(lián)表法進行原始分組和優(yōu)化分組之間或兩個優(yōu)化分組之間的比較。兩種優(yōu)化（基于物種的平均值和典型的初始對象）為我們選擇的4個“典型”對象產(chǎn)生相同的結果。注意這個選擇至關重要。與原來的Ward聚類4組結果比較顯示，只有一個對象第15樣方已從組2移動到組1。這樣移動可以改進輪廓圖并略微改變了沿著河流分布的地圖。請注意，樣方19的成員歸屬仍不清楚：它的輪廓寬度值仍為負值。
# 沿DOUBS河流優(yōu)化后Ward聚類4組分布圖
dev.new(title = "河流優(yōu)化后Ward聚類4組分布圖",
width = 9,
noRStudioGD = TRUE)
drawmap(xy = spa,
clusters = spech.ward.gk,
main = "沿DOUBS河流優(yōu)化后Ward聚類4組分布圖")
2. 圍繞中心點劃分（PAM）
圍繞中心點劃分：“從所有的數(shù)據(jù)觀測點尋找k個代表性的對象或形心點，這些代表性的對象應該反映數(shù)據(jù)的主體結構。k個形心點選定后，將每個觀測點分配給某個形心點構建k個聚類簇，不斷尋找最佳的k個代表性對象，使對象之間的相異性總和最小”。k-均值法是最小化組內(nèi)歐氏距離平方和。因此，k-均值法是傳統(tǒng)的最小二乘法，但PAM不是。
在R里，pam（）函數(shù)（cluster程序包）的輸入可以是原始數(shù)據(jù)，也可以是相異矩陣（與kmean（）相比，pam（）的優(yōu)勢是可以輸入更多類型的關聯(lián)測度），并且允許通過輪廓寬度值確定最佳的分組數(shù)量。
以下代碼兩的最后部分利用雙輪廓圖比較k-均值法和PAM法的分組結果。
基于弦距離的圍繞中心點劃分（PAM）
#聚類簇數(shù)量的選擇
#循環(huán)計算2至28組分類數(shù)下平均輪廓寬度
asw <- numeric(nrow(spe))
for (k in 2:(nrow(spe) - 1))
asw[k] <- pam(spe.ch, k, diss = TRUE)$silinfo$avg.width
k.best <- which.max(asw)
dev.new(title = "PAM", noRStudioGD = TRUE)
plot(
1:nrow(spe),
asw,
type = "h",
main = "聚類簇數(shù)量的選擇",
xlab = "k (群集數(shù))",
ylab = "平均輪廓寬度"
)
axis(
1,
k.best,
paste("optimum", k.best, sep = "\n"),
col = "red",
font = 2,
col.axis = "red"
)
points(k.best,
max(asw),
pch = 16,
col = "red",
cex = 1.5)
當k=2時，PAM具有最佳的方案（asw=0.3841），這并不是我們期望的結果。如果選擇常用的4組分法，從輪廓寬度角度分析結果并不好（asw=0.2736）。盡管如此，我們還是需要分析PAM分4組的情況。
# PAM分4組情況
spe.ch.pam <- pam(spe.ch, k = 4, diss = TRUE)
summary(spe.ch.pam)
spe.ch.pam.g <- spe.ch.pam$clustering
spe.ch.pam$silinfo$widths
# 將當前的分類結果與之前Ward聚類和k-均值劃分進行比較
table(spe.ch.pam.g, spech.ward.g)
table(spe.ch.pam.g, spe.kmeans.g)
PAM結果與Ward聚類和k-均值劃分結果顯著不同
#k=4組下與Ward聚類和k-均值劃分進行比較
dev.new(
title = "輪廓寬度值 - k-均值 and PAM",
width = 12,
height = 8,
noRStudioGD = TRUE
)
par(mfrow = c(1, 2))
k <- 4
sil <- silhouette(spe.kmeans.g, spe.ch)
rownames(sil) <- row.names(spe)
plot(sil,
main = "輪廓寬度值 - k-均值",
cex.names = 0.8,
col = 2:(k + 1))
plot(
silhouette(spe.ch.pam),
main = "輪廓寬度值 - PAM",
cex.names = 0.8,
col = 2:(k + 1)
)
基于此圖，可以分辨PAM和k-均值法哪個有更好的平均輪廓值。我們可以將當前結果與Ward聚類比較輪廓寬度圖比較。除了聚類簇成員不同外，在k-均值和最優(yōu)化的Ward聚類之間還有哪些不同?我們可以通過檢查兩個聚類比較的列聯(lián)表進行分析
#比較k-均值劃分和最優(yōu)化后的Ward聚類
table(spe.kmeans.g, spech.ward.gk)
PAM法表現(xiàn)更穩(wěn)定，一方面因為它是最小化相異系數(shù)總和，代替了最小化歐氏距離的平方和；另一方面在給定分組數(shù)k的條件下，PAM相異系數(shù)總和更容易收斂。但對于特定的研究目的，并不能保證PAM一定是最合適的聚類方法。
從前面的例子，可以看出即使是兩種目標相同且屬于同一大類（k-均值法和PAM同屬于非層次聚類）的聚類方法也可能產(chǎn)生完全不同的結果。我們需要根據(jù)哪種方法的結果能夠提供更多有用的信息，或能更好地被環(huán)境變量解釋來選擇合適的方法。
好了，第四章聚類分析的非層次聚類就這么多內(nèi)容，主要就是k-均值法和PAM這兩種方法，如果還不能好好理解，可以多回顧幾次，下次內(nèi)容是聚類分析中的使用環(huán)境變量來解釋來比較，敬請期待。
謝謝你的閱讀。
如有不足或錯誤之處，請批評指正。
有什么不明白的也歡迎留言討論。
歡迎關注微信公眾號：fafu 生信 小蘑菇
往期內(nèi)容：
《數(shù)量生態(tài)學：R語言的應用》第三章-R模式
《數(shù)量生態(tài)學：R語言的應用》第二版第三章-關聯(lián)測度與矩陣------Q模式
《數(shù)量生態(tài)學：R語言的應用》第二版筆記2
《數(shù)量生態(tài)學——R語言的應用》第二版閱讀筆記--緒論和第二章（一部分）
感謝你的閱讀?。?！你的點贊關注轉發(fā)是對我最大的鼓勵。
