从数学建模到药物发现：高维特征选择与机器学习模型实战解析-平芜编程栈

1. 项目概述：当数学建模遇上医学大数据

去年带队参加华为杯（现为“华为杯”中国研究生数学建模竞赛），我们组选的就是那道关于抗乳腺癌药物活性预测的C题。这道题在当时引起了不小的讨论，因为它完美地踩在了几个风口上：数学建模、大数据、机器学习，以及最前沿的生物医药计算。题目给了我们一堆化合物的分子描述符数据，目标是构建模型来预测这些化合物对雌激素受体α（ERα）的生物活性。ERα是治疗乳腺癌的一个重要靶点，能高效预测化合物活性，就意味着能大大加速新药研发的筛选过程，节省天文数字的研发成本。

这不仅仅是一道竞赛题，它几乎是当前“AI for Science”在药物发现领域的一个微型沙盘推演。你需要处理的是典型的高维、小样本生物医学数据（特征多，样本相对少），然后运用机器学习甚至深度学习的方法，从海量特征中挖掘出与生物活性相关的关键模式。对于参赛者来说，它考察的不仅仅是调包调用几个模型，更是对数据理解、特征工程、模型选择与解释、以及结果业务化落地的全链条思考。无论是数学、统计、计算机还是生物背景的同学，都能在这道题里找到发挥的空间，但也都会面临跨学科的挑战。接下来，我就结合我们当时的解题思路和后续的一些反思，拆解一下这道题的核心脉络与实操要点，希望能给未来参加类似赛题，或对“数据+AI”驱动药物研发感兴趣的朋友一些实在的参考。

2. 赛题核心与解题思路全拆解

2.1 问题本质：一个高维特征下的回归与分类混合预测任务

拿到题目和数据，第一步永远是穿透表象看本质。这道题提供的核心数据是一个包含1974个化合物样本的数据集，每个化合物用729个分子描述符（特征）进行表征。这些描述符可能包括分子的物理化学性质（如分子量、脂水分配系数logP）、拓扑结构指数、电子属性等，它们从不同维度数字化了一个分子的结构信息。我们的目标变量是这些化合物对ERα的pIC50值（生物活性指标，值越大通常表示活性越强）。

所以，从机器学习任务类型上看，首要任务是一个回归问题：用729维特征预测一个连续的pIC50值。但题目往往不止于此，它通常还会要求根据预测的活性进行化合物筛选（比如找出高活性的前N个化合物），这又引入了分类问题的思想（例如，设定一个pIC50阈值，如6.0或7.0，将化合物划分为“活性”与“非活性”）。因此，这是一个回归与分类思想交织的任务。

更关键的挑战在于数据特性：“维数灾难”。样本数（1974）远小于特征数（729），直接使用所有特征建模极易导致过拟合，模型会在训练集上表现完美，在未知数据上却一塌糊涂。因此，解题的核心主线必然围绕“特征工程”展开，目的是从729个特征中筛选出最相关、最不冗余、最具生物学解释性的一小部分子集。

2.2 解题总路线图：从数据清洗到模型集成

基于以上分析，一个稳健的解题流程可以概括为以下五个阶段，它构成了我们整个工作的骨架：

数据预处理与探索性分析：这是所有数据工作的基石。处理缺失值、异常值，进行初步的数据分布观察，了解特征与目标变量的基本关系。
特征工程与降维：这是本赛题的绝对核心。通过过滤法、包裹法、嵌入法等多种手段，大幅削减特征数量，保留精华。同时，可能需要进行特征构造或变换。
机器学习模型构建与训练：在降维后的特征子集上，选择合适的回归模型进行训练。这里会涉及模型选择、超参数调优、以及至关重要的交叉验证。
模型评估与解释：不仅仅看R²、RMSE等回归指标，还要从业务角度（药物筛选）评估模型。同时，利用SHAP、特征重要性等工具解释模型，让预测结果具有生物学意义。
结果整合与策略建议：根据模型预测结果，给出化合物筛选名单，并形成一套完整的计算机辅助药物设计流程建议。

这个路线图看似标准，但每一步在具体实施时，都有大量细节和技巧需要把握，尤其是在有限的竞赛时间内做出正确的权衡。

3. 特征工程：从729到核心特征的“瘦身”艺术

特征工程是本题成败的生命线。我们的目标是找到约20-50个核心特征，既能保证模型性能，又具备良好的解释性。我们采用了多轮次、多方法结合的“组合拳”策略。

3.1 第一轮：基于统计与相关性的快速过滤

首先进行粗筛，剔除明显无效的特征。

缺失值处理：检查每个特征的缺失值比例。对于缺失率过高的特征（例如>30%），直接剔除，因为填补大量缺失值会引入过多噪声。
方差过滤：使用VarianceThreshold。方差接近0的特征，意味着该特征在所有样本上基本取值相同，不携带信息，直接删除。
单变量相关性分析：计算每个特征与目标变量pIC50的相关系数（如皮尔逊相关系数）。剔除那些与目标变量绝对相关系数极低（例如<0.05）的特征。这一步可以快速去掉大量明显不相关的特征。

实操心得：相关系数阈值不宜设得太高，初期可以宽松一些（如0.03），目的是快速缩减特征规模到400-500左右，为后续更精细的方法减负。同时，要小心非线性关系，相关系数低不代表没关系，但作为第一轮过滤是高效的。

3.2 第二轮：处理多重共线性与包裹式选择

经过第一轮，特征数可能还剩400-500个，但特征之间可能存在高度的相关性（多重共线性），这会影响模型的稳定性和解释性。

相关性热图与聚类：计算剩余特征之间的相关系数矩阵，绘制热图。可以观察到哪些特征高度相关（例如相关系数>0.9）。对于高度相关的特征组，通常只保留其中一个（可以是与目标变量相关性最高的，也可以是业务上更易解释的）。
递归特征消除：这是一种强大的包裹式方法。我们选择了一个基础模型（如线性回归、随机森林），让RFE自动递归地剔除最不重要的特征，直到剩下指定数量的特征。RFE的结果依赖于选择的基础模型。

# 示例：使用线性回归作为基模型的RFE from sklearn.feature_selection import RFE from sklearn.linear_model import LinearRegression lr = LinearRegression() rfe = RFE(estimator=lr, n_features_to_select=50, step=10) # 目标选择50个特征，每次迭代剔除10个 rfe.fit(X_filtered, y) selected_features_rfe = X_filtered.columns[rfe.support_]

3.3 第三轮：基于树模型的嵌入法选择

嵌入法在模型训练过程中自动进行特征选择。树模型（如随机森林、XGBoost）能提供很好的特征重要性评分。

训练树模型并获取重要性：用全部或第二轮筛选后的特征训练一个随机森林回归模型。训练完成后，模型会输出每个特征的重要性得分（基于基尼不纯度减少或袋外误差）。
重要性排序与阈值选择：将特征按重要性降序排列。这里没有固定阈值，我们可以通过观察重要性得分的“拐点”（肘部法则）来决定保留前K个特征。也可以设定一个累积重要性贡献的阈值（如保留贡献了95%重要性的特征）。

from sklearn.ensemble import RandomForestRegressor import matplotlib.pyplot as plt rf = RandomForestRegressor(n_estimators=100, random_state=42, oob_score=True) rf.fit(X_train, y_train) importances = rf.feature_importances_ indices = np.argsort(importances)[::-1] # 绘制特征重要性排序图 plt.figure(figsize=(12,6)) plt.title("Feature Importances") plt.bar(range(X_train.shape[1]), importances[indices]) plt.xlabel("Feature Index") plt.ylabel("Importance") plt.show()

注意事项：随机森林的特征重要性倾向于偏向具有更多类别或数值范围更大的特征。对于高维数据，其重要性评估可能不够稳定。因此，最好将RFE和树模型重要性筛选的结果取交集或并集，再进行人工复审，这样得到的特征子集更可靠。我们最终通过三轮筛选，将特征数控制在了35个左右。

4. 模型构建、训练与超参数调优实战

特征准备好后，就进入模型构建阶段。我们面临多种模型选择，关键在于理解其特性并进行系统化比较。

4.1 模型选型：为什么是这些模型？

我们主要对比测试了以下几类模型，它们各有优劣：

线性模型：如岭回归、Lasso回归。优点是可解释性强，计算快。Lasso自带特征选择，可以进一步压缩特征。缺点是无法捕捉复杂的非线性关系。
支持向量机：特别是支持向量回归。在高维空间表现可能不错，但对参数和核函数选择敏感，训练速度相对慢。
树集成模型：如随机森林、梯度提升树。这是我们的重点候选。它们能自动处理非线性关系，对异常值不敏感，通常能取得较好的预测性能。XGBoost、LightGBM这类梯度提升框架更是竞赛常客。
多层感知机：即简单的神经网络。理论上可以拟合任何复杂关系，但在这种样本量不大、特征已降维的情况下，未必比树模型有优势，且调参更复杂，解释性差。

我们的策略是：先用线性模型和随机森林建立基线，再用梯度提升树进行精调。

4.2 交叉验证：防止过拟合的黄金准则

在竞赛和实际建模中，绝对不能只用一次训练集/测试集分割来评估模型。我们必须使用交叉验证。

我们选择了5折或10折交叉验证：将训练数据分成5或10份，轮流用其中4份或9份训练，1份验证，重复5或10次，取性能指标的平均值。这能更稳健地评估模型的泛化能力。
关键点：在进行任何基于验证集的操作（包括特征选择、超参数调优）时，都必须小心数据泄露。例如，特征缩放（StandardScaler）的fit操作只能在训练折上进行，然后transform训练折和验证折。我们使用Pipeline和GridSearchCV/RandomizedSearchCV来确保这个过程是干净的。

4.3 超参数调优：以LightGBM为例的实战

我们最终选择了LightGBM作为主力模型，因为它训练速度快，对类别特征友好，且性能强劲。以下是我们的调优步骤：

设定参数搜索空间：我们并不盲目网格搜索所有参数，而是有重点地调整。

param_grid = { 'learning_rate': [0.01, 0.05, 0.1], # 学习率，控制每棵树的影响力 'n_estimators': [100, 200, 500], # 树的数量 'max_depth': [3, 5, 7, -1], # 树的最大深度，-1表示不限制（需谨慎） 'num_leaves': [15, 31, 63], # 叶子节点数，与max_depth相关 'subsample': [0.8, 0.9, 1.0], # 样本采样比例 'colsample_bytree': [0.8, 0.9, 1.0], # 特征采样比例 'reg_alpha': [0, 0.1, 1], # L1正则化项 'reg_lambda': [0, 0.1, 1], # L2正则化项 'min_child_samples': [5, 10, 20] # 叶子节点最小样本数，防止过拟合 }

使用随机搜索：由于参数组合太多，网格搜索耗时太长。我们使用RandomizedSearchCV进行50-100轮的随机采样搜索，效率更高。

from sklearn.model_selection import RandomizedSearchCV import lightgbm as lgb lgb_model = lgb.LGBMRegressor(random_state=42, verbose=-1) random_search = RandomizedSearchCV( estimator=lgb_model, param_distributions=param_grid, n_iter=80, cv=5, scoring='neg_root_mean_squared_error', # 以负RMSE作为评分，越大越好 verbose=1, random_state=42, n_jobs=-1 ) random_search.fit(X_train_cv, y_train_cv)

锁定最佳参数并最终训练：随机搜索找到一组较优参数后，在其附近进行小范围的网格搜索微调，最终确定模型参数，并用全部训练数据重新训练最终模型。

踩坑实录：一开始我们贪图深度和叶子数，设置了很大的max_depth和num_leaves，结果模型在训练集上RMSE极低，但在交叉验证中波动很大，出现了明显的过拟合。后来通过加大reg_alpha、reg_lambda，限制max_depth在5-7之间，并增加min_child_samples，才使模型稳定下来。树模型调参的核心就是在偏差和方差之间做trade-off。

5. 模型评估、解释与业务化输出

模型训练好不是终点，如何评估其好坏，并让结果产生实际价值，才是关键。

5.1 多维度评估指标

我们不仅看一个R²，而是从多个角度评估：

回归指标：均方根误差、决定系数。这些是基础。
排序能力评估：由于最终目的是筛选高活性化合物，我们关心模型预测值的排序是否准确。可以计算预测值与真实值的斯皮尔曼等级相关系数。这个指标对我们来说甚至比RMSE更重要。
分类阈值评估：设定一个pIC50活性阈值（如6.3）。将模型预测值根据此阈值二分类，计算准确率、召回率、F1-score等。这直接模拟了药物筛选场景。

5.2 模型可解释性：SHAP值分析

“黑箱模型”在科研中是不被接受的。我们必须解释为什么模型认为某个化合物活性高。SHAP值是目前最强大的模型解释工具之一。

计算SHAP值：对测试集的样本，计算每个特征对该样本预测结果的贡献值。
可视化分析：
- 摘要图：可以看到哪些特征整体上最重要，以及特征值与SHAP值的关系（正向/负向影响）。
- 单个样本力图：可以针对某个具体的高活性预测化合物，可视化其各个特征是如何将预测值从基线（所有样本的平均预测）推高或拉低的。这能给出非常直观的“分子设计指导”，例如：“这个化合物活性高，主要是因为它的‘拓扑极性表面积’特征值很低，且‘碳原子数’特征值适中。”

import shap # 创建解释器 explainer = shap.TreeExplainer(best_lgb_model) shap_values = explainer.shap_values(X_test) # 绘制特征重要性摘要图 shap.summary_plot(shap_values, X_test, plot_type="dot")

5.3 结果输出与策略建议

最后，我们将所有流程整合，输出最终结果：

高活性化合物列表：使用最终模型对所有1974个化合物进行预测，按预测pIC50值降序排列，输出排名前50或100的化合物及其预测值、关键特征值。
关键分子描述符报告：结合特征重要性排序和SHAP分析，列出5-10个对ERα活性影响最大的分子描述符，并简要说明其物理化学意义（例如，MLogP表示预测的脂水分配系数，影响化合物透膜能力）。
虚拟筛选流程建议：在论文中，我们提出了一套完整的计算机辅助药物筛选流程建议：
- 初筛：使用我们构建的快速预测模型，对百万级虚拟化合物库进行第一轮粗筛，快速淘汰低活性化合物。
- 精筛：对粗筛出的几万个化合物，进行更精确的分子对接模拟或更复杂的QSAR模型预测。
- 实验验证：对精筛出的几十到几百个顶级候选化合物，进行湿实验验证。

6. 常见问题与避坑指南实录

在整个解题和后续复盘过程中，我们遇到了不少典型问题，这里总结出来，希望大家能避开这些坑。

6.1 数据预处理中的陷阱

缺失值处理不当：直接删除缺失值过多的特征是对的，但对于剩余特征的少量缺失值，采用中位数或众数填补是常用方法。但要注意，千万不要在划分训练集和测试集之前就对整个数据集进行填补！这会导致信息从“未来”（测试集）泄露到“过去”（训练集）。正确的做法是在交叉验证的每一折内，或者使用Pipeline时在训练折上fit填补器，再transform所有数据。
特征缩放的必要性：对于基于距离的模型，必须进行特征缩放。但对于树模型，理论上不需要。然而，如果使用了线性模型作为特征选择（如Lasso）或对比基线，或者后续要做特征重要性比较，进行标准化通常是好习惯。我们使用了StandardScaler，并在Pipeline中管理。

6.2 特征工程与模型训练中的误区

特征选择中的数据泄露：这是最致命的错误之一。如果在特征选择阶段（例如计算相关系数、进行RFE）时使用了全部数据（包括未来的测试集），那么选择出的特征已经“见过”测试集，评估结果会极度乐观，完全不具泛化性。必须将特征选择过程嵌入到交叉验证的循环内部，或者严格在训练集上进行特征选择，再将选择规则应用于测试集。Sklearn的SelectFromModel和RFE在Pipeline中与GridSearchCV结合可以很好地避免此问题。
盲目追求复杂模型：一开始我们尝试了复杂的深度学习模型，但效果并不比调优后的LightGBM好，且训练时间长，调参困难。在数据量不是特别巨大的情况下，梯度提升树模型往往是性价比最高的选择。
忽略模型集成：单一模型再好也有其局限性。我们后期尝试了将LightGBM、随机森林和岭回归的预测结果进行加权平均（Stacking的简单版），发现集成后的模型在交叉验证上的RMSE和稳定性均有小幅提升。这在竞赛最后阶段是提分的关键技巧。

6.3 比赛策略与时间管理

尽早确定基线：拿到数据后，用最简单的模型（如线性回归）和原始特征跑出一个基准分数。这个分数是所有改进的起点。
并行实验：特征工程和模型调优可以分头进行。一组同学专攻特征筛选和构造，另一组同学在筛选出的不同特征子集上测试不同模型和参数。
重视文档与可复现性：所有数据处理步骤、特征选择记录、模型参数、实验结果都必须即时记录。Jupyter Notebook的单元格编号有时会混乱，我们后来转向使用脚本配合配置文件的方式，确保每一步都可追溯、可复现。
论文写作与可视化同步：不要把所有分析做完再写论文。图表（如特征相关性热图、特征重要性图、SHAP图、模型性能对比图）在分析过程中就生成并保存好，并配上简要说明。这会让最后的论文撰写事半功倍。

这道华为杯的C题，是一个绝佳的将数学建模思想、大数据处理技术和机器学习算法应用于实际科学问题的案例。它教会我们的，远不止几个Sklearn的API调用，而是一套从问题定义、数据洞察、方法选择、实验验证到结果解释的完整数据科学工作流。在药物研发成本高企的今天，这种“干湿结合”的研究范式，正展现出越来越巨大的潜力。