高分生信SCI套路攻略！精选我最喜欢的3大套路！谁用谁高分！（附代码）
大家好，我是风风。单细胞系列的推文告一段路，我们把基本分析和常见的高级分析基本都走了一遍，剩下的就是实操进行排列组合了。今天我们来聊点新的内容——基于线性模型的特征筛选方法。特征筛选方法大家可能比较熟悉了，像：LASSO回归、随机森林、支持向量机等等。那这么多方法中，哪些是基于线性模型的特征筛选方法呢？又该怎么去实现呢？我们今天一起来看看经典的三种是基于线性模型的特征筛选方法：岭回归、LASSO回归、以及结合了岭回归和LASSO回归两者优势的弹性网络。
数据准备
rm(list = ls())# install.packages('~/推文整理55/ElemStatLearn_2015.6.26.tar.gz', # repos = NULL, # type = 'source') # 按照R包，R包压缩包已经在资料区提供library(ElemStatLearn) # 加载R包library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.3 v purrr 0.3.4## v tibble 3.1.2 v dplyr 1.0.6## v tidyr 1.1.3 v stringr 1.4.0## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --## x dplyr::filter() masks stats::filter()## x dplyr::lag() masks stats::lag()library(corrgram)# 查看数据data(prostate)str(prostate)## 'data.frame': 97 obs. of 10 variables:## $ lcavol : num -0.58 -0.994 -0.511 -1.204 0.751 ...## $ lweight: num 2.77 3.32 2.69 3.28 3.43 ...## $ age : int 50 58 74 58 62 50 64 58 47 63 ...## $ lbph : num -1.39 -1.39 -1.39 -1.39 -1.39 ...## $ svi : int 0 0 0 0 0 0 0 0 0 0 ...## $ lcp : num -1.39 -1.39 -1.39 -1.39 -1.39 ...## $ gleason: int 6 6 7 6 6 6 6 6 6 6 ...## $ pgg45 : int 0 0 20 0 0 0 0 0 0 0 ...## $ lpsa : num -0.431 -0.163 -0.163 -0.163 0.372 ...## $ train : logi TRUE TRUE TRUE TRUE TRUE TRUE ...# 按照官方解释，各个参数意义如下：# lcavol：肿瘤体积的对数值# lweight：前列腺重量的对数值# age：患者年龄# lbph：良性前列腺增生（BPH）量的对数值，非癌症性质的前列腺增生# svi：一个指标变量，表示癌细胞是否已经透过前列腺壁侵入精囊（1 = 是，0 = 否）# lcp：包膜穿透度的对数值，表示癌细胞扩散到前列腺包膜之外的程度# gleason：患者的Gleason评分，表示癌细胞的分化程度——评分越高越严重# pgg45：Gleason评分为4或5所占的百分比（晚期肿瘤）# lpsa：PSA值的对数值# train：用来区分训练数据和测试数据# 查看数据分布特征plot(prostate, col = 'steelblue')
# 我们发现变量svi和gleason图片异常，结合参数意义，发现这两个变量是等级变量，所以接着查看两个变量的分组情况table(prostate$svi)#### 0 1## 76 21table(prostate$gleason)#### 6 7 8 9## 35 56 1 5# 从结果看，gleason评分中，6和7占据大部分，而8和9合计只有6个样本，因此我们对数据进行处理，将6分作为单独一组，而大于6分的样本作为另外一组prostate$gleason <- if_else(prostate$gleason == 6, 0, 1)table(prostate$gleason)#### 0 1## 35 62# 查看各变量之间的相关性corrgram(prostate,lower.panel = panel.bar,upper.panel = panel.conf,col.regions = colorRampPalette(c('red', 'white', 'steelblue')),cor.method = 'pearson')
# 接下来我们区分下训练集和验证集train <- subset(prostate, train == TRUE)[, 1:9]str(train)## 'data.frame': 67 obs. of 9 variables:## $ lcavol : num -0.58 -0.994 -0.511 -1.204 0.751 ...## $ lweight: num 2.77 3.32 2.69 3.28 3.43 ...## $ age : int 50 58 74 58 62 50 58 65 63 63 ...## $ lbph : num -1.39 -1.39 -1.39 -1.39 -1.39 ...## $ svi : int 0 0 0 0 0 0 0 0 0 0 ...## $ lcp : num -1.39 -1.39 -1.39 -1.39 -1.39 ...## $ gleason: num 0 0 1 0 0 0 0 0 0 1 ...## $ pgg45 : int 0 0 20 0 0 0 0 0 0 30 ...## $ lpsa : num -0.431 -0.163 -0.163 -0.163 0.372 ...test <- subset(prostate, train == FALSE)[, 1:9]str(test)## 'data.frame': 30 obs. of 9 variables:## $ lcavol : num 0.737 -0.777 0.223 1.206 2.059 ...## $ lweight: num 3.47 3.54 3.24 3.44 3.5 ...## $ age : int 64 47 63 57 60 69 68 67 65 54 ...## $ lbph : num 0.615 -1.386 -1.386 -1.386 1.475 ...## $ svi : int 0 0 0 0 0 0 0 0 0 0 ...## $ lcp : num -1.386 -1.386 -1.386 -0.431 1.348 ...## $ gleason: num 0 0 0 1 1 0 0 1 0 0 ...## $ pgg45 : int 0 0 0 5 20 0 0 20 0 0 ...## $ lpsa : num 0.765 1.047 1.047 1.399 1.658 ...# 数据准备完成，接下来我们进行不同模型的构建岭回归
# 使用glmnet函数进行岭回归分析library(glmnet)## 载入需要的程辑包：Matrix#### 载入程辑包：'Matrix'## The following objects are masked from 'package:tidyr':#### expand, pack, unpack## Loaded glmnet 4.1-2# 设置自变量和因变量x <- as.matrix(train[, 1:8])y <- train[, 9]# glmnet包会在计算λ值之前首先对输入进行标准化，然后计算非标准化系数，这里我们的因变量lpsa是连续型变量，所以我们指定因变量的分布为gaussian。此外，还需要设置alpha = 0，表示进行岭回归ridge_model <- glmnet(x, # 自变量y, # 因变量（响应变量）family = 'gaussian', # 因变量数据类型alpha = 0 # 0表示进行岭回归)# 结果可视化plot(ridge_model,label = TRUE)
# 上图表示系数值和L1范数之间的关系，而如果我们想看系数值随解释偏差百分比变化的话，可以设置xvar = 'dev'plot(ridge_model, xvar = 'dev', label = TRUE)# 为了将解释偏差百分比与λ值进行对应，我们再画一个Lambda图plot(ridge_model,xvar = 'lambda',label = TRUE)
# 结合上面两张图片可以看出：当λ减小时，系数Coefficient会增大，解释偏差百分比也会增大,如果将λ值设为0，就会忽略掉模型的收缩惩罚。接下来使用coef函数计算模型各参数的系数值ridge_coef <- coef(ridge_model, s = ridge_model$lambda, exact = TRUE)# 指定lambda为0.1时的系数值ridge_coef1 <- coef(ridge_model, s = 0.1)ridge_coef1## 9 x 1 sparse Matrix of class 'dgCMatrix'## s1## (Intercept) 0.130475478## lcavol 0.457279371## lweight 0.645792042## age -0.017356156## lbph 0.122497573## svi 0.636779442## lcp -0.104712451## gleason 0.346022979## pgg45 0.004287179# 量化自变量和因变量之间的关系（类似于risk score值）ridge_train <- predict(ridge_model, newx = x, type = 'response', # 可以设置为“link”, “response”, “coefficients”, “nonzero”, “class”任意一个 s = 0.1)# 用验证集数据进行验证newx <- as.matrix(test[, 1:8])ridge_test <- predict(ridge_model, newx = newx, type = 'response', s = 0.1)# 绘制因变量与预测值之间的关系图plot(ridge_train, y, xlab = 'Predicted', ylab = 'Actual', main = 'Ridge Regression', col = 'steelblue', pch=20, cex = 3)abline(lm (y ~ ridge_train, data = train), col = 'darkgreen', lty = 2) # 绘制回归曲线LASSO回归
# 一样使用glmnet函数进行模型拟合# 设置自变量和因变量x <- as.matrix(train[, 1:8])y <- train[, 9]lasso_model <- glmnet(x,y,family = 'gaussian',alpha = 1)# lasso_model# 和岭回归一样，我们可以对结果进行可视化plot(lasso_model,xvar = 'dev',label = TRUE)
plot(lasso_model, xvar = 'lambda', label = TRUE)# 计算系数值，观察lasso_model的结果，如果尝试拟合7个自变量构成的模型，从lasso_model的7变量模型中观察到对应的Lambda最小为0.045，根据此参数进行系数计算lasso_coef <- coef(lasso_model, s = 0.045)lasso_coef## 9 x 1 sparse Matrix of class 'dgCMatrix'## s1## (Intercept) -0.1305900670## lcavol 0.4479592050## lweight 0.5910476764## age -0.0073162861## lbph 0.0974103575## svi 0.4746790830## lcp .## gleason 0.2968768129## pgg45 0.0009788059# 发现lcp的系数为0，表示lcp并未纳入模型，而其它自变量都有一一对应的系数值，这表示LASSO算法在λ值为0.045时，将lcp的系数归零，接下来一样使用predict函数lasso_train <- predict(lasso_model,newx = x,type = 'response',s = 0.045)# 用验证集数据进行验证newx <- as.matrix(test[, 1:8])lasso_test <- predict(lasso_model,newx = newx,type = 'response',s = 0.045)# 绘制因变量与预测值之间的关系图plot(lasso_train,y,xlab = 'Predicted',ylab = 'Actual',main = 'LASSO Regression',col = 'purple',pch=20,cex = 3)
# 我们发现在这里，岭回归和LASSO回归的统计图基本一致，说明两模型结果差别不大，接下来我们尝试进行弹性网络拟合弹性网络
library(caret)## 载入需要的程辑包：lattice#### 载入程辑包：'lattice'## The following object is masked from 'package:corrgram':#### panel.fill#### 载入程辑包：'caret'## The following object is masked from 'package:purrr':#### lift# caret包生成α值和λ值grid <- expand.grid(.alpha = seq(0, 1, by = .2), .lambda = seq(0.00, 0.2, by = 0.02))# 查看α值和λ值的组合table(grid) # α值在0和1之间，每次增加0.2，λ值在0和0.2之间，每次增加0.02## .lambda## .alpha 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2## 0 1 1 1 1 1 1 1 1 1 1 1## 0.2 1 1 1 1 1 1 1 1 1 1 1## 0.4 1 1 1 1 1 1 1 1 1 1 1## 0.6 1 1 1 1 1 1 1 1 1 1 1## 0.8 1 1 1 1 1 1 1 1 1 1 1## 1 1 1 1 1 1 1 1 1 1 1 1# 设置重取样方法LOOCVcontrol <- trainControl(method = 'LOOCV')# 使用train()函数确定最优的弹性网络参数enet_train <- train(lpsa ~ .,data = train,method = 'glmnet',trControl = control,tuneGrid = grid)head(enet_train[['results']])## alpha lambda RMSE Rsquared MAE## 1 0 0.00 0.7498311 0.6091434 0.5684548## 2 0 0.02 0.7498311 0.6091434 0.5684548## 3 0 0.04 0.7498311 0.6091434 0.5684548## 4 0 0.06 0.7498311 0.6091434 0.5684548## 5 0 0.08 0.7498311 0.6091434 0.5684548## 6 0 0.10 0.7508038 0.6079576 0.5691343# 由结果可以看到，RMSE最小值为0.7498，对应的最优组合为alpha = 0，lambda = 0.08，也就是相当于glmnet中s = 0.08的岭回归，接下来的过程就是岭回归的过程了x <- as.matrix(train[, 1:8])y <- train[, 9]enet <- glmnet(x, y, family = 'gaussian', nfolds = 10, alpha = 0, lambda = 0.08)enet_coef <- coef(enet, s = 0.08, exact = TRUE)enet_coef## 9 x 1 sparse Matrix of class 'dgCMatrix'## s1## (Intercept) 0.137811097## lcavol 0.470960525## lweight 0.652088157## age -0.018257308## lbph 0.123608113## svi 0.648209192## lcp -0.118214386## gleason 0.345480799## pgg45 0.00447826
交叉验证
# cv.glmnet进行交叉验证set.seed(565656)x <- as.matrix(train[, 1:8])y <- train[, 9]cv_res <- cv.glmnet(x, y, nfolds = 10)plot(cv_res)# 图片表示的是λ的对数值和均方误差MSE之间的关系，还带有模型中显著特征的参数数量。图中两条垂直的虚线表示取得MSE最小值的logλ值（左侧虚线）和距离最小值一个标准误差的logλcv_res$lambda.min # 取得MSE最小值的logλ值## [1] 0.003308928# 计算系数（这里用弹性网络的模型）coef_enet <- coef(enet, s = cv_res$lambda.1se)coef_enet## 9 x 1 sparse Matrix of class 'dgCMatrix'## s1## (Intercept) 0.137811097## lcavol 0.470960525## lweight 0.652088157## age -0.018257308## lbph 0.123608113## svi 0.648209192## lcp -0.118214386## gleason 0.345480799## pgg45 0.004478267coef_enet_cv <- predict(cv_res, newx = x, type = 'response', s = 'lambda.1se')plot(coef_enet_cv,y,xlab = 'Predicted',ylab = 'Actual',main = 'LASSO Regression',col = 'purple',pch=20,cex = 3)
《精通机器学习——基于R》
《R语言统计分析与应用》
《R语言实战》
《卫生统计学（第8版）》
赞 (0)