> library(pacman)
> p_load(dplyr, mclust, GGally)

聚类：在数据中识别相似行的技术。
常见聚类技术：k-means, DBSCAN, OPTICS

k-means 是一种基于划分的聚类算法，它以 k 为参数，把 n 个数据对象分成 k 个簇，使簇内具有较高的相似度，而簇间的相似度较低。常见算法有：Lloyd 算法；MacQueen 算法和Hartigan-Wong 算法。

1、k-means聚类思想

k-means 算法是根据给定的 n 个数据对象的数据集，构建 k 个划分聚类的方法，每个划分聚类即为一个簇。该方法将数据划分为 n 个簇，每个簇至少有一个数据对象，每个数据对象必须属于而且只能属于一个簇。同时要满足同一簇中的数据对象相似度高，不同簇中的数据对象相似度较小。聚类相似度是利用各簇中对象的均值来进行计算的。

k-means 算法的处理流程如下：
首先，随机地选择 k 个数据对象，每个数据对象代表一个簇中心，即选择 k 个初始中心；对剩余的每个对象，根据其与各簇中心的相似度（距离），将它赋给与其最相似的簇中心对应的簇；
然后重新计算每个簇中所有对象的平均值，作为新的簇中心；
不断重复以上过程，直到准则函数（通常采用均方差作为准则函数，即最小化每个点到最近簇中心的距离的平方和）收敛，也就是簇中心不发生明显的变化。

2、算法步骤

2.1 Lloyd算法步骤

1.选择k个中心。
2.在特征空间中随机初始化k个中心。
3.对每行计算其和每个中心之间的距离。
4.将各行分配到最近中心点的簇中。
5.将每个中心移动到其所在簇的平均值处。
6.重复步骤3和4，直到没有行改变簇或达到最大迭代次数。

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()函数计算距离。
dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)
x：数据矩阵或数据框
method：距离计算方法，包括"euclidean"-欧氏距离，"maximum"-切比雪夫距离，"manhattan"-曼哈顿距离（绝对值距离，马氏距离），"canberra"-兰氏距离，"binary"-定性变量距离，"minkowski"-明考斯基距离。
diag,upper：输出距离矩阵的式样，默认只输出下三角矩阵。
p：明氏距离的幂次。

在计算距离之前，一般需要将数据中心化或标准化处理，消除不同量纲或量级的影响。

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：聚类的个数或初识类重心
> # iter.max：最大迭代次数，默认10
> # nstart：随机集合个数，默认1
> # algorithm：动态聚类算法，默认为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：组内平方和。
(between_SS / total_SS = 50.0 %)：组间平方和/总平方和，用于衡量点聚集程度。
Available components：聚类后对象(gvhd.lloyd)的属性。比如要查看每个类的大小：

> gvhd.lloyd$size

## [1] 4818  433 1558

对象属性解读：
cluster，每个点的分组
centers，聚类的中心点坐标
totss，总平方和
withinss，组内平方和
tot.withinss，组内总平方和，sum(withinss)
betweenss，组间平方和，totss – tot.withinss
size，每个组中的数据点数量
iter，迭代次数
ifault，可能有问题的指标

对比三种算法的聚类结果，分别对比组内平方总和和组间平方和。

> 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.

要让组内平方和尽量小，组间平方和尽量大，所以"Lloyd"和"MacQueen"算法结果相差很小，都比"Hartigan-Wong"算法更差。

4 k-means聚类的优缺点

优点：原理简单，收敛快，效果优，参数少（只有Ｋ）。
缺点：Ｋ不好把握；不能自然地处理分类变量；对不同尺度的数据敏感；要求隐含类方差近似相等；对异常点敏感；采用迭代只是局部最优；由于初始中心点的随机性，每次聚类可能结果略有区别；优先寻找球形聚类，即使不符实情。

5 k-means聚类mlr3实现

真正最新理念、最新技术、最新一代的系统地实现机器学习算法的R包，比Python中的 sklearn 还先进。
任务（Task）：封装了数据及额外信息，如预测目标；
学习器（Learner）：封装了多种机器学习算法，能够训练模型并做预测，大多数学习器都有影响其操作的超参数；
度量（Measure）：基于预测值与真实值的差异计算数值得分；
重抽样（Resampling）：生成一系列的训练集和测试集。
另外，实际数据的机器学习建模工作流，还包括标准化、缺失值插补、特征选择等。

初识mlr3：https://www.jianshu.com/p/7766aa7f2ef7

5.1 创建任务

> p_load(mlr3, mlr3cluster)
> 
> # 创建任务
> task <- TaskClust$new(id = "gvhd", backend = gvhd)
> 
> # 探索性可视化
> mlr3viz::autoplot(task, type = "pairs")

探索性可视化

其实也是使用GGally包画的图。

5.2 创建学习器

> # 查看支持哪些学习器
> 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

> # 查看可以设置哪些超参数
> learner.kmeans$param_set$ids()

## [1] "centers"   "iter.max"  "algorithm" "nstart"    "trace"

跟kmeans函数的参数是一致的。
设置十折交叉验证：

> resampling <- rsmp("cv", folds = 10L)

5.3 训练和度量

度量函数：

> # 查看有哪些度量函数
> mlr_measures$keys("clust")

## [1] "clust.ch"         "clust.db"         "clust.dunn"       "clust.silhouette"

clust.ch：Calinski-Harabasz 伪 F 统计量，反映聚类间方差和聚类内方差的比率，CH越大代表着类自身越紧密，类与类之间越分散，即更优的聚类结果。
clust.db：Davies-Bouldin Index(戴维森堡丁指数)(分类适确性指标)(DB，DBI)，DB越小意味着类内距离越小，同时类间距离越大。缺点：因使用欧式距离，所以对于环状分布聚类评测很差。
clust.dunn：Dunn Validity Index (邓恩指数)(DVI)，DVI越大意味着类间距离越大，同时类内距离越小。缺点：对离散点的聚类测评很高、对环状分布测评效果差。
clust.silhouette：轮廓系数，是聚类效果好坏的一种评价方式。最早由 Peter J. Rousseeuw 在 1986 提出。它结合内聚度和分离度两种因素。可以用来在相同原始数据的基础上评价不同算法、或者算法不同运行方式对聚类结果所产生的影响。如果一个簇中的大多数样本具有比较高的轮廓系数（接近1），则簇会有较高的总轮廓系数，则整个数据集的平均轮廓系数越高，则聚类是合适的。如果许多样本点具有低轮廓系数甚至负值（接近-1），则聚类是不合适的，聚类的超参数K可能设定得太大或者太小。

> measure <- msrs(c("clust.ch", "clust.db", "clust.silhouette"))

根据不同的算法，设计训练：

> 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算法在该数据集上是最优的。
画个图看看：

> tibble(algorithm = c("Hartigan-Wong", "Lloyd", "MacQueen"),
+        # 数值同比例缩小，便于画图
+        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))

不同算法聚类结果对比

使用同样的方法可以找到最优的Ｋ值，然后使用最优的参数重新建模，并对新数据进行预测。