抗乳腺癌候选药物的优化建模
在研发治疗乳腺癌药物的过程中，能拮抗ERα活性的化合物是治疗乳腺癌的重要候选药物，同时也要考虑到化合物在人体内具备良好的药代动力学性质和安全性(ADMET 性质)，如果吸收性能、代谢速度、毒副作用等性质不佳，依然很难成为药物。本文对给定的1974个化合物的分子描述符、生物活性以及ADMET性质进行处理分析，探寻对生物活性有重要影响的分子描述符，构建化合物生物活性的定量预测模型和ADMET 性质的分类预测模型，并基于两个模型构建目标优化模型，找出具体的分子描述符范围。文章综合运用了内置随机森林重要性、基于排列的重要性、SHapley Additive exPlanation特征重要性排序、决策树、逻辑回归、Light Gradient Boosting Machine (LightGBM)、相关性分析、证据权重WOE 和信息值IV 筛选、过采样、机器学习、XGBoost分类、RandomForest分类、粒子群算法、主要目标法等经典机器学习算法和分析方法对相关问题进行量化分析和数学建模，使用了 Python、MATLAB等软件实现模型并得到问题答案。
针对问题一(变量选择) ，通过观察发现，数据中存在稀有变量和异常值的现象，首先，剔除270个稀有变量及使用均值法填充少量异常值，然后，通过 Spearman相关性分析发现729个分子描述符之间存在一些相关性很高的变量，对相关性较高的213个变量以及与生物活性相关性较低的50个变量进行初步筛选，最终得到196个分子描述符。随后，构建随机森林模型，利用基于内置随机森林重要性、基于排列的重要性、基于 SHapley Additive exPlanation 的重要性三种特征重要性计算方法进行特征重要性排序， 筛选出前20个对生物活性最具有显著影响的分子描述符, 发现MDEC-23、LipoaffinityIndex、maxHsOH在三种特征重要性计算算法下均排名前三，随后选择可解释性的SHAP 算法分析前20个分子描述符对生物活性的正负影响程度。
针对问题二(生物活性定值预测) ，首先，本文构建了基于决策树、逻辑回归、线性回归、Light Gradient Boosting Machine (LightGBM)等十二类算法的 ERα生物活性的定量回归预测模型，为了更好地筛选出重要的描述符，在次选择的变量为问题一中经过预处理和相关性筛选后的变量，随后采用MSE、RMSE等指标对各个模型的性能进行评估。结果发现，基于LightGBM算法的生物活性定值回归预测模型表现效果最好，MSE 值为最低0.4424。随后，本文计算基于LightGBM算法的分子描述符的SHAP 值，并与问题一得
到的前20个对生物活性最具有显著影响的分子描述符进行对比，选择出了交叉的15个分析描述符作为特征,对模型的参数(如 max depth和 num leaves) 进行调整, 得到性能最优的生物活性定值回归预测模型,最后,对文件“ERα activity. xlsx”的 test表中的50个化合物进行pIC50值预测, 并通过pIC50值计算对应的IC50值。
针对问题三(分类模型构建和预测)，首先，基于证据权重 WOE 和信息值Ⅳ方法，对影响化合物的 Caco-2、CYP3A4、hERG、HOB、MN五种ADMET性质的分子描述符进行筛选，确定用于预测不同ADMET 性质的变量类型。其次，对五种ADMET性质的分类变量数据分布进行分析，采用过采样方法均衡数据样本。再次，构建13种分类模型，包括11种机器学习模型和2种深度学习模型(LSTM、CNN)，通过对各个模型准确率、精度、召回率、F1值、ROC 曲线、AUC值及对数损失等指标的评价和比较，确定预测 Caco-2、CYP3A4、hERG 三种ADMET性质的最佳模型为XGBoost分类模型, 预测HOB、MN 两种ADMET 性质的最佳模型为随机森林(RandomForest)分类模型。然后，为了预测 test集合中的ADMET 性质， 一方面，评估了最佳模型的泛化能力， 结果表明各个模型泛化能力较强，在测试集上的学习能力分数都达到了0.9以上，其中，预测MN的RanomForest 分类模型和预测CYP3A4的XGBoost分类模型在测试集上学习能力分数超过了0.96；另一方面，为了提高模型预测的准确率，基于十折交叉验证方法，在训练集上对获得的最佳模型进行参数调优，以获得最优模型，并基于最优参数下的最优模型预测ADMET性质，结果表明在参数调优下,预测 Caco-2的XGBoost模型准确率最高可达到93.9%;预测CYP3A4的XGBoost模型准确率最高可达到96.8%; 预测hERG的XGBoost模型准确率最高可达到 92.6%; 预测 HOB 的 RandomForest 模型准确率最高可达到 92%; 预测 MN 的RandomForest 模型准确率最高可达到97.8%。最后，本文对 5个分类模型的特征重要性排序和重要变量的关系进行简要的可视化分析，结果表明不同分类模型的重要特征具有明显差异，不同重要变量的组合对分类结果的影响也有明显差异。
针对问题四(优化选取分子描述符并估计范围)，本文构建了以提高抑制Erα的生物活性和提高ADMET 性质为目标的多目标优化模型，根据前文进行数据清理得到的主要分子描述符，基于问题二、三所构造的回归模型和分类模型构造出符合题目要求的多目标优化模型，然后通过粒子群算法权衡两个目标函数之间的关系，求解 Pareto解集，之后再通过主要目标法将模型再次求解验证答案的有效性，并从中选出部分分子描述符展示目标函数最优的取值范围，并且通过散点图可视化部分分子描述符的重要影响和取值范围。
关键词：生物活性分析； 逻辑回归； 分类预测； 随机森林； XGBoost； 机器学习； 相关性分析； 粒子群算法； 多目标优化
一、 前言
1.1 问题背景
乳腺癌发病多见于女性且近年来发病率在年轻人群体中所占有的比率逐渐增大，致死率较高，对于乳腺癌的发病机制目前的研究并未完全透彻，但发现其与雌激素受体密切相关，实验表明ERα在乳腺癌治疗中起着十分重要的作用，因此在临床治疗时，通常选用能拮抗ERα活性的化合物作为治疗乳腺癌的候选药物，如他莫昔芬可以显著提高乳腺癌患者的生存几率。但是由于药物研发通常需要极高的时间、精力成本，为做好乳腺癌药物的研发工作, 通常会建立化合物的定量结构-活性关系 ( Quantitative Structure- Activity Relationship，QSAR)模型，即利用计算和统计的方法定量研究化合物的二维结构、三维结构等结构特征与化合物的药性、活性等生活效应特征之间的关系，该方法也是药物分子学设计中重要的方法基础和理论基础，能够指导药物分子结构的优化，通过该方法计算所得到的参数值可以指导生物活性预测、药物分子设计等[1]。
良好的乳腺癌生物活性是一个化合物成为候选药物必备的条件之一，除此之外，该化合物还需要具备稳定性、安全性、可代谢、低毒性等特征，以便能够适应人体内部化学生物反应，化合物对人体的药代动力学性质通常是由ADMET 性质描述( Absorption 吸收、 Distribution分布、 Metabolism代谢、 Excretion排泄、 Toxicity毒性),好的性质是该化合物成为候选药的另一个条件。
本题的命题宗旨是利用相关的数据建立影响生物活性的定量预测模型及构建5种ADMET 性质的分类预测模型，两个模型的精度和泛化能力越高越好，此外，还需要找出重要的分子描述符及其取值范围或值从而为化合物同时优化ERα拮抗剂的生物活性和ADMET 性质提供相应指导。
1.2 问题重述
基于前述研究背景, 该题目共提供了四个附件, 即“ERα activity. xlsx”、“ Molecular Descriptor. xlsx”、“ADMET. xlsx”和“分子描述符含义解释. xlsx”, 基于四个附件内容，拟要解决的研究问题如下：
(1) 寻找影响生物活性的重要变量：依据变量对生物活性影响的重要程度对729个分子描述符进行变量选择并进行排序，给出前20个对生物活性具有显著影响的分子描述符。
(2) 构建影响生物活性的定量预测模型：选择不超过20个分子描述符变量，构建化合物对ERα生物活性的定量预测回归模型,并预测“ERα activity. xlsx” test表中50个化合物的IC50值和对应的pIC50值。
(3) 构建5种 ADMET 性质的分类预测模型: 将“ Molecular Descriptor. xlsx”提供的729个分子描述符视为自变量,“ADMET. xlsx”的五种 ADMET 性质视为因变量,分别构建化合物的 Caco-2、CYP3A4、hERG、HOB、MN的二分类预测模型,并对“ADMET. xlsx”的 test表中的ADMET性质进行预测。
(4) 确定分子描述符及其取值范围来获得更好的生物活性和更好的ADMET性质：寻找并阐述化合物的哪些分子描述符，以及这些分子描述符在什么取值或者处于什么取值范围时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的ADMET 性质(给定的五个ADMET性质中，至少三个性质较好)。
1.3 基于 VOSviewer 的关键词共现时序分析
为了使研究方向更加精确，本文通过关键词共现来揭示有关药物的生物活性和ADMET 优化相关研究的关键内容，以CNKI数据库为数据源，通过主题词“抗肿瘤活性”和“ADMET”进行检索,设置检索时间为2006-2021年, 借助可视化分析工具VOSviewer 对文章的关键词进行共现与聚类，从关键词的角度对药物的生物活性和ADMET 优化的研究内容进行初探。VOSviewer作为知识图谱软件，能可视化揭示文献数据之间的关系。最终的关键词共现时序分析结果见图1.1和图1.2。关键词出现的频次越多，图中的圈越大。目前以及有不少有关抗肿瘤药物活性的研究，包括足校关系模型、有效成分提取、成分分析、分离鉴定、正交试验、响应面优化发等等，涉及到药物动力学、药代动力学、药效学的相关理论研究。与ADMET有关的研究则更多的应用于药物设计和药物发现，与分子对接、虚拟筛选(利用小分子化合物与药物靶标间的分子对接运算，快速地从大量分子中筛选出具有成药性的活性化合物) 等有关的研究较为丰富。
二、模型假设
2.1模型基本假设
根据化合药物研发情况和本题所给出的条件，本文作出如下假设：
(1) 所给数据中某一分子描述符在所有所给化合物中的值都为0时，即这一分子描述符对活性影响并不大。
(2) 除了所给的分子描述符之外，其他存在的因素对生物活性和ADMET性质的影响可以忽略。
(3) 分子描述符的所给值默认为连续变量。
2.2模型符号说明
本文所涉及的模型符号说明如表2.1所示
表2.1模型符号说明
| 序号 | 符号 | 说明 |
| 1 | ρ | Spearman 相关系数 |
| 2 | DR/₂ | 差异度 |
| 3 | V | 信息增益 |
| 4 | L | 损失函数 |
| 5 | Ω | 正则项 |
| 6 | Obj | 目标函数 |
| 7 | w | 权重函数 |
| 8 | σ | Sigmoid激活函数 |
| 9 | p₁ | 概率 |
| 10 | “ Caco-2” | Caco-2 性质对应的值 |
| 11 | z | 可行空间 |
| 12 | v₁ | 速度 |
| 13 | present | 当前位置 |
认为是冗余特征，从而影响模型的精度。
表 4.1部分分子描述符的统计信息列表
| 分子描述符 | mean | std | min | 0.25 | 0.5 | 0.75 | max |
| pIC50 | 6.586186 | 1.423052 | 2.456 | 5.38225 | 6.581 | 7.5685 | 10.337 |
| nAcid | 0.108409 | 0.3479 | 0 | 0 | 0 | 0 | 4 |
| ALogP | 1.110164 | 1.43425 | -23.105 | 0.3763 | 1.17095 | 1.9481 | 5.1817 |
| naAromAtom | 15.44681 | 5.155854 | 0 | 12 | 16 | 18 | 30 |
| nAromBond | 16.18946 | 5.635271 | 0 | 12 | 18 | 18 | 34 |
| nB | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| nC | 22.6079 | 6.631359 | 7 | 17 | 22 | 28 | 95 |
| nN | 1.508612 | 1.886457 | 0 | 1 | 1 | 2 | 46 |
| nS | 0.307497 | 0.562536 | 0 | 0 | 0 | 1 | 6 |
| nP | 0.001013 | 0.031822 | 0 | 0 | 0 | 0 | 1 |
| nCl nBr | 0.101317 0.061297 | 0.36283 0.254292 | 0 0 | 0 0 | 0 0 | 0 0 | 6 3 |
对照分子描述符含义解释掌握分子描述符的含义。查阅相关文献，为了后续数据的处理,对 Molecular Descriptor. xlsx中的数据进行简单处理。保留对数据分析处理和数据建模有价值的数据信息，剔除稀有变量270个，剔除的部分分子描述符如下表 4.2所示。
表 4.2部分筛选掉的描述符信息列表
| 编号 | 分子描述符 | 含义 |
| 1 | nBondsQ | 四重键数 |
| 2 | nssGeH2 | 原子级电子态计数：-GeH2- |
| 3 | nsssGeH | 原子级电子态计数：>GeH- |
| 4 | nssssGe | 原子级电子态计数：>Ge< |
| 5 | nHmisc | 原子级电子态计数:H bonded to |
| B, Si,P, Ge, As, Se, Sn or Pb | ||
| 6 | nsLi | 原子级电子态计数：-Li |
| 7 | nssBe | 原子级电子态计数：-Be- |
| 8 | nssssBem | 原子级电子态计数：>Be<-2 |
| 9 10 | nsBH2 nssSiH2 | 原子级电子态计数：-BH2 原子级电子态计数：-SiH2- |
4.2.1 去除异常值
本文采用基于箱线图法对异常值进行识别，以分子描述符ALogp2 为例，在第1563个样本中出现了异常值517.4294，如图4.2所示，对于这些异常值，使用该位点其它数值平均值代替。
ρ=∑x-πy-y∑ix1-x2∑iy1-y2 (4-1)
其中，x，y为两个待分析变量取值，x，y为两个变量的平均值。相关系数的绝对值越大，则代表两个变量之间的相关性越强。部分变量的相关性如图4.3所示，其中，对角线变量表示自相关，其他位置的变量表示互相之间的相关性，颜色越深，接近于红色表明相关性越强，反之变量之间的相关性越弱。相关性分析结果表明部分变量之间存在很强的相关性，而冗余特征会损害模型的可解释性，因此需要剔除相关性高的213个特征。
图4.3部分变量之间的相关性系数
通过相关性系数筛选出了相关性大于0.9的变量，其中部分变量之间的相关性结果如图4.4所示，从图中可以看出，这些变量之间呈现高度相关的特征。
4.3.2描述符与药物活性之间的相关性
为了证明分子描述符与药物活性之间具有相关性，对药物活性和分子描述符之间的相关性进行计算，筛选掉相关性较低的50个变量，剔除的部分分子描述符如下表 4.3所示。
表4.3描述符与药物活性之间的相关性系数
| 编号 | 分子描述符 | 含义 |
| 1 | VCH-3 | Valence chain, order 3 |
| 2 | VCH-4 | Valence chain, order 4 |
| 3 | VPC-6 | Valence path cluster, order 6 |
| 4 | SP-0 | Simple path, order 0 |
| 5 | SP-1 | Simple path, order 1 |
| 6 7 | SP-2 | Simple path, order 2 |
| VP-7 | Valence path, order 7 | |
| 8 | ndCH2 | Count of atom- type E- State: =CH2 |
| 9 10 | ntCH ndsCH | Count of atom- type E- State:#CH Count of atom- type E- State: =CH- |
4.4变量方法的构建与评估
4.4.1 变量方法调研
变量选择是指拟合线性和非线性的多元回归方程进行自变量选择、拟合判别函数的变量选择等。针对问题一的要求，需要对1974个化合物的729个分子描述符进行变量选择。目前主流的变量方法有子集选择、收缩方法、维数缩减。对这些变量选择方法的特点和优缺点进行调研，得到结果如下表4.4所示。
表4.4回归分析方法比较
| 方法描述 | 优缺点 | |
| 最优子集选择 | 将所有的特征组合进行建模， 然 | 适用于小数据量、简单的 |
| 后根据AIC、BIC 等准则选择最 | 关系， 但经常表现为高方 | |
| 优的模型， 如向前选择和向后选 | 差， 因此不容易降低全模 | |
| 择 | 型的预测误差 | |
| 压缩系数法 | 基于惩罚项变量选择方法，主要 | 用于强相关的自变量时， |
| 指岭回归和 Lasso回归 | 在进行参数估计时， 会导 | |
| 致解不可逆， 十分不稳定 | ||
| 降维法 | 将变量进行变换后的新变量进 | 对特征值的分解存在局限 |
| 行降维，主要是指主成分回归和 | 性 | |
| 偏最小二乘法 | ||
| 树结构的方法 | 如随机森林模型， 其本身可用于 | 可以用统一的方法处理数 |
| 预测的模型， 但在预测过程中， | 值型变量和分类型变量 | |
| 可以对变量重要性进行排序，然 | ||
| 后通过这种排序来进行变量筛 选 |
4.4.2基于随机森林的变量选择模型
本文采用 Gini指数作为分子描述符的重要性评判指标，针对一棵树中的每个节点k，计算其 Gini指数, 如公式(4-2) 所示。
Gk=2px1-pk
(4-2)
其中，ρₖ代表样本在节点k属于任何一类的概率估计值，一个节点的重要性程度由节点分裂前后 Gini指数的变化量来确定：
IΔk=Gk-Gk1-Gk2
(4-3)
其中，Gk₁和Gₖ₂分别表示节点k产生的子节点，针对每棵树进行递归，最终随机抽样样本和变量，产生包含T棵数的森林，如果变量X₁在第t棵树中出现N次，则变量X₁的重要性为：
Iu=1π∑t=1T∑j=1NIΔj
(4-4)
4.5特征重要性排序方法和结果
4.5.1 基于排列 ( Permutation- based) 的特征重要性
除了内置随机森林重要性( Gini指数) 计算，还可以采用基于排列的特征重要性计算公式。对于样本数据D中的特征j，随机打乱数据集D中的第j列取值，并将新的数据集
13
命名为Dₖj， 计算模型在数据集Dₖⱼ上的分数 skj， 则计算特征f₁的重要性I，如以公式 (4-5)所示：
Ij=s-1x∑k=1ksk,j
4.5.2 基于 SHAP 的药物优化特征重要性排序
特征重要性是一种为预测模型的输入特征进行评分的方法，可以揭示进行预测时每个特征的相对重要性。SHapley Additive exPlanation(SHAP) 属于模型事后解释的方法,可以增强复杂机器学习模型的可解释性。因此，可以采用SHAP算法对模型中对药物活性优化的影响因素进行解释分析。
SHAP 值来源于博弈论中的 Shapley value，主要是用于评估每个特征对模型预测的贡献值，其基本原理是计算每个特征对模型的边际贡献，然后计算该特征在所有特征序列中不同的边界贡献，最后该特征所有边际贡献的均值即为SHAP 值。
假设模型基准分(所有样本的目标变量的均值)为ybase，第i个样本为x₁，第i个个样本的第j个特征为xᵢⱼ，特征的边际共现为msᵢⱼ，边的权重为wₖ，模型对该样本的预测值为y₂，则第i个样本的第1个特征的SHAP值f(x₁)如(4-6)所示,同时SHAP 值要服从公式(4-7)。
fxij=∑k=1nmsi1wk
(4-6)
yi=ybase+∑s=1nfxis
(4-7)
4.5.3特征重要性排序结果
本文对比分析三种不同的特征重要性算法的排序结果，如表4.5所示，输入模型的变量为经过数据预处理得到的196个分子描述符，从表中发现，分子描述符MDEC-23、LipoaffinityIndex、maxHsOH在三种特征重要性计算算法下均排名前三, 尽管之后的分子描述符的特征重要性排序有所不同，但对定值预测模型的影响较小。
表4.5三种特征重要性排序算法的结果对比
| 分子描述符 | 特征重要 性(GINI) | 分子描述符 | 特征重要性 ( Permutation ) | 分子描述符 | 特征重要性 (SHAP) |
| MDEC-23 | 0.171 | MDEC-23 | 0.126 | MDEC-23 | 0.31269807 5 |
| LipoaffinityInde | 0.077 | LipoaffinityInde | 0.05 | LipoaffinityInde | 0.17819301 |
| x | x | x | 3 | ||
| maxHsOH | 0.05 | maxHsOH | 0.04 | maxHsOH | 0.11601104 2 |
| minsssN | 0.044 | minHsOH | 0.024 | BCUTc-11 | 0.09790052 |
| minHBint5 | 0.038 | ClSP2 | 0.021 | ClSP2 | 0.093620064 |
| BCUTc-11 | 0.038 | VC-5 | 0.017 | minHBint5 | 0.082440564 |
| ClSP2 | 0.032 | BCUTc-11 | 0.017 | minsssN | 0.073858329 |
| minHsOH nHBAcc | 0.03 0.02 | nHBAcc minHBint5 | 0.016 0.015 | minHsOH nHBAcc | 0.070860766 0.066484049 |
| minsOH | 0.014 | ATSc3 | 0.01 | VC-5 | 0.037619015 |
| ATSc3 | 0.014 | minsssN | 0.008 | ATSc3 | 0.031815638 |
| VC-5 | 0.013 | minsOH | 0.007 | minsOH | 0.026887697 |
| SHBint10 | 0.01 | TopoPSA | 0.005 | XLogP | 0.023957633 |
| MLFER_A | 0.01 | MLFER_A | 0.005 | MLFER_A | 0.023897194 |
| XLogP | 0.009 | MDEC-33 | 0.004 | ndssC | 0.02354074 |
| maxHBa | 0.009 | VCH-5 | 0.004 | CrippenLogP | 0.023160492 |
| MDEO-12 | 0.008 | MDEC-22 | 0.004 | MDEO-12 | 0.019767238 |
| hmin | 0.008 | CrippenLogP | 0.003 | mindssC | 0.019740989 |
| CrippenLogP | 0.008 | MDEO-12 | 0.003 | TopoPSA | 0.018172734 |
| TopoPSA | 0.008 | gmin | 0.003 | maxHBint5 | 0.01770153 |
由于传统的特征重要性排序无法反映出特征如何影响预测结果，而SHAP能反映出每一个样本中的特征的影响力，同时能展示出影响的正负性。从随机森林模型和SHapley Additive exPlanation算法中得到的特征重要性排序结果如图4.5所示, 从图中可以发现,MDEC-23对药物的生物活性影响最大，随着MDEC-23(所有二级和三级碳之间的分子距离边缘)、LipoaffinityIndex(脂亲和指数)、maxHsOH(最大原子类型HE态:-OH)、minHBint5(路径长度为5的潜在氢键强度的最小电子态描述符)、minsssN(最小原子类型电子态：>N-)等分子描述符值的增加，生物活性越大的可能性增加。但BCUTc-11( nhigh最低部分电荷加权BCUT)、C1SP2(与另一碳结合的双重追踪碳)、nHBAcc(氢键受体数量(使用CDK HbondAcceptor计数描述符算法) ) 等对生物活性都有负面影响。
4.5模型小结
针对药物的生物活性影响因素中的变量选择问题，本文首先对变量进行了预处理和相关性分析，通过分析分子描述符之间相关性，发现绝大多数变量不符合正态分布规律，因此采用了 Spearman相关系数进行相关性分析，初步筛选掉冗余变量和与生物活性相关性较低的变量。随后，采用基于随机森林模型的变量选择算法，结合内置随机森林重要性、基于排列的重要性、SHapley Additive exPlanation三种特征重要性计算方法进行特征重要性
排序, 最终采用SHapley Additive exPlanation的特征重要性结果, 得到20个对生物活性最具有显著影响的分子描述符分别为: MDEC-23、LipoaffinityIndex、maxHsOH、BCUTe-11、C1SP2、 minHBint5、 minsssN、 minHsOH、 nHBAcc、 VC-5、 ATSc3、 minsOH、 XLogP、MLFER A、 ndssC、 CrippenLogP、 MDEO-12、 mindssC、 TopoPSA、 maxHBint5, 同时分析了不同描述符对生物活性影响的正负影响程度。
LightGBM 是基于 Histogram的决策树算法，其原理是将连续的浮点型特征离散成k 个离散值，并构造宽度为k的直方图，遍历训练数据后，统计每个离散值在直方图的累计统计量，而在进行特征选择时，只需要根据直方图的离散值，遍历寻找最优的分割点，如图5.1所示。
相比传统的XGBoost、AdaBoost自适应提升模型, LightGBM 采用带有深度限制的按叶子生长算法( leaf- wise) 代替了按层生成的决策树生长策略，可以降低更多误差，如图5.2所示。
提升树是利用加模型与前向分布算法实现学习的优化过程，LightGBM算法能过很好的通过减少特征量而不影响精确度的问题，它主要包含两个算法单边梯度采样( Gradient- based One- Side Sampling) 和互斥特征绑定算法( Exclusive Feature Bundling) 。
②单边梯度采样
对于n个样本的训练集 x1x2x3⋯xn,每个样本x₁有m维特征，模型的每次梯度迭代，变量的损失函数的负梯度方向为 g1g2g3gn,决策树通过最大信息增益点将数据分配到各个节点，0表示某个固定节点的训练集，分割特征j的分割点p定义为：
nr|0id=∑Ixi∈O;x1d] (5-4)
遍历每个特征的分裂节点，找到 di*=argmaxaVjd并且计算最大的信息增益值，然后根据数据特征j'的分裂节点将数据分到左右子节点。在GOSS中，根据数据的梯度降序排列训练集，保留最高的a个数据示例，作为数据子集A，对于剩下的样本进行随机采样获得子样本B，最后，信息增益的公式如下：
③互斥特征绑定算法
在抗乳腺癌候选药物优化中，药物数据和描述符数据往往具有特征量多并且特征空间系数的特点，互斥特征绑定算法可以通过特征捆绑的方式减少特征维度，对于互斥特征，LightGBM 使用直方图算法对这些特征进行合并。在两个特征并不是完全互斥的情况下，可以用一个指标对特征不互斥的程度进行衡量，从而得到冲突比率，当冲突比率较小时，可在不影响最后精度的前提下，将不完全互斥的特征进行捆绑。
5.1.3多种回归模型的评估
采用平均绝对误差(MAE)、均方误差(MSE)、均方根误差(RMSE)和R2(R- Square)决定系数来确定最优回归模型，计算公式分别如X所示。
MAE=1N∑i=1N|yi-yi| (5-6)
MSE=1N∑i=1Nyi-yi2 (5-7)
RMSE=1N∑l=1Nyi-yi2 (5-8)
R2=1-∑13y-γ12∑i=13y2-yi2 (5-9)
12个回归模型的评估结果如图5.3所示，LightGBM 在四个指标上的表现均最好，其次是随机森林模型，ExtraTree模型则表现最差。
在LightGBM模型下，基于1974个样本数据的196个变量，绘制如图5.4所示的原值和预测值的对比情况。其中纵坐标对应预测结果，横坐标对应真实值，通过真实值和预测值拟合一条直线，可以发现散点
代码
| # coding= utf-8 import xgboost from numpy import loadtxt from xgboo st import XGBRegressor from sklearn. model _ selection import train _ test _ split import pandas as pd import matplotlib. pyplot as plt import numpy as np from sklearn. model _ selection import GridSearchCV,KFold import shap import matplotlib from sklearn. metrics import mean _ squared _ error, mean _ absolute _ error,r2 _ score from math import sqrt from sklearn import linear _ model, tree, sym, neighbors, ensemble, ExtraTreeRegressor,BaggingRegressor, Lasso from sklearn. neural _ network import MLPRegressor from sklearn. inspection import permutation _ importance # Linear Regression model _LinearRegression= linear _ model. LinearRegression() # Decision Tree Regressor model _DecisionTreeRegressor= tree. DecisionTreeRegressor() #SVM Regressor model _SVR= svm. SVR( cache _ size=300) #K Neighbors Regressor model _KNeighborsRegressor= neighbors. KNeighborsRegressor() # Random Forest Regressor model _RandomForestRegressor ensemble. RandomForestRegressor(n_ estimators=200, random _ state=1) # Adaboost Regressor model _AdaBoostRegressor= ensemble. AdaBoostRegressor(n_ estimators=50) # Gradient Boosting Random Forest Regressor model _GradientBoostingRegressor= ensemble. GradientBoostingRegressor(n_ estimators=150) # bagging Regressor model _BaggingRegressor=BaggingRegressor() #ExtraTree Regressor model _ExtraTreeRegressor= ExtraTreeRegressor() # 载入数据集 df = pd. read _ csv('D:/2021 年 中 国 研 究 生 数 学 建 模 竞 赛 赛 题 / Activity _ train _I_ nounique _nocorr3. csv', encoding='gb18030', header=0) X= df. drop(["pIC50"], axis=1)# df. drop([" highvalue"], axis=1) Y = df["pIC50"] |