
简介面向MATLAB环境下开展多变量统计分析与化学计量学研究的科研人员与研究生压缩包内提供PLS_Toolbox21偏最小二乘回归工具箱可高效处理高维数据与多重共线性问题覆盖标准PLS、正交PLS、PLSR等建模方式并包含中心化、标准化等预处理手段及VIP变量选择、交叉验证、Q2值计算等功能。rar压缩包共209个文件以189个m脚本和20个mat数据文件为主前者涵盖建模、预处理、绘图等模块后者便于直接运行验证与模型演示整包大小仅4.01MB轻量易部署。已有577人浏览学习适合化学计量、生物医学、材料科学等领域的实际应用。通过下载可取得完整工具源码、示例数据集和典型调用流程便于结合自身光谱或组学数据进行建模、预测与结果解释快速掌握PLS及其判别分析、主成分回归等进阶用法。1. 偏最小二乘与 PLS_Toolbox21为什么你调了一周预处理仍不如一套工具箱做近红外定量分析的人应该都有这种体会换一台仪器、换一批样本偏最小二乘模型的预测误差就涨一截。很多人第一反应是改预处理参数一调就是半天RMSECV 纹丝不动。我在反复折腾后换成 PLS_Toolbox21最大的感受不是它比 MATLAB 自带函数快多少而是它把建模流程里最容易翻车的环节——预处理、潜变量选择、交叉验证、离群点诊断——变成了看得见的闭环。这里从偏最小二乘的算法细节讲到工具箱的实际参数设置重点写我如何确定潜变量数、排查离群点。适合做光谱建模、过程分析、近红外和中红外定量模型的人直接照着复现。2. 选型与算法PLS 原理、工具箱核心模块和与内置函数的三个差距在开始建模前我希望先解释为何要选用 PLS_Toolbox21而不是自己写脚本。许多 MATLAB 用户在“工具箱”和“自带函数”之间摇摆。实际上建模项目大多要求快速交付三周后重跑模型还能稳定复现这一点比单次拟合优度更重要。内置的plsregress能做基础的偏最小二乘回归但当你需要交叉验证、预处理管道和离群点诊断时几乎每个环节都要自己写辅助代码。工具箱的价值是把建模中的共性问题变成可复用的对象而不是让你每次从零拼。偏最小二乘本身不是新方法这里只强调一个常被忽略的点PLS 对 X 和 Y 同时做双线性分解提取出的潜变量既要让 X 的方差尽量大又要让 X 得分与 Y 得分的协方差最大。这个区别决定了它比主成分回归更适合光谱这类“特征重叠、样本又少”的数据。2.1 普通最小二乘失效的场景共线性数据与潜变量的由来光谱数据有一个非常恼人的特点特征多样本少。一个中红外谱图动辄上千个波长点而标准样品可能只有三五十个。普通最小二乘要求解回归系数系数个数超过样本数时本身就是欠定系统解不唯一。即使把样本数堆到超过波长数相邻波长的强相关性也会让信息矩阵的条件数高得吓人回归系数对噪声极其敏感训练集拟合得很好第二天换个批次就崩。有人会说那就用主成分回归。PCR 先把 X 压缩成几个主成分再拿主成分当自变量做回归。但主成分只代表 X 自身方差最大的方向并不保证这个方向和 Y 真的有关系。实际项目中经常出现这样的情况前三个主成分累计解释了 X 方差的 98%可 Y 的主要变化却藏在一个方差很小的方向上PCR 会把这个方向直接丢掉。PLS 换个做法每提取一个潜变量同时分解 X 和 Y要求提取的得分向量 t 与 Y 得分向量 u 协方差最大。也就是说PLS 不只找 X 的主要结构而是找“和 Y 最相关”的 X 结构。这个机制带来的直接好处是同样的预测精度PLS 需要的潜变量数量通常比 PCR 少模型解释也更直接。工具箱里潜变量数的参数叫 latent variables缩写为 LV。新手容易把 LV 当成主成分个数来理解其实两者的提取准则差别不小。LV 太多会吸收噪声太少则欠拟合这个度在 4.1 里我会专门展开。2.2 工具箱的核心模块dataset 对象、预处理流、模型与验证的闭环PLS_Toolbox21 的建模链路一般是原始数据 → dataset 对象 → 预处理流 → pls 模型 → 交叉验证模型 → 诊断图。中间每一环都是独立对象可以单独拿出来检查。比如 dataset 对象不只能存 X 和 y还能存变量名、样本名和坐标轴信息。画图时自动显示波长位置不用每次手工 set。我一般会先把 Excel 读进来的数据整理成下面这样的结构% 读入 200 个样本的光谱数据前 120 列是波长最后 1 列是参考浓度 data xlsread(batch_spectra.xlsx); % 201 x 121 的矩阵 X data(:, 1:120); y data(:, 121); % 组织成 dataset 对象并标注维度名称 obj dataset(X, y); obj.name batch1; obj.label{1} {sample}; % 第一维样本 obj.label{2} {wavelength}; % 第二维波长通道dataset 构造函数把 X 和 y 捆成一个对象后续所有模型函数只要传这个对象不用反复传两个数组。更关键的是对象里保存了完整的坐标信息后续画载荷、画得分图的时候横轴就可以直接按波长显示。如果只拿两个普通数组丢进函数换一次仪器就得手工重新排列列非常容易串列。数据集建好以后下一步是定义预处理流。这一步只是登记处理方式真正执行发生在模型训练和预测阶段这样才能保证训练集和测试集用完全相同的预处理参数。比如这里登记一个自动缩放pre preprocess(default, autoscale);这个设计有一个很实际的好处当用测试集预测时模型会自动按训练集保存的均值和标准差做缩放。如果自己分开写标准化测试集漏用训练集统计量是家常便饭这种事我没少见。2.3 与内置 plsregress 的差距为什么我要额外装一套工具箱功能plsregressMATLAB 内置PLS_Toolbox21预处理无需要手工标准化Autoscale、SNV、导数、S-G 平滑支持串联交叉验证简单 k 折Venetian blinds、Contiguous block、随机、按区间离群点诊断无Q 残差与 Hotelling T² 图模型解释回归系数载荷、得分、VIP、变量筛选批次与模型迁移无支持多批次数据、累积建模拿近红外过程分析来说用内置函数往往要自己写一个很长的主流程标准化、建模、交叉验证、画残差图。换一个人做可能漏掉某个环节。工具箱把每步的参数和结果留在模型对象里回头排查时直接调属性。如果数据已经做过严格清洗且样本固定内置函数也能用但整套流程跑下来工具箱至少能省下三分之一的时间而且结果更容易复现。这张表不是纸上谈兵。我做玉米水分模型时拿同一份数据分别跑过内置plsregress和工具箱LV 同样取 8预测结果偏差不大。但换成测试集后内置函数的模型没有给出置信区间和 Q 残差我要排查异常样本只能把 X 重新抽出来算一遍马氏距离这一步非常费劲。工具箱把诊断信息直接挂在模型对象上一个命令就能调出图来这是它真正的价值。3. 建模实测从 Excel 原始光谱到交叉验证模型的可复现步骤第二章讲了为什么选它这一章直接上手。我以一份 200 个样本的近红外光谱为例120 列波长间隔 8 nm覆盖 900 到 1800 nm最后一列是实验室测得的参考浓度。建模前先扫一遍数据看看有没有全零列或空值然后再做模型训练。整个过程的重点不是把 LV 一步设对而是把交叉验证做扎实。3.1 单模型训练pls 命令与模型对象里存了些什么数据经过 2.2 的组织后可以直接调用pls命令。第一个模型不要想着一步到位随便挑一个 LV 数先把流程跑通后面再细化。我一般从 8 个潜变量开始% 用原始 X 和 y 直接建模潜变量数先设定为 8 model pls(X, y, 8, preprocessing, preprocess(autoscale)); % 查看模型概况确认没有 NaN、没有警告 model; % 画出训练集得分图看看有没有明显分群 plotscores(model);pls的前两个参数是光谱矩阵和浓度向量第三个参数 8 是潜变量数最后用名称-值对指定预处理。建模结束后model 对象里存了得分、载荷、回归系数、残差还有一个 detail 结构里面有每次潜变量迭代的方差解释率。训练集拟合效果再好看也不能说明泛化能力所以接下来必须做交叉验证。plotscores 是第一步检查。如果训练集样本明显分成几簇且这些簇对应不同批次或不同时间就要小心潜在的批次效应。这种效应对潜变量数选择影响很大后文会专门讲。3.2 交叉验证分组方式venetian blinds 与 contiguous block 的真实差异PLS_Toolbox21 里交叉验证命令写法比较直接关键在分组方式。默认我常用 ven即 venetian blinds把训练样本按顺序每隔 N 个取一个当验证集。这种分组适合样本顺序完全打乱的场景比如随机抽样的标准样品集。但当样本本身有批次顺序比如同一台仪器连续测量 50 分钟venetian blinds 会让验证样本在时间上离训练样本很近交叉验证误差会非常乐观。这种情况下我更推荐 con即 contiguous block按顺序把每块样本整体切出来验证。这才是模拟换一批新样本时的真实表现。% 从 1 到 10 个潜变量逐一跑交叉验证 rms_cv zeros(1, 10); for nLV 1:10 model_cv pls(X, y, nLV, preprocessing, preprocess(autoscale)); cv crossval(model_cv, cv, ven, 10); rms_cv(nLV) cv.rmsecv; % 交叉验证均方根误差 end % 画 RMSECV 随潜变量数变化的曲线 plot(1:10, rms_cv, -o); xlabel(潜变量数); ylabel(RMSECV);参数说明ven 后面的 10 表示把样本分成 10 份每份轮流当验证集。若改用分块交叉验证需要把 ven 换成 con后面接块大小或块数。在读曲线时注意RMSECV 是越小越好但它的最低点不一定是最佳模型点这个我在 3.3 里细说。再补充一个我在实际项目中常用的小技巧如果样本量少于 50我会直接用随机重复交叉验证比如 rand 加 10多跑几次取平均避免一次分组运气影响判断。后面 5.4 里的蒙特卡洛思路也是这么来的。3.3 RMSECV 曲线解读最低点不是终点折中是常态RMSECV 曲线是最常见的选 LV 依据。很多人看到第 8 个潜变量误差最小就直接取 8。我做过一个土壤有机质模型RMSECV 从第 6 个潜变量开始一路缓慢下降第 10 个潜变量时看起来最优可当我把 Q 残差图调出来从第 9 个开始有 3 个样本的残差突然跳到限外。这说明模型开始把噪声当信号吸收了再往后加 LV 只会过拟合。我把那组数据的误差趋势简化成下面这张表方便说明判断逻辑潜变量数RMSECVQ残差超限样本数60.82070.55080.49190.443100.405如果只看 RMSECV第 10 个潜变量最诱人。但 Q 残差超限样本数从 1 跳到 5说明模型已经在拟合个别样本的异常波动。这种情况下我会退回第 7 或第 8 个潜变量而不是选最低点。我的习惯是先看 RMSECV 是否有明显拐点再看 Q 残差超限比例不超过 5%两者交叉后才能定 LV。4. 参数调优潜变量数、预处理顺序与离群点诊断的实际边界建模最容易改动的是预处理器和潜变量数但最容易被带偏的也是这两个。工具箱把预处理设计成一条流水线顺序不同效果差很多。这一章我写一份实际项目中反复使用的参数检查单按“潜变量个数 → 预处理顺序 → 离群点诊断”的顺序走一遍。4.1 潜变量数定案不看最低点看 Y 方差解释率的边际收益RMSECV 曲线能看趋势但最终定案我还会看一个数字每增加一个潜变量Y 方差解释率提升了多少。工具箱把每次迭代的解释率存在模型对象的 detail 结构里。% 取出 Y 方差累计解释率第二列是累计值 explained model.detail.varY(:, 2); % 逐个 LV 打印解释率 for k 1:length(explained) fprintf(第%d个潜变量后累计Y方差解释率%.2f%%\n, ... k, explained(k) * 100); end我一般这样判断如果第 5 个 LV 后解释率已经到 98.2%第 6 个只增加 0.3%RMSECV 也没有显著下降就停在第 5 个。如果第 6 个虽然解释率增得少但 RMSECV 从 0.6 降到 0.45这个 LV 仍然值得保留。解释率和 RMSECV 要配合看单看一个都会出错。还有一个细节不同预处理组合下解释率曲线形态完全不同。同一份数据Autoscale 和 SNV 做出来的曲线可能一个陡峭一个平缓。所以 LV 数不能跨预处理方案比较每次改预处理后都要重新选 LV。4.2 预处理顺序SNV、导数、Autoscale 的排列次序与窗口参数光谱预处理不是越多越好而是顺序要对。我处理的固体颗粒光谱比如饲料、土壤颗粒散射会造成基线漂移和样本间光程差异第一步我会用 SNV 做标准正态变量变换。然后做一阶导数把基线漂移变成一次项减小低频背景干扰。最后再做 Autoscale把每个波长缩放到单位方差避免高吸收峰在计算中主导整个模型。% 定义预处理流水线SNV - 一阶导数(S-G平滑) - Autoscale pre preprocess(default, snv); pre preprocess(derivative, 1, pre, window, 15, order, 2); pre preprocess(autoscale, pre); % 用这条预处理流重新建模 model pls(X, y, 8, preprocessing, pre);这里window15是 Savitzky-Golay 平滑窗口order2是多项式阶数。窗口太大导数会把相邻波峰拉平窗口太小噪声被放大。土壤近红外我常用 9 到 21汽油近红外可以用 11具体要根据光谱分辨率自行调整没有绝对标准。顺序上有一条铁律散射校正和导数要放在 Autoscale 之前。如果先 AutoscaleSNV 的除法本来就带有缩放效果两者叠加会让一部分波长被过度放大后面再做导数时会把噪声一块放大。我自己第一次就掉进过这个坑后来养成了从后往前检查预处理列表的习惯。4.3 离群点诊断Q 残差与 Hotelling T² 的阈值到底看哪个PLS_Toolbox21 的诊断图里有两张常用图Hotelling T² 和 Q 残差。T² 衡量样本在主成分空间的距离Q 残差衡量样本在残差空间里的距离。理解这两者非常关键。T² 超限的样本说明它在模型覆盖范围之外可能是一个浓度很高的正常样本也可能是一个完全不同的物质带来的干扰。Q 残差超限的样本说明它没有被模型结构解释的部分太多常见原因包括仪器噪声突然变大、样品出现气泡或者光谱采集位置偏移。我建议把诊断图打开之后再删样本。判断准则可以写成一个简单的优先级% 画出 Q 残差对应 Hotelling T² 的诊断图 plotdiagnostics(model, Qhotellings);如果样本只在 T² 上超限且它的原始光谱形态正常多半是高杠杆点保留即可。如果样本 Q 残差超限且横向对应一个明显的毛刺峰先检查原始光谱再决定是否删除。删除前先想一下这个样本在量产时会不会出现如果不会删掉它再重新建模是合理的如果会那就应该把它当正常波动留下来。离群点删除后必须重新训练模型并重新做交叉验证不能把删除前后的性能指标混在一起对比。这个小细节经常被人忽略导致选出的 LV 数根本不可比。5. 常见问题与避坑记录PLS_Toolbox21 踩坑后的四个排查手段工具箱用久了会发现很多问题不是算法问题而是调用方式问题。我在几个项目里反复踩过同类坑按“现象 → 原因 → 解决”的方式记录在这里。把这些排查手段背下来能省不少调试时间。5.1 建模报 NaN最隐蔽的预处理除零现象预处理跑完之后调用 pls返回的模型对象里有 NaN画图时得分图一片空白。原因最常见的有两种情况。数据里有空值或者某一列的方差为 0。Autoscale 会把每个变量除以标准差零方差变量直接产生除零错误。光谱两端经常出现探测器坏点那里全是同一个数值正好命中。解决建模前先做一步删除操作把含 NaN 和零方差的列筛掉。% 删除含 NaN 的列同时保留方差大于 0 的列 keep ~any(isnan(X), 1) (var(X, 0, 1) eps); X X(:, keep);这一步做完再看 y确认 y 里也没有 NaN。大多数 NaN 都来自 Excel 表格中的空单元格导入后显示为 NaN工具箱自身不会帮你自动清理。5.2 交叉验证误差与外部测试集差距巨大现象交叉验证报 RMSECV 为 0.35但把模型放到新批次数据上预测RMSEP 达到 1.2差了三倍以上。原因分组方式不符合真实应用。样本按批次连续采集时我用了 ven同一批次的样本被分散到训练集和验证集里等于让模型“见过”同一批邻近样本。换到真正的异地批次性能自然掉下来。另一个常见原因是预处理统计量用了全集计算而不是只用训练集。解决改用 con 分组或手动按批次留出测试集。最稳妥的做法是直接抽出一个批次当测试集其余样本做训练train obj(1:160, :); test obj(161:200, :); model pls(train, preprocessing, pre); pred predict(test, model);验证集必须从建模流程的一开始就隔离出去。我见过有人把所有样本一起做 Autoscale 后再拆分这等于提前把测试集信息泄露给了模型外部测试结果自然虚高。5.3 工具箱版本接口不一致导致旧代码失效现象老代码用 pls(x, y, ncomp) 直接跑换到 PLS_Toolbox21 之后报 Undefined function 或者参数数量错误。原因新版本对输入参数做了更严格检查部分老函数改成了对象方法。工具箱更新后旧脚本没有同步修改接口自然对不上。解决先用ver查看版本再用help pls确认当前签名。通常需要把第三个参数从整数改成显式的潜变量数选项或者把两个输入数组改成 dataset 对象。改完之后重新跑一次交叉验证确认结果和旧版本基本一致再继续用。需要提醒的是版本更新可能让之前保存的模型对象无法直接被新版本读取备份模型时同时备份数据、预处理参数和 LV 数别只存一个 .mat。5.4 潜变量数只看最低点导致过拟合现象RMSECV 最低点在 10 个潜变量处模型训练集预测非常好但半年后用到新样本上误差比预期大很多。原因RMSECV 最低点往往是模型开始吸收噪声的位置。小样本场景下单个最低点很容易由偶然误差产生不具备稳定性。解决用随机重复交叉验证代替单次分组看多次分组的平均误差和波动范围。% 重复 50 次随机 10 倍交叉验证 rms_cv zeros(50, 1); for i 1:50 cv crossval(model_cv, cv, rand, 10); rms_cv(i) cv.rmsecv; end fprintf(平均 RMSECV%.3f标准差%.3f\n, ... mean(rms_cv), std(rms_cv));如果最优点附近的标准差特别大说明这些 LV 数对分组方式太敏感应该考虑降低一两个 LV。这个方法比单次交叉验证多花几十秒但换来的是长期稳定性。6. 进阶用法模型解释、系数导出与换仪器后的复验习惯模型建完不等于项目做完。交付时我固定做三件事第一用潜变量载荷找出对浓度变化贡献最高的波长区间第二把回归系数和 VIP 值连同样本号一起导出第三对测试集做一次严格的留存预测给工艺人员看结果。这三件事做完模型才不算黑匣子。提取重要变量时我常看 VIP 值。VIP 大于 1 的变量一般被认为对模型比较重要。在工具箱里可以用类似下面的代码把 VIP 值和回归系数一起导出% 提取变量重要性投影得分不同版本函数名可能略有差异 vip vipscore(model); coeffs model.regcoef; % 按 VIP 降序排列并打印前十个最重要变量 [~, idx] sort(vip, descend); fprintf(前10重要变量的原始索引%s\n, num2str(idx(1:10))); % 保存成带变量名的表格 T table((1:length(coeffs)), coeffs, vip, ... VariableNames, {wavelength_index, regcoef, vip}); writetable(T, pls_model_summary.csv);拿到索引后我会把它映射回波长坐标写进报告里“有效变量集中在 1100 到 1250 nm 区间”。这一步常被新手跳过但工艺人员最关心的是这个区间因为它对应 O-H 和 C-H 的合频振动。这样一写模型马上就变成可沟通的技术结论。换仪器或者换测量条件时我会强制走一遍复验流程先把新仪器采到的 10 个样本放进模型预测看 Q 残差有没有大范围超限。如果超限先检查数据预处理是否一致再检查波长范围和分辨率。不要一上来就重新建模。很多时候只是波长范围截取不一致模型被生生带偏。从那以后我每次交付模型都会强制导出这张汇总表无论对方是否要求。这个习惯在半年后重新复现模型时特别有用因为你能直接看到当时选了哪几个潜变量、哪几个波长最重要不用再把老代码翻一遍。希望这些踩坑记录能帮你少走一点弯路也希望这个导出习惯对你接下来的模型交付有实际帮助。本文还有配套的精品资源点击获取