【R语言】非监督学习之降维

非监督学习:不知道数据集中数据、特征之间的关系,而是要根据聚类或一定的模型得到数据之间的关系。回归(Regression)针对于连续型变量的;分类(Classification) 针对离散型的,输出的结果是有限
的。

为何要降维?方便可视化探索;减轻维度诅咒;缓解共线性。
降维方法:PCA(线性), t-SNE, UMAP, SOM, LLE

各类包汇总:
GGally包作为ggplot2扩展包之一,可以方便的显示配对图矩阵,散点图矩阵,平行坐标图,生存图,以及绘制网络。
factoextra包可以轻松提取和可视化探索性多变量数据分析的输出,包括PCA等。
Rtsne专门进行t-SNE降维分析。
umap进行UMAP降维可视化分析。

01 PCA主成分分析

PCA主成分分析视频讲解
主成分分析(Principal Component Analysis, 简写PCA)是将多个指标化为少数几个综合指标的一种统计分析方法,实质是一种降维技术, 利用正交变换把一系列可能线性相关的变量转换为一组线性不相关的新变量,也称为主成分,从而利用新变量在更小的维度下展示数据的特征。
原理:PCA的核心就是找坐标系,第一步通过投影到X,Y轴的距离,分别找到数据在X,Y轴投影的中心点,两个轴上的中心点所组成的坐标就是这组数据的中心,也是新坐标系的原点;第二步,穿过“新原点”的直线进行旋转,找到所有数据点投影在这条直线上的投影点到新原点的距离的平方之和最大,这条直线就是PC1,与之垂直的直线就是PC2,这就找到了新的坐标轴。投影点到原点之间距离的平方和称为PC1的特征值,PC1特征值的平方根为PC1奇异值,每个变量的比例称为载荷得分。
R语言中的PCA分析与可视化

pacman::p_load(dplyr, ggplot2, GGally, factoextra, mclust)

step1:构建PCA模型

note <- as_tibble(mclust::banknote)   #把banknote变成一个矩阵
note
#各自变量与因变量(Status)之间、各自变量之间的相关性、分布情况
ggpairs(note,aes(color = Status)) + theme_bw()
#【最重要】坐标轴变换,中心化,消除数据量纲不同的影响
pca <- dplyr::select(note,-Status) %>%    #去掉真假钞列
  prcomp(center = T, scale = T)   
pca
summary(pca)
pca$sdev     #特征值的开方,各主成分的标准差
pca$rotation    #特征向量,回归系数
pca$x   #样本得分score,亦即样本数据在新的坐标轴上的投影(此处为PC1-PC6共6根轴上的投影)
#变量荷载计算
map_dfc(1:6,~pca$rotation[,.]*sqrt(pca$sdev^2)[.])
# Map 旋转之后的轴和之前轴的关系, ”-0.9219364” PC2与length强烈的负相关;“0.847”显示,PC1与Diagonal强烈正相关。

> note
# A tibble: 200 x 7
   Status  Length  Left Right Bottom   Top Diagonal
   <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>
 1 genuine   215.  131   131.    9     9.7     141 
 2 genuine   215.  130.  130.    8.1   9.5     142.
 3 genuine   215.  130.  130.    8.7   9.6     142.
 4 genuine   215.  130.  130.    7.5  10.4     142 
 5 genuine   215   130.  130.   10.4   7.7     142.
 6 genuine   216.  131.  130.    9    10.1     141.
 7 genuine   216.  130.  130.    7.9   9.6     142.
 8 genuine   214.  130.  129.    7.2  10.7     142.
 9 genuine   215.  129.  130.    8.2  11       142.
10 genuine   215.  130.  130.    9.2  10       141.
# ... with 190 more rows

ggpairs()用法
ggpairs(note,aes(color = Status),
upper = list(continuous = "density"),
lower = list(continuous = "points"),
diag = list(continuous = "densityDiag")) +
theme_bw()
注:上三角plot图为相关系数和显著性检验结果图,下三角plot图为散点图,对角线plot图为密度图。
theme_bw() ggplot2::theme_grey中的一种主题,去掉灰色背景。
也可使用ggscatmat() ,ggscatmat() 画图快,可指定列,如:
ggscatmat(data,columns = 1:ncol(data),color = NULL,alpha = 1,corMethod = "pearson" )
ggscatmat(note,color = "Status")
ggscatmat(note,color = "Status",columns = 2:7)

ggpairs(note,aes(color = Status)) + theme_bw()

prcomp()是对变量矩阵(相关矩阵)采用SVD方法计算其奇异值(原理上是特征值的平方根),函数帮助中描述为函数结果中的sdev。
输入参数:变量矩阵(x),中心化(center,默认为true,当有比较多的异常点,均值受到偏离时不用center=T),标准化(scale,默认为false,建议改为true,scale数值差异化大,做方差归一化操作,防止量级大占主导),主成份个数(rank)。
函数输出:sdev(各主成分的标准差),rotation(特征向量,回归系数),x(score得分矩阵)。

> pca
Standard deviations (1, .., p=6):
[1] 1.7162629 1.1305237 0.9322192 0.6706480 0.5183405 0.4346031

Rotation (n x k) = (6 x 6):
                  PC1         PC2         PC3        PC4
Length    0.006987029 -0.81549497  0.01768066  0.5746173
Left     -0.467758161 -0.34196711 -0.10338286 -0.3949225
Right    -0.486678705 -0.25245860 -0.12347472 -0.4302783
Bottom   -0.406758327  0.26622878 -0.58353831  0.4036735
Top      -0.367891118  0.09148667  0.78757147  0.1102267
Diagonal  0.493458317 -0.27394074 -0.11387536 -0.3919305
                PC5         PC6
Length   -0.0587961  0.03105698
Left      0.6394961 -0.29774768
Right    -0.6140972  0.34915294
Bottom   -0.2154756 -0.46235361
Top      -0.2198494 -0.41896754
Diagonal -0.3401601 -0.63179849


> summary(pca)       #清楚看出降到3维
Importance of components:   #组成元素的重要性
                          PC1    PC2    PC3     PC4     PC5
Standard deviation     1.7163 1.1305 0.9322 0.67065 0.51834   #标准差
Proportion of Variance 0.4909 0.2130 0.1448 0.07496 0.04478   #方差比例
Cumulative Proportion  0.4909 0.7039 0.8488 0.92374 0.96852   #累计比例
                           PC6
Standard deviation     0.43460
Proportion of Variance 0.03148
Cumulative Proportion  1.00000

不同主成分对数据差异的贡献和主成分与原始变量的关系:

  1. 主成分标准差的平方为特征值,其含义为每个主成分可以解释的数据差异,计算方式为 eigenvalues = (pca$sdev)^2;
  2. 每个主成分可以解释的数据差异的比例为percent_var = eigenvalues*100/sum(eigenvalues)
  3. 可以使用summary(pca)获取以上两条信息。

通过以下信息可以判断主成分分析的质量:成功的降维需要保证在前几个为数不多的主成分对数据差异的解释可以达到80-90%。

指导选择主成分的数目:

  1. 选择的主成分足以解释的总方差大于80% (方差比例碎石图)
  2. 从前面的协方差矩阵可以看到,自动定标(scale)的变量的方差为1 (协方差矩阵对角线的值)。待选择的主成分应该是那些方差大于1的主成分,即其解释的方差大于原始变量(特征值碎石图,方差大于1,特征值也会大于1,反之亦然)。

经过主成分分析后,六个自变量变为了PC1到PC6,累积方差解释比(Cumulative Proportion)依次为0.4909、0.7039、0.8488、0.92374、0.96852、1.00000,意思是如果只选择PC1到PC4四个变量,但是却能解释整个数据方差的92.374%,从而达到降维的目的。
成功的降维需要保证在前几个为数不多的主成分对数据差异的解释可以达到80-90%。后面的建模直接使用PC1到PC4的数据,而不再使用原始数据。

> map_dfc(1:6,~pca$rotation[,.]*sqrt(pca$sdev^2)[.])
New names:
* NA -> ...1
* NA -> ...2
* NA -> ...3
* NA -> ...4
* NA -> ...5
* ...
# A tibble: 6 x 6
     ...1   ...2    ...3    ...4    ...5    ...6
    <dbl>  <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
1  0.0120 -0.922  0.0165  0.385  -0.0305  0.0135
2 -0.803  -0.387 -0.0964 -0.265   0.331  -0.129 
3 -0.835  -0.285 -0.115  -0.289  -0.318   0.152 
4 -0.698   0.301 -0.544   0.271  -0.112  -0.201 
5 -0.631   0.103  0.734   0.0739 -0.114  -0.182 
6  0.847  -0.310 -0.106  -0.263  -0.176  -0.275

# 主成分构成的新矩阵
> pca$rotation   
                  PC1         PC2         PC3        PC4
Length    0.006987029 -0.81549497  0.01768066  0.5746173
Left     -0.467758161 -0.34196711 -0.10338286 -0.3949225
Right    -0.486678705 -0.25245860 -0.12347472 -0.4302783
Bottom   -0.406758327  0.26622878 -0.58353831  0.4036735
Top      -0.367891118  0.09148667  0.78757147  0.1102267
Diagonal  0.493458317 -0.27394074 -0.11387536 -0.3919305
                PC5         PC6
Length   -0.0587961  0.03105698
Left      0.6394961 -0.29774768
Right    -0.6140972  0.34915294
Bottom   -0.2154756 -0.46235361
Top      -0.2198494 -0.41896754
Diagonal -0.3401601 -0.63179849

#每个主成分可以解释的数据差异,特征值eigenvalues = (pca$sdev)^2
> pca$sdev   
[1] 1.7162629 1.1305237 0.9322192 0.6706480 0.5183405 0.4346031

prcomp函数输出有sdev(各主成份的奇异值:原理上是特征值的平方根),rotation(特征向量,回归系数),x(score得分矩阵)。


step2:可视化

#双标图,看新坐标与原变量关系
fviz_pca_biplot(pca, label = "var")
#变量负载图
fviz_pca_var(pca)
# 特征值碎石图,一般选择方差大于1的主成分
fviz_screeplot(pca, addlabels = T, choice = "eigenvalue")
# 方差比例碎石图,一般选择总方差大于80%
fviz_screeplot(pca,addlabels = T, choice = "variance")
# 选择前两个主成分,加入因变量
notePCA <- note %>%
  mutate(PCA1 = pca$x[,1],PCA2 = pca$x[,2])
# 画图
ggplot(notePCA,aes(PCA1,PCA2,color = Status))+
  geom_point()
双标图

在空间上,PCA可以理解为把原始数据投射到一个新的坐标系统,第一主成分为第一坐标轴(Dim1),它的含义代表了原始数据中多个变量经过某种变换得到的新变量的变化区间;第二成分为第二坐标轴(Dim2),代表了原始数据中多个变量经过某种变换得到的第二个新变量的变化区间。这样把利用原始数据解释样品的差异转变为利用新变量解释样品的差异。
这种投射方式会有很多,为了最大限度保留对原始数据的解释,一般会用最大方差理论或最小损失理论,使得第一主成分有着最大的方差或变异数 (就是说其能尽量多的解释原始数据的差异);随后的每一个主成分都与前面的主成分正交,且有着仅次于前一主成分的最大方差 (正交简单的理解就是两个主成分空间夹角为90°,两者之间无线性关联,从而完成去冗余操作)。

变量负载图

在合并冗余原始变量得到主成分过程中,会发现某些原始变量对同一主成分有着相似的贡献,也就是说这些变量之间存在着某种相关性,为相关变量。同时也可以获得这些变量对主成分的贡献程度。

特征值碎石图,一般选择方差大于1的主成分
方差比例碎石图,一般选择总方差大于80%
> notePCA
# A tibble: 200 x 9
   Status  Length  Left Right Bottom   Top Diagonal   PCA1    PCA2
   <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl>    <dbl>  <dbl>   <dbl>
 1 genuine   215.  131   131.    9     9.7     141  -1.74  -1.65  
 2 genuine   215.  130.  130.    8.1   9.5     142.  2.27   0.537 
 3 genuine   215.  130.  130.    8.7   9.6     142.  2.27   0.107 
 4 genuine   215.  130.  130.    7.5  10.4     142   2.28   0.0874
 5 genuine   215   130.  130.   10.4   7.7     142.  2.63  -0.0391
 6 genuine   216.  131.  130.    9    10.1     141. -0.757 -3.08  
 7 genuine   216.  130.  130.    7.9   9.6     142.  2.51  -1.22  
 8 genuine   214.  130.  129.    7.2  10.7     142.  2.70   1.13  
 9 genuine   215.  129.  130.    8.2  11       142.  2.03   0.314 
10 genuine   215.  130.  130.    9.2  10       141. -0.317 -1.30  
# ... with 190 more rows
主成分分类(红色表示假钞,蓝色表示真钞)

step3:应用PCA检验效果

如果新来了两张钞票,可以通过可视化辨别是真钞还是假钞。

#增加2个新数据
newNotes <- tibble(
  Length = c(214, 216),
  Left = c(130, 128),
  Right = c(132, 129),
  Bottom = c(12, 7),
  Top = c(12, 8),
  Diagonal = c(138, 142)
)
#预测,转化成tibble
newPCA <- predict(pca,newNotes) %>% as_tibble()
#原始图真钞伪钞+新加入两个点
ggplot(notePCA,aes(PCA1,PCA2,color = Status))+
  geom_point()+
  stat_ellipse(level = 0.90)+
  geom_point(data = newPCA,aes(PC1,PC2),col="black",size=4)

> newPCA
# A tibble: 2 x 6
    PC1    PC2    PC3   PC4   PC5   PC6
  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
1 -4.73  2.00  -0.106 -1.66 -3.20  1.62
2  6.47 -0.892 -0.821  3.47 -1.84  2.34
预测

PCA的优劣
优势:
1.直接利用原有变量建立新轴;
2.新数据可以直接投影;
3.纯数学变换,计算量不大。
劣势:
1.高维到低维可能非线性;
2.无法处理分类变量;
3.降低维度需要人工判断。




02 t-SNE非线性降维

t-SNE:t-distributed stochastic neighbor embedding:t分布随机邻域嵌入是一种用于探索高维数据的非线性降维算法。它将多维数据映射到适合于人类观察的两个或多个维度。
t-SNE非线性降维算法通过基于具有多个特征的数据点的相似性识别观察到的簇来在数据中找到模式。本质上是一种降维和可视化技术。另外t-SNE的输出可以作为其他分类算法的输入特征。

1、降维步骤
t-SNE:将数据点之间的相似度转换为概率。原始空间中的相似度由高斯联合概率表示,嵌入空间的相似度由“学生t分布”表示。
第一步:计算数据集中每行与其他行的距离(默认为欧式距离),转换为概率向量;
第二步:对每一行重复操作,得到概率矩阵;
第三步:沿两条新轴用学生t分布对数据随机化;
第四步:逐渐迭代,通过最小化KL散度,使得二维空间的新概率矩阵尽可能接近原高维空间的。

2、t-SNE特点
相较于正态分布,使用t分布能更好地分散可能的数据簇,更易识别;基于所实现的精度,将t-SNE与PCA和其他线性降维模型相比,结果表明t-SNE能够提供更好的结果,这是因为算法定义了数据的局部和全局结构之间的软边界。
缺点:1.t-SNE在低维容易保持局部结构,但全局结构未明确保留;2.计算费时,计算时间随簇数显著增加,在数百万个样本数据集中可能需要几个小时,而PCA可以在几秒钟或几分钟内完成;3.无法像PCA一样投影新数据;4.簇间距离意义不大。

3、超参数
perplexity:控制距离转化为概率的分布:局部结构 5-30-50 全局结构,数据集越大,需要参数值越大,取值必须小于 (nrow(data) - 1 )/ 3 ;
theta:权衡速度与精度:精确 0-0.5-1 最快;
dims :设置降维之后的维度,默认值为2;
eta:学习率,越低越精确,越大更少迭代,默认值200;
max_iter:最多迭代次数:默认值1000。

pacman::p_load(tidyr, dplyr, Rtsne, ggplot2, mclust)

noteTsne <- mclust::banknote %>% 
  dplyr::select(-Status) %>%
Rtsne(perplexity = 30, theta = 0, max_iter = 5000, verbose = T)
# 降维后的数据
head(noteTsne$Y)

## 将数值型变量标准化,将降维后的数据加入数据框
noteTsne1 <- banknote %>% 
  mutate_if(.funs = scale, .predicate = is.numeric, scale = F) %>%   
  mutate(tSNE1 = noteTsne$Y[,1], tSNE2 = noteTsne$Y[,2]) %>%
  gather(key = "Variable", value = "Value", c(-tSNE1, -tSNE2, -Status))
# 保留tsne1, tsne2, Status列,将其他列宽表变长表,便于画图
#也可写成:pivot_longer(2:7,names_to ="Variable",values_to = "Value")
#查看每个特征的降维效果图
ggplot(noteTsne1, aes(tSNE1, tSNE2, col = Value, shape = Status)) +
  facet_wrap(~ Variable) +
  geom_point(size = 2) +
  scale_color_gradient(low = "dark blue", high = "cyan") +   # 梯度填充颜色
  theme_minimal()
> head(noteTsne$Y)
           [,1]      [,2]
[1,] -0.7772086 -14.33556
[2,] -5.1238226 -19.98535
[3,] -3.4972041 -20.19085
[4,] -9.1770545 -15.85892
[5,] -0.4497331 -21.11881
[6,] -1.8790969 -14.02269

gather()用于长数据和宽数据之间的转换,参考R语言数据格式 长数据 和 宽数据 之间的转换R语言之数据处理常用包
使用gather()将宽数据转换为长数据:
gather(data, key, value, …, na.rm = FALSE, convert = FALSE)
data:需要被转换的宽形表
key:将原数据框中的所有列赋给一个新变量key
value:将原数据框中的所有值赋给一个新变量value
…:可以指定哪些列聚到同一列中
na.rm:是否删除缺失值

> head(noteTsne1)
   Status      tSNE1     tSNE2 Variable  Value
1 genuine -0.7772086 -14.33556   Length -0.096
2 genuine -5.1238226 -19.98535   Length -0.296
3 genuine -3.4972041 -20.19085   Length -0.096
4 genuine -9.1770545 -15.85892   Length -0.096
5 genuine -0.4497331 -21.11881   Length  0.104
6 genuine -1.8790969 -14.02269   Length  0.804
特征分面图

查看每个特征的降维效果图,其中Diagonal特征比较容易区分两种钞票,其他特征分辨能力相对较差,因为颜色梯度填充比较混杂。




03 UMAP非线性降维

UMAP(Uniform Manifold Approximation and Projection,统一流形逼近与投影)是一种非线性降维的算法,相对于t-SNE,UMAP算法更加快速,该方法的原理是利用流形学和投影技术,达到降维目的,首先计算高维空间中的点之间的距离,将它们投影到低维空间,并计算该低维空间中的点之间的距离。然后,它使用随机梯度下降来最小化这些距离之间的差异。
相比 t-SNE 的优势:
1.速度快得多
2.确定性算法
3.保留双结构

UMAP的超参数
n_neighbors:控制模糊搜索区域的半径:更少邻域到更多邻域;通常其设置在2-100之间;
n_components降维的维数大小,默认是2,其范围最好也在2-100之间;
min_dist:低维下允许的行间最小距离:更集中到更分散;
metric:选择距离的测度方法:欧氏距离euclidean、曼哈顿距离manhattan等;
n_epochs:优化步骤的迭代次数;
verbose = T能看到过程;

pacman::p_load(dplyr, umap, ggplot2)

note <- as_tibble(mclust::banknote)
noteUmap <- note %>% 
  dplyr::select(-Status) %>% 
  as.matrix() %>% 
  umap(n_neibours = 7, min_dist = 0.1, 
       metric = "manhattan", n_components = 2,n_epochs = 200, verbose = T)
# 查看降维后的数据
noteUmap$layout %>% head()
# 将数值型变量标准化,将降维后的数据加入数据框
noteUmap1 <- note %>% 
  mutate_if(.funs = scale, .predicate = is.numeric, scale = F) %>% 
  mutate(UMAP1 = noteUmap$layout[,1], UMAP2 = noteUmap$layout[,2]) %>% 
  gather(key = "Variable", value = "Value", c(-UMAP1, -UMAP2, -Status))
#画特征分面图
ggplot(noteUmap1, aes(UMAP1, UMAP2, col = Value, shape = Status)) +
  facet_wrap(~ Variable) +
  geom_point(size = 2) +
  scale_color_gradient(low = "dark blue", high = "cyan")
#增加2个新数据,预测模型在新数据情况下的值
newNotes <- tibble(
  Length = c(214, 216),
  Left = c(130, 128),
  Right = c(132, 129),
  Bottom = c(12, 7),
  Top = c(12, 8),
  Diagonal = c(138, 142)
)
predict(noteUmap, newNotes)
> noteUmap$layout %>% head()
         [,1]      [,2]
[1,] 3.480542 -4.748042
[2,] 3.763205 -6.797228
[3,] 4.385722 -6.779360
[4,] 0.670081 -6.599834
[5,] 5.057274 -6.488332
[6,] 3.257957 -4.742613

> noteUmap1
# A tibble: 1,200 x 5
   Status  UMAP1 UMAP2 Variable    Value
   <fct>   <dbl> <dbl> <chr>       <dbl>
 1 genuine 3.48  -4.75 Length   -0.0960 
 2 genuine 3.76  -6.80 Length   -0.296  
 3 genuine 4.39  -6.78 Length   -0.0960 
 4 genuine 0.670 -6.60 Length   -0.0960 
 5 genuine 5.06  -6.49 Length    0.104  
 6 genuine 3.26  -4.74 Length    0.804  
 7 genuine 3.58  -6.62 Length    0.604  
 8 genuine 1.10  -7.00 Length   -0.396  
 9 genuine 0.501 -6.31 Length    0.00400
10 genuine 3.32  -4.62 Length    0.304  
# ... with 1,190 more rows
特征分面图
> predict(noteUmap, newNotes)
[2021-06-17 23:31:43]  creating graph of nearest neighbors
[2021-06-17 23:31:43]  creating initial embedding
[2021-06-17 23:31:43]  optimizing embedding
[2021-06-17 23:31:44]  done
        [,1]      [,2]
1 -0.4107805  6.933556
2  3.9133754 -7.397181

t-SNE 及 UMAP 的优点:
1.可学习非线性规律;
2.能更好地分簇;
3.UMAP 可预测新数据;
4.UMAP 计算量不太大;
5.UMAP 保持双结构。
t-SNE 及 UMAP 的缺点:
1.新轴线难解释;
2.t-SNE 不可预测新数据;
3.t-SNE 计算量太大;
4.t-SNE 不保持全局结构;
5.不能直接处理分类变量。




04 SOM非线性降维

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

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

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

pacman:: p_load(dplyr, kohonen, GGally)

data(flea)
ggpairs(flea, aes(color = species))
#查看6个变量之间的关系,相关性较大可以降维
#或者使用ggscatmat(flea, column =2:6, color = "species")
#建立SOM网络(碗阵)
somGrid <- somgrid(xdim = 5, ydim = 5, topo = "hexagonal",
                   neighbourhood.fct = "bubble", toroidal = F)
#标准化,有助于使每个特征对于计算相似度(距离)的贡献相同
flea1 <- as_tibble(flea)
flea2 <- dplyr::select(flea1, -species) %>% 
  scale()
# rlen:迭代次数
# alpha:学习速率,默认从0.05下降到0.01
fleaSom <- som(flea2, grid = somGrid, rlen = 5000, alpha = c(0.05, 0.01))
#fleaSom$data 中心化后的数据
#fleaSom$unit.classif  查看降维后的结果,碗阵
#根据plot.type中的类型,依次画图,并按2×3排列
par(mfrow = c(2,3))   
plotType <- c("codes", "changes", "counts", "quality", "dist.neighbours", "mapping")
walk(plotType, ~plot(fleaSom, type = ., shape = "straight"))
# property:属性值,main:标题
getCodes(fleaSom) %>% 
  as_tibble() %>% 
 purrr:: iwalk(~plot(fleaSom, type = "property", property = .,
              main = .y, shape = "straight"))
#"property":每个单元的属性可以计算并显示在颜色代码中,用来可视化一个特定对象与映射中所有单元的相似性。

#拓展:clusters层级聚类
clusters<-cutree(hclust(dist(fleaSom$codes[[1]],method = "manhattan")),3)
comClusters<-map_dbl(clusters,~{if(.==1)3  else if(.==2)2 else 1})
plot(fleaSom,type="mapping",pch=21,bg=nodeCol[as.numeric(flea1$sp
ecies)],shape = "round",bgcol=nodeCol[as.integer(comClusters)])
add.cluster.boundaries(fleaSom,comClusters)

#新数据上应用SOM
#给新数据
newDataFlea <- tibble::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)) 
# 使用之前的属性值标准化
newDataFlea <- newDataFlea %>% 
  scale(center = attr(flea2, "scaled:center"),
        scale = attr(flea2, "scaled:scale"))
#预测
predicted <- predict(fleaSom, newDataFlea)
par(mfrow=c(1,1))   #恢复1个排列
plot(fleaSom, type = "mapping", classif = predicted, shape = "round")

> 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 ...
变量相关性

用法:somgrid(xdim = 8, ydim = 6, topo = c("rectangular", "hexagonal"),neighbourhood.fct = c("bubble", "gaussian"), toroidal = FALSE)
xdim,ydim和topo是SOM网络的最基本拓扑信息,包括网络大小和拓扑形状,topo ="hexagonal"六边形,"rectangular"长方形
neighbourhood.fct是邻域函数,可以选高斯函数和“泡泡”bubble函数;
toroidal = T把默认的SOM拓扑图变成环形,此时SOM拓扑相当于是把环形切割后并展开;toroidal的作用看情况,主要改变的是SOM的拓扑关系,对训练结果本身并没有影响。

> somGrid
$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"
> str(flea2)
 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" ...
> fleaSom
SOM of size 5x5 with a hexagonal topology.
Training data included.
依次画图

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

针对每个特征会有一幅图
重新聚类
新数据映射



05 LLE非线性降维

LLE: Locally Linear Embedding,局部线性嵌入,属于流形学习(Manifold Learning)的一种,其假设数据在较小的局部是线性的,也就是说,某一个数据可以由它邻域中的几个样本来线性表示。LLE 非常适合处理卷起或扭曲状的数据。

LLE步骤总结
1.计算行间距,设定超参数k。
2.对一行选出其最近的k行,表示为其线性组合,该线性组合系数为权重。
3.对每行重复操作,使得数据在2或3维空间中(近乎)保持该线性组合关系。

pacman:: p_load(dplyr, lle, ggplot2, plot3D, plot3Drgl)

#LLE 数据展示
data(lle_scurve_data)
colnames(lle_scurve_data) <- c("x","y","z")
scurve <- as_tibble(lle_scurve_data)
scatter3D(x = scurve$x, y = scurve$y, z = scurve$z, pch = 19,
          bty = "b2", colkey = F, theta = 35, phi = 10,
          col = ramp.col(c("red", "cyan")))
# 让3D图像可以用鼠标转动
plot3Drgl::plotrgl()   

#LLE 超参数选择
llek <- calc_k(scurve, m = 2, kmin = 1, kmax = 20, 
               parallel = T, cpus = parallel::detectCores())
# 找出使rho最小的K值
llekStar <- dplyr::filter(llek, rho == min(llek$rho))
#使用最优的K值降维
lleCurve <- lle(scurve, m = 2, k = llekStar$k)
#使用降维后的数据画图
scurve <- scurve %>% 
  mutate(
    LLE1 = lleCurve$Y[,1],
    LLE2 = lleCurve$Y[,2]
  )
ggplot(scurve, aes(LLE1,LLE2, col = z)) +   # 按z值大小填充渐变色
  geom_point() +
  scale_color_gradient(low = "cyan", high = "red")
数据三维图

除了维数,k (近邻数量)是唯一需要确定的超参数。
K可以通过函数计算出来:calc_k()
m 表示维数,通常2 或 3;
kmin,kmax 决定 k 取值域;
parallel,是否多核运行,默认为否;
cpus 指定使用 cpu 核数;

llek:红线代表最优解,rho越小越好
> llekStar
   k       rho
1 17 0.1468803
二维图

练习

1.给 flea 的LLE 二维散点图上加上代表种类的椭圆。

data(flea, package = "GGally")
# 去掉species变量,标准化
flea.scale <- flea %>% 
  as_tibble() %>%
  dplyr::select(-species) %>% 
  scale()

flea.k <- calc_k(flea.scale, m = 2, kmin = 1, kmax = 20, plotres = F,
                 parallel = T, cpus = parallel::detectCores())

flea.lle <- lle(flea.scale, m = 2, k = filter(flea.k, rho == min(flea.k$rho))$k)
#画图
p1 <- tibble(LLE1 = flea.lle$Y[, 1],
             LLE2 = flea.lle$Y[, 2],
             species = flea$species) %>% 
    ggplot(aes(LLE1, LLE2, col = species)) +
    geom_point(size = 2) +
    stat_ellipse(level = 0.93) +
    theme_bw() +
    theme(legend.position = "top")
p1
标准化后的数据

2.设定三维来做LLE,并用 scatter3D() 函数来作图。

scatter3D(x = flea.lle$Y[, 1], y = flea.lle$Y[, 2], z = flea$head, pch = 19,
          bty = "b2", colkey = F, theta = 35, phi = 10,
          col = ramp.col(c("red", "cyan")))
三维图

3.不使用scale函数直接做二维LLE,并绘制散点图。

flea2 <- flea %>% as_tibble() %>% dplyr::select(-species)
 
flea.k2 <- calc_k(flea2, m = 2, kmin = 1, kmax = 20, plotres = F,
                  parallel = T, cpus = parallel::detectCores())
flea.lle2 <- lle(flea2, m = 2, k = filter(flea.k, rho == min(flea.k$rho))$k)

p2 <- tibble(LLE1 = flea.lle2$Y[, 1],
             LLE2 = flea.lle2$Y[, 2],
             species = flea$species) %>% 
         ggplot(aes(LLE1, LLE2, col = species)) +
         geom_point(size = 2) +
         stat_ellipse(level = 0.93) +
         theme_bw() +
         theme(legend.position = "none")
#与标准化的图拼在一起对比
p_load(patchwork)
p1 + p2
标准化与非标准化数据对比

SOM与LLE
1.优点:
1.非线性还原算法。
2.新数据可以映射到SOM上。
3.训练成本相当不高。
4.LLE算法可重复。
2.缺点:
1.不能处理分类变量。
2.不能直接用原始变量解释。
3.对不同尺度的数据很敏感。
4.新数据不能映射到LLE上。
5.不一定保留数据的全局结构。
6.SOM算法每次都会产生不同的结果。
7.SOM在大型数据集上效果更好。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 193,812评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,626评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,144评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,052评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 60,925评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,035评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,461评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,150评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,413评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,501评论 2 307
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,277评论 1 325
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,159评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,528评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,868评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,143评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,407评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,615评论 2 335

推荐阅读更多精彩内容