R语言实现非负矩阵分析

著名的科学杂志《Nature》于1999年刊登了两位科学家D.D.Lee和H.S.Seung对数学中非负矩阵研究的突出成果。该文提出了一种新的矩阵分解思想――非负矩阵分解(Non-negative Matrix Factorization，NMF)算法，即NMF是在矩阵中所有元素均为非负数约束条件之下的矩阵分解方法。

上图引自网络（有出处请告知），NMF的思想：V=WH（W权重矩阵、H特征矩阵、V原矩阵），通过计算从原矩阵提取权重和特征两个不同的矩阵出来。属于一个无监督学习的算法，其中限制条件就是W和H中的所有元素都要大于0。

今天我们给大家讲下在R语言中是如何实现的。先来看下NMF包的安装。这个有点麻烦，我们首先必须要把我们的R版本升级到3.6及以上，因为有个rngtools的依赖包所需要的环境是3.6及以上。

另外，Biobase包需要从Bioconductor进行安装需要在安装NMF前安装好，如果没有安装则会显示下面信息：

以上都按照要求安装好后则会出现下面的成果信息：

那么接下来我们运行下NMF自带的例子：

meth <- nmfAlgorithm(version='R')
meth <- c(names(meth), meth)
meth
data(esGolub)
res <- nmf(esGolub, 3, meth,seed=123456)

接下来我们看下nmf函数的主要参数：

Rank：就是因式分解的级别。其中自带了计算最优等级的函数nmfEstimateRank：

其中重要参数：Range：等级的分布范围。例seq(1,5)；Method:就是NMF算法的选择；Nrun：范围内每个值的迭代次数；Model：如果有已经构建好的模型，那就是直接把模型写在这里，优化等级计算。构建模型的函数是nmfModel（rank,c(features,samples)）或者是nmfModel（rank,data,W,H）。

Methods：就是对应的NMF中的算法。

.options 中可以设置是否保留每次的运算结果：keep.all=T。例：.options=list(keep.all=TRUE);如果设置parallel – p则可以让函数进行并行运算（例p1，p2）。如果在4个以内可以自行调内部函数解决：

res1 <- nmf(esGolub, 3, nrun=5,.opt='vp2', seed=123)

当大于4的时候就需要载入并行包doMPI。在这我们就不去细讲了。如果需要可以参考下提供的实例：

library(doMPI)
cl <- startMPIcluster()
registerDoMPI(cl)
library(NMF)
# run on all workers using the currentparallel backend
data(esGolub)
res <- nmf(esGolub, 3, 'brunet', nrun=n,.opt='p', .pbackend=NULL)
17http://cran.r-project.org/package=doMPI
17
# save result
save(res, file='result.RData')
## 4. Shutdown the cluster and quit MPI
closeCluster(cl)
mpi.quit()

接下来是结果的可视化展示：

评估结果的绘制：

首先构建需要评估的模型：

estim.r <- nmf(esGolub, 2:6, nrun=10,seed=123456)
plot(estim.r)

当然我们也可以只提取里面的任何一个数据然后绘制图像：

coph <- estim.r$measures$cophenetic
plot(2:6,coph, type="b",col="purple")

判断最佳rank值的准则就是，cophenetic值随K变化的最大变动的前点，如4-5变化最大，所以选择最佳rank值为4。

有时候我们可能需要对多个等级进行对比然后选择适合我们模型的等级，那么可视化将会帮助我们大忙。例：

V.random <- randomize(esGolub)
# estimate quality measures from theshuffled data (use default NMF algorithm)
estim.r.random <- nmf(V.random, 2:6,nrun=10, seed=123456)
# plot measures on same graph
plot(estim.r, estim.r.random)

热图的绘制：

consensusmap(estim.r, annCol=esGolub,labCol=NA, labRow=NA)

热图反映了对应的等级的选择最后所展示的聚类效果，从图中rank=4我们可以看出其聚类的情况最接近实际情况，当然这个图你是看不出来了，自己绘制下，得放大了看。

误差追踪主要是为了帮我们筛选适合我们的模型：

res <- nmf(esGolub, 3, .options='t')
 
res.multi.method <- nmf(esGolub, 3,list('brunet', 'lee', 'ns'), seed=123456, .options='t')

compare(res.multi.method)#对比各模型计算的参数：
 
compare(res.multi.method,$Cell)#当然也可以选择其中子数据进行评估。

接下来就是误差追踪：

plot(res)
plot(res.multi.method)

那么，我们的W，H在哪呢，我们最后还是需要我们的特征向量。如下参数可以获取：

res@fit@H

那接下来就可以去拿特征向量去做自己的研究了。

至此我们就可以筛选到属于我们的最优模型以及最优rank。那么我们总结下具体的一个流程：