前言
XGBoost（eXtreme Gradient Boosting）和LightGBM（Light Gradient Boosting Machine）是两种流行的机器学习框架，特别适用于梯度提升算法。目前，在很多文章中，尤其是非肿瘤的文章中应用广泛。XGBoost和LightGBM生成的模型被视为“黑箱模型”，即通过常规的重要性评估，我们只能确定对模型最关键的特征。然而，对于这些关键特征如何影响模型的最终结果指标，我们却一筹莫展。因此，2017年，Lundberg和Lee的论文提出了SHAP值这一通用方法，用于解释各种模型（包括分类和回归）。其中最大的受益者是之前难以理解的黑箱模型，例如XGBoost和LightGBM模型。下面跟着果叔来复现下面文章中的机器学习过程吧！
今天果叔和大家分享的是2022年12月发表在 “Frontiers in Immunology” 上的分的非肿瘤免疫分型结合多种机器学习算法文章——《Identification of immune microenvironment subtypes and signature genes for Alzheimer's disease diagnosis and risk prediction based on explainable machine learning》。
原文中机器学习模型的应用和解释
1.通过SHAP值阐明每个特征变量对预测模型的影响。基于 XGBoost 模型的 SHAP 汇总图对特征变量的重要性排序表明，对 XGBoost 模型贡献最大的前 5 个变量是 CXCR4、PPP3R1、HSP90AB1、CXCL10 和 S100A12（图6A）。
2. 采用 SHAP 依赖分析来描述单个特征变量如何影响 XGBoost 预测模型的结果（图 6B）。特征变量的 SHAP 值越高，AD 的可能性就越大。
3.LightGBM模型中最重要的5个变量与XGBoost模型中的变量一致（图6C）。SHAP 依赖分析也用于解释每个特征变量对 LightGBM 模型输出的影响（图6D）。
下面开始我们的实战了，果叔提醒大家本次复现的内容运算起来比较费时、占用内存比较大，建议各位小伙伴可以使用服务器运行，这样就不会出现跑不动或者黑屏的情况了！！如果没有自己的服务器欢迎联系我们使用服务器租赁~价格便宜到超乎你想象！！真的是花小价钱办大事，性价比超级高哦！！！
跑代码时卡顿、电脑不给力让人抓狂！找大侠试用稳定高速的服务器，让分析顺畅无比！
代码学不会？bug 频繁出现，束手无策？实操生信分析课程赶快学起来！滴滴大侠领取体验课程哦~
线上课程教学
课题设计、定制生信分析
云服务器租赁
加微信备注99领取使用
代码复现
首先，我们来看看如何构建XGBoost机器学习模型并通过模型进行特征选择。
XGBoost 机器学习模型是一种高效且可扩的展的机器学习分类器，由 Chen 和 Guestrin 在 2016 年推广。XGBoost原理是是在决策树的基础上产生迭代，它以 boosting 的方式结合了多个决策树。通常创建每棵新树是为了通过梯度提升来减少先前模型的误差，误差指的是实际值和预测值之间的差异。把误差作为协变量参与下一个模型的预测，反复执行这个过程，降低出错率，直到决策树指定阈值，模型已经被训练成功。XGBoost 具有与梯度提升相同的原理，它使用提升次数、学习率、二次采样率和最大树深度来控制过度拟合并增强更好的性能。
Xgboost相对于线性模型在进行预测时往往有更好的精度，但是同时也失去了线性模型的可解释性。所以Xgboost通常被认为是黑箱模型。
2017年，Lundberg和Lee的论文提出了SHAP值这一广泛适用的方法用来解释各种模型（分类以及回归），其中最大的受益者莫过于之前难以被理解的黑箱模型，如boosting和神经网络模型。
本教程中，我们在真实数据集上进行实操，利用SHAP来解释Xgboost模型。
XGBoost
shapviz: 提供了用于SHAP值可视化的函数。
ggplot2: 用于绘制图形。
xgboost: XGBoost模型。
##加载相关的R包= ls())library(ggplot2)library(xgboost)inputFile="SHAP.csv" #输入文件#读取输入文件data=data.table::fread(inputFile)data=data[,-1]data=t(data)data=t(data)dat <-data#对数据进行一些预处理，包括转置（t）和删除缺失值（na.omit）。dat <- na.omit(dat)vfactor <- c("group")set.seed(123)dtrain <- xgb.DMatrix(data.matrix(dat[,-1]),label = dat[,1])# 构建模型#使用xgb.DMatrix将数据转换为XGBoost的数据格式。#使用xgboost函数构建XGBoost模型，指定一些参数，如eta（学习率）、eval_metric（评估指标）、nrounds（迭代次数）等。fit <- xgboost(params = list(eta = 0.1,eval_metric = "logloss",nthread=2,eval_metric="auc"),data = dtrain,prediction=TRUE,nrounds = 100)#使用shapviz中的函数计算SHAP值，并通过sv_waterfall、sv_force、sv_importance、sv_dependence等函数对SHAP值进行可视化。shap <- shapviz(fit,X_pred =data.matrix(dat[,-1]))row_id = 3,max_display = 30,fill_colors=c("#FF0000", "#0085FF"))sv_force(shap,row_id = 2,max_display = 30,fill_colors=c("#FF0000", "#0085FF"))kind = "beeswarm",max_display = 30)# 绘制条形图= "#0085FF",max_display = 30)
以上是对模型重要的特征。然而，如何去解释以上特征对模型的贡献呢？这个时候就需要计算SHAP值。
#SHAP dependence plots（SHAP依赖图）v = "SERPINA3", ##这个地方注意要换成自己的基因哦！！color_var = NULL)#计算并可视化SHAP值的交互作用，使用sv_interaction函数。shp_i <- shapviz(fit,X_pred =data.matrix(dat[,-1]),interactions = TRUE)= element_text(angle = 45, vjust = 1, hjust = 1))
图解：变量名沿纵坐标表示，模型的重要性从上到下递减。紧接着的数值是所有样本SHAP值的平均值。图中每个点代表一个样本的SHAP值，颜色接近黄色表示值较大，接近紫色表示值较小。
LightGBM
使用LightGBM模型计算SHAP（Shapley Additive exPlanations）值，并使用shapviz包中的可视化函数对模型的解释结果进行展示。SHAP值提供了对模型预测输出的解释，以及每个特征对于预测结果的贡献程度。
LightGBM在XGBoost上主要有3方面的优化。
1，Histogram算法:直方图算法。
2，GOSS算法:基于梯度的单边采样算法。
3，EFB算法:互斥特征捆绑算法。
XGBoost生成一片叶子的复杂度可以粗略估计为 = 特征数量*候选分裂点数量*样本的数量。
Histogram算法的主要作用是减少候选分裂点数量，
GOSS算法的作用是减少样本的数量，
EFB算法的作用是减少特征的数量。
通过这3个算法的引入，LightGBM生成一片叶子需要的复杂度大大降低了，从而极大节约了计算时间。同时Histogram算法还将特征由浮点数转换成0~255位的整数进行存储，从而极大节约了内存存储。
加载R包
#shapviz: 提供了用于SHAP值可视化的函数。#ggplot2: 用于绘制图形。#lightgbm: LightGBM模型。library(shapviz)library(ggplot2)library(lightgbm)inputFile="SHAP.csv" #输入文件#读取输入文件data=data.table::fread(inputFile)data=data[,-1]data=t(data)data=t(data)dat <-datadat <- na.omit(dat)vfactor <- c("group")<- lapply(df[vfactor], factor)set.seed(123)#模型构建#使用LightGBM的lgb.Dataset将数据转换为LightGBM的数据格式。#使用lgb.train函数构建LightGBM模型，指定一些参数，如objective（目标函数）、learning_rate（学习率）等。dftrain <- lgb.Dataset(data=data.matrix(dat[,-1]),label = dat[,1])lgb <- lgb.train(params = list(objective="binary",learning_rate=0.3),data=dftrain,nrounds=1000,verbose = -1)# 计算SHAP值#使用shapviz中的函数计算SHAP值，并通过sv_waterfall、sv_importance等函数对SHAP值进行可视化。shp <- shapviz(lgb,X_pred =data.matrix(dat[,-1]))#waterfall plot= 2,max_display = 30)
#SHAP summary plot（变量重要性）sv_importance(shp, kind = "beeswarm",max_display = 30)# 条形图sv_importance(shp,fill="#0085FF",max_display = 30)
结语
好啦，今天的学习就到此结束了，小伙伴们快去试试吧！想要今天示例数据的小伙伴们，赶快去后台领取吧！！最后，大家如果对生信分析感兴趣但还不熟悉，又想尝试一下处理自己的数据，不妨试一下果叔开发的生信云平台哦，一键出图，一键导出CNS级别的Figture！！赶快去试试吧！！点击 http://www.biocloudservice.com/home.html。
不会分析还想用生信工具助力发文咋办？有这顾虑的朋友，想一步到位就带着想法来，不论是代码实操还是在线文章结果复现，果叔照样能提供，还有大家都想要的服务器，找果叔获取就对了！
定制生信分析
服务器租赁
扫码咨询果叔
往期回顾
01 “顶刊中的小弟，水刊中的大哥”说的就是Nature Communications！杭师大团队靠UKB数据库，简简单单5张图搞定！ |
02 |
03 |
04 IF:18.8!这套思路转化价值极高！深度学习+多模态数据+生物标志物，正成为国自然临床口子新宠，原来可以做这么多方向！ |
