制剂 CQA 预测模型开发教程12非线性建模——随机森林、梯度提升与核方法版本声明块工具/软件Python 3.10.11 scikit-learn 1.7.2 numpy 2.2.6 scipy 1.15.3数据NIR Shootout 2002 片剂数据集655 片 × 650 波长本文目标能判断什么时候该换掉 PLS并且不会被树模型漂亮的 R²cal 骗到。一句话结论在 NIR Shootout 2002 的 assay 数据SNV 预处理上SVRrbf、C100、ε1.0是唯一在测试集上超过偏最小二乘回归PLS的模型R²test 0.9144 vs 0.8979、RPD 3.42 vs 3.13梯度提升的 R²cal 0.9994、RPDcal 39.85 是典型过拟合测试集仅 0.8697而随机森林在未预处理的原始光谱上退化到 R²test 0.5455、RPD 1.48落入 RPD1.5 的不可用区。〇、本篇要解决的认知问题光谱与 CQA 的关系一定是线性的吗PLS 的线性假设什么时候会失效树模型随机森林、梯度提升为什么在高维共线光谱上不一定比 PLS 强为什么梯度提升的 R²cal 能到 0.9994测试集却只有 0.8697SVR 为什么能在小样本校正集仅 155 片上超过 PLS什么情况下才值得放弃 PLS 去换非线性模型一、机制解析1.1 树模型在高维共线光谱上的三个固有弱点随机森林与梯度提升都以决策树为基学习器而决策树的分裂方式是选一个特征、找一个阈值、把样本切开——这种轴平行切分axis-aligned split在 NIR 场景下有三个天生的短板弱点一无法表达多个波长的线性组合这种斜超平面。PLS 的每一个潜变量都是 650 个波长的加权和天然是一条斜线方向的投影而单棵树只能用波长 1218 nm 某值这样的方盒子去逼近斜方向。要逼近一条斜线树必须切出大量小方块代价是深度与样本量。弱点二共线特征互相抢分裂机会重要性被稀释。NIR 相邻波长的相关系数极高1208–1236 nm 一带尤其如此见第 09 篇的 VIP 表。在 650 个高度相关的候选里挑分裂点等价于在大量看起来偶然更好的切分中取最大值——这是典型的选择偏差放大。同时真正有意义的信息带被摊薄到几十个几乎等价的波长上单棵树只随机用到其中几个森林平均后信号被稀释。弱点二的定量解释箱型划分如何把相邻波长的信息打散。一次分裂只用一个波长、一个阈值把样本空间切成两个轴平行的箱子而 NIR 吸收峰是横跨多个波长点的连续形状。本系列实测的 VIP 中最靠前的 15 个高权波长连续排布于 1208、1210 … 1236 nm步长正是数据集的 2 nmVIP 从 2.162 递减到 1.614——它们是同一份 C-H 键二级倍频信息的 15 次重复采样彼此高度共线。对 PLS 而言这种重复是资产一个潜变量把 650 个波长按载荷一次性线性组合15 个相邻点的信息被打包进同一方向不额外消耗自由度这正是SNVSG(21,2,d2)只用 3 个潜变量拿到 R²test 0.9200 的原因。对树而言它是负债某个节点一旦选中 1218 nm1216 nm 与 1220 nm 在该节点上就再没有发言权要补齐峰形只能靠后续节点或其它树再切而每多切一层样本就被再切一次——树深 d 的理论叶子数为 2 的 d 次幂落到叶子的平均样本数随之衰减。只有 155 个校正样本时叶子很快变稀疏模型只能靠记住样本来填满它这正是 1.3 节 GradientBoosting 出现 R²cal 0.9994 却只有 R²test 0.8697 的由来。自助采样bootstrap让每棵树只见约六成唯一样本进一步压薄了样本300 棵树平均后相邻波长的贡献被摊平信号被稀释、噪声被保留。弱点三树的预测值无法外推。树的输出是叶子内样本 y 的平均因此预测值永远被夹在训练集 y 的取值范围内是分段常数。如果测试集里有参考值超出手工训练范围外推的样本树模型在结构上就给不出合理答案。一句总结PLS 用全局线性组合换来了对高维共线的稳定性树模型用局部方盒子换来了非线性拟合能力——在 155 个校正样本、650 个共线波长的场景里前者的交易更划算。1.2 核方法与集成方法的适用边界模型族代表核心机制优势场景主要风险潜变量线性PLSNIPALS 迭代X/Y 协方差最大化小样本、高维、强共线本数据集无法表达非线性核方法SVRrbf核映射到高维 ε 不敏感损失 支持向量稀疏化小样本非线性对高维不敏感超参数C、ε、γ敏感无不确定度输出核方法贝叶斯GPR高斯过程先验直接输出预测方差需要不确定度、样本极少O(n³) 复杂度核函数选择难Bagging 集成RandomForest多棵去相关树平均降方差特征交互多、样本充足分段常数、无法外推、高维共线时退化Boosting 集成GradientBoosting逐轮拟合残差降偏差表格数据、中等维度极易过拟合对噪声敏感Boosting加速HistGradientBoosting直方图分箱加速分裂点搜索样本量大万级仍属 boosting同样需要早停适用边界一句话变量维度高、样本少、共线强时首选 PLS确认存在非线性且样本量还能支撑时先试 SVR只有样本量足够经验上校正集 200且需要建模复杂交互时树模型才值得上场。1.3 六组模型的真实对比以下是同一份 assay 数据SNV 预处理上的实测结果全部来自本系列的真实运行未经修饰模型R²calRPDcalR²testRMSEPRPDtestPLSnLV50.96165.120.89795.03003.13SVRrbf, C100, ε1.00.94264.190.91444.60743.42RandomForest300 树0.98608.470.85915.91082.67GradientBoosting0.999439.850.86975.68402.77HistGradientBoosting0.997519.910.86065.87832.68RandomForest原始光谱未预处理0.96775.580.545510.61551.48口径说明本表 PLS 基线为 R²test 0.8979 / RMSEP 5.0300 / RPDtest 3.13与第 11 篇锚点值 0.8981 / 5.0272 / 3.14 在小数第 4 位略有差异属不同对比实验流程的正常波动全系列以第 11 篇的 0.8981 / 5.0272 / 3.14 为锚点。三条必须写出来的关键结论① 梯度提升的 R²cal 0.9994、RPDcal 39.85 是典型过拟合不是本事。R²cal 0.9994 意味着模型几乎完美复现了校正集的每一个点——这通常不该发生因为 NIR 光谱本身含有噪声。它的测试集 R² 只有 0.8697低于 PLS 的 0.8979。这类训练集完美、测试集反而更差的模式是 boosting 在小样本高维数据上的标准失败姿势每轮拟合残差会把噪声也当成信号学进去。判据很简单——R²cal 与 R²test 的落差0.9994 − 0.8697 0.1297远大于 PLS 的落差0.9616 − 0.8979 0.0637。HistGradientBoosting 同理0.9975 → 0.8606。② SVR 是唯一在测试集上超过 PLS 的模型。SVRrbf 核取得 R²test 0.9144、RMSEP 4.6074、RPD 3.42三项全面优于 PLS 的 0.8979 / 5.0300 / 3.13而它的 R²cal 0.9426 还低于PLS 的 0.9616——这正是健康模型的样子校正集不追求满分泛化不退化。原因在于 rbf 核能在 155 个样本上用少量支持向量刻画非线性流形且 ε 不敏感带天然忽略了部分噪声。③ 随机森林在未预处理的原始光谱上大幅退化说明预处理与模型同等重要。同样是 RandomForest300 棵树在 SNV 预处理后 R²test 0.8591、RPD 2.67换成未预处理的原始光谱R²test 直接掉到 0.5455、RPD 1.48——跌入 RPD 1.5 的不可用区第 10 篇的分级。原因有两层树模型的分裂阈值直接建立在原始反射率上散射与基线漂移被当成真实信号去拟合同时未标准化让不同波段的量级差异主导了分裂点选择。结论换模型之前先换预处理。可选提及本环境无法安装故不给任何性能数字LightGBM与XGBoost属于同一族梯度提升实现直方图分箱 正则化 并行化在工程上常比 scikit-learn 的GradientBoostingRegressor更快。但本系列写作环境实测无法访问 PyPI因此不为它们提供任何可运行代码与性能数字——任何声称LightGBM 在本数据上 R²0.9x的说法在本系列中都不成立也不应被引用。1.4 何时值得换模型用验证信号决定的判定流程要不要换模型必须由校正集内部可得的信号回答而不是拿测试集反复试错铁律 2/5。按下面顺序逐级排查任一步为否就先解决该步不要跳到换模型① 预处理到位了吗 判据同一模型300 树 RandomForest在原始光谱与 SNV 光谱上的测试集 RPD —— 1.48 对 2.67 跨越了 RPD1.5 的不可用线。 否 → 先修预处理第 03、04 篇不要换模型。 │是 ▼ ② 线性模型是否已经把线性结构榨干 判据RMSECV 曲线是否在极小值后立刻回升第 08 篇的肘点。assay SNV 实测在 nLV4 取极小RMSECV 4.7564、Q² 0.9529后单调上升说明线性模型已经够用 若曲线一路平到底、没有肘点才说明线性假设不足、值得考虑核方法。 │是 ▼ ③ 残差还有系统结构吗 判据校正集内部交叉验证的残差对预测值作散点是否出现弯曲、漏斗形或分段偏移。健康模型的 残差应围绕 0 呈带状SVR 的 R²cal 0.9426 低于 PLS 的 0.9616但 R²test 反高0.9144 对 0.8979。 │有 ▼ ④ 过硬度门槛R²test 至少提升 0.01 且 RPD 至少提升 0.1否则不值得承担额外复杂度与验证成本。 实测 SVR 对 PLS0.8979 升到 0.9144、3.13 升到 3.42过门槛RandomForest0.8591与 GradientBoosting0.8697均低于 PLS 的 0.8979直接淘汰。 │过 ▼ ⑤ 换完必须补三件事重算适用域第 15 篇h* 里的 p 随模型而变、重出模型卡与版本号第 19 篇、 重跑一次预处理对照以确认增益不是预处理的功劳。第 ②、③ 步只花校正集就能算故须排在 ④ 之前一上来就跑测试集看谁高等于提前用掉它唯一的一次机会。二、完整代码与逐行剖析2.1 公平对比脚本第 12 篇PLS / SVR / 树模型在 assay 上的公平对比同一预处理、同一划分、同一指标importnumpyasnpimportpandasaspdfromscipy.ioimportloadmatfromsklearn.cross_decompositionimportPLSRegressionfromsklearn.svmimportSVRfromsklearn.ensembleimport(RandomForestRegressor,GradientBoostingRegressor,HistGradientBoostingRegressor)fromsklearn.metricsimportr2_score,root_mean_squared_errordefsnv(X):逐样本标准正态变量变换铁律 10不许用全局均值/方差做标准化。return(X-X.mean(axis1,keepdimsTrue))/X.std(axis1,ddof1,keepdimsTrue)defevaluate(name,model,Xcal,ycal,Xtest,ytest):统一的评价入口所有模型走同一套指标避免口径不一致导致的假比较。model.fit(Xcal,ycal)# 全部用默认随机种子之外的显式设置yp_calmodel.predict(Xcal).ravel()# ravel树模型的输出恒为 (n,)统一形状yp_testmodel.predict(Xtest).ravel()rmse_croot_mean_squared_error(ycal,yp_cal)# 铁律 1rmse_proot_mean_squared_error(ytest,yp_test)return{model:name,R2cal:r2_score(ycal,yp_cal),RPDcal:np.std(ycal,ddof1)/rmse_c,# 校正集 SD 口径R2test:r2_score(ytest,yp_test),RMSEP:rmse_p,RPDtest:np.std(ytest,ddof1)/rmse_p,# 测试集自身 SD 口径}if__name____main__:mloadmat(nir_shootout_2002.mat)uplambdak:np.asarray(m[k][data][0,0],dtypenp.float64)# 解包 转 float64Xcal,Xtestsnv(up(calibrate_1)),snv(up(test_1))ycal,ytestup(calibrate_Y)[:,2],up(test_Y)[:,2]# 第 3 列 assaymodels[(PLS (nLV5),PLSRegression(n_components5,scaleTrue)),(SVR (rbf),SVR(kernelrbf,C100.0,epsilon1.0)),(RandomForest,RandomForestRegressor(n_estimators300,random_state42,n_jobs-1)),(GradientBoosting,GradientBoostingRegressor(random_state42)),(HistGradientBoosting,HistGradientBoostingRegressor(random_state42)),]rows[evaluate(n,mo,Xcal,ycal,Xtest,ytest)forn,moinmodels]print(pd.DataFrame(rows).round(4).to_string(indexFalse))2.2 对照组同一模型、只换预处理# 对照组RandomForest 在原始光谱与SNV 后光谱上的差异raw_cal,raw_testup(calibrate_1),up(test_1)rflambda:RandomForestRegressor(n_estimators300,random_state42,n_jobs-1)fortag,Xc,Xtin[(SNV 预处理,Xcal,Xtest),(原始光谱,raw_cal,raw_test)]:revaluate(fRF {tag},rf(),Xc,ycal,Xt,ytest)print(f{r[model]:16s}R2test{r[R2test]:.4f}RMSEP{r[RMSEP]:.4f}RPD{r[RPDtest]:.2f})实测输出SNV 预处理下行 R²test 0.8591 / RMSEP 5.9108 / RPD 2.67原始光谱下行 R²test 0.5455 / RMSEP 10.6155 / RPD 1.48。同一个模型、同一批数据只因为少了 SNVRPD 从 2.67 掉到 1.48跨越了可定量与不可用的分界线。2.3 用排列重要性看树模型学到了哪些波长fromsklearn.inspectionimportpermutation_importancefromsklearn.model_selectionimporttrain_test_split# 在校正集内部再切一小块做重要性评估避免用测试集铁律 2Xa,Xb,ya,ybtrain_test_split(Xcal,ycal,test_size0.3,random_state42)rf_modelRandomForestRegressor(n_estimators300,random_state42,n_jobs-1).fit(Xa,ya)imppermutation_importance(rf_model,Xb,yb,n_repeats10,random_state42,scoringneg_root_mean_squared_error)waveup(axisscale).ravel()topnp.argsort(imp.importances_mean)[::-1][:10]print([f{wave[j]:.0f}nmforjintop])读法若排列重要性给出的高分波长同样落在 1200–1300 nm 一带说明树模型与 PLS 的 VIP 在哪些波段有用上达成了一致——这是模型可信的正向信号若高分波长散落在 600–750 nm 或 1770–1898 nm第 09 篇已证明这两个区间几乎不含 assay 信息则说明树模型在拟合噪声应立刻回到 PLS。2.4 反直觉的默认值与陷阱项直觉实际GradientBoostingRegressor默认参数默认就是个好模型默认n_estimators100、learning_rate0.1在小样本上很容易过拟合本数据集实测 R²cal 0.9994RandomForestRegressor.n_jobs不影响结果只影响速度但若与random_state配合使用不当如未设种子结果不可复现铁律 9SVR的C/epsilon越大越好epsilon的量纲与 y 相同assay 的 y 波动SD 21.9803 mg远大于 weightSD 5.5842 mg同一组超参数不能跨 CQA 复用树模型的predict输出与 PLS 同形状PLS 输出(n, 1)树模型输出(n,)混用会让r2_score报形状错三、常见报错与排查1.ValueError: Found input variables with inconsistent numbers of samples: [460, 460, 1]现象评价树模型时r2_score(ytest, yp)报错。根因PLSRegression.predict()返回(n, 1)而树模型返回(n,)混用时形状约定不一致。解法统一在预测后加.ravel()。2. 模型性能极好R²cal 0.999现象GradientBoosting 的 R²cal 0.9994、RPDcal 39.85。根因boosting 在小样本高维数据上把噪声也学进去了。解法先算 Q²第 10 篇的q2_and_oof交叉核对若 Q² 远低于 R²cal就是过拟合应降低n_estimators、调小learning_rate或加subsample而不是庆祝高 R²。3. SVR 预测值几乎是一条水平线现象SVR 预测结果方差很小R²test 接近 0。根因epsilon相对 y 的量纲设得过大宽容带吃掉全部残差或C过小导致欠拟合。解法把epsilon按 y 的噪声水平设定assay 的 SD 是 21.9803 mgε1.0 是合理量级若把它设成 20 就会失效并保证输入已做 SNV。4. RandomForest 的 R²test 只有 0.55 左右现象SNV 下 R²test 0.8591换成原始光谱后掉到 0.5455、RPD 1.48。根因未预处理的散射与基线漂移被树当成信号。解法先加 SNV 或 MSC第 03 篇必要时再加 S-G 导数第 04 篇再谈换模型。5.ImportError: No module named lightgbm/xgboost现象按网上教程复制import lightgbm直接失败。根因本系列运行环境无法访问 PyPI这两个包不可安装。解法改用 scikit-learn 自带的GradientBoostingRegressor或HistGradientBoostingRegressor完成同样的梯度提升实验LightGBM/XGBoost在本系列中只能作为概念提及。6. 排列重要性把 600–750 nm 或 1770–1898 nm 排到前列现象2.3 节输出的 top-10 波长落在 600–750 nm 或 1770–1898 nm而第 09 篇的 VIP 与 iPLS 都指向 1120–1248 nm。根因预处理未到位时树会把散射与基线漂移当成信号——iPLS 实测 10 号区间1770–1898 nmR²test 为 −0.1759该波段几乎不含 assay 信息树却给它高分。解法先补 SNV 与 S-G 导数再重算若高分波段仍不在 1200 nm 一带判定该树模型不可信回到 PLS。7.HistGradientBoostingRegressor反而比GradientBoostingRegressor更差现象直方图梯度提升的 R²test 0.8606低于普通梯度提升的 0.8697也远低于 SVR 的 0.9144。根因直方图分箱是为大样本万级以上设计的加速手段155 个校正样本被固定桶数分箱后每桶样本极少分裂点搜索精度反而下降。解法小样本优先用GradientBoostingRegressor或核方法HistGradientBoostingRegressor留给样本充足的场景。8. 改了n_jobs之后结果跟着变现象把n_jobs从 1 改成 −1R²test 在小数第 4 位漂移。根因n_jobs只影响速度、不影响结果结果变了必是别处引入了未固定的随机性铁律 9。解法全流程显式传random_state42。四、动手练习练习 1复现对比表跑通 2.1 的脚本记录五个模型的双集指标。判定标准SVR 的 R²test 应不低于 0.91GradientBoosting 的 R²cal 应不低于 0.99 而 R²test 低于 0.88——这两个特征同时出现才说明复现成功。练习 2预处理对照实验把 RandomForest 分别跑在原始光谱与 SNV 光谱上。判定标准原始光谱的 R²test 落在 0.53 ~ 0.57 之间、RPD 落在 1.45 ~ 1.51 之间SNV 后 R²test 应升到 0.85 以上。练习 3过拟合判据为 GradientBoosting 计算 Q²10 折random_state42并与它的 R²cal 0.9994 对比。判定标准Q² 应明显低于 R²cal若你能观察到 R²cal − R²test 0.12就应判定该模型过拟合并写出至少两条缓解措施。练习 4验证箱型划分打散信息的机理把 RandomForest 的max_depth设为 2、5、10 三档在同一份 SNV 光谱上各跑一次。判定标准随深度增加 R²cal 单调上升而 R²test 不上升深度为 10 时若 R²cal − R²test 0.12即判定树在靠深度记忆噪声须回退浅树或改用 SVR。练习 5树模型与 PLS 的波长对表运行 2.3 节的排列重要性脚本把 top-10 波长与第 09 篇 VIP top-151208–1236 nm对表。判定标准top-10 中落在 1200–1300 nm 的波长不少于 5 个则树模型学到了化学信号若落在 600–750 nm 或 1770–1898 nm 的超过 5 个则其在拟合噪声须先修预处理、再谈换模型。五、小结与下一篇预告树模型不是更先进的 PLS在高维共线、小样本的光谱场景里轴平行切分、共线特征抢分裂点、无法外推这三个弱点会让它落在 PLS 后面而 boosting 更容易用 R²cal 0.9994 这种假象骗人。真正的赢家是 SVR——它在测试集上以 0.9144 / 4.6074 / 3.42 全面超过 PLS 的 0.8979 / 5.0300 / 3.13且校正集 R² 并不夸张0.9426。同时要记住RandomForest 在未预处理原始光谱上会掉到 R²test 0.5455 / RPD 1.48预处理与模型同等重要换模型之前先换预处理。换模型的判据只有一条硬底线R²test 提升至少 0.01 且 RPD 提升至少 0.1才值得承担额外的复杂度与验证成本。SVR 既然赢了就值得专门讲透。**第 13 篇《支持向量回归与高斯过程回归小样本与不确定度》**会拆解 SVR 的 ε 不敏感损失与核技巧并用高斯过程回归的return_stdTrue输出预测带检验它的置信区间是否真的覆盖了真实值。本篇认知问题回显FAQQ1光谱与 CQA 的关系一定是线性的吗A不一定但对 NIR 片剂 assay 这类朗伯-比尔近似成立的场景线性 PLS 已能取得 R²test 0.8979换非线性模型只带来小幅提升。Q2随机森林为什么在高维共线光谱上不一定强于 PLSA树只做轴平行的单特征切分无法表达多波长线性组合且共线特征互相抢分裂机会使重要性被稀释、噪声被放大。Q3梯度提升的 R²cal 达到 0.9994 正常吗A不正常这是典型过拟合。它测试集 R² 仅 0.8697、低于 PLS 的 0.8979说明 boosting 在小样本高维数据上把噪声也学进去了。Q4SVR 为什么能在校正集只有 155 片时超过 PLSArbf 核映射配合 ε 不敏感损失用少量支持向量刻画非线性并忽略部分噪声实测 R²test 0.9144、RPD 3.42均高于 PLS。Q5什么时候才值得放弃 PLS 改用非线性模型A当 R²test 至少提升 0.01 且 RPD 至少提升 0.1 时才值得同时必须先确认预处理已经到位否则提升可能只是预处理的功劳。