【R数据处理】随机森林机器学习法
在生态学中随机森林机器学习法的应用有很多，其中一种用法是用来寻找生物标记（志）物（biomarkers）。什么叫生物标志物，简单来说生物标志物就是能用来表征整体生物群落特征的某一种或几种生物。对于单一生物体来说，生物标志物指生物体包含的能够标志该生物体某一类特征的性状。类似于植物的树冠大小、叶片大小、根长、株高等生物特征。
生物标志物有什么用处呢？举例来说，若研究者调查生长在不同海拔高度的植物性状如何随温度变化而改变的，其中一项比较重要的标志物就是叶片的大小。现在，我们很清楚叶片大小与温度的关系。因此，大叶片的植物往往意味着该植物生长在温暖的低海拔地区，而小叶片的植物往往意味着该植物生长在寒冷的高海拔地区。由此可见，生物标志物的确定对研究某些潜在的关系以及科学预测至关重要。（哇，我竟然写了这么多文字。。。）
计算分组对数据的解释量
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%以上。
寻找生物标志物
生物标志物的寻找原理其实比较简单，类似于回归分析中的逐步回归分析；也就是说当存在该物种时，判断模型的准确度是增加还是降低，从而找到一组物种的组合，使该组合满足数量相对较少，而模型的准确度又相对较高。
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])
画图
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提取生物标志物和数据输出
经过上述一统计算我们终于知道了多少生物标志物的组合是最优了，那么我们把这些生物标志物在所有（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)