RandomForest：随机森林预测生物标记biomarker——回归

关于随机森林的简介和应用理论，请阅读之前分享的文章：

关于随机森林进行分类的入门实战，请阅读之前分享的

大家可以学习此文，实现分组挖掘两组或多组的特异Features，也可以展示特征的贡献度，获得分类评估的准确度，以及使用新数据进行预测，无监督的随机森林等基础技能。

今天我们讲使用randomForest实现回归分析的实战代码。回归的应用主要包括时间序列预测模式，如预测股票、尸体死亡时间等。

本节不需要先难知识可也直接学习使用。

RandomForest安装与加载

# 安装 install.packages("randomForest") # 加载 library(randomForest)

回归Classification

先了解一下输入数据格式，方便准备

使用R内置按天记录的空气质量数据

data(airquality) head(airquality)

数据包括157天中，臭氧、太阳强度、风和温度，部分有缺失。前4列属性数据，后2列时间月和日为分组数据。

 Ozone Solar.R Wind Temp Month Day 1    41     190  7.4   67     5   1 2    36     118  8.0   72     5   2 3    12     149 12.6   74     5   3 4    18     313 11.5   62     5   4 5    NA      NA 14.3   56     5   5 6    28      NA 14.9   66     5   6

设置随机数种子保证结果可重复

set.seed(315)

随机森林回归臭氧与其它所有属性

ozone.rf= randomForest(Ozone ~ ., data=airquality, mtry=3,                     importance=TRUE, na.action=na.omit) print(ozone.rf)

结果如下：包括分析的命令，分析类型，树数量，重要的变量(Feature)个数，平均残差平方，解析率。

   Call: randomForest(formula = Ozone ~ ., data = airquality, mtry = 3,      importance = TRUE, na.action = na.omit)               Type of random forest: regression                     Number of trees: 500 No. of variables tried at each split: 3          Mean of squared residuals: 304.4269                    % Var explained: 72.26

查看每个变量的分类贡献度，显示仅保留两位小数可读性更好

round(importance(ozone.rf), 2)        %IncMSE IncNodePurity Solar.R   10.40      10833.08 Wind      23.55      43838.66 Temp      47.13      53731.95 Month      2.04       1504.72 Day        0.91       6306.42

结果为每个相关变量对应两列值。%IncMSE是Increased in mean squared error (%)，直译为增长的错误率平方均值，即去除该变量后，对目标预测的准确度下降的低，可理解为对目标变量预测准确的贡献度。IncNodePurity是Increased node purity，是另一种评估的方法。这里我们只关注%IncMSE就够了。

varImpPlot(ozone.rf)  

交叉验证cross-validation

# 先清空NA的样本，验证不允许有NA airquality = na.omit(airquality) myairquality= cbind(airquality[1:6], matrix(runif(96 * nrow(airquality)), nrow(airquality), 96)) # 交驻验证添加了随机数的训练集，分组，交叉验证的次数 result= rfcv(myairquality, airquality$Ozone, cv.fold=3) # 绘制错误率曲线，观查错误率与使用Markers数量的变化 with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))

我们看到一个现象，不是feature越多越好，无关的feature如果多了，反而错误率上升，会影响预测的准确度。

# 使用replicate进行多次交叉验证，可选 result= replicate(5, rfcv(myairquality, airquality$Ozone), simplify=FALSE) error.cv= sapply(result, "[[", "error.cv") matplot(result[[1]]$n.var, cbind(rowMeans(error.cv), error.cv), type="l",        lwd=c(2, rep(1, ncol(error.cv))), col=1, lty=1, log="x",        xlab="Number of variables", ylab="CV Error")

多次验证结果类型，更能说明结果的可信度。

严谨总没有坏处，好的结果都是多角度证明的。

大家学习随机森林的分类、和回归。将来有时间，将带大家上手重复一些高水平文章中的分析，结合具体生物学问题会更有意思。

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：生信宝典 学术图表 高分文章 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外150+ PI，1500+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读