本期由吳濤師弟講解聚類，內容很詳實，推薦感興趣的讀者通過原文連結觀看介紹視頻。
❞聚類分析的思想:對於有p個變量的數據集來說,每個觀測值都是p維空間中的一個點,所以屬於同一類的點在空間中的距離應該顯著小於屬於不同類的點之間的距離
聚類距離測度1.歐氏(Euclidean)距離：
2.曼哈頓(Manhattan)距離：
x,y都是n維的向量
還有一些不相似性度量也可以來表示距離：
3.Pearson線性相關距離：也就是1減去Pearson相關係數，Pearson相關是衡量兩個變量的線性相關性的，對數據的假設：Pearson相關係數應用於連續變量，假定兩組變量均為正態分布、存在線性關係且等方差
4.cosine相關距離:是Pearson相關的一個特例，就是將
5.Spearman相關距離：spearman相關計算的是變量秩之間的相關性,也是1減去Spearman相關係數
6.Kendall 相關距離：
Kendall相關方法是衡量變量秩的correspondence
對於大小是n的變量x和y,可能的匹配對數是
距離的選擇如果我們關注的是變量的變化趨勢而不是變量的具體的值的時候，比如基因的表達，我們就可以使用基於相關性的距離測度，另外要注意的是pearson相關對outliers比較敏感,這個時候可以使用spearman相關
當我們關注的是變量的值的大小，可以使用歐氏距離來聚類
數據標準化當變量是由不同的標度測量的時候，最好要對數據進行標準化使之可以進行比較；一般情況在下對變量進行縮放使之：標準差是1，均值是0；當變量的均值或者標準差相差較大的時候也可以對數據進行scale：
center(x)可以是均值或者中位數；scale(x)可以是標準差,四分位間距,或者絕對中位差(median absolute deviation,MAD)，R裡面可以使用scale()函數進行標準化
❝MAD的定義：數據點到中位數的絕對偏差的中位數
❞計算距離矩陣使用的數據集為USArrests：
df <- USArrestsdf_scaled <- scale(df)##標準化計算距離的R函數有很多，如：
get_dist() factoextra包裡面的,可以計算基於相關性的距離，包括「pearson」, 「kendall」 「spearman」daisy()cluster包裡面的，可以處理除了數值變量以外的其他變量類型(如分類變量,定序變量等)注意：這些計算距離的函數都是計算行之間的距離，所以如果我們要計算變量(列)之間的距離，先要將其轉置
dist.eucl <- dist(df_scaled,method = "euclidean")##歐氏距離class(dist.eucl)#[1] "dist"as.matrix(dist.eucl)[1:3,1:3]# Alabama Alaska Arizona# Alabama 0.000000 2.703754 2.293520# Alaska 2.703754 0.000000 2.700643# Arizona 2.293520 2.700643 0.000000library(factoextra)dist_cor <- get_dist(df_scaled,method = "pearson")##相關性距離as.matrix(dist_cor)[1:3,1:3]# Alabama Alaska Arizona# Alabama 0.0000000 0.7138308 1.4465948# Alaska 0.7138308 0.0000000 0.8307246# Arizona 1.4465948 0.8307246 0.0000000
library(cluster)data("flower")head(flower,3)# V1 V2 V3 V4 V5 V6 V7 V8# 1 0 1 1 4 3 15 25 15# 2 1 0 0 2 1 3 150 50# 3 0 1 0 3 3 1 150 50str(flower)# 'data.frame': 18 obs. of 8 variables:# $ V1: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 2 2 ...# $ V2: Factor w/ 2 levels "0","1": 2 1 2 1 2 2 1 1 2 2 ...# $ V3: Factor w/ 2 levels "0","1": 2 1 1 2 1 1 1 2 1 1 ...# $ V4: Factor w/ 5 levels "1","2","3","4",..: 4 2 3 4 5 4 4 2 3 5 ...# $ V5: Ord.factor w/ 3 levels "1"<"2"<"3": 3 1 3 2 2 3 3 2 1 2 ...# $ V6: Ord.factor w/ 18 levels "1"<"2"<"3"<"4"<..: 15 3 1 16 2 12 13 7 4 14 ...# $ V7: num 25 150 150 125 20 50 40 100 25 100 ...# $ V8: num 15 50 50 50 15 40 20 15 15 60 ...dd <- daisy(flower)##計算含有分類變量，定序變量的距離as.matrix(dd)[1:3,1:3]# 1 2 3# 1 0.0000000 0.8875408 0.5272467# 2 0.8875408 0.0000000 0.5147059# 3 0.5272467 0.5147059 0.0000000劃分聚類(Partitioning clustering)劃分聚類需要我們指定類別的數量
最常用的有：
K-medoids clustering (PAM)K均值聚類k表示我們想要數據聚成的類數，最終的結果是實現高的類內相似性和低的類間相似性
我們的目的就是使上式最小化
算法隨機選取k個點作為起始聚類中心(initial cluster centers)更新聚類中心：計算每個類的數據點的平均值作為新的聚類中心迭代3,4步，直到聚類狀態不再變化或者達到最大的迭代數目(R中默認是10)Rkmeans(x, centers, iter.max = 10, nstart = 1)centers: 類數或者起始的距離中心,如果輸入的是一個數值的話則隨機選取x的行作為初始聚類中心nstart: 開始選擇隨機聚類中心的次數,比如nstart=5,則是開始隨機選擇5次k個聚類中心,最後選擇結果最好的如何選擇最佳聚類數？
一個簡單的方法就是嘗試不同的聚類數目k,計算上面的total within sum of square；隨著聚類數目的增加WSS的趨勢一定會下降(最極端的情況就是每個點都是一個類),當k小於真實聚類數時WSS下降幅度會很大,而當k大於真實聚類數時，再增加k WSS的下降幅度會驟減,所以會存在一個拐點
library(factoextra)fviz_nbclust(df_scaled,kmeans,method = "wss")+ geom_vline(xintercept = 4,linetype=2)image-20200720233852847最佳聚類數為4：
set.seed(2020720)km_res <- kmeans(df_scaled,4,nstart = 25)print(km_res)
names(km_res)# [1] "cluster" "centers" "totss" "withinss" "tot.withinss"# [6] "betweenss" "size" "iter" "ifault"結果包括聚類的中心點，每個觀測值所屬的類，wss
對於高維的數據我們可以先降維再可視化聚類的結果：
fviz_cluster(km_res, data = df_scaled)image-20200720235320015K-Medoids在k-medoids聚類中每個類由類內的某個點來代替，這些點就叫聚類中心(cluster medoids)
在 K-means 算法中，我們每次選簇的平均值作為新的中心，迭代直到簇中對象分布不再變化。因此一個具有很大極端值的對象會扭曲數據分布，造成算法對極端值敏感; K-Medoids算法不選用平均值而是用中心點作為參照點
最常用的k-medoids聚類方法是PAM算法(Partitioning Around Medoids)
PAM 算法隨機選擇k個點作為medoids(或者指定k個點)在每一類裡面，對除初始的medoids點外的所有其他點，按順序計算當其為新的medoids時準則函數的值是否下降,選擇使其下降最多的點作為新的中心點(準則函數為所有點到其最近中心點的距離的和)迭代3，4直到準則函數不再下降(medoids不再變化)Rcluster::pam(x, k, metric = "euclidean", stand = FALSE)x : 可以是數值矩陣或者數據框,行是觀測,列是變量；也可以是距離矩陣metric : 計算距離的方法,可以是euclidean或者manhattan首先需要估計最佳聚類數，可以使用平均輪廓法(average silhouette method),平均輪廓值越高說明聚類質量越好
可以使用factoextra包中的fviz_nbclust函數來計算：
fviz_nbclust(df_scaled,pam,method = "silhouette")+ theme_classic()image-20200721222804298最佳聚類數為2：
pam_res <- pam(df_scaled,2)names(pam_res)# [1] "medoids" "id.med" "clustering" "objective" "isolation" "clusinfo" # [7] "silinfo" "diss" "call" "data"繼續使用fviz_cluster進行PCA降維後可視化聚類的效果：
fviz_cluster(pam_res, data = df_scaled)image-20200721223407152CLARA聚類CLARA (Clustering Large Applications)是k-medoids聚類的延伸，用來處理比較大的數據集
算法應用PAM算法找出每個亞數據集的中心點，分別將每個亞數據集的中心點應用到整個數據集計算所有數據點到最近中心點的距離和，保留最小距離和的亞數據集的中心點重複1，2步如果計算的距離和小於上次最小的距離和則用新的中心點代替原來的中心點直至中心點不再變化R##generate dataset.seed(2020721)df <- rbind(cbind(rnorm(200,0,8), rnorm(200,0,8)), cbind(rnorm(300,50,8), rnorm(300,50,8)))colnames(df) <- c("x","y")rownames(df) <- paste("S",1:nrow(df),sep = "")cluster::clara(x, k, metric = "euclidean", stand = FALSE, samples = 5, pamLike = FALSE)metric : 距離計算方法可以選擇euclidean或者manhattanpamLike：是否和pam()函數使用相同的算法首先使用silhouette方法來估計最佳聚類數：
fviz_nbclust(df, clara, method = "silhouette")+ theme_classic()image-20200721233550297最佳聚類數為2：
clara_res <- clara(df,2,samples = 50,pamLike = T)print(clara_res)names(clara_res)# [1] "sample" "medoids" "i.med" "clustering" "objective" "clusinfo" # [7] "diss" "call" "silinfo" "data"fviz_cluster(clara_res,data = df)image-20200721234056740層次聚類(Hierarchical clustering)層次聚類和劃分聚類一個顯著不同就是層次聚類不需要預先規定聚類數目
凝聚方法（agglomerative hierarchical clustering):自底向上，每個觀察值最初都被視為一類(葉),然後將最相似的類連續合併，直到只有一個大類(根)為止分裂方法（divisive hierarchical clustering):自上向下，是凝聚聚類的逆過程，從根開始，所有觀測值都包含在一個類中然後將最不均一的聚類相繼劃分直到所有觀測值都在它們自己的類中(葉)image-20200722083259840凝聚聚類使用連接函數(linkage function)基於距離信息將對象連接成層次聚類樹連接函數獲取由函數dist()返回的距離信息,並根據對象的相似性將對象對分組;重複此過程,直到原始數據集中的所有對象在層次樹中連結在一起為止
res_hc <- stats::hclust(d = dist.eucl, method = "ward.D2")method: 連接函數,可以是「ward.D」, 「ward.D2」, 「single」,「complete」, 「average」, 「mcquitty」, 「median」 「centroid」主要使用的連接函數(也就是類間距離)有:
最長距離法(complete-linkage):兩個類的距離定義為兩個類的元素的所有成對距離的最大值最短距離法(single-linkage): 兩個類的距離定義為兩個類的元素的所有成對距離的最小值類平均法(mean or average linkage,UPGMA): 兩個類的距離定義為兩個類的元素的所有成對距離的平均值中心法(centroid linkage,UPGMC): 兩個聚類之間的距離定義為兩個類的質心(centroid,該類所有點的均值)的距離Ward法: 最小化總的within-cluster variationplot(res_hc)fviz_dend(res_hc, cex = 0.5)image-20200722133701585連接兩個對象的豎線的高度衡量了這兩個對象的距離,越長距離越大,這個高度也叫這兩個對象的共同距離cophenetic distance
兩個點的共同距離是這兩個點第一次被聚在一起時的節點的高度，截取一小部分說明：
image-20200725182546137我們可以看聚類樹的共同距離和原始的距離矩陣的相似性來衡量聚類的好壞：
res_coph <- cophenetic(res_hc)cor(dist_eucl, res_coph)#[1] 0.6975266
res_hc2 <- hclust(dist_eucl, method = "average") ##類平均法cor(dist_eucl, cophenetic(res_hc2))#[1] 0.7180382可以使用cutree函數來分割樹進行聚類結果的展示：
cluster <- cutree(res_hc, k = 4)table(cluster)# cluster# 1 2 3 4 # 7 12 19 12
factoextra::fviz_dend(res_hc, k = 4,rect = TRUE)##可視化
fviz_cluster(list(data = df_scaled, cluster = cluster))##先PCA再展現聚類結果image-20200722133714891image-20200722133728396另外cluster包也可以很方便的進行凝聚或者分裂聚類:
library("cluster")# Agglomerative res.agnes <- agnes(x = USArrests, # data matrix stand = TRUE, # Standardize the data metric = "euclidean", # metric for distance matrix method = "ward" # Linkage method)# Divisiveres.diana <- diana(x = USArrests, # data matrix stand = TRUE, # standardize the data metric = "euclidean" # metric for distance matrix )比較樹狀圖使用dendextend包
首先創建兩個不同的樹狀圖：
dend1 <- stats::as.dendrogram(res_hc) dend2 <- stats::as.dendrogram(res_hc2)dend_list <- dendextend::dendlist(dend1, dend2)dendextend::tanglegram(dend1, dend2)image-20200722135343171比對的質量可以使用entanglement()函數來計算，這個值是0到1之間的，越小說明比對越好
dendextend::entanglement(dend1, dend2)#[1] 0.8342094也可以使用cor.dendlist()函數來計算兩個樹的相關性,有兩種方法cophenetic和baker
# Cophenetic correlation matrixdendextend::cor.dendlist(dend_list, method = "cophenetic")# [,1] [,2]# [1,] 1.000000 0.843143# [2,] 0.843143 1.000000# Baker correlation matrixdendextend::cor.dendlist(dend_list, method = "baker")# [,1] [,2]# [1,] 1.0000000 0.8400675# [2,] 0.8400675 1.0000000選擇最佳聚類數直接法：最優化一個準則，比如within cluster variation(肘方法,elbow method)或者平均輪廓(輪廓法,silhouette method)將總的WSS看做是聚類數的函數，當增加聚類數不會大幅度降低WSS時會出現一個拐點，選擇該點作為最佳聚類數
Average silhouette method(平均輪廓法)該方法需要計算輪廓係數：
計算對象i到同類其他對象的平均距離
image-20200722142645591「所有樣本的
和肘方法相似，計算不同聚類數目的輪廓係數，輪廓係數最大的聚類數為最佳聚類數
一般選擇B=500，結果就比較穩健
Rfactoextra::fviz_nbclust(x, FUNcluster, method = c("silhouette", "wss", "gap_stat"))FUNcluster: 聚類方法包括kmeans, pam, clara 和 hcut# Elbow methodp1 <- factoextra::fviz_nbclust(df_scaled, kmeans, method = "wss") + geom_vline(xintercept = 4, linetype = 2)+ labs(subtitle = "Elbow method")# Silhouette methodp2 <- factoextra::fviz_nbclust(df_scaled, kmeans, method = "silhouette")+ labs(subtitle = "Silhouette method")# Gap statistic# nboot = 50 to keep the function speedy.# recommended value: nboot= 500 for your analysis.# Use verbose = FALSE to hide computing progression.set.seed(123)p3 <- factoextra::fviz_nbclust(df_scaled, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+ labs(subtitle = "Gap statistic method")cowplot::plot_grid(p1,p2,p3)image-20200722144751617NbClust包的NbClust()函數可以提供30種指標來計算最佳聚類數
NbClust(data = NULL, diss = NULL, distance = "euclidean", min.nc = 2, max.nc = 15, method = NULL)diss：不相似矩陣，如果提供diss，後面的distance應該為NULLdistance: 計算距離矩陣的方法包括euclidean manhattan和NULLmethod: 聚類方法"ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median", "centroid", "kmeans"nb <- NbClust::NbClust(df_scaled, distance = "euclidean", min.nc = 2, max.nc = 10, method = "kmeans")factoextra::fviz_nbclust(nb)image-20200722145610875還有一個增強版的hclust:fastcluster::hclust:更快，能處理更大的數據
image-20200722150942964參考：
https://cran.r-project.org/web/packages/flashClust/index.html
https://www.biostars.org/p/85150/
https://www.cnblogs.com/lexus/archive/2012/12/13/2815769.html
Kassambara A. Practical guide to cluster analysis in R: Unsupervised machine learning[M]. Sthda, 2017.
Reference[1]R 聚類圖書: r-cluster-book.pdf