【R数据处理】随机森林机器学习法

 知其然,也要知其所以然。”   --科白君
"R数据分析"专题·第14篇
  编辑 | free傻孩子
  5021字 |13分钟阅读
本期推送内容
机器学习研究和构建的是一种特殊算法(而非某一个特定的算法),能够让计算机自己在数据中学习从而进行预测。所以,机器学习不是某种具体的算法,而是很多算法的统称。机器学习包含了很多种不同的算法,深度学习就是其中之一,其他方法包括决策树,聚类,贝叶斯等。随机森林机器学习法就是利用决策树对数据进行分析的一种方法。

在生态学中随机森林机器学习法的应用有很多,其中一种用法是用来寻找生物标记(志)物(biomarkers)。什么叫生物标志物,简单来说生物标志物就是能用来表征整体生物群落特征的某一种或几种生物。对于单一生物体来说,生物标志物指生物体包含的能够标志该生物体某一类特征的性状。类似于植物的树冠大小、叶片大小、根长、株高等生物特征。

生物标志物有什么用处呢?举例来说,若研究者调查生长在不同海拔高度的植物性状如何随温度变化而改变的,其中一项比较重要的标志物就是叶片的大小。现在,我们很清楚叶片大小与温度的关系。因此,大叶片的植物往往意味着该植物生长在温暖的低海拔地区,而小叶片的植物往往意味着该植物生长在寒冷的高海拔地区。由此可见,生物标志物的确定对研究某些潜在的关系以及科学预测至关重要。(哇,我竟然写了这么多文字。。。)

01

计算分组对数据的解释量

1)建立分组

在随机森林分析中,分组必须为因子变量,因此在计算之前需把分组转变为因子变量。

首先,加载需要用到的包和数据,代码如下:

library(tidyverse)
library(randomForest)

load("testData.RData")
head(testData[,1:6])#数据格式如下

注意:数据为数据框格式,包含行名和列名。

testData2 <- testData %>%
mutate(treatments = as.factor(treatments))#把分组转变为因子变量

2)计算分组预测out分布的准确性

在本测试数据中分组有3种处理,分别是in_situ, upward和downward。为了检测该分组对OTU水平数据分布预测的准确性,我们使用了随机森林分析。代码如下:

set.seed(123)
otu_rf= randomForest(treatments ~ ., data = testData2, importance=TRUE, proximity=TRUE)
otu_rf

我们的结果表明,在OTU水平,分组对数据分布预测的准确性为69.57%(100 - OOB)。因为我们的数据是测试数据,分组是相对随便的,所以预测的准确性不太高,如果是真实数据,准确性一般在90%以上。

02

寻找生物标志物

生物标志物的寻找原理其实比较简单,类似于回归分析中的逐步回归分析;也就是说当存在该物种时,判断模型的准确度是增加还是降低,从而找到一组物种的组合,使该组合满足数量相对较少,而模型的准确度又相对较高。

1)计算otu对模型的重要性

imp_otu <- as_tibble(round(importance(otu_rf), 2),rownames = "OTUid") %>%
  arrange(desc(MeanDecreaseAccuracy))
imp_otu

结果已经按照“MeanDecreaseAccuracy”进行了排序(降序)。由结果表明,当otu655不存在时,模型的准确度降低了5.06%,相似的当otu657不存在时模型的准确度降低了4.87%,这意味着这几个out对模型很重要。

2)5次重复的10倍交叉验证

ncol(testData2)#列的数量
myotu= testData2[,-1]#删除因子(第几列为因子这里的数字就是几)
set.seed(123)

result<-rfcv(myotu, testData2$treatments, cv.fold=5, scale = "log", step = 0.9)
result1<-result
result1$n.var # 723

with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))#绘图,查看图形走势

结果表明,在n(选取的otu数量)为22时,误差CV最低,随着数量的增加误差CV先增加后保持平稳,因此,对于本测试数据来说,最优的生物标志物数量为22。

3)再次对该结果进行验证

#5次重复计算

#5次重复计算
for (i in 124:127){
  print(i)
  set.seed(i)
  result= rfcv(myotu, testData2$treatments, cv.fold=5, scale = "log", step = 0.9)
  error.cv = cbind(error.cv, result$error.cv)
}

n.var <- error.cv$num
error.cv <- error.cv[,-1]
colnames(error.cv)<- paste('err',1:5,sep='.')#添加列名
err.mean <-apply(error.cv,1,mean)#求5列的平均值
allerr<-data.frame(num=n.var,err.mean=err.mean,error.cv)
head(allerr[,1:6])

03

画图

optimal = 22#上述计算获得的最佳生物标志物的数量
main_theme = theme(panel.background=element_blank(),
panel.grid=element_blank(),
axis.line.x=element_line(size=0.5, colour="black"),
axis.line.y=element_line(size=0.5, colour="black"),
axis.ticks=element_line(color="black"),
axis.text=element_text(color="black", size=12),
legend.position="right",
legend.background=element_blank(),
legend.key=element_blank(),
legend.text= element_text(size=12),
text=element_text(family="sans", size=12))
max(allerr$num)#查看要绘图的坐标轴最大值
ggplot() +
  geom_line(aes(x = allerr$num, y = allerr$err.1), colour = 'grey',size=0.5) +
  geom_line(aes(x = allerr$num, y = allerr$err.2), colour = 'grey',size=0.5) +
  geom_line(aes(x = allerr$num, y = allerr$err.3), colour = 'grey',size=0.5) +
  geom_line(aes(x = allerr$num, y = allerr$err.4), colour = 'grey',size=0.5) +
  geom_line(aes(x = allerr$num, y = allerr$err.5), colour = 'grey',size=0.5) +
  geom_line(aes(x = allerr$num, y = allerr$err.mean), colour = 'black',size=0.5) +
  geom_vline(xintercept = optimal, colour='black', lwd=0.36, linetype=2) +
  coord_trans(x = "log2") +
  scale_x_continuous(breaks = c(1, 5, 10, 20, 30, 50, 100, 300, 723)) + # , max(allerr$num)
  labs( x='Number of OTUs ', y='Cross-validation error rate') +
  annotate("text", x = optimal, y = max(allerr$err.mean), label=paste("Optimal = ", optimal, sep=""))+
  main_theme

由图,我们发现当选择22作为生物标志物的数量时,似乎5次重复的平均误差值并不是最低的,这时我们可以调整要选择的生物标志物的数量以使平均误差值最低。在本图中,我们可以看到最低值出现在30-50之间,根据上述我们计算的结果如下:

我们看到在30-50之间有5个数,分别是31,34,38,42和47,分别验证后我们发现Optimal为38时平均误差最低,结果如下:

optimal = 38
ggplot() +
geom_line(aes(x = allerr$num, y = allerr$err.1), colour = 'grey',size=0.5) +
geom_line(aes(x = allerr$num, y = allerr$err.2), colour = 'grey',size=0.5) +
geom_line(aes(x = allerr$num, y = allerr$err.3), colour = 'grey',size=0.5) +
geom_line(aes(x = allerr$num, y = allerr$err.4), colour = 'grey',size=0.5) +
geom_line(aes(x = allerr$num, y = allerr$err.5), colour = 'grey',size=0.5) +
geom_line(aes(x = allerr$num, y = allerr$err.mean), colour = 'black',size=0.5) +
geom_vline(xintercept = optimal, colour='black', lwd=0.36, linetype=2) +
coord_trans(x = "log2") +
scale_x_continuous(breaks = c(1, 5, 10, 20, 30, 50, 100, 300, 723)) + # , max(allerr$num)
labs( x='Number of OTUs ', y='Cross-validation error rate') +
annotate("text", x = optimal, y = max(allerr$err.mean), label=paste("Optimal = ", optimal, sep=""))+
main_theme ->p1
p1
04

提取生物标志物和数据输出

经过上述一统计算我们终于知道了多少生物标志物的组合是最优了,那么我们把这些生物标志物在所有(723个OTUs)中拿出来吧。

#抽取最重要的前38个otu

write.csv(imp_otu[1:38,], "biomarks.csv")
imp_otu[1:38,] %>%
select(OTUid,MeanDecreaseAccuracy) %>%
arrange(MeanDecreaseAccuracy) %>%
mutate(OTUid = forcats::fct_inorder(OTUid)) %>%
ggplot(aes(x = OTUid, y = MeanDecreaseAccuracy))+
geom_bar(stat = "identity")+
labs(x = "", y = "Mean decrease accuracy")+
coord_flip()+
main_theme ->p2
p2

组图

library(customLayout)
library(magrittr)
library(export)

lay1 = lay_new(mat = matrix(1, ncol = 1))
lay2 =lay_bind_col(lay1, lay1, widths = c(1, 1))

lay_show(lay2)
lay_grid(list(p1,p2), lay2)
graph2ppt(file = "random.pptx",width = 12,height = 6, append =FALSE)

(0)

相关推荐

  • ggClusterNet:手写算法在R语言中实现Gephi布局之一

    写在前面 像不像Gephi出图,Gephi布局还是有很多人喜欢,R语言确实不好实现,这个算法需要大量的迭代,让模块化点之间连线最短,不断修正,很耗费时间,我将聚类算法和修改后的距离算法合并起来并结合节 ...

  • 微生信生物新年放大招:一条代码完成堆叠柱状图-冲击图的操作-终结版

    00. 写在前面 眼下2019年的余额不足一天了,2020年最终还是要到来,首先祝愿大家元旦快乐,其次新的一年祝大家心想事成.说实在的最近老是想出去玩,哈哈哈,文章还没写完,不可以!不可以! 这只是大 ...

  • RDA_环境因子_群落结构_统计检验_可视化

    RDA环境因子群落结构统计检验可视化 环境因子的筛选及数据的转化方面请参阅宏基因组公众号之前的推文,本文主要侧重统计分析与可视化 看到师兄文章里的图自己可能用到,想复现一下,于是就尝试了一下,顺便写个 ...

  • Ven图-Upset-花瓣图-三个NO-二分网络-YES

    写在前面 让点随机排布在一个区域,保证点之间不重叠,并且将点的图层放到最上层,保证节点最清晰,然后边可以进行透明化,更加突出节点的位置.这里我新构建了布局函数PolyRdmNotdCirG 来做这个随 ...

  • 对邓老师组开发不同域之间二分网络的思考和ggCLusterNet的进一步完善

    不同域之间的相互关系 例如:细菌和真菌,不考虑内部相关,只考虑细菌和真菌之间的相关. 邓老师在前些日子的讲座这提到的这种网络,称之为"二分网络",其实二分网络是我17年听到的一个词 ...

  • 在ggClusterNet中仿造cytosccape添加多行的聚类布局

    写在前面 算法肯定不同,功能相近. 时间戳:现在2020年8月8日,我在上周便构造了这个函数,用于模仿cytoscape网络的矩阵布局.类似下图的样式. 但是我觉得这个算法我写的不够完善,后面应该还会 ...

  • 排序出图三个重复带来的一系列问题的解决方案

    wentao 2021/8/8 写在前面 目前patchwork包拼图很不错,我尝试将群落两两比较结果和排序图放到一起看看: 三个样本的排序图置信椭圆绘制这个问题,我来做两个方案: 三个重复的置换检验 ...

  • ggClusterNet-多边形模块化排布网络(PolygonClusterG)

    写在前面 专辑ggClusterNet包今日开始连载,虽然本包中函数不多,但都是精巧实用之辈.相信可以帮助到大家. devtools::install_github("taowenmicro ...

  • 16s分析之进化树+差异分析(一)

    之前我便说过,16s出图主要是用R,自这两年来Y叔的ggtree出来,好评不断,引用不断攀升 这里当然借用Y树的包来做我们的进化树,当然Y叔公众号讲的更加深入和详细,有需要请关注:biobabble ...

  • ggClusterNet-多边形可变变径布局网络(randomRrClusterG)

    写在前面 多边形排布,将不同的模块按照多边形排布进行.模块越大的部分,其半径越大,是按照模块内节点数量判断到多边形中心的距离.说通俗一点,就是这个函数可以根据模块节点多少设置合适的大小展示,并且不会重 ...

  • 把玫瑰图大胆用起来

    写在前面 玫瑰图虽然在疫情数据上有了应用,但是在我们微生物生态领域可以完好的使用,除非你解决一下两个问题: 1 如何让这个圈比较小的部分可以清楚的被标记和识别 2 我不想要一整个圈,来半个圈行不行. ...

  • 人工设置网络模块位置2-大尺度空间网络布局应用

    写在前面 多变形排布,将不同的模块按照多变形排布进行.模块越大的部分,其半径越大,是按照模块内节点数量判断到多变形中心的距离. 微生物网络 输入文件 #--导入所需R包#------- library ...

  • 随机森林和特征可视化-(火柴图和玫瑰图)

    写在前面2 今天补充几句话,这部分机器学习的完整内容也会整理到刘老师主编的.实验和分析方法的项目中,这里只是一个引子,供大家理解和学习. 写在前面1 现在机器学习,尤其是随机森林在微生物组学的应用越来 ...