github-actions[bot]
commented
2 years ago 机器学习系列|R语言实践机器学习案例分享-Bootstrap的R语言实现 by 科研生信充电宝
R语言实践
机器学习案例
R
Bootstrap的R语言实现
- 05·21-
宜分享
R语言实践机器学习案例迎来第五期啦~
昨天python实现Bootstrap的推送大受欢迎。有小伙伴留言说，能否用R语言实现Bootstrap呢？
本期，为大家带来《Bootstrap的R语言实现》
文献示例1：
图. 科研文章中常见使用Bootstrap计算置信区间
文献示例2：
图. 科研文章中常见使用Bootstrap来重构模型
至于该文献中的这种方法是否能够评价模型的内部稳健性呢？暂且存疑、、、、
在之前的文章《一文理清重采样（Resampling）》总结系列|一文理清重采样（Resampling）中，我们已经详细介绍了Bootstrap的基本原理及其Python实现。
下面简要回顾下：
模型选择与评价的步骤 .
① 抽样（sampling）：将数据分为训练数据（train data）、验证数据（validate data）、测试数据（test data），抽样方法有保留方法抽样、分层抽样、过采样、自助法抽样、袋外抽样；
② 交叉验证（cross-validation）：交叉验证是将数据分成k等份，每次将取出k-1部分作训练数据，另为一部分作验证数据。这样可以选择更正确的模型；
③ 模型选择：模型的复杂度、偏差与方差、模型调参；
④ 测试数据预测：模型选择出来，要利用测试数据预测结果；
⑤ 模型评价：将预测结果和测试数据的目标变量值做比较。目标变量是分类因子，评价用混淆矩阵；目标变量是数值，评价用误差均方和；
图. 模型选择与评价步骤
大数据抽样方法 .
目前大数据的抽样方法很多，具体可以见上一篇推送。这里概括如下图：
图. 大数据的抽样方法
而在R语言的使用中，常常会用到如下几种抽样方法：
图. R语言抽样方法的名称和参数
常见的有:
① 留出法抽样（Holdout）：
# 留出法抽样；data <- mtcarsstr(data) # 11个变量，32个观察set.seed(2022)rand <- order(runif(32)) # 1到32不重置随机抽样train <- data[rand[1:20],] # 20个训练数据validate <- data[rand[21:25],] # 5个验证数据test <- data[rand[26,32],] # 7个测试数据str(train);str(validate);str(test)table(data$mpg) # 查看目标类别
② 分层抽样：十分常用
# 分层抽样library(caret)data <- mtcarsstr(data) # 11个变量，32个观察set.seed(2022)# 根据变量mpg进行分层抽样75%为训练集trval <- createDataPartition(data$mpg,p=0.75,list=FALSE)train <- data[trval,] # 训练集test <- data[-trval,] # 测试集str(train)table(train$mpg)table(test$mpg)
Bootstrap重抽样思维 .
Bootstrap是一种统计学思维逻辑：
随机森林算法就用到了这种思维：
今天，我们继续介绍Bootsrap的R语言实现，分为以下几个部分：
赞赏本文，即可解锁资源哦~
Boot包中的自助法
Brief Introduction
一、概述：
Boot包扩展了自助法和重抽样的相关用途。可以实现对一个统计量或一个统计量向量使用自助法。
二、步骤：
① 写一个能返回待研究统计量值的函数；
② 使用boot()函数生成有效的统计量重复样本；
③ 使用boot.ci()函数生成统计量的置信区间；
三、主要函数：
boot()函数：
# boot()函数格式bootobject <- boot(data=,statistic=,R=,……)
参数描述：
函数意义：
boot()函数可调用统计量函数R次，每次都从整数1：nrow(data)中生成一列有放回的随机指标，这些指标被统计量函数用来选择样本。
函数返回值：
bootobject$t0:从原始数据中得到的k个统计量的观测值；
bootobject$t:一个R*k矩阵，每行即k个统计量的自助重复值；
查看结果：
可以使用print()函数和plot()函数检查结果。
boot.ci()函数：
# boot.ci()函数格式boot.ci(bootobject,conf=,type=)
参数描述：
bootobject：boot()函数返回的对象；
conf：预期的置信区间（默认conf=0.95）
type：返回置信区间的类型，可为：norm、basic、stud、perc、bca和all（默认type="all"）
type参数注释：
perc方法，即分位数，展示的是样本均值；
bca方法，根据偏差对区间做简单调整，并且在大部分情况中都是更可取的；
四、对单个统计量使用自助法：
例：针对mtcars数据集，使用车重（wt）和发动机排量（disp）来构建多元回归模型预测汽车每加仑行驶的英里数（mpg）。此时，审稿人希望你给出R2的95%置信区间。
① 构建获取R2的函数：
# ① 构建获取R2的函数：rsq <- function(formula,data,indices){d <- data[indices,]fit <- lm(formula,data=d)return(summary(fit)$r.square)}
② 进行1000次bootstrap：
# ② 进行1000次bootstraplibrary(boot)set.seed(2022)results <- boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)
③ 输出结果：
# ③ 输出结果print(results)
④ 可视化结果：
# ④ 可视化结果plot(results)
解读：自助法求得的R2值不呈正态分布。
⑤ 计算95%置信区间：
# ⑤ 使用两种方法计算95%置信区间：boot.ci(results,type=c("perc","bca"))
解读：不同方法生成的置信区间不同。
五、对多个统计量使用自助法：
例：在上例的基础上，计算三个回归系数（截距项、车重、发动机排量）的95%置信区间。
① 创建返回回归系数向量的函数：
# 对多个统计量使用自助法：# ① 创建返回回归系数向量的函数：bs <- function(formula,data,indices){d <- data[indices,]fit <- lm(formula,data=d)return(coef(fit))}
② 自助抽样1000次：
# ② 自助抽样1000次：library(boot)set.seed(2022)results <- boot(data=mtcars,statistic = bs,R=10000,formula=mpg~wt+disp)
③ 输出结果：
# ③ 输出结果：print(results)
解读：在本例中t1指截距项；t2指wt；t3指disp；
④ 可视化结果：
针对多个统计量使用自助法时，在可视化和计算置信区间时，需添加一个索引参数来指明t的列。
# ④ 可视化结果：plot(results,index=2)# wt
⑤ 计算95%置信区间：
# ⑤ 计算95%置信区间：boot.ci(results,type="bca",index=2) # wt
回归预测模型的自举重采样验证
Brief Introduction
例：针对mtcars数据集，使用所有自变量来构建多元回归模型预测汽车每加仑行驶的英里数（mpg），在模型构建后使用100个重采样的自举法来测试线性回归模型。
① 加载包：
# 回归预测模型的自举重采样验证# ① 加载包：library(tidyverse)library(caret)
② 加载数据：
# ② 加载数据data=mtcars# 查看数据（随机选取3行数据来查看）sample_n(data,3)
③ 构建回归模型：
# ③ 构建回归模型：model <- lm(mpg~.,data=data)print(model)
计算模型的R2和RMSE：
# 计算模型的R2和RMSEpredictions <- model%>%predict(data)# 模型评价data.frame(RMSE=RMSE(predictions,data$mpg),R2=R2(predictions,data$mpg))
④ 返回回归模型系数：
该实现步骤同前
# ④ 返回回归模型系数# 构建返回的函数model_coef <- function(data,index){ coef(lm(mpg ~.,data =data,subset =index))}model_coef(data,1:32) # 32为样本数
⑤ 计算100个自举数据集的各个截距和预测变量的系数的标准误差：
# 计算100个自举数据集的各个截距和预测变量的系数的标准误差：library(boot)results <- boot(data,model_coef,100)print(results)
⑤ 返回回归系数的置信区间：
# ⑥ 计算系数的置信区间boot.ci(results,type="bca",index=2)
⑥ 对比：
传统lm()方法有四个假设：
Linearity 线性：应变量和每个自变量都是线性关系；
Indpendence 独立性：对于所有的观测值，它们的误差项相互之间是独立的；
Normality 正态性：误差项服从正态分布；
Equal-variance 等方差：所有的误差项具有同样方差；
传统lm()参数查看：
# 传统lm()参数查看summary(lm(mpg ~.,data =data))$coef
和bootstrap结果对比：
解读：回归系数基本相同，但标准差相对来说要大一些。自举重采样验证法不依赖于线性模型做出的任何假设，因此，它可能提供比summary（）函数更准确的系数标准误差估计。
Cox回归模型的的自举重采样验证
Brief Introduction
一般使用Harrell’s C-index来评价Cox回归模型的拟合效果。
下面使用之前的示例数据集来进行Cox回归模型的的自举重采样验证。
① 加载包：
library(caret) # 划分训练集和测试集library(Hmisc) # 计算C统计量library(survival)library(rms) # 拟合模型
② 加载数据集：
# boot()函数格式# ② 加载数据surdata <- read.csv("F:\\05_15rcs.csv")# 查看数据（随机选取3行数据来查看）sample_n(surdata,3)
③ 处理数据集并划分：
# 构建数据集surdatanew <- surdata[,2:5]sample_n(surdatanew,3)# 使用caret包简单划分训练集和测试集set.seed(2022)ind <- createDataPartition(surdatanew$status,p = 0.8,list = F,times = 1)train <- surdatanew[ind,]table(train$status)test <- surdatanew[-ind,]table(test$status)
④ 使用训练集构建Cox模型：
# ④ 使用训练集构建Cox模型：library(Hmisc) # 计算C统计量library(survival)library(rms)fit <- cph(Surv(time,status)~gxr+gender,x=TRUE,y=TRUE,surv=TRUE,data=train)# 查看拟合结果fit
⑤ 计算训练集的C统计量：
# ⑤ 计算训练集的C统计量：cindex <- as.data.frame(1-rcorr.cens(predict(fit,train),Surv(train$time,train$status)))[1,1]cindex
输出：0.7030315
⑥ 计算测试集的C统计量：
# 计算测试集的C统计量cindextest <- as.data.frame(1-rcorr.cens(predict(fit,test),Surv(test$time,test$status)))[1,1]cindextest
输出：0.6658172
⑦ 使用Bootstrap自助法计算训练集的Harrell’s C-index的置信区间：Bootstrap200次
# 使用Bootstrap自助法计算Harrell’s C-index的置信区间library(boot)cindexbootstrap <- function(data,indices){d <- data[indices,]cindex <- as.data.frame(1-rcorr.cens(predict(fit,d),Surv(d$time,d$status)))[1,1]return(cindex)}results <- boot(data=train,statistic = cindexbootstrap,R=200)print(results)
解读：这里显示了C-index的均值和标准差；
# 输出置信区间boot.ci(results,type=c("norm","basic","perc"))
解读：这里分别使用了三种方法计算了置信区间。
⑧ 使用Bootstrap自助法计算测试集的Harrell’s C-index的置信区间：Bootstrap200次
# 使用Bootstrap自助法计算测试集Harrell’s C-index的置信区间library(boot)cindexbootstrap <- function(data,indices){d <- data[indices,]cindex <- as.data.frame(1-rcorr.cens(predict(fit,d),Surv(d$time,d$status)))[1,1]return(cindex)}results <- boot(data=test,statistic = cindexbootstrap,R=200)print(results)plot(results)# 输出置信区间boot.ci(results,type=c("norm","basic","perc"))
使用caret包构建基于bootstrap的机器学习模型
Brief Introduction
在以上这几个板块中，我们都是在已经构建好的模型（线性回归或Cox回归）上，使用bootstrap方法来计算一个或几个统计量的置信区间。而在模型构建阶段，我们都没有采用bootstrap技术来提升模型性能。
借助于caret包，我们可以实现基于bootstrap的机器学习模型构建。
caret包 .
全称：Classification and Regression Training，其是一系列函数的集合。
它试图对创建预测模型的过程进行流程化。主要包括图形可视化，数据预处理，特征选择，数据分割，模型的构建，抽样、模型调参与模型的预测等方面，caret包为 R 中超过 200 个包提供了统一的接口，方便了我们的使用。
在进行建模时，caret包的常用函数有：
① modelLookup()：访问模型的参数；
② train()：训练模型；
③ trainControl()：调整模型；
④ expand.grid()：调整模型参数；
⑤ predict()：预测测试数据；
在caret包中，有两种方法实现bootstrap。
① 一种是在数据分割的时候：
传统的分层分割方法是使用createDataPartition()函数：
# 加载数据data <- mtcars# 传统数据分割set.seed(2022)ind <- createDataPartition(data$mpg,# 指定目标列别p = 0.8,# 训练集和测试集的比例list = F,# FALSE不将数据以列表方式返回times = 1)# 创建拆分的数量train <- data[ind,] # 训练集table(train$mpg) # 训练集中目标的类别统计test <- data[-ind,] # 测试集table(test$mpg) # 测试集中目标的类别统计
分层k折交叉验证的分割方法是使用createfold()函数：
# 分层k折交叉验证set.seed(2022)indcv <- createFolds(data$mpg,# 指定目标列别k = 10,# 折数，默认为10list = T,# 是否返回抽样的真实值，默认返回索引值returnTrain = F)indcv
bootstrap的分割方法是使用createResample()函数：
# bootstrap数据分割set.seed(2022)indbootstrap <- createResample(data$mpg,# 指定目标列别times = 10, # 指定抽样组数，默认为10组list = T)# TRUE将数据以列表方式返回indbootstrap
这种方法可以实现在数据划分时实现bootstrap，并且返回每一个bootstrap样本集都选用了哪些样本。但这个方法直接用于后续的模型构建的略微复杂。
② 在训练模型的时候：
在进行建模时，caret包中主要使用：train（）函数。
train（）函数格式如下：
# train（）函数格式train(x, y, method = "rf", preProcess = NULL, ...,weights = NULL, metric = ifelse(is.factor(y), "Accuracy", "RMSE"),maximize = ifelse(metric %in% c("RMSE", "logLoss", "MAE"), FALSE, TRUE),trControl = trainControl(), tuneGrid = NULL,tuneLength = ifelse(trControl$method == "none", 1, 3))
其中：
x：行为样本，列为特征的矩阵或数据框。列必须有名字；
y：每个样本的结果，数值或因子型；
method：指定具体的模型形式，支持大量训练模型，可在此查询：http://topepo.github.io/caret/available-models.html；
preProcess ：代表自变量预处理方法的字符向量。默认为空，可以是 "BoxCox", "YeoJohnson", "expoTrans", "center", "scale", "range", "knnImpute", "bagImpute", "medianImpute", "pca", "ica" and "spatialSign"；
weights：加权的数值向量。仅作用于允许加权的模型
metric：指定将使用什么汇总度量来选择最优模型。默认情况下，"RMSE" and "Rsquared" for regression and "Accuracy" and "Kappa" for classification；
maximize：逻辑值，metric是否最大化；
trControl：定义函数运行参数的列表。具体见下；
tuneGrid：可能的调整值的数据框，列名与调整参数一致；
tuneLength：调整参数网格中的粒度数量，默认时每个调整参数的level的数量；
接下来，再具体介绍下trainControl()函数：
# trainControl()函数：trainControl(method = "boot", number = ifelse(grepl("cv", method), 10, 25),repeats = ifelse(grepl("[d_]cv$", method), 1, NA), p = 0.75,search = "grid", initialWindow = NULL, horizon = 1,fixedWindow = TRUE, skip = 0, verboseIter = FALSE, returnData = TRUE,returnResamp = "final",.....)
其中：
method：重抽样方法："boot", "boot632", "optimism_boot", "boot_all", "cv", "repeatedcv", "LOOCV", "LGOCV" (for repeated training/test splits), "none" (only fits one model to the entire training set), "oob" (only for random forest, bagged trees, bagged earth, bagged flexible discriminant analysis, or conditional tree forest models), timeslice, "adaptive_cv", "adaptive_boot" or "adaptive_LGOCV"
number：folds的数量或重抽样的迭代次数；
repeats：仅作用于k折交叉验证：代表要计算的完整折叠集的数量；
p：仅作用于分组交叉验证：代表训练集的百分比；
search：Either "grid" or "random"，表示如何确定调整参数网格；
由此，我们可以根据trainControl()函数中的method来选择bootstrap方法。
下面以构建线性回归模型为例：
对比几种不同的抽样方法：
① 加载数据-随机排序-分层抽样
# 分层抽样library(caret)data <- mtcars# 随机排序set.seed(2022)rows <- sample(nrow(data))data <- diamonds[data, ]str(data) # 11个变量，32个观察# 根据变量mpg进行分层抽样75%为训练集set.seed(2022)trval <- createDataPartition(data$mpg,p=0.75,list=FALSE)trainset <- data[trval,] # 训练集testset <- data[-trval,] # 测试集str(trainset)table(trainset$mpg)table(testset$mpg)
② 构建传统线性回归模型：
# 在训练集构建传统模型model <- lm(mpg~.,data=trainset)print(model)# 模型评价# 计算模型在测试集的R2和RMSEpredictions <- model%>%predict(testset)data.frame(RMSE=RMSE(predictions,testset$mpg),R2=R2(predictions,testset$mpg))# 计算模型在训练集的R2和RMSEpredictions <- model%>%predict(trainset)data.frame(RMSE=RMSE(predictions,trainset$mpg),R2=R2(predictions,trainset$mpg))
③ 构建caret线性回归模型：
# 在训练集构建caret线性回归模型library(tidyverse)library(caret)model<-train(mpg~.,data=trainset,method= "lm")print(model)# 模型评价# 计算模型在测试集的R2和RMSEpredictions <- model%>%predict(testset)data.frame(RMSE=RMSE(predictions,testset$mpg),R2=R2(predictions,testset$mpg))# 计算模型在训练集的R2和RMSEpredictions <- model%>%predict(trainset)data.frame(RMSE=RMSE(predictions,trainset$mpg),R2=R2(predictions,trainset$mpg))
④ 构建基于bootstrap的线性回归模型：
# 在训练集构建基于bootstrap的回归模型：set.seed(2022)train.control<-trainControl(method= "optimism_boot",number= 500)model<-train(mpg~.,data=trainset,method= "lm", trControl=train.control,metric = "R2")print(model)# 计算模型在测试集的R2和RMSEpredictions <- model%>%predict(testset)data.frame(RMSE=RMSE(predictions,testset$mpg),R2=R2(predictions,testset$mpg))# 计算模型在训练集的R2和RMSEpredictions <- model%>%predict(trainset)data.frame(RMSE=RMSE(predictions,trainset$mpg),R2=R2(predictions,trainset$mpg))
注意：这里给出的结果是Resampling后的平均结果。
⑤ 构建基于多重K折交叉验证的线性回归模型：
# 在训练集构建基于repeatedcv的回归模型set.seed(2022)train.control<-trainControl(method= "repeatedcv",number= 5,repeats=10)model<-train(mpg~.,data=data,method= "lm", trControl=train.control)print(model)# 计算模型在测试集的R2和RMSEpredictions <- model%>%predict(testset)data.frame(RMSE=RMSE(predictions,testset$mpg),R2=R2(predictions,testset$mpg))# 计算模型在训练集的R2和RMSEpredictions <- model%>%predict(trainset)data.frame(RMSE=RMSE(predictions,trainset$mpg),R2=R2(predictions,trainset$mpg))
注意：这里给出的结果是Resampling后的平均结果。
解读：通过对比不同模型在测试集中的表现，可以看到多重K折交叉验证的model性能好于bootstrap和原始lm()模型。
附
Brief Introduction
下面利用Boston数据集，简要对比下10次随机留出法和10次10折交叉验证的区别：
① 10次随机留出法：
重复运用验证集方法10次，每次用一种不同的随机分割把观测分为一个训练集和一个测试集
# 10次留出法library(MASS)set.seed ( 1 )train1 <- sample ( 506 , 506 / 2 )lmfit1 <- lm ( medv ~. , data = Boston , subset = train1 )attach ( Boston )pred1 <- predict ( lmfit1 , Boston [ - train1 , ] )mean ( ( medv [ - train1 ] - pred1 ) ^ 2 )err1 <- rep ( 0 , 10 )for ( i in 1 : 10 ) {train2 <- sample ( 506 , 506 / 2 )lmfit2 <- lm ( medv ~. , data = Boston , subset = train2 )pred2 <- predict ( lmfit2 , Boston [ - train2 , ] )err1 [ i ] <- mean ( ( medv [ - train2 ] - pred2 ) ^ 2 )}plot ( 1 : 10 , err1 , xlab = "" , ylim = c ( 20 , 30 ) , type = "l" ,main = "选取10个不同的训练集对应的测试误差" )detach ( Boston )
② 10次10折交叉验证：
运用10折CV方法10次，每次用一种不同给的随机分割把数据分为10个部分。
# 10折交叉验证library(boot)set.seed ( 3 )glmfit2 <- glm ( medv ~. , data = Boston )cv.err2 <- cv.glm ( Boston , glmfit2 , K = 10 )cv.err2 $ deltaerr2 <- rep ( 0 , 10 )for ( i in 1 : 10 ) {glmfit3 <- glm ( medv ~. , data = Boston )cv.err3 <- cv.glm ( Boston , glmfit3 , K = 10 )err2 [ i ] <- cv.err3 $ delta [ 1 ]}plot ( 1 : 10 , err2 , xlab = "" , ylim = c ( 20 , 30 ) , type = "l" ,main = "10次不同的CV误差" )
解读：可以看出，使用验证集方法所产生的测试均方误差具有较大的波动性。
总结
Brief Introduction
在今天的学习中，我们回顾了模型选择与评价的步骤：
复习了大数据抽样的方法：
重点介绍了R语言中的Bootstrap方法实现：
最后通过模型比较，我们发现，采取多次K-折交叉验证的效果最好。
附：caret包使用指南
现在：
长按扫码关注：科研生信充电宝
10元赞赏本文，即喜欢作者~
即可直接解锁：
《Bootstrap的R语言实现》全套代码
看到这里你还不心动吗？
赶紧关注、转发、点赞、分享，领取你的专属福利吧~
好啦，以上就是今天推文的全部内容啦！
今天的分享就到这里啦~
另外，本公众号建立了R语言和Python学习分享群，联系后台，拉您入群。
版权声明：本文内容由互联网用户自发贡献，版权归作者所有，本公众号不拥有所有权，也不承担相关法律责任。
如果您发现本公众号中有涉嫌抄袭的内容，欢迎发送邮件至：kysxcdb@163.com 进行举报，一经查实，本公众号将立刻删除涉嫌侵权内容。
BIOINFOR
微信号｜科研生信充电宝
· BIOINFOR ·
永远相信美好的事情
即将发生
https://mp.weixin.qq.com/s/12_IsC599stYeTWJ4yWY4w