library(pacman)
p_load(dbscan, ggplot2, dplyr)

DBSCAN在输入参数的选取上比较困难，即DBSCAN对输入参数比较敏感。当给定全局参数eps和minPts时，会存在问题：不同的全局参数会得到不同的聚类结果。

OPTICS(Ordering Points To Identify the Clustering Structure, OPTICS)实际上是DBSCAN算法的一种有效扩展，主要解决对输入参数敏感的问题。它使得基于密度的聚类结构能够呈现出一种特殊的顺序，该顺序所对应的聚类结构包含了每个层级的聚类的信息，并且便于分析。
OPTICS并不显示产生的结果类簇，而是为聚类分析生成一个排序，这个排序代表了各样本点基于密度的聚类结构。换句话说，从这个排序中可以得到基于任何参数eps和minPts的DBSCAN算法的聚类结果。

1 基本概念

核心距离：假定P是核心对象，人为给定一个阈值Ａ，然后计算关于Ｐ点满足阈值Ａ的最小的半径Ｒ，即在Ｒ内，Ｐ最少有给定Ａ个点数。
可达距离：对象q到对象p的可达距离是指p的核心距离和p与q之间欧几里得距离之间的较大值。如果p不是核心对象，p和q之间的可达距离没有意义。

2 聚类实例

> data(banknote, package = "mclust")
> bn <- as_tibble(banknote)
> str(bn)

## tibble [200 × 7] (S3: tbl_df/tbl/data.frame)
##  $ Status  : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Length  : num [1:200] 215 215 215 215 215 ...
##  $ Left    : num [1:200] 131 130 130 130 130 ...
##  $ Right   : num [1:200] 131 130 130 130 130 ...
##  $ Bottom  : num [1:200] 9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
##  $ Top     : num [1:200] 9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 11 10 ...
##  $ Diagonal: num [1:200] 141 142 142 142 142 ...

> DataExplorer::profile_missing(bn)

## # A tibble: 7 x 3
##   feature  num_missing pct_missing
##   <fct>          <int>       <dbl>
## 1 Status             0           0
## 2 Length             0           0
## 3 Left               0           0
## 4 Right              0           0
## 5 Bottom             0           0
## 6 Top                0           0
## 7 Diagonal           0           0

> # 因为做聚类，所以去掉类别列
> bn <- bn[, -1] %>% 
+   # 标准化
+   mutate(across(everything(), scale))
> str(bn)

## tibble [200 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Length  : num [1:200, 1] -0.255 -0.786 -0.255 -0.255 0.276 ...
##   ..- attr(*, "scaled:center")= num 215
##   ..- attr(*, "scaled:scale")= num 0.377
##  $ Left    : num [1:200, 1] 2.43 -1.17 -1.17 -1.17 -1.44 ...
##   ..- attr(*, "scaled:center")= num 130
##   ..- attr(*, "scaled:scale")= num 0.361
##  $ Right   : num [1:200, 1] 2.83 -0.635 -0.635 -0.882 -0.635 ...
##   ..- attr(*, "scaled:center")= num 130
##   ..- attr(*, "scaled:scale")= num 0.404
##  $ Bottom  : num [1:200, 1] -0.289 -0.912 -0.497 -1.327 0.68 ...
##   ..- attr(*, "scaled:center")= num 9.42
##   ..- attr(*, "scaled:scale")= num 1.44
##  $ Top     : num [1:200, 1] -1.184 -1.433 -1.308 -0.312 -3.675 ...
##   ..- attr(*, "scaled:center")= num 10.7
##   ..- attr(*, "scaled:scale")= num 0.803
##  $ Diagonal: num [1:200, 1] 0.448 1.056 1.49 1.316 1.143 ...
##   ..- attr(*, "scaled:center")= num 140
##   ..- attr(*, "scaled:scale")= num 1.15

> bn.clust <- optics(bn, minPts = 9)
> plot(bn.clust)

OPTICS聚类

> # 使用xi方法提取不同密度的层次聚类
> bn.xi <- extractXi(bn.clust, xi = 0.05)
> str(bn.xi)

## List of 10
##  $ order      : int [1:200] 1 199 198 193 144 147 133 124 119 112 ...
##  $ reachdist  : num [1:200] Inf 1.117 1.121 0.886 1.935 ...
##  $ coredist   : num [1:200] 3.01 1.2 1.2 1.09 2.69 ...
##  $ predecessor: int [1:200] NA 91 95 90 16 34 69 78 14 34 ...
##  $ minPts     : num 9
##  $ eps        : num 3.49
##  $ eps_cl     : logi NA
##  $ xi         : num 0.05
##  $ clusters_xi:Classes 'xics' and 'data.frame':  3 obs. of  3 variables:
##   ..$ start     : num [1:3] 1 3 94
##   ..$ end       : num [1:3] 199 90 185
##   ..$ cluster_id: int [1:3] 1 2 3
##  $ cluster    : num [1:200] 1 3 3 3 1 1 3 3 3 3 ...
##  - attr(*, "class")= chr "optics"

画图：

> plot(bn.xi)

xi=0.05

当密度为0.05时，数据集被聚为了3类。
当密度为0.06时，数据集被聚为了2类：

> # 使用xi方法提取不同密度的层次聚类
> bn.xi <- extractXi(bn.clust, xi = 0.06)
> plot(bn.xi)

xi=0.06

> bn %>% 
+   # 添加类别列
+   mutate(cluster = bn.xi$cluster) %>% 
+   # 转换为因子
+   mutate(cluster = as.factor(cluster)) %>% 
+   # 去掉噪声点
+   filter(cluster != 0) %>% 
+   # 画图
+   pairs(col = .$cluster)

聚类后变量相关图