查看: 937|回复: 3

[技术前沿] 如何构建随机森林及筛选重要特征变量?

[复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
907
注册时间
2020.6.16
在线时间
485 小时

发表于 2022.6.24 10:54:43 | 显示全部楼层 |阅读模式
基于TCGA临床数据整理、生存分析、生存模型构建和模型预测等内容可戳往期:
《一文拿捏住TCGA临床和基因表达量数据处理!》
《TCGA更新了?!最新数据处理拿捏住!》
《重要的生存分析前,你必须做好这个准备!》
《真香!被自己做的生存分析美丽住了!》
《真香!被自己炫酷的批量生存分析操作厉害住了!》
《嘿!高端的单因素Cox回归,往往只需要朴素的烹饪!!》
《如何使用lasso回归筛选重要特征变量?》
《如何进行数据切割并使用ROC评估lasso建模好坏?》
《如何进行多因素cox回归分析及可视化?》
《如何构建最优多因素cox回归模型?》

除了使用cox或lasso回归进行生存模型构建和预测,今天我们继续学习一种基于决策树的集成学习算法——随机森林(random forest),它使用引导聚合(Bagging)和预测变量随机化来实现高度的预测准确性。


简单来说,随机森林通过自主抽样(bootstrap),从原数据集中有放回地抽取多个样本,对抽取的样本用决策树进行训练,生成多棵决策树。每棵决策树都是一个分类或预测器,只要我们输入一个样本,n棵决策树就会有n个分类或预测的结果。随机森林集成了所有决策树的投票,将投票次数最多的结果指定为最终输出。

(图片源自网络)

每棵决策树使用的样本和变量特征各不相同,因为我们在训练之初并不知道哪些特征对训练的结果影响更大,随机性降低了对结果的影响。

我们同样可以使用随机森林模型来筛选重要的特征变量。

1. 将数据切割为训练集和测试集

  1. #randomforest:
  2. #加载TCGA-KIRC表达矩阵和临床信息(往期整理):
  3. load("KIRC_exp_cl.Rdata")
  4. View(exp2)
复制代码


  1. View(cl)
复制代码


  1. #加载所需R包:
  2. library(randomForest)#随机森林建模
  3. library(caret)#切割数据
  4. library(ggplot2)
  5. library(pROC)#ROC曲线

  6. #将总数据集切割为训练集(70%)和测试集(30%):
  7. set.seed(233)
  8. inA <- createDataPartition(cl$event, p = 0.7,times = 1,list = F)
  9. #训练集:
  10. train_exp <- exp2[,inA]
  11. train_cl <- cl[inA,]
  12. #测试集:
  13. test_exp <- exp2[,-inA]
  14. test_cl <- cl[-inA,]

  15. #切割比例:
  16. prop.table(table(train_cl$event))
  17. prop.table(table(test_cl$event))
复制代码


2. 训练集构建随机森林

  1. #模型拟合(会拟合较长时间,需耐心等待):
  2. x = t(train_exp) #转置
  3. y = train_cl$event
  4. train_forest <- randomForest(x = x,
  5. y = y,
  6. importance = TRUE, #计算每个特征的重要性
  7. ntree = 500, #决策树棵数(默认500)
  8. mtry = if (!is.null(y) && !is.factor(y))
  9. max(floor(ncol(x)/3), 1) else floor(sqrt(ncol(x))), #mytry指定节点中用于二叉树的变量个数,默认为数据集变量个数的二次方根(分类模型)或三分之一(预测模型),也可以人为逐次挑选最佳变量数;
  10. proximity = TRUE )
  11. train_forest #模型概要
复制代码


  1. plot(train_forest)
复制代码


可以看到在100棵决策树后error基本趋于稳定。

  1. #模型预测:训练集模型预测自身:
  2. train_predict <- predict(train_forest, x)
  3. compare_train <- as.data.frame(cbind(y ,train_predict))
  4. colnames(compare_train) <- c("event","train_predict")
  5. head(compare_train)
复制代码


  1. #箱线图展示预测结果:
  2. compare_train$event <- as.factor(compare_train$event)
  3. p <- ggplot(compare_train,
  4. aes(x = event,y = train_predict,group = event))+
  5. geom_boxplot(compare_train,
  6. mapping = aes(color = event,fill = event,group = event),
  7. alpha = 0.4,
  8. show.legend = F)+
  9. geom_jitter(aes(fill = event,color = event),
  10. alpha = 0.3,width = 0.2,shape = 21,size = 1.5)+
  11. scale_y_continuous(limits = c(0,1))+
  12. theme_bw()
  13. p
复制代码


肉眼可见预测结果的准确性非常高。我们再用ROC曲线判断一下:

  1. #ROC曲线:
  2. roc1 <- roc(event ~ train_predict, compare_train,
  3. direction =
  4. "<")
  5. print(roc1)
复制代码


  1. plot.roc(roc1,
  2. print.auc=TRUE,
  3. auc.polygon=T,
  4. auc.polygon.col="skyblue")
复制代码


训练集建模预测训练集数据,AUC高达1,完美分类器!下面我们筛选影响模型结果的重要基因(变量/特征):

3. 变量重要性

  1. #展示重要基因及评价指标:
  2. varImpPlot(train_forest,
  3. n.var=30, #显示解释模型top30的重要基因
  4. scale=FALSE, #是否将度量划分为标准误
  5. main="Top 30 genes")
复制代码


两张图,分别为影响模型结果的两个重要评价指标,各显示了top30;

%IncMSE(percentage of increase of mean square error):通过对每一个预测变量随机赋值,如果该预测变量更为重要,那么其值被随机替换后模型预测的误差会增大;该值越大表示该变量的重要性越大;

IncNodePurity(Increase in Node Purity):通过残差平方和来度量每个变量对分类树每个节点上观测值的异质性的影响,该值越大表示该变量的重要性越大。

  1. #提取重要基因/变量:
  2. importance <-importance(train_forest) #提取基因和对应重要性评价指标取值
  3. head(importance)#此时还未进行排序;
复制代码


  1. importance_order <- importance[order(importance[,2],decreasing = TRUE),]#按照IncNodePurity值从高到低排序;
  2. head(importance_order)
复制代码


  1. randomforest_genes <- rownames(importance_order[1:30,])#提取前30的基因;
  2. randomforest_genes
复制代码


4. 测试集数据预测

  1. #训练集模型预测测试集:
  2. x = t(test_exp)
  3. y = test_cl$event

  4. test_predict <- predict(train_forest, x)
  5. compare_test <- as.data.frame(cbind(y ,test_predict))
  6. colnames(compare_test) <- c("event","test_predict")
  7. head(compare_test)
复制代码


  1. #箱线图展示预测结果:
  2. compare_test$event <- as.factor(compare_test$event)
  3. p1 <- ggplot(compare_test,
  4. aes(x = event,y = test_predict,group = event))+
  5. geom_boxplot(compare_test,
  6. mapping = aes(color = event,fill = event,group = event),
  7. alpha = 0.4,
  8. show.legend = F)+
  9. geom_jitter(aes(fill = event,color = event),
  10. alpha = 0.3,width = 0.2,shape = 21,size = 1.5)+
  11. scale_y_continuous(limits = c(0,1))+
  12. theme_bw()
  13. p1
复制代码


用来预测新数据准确度肉眼变低。

  1. #ROC曲线:
  2. roc2 <- roc(event ~ test_predict, compare_test,
  3. direction = "<")
  4. print(roc2)
复制代码


  1. plot.roc(roc2,
  2. print.auc=TRUE,
  3. auc.polygon=T,
  4. auc.polygon.col="skyblue")
复制代码


AUC值还是有0.738,可见模型可是不错的。

好啦,今天的分享就到这里!

参考资料:https://www.stat.berkeley.edu/~breiman/RandomForests/

本文作者:基迪奥-喵酱

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2016.10.5
在线时间
47 小时

发表于 2022.6.24 14:22:44 | 显示全部楼层
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
63 小时

发表于 2022.6.25 19:18:10 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2021.12.31
在线时间
6 小时

发表于 2022.11.24 14:10:28 | 显示全部楼层
云平台上怎么找不到随机森林分析这个工具了呀
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表