> library(pacman)
> p_load(dplyr, mclust, GGally)
聚類：在數(shù)據(jù)中識別相似行的技術。
常見聚類技術：k-means, DBSCAN, OPTICS
k-means 是一種基于劃分的聚類算法，它以 k 為參數(shù)，把 n 個數(shù)據(jù)對象分成 k 個簇，使簇內(nèi)具有較高的相似度，而簇間的相似度較低。常見算法有：Lloyd 算法；MacQueen 算法和Hartigan-Wong 算法。
1、k-means聚類思想
k-means 算法是根據(jù)給定的 n 個數(shù)據(jù)對象的數(shù)據(jù)集，構建 k 個劃分聚類的方法，每個劃分聚類即為一個簇。該方法將數(shù)據(jù)劃分為 n 個簇，每個簇至少有一個數(shù)據(jù)對象，每個數(shù)據(jù)對象必須屬于而且只能屬于一個簇。同時要滿足同一簇中的數(shù)據(jù)對象相似度高，不同簇中的數(shù)據(jù)對象相似度較小。聚類相似度是利用各簇中對象的均值來進行計算的。
k-means 算法的處理流程如下：
首先，隨機地選擇 k 個數(shù)據(jù)對象，每個數(shù)據(jù)對象代表一個簇中心，即選擇 k 個初始中心；對剩余的每個對象，根據(jù)其與各簇中心的相似度（距離），將它賦給與其最相似的簇中心對應的簇；
然后重新計算每個簇中所有對象的平均值，作為新的簇中心；
不斷重復以上過程，直到準則函數(shù)（通常采用均方差作為準則函數(shù)，即最小化每個點到最近簇中心的距離的平方和）收斂，也就是簇中心不發(fā)生明顯的變化。
2、算法步驟
2.1 Lloyd算法步驟
1.選擇k個中心。
2.在特征空間中隨機初始化k個中心。
3.對每行計算其和每個中心之間的距離。
4.將各行分配到最近中心點的簇中。
5.將每個中心移動到其所在簇的平均值處。
6.重復步驟3和4，直到?jīng)]有行改變簇或達到最大迭代次數(shù)。
2.2 MacQueen算法步驟
Lloyd的算法:批處理或離線算法。
MacQueen的算法：增量或在線算法。
1.選擇k個中心。
2.在特征空間中隨機初始化k個中心。
3.將各行按距離分配到最近中心點的簇中。
4.將每個中心移動到其所在簇的平均值處。
5.對每行，與各中心計算距離并分配到最近中心點的簇中。若有行換簇，則更新中心。
6.每行處理完后，更新所有中心點。
7.若沒有行改變簇則停止，否則重復步驟5。
2.3 Hartigan-Wong算法步驟
1.選擇k個中心。
2.在特征空間中隨機初始化k個中心。
3.將各行按距離分配到最近中心點的簇中。
4.將每個中心移動到其所在簇的平均值處。
5.對每行，計算其在每簇時的SSE，并分到SSE最小的簇中。若有行換簇，則更新中心。
6.若沒有行改變簇則停止，否則重復步驟5。
2.4 距離計算
R中dist()函數(shù)計算距離。
dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
x：數(shù)據(jù)矩陣或數(shù)據(jù)框
method：距離計算方法，包括"euclidean"-歐氏距離，"maximum"-切比雪夫距離，"manhattan"-曼哈頓距離（絕對值距離，馬氏距離），"canberra"-蘭氏距離，"binary"-定性變量距離，"minkowski"-明考斯基距離。
diag,upper：輸出距離矩陣的式樣，默認只輸出下三角矩陣。
p：明氏距離的冪次。
在計算距離之前，一般需要將數(shù)據(jù)中心化或標準化處理，消除不同量綱或量級的影響。
3 k-means聚類實例
> data(GvHD, package = "mclust")
> gvhd <- as_tibble(GvHD.control)
> str(gvhd)
## tibble [6,809 × 4] (S3: tbl_df/tbl/data.frame)
## $ CD4 : num [1:6809] 199 294 85 19 35 376 97 200 422 391 ...
## $ CD8b: num [1:6809] 420 311 79 1 29 346 329 342 433 390 ...
## $ CD3 : num [1:6809] 132 241 14 141 6 138 527 145 163 147 ...
## $ CD8 : num [1:6809] 226 164 218 130 135 176 406 189 47 190 ...
> DataExplorer::profile_missing(gvhd)
## # A tibble: 4 x 3
## feature num_missing pct_missing
## <fct> <int> <dbl>
## 1 CD4 0 0
## 2 CD8b 0 0
## 3 CD3 0 0
## 4 CD8 0 0
探索性可視化：
> ggpairs(gvhd, upper = list(continuous = "density"),
+ lower = list(continuous = wrap("points", size = 0.4)),
+ diag = list(continuous = "densityDiag"))
從圖中隱約可以看到分成３類較為合適。
k-means聚類：
> # 標準化
> gvhd.scale <- scale(gvhd)
> # 分別使用三種算法聚類
> # centers：聚類的個數(shù)或初識類重心
> # iter.max：最大迭代次數(shù)，默認10
> # nstart：隨機集合個數(shù)，默認1
> # algorithm：動態(tài)聚類算法，默認為Hartigan-Wong
> gvhd.lloyd <- kmeans(gvhd.scale, centers = 3, algorithm = "Lloyd")
> gvhd.hartiganwong <- kmeans(gvhd.scale, centers = 3, algorithm = "Hartigan-Wong")
> gvhd.macqueen <- kmeans(gvhd.scale, centers = 3, algorithm = "MacQueen")
結果解讀：
> gvhd.lloyd
## K-means clustering with 3 clusters of sizes 4818, 433, 1558
##
## Cluster means:
## CD4 CD8b CD3 CD8
## 1 0.4891758 0.3349026 -0.02166058 -0.24424752
## 2 -0.9043872 0.5886317 1.68486535 2.77669629
## 3 -1.2613925 -1.1992544 -0.40127474 -0.01638313
##
## Clustering vector:
## [1] 1 1 3 3 3 1 2 1 1 1 1 1 1 3 3 2 3 3 3 3 1 3 1 1 1 1 3 1 1 1 1 1 1 1 1 3 1 2 3 1 1 2 1 3 1 3 1 1 2 1 1 3
## [53] 1 1 1 1 3 1 1 1 3 1 3 1 3 1 1 1 1 1 1 1 2 1 3 1 1 1 1 1 1 1 3 3 2 1 1 3 1 3 1 1 1 2 1 1 1 1 1 3 1 1 1 3
## [105] 3 1 1 1 1 1 1 1 1 1 1 3 1 1 1 3 3 1 1 3 1 1 3 1 1 1 1 1 1 1 1 1 3 1 3 1 1 1 3 1 1 3 1 3 1 1 1 1 1 3 3 3
## [157] 1 3 1 1 3 2 1 1 2 1 1 1 3 1 1 1 1 1 1 1 3 1 3 3 1 1 3 1 1 3 1 1 1 3 3 3 1 3 1 1 2 1 1 1 1 1 1 1 3 3 3 1
## [209] 1 3 1 1 3 3 3 1 1 3 3 1 1 3 1 1 3 1 1 3 1 1 3 1 1 1 3 3 1 1 1 3 1 1 1 1 2 1 1 1 1 1 3 1 3 3 3 3 3 1 2 3
## [261] 1 3 3 1 2 3 1 3 1 1 1 3 1 3 1 1 1 3 3 1 1 1 1 3 3 1 1 1 3 2 3 1 1 1 3 1 1 3 1 3 1 2 1 3 3 1 1 1 3 3 2 1
## [313] 1 3 3 3 1 1 1 1 1 3 1 1 1 1 1 1 3 1 1 1 1 1 1 1 3 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1
## [365] 3 1 1 1 1 2 1 3 1 1 1 1 3 3 1 3 3 1 3 1 1 1 1 1 3 3 3 1 1 1 1 1 1 3 1 1 2 3 1 3 1 1 3 3 1 1 1 3 1 1 1 1
## [417] 1 1 1 1 3 1 1 1 1 1 1 1 3 1 3 1 1 1 1 3 1 1 1 1 1 1 1 2 3 1 2 2 1 2 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 3 1
## [469] 1 1 3 1 3 3 1 1 1 1 1 1 1 1 1 3 1 3 1 1 3 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 3 3 1 1 3 1 3 1 1 1 1 1 2 1 1 1
## [521] 1 1 1 1 1 3 3 3 2 1 1 1 1 1 1 1 1 2 1 1 3 3 3 1 3 3 1 1 1 1 1 1 3 1 1 1 1 3 1 1 2 1 1 3 3 1 3 1 2 3 1 1
## [573] 2 1 3 2 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 1 1 1 3 2 3 1 1 1 1 1 2 1 1
## [625] 1 1 1 3 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 3 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 3 3 3 1 3 1 1 1 3 1 2
## [677] 1 3 1 3 1 2 3 3 3 1 2 3 1 1 1 1 3 3 1 3 1 3 1 1 3 1 3 2 1 1 3 1 2 1 1 1 1 1 3 2 1 3 1 1 1 1 1 1 1 3 1 1
## [729] 1 1 1 1 3 1 2 1 1 1 3 3 1 1 3 3 1 3 1 1 1 1 1 1 2 3 1 1 1 1 1 1 1 1 3 3 2 1 3 1 1 1 1 1 3 1 3 1 1 1 1 1
## [781] 1 2 1 3 1 1 1 1 1 1 1 3 1 1 3 1 1 1 3 3 1 1 3 1 1 1 1 1 1 1 1 1 1 3 3 1 1 3 2 3 1 1 3 1 1 1 3 3 1 1 1 1
## [833] 1 1 1 1 3 1 3 1 3 2 1 1 1 3 3 1 1 1 1 3 1 3 1 3 1 1 3 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 3 3 1
## [885] 3 2 3 1 1 1 1 3 3 1 1 1 3 3 3 1 1 1 3 1 1 3 3 3 3 1 2 1 1 1 3 1 1 2 1 1 1 3 1 1 3 1 1 3 1 2 3 2 1 1 1 2
## [937] 1 1 1 1 1 1 2 1 3 1 1 1 1 1 1 1 1 1 3 1 1 3 1 1 1 1 2 1 1 1 1 1 1 1 1 3 1 3 1 1 1 1 1 1 1 2 1 1 1 1 1 1
## [989] 1 1 3 1 2 1 1 1 1 1 1 1
## [ reached getOption("max.print") -- omitted 5809 entries ]
##
## Within cluster sum of squares by cluster:
## [1] 11329.51 1451.36 2425.35
## (between_SS / total_SS = 44.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss" "betweenss" "size"
## [8] "iter" "ifault"
K-means clustering with 3 clusters of sizes 4293, 1847, 669：３個類的大小分別為4293, 1847, 669。
Cluster means：中心點的值。
Clustering vector：分組索引。
Within cluster sum of squares by cluster：組內(nèi)平方和。
(between_SS / total_SS = 50.0 %)：組間平方和/總平方和，用于衡量點聚集程度。
Available components：聚類后對象(gvhd.lloyd)的屬性。比如要查看每個類的大小：
> gvhd.lloyd$size
## [1] 4818 433 1558
對象屬性解讀：
cluster，每個點的分組
centers，聚類的中心點坐標
totss，總平方和
withinss，組內(nèi)平方和
tot.withinss，組內(nèi)總平方和，sum(withinss)
betweenss，組間平方和，totss – tot.withinss
size，每個組中的數(shù)據(jù)點數(shù)量
iter，迭代次數(shù)
ifault，可能有問題的指標
對比三種算法的聚類結果，分別對比組內(nèi)平方總和和組間平方和。
> tibble(algorithm = c("Lloyd", "Hartigan-Wong","MacQueen"),
+ withinss = c(gvhd.lloyd$totss, gvhd.hartiganwong$totss,
+ gvhd.macqueen$totss),
+ betweenss = c(gvhd.lloyd$betweenss, gvhd.hartiganwong$betweenss,
+ gvhd.macqueen$betweenss))
## # A tibble: 3 x 3
## algorithm withinss betweenss
## <chr> <dbl> <dbl>
## 1 Lloyd 27232. 12026.
## 2 Hartigan-Wong 27232. 13396.
## 3 MacQueen 27232. 12258.
要讓組內(nèi)平方和盡量小，組間平方和盡量大，所以"Lloyd"和"MacQueen"算法結果相差很小，都比"Hartigan-Wong"算法更差。
4 k-means聚類的優(yōu)缺點
優(yōu)點：原理簡單，收斂快，效果優(yōu)，參數(shù)少（只有Ｋ）。
缺點：Ｋ不好把握；不能自然地處理分類變量；對不同尺度的數(shù)據(jù)敏感；要求隱含類方差近似相等；對異常點敏感；采用迭代只是局部最優(yōu)；由于初始中心點的隨機性，每次聚類可能結果略有區(qū)別；優(yōu)先尋找球形聚類，即使不符實情。
5 k-means聚類mlr3實現(xiàn)
真正最新理念、最新技術、最新一代的系統(tǒng)地實現(xiàn)機器學習算法的R包，比Python中的 sklearn 還先進。
任務（Task）：封裝了數(shù)據(jù)及額外信息，如預測目標；
學習器（Learner）：封裝了多種機器學習算法，能夠訓練模型并做預測，大多數(shù)學習器都有影響其操作的超參數(shù)；
度量（Measure）：基于預測值與真實值的差異計算數(shù)值得分；
重抽樣（Resampling）：生成一系列的訓練集和測試集。
另外，實際數(shù)據(jù)的機器學習建模工作流，還包括標準化、缺失值插補、特征選擇等。
初識mlr3：http://www.itdecent.cn/p/7766aa7f2ef7
5.1 創(chuàng)建任務
> p_load(mlr3, mlr3cluster)
>
> # 創(chuàng)建任務
> task <- TaskClust$new(id = "gvhd", backend = gvhd)
>
> # 探索性可視化
> mlr3viz::autoplot(task, type = "pairs")
其實也是使用GGally包畫的圖。
5.2 創(chuàng)建學習器
> # 查看支持哪些學習器
> mlr_learners$keys("clust")
## [1] "clust.agnes" "clust.cmeans" "clust.dbscan" "clust.diana" "clust.fanny"
## [6] "clust.featureless" "clust.kmeans" "clust.pam" "clust.xmeans"
可以看到，支持的聚類算法有：clust.agnes, clust.cmeans, clust.dbscan, clust.diana, clust.fanny, clust.featureless, clust.kmeans, clust.pam, clust.xmeans。
我們使用k-means聚類：
> learner.kmeans <- mlr_learners$get("clust.kmeans")
> # 查看信息
> print(learner.kmeans)
## <LearnerClustKMeans:clust.kmeans>
## * Model: -
## * Parameters: centers=2
## * Packages: stats, clue
## * Predict Type: partition
## * Feature types: logical, integer, numeric
## * Properties: complete, exclusive, partitional
> # 查看可以設置哪些超參數(shù)
> learner.kmeans$param_set$ids()
## [1] "centers" "iter.max" "algorithm" "nstart" "trace"
跟kmeans函數(shù)的參數(shù)是一致的。
設置十折交叉驗證：
> resampling <- rsmp("cv", folds = 10L)
5.3 訓練和度量
度量函數(shù)：
> # 查看有哪些度量函數(shù)
> mlr_measures$keys("clust")
## [1] "clust.ch" "clust.db" "clust.dunn" "clust.silhouette"
clust.ch：Calinski-Harabasz 偽 F 統(tǒng)計量，反映聚類間方差和聚類內(nèi)方差的比率，CH越大代表著類自身越緊密，類與類之間越分散，即更優(yōu)的聚類結果。
clust.db：Davies-Bouldin Index(戴維森堡丁指數(shù))(分類適確性指標)(DB，DBI)，DB越小意味著類內(nèi)距離越小，同時類間距離越大。缺點：因使用歐式距離，所以對于環(huán)狀分布聚類評測很差。
clust.dunn：Dunn Validity Index (鄧恩指數(shù))(DVI)，DVI越大意味著類間距離越大，同時類內(nèi)距離越小。缺點：對離散點的聚類測評很高、對環(huán)狀分布測評效果差。
clust.silhouette：輪廓系數(shù)，是聚類效果好壞的一種評價方式。最早由 Peter J. Rousseeuw 在 1986 提出。它結合內(nèi)聚度和分離度兩種因素?？梢杂脕碓谙嗤紨?shù)據(jù)的基礎上評價不同算法、或者算法不同運行方式對聚類結果所產(chǎn)生的影響。如果一個簇中的大多數(shù)樣本具有比較高的輪廓系數(shù)（接近1），則簇會有較高的總輪廓系數(shù)，則整個數(shù)據(jù)集的平均輪廓系數(shù)越高，則聚類是合適的。如果許多樣本點具有低輪廓系數(shù)甚至負值（接近-1），則聚類是不合適的，聚類的超參數(shù)K可能設定得太大或者太小。
> measure <- msrs(c("clust.ch", "clust.db", "clust.silhouette"))
根據(jù)不同的算法，設計訓練：
> design <- benchmark_grid(
+ task = task,
+ learner = list(
+ lrn("clust.kmeans", algorithm = "Hartigan-Wong"),
+ lrn("clust.kmeans", algorithm = "Lloyd"),
+ lrn("clust.kmeans", algorithm = "MacQueen")
+ ),
+ resampling = resampling)
> print(design)
## task learner resampling
## 1: <TaskClust[41]> <LearnerClustKMeans[31]> <ResamplingCV[19]>
## 2: <TaskClust[41]> <LearnerClustKMeans[31]> <ResamplingCV[19]>
## 3: <TaskClust[41]> <LearnerClustKMeans[31]> <ResamplingCV[19]>
> bmr <- benchmark(design = design)
> results <- bmr$aggregate(measures = measure)
結果表明Hartigan-Wong算法在該數(shù)據(jù)集上是最優(yōu)的。
畫個圖看看：
> tibble(algorithm = c("Hartigan-Wong", "Lloyd", "MacQueen"),
+ # 數(shù)值同比例縮小，便于畫圖
+ clust.ch = results$clust.ch / 200,
+ clust.db = results$clust.db,
+ clust.silhouette = results$clust.silhouette) %>%
+ # 寬表變長表
+ tidyr::pivot_longer(names_to = "metric", values_to = "value", -c("algorithm")) %>%
+ ggplot(aes(metric, value, fill = algorithm)) +
+ geom_bar(stat = "identity", position = "dodge", width = 0.5) +
+ labs(x = "", y = "", title = "不同聚類算法對比") +
+ theme_bw() +
+ theme(legend.position = "top",
+ plot.title = element_text(hjust = 0.5))
使用同樣的方法可以找到最優(yōu)的Ｋ值，然后使用最優(yōu)的參數(shù)重新建模，并對新數(shù)據(jù)進行預測。
