模糊c均值聚類（FCM）方法概述
模糊c均值聚類（Fuzzy C-Means Clustering，FCM）是最廣泛使用的模糊聚類演算法之一，與k均值劃分聚類較為相似，透過不斷迭代使組內平方和最小化。二者最主要的區別在於，k均值劃分聚類中每個物件只由一個聚類中心約束，模糊c均值聚類中物件與所有聚類中心都有關。
組內平方和計算如下：
其中uijm透過下式計算：
uij是物件xi屬於簇cj的程度，與物件xi離cj質心（聚類中心）的距離成反比；uj是簇j的質心。模糊c均值聚類中，一個物件理論上可以被分配為所有組，並透過給定物件在各組中的成員值（membership value）描述物件屬於某一類的程度（可以理解為機率）。成員值介於0和1之間，接近0代表了物件遠離該聚類中心，接近1代表物件靠近該聚類中心，每個物件在各類中成員值的總和為1。
聚類中心計算為組內所有物件的質心，並根據物件的成員度進行加權。
m是模糊引數，理論上可以取大於或等於1的任何值，當接近於1時聚類結果將與k均值劃分等硬聚類相似，值接近無限將導致完全模糊，通常使用m=2。
模糊c均值聚類過程可以總結如下：
（1）指定聚類數k；
（2）為每個物件任意分配一組初始成員值；
（3）使用上述公式計算每個聚類簇的質心；
（4）使用上述公式重新計算每個物件在各聚類簇中的成員值；
（5）重複（3）、（4）過程，直到成員值收斂（不再變動）或達到最大迭代次數為止。
R包cluster的模糊c均值聚類
接下來展示R語言執行模糊c均值聚類的方法。
示例資料和R程式碼的百度盤連結：
https://pan.baidu.com/s/1KjnCaP4Gfdys9lSrw5mUUA
cluster包提供了模糊c均值聚類方法，透過函式fanny()實現。
library(cluster)
#示例資料包含 15 個樣本（物件），20 個變數
dat <- read.delim('data.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
#推薦變數標準化，主要用於消除量綱差異，降低極端值比重
dat <- scale(dat)
#已知這些樣本來自於 3 組資料，但不明確具體哪些樣本是一組的，用來演示模糊 c 均值聚類過程
#模糊 c 均值聚類，詳情 ?fanny
set.seed(123)
cm.fanny <- fanny(dat, k = 3, metric = 'euclidean', maxit = 100)
cm.fanny
#或者也可計算物件間距離，將距離作為模糊 c 均值聚類的輸入，不使用原始資料
dis_euc <- vegan::vegdist(dat, method = 'euclidean')
set.seed(123)
cm.fanny <- fanny(dis_euc, k = 3, diss = TRUE, maxit = 100)
cm.fanny
#檢視主要結果
names(cm.fanny)
head(cm.fanny$membership) #物件屬於某一類的程度，即成員值，行和為 1
cm.fanny$clustering #最佳分類，即各物件最接近的聚類簇
#輪廓圖
plot(silhouette(cm.fanny))
輪廓圖顯示物件分類為3組的情況良好。
#以 PCoA 為例降維，並將聚類結果標註在排序圖中
#以上述所得物件間歐幾里得距離為例，計算 PCoA
pcoa <- cmdscale(dis_euc, k = (nrow(dat) – 1), eig = TRUE)
eig_prop <- 100*pcoa$eig/sum(pcoa$eig)
pcoa_site <- pcoa$point[ ,1:2]
plot(pcoa_site,
xlab = paste('PCoA axis1:', round(eig_prop[1], 2), '%'),
ylab = paste('PCoA axis2:', round(eig_prop[2], 2), '%'))
#標註各物件最接近的聚類簇
for (i in 1:3) {
gg <- pcoa_site[cm.fanny$clustering == i, ]
hpts <- chull(gg)
hpts <- c(hpts, hpts[1])
lines(gg[hpts, ], col = 1+1)
}
#星圖展示了各物件的成員值
stars(cm.fanny$membership, location = pcoa_site, draw.segments = TRUE,
add = TRUE, scale = FALSE, col.segments = 2:4, len = 0.5)
星圖展示了物件屬於各聚類簇的成員值，面積越大代表成員值越高，並將各物件所屬的最佳聚類簇用紅線連線。資料集整體特徵一目瞭然。
相比上述plot()作圖，factoextra包提供了更便捷的方法，以下是一些參考。
#factoextra 包的視覺化方案
library(factoextra)
#聚類和 PCA 結合，點的顏色和形狀代表分組，分組橢圓展示為 95% 置信區間
fviz_cluster(cm.fanny, ellipse.type = 'norm', repel = TRUE, palette = 'jco',
ellipse.level = 0.95, ggtheme = theme_minimal(), legend = 'right')
#輪廓圖
fviz_silhouette(cm.fanny, palette = 'jco', ggtheme = theme_minimal())
R包e1071的模糊c均值聚類
e1071包也是可用於模糊c均值聚類的常用R包，繼續使用上文的示例展示e1071包中的方法。
做過表達譜時間趨勢分析的同學們可能用過一個R包，Mfuzz，它就是以e1071包中的模糊c均值方法為基礎進行時間表達模式聚類的。
library(e1071)
#讀取資料
dat <- read.delim('data.txt', sep = '\t', row.names = 1, stringsAsFactors = FALSE, check.names = FALSE)
#推薦標準化，主要用於消除量綱差異，降低極端值比重
dat <- scale(dat)
#模糊 c 均值聚類，詳情 ?cmeans
#期望分為 3 類，最大迭代 100 次
set.seed(123)
cm.cmeans <- cmeans(dat, centers = 3, iter.max = 100, dist = 'euclidean', method = 'cmeans')
cm.cmeans
#檢視主要結果
summary(cm.cmeans)
head(cm.cmeans$membership) #物件屬於某一類的程度，即成員值，行和為 1
cm.cmeans$cluster #最佳分類，即各物件最接近的聚類簇
可以結合一些其它包的視覺化方法，觀測物件分類詳情。
#使用 corrplot 包描述物件歸屬的成員值
library(corrplot)
corrplot(cm.cmeans$membership, is.corr = FALSE)
該圖代表了模糊c均值聚類計算的各物件在各類別中的歸屬程度，明確屬於某一聚類簇的物件在該簇中具有較高的成員值。
同樣地，再結合排序圖觀測各類別的邊界，將聚類結果標註在圖中。
#factoextra 包的視覺化方案，和 PCA 結合
#點的顏色和形狀代表分組，分組橢圓展示為 95% 置信區間
library(factoextra)
fviz_cluster(list(data = dat, cluster = cm.cmeans$cluster), palette = 'jco',
ellipse.type = 'norm', ellipse.level = 0.95, ggtheme = theme_minimal())
不同顏色和形狀的點代表物件屬於不同的類別，各組中心最大的點代表了聚類中心，顯示該示例資料的樣本分類還是很好的。
參考資料
