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

SOM（Self Organizing Maps，自组织映射）本质上是一种只有输入层--隐藏层的神经网络。输入层神经元的数量是由输入数据的维度决定的，一个神经元对应一个特征，隐藏层中的一个节点代表一个需要聚成的类。训练时采用“竞争学习”的方式，每个输入的样例在隐藏层中找到一个和它最匹配的节点，称为它的激活节点，也叫“winning neuron”。紧接着用随机梯度下降法更新激活节点的参数。同时，和激活节点临近的点也根据它们距离激活节点的远近而适当地更新参数。两种常见邻域函数：bubble function和Gaussian function。

1、SOM步骤总结

1.创建节点网格。
2.随机给每个节点分配权重（数据集中每个变量一个权重（很小的随机数））。
3.随机选择一行，并计算其与网格中每个节点权重的距离（相似度，通常为欧式距离）。
4.把此行放到权重与该行距离最小的节点中（BMU，best matching unit）。
5.更新BMU（基本思想是：越靠近优胜节点，更新幅度越大；越远离优胜节点，更新幅度越小）及其邻域内节点的权重（取决于邻域函数）。
6.重复步骤3-5，迭代指定次数。

2、kohonen包最重要的四个函数

som()、xyf()、supersom()、somgrid()

简单说，som()和xyf()是supersom()的封装版本，分别对应单层SOM和双层SOM，如果是两层以上的多层SOM，必须使用supersom()。somgrid()函数用于建立SOM网络。

3、实例

> data(flea, package = "GGally")
> str(flea)

## 'data.frame':    74 obs. of  7 variables:
##  $ species: Factor w/ 3 levels "Concinna","Heikert.",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ tars1  : int  191 185 200 173 171 160 188 186 174 163 ...
##  $ tars2  : int  131 134 137 127 118 118 134 129 131 115 ...
##  $ head   : int  53 50 52 50 49 47 54 51 52 47 ...
##  $ aede1  : int  150 147 144 144 153 140 151 143 144 142 ...
##  $ aede2  : int  15 13 14 16 13 15 14 14 14 15 ...
##  $ aede3  : int  104 105 102 97 106 99 98 110 116 95 ...

74行7列，1列为因子型，其他全为整数型。

> DataExplorer::profile_missing(flea)

##   feature num_missing pct_missing
## 1 species           0           0
## 2   tars1           0           0
## 3   tars2           0           0
## 4    head           0           0
## 5   aede1           0           0
## 6   aede2           0           0
## 7   aede3           0           0

无缺失值。

> ggpairs(flea, axisLabels = "none", aes(col = species), 
+         upper = list(continuous = ggally_density, 
+                      combo = ggally_box_no_facet)) +
+   theme_bw()

变量相关性

> # xdim/ydim：网格尺寸
> # topo：六边形或矩形，Hexagonal或Rectangular
> # neighbourhood.fct：邻近函数，bubble或gaussian
> # toroidal：是否为环形
> som.grid <- somgrid(xdim = 5, ydim = 5, topo = "hexagonal",
+                     neighbourhood.fct = "bubble", toroidal = F)
> som.grid

## $pts
##         x         y
##  [1,] 1.5 0.8660254
##  [2,] 2.5 0.8660254
##  [3,] 3.5 0.8660254
##  [4,] 4.5 0.8660254
##  [5,] 5.5 0.8660254
##  [6,] 1.0 1.7320508
##  [7,] 2.0 1.7320508
##  [8,] 3.0 1.7320508
##  [9,] 4.0 1.7320508
## [10,] 5.0 1.7320508
## [11,] 1.5 2.5980762
## [12,] 2.5 2.5980762
## [13,] 3.5 2.5980762
## [14,] 4.5 2.5980762
## [15,] 5.5 2.5980762
## [16,] 1.0 3.4641016
## [17,] 2.0 3.4641016
## [18,] 3.0 3.4641016
## [19,] 4.0 3.4641016
## [20,] 5.0 3.4641016
## [21,] 1.5 4.3301270
## [22,] 2.5 4.3301270
## [23,] 3.5 4.3301270
## [24,] 4.5 4.3301270
## [25,] 5.5 4.3301270
## 
## $xdim
## [1] 5
## 
## $ydim
## [1] 5
## 
## $topo
## [1] "hexagonal"
## 
## $neighbourhood.fct
## [1] bubble
## Levels: bubble gaussian
## 
## $toroidal
## [1] FALSE
## 
## attr(,"class")
## [1] "somgrid"

标准化，有助于使每个特征对于计算相似度（距离）的贡献相同。

> flea.scale <- flea %>% 
+   as_tibble() %>%
+   # 去掉species变量
+   select(-species) %>% 
+   # 标准化
+   scale()
> str(flea.scale)

##  num [1:74, 1:6] 0.467 0.263 0.773 -0.145 -0.213 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:6] "tars1" "tars2" "head" "aede1" ...
##  - attr(*, "scaled:center")= Named num [1:6] 177.3 124 50.4 134.8 13 ...
##   ..- attr(*, "names")= chr [1:6] "tars1" "tars2" "head" "aede1" ...
##  - attr(*, "scaled:scale")= Named num [1:6] 29.41 8.48 2.75 10.35 2.14 ...
##   ..- attr(*, "names")= chr [1:6] "tars1" "tars2" "head" "aede1" ...

> # 标准化的两个属性
> # 新数据也需要使用这个属性，否则结果将不一致
> attr(flea.scale, "scaled:center")

##     tars1     tars2      head     aede1     aede2     aede3 
## 177.25676 123.95946  50.35135 134.81081  12.98649  95.37838

> attr(flea.scale, "scaled:scale")

##     tars1     tars2      head     aede1     aede2     aede3 
## 29.412541  8.481146  2.751998 10.350932  2.142162 14.304614

> # rlen：迭代次数
> # alpha：学习速率，默认从0.05下降到0.01
> flea.som <- som(flea.scale, grid = som.grid, rlen = 5000, alpha = c(0.05, 0.01))
> flea.som

## SOM of size 5x5 with a hexagonal topology.
## Training data included.

画图：

> par(mfrow = c(2, 3))
> plot.type <- c("codes", "changes", "counts", "quality", "dist.neighbours", "mapping")
> # 根据plot.type中的类型，依次画图，并按2×3排列
> purrr::walk(plot.type, ~ plot(flea.som, type = ., shape = "straight"))

依次画图

"changes" - Training progress：展示训练过程，距离随着迭代减少的趋势，判断迭代是否足够，最后趋于平稳比较好。
"codes" - Codes plot：查看SOM中心点的变化趋势。
"counts" - Counts plot：展示每个SOM中心点包含的样本数目。可以跟“mapping”一起看，“counts”颜色越浅，对应的“mapping”数量越多。
"dist.neighbours" - Neighbours distance plot：邻近距离，查看潜在边界点，颜色越深表示与周边点差别越大，越可能是边界点。
"mapping" - Mapping plot：展示每个样本的映射。
"quality" - Quality plot：计量SOM中心点的内敛性和质量，距离越小展示得越好。

"property"：每个单元的属性可以计算并显示在颜色代码中。用来可视化一个特定对象与映射中所有单元的相似性，以显示所有单元和映射到它们的对象的平均相似性。

> getCodes(flea.som) %>% 
+   as_tibble() %>% 
+   # property：属性值
+   # main：标题
+   purrr::iwalk(~ plot(flea.som, type = "property", property = .,
+                main = .y, shape = "straight"))

tars1

针对每个特征会有一幅图，就不一一展示了。

4、SOM中心点相关的样本数量

> table(flea.som$unit.classif)

## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 23 24 25 
##  3  3  4  2  2  3  3  7  1  3  3  2  2  2  4  3  3  2  4  5  4  2  4  3

因为定义的网格为5×5，所以一共25个中心点。

> # 不能用flea.scale，因为tibble没有行名
> code.class <- tibble(name = rownames(flea),
+                      class = flea.som$unit.classif)
> head(code.class)

## # A tibble: 6 x 2
##   name  class
##   <chr> <dbl>
## 1 1         3
## 2 2         8
## 3 3         2
## 4 4        13
## 5 5        12
## 6 6        18

5、SOM结果进一步聚类

> # 转换为数据框再转换为矩阵
> mydata <- as.matrix(as.data.frame(flea.som$codes))
> wss <- (nrow(mydata) - 1) * sum(apply(mydata, 2, var)) 
> for (i in 2:15) wss[i] <- sum(kmeans(mydata, centers = i)$withinss)
> 
> plot(1:15, wss, type = "b", xlab = "聚类数量",
+      ylab = "类内平方和", main = "类内平方和 (WCSS)")

选择聚类数量

选择曲线逐渐开始平缓的第一个拐点，本例选择3类。

> # hclust聚类后剪枝为3类
> som.cluster <- cutree(hclust(dist(mydata)), 3)
> 
> # 定义色块颜色
> cluster.palette <- function(x, alpha = 0.6) {
+   n = length(unique(x)) * 2
+   rainbow(n, start = 1/3, end = 3/3, alpha = alpha)[seq(n, 0, -2)]
+ }
> cluster.palette.init <- cluster.palette(som.cluster)
> bgcol <- cluster.palette.init[som.cluster]
> 
> plot(flea.som, type="codes", bgcol = bgcol, main = "Clusters", codeRendering="lines")
> add.cluster.boundaries(flea.som, som.cluster)

重新聚类

查看数据所在新聚的类：

> code.class <- code.class %>% 
+   bind_cols(new_class = som.cluster[code.class$class])
> head(code.class)

## # A tibble: 6 x 3
##   name  class new_class
##   <chr> <dbl>     <int>
## 1 1         3         1
## 2 2         8         1
## 3 3         2         1
## 4 4        13         1
## 5 5        12         2
## 6 6        18         2

6、新数据上应用SOM

> new.data <- tibble(tars1 = c(120, 200),
+                    tars2 = c(125, 120),
+                    head = c(52, 48),
+                    aede1 = c(140, 128),
+                    aede2 = c(12, 14),
+                    aede3 = c(100, 85))
> 
> new.flea <- new.data %>% 
+   # 使用之前的属性值标准化
+   scale(center = attr(flea.scale, "scaled:center"),
+         scale = attr(flea.scale, "scaled:scale")) %>% 
+   # 预测
+   predict(flea.som, newdata = .)
> 
> plot(flea.som, type = "counts", classif = new.flea, shape = "round")

新数据映射

新数据所在新聚的类：

> # 预测的聚类
> new.flea$unit.classif

## [1] 11 14

> # 新聚类
> som.cluster[new.flea$unit.classif]

## V11 V14 
##   2   3

7、与K-Means的比较

（1）K-Means需要事先定下类的个数，也就是K的值。SOM则不用，隐藏层中的某些节点可以没有任何输入数据属于它。所以，K-Means受初始化的影响要比较大。
（2）K-means为每个输入数据找到一个最相似的类后，只更新这个类的参数。SOM则会更新临近的节点。所以K-mean受noise data的影响比较大，SOM的准确性可能会比k-means低（因为也更新了临近节点）。
（3）SOM的可视化比较好。

8、练习

1.设定 topo=rectangular,toroidal=T，重新运行 SOM 比较。

> som.grid2 <- somgrid(xdim = 5, ydim = 5, topo = "rectangular",
+                      neighbourhood.fct = "bubble", toroidal = T)
> flea.som2 <- som(flea.scale, grid = som.grid2, 
+                  rlen = 5000, alpha = c(0.05, 0.01))
> 
> par(mfrow = c(1, 2))
> plot(flea.som, type = "mapping", shape = "straight", main = "Hexagonal")
> plot(flea.som2, type = "mapping", shape = "straight", main = "Rectangular")

比较

2.利用同一个 somgrid，迭代次数改为 10000 次，alpha 设为 c(0.1, 0.001) 来做SOM。

> flea.som3 <- som(flea.scale, grid = som.grid, 
+                  rlen = 10000, alpha = c(0.1, 0.001))
> 
> par(mfrow = c(1, 2))
> plot(flea.som, type = "mapping", shape = "straight", main = "5000次")
> plot(flea.som3, type = "mapping", shape = "straight", main = "10000次")

比较

96-非监督学习之SOM非线性降维