k均值劃分（k-means）
k均值劃分（k-means partitioning，k-means）從任意分配的k組資料開始，經迭代最佳化產生最終結果。
k-means演算法簡述
對於給定資料集（a），透過k-means聚類的計算方法可概述如下。
（b）首先隨機分配k個聚類中心，這裡設定k=3；
（c）分配每個物件到離它最近的聚類中心，得到k個類；每個物件被分配到使下式得到最小值的那一類中：
xij表示第i個物件中第j個變數的值，`xkj表示第k個類中第j個變數的均值，n是物件個數，p是變數個數。
（d）計算每個類的質心，將該質心作為新的聚類中心；
（e）分配每個物件到離它最近的聚類中心，得到k個類，分類標準同（c）；
（f）重複（d）（e）過程，直到物件分類不再變動或達到最大迭代次數為止。
可以看到，k-means是一種線性方法，使用資料區域性結構構建聚類簇，把物件分成k組並使得物件到指定的聚類中心的平方總和為最小。它基於物件之間的歐幾里得距離，透過不斷迭代將總誤差平方和（TESS，或稱組內平方和）最小化，分組標準和Ward最小方差聚類相似。
k-means聚類在處理較大樣本量的資料集時，通常比層次聚類更有優勢。但是均值的使用意味著所有的變數必須是連續的，並且受異常值影響較大。物件聚類與初始指定k的數量密切相關，且對初始中心值的選擇也很敏感。
R語言執行k-means
在R中，可使用kmeans()執行k-means聚類。
示例資料和R程式碼的獲取連結：
https://pan.baidu.com/s/1WEHBd2gSzA9Q3s99qhEBFg
一個簡單示例
示例資料集GeneExpExample5000，記錄了某項研究中獲得的腎臟（Kidney）和肝臟（Liver）RNA樣本的基因表達值資訊，從中取了一部分資料作為示例資料集演示，共5000行（5000個基因的表達值）。
使用k-means的方法劃分樣本歸類。
#示例資料集，該資料集實際上來自 DEGseq 包，詳情可在載入 DEGseq 包後 ?GeneExpExample5000 檢視
#已知的先驗分組，Kidney_group=c(7, 9, 12, 15, 17, 18, 20)；Liver_group=c(8, 10, 11, 13, 14, 16, 19)
gene_express <- read.delim('GeneExpExample5000.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
reads_count <- t(gene_express[7:20])
#當變數量綱不同時，需要標準化資料使方差有意義；變數標準化還能降低極大值的影響
#本示例變數量綱一致，但如果覺得基因表達值的數量級相差過大可能導致較大離散度，可選標準化
#reads_count <- scale(reads_count)
#k-means 聚類，詳情 ?kmeans
#由於 k-means 在開始時隨機指定中心點，每次呼叫函式時可能獲得不同的方案，使用 set.seed() 保證可重複性
#設定 k=2，並生成 25 種初始配置（初始中心點）並輸出最好的一個
set.seed(123)
gene_kmeans <- kmeans(reads_count, centers = 2, nstart = 25)
summary(gene_kmeans)
gene_kmeans$cluster
#輪廓圖評估
plot(silhouette(gene_kmeans$cluster, vegan::vegdist(reads_count, method = 'euclidean')))
對於最佳聚類數量的選擇
上述示例是已知資料先驗分組的情況。
但有時，我們可能並不知道樣本來源，不知道聚為多少類是合適的。此時對於k值的選擇，R中也給出了一些參考標準。
#隨機生成模擬資料
set.seed(123)
rand <- rbind(matrix(rnorm(150, mean = 10, sd = 0.3), ncol = 5),
matrix(rnorm(150, mean = 3, sd = 0.2), ncol = 5),
matrix(rnorm(150, mean = 1, sd = 0.1), ncol = 5),
matrix(rnorm(150, mean = 6, sd = 0.3), ncol = 5),
matrix(rnorm(150, mean = 9, sd = 0.3), ncol = 5))
#NbClust 包的最佳聚類數量評估，詳情 ?NbClust
library(NbClust)
nc <- NbClust(rand, distance = 'euclidean', min.nc = 2, max.nc = 10, method = 'kmeans', index = 'hubert')
plot(nc)
模擬生成了5組均值和方差不等的資料，NbClust的方法給出評估，聚類為5組是最佳的（在cluster=5的位置出現明顯的拐點）。
#vegan 包的最佳聚類數量評估，詳情 ?cascadeKM
library(vegan)
cc <- cascadeKM(rand, inf.gr = 2, sup.gr = 10, iter = 100, criterion = 'ssi')
plot(cc)
同樣顯示聚類為5組是最佳的（k=5時，變數的結構最清晰，ssi值最大）。
至於最終聚類是否和模擬生成的5組資料完全對應，執行k=5時的k-means對比觀察即可，就不再演示了。
非歐式距離的應用
如果不期望使用歐式距離，或者只提供了非歐式距離的距離矩陣而缺少原始資料時，該方法可供參考。
例如期望透過物件間的Bray-curtis距離完成k-means。由於Bray-curtis距離並非歐氏距離，可以首先進行平方根轉化，使其具有歐式幾何屬性，然後將轉化後的距離應用主座標分析（PCoA），並將所有PCoA軸作為k-means的輸入，完成聚類。
#非歐式距離的應用，如 Bray-curtis 距離
#讀取距離矩陣
dis_bray <- read.delim('bray_distance.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
#平方根轉化
dis <- (as.dist(dis_bray))^0.5
#PCoA
pcoa <- cmdscale(dis, k = (nrow(dis_bray) – 1), eig = TRUE)
#獲取物件排序軸
site <- pcoa$point
#k-means，聚為 3 類
bray_kmeans <- kmeans(site, centers = 3, nstart = 25)
bray_kmeans$cluster
圍繞中心點劃分（PAM）
由於k-means聚類基於最小化組內歐式距離平方和，所以它對異常值敏感。相比之下，圍繞中心點劃分（Partitioning Around Mediods，PAM）是一個更穩健的方法。
PAM演算法簡述
（1）隨機選擇k個觀測物件作為聚類中心；
（2）計算每個聚類中心到每個物件的距離，可以為任意的距離測度；並分配每個物件到離它最近的聚類中心，得到k個類；
（3）計算每個聚類中心到該類中每個物件的距離的總和（Dsum）；
（4）在各類中，重新選擇一個物件作為新的聚類中心，並計算新的聚類中心到每個物件的距離；並分配每個物件到離它最近的聚類中心，得到新的k個類；
（5）計算新的每個聚類中心到該類中每個物件的距離的總和（D’sum）；
（6）如果D’sum<Dsum，則使用新的聚類中心；
（7）重複（4）–（6）過程，直到物件分類不再變動為止。
可以看到，PAM的目的是使組內物件之間的相異性總和最小，這種相異性可以為任意的距離測度，不再侷限於歐幾里得距離。這個特徵允許PAM容納混合資料型別，不僅限於連續變數，因此PAM在方法上比k-means更靈活。
R語言執行PAM
在R中，可使用cluster包pam()執行PAM聚類。
以上文示例為例繼續。
library(cluster)
#PAM 聚類，詳情 ?pam
#可以輸入原始資料，指定距離測度，同樣以歐幾里得距離為例
gene_pam <- pam(reads_count, k = 2, metric = 'euclidean')
summary(gene_pam)
#也可以首先計算距離矩陣，將距離矩陣作為輸入，同樣以歐幾里得距離為例
dis_euc <- vegan::vegdist(reads_count, method = 'euclidean')
gene_pam <- pam(dis_euc, k = 2, diss = TRUE)
summary(gene_pam)
#可透過輪廓圖評估上述分類是否可靠
plot(gene_pam)
