制剂 CQA 预测模型开发教程（9）：变量选择——VIP、iPLS 与 CARS 的波长筛选实战

制剂 CQA 预测模型开发教程9变量选择——VIP、iPLS 与 CARS 的波长筛选实战版本声明块工具/软件Python 3.10.11 scikit-learn 1.7.2 numpy 2.2.6 scipy 1.15.3数据NIR Shootout 2002 片剂数据集655 片 × 650 波长本文目标用十行代码从任意PLSRegression模型算出 VIP 分数并用 iPLS 找出信息最密集的波段同时看清变量选择不等于精度提升。一句话结论在 NIR Shootout 2002 的 assay 数据SNV 预处理上VIP1 的波长共 193/650、集中在 1208–1236 nm 的 C-H 键二级倍频区筛选后 R²test 仅由 0.9074 微升至 0.90980.0024变量数却降到 29.7%真正靠少变量取胜的是 iPLS 第 5 区间1120–1248 nm仅 65 个波长即取得 R²test 0.9263 / RMSEP 4.2745 / RPD 3.69。〇、本篇要解决的认知问题650 个波长点全都同等重要吗VIP 分数凭什么能量化哪个波长重要VIP 1 这个阈值是拍脑袋定的还是有数学来源为什么sklearn.cross_decomposition.PLSRegression没有vip_属性必须自己算iPLS 与 CARS 都在做波长选择它们的搜索粒度差在哪删掉无信息变量一定能让模型变准吗一、机制解析1.1 从 PLS 的三个属性推导 VIP偏最小二乘回归PLS在建模时已经把方差结构拆成了三块只是 scikit-learn 没有把它们组合成一个重要性指标X (n×p) ──NIPALS──► x_scores_ t : (n × A) 样本在潜空间的坐标 x_weights_ w : (p × A) 每个波长对潜变量的贡献方向 y_loadings_ q : (m × A) 每个潜变量对 Y 的解释强度变量投影重要性VIPVariable Importance in Projection要做的事就是把方向贡献与解释强度相乘再按潜变量加权求和。第 a 个潜变量对 Y 的解释平方和记为SSY_a q_a² · (t_aᵀ t_a) # q_a 是 y_loadings_ 第 a 列t_a 是 x_scores_ 第 a 列于是第 j 个波长的 VIP 分数为VIP_j sqrt( p · Σ_a [ SSY_a · (w_ja / ‖w_a‖)² ] / Σ_a SSY_a )其中p是变量总数650‖w_a‖是第 a 个权向量的欧氏范数。物理含义很清楚SSY_a / ΣSSY_a 是该潜变量的权重w_ja/‖w_a‖ 是归一化后该波长在该潜变量上的贡献。一个波长若在多个解释力强的潜变量上都贡献大VIP 就高。这里的SSY_a不是随手定义的NIPALS 把 Y 载荷取为q_a (yᵀt_a)/(t_aᵀt_a)于是把 y 对单个潜变量t_a做一元回归时它的回归平方和恰好就是q_a²·(t_aᵀt_a) SSY_a。也就是说SSY_a 就是第 a 个潜变量单独能解释多少 y 的方差。又因为 NIPALS 逐轮正交化各潜变量得分满足t_aᵀt_b 0 (a ≠ b)这些解释量可以直接相加而不产生交叉项——这正是 VIP 定义式能写成求和式的几何前提。1.2 VIP 1 阈值的数学来源阈值 1 不是经验拍脑袋而是一个归一化恒等式的结果。对任一潜变量 a权向量按列归一化后满足Σ_j (w_ja / ‖w_a‖)² 1把它代入 VIP 定义并求全体波长的平方和Σ_j VIP_j² p · Σ_a [ SSY_a · Σ_j(w_ja/‖w_a‖)² ] / Σ_a SSY_a p · Σ_a SSY_a / Σ_a SSY_a p也就是说VIP² 的全体平均恰好等于 1。所以VIP 1等价于该波长的重要性高于平均水平——阈值来自mean(VIP²) 1这个恒等式而不是人为规定。写代码时可以用它自检算完 VIP 后np.mean(vip**2)应当等于 1。1.3 VIP 的两个隐含前提VIP 的公式很短但它有两个容易被忽略的前提直接决定了 VIP 什么时候能信、什么时候不能信。前提一VIP 只有在潜变量数固定之后才可比较。权重SSY_a / ΣSSY_a依赖建模用了几个潜变量nLV 一变例如从 5 到 19新增潜变量会分走一部分 SSY 权重VIP 排序与数值随之改变。在两个不同 nLV 的模型之间直接比 VIP 没有意义——必须先按第 08 篇把 nLV 定下来铁律 5。本篇 VIP 数字均绑定在 nLV5 的模型上。前提二VIP 只回答重要不回答方向。公式中的(w_ja/‖w_a‖)²把符号平方掉了所以一个与 y 强烈负相关的波长其 VIP 同样很高。VIP 高只说明该波长参与了区分样本不能读成该波长升高时 assay 也升高。判断方向要回到 PLS 回归系数coef_或载荷的符号第 14 篇展开。1.4 iPLS把 650 个点切成整段来搜索iPLSinterval PLS区间偏最小二乘的思路与 VIP 完全不同它不做逐变量打分而是把光谱等分成若干连续区间每个区间单独建模用性能指标通常是交叉验证 RMSECV工程上也可用测试集 RMSEP 做教学演示挑出最好的区间。波长轴 600 ────────────────────────────────────────────► 1898 nm | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 600-728 ... 1120-1248 ... 1770-1898 65 点/段共 10 段 × 65 650 点 每个区间 → 独立 PLS → 比较 R²test / RMSEP / RPD → 选最优段iPLS 的隐含假设是携带同一化学信息的吸收带是连续的。这与 NIR 的物理事实吻合——C-H、N-H、O-H 的倍频与合频吸收本来就落在连续波段里。区间宽度与潜变量数之间还藏着一层耦合这是 iPLS 最容易被忽略的地方区间内变量数 p_int 650 / k 该区间 nLV 上限 min(n_samples − 1, p_int)k 越大区间越窄p_int 越小该区间能支撑的潜变量越少一旦 p_int 小于所需的 nLV模型连一个完整潜变量都装不下。反过来k 太小区间太宽又会把有用信息带和噪声带混在同一段里。本篇实测 k10每段 65 点时区间 9 在 65 个变量上选了 nLV20这本身就是一个过拟合信号它对应 R²test 只有 0.7871。换句话说表里哪一段最好同时受两个因素牵引这段波段含有多少信息以及这段波段能撑起多复杂的模型。要让区间之间公平可比应固定 nLV 后再比较或对每段各自用交叉验证选 nLV本篇采用后者并在报告里写明用的是哪一种。1.5 六种变量选择策略的对照把常见的几类波长选择策略放在一张表里对照选型时按先看数据规模、再看化学先验的顺序判断策略粒度是否用到 Y随机性适用场景主要风险全波段不选无否无基线模型、变量数与样本数同量级共线性与噪声变量稀释信号相关系数阈值单变量是逐点无快速探查、给出波段大致轮廓忽略共线性删掉单独弱、组合强的协同波长VIP 阈值本篇单变量是经潜变量无需要可解释的波长清单阈值处排序敏感换 nLV 即变不区分正负iPLS 单区间连续区间是区间级无信息带连续、在线端要求少变量区间边界人为划定可能切断信息带iPLS 前向/后向组合多区间是区间级无单区间不足、需跨带组合信息组合搜索次数随区间数增长易过拟合验证集CARS单变量是多轮系数有变量极多、需要自动降维采样随机性大必须固定种子铁律 9读表要点只有一句用不用到 Y与有没有随机性是两把标尺。VIP 与 iPLS 都不含随机步骤天然可复现适合受监管场景CARS 引入蒙特卡洛采样可复现性完全依赖种子管理而相关系数阈值虽然最简单却是唯一会系统性漏掉协同变量的做法。1.6 CARS指数衰减的逐变量淘汰赛CARSCompetitive Adaptive Reweighted Sampling竞争性自适应重加权采样是更细粒度的搜索机制分三步蒙特卡洛采样每次随机抽取一部分样本建 PLS 模型用回归系数的绝对值 |b_j| 作为该轮每个变量的重要性得分。指数衰减函数EDF强制降维第 i 轮保留比例r_i a·exp(−k·i)其中 a 与 k 由首轮与末轮的保留比例解出例如首轮保留全部、末轮保留 0.001 时k ln(1/0.001)/(N−1)。变量按 |b_j| 排序只保留|b_j| r_i的前若干个因此变量数逐轮下降。自适应重加权采样 交叉验证定优在保留变量中按权重再采样建模型跑 N 轮后取 RMSECV 最小的那一轮变量子集。iPLS 与 CARS 的差异可以用一张表说清维度iPLSCARS搜索粒度区间整段保留或丢弃单个波长搜索方式穷举所有候选区间蒙特卡洛采样 指数衰减淘汰重要性来源区间级模型性能变量级 |PLS 回归系数|迭代轮次单轮可扩展为前向/反向多轮多轮变量数单调下降主要风险区间边界人为划定可能切断信息带采样随机性大必须固定随机种子铁律 9二、完整代码与逐行剖析2.1 公共部分数据加载与 SNVimportnumpyasnpfromscipy.ioimportloadmatfromsklearn.cross_decompositionimportPLSRegressionfromsklearn.metricsimportr2_score,root_mean_squared_errordefload_shootout(path):加载 NIR Shootout 2002 片剂数据集返回三集光谱与三个 CQA。mloadmat(path)defunpack(name):# MATLAB DataSet Object 外层是 (1,1) 结构体必须 [data][0, 0] 逐层解包arrm[name][data][0,0]# 原始 dtype 为大端 f8显式转 float64否则后续矩阵运算类型不匹配returnnp.asarray(arr,dtypenp.float64)return{Xcal:unpack(calibrate_1),Xval:unpack(validate_1),Xtest:unpack(test_1),Ycal:unpack(calibrate_Y),Yval:unpack(validate_Y),Ytest:unpack(test_Y),wave:unpack(axisscale).ravel(),# 600 → 1898 nm步长 2 nmcqa:[weight,hardness,assay],# 列顺序已实测}defsnv(X):逐样本标准正态变量变换每个样本减自身均值、除自身标准差。muX.mean(axis1,keepdimsTrue)# keepdims 保证 (n,1) 可广播不能用全局均值sdX.std(axis1,ddof1,keepdimsTrue)return(X-mu)/sd2.2 从零实现 VIP十行核心代码defvip_scores(model,p):从已拟合的 PLSRegression 模型计算 VIP 分数。 model.x_scores_ : (n, A) 潜变量得分 t model.x_weights_ : (p, A) 权重 w model.y_loadings_: (m, A) Y 载荷 q p : 变量总数 光谱波长数 tmodel.x_scores_# (n, A)wmodel.x_weights_# (p, A)qmodel.y_loadings_.ravel()# (A,) —— 注意不是 (A,1)ssy(q**2)*(t**2).sum(axis0)# SSY_a q_a² · Σt_a²wn(w/np.linalg.norm(w,axis0))**2# 按【列】归一化后再平方returnnp.sqrt(p*(wn*ssy).sum(axis1)/ssy.sum())# (p,)if__name____main__:dload_shootout(nir_shootout_2002.mat)Xsnv(d[Xcal])# 铁律 3SNV 逐样本无需在训练集 fityd[Ycal][:,2]# 第 3 列是 assaymgwaved[wave]plsPLSRegression(n_components5,scaleTrue)# 与 SPEC 基线一致的 nLV5pls.fit(X,y)vipvip_scores(pls,X.shape[1])print(mean(VIP²) ,round(float(np.mean(vip**2)),6))# 应恒等于 1.0用于自检print(VIP 1 的波长数,int((vip1).sum()),/,X.shape[1])ordernp.argsort(vip)[::-1][:15]# 取 VIP 最高的 15 个波长forrank,jinenumerate(order,1):print(f{rank:2d}{wave[j]:6.0f}nm VIP{vip[j]:.3f})这段代码在真实数据上的输出已实测为mean(VIP²) 1.0VIP 1 的波长数 193 / 650VIP 取值范围0.329 ~ 2.162。VIP 最高的 15 个波长与分数如下本系列实测值未经修饰排名波长 (nm)VIP排名波长 (nm)VIP排名波长 (nm)VIP112182.162612242.0691112101.850212202.153712262.0111212321.805312162.145812122.0091312341.732412222.120912281.9451412361.661512142.0981012301.8761512081.614前 15 名全部落在1208–1236 nm这是一个可以直接用化学解释的现象该区间是C-H 键二级倍频second overtone区。片剂 API 含量变化直接改变 C-H 基团的浓度NIR 光谱在这个波段对此最敏感——统计重要性最高的波长恰好落在化学上说得通的吸收带上这才是变量选择可信的前提。2.3 变量选择的效果必须如实报告的结论把 VIP1 的 193 个波长留下来重新建模与全波段 650 个波长对比assaySNV变量集变量数R²testRMSEPRPDtest全波段6500.90744.79253.29VIP 1 筛选19329.7%0.90984.72973.33真实结论VIP 筛选后性能仅微弱提升R² 0.0024它的真正价值是把变量从 650 个削减到 193 个29.7%。这一点必须讲清楚——很多教程会把变量选择包装成精度利器但在本数据集上它不是。变量削减的真实收益在于模型存储更小、在线预测更快、维护时可解释的波长更少而精度上几乎没有代价。把 0.0024 说成显著提升就是失真。2.4 iPLS 十区间扫描defipls_scan(X,y,Xtest,ytest,wave,n_intervals10):把波长轴等分成 n_intervals 段每段独立建模返回逐区间性能表。results[]foridx,subinenumerate(np.array_split(np.arange(X.shape[1]),n_intervals),1):Xi,XtX[:,sub],Xtest[:,sub]# 只取本区间的列# 每段的最优潜变量数独立搜索这里用测试集 R² 演示区间差异正式项目须用 CV铁律 5bestNonefornlvinrange(1,21):mPLSRegression(n_componentsnlv,scaleTrue).fit(Xi,y)r2r2_score(ytest,m.predict(Xt).ravel())ifbestisNoneorr2best[1]:best(nlv,r2,m)nlv,r2,mbest rmseroot_mean_squared_error(ytest,m.predict(Xt).ravel())rpdnp.std(y,ddof1)/rmse results.append({interval:idx,n_wave:len(sub),nLV:nlv,wave_lo:wave[sub[0]],wave_hi:wave[sub[-1]],R2test:r2,RMSEP:rmse,RPDtest:rpd,})returnresults在 650 点等分 10 段每段 65 点的实测结果区间波长范围 (nm)nLVR²testRMSEPRPDtest1600–72830.437611.80831.332730–85810.012915.64341.013860–98860.87595.54692.844990–111830.91004.72343.3451120–124880.92634.27453.6961250–137860.91234.66253.3871380–150840.89795.03243.1381510–163880.91714.53383.4891640–1768200.78717.26552.17101770–18981−0.175917.07410.92这张表有三条值得记住的信息最佳区间 51120–1248 nm仅用 65 个波长就超过全波段 650 个波长R²test 0.9263 vs 0.9074RPD 3.69 vs 3.29。它与 VIP 最高峰的 1208–1236 nm 高度重合——两种完全不同的方法指向了同一片化学信息带这种交叉验证才是变量选择值得相信的理由。最差区间 101770–1898 nmR²test 为 −0.1759比直接用均值预测还差说明该波段几乎不含 assay 信息这一带更多对应 O-H、N-H 合频与其他组分的强吸收。区间 1600–728 nm与区间 2730–858 nm性能同样很差RPD 1.33、1.01属于典型的无信息噪声区。2.5 容易踩的默认值API / 参数常见误解实际行为PLSRegression.scale光谱量级相近所以缩放无害默认True会逐列自动缩放反而破坏原始方差结构第 11 篇详述PLSRegression.y_loadings_以为是 (A,) 一维数组实际是 (1, A)直接乘 (p, A) 会广播错误x_weights_归一化方向以为按行归一化必须按【列】归一化axis0np.array_split用len() // k手动切650 不能整除时手切会丢尾巴array_split自动分配余数2.6 工程视角变量选择的三条经验法则非铁律以下三条为经验法则不占用十条铁律的额度收益优先看变量数而不是 R²。本数据集 VIP 筛选把 R²test 从 0.9074 抬到 0.90980.0024几乎落在噪声内真正的收益是变量数降到 29.7%。若目标是在线端PAT 场景的预测速度与维护成本这就够了若目标是继续涨 R²应换预处理或换模型族第 12、13 篇而不是抠波长。比较必须在交叉验证保护下做。用同一套折固定random_state42铁律 9对全波段与筛选后跑 CV 比 RMSECV 或 Q²测试集只出现一次铁律 2。筛选结果要过化学一致性这一关。若选中的波长落在噪声区或与目标 CQA 无化学关联的波段即便 R² 略升也应怀疑过拟合本篇 VIP 前 15 名全部落在 1208–1236 nm 的 C-H 二级倍频区才是可信的证据。三、常见报错与排查1.ValueError: operands could not be broadcast together with shapes (650,5) (1,5)现象算 VIP 时形状对不上。根因model.y_loadings_的形状是(n_targets, n_components) (1, A)被当成(A,)参与了运算。解法先q model.y_loadings_.ravel()压成一维再算q**2。2.mean(VIP²)不等于 1现象算出来的 VIP 全部大于 1或全部小于 1阈值判断失效。根因x_weights_归一化时写成了w / np.linalg.norm(w)全局范数而不是w / np.linalg.norm(w, axis0)逐列范数。解法加上axis0并用np.mean(vip**2) 1做断言自检。3. VIP 筛选后测试集指标暴涨现象筛选后 R²test 从 0.9 跳到 0.99。根因把校正集与测试集拼在一起做变量选择或 PLS 拟合测试集信息泄漏进了变量选择过程违反铁律 2测试集只能用一次。解法变量选择的所有拟合都只在校正集或校正集内部的交叉验证上完成测试集只在最后评估时出现。4. iPLS 切分报错或最后一段样本数异常现象ValueError: Found array with 0 feature(s)或第 10 段只有 5 个波长。根因用650 // 10 65手写切片若波长数不能被区间数整除如 651尾部会被丢弃或切出空数组。解法统一用np.array_split(np.arange(p), n_intervals)它会自动把余数分配到前几段。5. iPLS 结果每次运行都不一样现象同一个数据集两次运行的最佳区间不同。根因区间内部用了带随机性的交叉验证却没有固定种子。解法所有KFold/ShuffleSplit显式传random_state42铁律 9。四、动手练习练习 1阈值自检在自己的 PLS 模型上运行vip_scores()断言abs(np.mean(vip**2) - 1.0) 1e-9。若断言失败检查x_weights_的归一化轴。判定标准np.mean(vip**2)与 1 的偏差小于 1e-9。练习 2VIP 阈值扫描分别取 VIP 阈值 0.8、1.0、1.2、1.5记录每个阈值留下的变量数与测试集 R²test。判定标准能观察到变量数单调下降且 R²test 的变化幅度小于 0.02即在本数据集上变量选择不带来显著精度变化。练习 3改变区间数把 iPLS 的n_intervals从 10 依次改为 5、13、26 重跑。判定标准当区间数达到 26每段 25 点时最佳区间的 R²test 应仍不低于 0.90若某次出现 R²test 0记录是哪个波长区间并对照 1770–1898 nm 的解释。五、小结与下一篇预告VIP 与 iPLS 是两种粒度不同的波长选择方法前者从 PLS 的权重与载荷推导出逐变量分数阈值 1 来自mean(VIP²) 1的归一化恒等式后者把光谱切成连续区间整体筛选。在 NIR Shootout 2002 的 assay 上两者都指向 1100–1250 nm 的 C-H 二级倍频带但只有 iPLS 真正以小博大——65 个波长取得 R²test 0.9263 / RPD 3.69而 VIP 筛选 193 个波长只把 R²test 从 0.9074 微升到 0.9098。变量选择的核心价值在精简与可解释而非精度。CARS 的指数衰减权重机制比 iPLS 更细但随机性更大必须固定随机种子。没有评价指标就无法比较任何两个模型——既然本篇反复在用 R²test 与 RPD下一篇就一次性把它们讲透。**第 10 篇《模型评价体系RMSEC/RMSEP/R²/Q²/RPD/RER》**会给出六个指标的公式与关系、RPD 分级经验法则并复现 scikit-learn 1.7.2 下mean_squared_error(squaredFalse)的TypeError与修复方案。本篇认知问题回显FAQQ1PLS 的 VIP 分数怎么从 scikit-learn 的 PLSRegression 属性里算出来A用x_scores_、x_weights_、y_loadings_三个属性按VIP_j sqrt(p·Σ_a SSY_a·(w_ja/‖w_a‖)²/Σ_a SSY_a)计算其中SSY_a q_a²·Σt_a²。Q2VIP 大于 1 的阈值是怎么来的是经验法则吗A不是拍脑袋。权向量逐列归一化后Σ_j VIP_j² p故 VIP² 的全体均值恒为 1VIP1 即重要性高于平均。Q3为什么 sklearn 的 PLSRegression 没有 vip_ 属性Ascikit-learn 只暴露 NIPALS 分解结果VIP 是化学计量学的后处理指标需用户用x_weights_等属性自行组合官方不为它提供属性。Q4iPLS 与 CARS 做波长选择有什么本质区别AiPLS 以连续区间为单位穷举比较、整段保留或丢弃CARS 以单个波长为单位靠指数衰减函数逐轮淘汰低|PLS 回归系数|变量。Q5删掉无信息波长一定能让光谱模型变得更好吗A不一定。本系列在 assay 数据上实测 VIP1 筛选后 R²test 仅由 0.9074 升至 0.90980.0024收益主要是变量数降到 29.7% 的精简。