60-R语言中的神经网络

Posted 2023-04-28

参考技术A

《深度学习精要（基于R语言）》学习笔记

机器学习主要用于开发和使用那些从原始数据中学习、总结出来的用于进行预测的算法。
深度学习是一种强大的多层架构，可以用于模式识别、信号检测以及分类或预测等多个领域。
神经网络包括一系列的神经元，或者叫作节点，它们彼此连结并处理输入。神经元之间的连结经过加权处理，权重取决于从数据中学习、总结出的使用函数。一组神经元的激活和权重（从数据中自适应地学习）可以提供给其他的神经元，其中一些最终神经元的激活就是预测。
经常选择的激活函数是sigmoid函数以及双曲正切函数tanh，因为径向基函数是有效的函数逼近，所以有时也会用到它们。
权重是从每个隐藏单元到每个输出的路径，对第i个的输出通过(w_i)表示。如创建隐藏层的权重，这些权重也是从数据中学习得到的。分类会经常使用一种最终变换，softmax函数。线性回归经常使用恒等（identity）函数，它返回输入值。权重必须从数据中学习得到，权重为零或接近零基本上等同于放弃不必要的关系。

R中神经网络相关包：

一旦集群完成初始化，可以使用R或本地主机（127.0.0.1:54321）提供的Web接口与它连接。

如果数据集已经加载到R，使用as.h2o()函数：

如果数据没有载入R，可以直接导入到h2o中：

也可以直接导入网络上的文件：

导入基于图片识别手写体数字，数据集的每一列（即特征），表示图像的一个像素。每张图像都经过标准化处理，转化成同样的大小，所以所有图像的像素个数都相同。第一列包含真实的数据标签，其余各列是黑暗像素的值，它用于分类。

使用caret包训练模型：

生成数据的一组预测，查看柱状图：

跟训练集数据柱状图对比，很明显模型不是最优的。
通过混淆矩阵检查模型性能：

No Information Rate（无信息率）指不考虑任何信息而仅仅通过猜测来决定最频繁的类的准确度期望。在情形“1”中，它在11.16%的时间中发生。P值（P-Value [Acc > NIR]）检验了观测准确度（Accuracy : 0.3674）是否显著不同于无信息率（11.16%）。
Class: 0的灵敏度（Sensitivity）可以解释为：89.07%的数字0被正确地预测为0。特异度（Specificity）可以解释为：95.14%的预测为非数字0被预测为不是数字0。
检出率（Detection Rate）是真阳性的百分比，而最后的检出预防度（detection prevalence）是预测为阳性的实例比例，不管它们是否真的为阳性。
平衡准确度（balanced accuracy）是灵敏度和特异度的平均值。

接下来我们通过增加神经元的个数来提升模型的性能，其代价是模型的复杂性会显著增加：

隐藏神经元的数量从5个增加到10个，样本内性能的总准确度从36.74% 提升到了 65.4%。我们继续增加隐藏神经元的数量：

增加到40个神经元后准确度跟10个神经元的一样，还是65.4%。如果是商业问题，还需要继续调节神经元的数量和衰变率。但是作为学习，模型对数字9的表现比较差，对其他数字都还行。

RSNNS包提供了使用斯图加特神经网络仿真器（Stuttgart Neural Network Simulator ， SNNS）模型的接口，但是，对基本的、单隐藏层的、前馈的神经网络，我们可以使用mlp()这个更为方便的封装函数，它的名称表示多层感知器（multi-layer perceptron）。
RSNNS包要求输入为矩阵、响应变量为一个哑变量的 矩阵 ，因此每个可能的类表示成矩阵列中的 0/1 编码。

通过decodeClassLabels()函数可以很方便的将数据转换为哑变量矩阵。

预测结果的值为1-10，但是实际值为0-9，所以在生成混淆矩阵时，需要先减去1：

RSNNS包的学习算法使用了相同数目的隐藏神经元，计算结果的性能却有极大提高。
函数I()有两个作用：
1.在对data.frame的调用中将对象包含在I()中来保护它，防止字符向量到factor的转换和名称的删除，并确保矩阵作为单列插入。
2.在formula函数中，它被用来禁止将“+”、“-”、“*”和“^”等运算符解释为公式运算符，因此它们被用作算术运算符。

从RSNNS包返回的预测值（pred.ml4）中可以看到，一个观测可能有40%的概率成为“5”，20%的概率成为“6”，等等。最简单的方法就是基于高预测概率来对观测进行分类。RSNNS包有一种称为赢者通吃（winner takes all，WTA）的方法，只要没有关系就选择概率最高的类，最高的概率高于用户定义的阈值（这个阈值可以是0），而其他类的预测概率都低于最大值减去另一个用户定义的阈值，否则观测的分类就不明了。如果这两个阈值都是0（缺省），那么最大值必然存在并且唯一。这种方法的优点是它提供了某种质量控制。
但是在实际应用中，比如一个医学背景下，我们收集了病人的多种生物指标和基因信息，用来分类确定他们是否健康，是否有患癌症的风险，是否有患心脏病的风险，即使有40%的患癌概率也需要病人进一步做检查，即便他健康的概率是60%。RSNNS包中还提供一种分类方法称为“402040”，如果一个值高于用户定义的阈值，而所有的其他值低于用户定义的另一个阈值。如果多个值都高于第一个阈值，或者任何值都不低于第二个阈值，我们就把观测定性为未知的。这样做的目的是再次给出了某种质量控制。

“0”分类表示未知的预测。

通常来说，过拟合指模型在训练集上的性能优于测试集。过拟合发生在模型正好拟合了训练数据的噪声部分的时候。因为考虑了噪声，它似乎更准确，但一个数据集和下一个数据集的噪声不同，这种准确度不能运用于除了训练数据之外的任何数据 — 它没有一般化。
使用RSNNS模型对样本外数据预测：

模型在第一个5000行上的准确度为85.1%，在第二个5000行上的准确度减少为80%，损失超过5%，换句话说，使用训练数据来评价模型性能导致了过度乐观的准确度估计，过度估计是5%。
这个问题我们后面再处理。

R语言中的网络可视化

网络分析适合用来研究多样本或特性间的关系，这类关系通常用互相连接的节点来表示，在可视化中节点一般指代一个样本或特性，连线则代表了样本间或特性间的关系。也就是说，网络的最小单元就是一个两点连线，虽然描述一个网络很直观，但具体到数据结构上就存在一些问题。常规样本数据一般是每一行代表一个样本，每一列代表一个描述样本的维度或特性，样本或特性间的关系并不能展示在原始数据结构里而需要额外计算，所以我们需要将样本-维度的数据结构转成描述网络的点对点数据结构才好可视化。但事实上，更通用的方法是定义一种描述网络的数据类型，然后根据类型定义可视化方法，也就是将原始数据转为网络数据类型，这种类型定义也方便了除可视化外其他针对网络的分析方法开发与使用。

在 R 中，有两个包提供了描述网络的基础类型定义，一个是 network 包，另一个是 igraph 包。这两个包允许用户生成一个专门描述网络的对象，也定义了该对象类型的绘图方法，也就是说，如果你可以直接 plot 一个网络对象，实现快速可视化。很多网络可视化的工具，例如 ggnetwork 包或 ggraph 包或 GGally 包的 ggnet2 函数，都支持输入的对象为 network 包或 igraph 包里定义的网络类型。不过，这里面 ggraph 包还可以基于 tidygraph 包使用 tbl_graph 对象来描述网络关系，几乎完全覆盖了 igraph 包的内容，当然，装这个就得装 tidyverse 全家桶。

不论哪一种对象类型，网络对象一定可以抽象出两张表，一张表保存节点的属性，另一张表保存节点间连线的属性。而对于可视化而言，节点属性表其实就是原始数据框，而连线属性表则要保存我们计算出的节点间关系，例如节点间相关性、距离等。然后很自然每一种对象类型都设计了独立的针对节点与边的赋值方法，有了这些方法就可以自定义一些节点或边的属性方便可视化。

不过，这里很多人手头只有数据框，也就是只有节点那张表，表示两两关系的表还是要自己计算生成的，在图论里这叫做邻接矩阵（ adjacency matrix ）。最简单就是一个相关矩阵或距离矩阵，不过也不一定就是方阵或三角阵，这里面涉及不同的网络类型。多数可视化包都支持你导入一个邻接矩阵或者更原始的点对点数据框来生成网络对象，然后你可以对节点与边进行属性定义，在可视化时指定需要可视化的属性就可以了。很多人（其实就是我）卡在第一步数据导入上就放弃了，但过了第一步能生成网络对象后后面就特别容易进行后续分析。当然，至于说网络可视化的具体图形设计，其实是存在一些预先设定好的美观布局设计方案的例如 Kamada Kawai、Frucherman-Reingold、树形等，可以根据需求进行调整，但不要误导读者，做好图例。

这里特别提一下 qgraph 包，这个包几乎依赖了上面所有提过的包，但这种拿来主义在应用学科的开发者身上并不少见，例如这个包主要是为心理测量学设计的。但很有意思的是，如果你去搜索 R 语言的网络可视化教程，基本都会找到心理测量学或社会科学背景的人写的东西，而且质量很高，例如 Katya Ognyanova 的博客， Sacha Epskamp 的博客， Claudia van Borkulo  的博客还有 psych-networks ，特别最后一个总结并追踪相当多近些年网络科学的研究进展。打个比方，通常我们说节点间的关系，一般就是想到相关性，但两个节点间也可以用是否独立来构建联系而相关只是独立与否的一种，偏相关行不行？或者如果计算二分类变量（社会科学里定量研究常用）间的独立性就需要用到 Ising 模型。同时构建出的网络是不是稳定也需要正则化例如 lasso 或重采样来对变量间关系进行调整，去掉不稳定的联系。另外，如何检测一个网络中的社群？有哪些算法？其实背后也是潜在变量分析的影子。

如果跳到生物信息学领域，有一个 WGCNA 包用的特别多，但很多写教程的人都没搞清楚原理与模型假设。WGCNA 包是基于巴拉巴西的无尺度网络构建的，基本原理就是先对所有基因构建两两相关性矩阵，然后从相关性矩阵中探索出共表达的基因模块，相当于把几千维的基因降维到不到十维且最好能联系上生物学意义例如某个通路啥的。具体计算则是每个模块进行主成分分析（其实是 SVD 分解），然后用第一个主成分作为这个模块的代表对你的研究分组进行差异分析，找出哪个模块有影响然后解释。

这里面核心步骤里有两个坑，第一个在相关性矩阵到基因模块这里，第二个在主成分分析那边。第一个坑是因为其探索模块用了无尺度网络的假设，首先得去选一个幂级数来计算邻接矩阵，这个幂级数是拟合无尺度网络的度分布搞出来的，很多数据本身不符合无尺度网络的度分布假设，所以硬套这个假设是不合适的。第二个坑跟第一个有关系，只有模块内部第一个主成分可解释方差很高才能这么用，但由于第一个坑很多人用了默认值，第二个坑也就只用了第一个主成分，很多时候方差解释连三分之一都不到，虽然能讲故事，但明显是有偏的。当然这个包里也是涵盖了很多对于用户而言天书级的概念，很多人不求甚解套默认值也把文章给发了，完全就当神奇降维盒子在用了。

其实说白了网络分析是另一层意义上的因子分析，起一个降维作用，只是降维方式不是简单的线性组合而是引入了图论的一些统计量罢了，但我看到很多人用起来套代码，解释上完全就是意会，特别是代谢组学分析里会出现套基因组学的分析方法而不验证假设盲目追求自动化的默认值。不过我也看到了很多基于图论的生物信息学文章，很多想法非常超前但引用非常少且真正生产实验数据的人基本看不懂或不会去读，这就还不如心理测量学那边研究人员的学科内科普做得好。

这里说个题外话，非统计与计算机科学背景开发者写的包其实非常多，多数都不在 CRAN 上。他们一般只会用基础R包函数来实现自己想要的功能，如果没有就会去依赖其他包来调用函数，很少用 Rcpp ，不开并行计算，基本不关心速度，用 S3 对象而不是 S4 ，这倒是科研编程的日常状态。有时候读他们的源码有种见字如面的感觉：有的人明显是其他语言转过来的，有次读一个包的代码怎么看怎么别扭，后来发现这个开发者的母语是 java ，很多定义方式都是那边继承过来的。有的人注释掉的代码中其实有另一重意思，跟最终版的实现方式完全不同，源码里保留了很多进化遗迹，类似化石。有的人严谨，每个函数都写测试，文档明显打磨过语言。有的人飘逸，通篇找不到注释，很多编码风格都不一致，感觉是爆栈网复制过来的。有的人很明显是 tidyverse 风格出现后才开始学的 R ，对基础函数用法非常不熟。非统计与计算机科学开发者的代码通常存在很多不严谨的地方，没有经过软件工程的训练，更多是为了解决特定目的而快速实现的，不过很多代码展示出了丰富多彩的想象力。

下面生成几个网络并做下基础可视化，这里不会用最常见的那种点对点数据结构，因为这个东西是需要从原始数据生成的，很少有原始数据本身就是这种关系结构，而且此处也不涉及 ggplot2 风格的绘图包，用基础绘图系统来做，函数统一为 plot ，当然不同的对象类型会有不同的绘图参数，此处也不会区分有向图与无向图，按默认值来。理解清楚了这些最基础的代码实现过程，自定义出图上就会很轻松了。

network 版

set.seed(110)library(network)# 生成一个3节点网络net <- network.initialize(3)# 画出来plot(net,vertex.cex=10)

# 添加一条边add.edge(net,2,3) # 画出来plot(net,vertex.cex=10, displaylabels=T)

# 添加两个点add.vertices(net,2)# 画出来plot(net,vertex.cex=10, displaylabels=T)

# 模拟一个5*12的数据框df <- matrix(rnorm(60),5)# 用邻接矩阵直接生成网络dfcor <- cor(df)# 去掉低相关性边dfcor[dfcor<0.5] <- 0netcor <- as.network(dfcor,matrix.type = 'adjacency')plot(netcor)

# 增加节点/边属性set.vertex.attribute(netcor, "class", length(netcor$val):1)set.edge.attribute(netcor,"color",length(netcor$mel):1)# 可视化属性plot(netcor,vertex.cex=5,vertex.col=get.vertex.attribute(netcor,"class"),edge.col=get.edge.attribute(netcor,'color'))

igraph 版

set.seed(110)library(igraph)# 生成一个3节点网络net <- graph.empty(n=3, directed=TRUE)# 画出来plot(net)

# 添加两条边new_edges <- c(1,3, 2,3)net <- add.edges(net, new_edges)# 画出来plot(net)

# 添加两个点net <- add.vertices(net, 2)# 画出来plot(net)

# 模拟一个5*12的数据框df <- matrix(rnorm(60),5)# 用邻接矩阵直接生成网络dfcor <- cor(df)# 去掉低相关性边dfcor[dfcor<0.5] <- 0net <- graph.adjacency(dfcor,weighted=TRUE,diag=FALSE)plot(net)

# 增加节点/边属性V(net)$name <- letters[1:vcount(net)]E(net)$color <- "red"E(net)[ weight < 0.7 ]$width <- 2E(net)[ weight < 0.7 ]$color <- "green"# 可视化属性plot(net)

网络可视化只是网络分析的基础，也只有理解了其基础才能更好进行下一步的分析，很多基于网络稳定性分析还有网络群组分析都是可以基于更基础的概率图模型来进行。这些分析都有明确的背景问题来源，但涉及的知识点非常多，从统计物理到图论到随机过程，不过如果带着自己的问题去探索，总会有新的发现。

关于作者

于淼，现任编辑部主编，中科院生态环境研究中心博士，美国西奈山伊坎医学院博后，研究方向为环境化学与暴露组学数据分析，博客https://yufree.cn

审稿：夏骁凯

编辑：向悦

统计之都：专业、人本、正直的中国统计学社区。

往期推送：进入统计之都会话窗口，点击右上角小人图标，查看历史消息即可。