> 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))
使用同样的方法可以找到最优的Ｋ值，然后使用最优的参数重新建模，并对新数据进行预测。
