一套完整的基于随机森林的机器学习流程(特征选择、交叉验证、模型评估))...
机器学习实操(以随机森林为例)
为了展示随机森林的操作,我们用一套早期的前列腺癌和癌旁基因表达芯片数据集,包含102
个样品(50个正常,52个肿瘤),2
个分组和9021
个变量 (基因)。(https://file.biolab.si/biolab/supp/bi-cancer/projections/info/prostata.html)
数据格式和读入数据
输入数据为标准化之后的表达矩阵,基因在行,样本在列。随机森林对数值分布没有假设。每个基因表达值用于分类时是基因内部在不同样品直接比较,只要是样品之间标准化的数据即可,其他任何线性转换如log2
,scale
等都没有影响 (数据在:https://gitee.com/ct5869/shengxin-baodian/tree/master/machinelearning)。
样品表达数据:
prostat.expr.txt
样品分组信息:
prostat.metadata.txt
expr_file <- "data/prostat.expr.symbol.txt"
metadata_file <- "data/prostat.metadata.txt"# 每个基因表达值是内部比较,只要是样品之间标准化的数据即可,其它什么转换都关系不大
# 机器学习时,字符串还是默认为因子类型的好
expr_mat <- read.table(expr_file, row.names = 1, header = T, sep="\t", stringsAsFactors =T)# 处理异常的基因名字
rownames(expr_mat) <- make.names(rownames(expr_mat))metadata <- read.table(metadata_file, row.names=1, header=T, sep="\t", stringsAsFactors =T)dim(expr_mat)## [1] 9021 102
基因表达表示例如下:
expr_mat[1:4,1:5]## normal_1 normal_2 normal_3 normal_4 normal_5
## AADAC 1.3 -1 -7 -4 5
## AAK1 0.4 0 10 11 8
## AAMP -0.4 20 -7 -14 12
## AANAT 143.3 19 397 245 328
Metadata表示例如下
head(metadata)## class
## normal_1 normal
## normal_2 normal
## normal_3 normal
## normal_4 normal
## normal_5 normal
## normal_6 normaltable(metadata)## metadata
## normal tumor
## 50 52
样品筛选和排序
对读入的表达数据进行转置。通常我们是一行一个基因,一列一个样品。在构建模型时,数据通常是反过来的,一列一个基因,一行一个样品。每一列代表一个变量 (variable
),每一行代表一个案例 (case
)。这样更方便提取每个变量,且易于把模型中的x,y
放到一个矩阵中。
样本表和表达表中的样本顺序对齐一致也是需要确保的一个操作。
# 表达数据转置
# 习惯上我们是一行一个基因,一列一个样品
# 做机器学习时,大部分数据都是反过来的,一列一个基因,一行一个样品
# 每一列代表一个变量
expr_mat <- t(expr_mat)
expr_mat_sampleL <- rownames(expr_mat)
metadata_sampleL <- rownames(metadata)common_sampleL <- intersect(expr_mat_sampleL, metadata_sampleL)# 保证表达表样品与METAdata样品顺序和数目完全一致
expr_mat <- expr_mat[common_sampleL,,drop=F]
metadata <- metadata[common_sampleL,,drop=F]
判断是分类还是回归
前面读数据时已经给定了参数stringsAsFactors =T
,这一步可以忽略了。
如果group对应的列为数字,转换为数值型 - 做回归
如果group对应的列为分组,转换为因子型 - 做分类
# R4.0之后默认读入的不是factor,需要做一个转换
# devtools::install_github("Tong-Chen/ImageGP")
library(ImageGP)# 此处的class根据需要修改
group = "class"# 如果group对应的列为数字,转换为数值型 - 做回归
# 如果group对应的列为分组,转换为因子型 - 做分类
if(numCheck(metadata[[group]])){if (!is.numeric(metadata[[group]])) {metadata[[group]] <- mixedToFloat(metadata[[group]])}
} else{metadata[[group]] <- as.factor(metadata[[group]])
}
随机森林一般分析
library(randomForest)# 查看参数是个好习惯
# 有了前面的基础概述,再看每个参数的含义就明确了很多
# 也知道该怎么调了
# 每个人要解决的问题不同,通常不是别人用什么参数,自己就跟着用什么参数
# 尤其是到下游分析时
# ?randomForest# 查看源码
# randomForest:::randomForest.default
加载包之后,直接分析一下,看到结果再调参。
# 设置随机数种子,具体含义见 https://mp.weixin.qq.com/s/6plxo-E8qCdlzCgN8E90zg
set.seed(304)# 直接使用默认参数
rf <- randomForest(expr_mat, metadata[[group]])
查看下初步结果, 随机森林类型判断为分类
,构建了500
棵树,每次决策时从随机选择的94
个基因中做最优决策 (mtry
),OOB
估计的错误率是9.8%
,挺高的。
分类效果评估矩阵Confusion matrix
,显示normal
组的分类错误率为0.06
,tumor
组的分类错误率为0.13
。
rf##
## Call:
## randomForest(x = expr_mat, y = metadata[[group]])
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 94
##
## OOB estimate of error rate: 9.8%
## Confusion matrix:
## normal tumor class.error
## normal 47 3 0.0600000
## tumor 7 45 0.1346154
随机森林标准操作流程 (适用于其他机器学习模型)
拆分训练集和测试集
library(caret)
seed <- 1
set.seed(seed)
train_index <- createDataPartition(metadata[[group]], p=0.75, list=F)
train_data <- expr_mat[train_index,]
train_data_group <- metadata[[group]][train_index]test_data <- expr_mat[-train_index,]
test_data_group <- metadata[[group]][-train_index]dim(train_data)## [1] 77 9021dim(test_data)## [1] 25 9021
Boruta特征选择鉴定关键分类变量
# install.packages("Boruta")
library(Boruta)
set.seed(1)boruta <- Boruta(x=train_data, y=train_data_group, pValue=0.01, mcAdj=T, maxRuns=300)boruta## Boruta performed 299 iterations in 1.937513 mins.
## 46 attributes confirmed important: ADH5, AGR2, AKR1B1, ANGPT1,
## ANXA2.....ANXA2P1.....ANXA2P3 and 41 more;
## 8943 attributes confirmed unimportant: AADAC, AAK1, AAMP, AANAT, AARS
## and 8938 more;
## 32 tentative attributes left: ALDH2, ATP6V1G1, C16orf45, CDC42BPA,
## COL4A6 and 27 more;
查看下变量重要性鉴定结果(实际上面的输出中也已经有体现了),54
个重要的变量,36
个可能重要的变量 (tentative variable
, 重要性得分与最好的影子变量得分无统计差异),6,980
个不重要的变量。
table(boruta$finalDecision)##
## Tentative Confirmed Rejected
## 32 46 8943
绘制鉴定出的变量的重要性。变量少了可以用默认绘图,变量多时绘制的图看不清,需要自己整理数据绘图。
定义一个函数提取每个变量对应的重要性值。
library(dplyr)
boruta.imp <- function(x){imp <- reshape2::melt(x$ImpHistory, na.rm=T)[,-1]colnames(imp) <- c("Variable","Importance")imp <- imp[is.finite(imp$Importance),]variableGrp <- data.frame(Variable=names(x$finalDecision), finalDecision=x$finalDecision)showGrp <- data.frame(Variable=c("shadowMax", "shadowMean", "shadowMin"),finalDecision=c("shadowMax", "shadowMean", "shadowMin"))variableGrp <- rbind(variableGrp, showGrp)boruta.variable.imp <- merge(imp, variableGrp, all.x=T)sortedVariable <- boruta.variable.imp %>% group_by(Variable) %>% summarise(median=median(Importance)) %>% arrange(median)sortedVariable <- as.vector(sortedVariable$Variable)boruta.variable.imp$Variable <- factor(boruta.variable.imp$Variable, levels=sortedVariable)invisible(boruta.variable.imp)
}boruta.variable.imp <- boruta.imp(boruta)head(boruta.variable.imp)## Variable Importance finalDecision
## 1 AADAC 0 Rejected
## 2 AADAC 0 Rejected
## 3 AADAC 0 Rejected
## 4 AADAC 0 Rejected
## 5 AADAC 0 Rejected
## 6 AADAC 0 Rejected
只绘制Confirmed
变量。
library(ImageGP)sp_boxplot(boruta.variable.imp, melted=T, xvariable = "Variable", yvariable = "Importance",legend_variable = "finalDecision", legend_variable_order = c("shadowMax", "shadowMean", "shadowMin", "Confirmed"),xtics_angle = 90)
提取重要的变量和可能重要的变量
boruta.finalVarsWithTentative <- data.frame(Item=getSelectedAttributes(boruta, withTentative = T), Type="Boruta_with_tentative")
看下这些变量的值的分布
caret::featurePlot(train_data[,boruta.finalVarsWithTentative$Item], train_data_group, plot="box")
交叉验证选择参数并拟合模型
定义一个函数生成一些列用来测试的mtry
(一系列不大于总变量数的数值)。
generateTestVariableSet <- function(num_toal_variable){max_power <- ceiling(log10(num_toal_variable))tmp_subset <- c(unlist(sapply(1:max_power, function(x) (1:10)^x, simplify = F)), ceiling(max_power/3))#return(tmp_subset)base::unique(sort(tmp_subset[tmp_subset<num_toal_variable]))
}
# generateTestVariableSet(78)
选择关键特征变量相关的数据
# 提取训练集的特征变量子集
boruta_train_data <- train_data[, boruta.finalVarsWithTentative$Item]
boruta_mtry <- generateTestVariableSet(ncol(boruta_train_data))
使用 Caret 进行调参和建模
library(caret)
# Create model with default parameters
trControl <- trainControl(method="repeatedcv", number=10, repeats=5)seed <- 1
set.seed(seed)
# 根据经验或感觉设置一些待查询的参数和参数值
tuneGrid <- expand.grid(mtry=boruta_mtry)borutaConfirmed_rf_default <- train(x=boruta_train_data, y=train_data_group, method="rf", tuneGrid = tuneGrid, # metric="Accuracy", #metric='Kappa'trControl=trControl)
borutaConfirmed_rf_default## Random Forest
##
## 77 samples
## 78 predictors
## 2 classes: 'normal', 'tumor'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 71, 69, 69, 69, 69, 69, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 1 0.9352381 0.8708771
## 2 0.9352381 0.8708771
## 3 0.9352381 0.8708771
## 4 0.9377381 0.8758771
## 5 0.9377381 0.8758771
## 6 0.9402381 0.8808771
## 7 0.9402381 0.8808771
## 8 0.9452381 0.8908771
## 9 0.9402381 0.8808771
## 10 0.9452381 0.8908771
## 16 0.9452381 0.8908771
## 25 0.9477381 0.8958771
## 36 0.9452381 0.8908771
## 49 0.9402381 0.8808771
## 64 0.9327381 0.8658771
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 25.
可视化不同参数的准确性分布
plot(borutaConfirmed_rf_default)
可视化Top20重要的变量
dotPlot(varImp(borutaConfirmed_rf_default))
提取最终选择的模型,并绘制 ROC 曲线评估模型
borutaConfirmed_rf_default_finalmodel <- borutaConfirmed_rf_default$finalModel
先自评,评估模型对训练集的分类效果
采用训练数据集评估构建的模型,Accuracy=1; Kappa=1
,非常完美。
模型的预测显著性P-Value [Acc > NIR] : 2.2e-16
。其中NIR
是No Information Rate
,其计算方式为数据集中最大的类包含的数据占总数据集的比例。如某套数据中,分组A
有80
个样品,分组B
有20
个样品,我们只要猜A
,正确率就会有80%
,这就是NIR
。如果基于这套数据构建的模型准确率也是80%
,那么这个看上去准确率较高的模型也没有意义。confusionMatrix
使用binom.test
函数检验模型的准确性Accuracy
是否显著优于NIR
,若P-value<0.05
,则表示模型预测准确率显著高于随便猜测。
# 获得模型结果评估矩阵(`confusion matrix`)predictions_train <- predict(borutaConfirmed_rf_default_finalmodel, newdata=train_data)
confusionMatrix(predictions_train, train_data_group)## Confusion Matrix and Statistics
##
## Reference
## Prediction normal tumor
## normal 38 0
## tumor 0 39
##
## Accuracy : 1
## 95% CI : (0.9532, 1)
## No Information Rate : 0.5065
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 1
##
## Mcnemar's Test P-Value : NA
##
## Sensitivity : 1.0000
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 1.0000
## Prevalence : 0.4935
## Detection Rate : 0.4935
## Detection Prevalence : 0.4935
## Balanced Accuracy : 1.0000
##
## 'Positive' Class : normal
##
盲评,评估模型应用于测试集时的效果
绘制ROC
曲线,计算模型整体的AUC
值,并选择最佳模型。
# 绘制ROC曲线prediction_prob <- predict(borutaConfirmed_rf_default_finalmodel, newdata=test_data, type="prob")
library(pROC)
roc_curve <- roc(test_data_group, prediction_prob[,1])roc_curve##
## Call:
## roc.default(response = test_data_group, predictor = prediction_prob[, 1])
##
## Data: prediction_prob[, 1] in 12 controls (test_data_group normal) > 13 cases (test_data_group tumor).
## Area under the curve: 0.9872# roc <- roc(test_data_group, factor(predictions, ordered=T))
# plot(roc)
基于默认阈值的盲评
基于默认阈值绘制混淆矩阵并评估模型预测准确度显著性,结果显著P-Value [Acc > NIR]<0.05
。
# 获得模型结果评估矩阵(`confusion matrix`)predictions <- predict(borutaConfirmed_rf_default_finalmodel, newdata=test_data)
confusionMatrix(predictions, test_data_group)## Confusion Matrix and Statistics
##
## Reference
## Prediction normal tumor
## normal 12 2
## tumor 0 11
##
## Accuracy : 0.92
## 95% CI : (0.7397, 0.9902)
## No Information Rate : 0.52
## P-Value [Acc > NIR] : 2.222e-05
##
## Kappa : 0.8408
##
## Mcnemar's Test P-Value : 0.4795
##
## Sensitivity : 1.0000
## Specificity : 0.8462
## Pos Pred Value : 0.8571
## Neg Pred Value : 1.0000
## Prevalence : 0.4800
## Detection Rate : 0.4800
## Detection Prevalence : 0.5600
## Balanced Accuracy : 0.9231
##
## 'Positive' Class : normal
##
选择模型分类最佳阈值再盲评
youden:max(sensitivities + r × specificities)
closest.topleft:min((1 − sensitivities)2 + r × (1 − specificities)2)
r
是加权系数,默认是1
,其计算方式为r = (1 − prevalenc**e)/(cos**t * prevalenc**e).
best.weights
控制加权方式:(cost
, prevalence
)默认是(1
,0.5
),据此算出的r
为1
。
cost: 假阴性率占假阳性率的比例,容忍更高的假阳性率还是假阴性率
prevalence: 关注的类中的个体所占的比例 (
n.cases/(n.controls+n.cases)
).
best_thresh <- data.frame(coords(roc=roc_curve, x = "best", input="threshold", transpose = F, best.method = "youden"))best_thresh$best <- apply(best_thresh, 1, function (x) paste0('threshold: ', x[1], ' (', round(1-x[2],3), ", ", round(x[3],3), ")"))best_thresh## threshold specificity sensitivity best
## 1 0.672 0.9166667 1 threshold: 0.672 (0.083, 1)
准备数据绘制ROC曲线
library(ggrepel)
ROC_data <- data.frame(FPR = 1- roc_curve$specificities, TPR=roc_curve$sensitivities)
ROC_data <- ROC_data[with(ROC_data, order(FPR,TPR)),]p <- ggplot(data=ROC_data, mapping=aes(x=FPR, y=TPR)) +geom_step(color="red", size=1, direction = "vh") +geom_segment(aes(x=0, xend=1, y=0, yend=1)) + theme_classic() + xlab("False positive rate") + ylab("True positive rate") + coord_fixed(1) + xlim(0,1) + ylim(0,1) +annotate('text', x=0.5, y=0.25, label=paste('AUC=', round(roc_curve$auc,2))) +geom_point(data=best_thresh, mapping=aes(x=1-specificity, y=sensitivity), color='blue', size=2) + geom_text_repel(data=best_thresh, mapping=aes(x=1.05-specificity, y=sensitivity ,label=best))
p
基于选定的最优阈值制作混淆矩阵并评估模型预测准确度显著性,结果显著P-Value [Acc > NIR]<0.05
。
predict_result <- data.frame(Predict_status=c(T,F), Predict_class=colnames(prediction_prob))head(predict_result)## Predict_status Predict_class
## 1 TRUE normal
## 2 FALSE tumorpredictions2 <- plyr::join(data.frame(Predict_status=prediction_prob[,1] > best_thresh[1,1]), predict_result)predictions2 <- as.factor(predictions2$Predict_class)confusionMatrix(predictions2, test_data_group)## Confusion Matrix and Statistics
##
## Reference
## Prediction normal tumor
## normal 11 0
## tumor 1 13
##
## Accuracy : 0.96
## 95% CI : (0.7965, 0.999)
## No Information Rate : 0.52
## P-Value [Acc > NIR] : 1.913e-06
##
## Kappa : 0.9196
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9167
## Specificity : 1.0000
## Pos Pred Value : 1.0000
## Neg Pred Value : 0.9286
## Prevalence : 0.4800
## Detection Rate : 0.4400
## Detection Prevalence : 0.4400
## Balanced Accuracy : 0.9583
##
## 'Positive' Class : normal
##
机器学习系列教程
从随机森林开始,一步步理解决策树、随机森林、ROC/AUC、数据集、交叉验证的概念和实践。
文字能说清的用文字、图片能展示的用、描述不清的用公式、公式还不清楚的写个简单代码,一步步理清各个环节和概念。
再到成熟代码应用、模型调参、模型比较、模型评估,学习整个机器学习需要用到的知识和技能。
一图感受各种机器学习算法
机器学习算法 - 随机森林之决策树初探(1)
机器学习算法-随机森林之决策树R 代码从头暴力实现(2)
机器学习算法-随机森林之决策树R 代码从头暴力实现(3)
机器学习算法-随机森林之理论概述
机器学习算法-随机森林初探(1)
机器学习 - 随机森林手动10 折交叉验证
机器学习 模型评估指标 - ROC曲线和AUC值
机器学习 - 训练集、验证集、测试集
一个函数统一238个机器学习R包,这也太赞了吧
基于Caret和RandomForest包进行随机森林分析的一般步骤 (1)
Caret模型训练和调参更多参数解读(2)
基于Caret进行随机森林随机调参的4种方式
机器学习第17篇 - 特征变量筛选(1)
机器学习第18篇 - Boruta特征变量筛选(2)
机器学习第19篇 - 机器学习系列补充:数据集准备和更正YSX包
机器学习第20篇 - 基于Boruta选择的特征变量构建随机森林
机器学习第21篇 - 特征递归消除RFE算法 理论
机器学习第22篇 - RFE筛选出的特征变量竟然是Boruta的4倍之多
机器学习第23篇 - 更多特征变量却未能带来随机森林分类效果的提升
机器学习相关书籍分享
UCI机器学习数据集
送你一个在线机器学习网站,真香!
多套用于机器学习的多种癌症表达数据集
这个统一了238个机器学习模型R包的参考手册推荐给你
莫烦Python机器学习
机器学习与人工智能、深度学习有什么关系?终于有人讲明白了
往期精品(点击图片直达文字对应教程)
机器学习
后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集
一套完整的基于随机森林的机器学习流程(特征选择、交叉验证、模型评估))...相关推荐
- 机器学习项目实战:基于随机森林进行心脏病分类(含多种模型解释方法)
本项目是Kaggle上面的一个经典竞赛题,心脏病分类问题,题目链接在这里. 主要基于随机森林的bagging集成学习框架,通过13个生理特征数据,实现对心脏病分类的预测. 由于自己想要在这个项 ...
- 数据分享|R语言逻辑回归、线性判别分析LDA、GAM、MARS、KNN、QDA、决策树、随机森林、SVM分类葡萄酒交叉验证ROC...
全文链接:http://tecdat.cn/?p=27384 在本文中,数据包含有关葡萄牙"Vinho Verde"葡萄酒的信息(点击文末"阅读原文"获取完整代 ...
- 随机森林c语言编程,一种基于随机森林的C语言源代码静态评分方法与流程
本发明涉及评分领域及机器学习领域,它特别涉及一种基于随机森林的c语言源代码静态评分的构建方法. 背景技术: c语言源代码静态评分是指对根据一定的题目描述做出解答的c语言源代码进行评分.目前,主要流行的 ...
- 基于随机森林算法的贷款违约预测模型研究(Give me some credit)
本文实验代码github:https://github.com/dengsiying/give-me-some-credit/tree/master 参考链接:从决策树到随机森林 http://www ...
- 基于随机森林算法的糖尿病数据集回归
基于随机森林算法的糖尿病数据集回归 作者介绍 1. 随机森林算法原理 1.1决策树与Bagging 1.2 随机森林算法 2. 实验过程 2.1 糖尿病数据集 2.2 实验过程 2.3 实验结果展示 ...
- 基于随机森林算法的多因子选股方法分析与实现(2,代码实现)
摘要 量化投资中经常听到的"多因子模型"是个什么鬼?因子是影响因素的简称,或简单理解成指标.我们都知道股票收益受到多重因素的影响,比如宏观.行业.流动性.公司基本面.交易情绪等等. ...
- kaggle项目:基于随机森林模型的心脏病患者预测分类!
公众号:尤而小屋 作者:Peter 编辑:Peter 大家好,我是Peter~ 新年的第一个项目实践~给大家分享一个新的kaggle案例:基于随机森林模型(RandomForest)的心脏病人预测分类 ...
- RDKit | 基于随机森林的化合物活性二分类模型
基于随机森林算法的化合物二分类机器学习模型 代码示例 #导入依赖包 import pandas as pd import numpy as np from rdkit import Chem, Dat ...
- R语言基于随机森林进行特征选择(feature selection)
R语言基于随机森林进行特征选择(feature selection) 目录 R语言基于随机森林进行特征选择(feature selection)
最新文章
- mysql memcache搭建_Memcached 搭建过程
- mysql支持的并发数_重学MySQL系列(五):谈谈对MySQL的存储引擎的理解
- 李飞飞、邓中翰当选美国国家工程院院士
- mysql安装图形化管理界面phpMyAdmin
- 为什么WebUI里新建任务的状态字段下拉框里没有released选项
- maven远程发布jar
- [财务知识] debt debit credit 的区别于联系
- 几位无人驾驶领域的杰出科学家
- LoftLoader Pro–功能丰富预加载WordPress插件
- Eclipse集成ijkplayer并实现本地和网络视频播放等
- 云服务器学习linux_云服务器怎么选linux系统
- 计算机系统的工作方式,计算机的操作系统与工作方式
- ARP 地址分类 NAT技术
- h5聊天页面 jquery_基于Jquery WeUI的微信开发H5页面控件的经验总结(1)
- 严蔚敏数据结构习题第七章
- 单片机控制步进电机程序c语言正反转停止,51单片机控制步进电机的启动、停止、正转、反转...
- 基于ssm的图书馆预约占座系统 java mysql
- 维基百科中文语料库训练word2vec模型和使用总结
- 性能分析之系统资源饱和度
- 第十六届全国大学生智能车竞赛三个创意组别总决赛比赛成绩
热门文章
- 基于炼铁大数据智能互联平台推动传统工业转型升级
- 《大数据》第1期“专题”——大数据机器学习系统研究进展(下)
- C++——《算法分析》实验壹——二分搜索算法
- 求第k小的数(洛谷P1923题题解,Java/C++语言描述)
- 【Java】命令行下的Java包结构编译与执行
- 【面向对象】面向对象程序设计测试题2-Java基本语法测试题
- Anacoda 的 Python 包管理器 Conda 的常用使用命令
- 香港政府活用无人机,正式应用到调研检测领域
- Python2.7本地安装numpy包
- Selenium Web 自动化 - 项目实战(三)