數(shù)量生態(tài)學：R語言的應(yīng)用 第四章 聚類分析3—非層次聚類
在聚類分析中層次聚類被經(jīng)常使用溪王，層次聚類通過某種相似性測度計算節(jié)點之間的相似性，并按相似度由高到低排序值骇，逐步重新連接個節(jié)點。而另一類非層次聚類是對一組對象進行簡單分組的方法移国，盡量使組內(nèi)對象之間比組間對象之間的相似度更高吱瘩。但是需要我們自己來決定分組的組數(shù)k。
通常是設(shè)定不同的初始結(jié)構(gòu)迹缀，然后通過大量的迭代以找到最佳的解決方法使碾。本次介紹兩種非層次聚類算法：k-均值劃分和圍繞中心點劃分（PAM）。這些都是歐氏距離空間屬性的算法祝懂。
注意：不同綱量的變量在進行非層次聚類之前應(yīng)該進行標準化票摇，因為方差的綱量等于變量綱量的平方，如果不標準化砚蓬，不同量綱的方差相加無意義矢门。
1. k-均值劃分
k-均值劃分使用數(shù)據(jù)局部結(jié)構(gòu)構(gòu)建聚類簇：通過確認數(shù)據(jù)高密度區(qū)構(gòu)建分類組。
其步驟是：預(yù)估將數(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")
#從聚類結(jié)果獲得二元差異矩陣的函數(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è)定的組數(shù)后，使用stats包的 kmeans（）函數(shù)運行k-均值劃分物延。kmeans（）函數(shù)會以隨機設(shè)定的初始結(jié)構(gòu)為基礎(chǔ)宣旱，不斷的自動重復迭代到設(shè)定次數(shù)（參數(shù)nstart）為止，并獲得當前次數(shù)下最佳的分類方案（即最小SSE值）叛薯。
k-均值劃分是一種線性模型的方法浑吟，因此不適合含很多零值的原始數(shù)據(jù)。遇到這樣的數(shù)據(jù)耗溜，有兩種方法進行處理组力。
- 使用非歐氏相異矩陣，如百分數(shù)差異（又名Bray-Curtis）强霎；這樣的非歐氏相異矩陣應(yīng)進行平方根變換并提交主坐標分析（PCoA）以獲得歐氏屬性的主坐標用于k-均值劃分侮措。
- 對數(shù)據(jù)進行預(yù)轉(zhuǎn)化。為了與前面聚類分析的對象保持一致瞳别，可以使用之前已經(jīng)范數(shù)標準化后（弦轉(zhuǎn)化）的物種數(shù)據(jù)作為案例忙灼。為了與之前的Ward聚類結(jié)果比較，設(shè)定k=4組進行k-均值劃分家夺，并將結(jié)果與Ward層次聚類4組樣方的結(jié)果進行比較脱柱。
為了與之前的Ward聚類結(jié)果比較，設(shè)定k=4組進行k-均值劃分拉馋，并將其結(jié)果與ward結(jié)果進行比較榨为。
#預(yù)轉(zhuǎn)化后物種數(shù)據(jù)k-均值劃分
spe.kmeans <- kmeans(spe.norm, centers = 4, nstart = 100)
spe.kmeans
注意：即使給定的nstart相同，每次運行上述命令煌茴，所產(chǎn)生的結(jié)果也不一定完全相同随闺，因為每次運算設(shè)定的初始結(jié)構(gòu)是隨機的。
# 比較當前分4組的分類結(jié)果與之前Ward聚類的結(jié)果蔓腐。
table(spe.kmeans$cluster, spech.ward.g)
#這兩個聚類結(jié)果是否非常相似矩乐？哪個（或哪些）對象有差別？
spe.kmeans$cluster
spech.ward.g
對于那些無法從原始數(shù)據(jù)轉(zhuǎn)化后進行歐氏距離計算的距離指數(shù)（例如百分數(shù)差異指數(shù)回论，又名Bray-Curtis相異指數(shù)）散罕，必須先將原始的距離矩陣通過主坐標分析（PCoA）獲得一個n行的矩形數(shù)據(jù)表格，然后再進行k-均值劃分傀蓉。
對于百分數(shù)差異相異指數(shù)欧漱，必須先取平方根進行PCoA分析，或使用校對負特征根的PCoA函數(shù)葬燎，才能獲得完全的歐氏屬性數(shù)據(jù)表格误甚。
kmeans（）函數(shù)每次分析只能產(chǎn)生一個簡單的預(yù)先設(shè)定組數(shù)的分組結(jié)果缚甩。如果需要嘗試不同組數(shù)的結(jié)果，需要重新運行命令靶草。但到底多少組是最好的分類方案呢蹄胰？
最好方案的標準有很多，cclust程序包clustIndex（）函數(shù)內(nèi)有一部分標準奕翔。推薦使用最大化Calinski-Harabasz指數(shù)（比較分類結(jié)果組間與組內(nèi)平方和的F統(tǒng)計量）裕寨，盡管在分類組所含對象數(shù)量差異比較大情況下Calinski-Harabasz指數(shù)的值可能比較低∨杉蹋“ssi”（簡單結(jié)構(gòu)指數(shù)）也是衡量分類結(jié)果的另一個不錯的指標宾袜。
我們不必手動多次運行kmeans（）函數(shù)來獲得不同組數(shù)的分類結(jié)果，可以使用vegan包里的函數(shù)cascadeKM（）驾窟。cascadeKM（）函數(shù)可以一次運行生成組數(shù)k從少（參數(shù)Inf.gr）到多（參數(shù)sup.gr）的聚類結(jié)果庆猫。
以下使用cascadeKM（）函數(shù)，組數(shù)設(shè)定從2組到10組绅络，并使用ssi評估聚類的質(zhì)量月培，同時繪制聚類的結(jié)果。
#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ù)結(jié)果举娩。其中“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}召边，每次都是獨立運行铺呵。因此上圖從下到上結(jié)構(gòu)互相獨立，與層次聚類樹的嵌套結(jié)構(gòu)不同隧熙。
確定樣方聚類簇后片挂，可以查看聚類簇的內(nèi)容。最簡單的方法是基于分組圖形定義樣方的亞組合計算基本統(tǒng)計量贞盯。以下對4組k-均值劃分結(jié)果進行分析音念。
# 按照k-均值劃分結(jié)果重新排列樣方
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）”件余。
讓我們應(yīng)用這個想法來驗證由魚類數(shù)據(jù)Ward聚類獲得的4個組讥脐，然后被k-均值修改分組。
以Ward 聚類獲得4組中各個物種的平均多度矩陣作為初始結(jié)構(gòu)進行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)
與前面的方法有點不同啼器，這里需要返回看一下聚類（層次結(jié)構(gòu)聚類）的結(jié)果旬渠，確認每組最典型的對象（輪廓寬度值最大），將這些對象作為k-均值劃分的初始結(jié)構(gòu)（參數(shù)centers）
startobjects <- spe.norm[c(2, 17, 21, 23), ]
spe.kmeans3 <- kmeans(spe.norm, centers = startobjects)
#將k-均值劃分結(jié)果與Ward聚類分析結(jié)果進行比對
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)生相同的結(jié)果。注意這個選擇至關(guān)重要更哄。與原來的Ward聚類4組結(jié)果比較顯示芋齿，只有一個對象第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個代表性的對象或形心點麻敌，這些代表性的對象應(yīng)該反映數(shù)據(jù)的主體結(jié)構(gòu)栅炒。k個形心點選定后，將每個觀測點分配給某個形心點構(gòu)建k個聚類簇术羔，不斷尋找最佳的k個代表性對象赢赊，使對象之間的相異性總和最小”。k-均值法是最小化組內(nèi)歐氏距離平方和级历。因此释移，k-均值法是傳統(tǒng)的最小二乘法，但PAM不是寥殖。
在R里玩讳，pam（）函數(shù)（cluster程序包）的輸入可以是原始數(shù)據(jù)涩蜘，也可以是相異矩陣（與kmean（）相比，pam（）的優(yōu)勢是可以輸入更多類型的關(guān)聯(lián)測度）熏纯，并且允許通過輪廓寬度值確定最佳的分組數(shù)量同诫。
以下代碼兩的最后部分利用雙輪廓圖比較k-均值法和PAM法的分組結(jié)果。
基于弦距離的圍繞中心點劃分（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）误窖，這并不是我們期望的結(jié)果。如果選擇常用的4組分法往扔，從輪廓寬度角度分析結(jié)果并不好（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
# 將當前的分類結(jié)果與之前Ward聚類和k-均值劃分進行比較
table(spe.ch.pam.g, spech.ward.g)
table(spe.ch.pam.g, spe.kmeans.g)
PAM結(jié)果與Ward聚類和k-均值劃分結(jié)果顯著不同
#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-均值法哪個有更好的平均輪廓值。我們可以將當前結(jié)果與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)生完全不同的結(jié)果闺骚。我們需要根據(jù)哪種方法的結(jié)果能夠提供更多有用的信息，或能更好地被環(huán)境變量解釋來選擇合適的方法妆档。
好了僻爽，第四章聚類分析的非層次聚類就這么多內(nèi)容，主要就是k-均值法和PAM這兩種方法贾惦，如果還不能好好理解胸梆，可以多回顧幾次，下次內(nèi)容是聚類分析中的使用環(huán)境變量來解釋來比較须板，敬請期待碰镜。
謝謝你的閱讀。
如有不足或錯誤之處习瑰，請批評指正绪颖。
有什么不明白的也歡迎留言討論。
歡迎關(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ā)是對我最大的鼓勵滓鸠。
