1. 項(xiàng)目概述當(dāng)數(shù)學(xué)建模遇上醫(yī)學(xué)大數(shù)據(jù)去年帶隊(duì)參加華為杯現(xiàn)為“華為杯”中國研究生數(shù)學(xué)建模競賽我們組選的就是那道關(guān)于抗乳腺癌藥物活性預(yù)測的C題。這道題在當(dāng)時(shí)引起了不小的討論因?yàn)樗昝赖夭仍诹藥讉€(gè)風(fēng)口上數(shù)學(xué)建模、大數(shù)據(jù)、機(jī)器學(xué)習(xí)以及最前沿的生物醫(yī)藥計(jì)算。題目給了我們一堆化合物的分子描述符數(shù)據(jù)目標(biāo)是構(gòu)建模型來預(yù)測這些化合物對雌激素受體αERα的生物活性。ERα是治療乳腺癌的一個(gè)重要靶點(diǎn)能高效預(yù)測化合物活性就意味著能大大加速新藥研發(fā)的篩選過程節(jié)省天文數(shù)字的研發(fā)成本。這不僅僅是一道競賽題它幾乎是當(dāng)前“AI for Science”在藥物發(fā)現(xiàn)領(lǐng)域的一個(gè)微型沙盤推演。你需要處理的是典型的高維、小樣本生物醫(yī)學(xué)數(shù)據(jù)特征多樣本相對少然后運(yùn)用機(jī)器學(xué)習(xí)甚至深度學(xué)習(xí)的方法從海量特征中挖掘出與生物活性相關(guān)的關(guān)鍵模式。對于參賽者來說它考察的不僅僅是調(diào)包調(diào)用幾個(gè)模型更是對數(shù)據(jù)理解、特征工程、模型選擇與解釋、以及結(jié)果業(yè)務(wù)化落地的全鏈條思考。無論是數(shù)學(xué)、統(tǒng)計(jì)、計(jì)算機(jī)還是生物背景的同學(xué)都能在這道題里找到發(fā)揮的空間但也都會(huì)面臨跨學(xué)科的挑戰(zhàn)。接下來我就結(jié)合我們當(dāng)時(shí)的解題思路和后續(xù)的一些反思拆解一下這道題的核心脈絡(luò)與實(shí)操要點(diǎn)希望能給未來參加類似賽題或?qū)Α皵?shù)據(jù)AI”驅(qū)動(dòng)藥物研發(fā)感興趣的朋友一些實(shí)在的參考。2. 賽題核心與解題思路全拆解2.1 問題本質(zhì)一個(gè)高維特征下的回歸與分類混合預(yù)測任務(wù)拿到題目和數(shù)據(jù)第一步永遠(yuǎn)是穿透表象看本質(zhì)。這道題提供的核心數(shù)據(jù)是一個(gè)包含1974個(gè)化合物樣本的數(shù)據(jù)集每個(gè)化合物用729個(gè)分子描述符特征進(jìn)行表征。這些描述符可能包括分子的物理化學(xué)性質(zhì)如分子量、脂水分配系數(shù)logP、拓?fù)浣Y(jié)構(gòu)指數(shù)、電子屬性等它們從不同維度數(shù)字化了一個(gè)分子的結(jié)構(gòu)信息。我們的目標(biāo)變量是這些化合物對ERα的pIC50值生物活性指標(biāo)值越大通常表示活性越強(qiáng)。所以從機(jī)器學(xué)習(xí)任務(wù)類型上看首要任務(wù)是一個(gè)回歸問題用729維特征預(yù)測一個(gè)連續(xù)的pIC50值。但題目往往不止于此它通常還會(huì)要求根據(jù)預(yù)測的活性進(jìn)行化合物篩選比如找出高活性的前N個(gè)化合物這又引入了分類問題的思想例如設(shè)定一個(gè)pIC50閾值如6.0或7.0將化合物劃分為“活性”與“非活性”。因此這是一個(gè)回歸與分類思想交織的任務(wù)。更關(guān)鍵的挑戰(zhàn)在于數(shù)據(jù)特性“維數(shù)災(zāi)難”。樣本數(shù)1974遠(yuǎn)小于特征數(shù)729直接使用所有特征建模極易導(dǎo)致過擬合模型會(huì)在訓(xùn)練集上表現(xiàn)完美在未知數(shù)據(jù)上卻一塌糊涂。因此解題的核心主線必然圍繞“特征工程”展開目的是從729個(gè)特征中篩選出最相關(guān)、最不冗余、最具生物學(xué)解釋性的一小部分子集。2.2 解題總路線圖從數(shù)據(jù)清洗到模型集成基于以上分析一個(gè)穩(wěn)健的解題流程可以概括為以下五個(gè)階段它構(gòu)成了我們整個(gè)工作的骨架數(shù)據(jù)預(yù)處理與探索性分析這是所有數(shù)據(jù)工作的基石。處理缺失值、異常值進(jìn)行初步的數(shù)據(jù)分布觀察了解特征與目標(biāo)變量的基本關(guān)系。特征工程與降維這是本賽題的絕對核心。通過過濾法、包裹法、嵌入法等多種手段大幅削減特征數(shù)量保留精華。同時(shí)可能需要進(jìn)行特征構(gòu)造或變換。機(jī)器學(xué)習(xí)模型構(gòu)建與訓(xùn)練在降維后的特征子集上選擇合適的回歸模型進(jìn)行訓(xùn)練。這里會(huì)涉及模型選擇、超參數(shù)調(diào)優(yōu)、以及至關(guān)重要的交叉驗(yàn)證。模型評估與解釋不僅僅看R2、RMSE等回歸指標(biāo)還要從業(yè)務(wù)角度藥物篩選評估模型。同時(shí)利用SHAP、特征重要性等工具解釋模型讓預(yù)測結(jié)果具有生物學(xué)意義。結(jié)果整合與策略建議根據(jù)模型預(yù)測結(jié)果給出化合物篩選名單并形成一套完整的計(jì)算機(jī)輔助藥物設(shè)計(jì)流程建議。這個(gè)路線圖看似標(biāo)準(zhǔn)但每一步在具體實(shí)施時(shí)都有大量細(xì)節(jié)和技巧需要把握尤其是在有限的競賽時(shí)間內(nèi)做出正確的權(quán)衡。3. 特征工程從729到核心特征的“瘦身”藝術(shù)特征工程是本題成敗的生命線。我們的目標(biāo)是找到約20-50個(gè)核心特征既能保證模型性能又具備良好的解釋性。我們采用了多輪次、多方法結(jié)合的“組合拳”策略。3.1 第一輪基于統(tǒng)計(jì)與相關(guān)性的快速過濾首先進(jìn)行粗篩剔除明顯無效的特征。缺失值處理檢查每個(gè)特征的缺失值比例。對于缺失率過高的特征例如30%直接剔除因?yàn)樘钛a(bǔ)大量缺失值會(huì)引入過多噪聲。方差過濾使用VarianceThreshold。方差接近0的特征意味著該特征在所有樣本上基本取值相同不攜帶信息直接刪除。單變量相關(guān)性分析計(jì)算每個(gè)特征與目標(biāo)變量pIC50的相關(guān)系數(shù)如皮爾遜相關(guān)系數(shù)。剔除那些與目標(biāo)變量絕對相關(guān)系數(shù)極低例如0.05的特征。這一步可以快速去掉大量明顯不相關(guān)的特征。實(shí)操心得相關(guān)系數(shù)閾值不宜設(shè)得太高初期可以寬松一些如0.03目的是快速縮減特征規(guī)模到400-500左右為后續(xù)更精細(xì)的方法減負(fù)。同時(shí)要小心非線性關(guān)系相關(guān)系數(shù)低不代表沒關(guān)系但作為第一輪過濾是高效的。3.2 第二輪處理多重共線性與包裹式選擇經(jīng)過第一輪特征數(shù)可能還剩400-500個(gè)但特征之間可能存在高度的相關(guān)性多重共線性這會(huì)影響模型的穩(wěn)定性和解釋性。相關(guān)性熱圖與聚類計(jì)算剩余特征之間的相關(guān)系數(shù)矩陣?yán)L制熱圖?？梢杂^察到哪些特征高度相關(guān)例如相關(guān)系數(shù)0.9。對于高度相關(guān)的特征組通常只保留其中一個(gè)可以是與目標(biāo)變量相關(guān)性最高的也可以是業(yè)務(wù)上更易解釋的。遞歸特征消除這是一種強(qiáng)大的包裹式方法。我們選擇了一個(gè)基礎(chǔ)模型如線性回歸、隨機(jī)森林讓RFE自動(dòng)遞歸地剔除最不重要的特征直到剩下指定數(shù)量的特征。RFE的結(jié)果依賴于選擇的基礎(chǔ)模型。# 示例使用線性回歸作為基模型的RFE from sklearn.feature_selection import RFE from sklearn.linear_model import LinearRegression lr LinearRegression() rfe RFE(estimatorlr, n_features_to_select50, step10) # 目標(biāo)選擇50個(gè)特征每次迭代剔除10個(gè) rfe.fit(X_filtered, y) selected_features_rfe X_filtered.columns[rfe.support_]3.3 第三輪基于樹模型的嵌入法選擇嵌入法在模型訓(xùn)練過程中自動(dòng)進(jìn)行特征選擇。樹模型如隨機(jī)森林、XGBoost能提供很好的特征重要性評分。訓(xùn)練樹模型并獲取重要性用全部或第二輪篩選后的特征訓(xùn)練一個(gè)隨機(jī)森林回歸模型。訓(xùn)練完成后模型會(huì)輸出每個(gè)特征的重要性得分基于基尼不純度減少或袋外誤差。重要性排序與閾值選擇將特征按重要性降序排列。這里沒有固定閾值我們可以通過觀察重要性得分的“拐點(diǎn)”肘部法則來決定保留前K個(gè)特征。也可以設(shè)定一個(gè)累積重要性貢獻(xiàn)的閾值如保留貢獻(xiàn)了95%重要性的特征。from sklearn.ensemble import RandomForestRegressor import matplotlib.pyplot as plt rf RandomForestRegressor(n_estimators100, random_state42, oob_scoreTrue) 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()注意事項(xiàng)隨機(jī)森林的特征重要性傾向于偏向具有更多類別或數(shù)值范圍更大的特征。對于高維數(shù)據(jù)其重要性評估可能不夠穩(wěn)定。因此最好將RFE和樹模型重要性篩選的結(jié)果取交集或并集再進(jìn)行人工復(fù)審這樣得到的特征子集更可靠。我們最終通過三輪篩選將特征數(shù)控制在了35個(gè)左右。4. 模型構(gòu)建、訓(xùn)練與超參數(shù)調(diào)優(yōu)實(shí)戰(zhàn)特征準(zhǔn)備好后就進(jìn)入模型構(gòu)建階段。我們面臨多種模型選擇關(guān)鍵在于理解其特性并進(jìn)行系統(tǒng)化比較。4.1 模型選型為什么是這些模型我們主要對比測試了以下幾類模型它們各有優(yōu)劣線性模型如嶺回歸、Lasso回歸。優(yōu)點(diǎn)是可解釋性強(qiáng)計(jì)算快。Lasso自帶特征選擇可以進(jìn)一步壓縮特征。缺點(diǎn)是無法捕捉復(fù)雜的非線性關(guān)系。支持向量機(jī)特別是支持向量回歸。在高維空間表現(xiàn)可能不錯(cuò)但對參數(shù)和核函數(shù)選擇敏感訓(xùn)練速度相對慢。樹集成模型如隨機(jī)森林、梯度提升樹。這是我們的重點(diǎn)候選。它們能自動(dòng)處理非線性關(guān)系對異常值不敏感通常能取得較好的預(yù)測性能。XGBoost、LightGBM這類梯度提升框架更是競賽?？?。多層感知機(jī)即簡單的神經(jīng)網(wǎng)絡(luò)。理論上可以擬合任何復(fù)雜關(guān)系但在這種樣本量不大、特征已降維的情況下未必比樹模型有優(yōu)勢且調(diào)參更復(fù)雜解釋性差。我們的策略是先用線性模型和隨機(jī)森林建立基線再用梯度提升樹進(jìn)行精調(diào)。4.2 交叉驗(yàn)證防止過擬合的黃金準(zhǔn)則在競賽和實(shí)際建模中絕對不能只用一次訓(xùn)練集/測試集分割來評估模型。我們必須使用交叉驗(yàn)證。我們選擇了5折或10折交叉驗(yàn)證將訓(xùn)練數(shù)據(jù)分成5或10份輪流用其中4份或9份訓(xùn)練1份驗(yàn)證重復(fù)5或10次取性能指標(biāo)的平均值。這能更穩(wěn)健地評估模型的泛化能力。關(guān)鍵點(diǎn)在進(jìn)行任何基于驗(yàn)證集的操作包括特征選擇、超參數(shù)調(diào)優(yōu)時(shí)都必須小心數(shù)據(jù)泄露。例如特征縮放StandardScaler的fit操作只能在訓(xùn)練折上進(jìn)行然后transform訓(xùn)練折和驗(yàn)證折。我們使用Pipeline和GridSearchCV/RandomizedSearchCV來確保這個(gè)過程是干凈的。4.3 超參數(shù)調(diào)優(yōu)以LightGBM為例的實(shí)戰(zhàn)我們最終選擇了LightGBM作為主力模型因?yàn)樗?xùn)練速度快對類別特征友好且性能強(qiáng)勁。以下是我們的調(diào)優(yōu)步驟設(shè)定參數(shù)搜索空間我們并不盲目網(wǎng)格搜索所有參數(shù)而是有重點(diǎn)地調(diào)整。param_grid { learning_rate: [0.01, 0.05, 0.1], # 學(xué)習(xí)率控制每棵樹的影響力 n_estimators: [100, 200, 500], # 樹的數(shù)量 max_depth: [3, 5, 7, -1], # 樹的最大深度-1表示不限制需謹(jǐn)慎 num_leaves: [15, 31, 63], # 葉子節(jié)點(diǎn)數(shù)與max_depth相關(guān) subsample: [0.8, 0.9, 1.0], # 樣本采樣比例 colsample_bytree: [0.8, 0.9, 1.0], # 特征采樣比例 reg_alpha: [0, 0.1, 1], # L1正則化項(xiàng) reg_lambda: [0, 0.1, 1], # L2正則化項(xiàng) min_child_samples: [5, 10, 20] # 葉子節(jié)點(diǎn)最小樣本數(shù)防止過擬合 }使用隨機(jī)搜索由于參數(shù)組合太多網(wǎng)格搜索耗時(shí)太長。我們使用RandomizedSearchCV進(jìn)行50-100輪的隨機(jī)采樣搜索效率更高。from sklearn.model_selection import RandomizedSearchCV import lightgbm as lgb lgb_model lgb.LGBMRegressor(random_state42, verbose-1) random_search RandomizedSearchCV( estimatorlgb_model, param_distributionsparam_grid, n_iter80, cv5, scoringneg_root_mean_squared_error, # 以負(fù)RMSE作為評分越大越好 verbose1, random_state42, n_jobs-1 ) random_search.fit(X_train_cv, y_train_cv)鎖定最佳參數(shù)并最終訓(xùn)練隨機(jī)搜索找到一組較優(yōu)參數(shù)后在其附近進(jìn)行小范圍的網(wǎng)格搜索微調(diào)最終確定模型參數(shù)并用全部訓(xùn)練數(shù)據(jù)重新訓(xùn)練最終模型。踩坑實(shí)錄一開始我們貪圖深度和葉子數(shù)設(shè)置了很大的max_depth和num_leaves結(jié)果模型在訓(xùn)練集上RMSE極低但在交叉驗(yàn)證中波動(dòng)很大出現(xiàn)了明顯的過擬合。后來通過加大reg_alpha、reg_lambda限制max_depth在5-7之間并增加min_child_samples才使模型穩(wěn)定下來。樹模型調(diào)參的核心就是在偏差和方差之間做trade-off。5. 模型評估、解釋與業(yè)務(wù)化輸出模型訓(xùn)練好不是終點(diǎn)如何評估其好壞并讓結(jié)果產(chǎn)生實(shí)際價(jià)值才是關(guān)鍵。5.1 多維度評估指標(biāo)我們不僅看一個(gè)R2而是從多個(gè)角度評估回歸指標(biāo)均方根誤差、決定系數(shù)。這些是基礎(chǔ)。排序能力評估由于最終目的是篩選高活性化合物我們關(guān)心模型預(yù)測值的排序是否準(zhǔn)確?？梢杂?jì)算預(yù)測值與真實(shí)值的斯皮爾曼等級相關(guān)系數(shù)。這個(gè)指標(biāo)對我們來說甚至比RMSE更重要。分類閾值評估設(shè)定一個(gè)pIC50活性閾值如6.3。將模型預(yù)測值根據(jù)此閾值二分類計(jì)算準(zhǔn)確率、召回率、F1-score等。這直接模擬了藥物篩選場景。5.2 模型可解釋性SHAP值分析“黑箱模型”在科研中是不被接受的。我們必須解釋為什么模型認(rèn)為某個(gè)化合物活性高。SHAP值是目前最強(qiáng)大的模型解釋工具之一。計(jì)算SHAP值對測試集的樣本計(jì)算每個(gè)特征對該樣本預(yù)測結(jié)果的貢獻(xiàn)值?？梢暬治稣獔D可以看到哪些特征整體上最重要以及特征值與SHAP值的關(guān)系正向/負(fù)向影響。單個(gè)樣本力圖可以針對某個(gè)具體的高活性預(yù)測化合物可視化其各個(gè)特征是如何將預(yù)測值從基線所有樣本的平均預(yù)測推高或拉低的。這能給出非常直觀的“分子設(shè)計(jì)指導(dǎo)”例如“這個(gè)化合物活性高主要是因?yàn)樗摹負(fù)錁O性表面積’特征值很低且‘碳原子數(shù)’特征值適中?！眎mport shap # 創(chuàng)建解釋器 explainer shap.TreeExplainer(best_lgb_model) shap_values explainer.shap_values(X_test) # 繪制特征重要性摘要圖 shap.summary_plot(shap_values, X_test, plot_typedot)5.3 結(jié)果輸出與策略建議最后我們將所有流程整合輸出最終結(jié)果高活性化合物列表使用最終模型對所有1974個(gè)化合物進(jìn)行預(yù)測按預(yù)測pIC50值降序排列輸出排名前50或100的化合物及其預(yù)測值、關(guān)鍵特征值。關(guān)鍵分子描述符報(bào)告結(jié)合特征重要性排序和SHAP分析列出5-10個(gè)對ERα活性影響最大的分子描述符并簡要說明其物理化學(xué)意義例如MLogP表示預(yù)測的脂水分配系數(shù)影響化合物透膜能力。虛擬篩選流程建議在論文中我們提出了一套完整的計(jì)算機(jī)輔助藥物篩選流程建議初篩使用我們構(gòu)建的快速預(yù)測模型對百萬級虛擬化合物庫進(jìn)行第一輪粗篩快速淘汰低活性化合物。精篩對粗篩出的幾萬個(gè)化合物進(jìn)行更精確的分子對接模擬或更復(fù)雜的QSAR模型預(yù)測。實(shí)驗(yàn)驗(yàn)證對精篩出的幾十到幾百個(gè)頂級候選化合物進(jìn)行濕實(shí)驗(yàn)驗(yàn)證。6. 常見問題與避坑指南實(shí)錄在整個(gè)解題和后續(xù)復(fù)盤過程中我們遇到了不少典型問題這里總結(jié)出來希望大家能避開這些坑。6.1 數(shù)據(jù)預(yù)處理中的陷阱缺失值處理不當(dāng)直接刪除缺失值過多的特征是對的但對于剩余特征的少量缺失值采用中位數(shù)或眾數(shù)填補(bǔ)是常用方法。但要注意千萬不要在劃分訓(xùn)練集和測試集之前就對整個(gè)數(shù)據(jù)集進(jìn)行填補(bǔ)這會(huì)導(dǎo)致信息從“未來”測試集泄露到“過去”訓(xùn)練集。正確的做法是在交叉驗(yàn)證的每一折內(nèi)或者使用Pipeline時(shí)在訓(xùn)練折上fit填補(bǔ)器再transform所有數(shù)據(jù)。特征縮放的必要性對于基于距離的模型必須進(jìn)行特征縮放。但對于樹模型理論上不需要。然而如果使用了線性模型作為特征選擇如Lasso或?qū)Ρ然€或者后續(xù)要做特征重要性比較進(jìn)行標(biāo)準(zhǔn)化通常是好習(xí)慣。我們使用了StandardScaler并在Pipeline中管理。6.2 特征工程與模型訓(xùn)練中的誤區(qū)特征選擇中的數(shù)據(jù)泄露這是最致命的錯(cuò)誤之一。如果在特征選擇階段例如計(jì)算相關(guān)系數(shù)、進(jìn)行RFE時(shí)使用了全部數(shù)據(jù)包括未來的測試集那么選擇出的特征已經(jīng)“見過”測試集評估結(jié)果會(huì)極度樂觀完全不具泛化性。必須將特征選擇過程嵌入到交叉驗(yàn)證的循環(huán)內(nèi)部或者嚴(yán)格在訓(xùn)練集上進(jìn)行特征選擇再將選擇規(guī)則應(yīng)用于測試集。Sklearn的SelectFromModel和RFE在Pipeline中與GridSearchCV結(jié)合可以很好地避免此問題。盲目追求復(fù)雜模型一開始我們嘗試了復(fù)雜的深度學(xué)習(xí)模型但效果并不比調(diào)優(yōu)后的LightGBM好且訓(xùn)練時(shí)間長調(diào)參困難。在數(shù)據(jù)量不是特別巨大的情況下梯度提升樹模型往往是性價(jià)比最高的選擇。忽略模型集成單一模型再好也有其局限性。我們后期嘗試了將LightGBM、隨機(jī)森林和嶺回歸的預(yù)測結(jié)果進(jìn)行加權(quán)平均Stacking的簡單版發(fā)現(xiàn)集成后的模型在交叉驗(yàn)證上的RMSE和穩(wěn)定性均有小幅提升。這在競賽最后階段是提分的關(guān)鍵技巧。6.3 比賽策略與時(shí)間管理盡早確定基線拿到數(shù)據(jù)后用最簡單的模型如線性回歸和原始特征跑出一個(gè)基準(zhǔn)分?jǐn)?shù)。這個(gè)分?jǐn)?shù)是所有改進(jìn)的起點(diǎn)。并行實(shí)驗(yàn)特征工程和模型調(diào)優(yōu)可以分頭進(jìn)行。一組同學(xué)專攻特征篩選和構(gòu)造另一組同學(xué)在篩選出的不同特征子集上測試不同模型和參數(shù)。重視文檔與可復(fù)現(xiàn)性所有數(shù)據(jù)處理步驟、特征選擇記錄、模型參數(shù)、實(shí)驗(yàn)結(jié)果都必須即時(shí)記錄。Jupyter Notebook的單元格編號有時(shí)會(huì)混亂我們后來轉(zhuǎn)向使用腳本配合配置文件的方式確保每一步都可追溯、可復(fù)現(xiàn)。論文寫作與可視化同步不要把所有分析做完再寫論文。圖表如特征相關(guān)性熱圖、特征重要性圖、SHAP圖、模型性能對比圖在分析過程中就生成并保存好并配上簡要說明。這會(huì)讓最后的論文撰寫事半功倍。這道華為杯的C題是一個(gè)絕佳的將數(shù)學(xué)建模思想、大數(shù)據(jù)處理技術(shù)和機(jī)器學(xué)習(xí)算法應(yīng)用于實(shí)際科學(xué)問題的案例。它教會(huì)我們的遠(yuǎn)不止幾個(gè)Sklearn的API調(diào)用而是一套從問題定義、數(shù)據(jù)洞察、方法選擇、實(shí)驗(yàn)驗(yàn)證到結(jié)果解釋的完整數(shù)據(jù)科學(xué)工作流。在藥物研發(fā)成本高企的今天這種“干濕結(jié)合”的研究范式正展現(xiàn)出越來越巨大的潛力。