正交实验设计与数据处理：L9(3^4)、极差分析和方差分析实战

我带的第一个配方项目里有位做了十几年实验的老师傅看到我摆出来的那张 L9(3^4) 表格第一句话是这不就是每次少做几次实验嘛。这话对了一半。正交表确实把 4 因素 3 水平本来要做的 81 组压缩成了 9 组但它真正值钱的地方不在少做而在于这 9 组做完之后你能对那没做的 72 组给出有统计依据的推断——哪几个因素真的在起作用哪个水平组合最优误差有多大预测值能落在什么区间。少了后面这半截少做实验就只是偷懒。正交实验设计和数据处理这两件事在很多人手里是割裂的设计的时候照着表格抄一组试验号数据处理的时候在 Excel 里拉几个平均数挑个最高的就宣布找到了最优配方。中间那些 k 值、R 值、平方和、自由度、F 检验要么被跳过要么被当成走流程的形式。我写这篇是想把这条链路完整走一遍——从因素水平的确定到正交表的选型到极差分析和方差分析再到 Python、Excel、MATLAB 三种工具下的具体实现最后把我这些年踩过的坑挨个摊开讲。适合谁看如果你正在准备做一组多因素实验或者手里已经有一批正交实验数据但不太确定怎么解读或者你被要求用代码把整套流程自动化跑起来这篇应该都能用上。1. 先算清正交表替你省了什么L9(3^4) 的均衡性与自由度约束1.1 从 81 组到 9 组省掉的到底是哪些组合假设你要考察 4 个因素每个因素取 3 个水平。全因子实验的组合数是 3×3×3×381 组。做 81 组是什么概念如果每组实验从准备到出结果要 4 小时那就是 324 小时接近连续干两周不休息。多因素实验的规模是乘法级增长的加到 5 个因素就是 243 组到 6 个因素就是 729 组任何实验室都扛不住。正交表 L9(3^4) 只要求你做 9 组。理解它省了什么关键看省的方式。它不是随机砍掉 72 组而是让 9 组样本满足一条特殊性质对任意两个因素它们的 9 种水平组合中每一种都恰好出现一次。这就是所谓的均衡分散。你可以拿标准的 L9(3^4) 表验证一下第 1 列和第 2 列交叉看(1,1)、(1,2)、(1,3)、(2,1)……一直到 (3,3)9 种组合一个不多一个不少各出现 1 次。这条性质带来的直接结果是虽然每个因素只做了 9 次实验但每个水平1、2、3都出现了 3 次而且这 3 次是在其他因素所有水平的均匀搭配下完成的。也就是说当你比较因素 A 的第 1 水平和第 2 水平时B、C、D 的干扰被自动抵消掉了——这就叫整齐可比。这一条才是极差分析能成立的数学基础。1.2 正交性的两条硬约束均衡分散与整齐可比我见过不少人在这一步卡住为什么正交表上的行可以随便排、列可以随便分配但结论却对得上因为正交表的构造本身满足两条约束只要你不破坏每列水平出现次数相同这个前提统计推断就成立。第一条是每列中每个水平出现次数相同。L9(3^4) 每列各水平出现 3 次L8(2^7) 每列各水平出现 4 次。第二条是任意两列的水平组合出现次数相同。这两条合起来保证了正交表上的数据可以用分离各因素效应的方式去分析。注意一旦你在试验过程中临时改了某个因素的某个水平比如某一组温度设错了正交性就被破坏了后面所有极差分析和方差分析的结果都会失真。这种情况的补救办法是重新补做那一组而不是反正误差不大改了算了。我吃过这个亏。有次做 L16(4^5)第 5 列是一个空列操作员觉得空列没意义就随便填了水平结果最后算误差平方和的时候发现那一列的离差比所有因素都大——因为空列本应反映纯误差被人为改动之后它变成本身就带效应的假误差F 检验全乱了。1.3 列的自由度怎么算以及 L9(3^4) 为什么最多塞 4 个三水平因素正交表括号里的数字含义要拆开看L9 是试验次数3 是水平数4 是列数。选表时有两条硬性约束必须同时满足。约束一表的列数 ≥ 你要安排的因素个数加上交互作用占用的列。约束二表的总自由度 ≥ 试验要考察的总自由度。自由度的算法是一个因素的自由度 水平数 - 1。三水平因素自由度是 2二水平因素自由度是 1。L9(3^4) 的总自由度 9 - 1 8而 4 个三水平因素的自由度之和 4×2 8刚好占满。这就是为什么 L9(3^4) 上排 4 个三水平因素时一点误差自由度都不剩只能靠重复实验来估误差。常用正交表的基本参数我整理成了下面这张表选表的时候直接对号入座| 正交表 | 试验次数 | 列数与水平 | 可安排的三水平因素上限 | 剩余误差自由度 | | L9(3^4) | 9 | 4 列全三水平 | 4 | 0不重复时 | | L18(2^1×3^7) | 18 | 8 列混合水平 | 7 | 2 | | L27(3^13) | 27 | 13 列全三水平 | 13 | 0 | | L8(2^7) | 8 | 7 列全二水平 | — | 0 | | L16(4^5) | 16 | 5 列全四水平 | — | 0 |看这张表能得出一个很实用的结论试验次数正好等于自由度数加 1 的那些表都是饱和的。饱和表只适合主效应分析一旦要考虑交互作用或者要估误差就得往上跳一档改用 L18 或 L27。这个取舍我在第 2.5 节还会展开。2. 因素与水平的确定这一步草率后面所有分析都白做2.1 因素筛选别把能想到的都放进去正交表的槽位是有限的但更大的限制是实验成本。我经常看到有人列了七八个因素一股脑全塞进表里理由是反正正交表能排下。问题在于因素越多每个因素的效应水平被其他因素稀释的可能性越大而且你最终拿到一份主次顺序表很可能前五个因素的 R 值都很接近谁主谁次根本分不出来。我的做法是先做一轮因素筛选。跟工艺、配方、设备相关的一线人员聊让他们列出所有可能影响指标的变量然后按两个维度打分一是理论上有没有明确的因果链二是现有条件下能不能方便地调节并在试验中稳定控制。因果链不清楚、或者调起来很费劲的先剔除。通常一轮筛下来能砍掉一半。剩下 3 到 5 个因素才是正交实验的合理规模。如果你实在放不下某个因素可以把它先固定在最有把握的水平上等主实验做完再单独用单因素实验去微调。这比一次性铺大摊子要稳得多。2.2 水平间隔怎么定等间距、非等间距与边界效应水平取几个二水平只能看出线性趋势三水平能看拐点四水平能看更细的形状。工程上三水平是性价比最高的选择——既能判断是否存在极值点实验量又不至于失控。水平间隔的确定有个循序渐进的办法先根据文献或经验确定一个大致范围然后按等间距取三个值。这里有个容易忽略的点——水平的边界不要顶到设备的物理极限或工艺的安全边界。我见过有人把温度水平定在 200、250、300℃结果 300℃ 那一组样品直接碳化了指标完全不可用整张表的分析都要重来。另一个考虑是非等间距。如果初步判断响应在某个区间变化剧烈、在另一个区间很平缓可以有意加密剧烈区间的取点。这种非等间距设计在正交表里是允许的因为正交性只要求每个水平出现次数相同并不要求水平值等间距。代价是极差分析时 R 值的解释会有轻微偏差——R 值反映的是你所选这几个水平之间的响应变化幅度不等间距的时候R 值大的因素未必灵敏度就真的更高它可能只是因为水平跨度定得更宽。2.3 响应指标必须是标量从曲线、光谱、点云里把特征提出来这一节是很多人做正交实验时真正卡住的地方。极差分析、方差分析处理的都是一个数对应一组实验的数据结构但现实里一次实验的输出往往是一条曲线、一张光谱图、一组三维点云。如果你不去处理直接把原始数据塞进分析流程结果一定是乱的。正确的做法是在进入正交分析之前先做一步特征提取把每条曲线、每张图压缩成一个或少数几个标量指标。| 原始数据形态 | 可提取的标量指标 | 提取时要注意的点 | | 时序曲线力、电流、温度 | 峰值、峰位、上升时间、积分面积 | 采样率和基线要统一否则组间不可比 | | 光谱曲线 | 特征峰强度、峰位偏移、峰面积比 | 先做基线校正和归一化 | | 三维点云 | 表面粗糙度 Ra、平面度、体积、缺陷占比 | 配准和坐标系要先统一 | | 图像 | 灰度均值、缺陷面积比、边缘锐度 | 光照条件必须一致 |我个人的经验是特征提取这一步的透明性比技巧性更重要。你提取了什么指标、用什么公式算的、参数怎么设的一定要写清楚。因为后面一旦发现某个因素不显著你得能分清是真的没有效应还是被特征提取过程抹掉了。我自己就遇到过一次提取峰面积时用了固定的积分窗口而某个因素的效应恰好是让峰发生位移位移之后峰跑出了窗口面积反而变小——最后得出了完全相反的结论白白浪费了一个 L18 的实验量。2.4 望目、望大、望小三种特性下的 S/N 比换算如果你用的是田口方法体系还要处理信噪比的问题。指标的期望方向不同处理方式也不同。望大特性越大越好比如强度、附着率信噪比 -10 × log10( (1/n) × Σ(1/y²) )望小特性越小越好比如磨损量、杂质含量信噪比 -10 × log10( (1/n) × Σ y² )望目特性越接近目标值越好比如尺寸、电阻信噪比 10 × log10( ȳ² / S² )其中 S² 是这批数据的样本方差这三个公式本质上都是在把均值和波动揉成一个数。我一般会把信噪比和原始均值分别做一轮极差分析两边对照着看。如果某个因素在信噪比上很显著、在均值上不显著说明它主要是影响稳定性而不是影响水平这种因素在新工艺导入阶段价值很大。2.5 选表与表头设计交互作用列怎么摆这一步是设计阶段技术含量最高的地方。交互作用不是随便忽略的得先判断两个因素之间物理上有没有理由互相影响。温度和时间在固化反应里几乎必然有交互两个毫不相干的助剂之间通常可以假定无交互。以 L9(3^4) 为例它的交互作用表是这样的| 因素所在列 | 交互作用占用列 | | 1 × 2 | 3、4 | | 1 × 3 | 2、4 | | 1 × 4 | 2、3 | | 2 × 3 | 1、4 | | 2 × 4 | 1、3 | | 3 × 4 | 1、2 |这张表藏着一个非常关键的结论在 L9(3^4) 上考察一对交互作用会直接吃掉两列。也就是说如果你把因素 A 放在第 1 列、因素 B 放在第 2 列A×B 的交互作用会占据第 3 列和第 4 列这时候你根本没有地方再放因素 C。这是一个设计上的硬约束。想要在考察交互作用的同时安排 3 个以上的因素必须跳表。L27(3^13) 有 13 列总自由度 26安排 4 个因素8 自由度加上两对交互作用每对占 4 个自由度还有富余的误差自由度这才是正确的选择。很多人硬在 L9(3^4) 上塞 3 个因素然后拿第 4 列当空列估误差这在不考虑交互作用的前提下勉强说得过去但一旦 A 和 B 之间存在交互A×B 的效应会部分混进第 3 列也就是因素 C 那一列C 的显著性判断就是错的。3. 极差分析k 值、R 值和那张最容易读错的趋势图3.1 k 值背后的整齐可比一个手算实例极差分析的运算本身很简单但每一步的含义要清楚。我用一个涂层附着力的例子来走一遍。因素和水平是这样的| 水平 | A 固化温度℃ | B 固化时间min | C 助剂用量份 | | 1 | 150 | 20 | 1.0 | | 2 | 170 | 30 | 1.5 | | 3 | 190 | 40 | 2.0 |按 L9(3^4) 安排 A、B、C 三个因素第 4 列留空。实验数据附着力越大越好如下| 试验号 | A | B | C | 空列 D | 附着力 y | | 1 | 1 | 1 | 1 | 1 | 62 | | 2 | 1 | 2 | 2 | 2 | 71 | | 3 | 1 | 3 | 3 | 3 | 68 | | 4 | 2 | 1 | 2 | 3 | 75 | | 5 | 2 | 2 | 3 | 1 | 82 | | 6 | 2 | 3 | 1 | 2 | 70 | | 7 | 3 | 1 | 3 | 2 | 80 | | 8 | 3 | 2 | 1 | 3 | 73 | | 9 | 3 | 3 | 2 | 1 | 79 |对因素 A 的第 1 水平对应试验号 1、2、3指标之和是 627168201平均值 k167.00。第 2 水平对应试验号 4、5、6和是 758270227k275.67。第 3 水平对应 7、8、9和是 807379232k377.33。这三个平均值反映的是在 B、C 的水平均匀搭配的情况下A 取不同水平时指标的平均表现。因为正交性保证了这一点所以 k1、k2、k3 之间的差异可以归因于 A 本身而不是 B、C 的偏态分布。这就是整齐可比的实际含义。同法算出 B 的 k 值是 72.33、75.33、72.33C 的 k 值是 68.33、75.00、76.67空列 D 的 k 值是 74.33、73.67、72.00。3.2 R 值排序定主次它只是个粗筛别当成判决书R 值就是同一因素各水平平均值中的最大值减最小值。算出来| 因素 | k1 | k2 | k3 | R | 主次 | | A 固化温度 | 67.00 | 75.67 | 77.33 | 10.33 | 1 | | B 固化时间 | 72.33 | 75.33 | 72.33 | 3.00 | 3 | | C 助剂用量 | 68.33 | 75.00 | 76.67 | 8.33 | 2 | | D 空列 | 74.33 | 73.67 | 72.00 | 2.33 | — |R 值排序给出主次顺序A C B。注意空列的 R 值是 2.33比 B 略小。这说明什么说明因素 B 的效应水平和纯随机波动差不多大——如果空列的 R 值超过了某个真实因素的 R 值那这个因素的效应基本上可以判定为噪声。提示空列的 R 值是一个非常重要的参照系。它是你在这张表上的噪声水位线。任何 R 值低于空列 R 值的因素都不要急着解释它的主次先怀疑是误差。但极差分析有一个致命的短板它给不出显著性判断。R 值大的因素就一定显著吗不一定。因素的 R 值大小和试验次数、误差大小直接相关同样的效应量在误差大的体系里会被淹没。要判断这个效应是不是真的存在必须转到方差分析。另外极差分析默认所有因素之间的效应是可加的、没有交互作用。如果存在交互R 值排序可能完全是错的——比如 A 单独看效应很小但它和 B 组合起来效应巨大极差分析会把 A 排在后面。这是它和方差分析最本质的区别。3.3 趋势图的三种形态单调、拐点与水平顶到边界的警告把每个因素的 k1、k2、k3 画成折线图横轴是水平纵轴是指标平均值能看出三种典型形态。第一种是单调递增或递减。这说明你所选的水平区间还没有覆盖到极值点真正的优水平可能在你选的最高水平之外。这时候要么直接选边界水平要么补一组区间外的实验确认。第二种是中间高两头低。这就是存在极值点的信号此因素的优水平应该取中间那个值或者进一步在中间值附近加密水平做精细化实验。第三种是变化平缓。R 值小、趋势平这类因素通常在方差分析里也不显著可以考虑固定在一个成本更低或操作更方便的水平上不必在它身上花精力。我最想强调的是第一条的警示作用。我见过一个案例因素水平的第三个值是 190℃趋势还在往上走团队直接选了 190℃ 作为最优。结果中试放大的时候温度控制精度不够实际波动到 195℃涂层开始出现微裂纹。如果当时在 190℃ 和 200℃ 之间补一组实验就能提前看到拐点不会在设计阶段埋这个雷。趋势图不平坦地顶到边界永远值得再补一组验证。4. 方差分析把误差自由度从哪里挤出来4.1 空列法与重复实验法两种截然不同的误差来源方差分析的第一步是把总平方和拆开一部分归因于各因素的效应剩下那部分归因于误差。麻烦在于误差平方和不会自己跳出来你得想办法把它挤出来。空列法是在正交表上留一两列不安排任何因素默认这些列的离差完全由随机误差贡献把它当作误差项。上面的 L9(3^4) 例子中第 4 列 D 就是空列它的平方和就是误差平方和。这个办法的优点是零额外成本缺点是误差自由度非常小——三水平因素的单个空列只提供 2 个自由度。重复实验法是每个试验号做多次通常 2 到 3 次用组内变异作为误差。它的误差自由度 试验号数 ×重复次数 - 1信息量大得多而且能捕捉到试验顺序带来的漂移。代价是实验量翻倍。我个人的取舍原则是能重复就重复尤其是在效应量预期比较小的时候。空列法的误差自由度只有 2 时F 检验的临界值高得离谱很多真实效应会被判为不显著。这个坑我在第 4.3 节会用一个具体数字说明。如果实在不能重复至少留两个空列并把它们合并成误差项前提是这两列的平方和差异不大可以用 F 检验先检验一下它们是否同质这样误差自由度能到 4检验功效会好很多。4.2 平方和、自由度、均方、F 值完整手算流程把公式列清楚这样你用任何工具实现的时候都能对得上。总平方和SS_T Σ(y_i - ȳ)²其中 ȳ 是所有数据的平均值自由度 df_T N - 1N 是试验总次数。因素平方和SS_j (Σ T_ij²) / n_i - (Σ y)² / N。这里 T_ij 是第 j 个因素第 i 个水平下所有指标之和n_i 是该水平出现的次数均衡设计时每个水平次数相同最后那一项叫校正项 CT (Σy)²/N。因素自由度df_j 水平数 - 1。误差平方和SS_e SS_T - Σ SS_j误差自由度df_e df_T - Σ df_j。均方MS SS / dfF 值F_j MS_j / MS_e。判断显著性时查 F 分布表的临界值。三水平因素、误差自由度不同的情况下临界值是| 误差自由度 df_e | F0.05(2, df_e) | F0.01(2, df_e) | | 2 | 19.00 | 99.00 | | 4 | 6.94 | 18.00 | | 6 | 5.14 | 10.92 | | 8 | 4.46 | 8.65 | | 12 | 3.89 | 6.93 |这张表最值得盯的是第一行和第二行的差距。误差自由度从 2 提到 40.05 水平的临界值从 19.00 掉到 6.94接近三倍。也就是说同样的因素均方用空列法可能判不显著用重复实验法就是显著的。这不是数据变了是你手里的尺子精度变了。4.3 用我那个涂层数据走一遍A 显著而 C 卡在临界值下面说明了什么继续用第 3 节的数据。总指标和是 660N9CT 660²/9 48400总平均 ȳ 73.33。各因素的平方和A(201² 227² 232²)/3 - 48400 48584.67 - 48400 184.67B(217² 226² 217²)/3 - 48400 48418.00 - 48400 18.00C(205² 225² 230²)/3 - 48400 48516.67 - 48400 116.67空列 D(223² 221² 216²)/3 - 48400 48408.67 - 48400 8.67总平方和等于 Σ(y - 73.33)²算出来约 328.00。校验一下184.67 18.00 116.67 8.67 328.01对得上。自由度方面每个因素 df2总 df8误差 df 8 - 6 2。误差平方和取空列 D 的 8.67均方 MS_e 4.335。| 来源 | SS | df | MS | F 值 | 临界值 F0.05(2,2) | 结论 | | A | 184.67 | 2 | 92.34 | 21.30 | 19.00 | 显著 | | C | 116.67 | 2 | 58.34 | 13.46 | 19.00 | 不显著 | | B | 18.00 | 2 | 9.00 | 2.08 | 19.00 | 不显著 | | 误差 e | 8.67 | 2 | 4.34 | — | — | — |这张表里最值得琢磨的是 C。它的 R 值排第二8.33比 B 大得多趋势也很明确68.33 → 75.00 → 76.67 单调上升但 F 值 13.46 卡在临界值 19.00 下面按 0.05 水平判不显著。这是数据不显著还是尺子不够好我认为是后者。误差自由度只有 2MS_e 本身估计得不稳F 检验的功效被压得很低。如果这个因素是工艺上必须调的一个参数我会选择补做重复实验每个试验号重复 3 次误差自由度变成 9×218这时候 F0.05(2,18)3.55C 的显著性判断就完全不一样了。这也是我给所有做正交实验的人的第一个建议在算方差分析之前先看看你的误差自由度有多少。低于 4 的话不要太相信 F 检验的结论。极差分析在这里反而更有参考价值——它虽然不给显著性但对效应大小的排序是稳健的。4.4 预测值、预测区间与验证实验确定了显著因素和各自的优水平之后可以预测最优组合下的指标值。公式是预测值 ȳ Σ(k_优 - ȳ)拿上面的例子优水平是 A3、B2、C3这三个水平的 k 值分别最大。预测值 73.33 (77.33-73.33) (75.33-73.33) (76.67-73.33) 73.33 4.00 2.00 3.34 82.67。注意如果某个因素经检验不显著原则上可以不纳入预测计算直接把它的效应当作噪声如果出于成本或工艺原因必须固定某个水平那就用实际固定的水平去算。预测区间会更宽一些大约是预测值 ± t(α/2, df_e) × sqrt(MS_e × (1 因素效应所占的自由度数)/N 之类的修正系数)。在这个例子里df_e2t(0.025, 2)4.303算出来的半宽接近 6预测区间大约落在 77 到 89 之间。这个区间宽得几乎没什么实用价值——但它本身就是一条很重要的信息你的误差自由度太小导致你的预测精度差。这时候不要急着下结论说最优组合的指标就是 82.67而要去做验证实验。验证实验做 3 次看实测均值是否落在这个区间内。如果落在里面说明模型可信如果明显超出说明有未考虑到的交互作用或系统误差需要回头检查表头设计。5. Python 全流程实现建表、极差分析、方差分析一气呵成5.1 依赖与选型取舍为什么我不用 pyDOE2 生成三水平表做正交实验的 Python 生态里pyDOE2 是绕不开的一个库。它提供了fullfact全因子、fracfact两水平部分因子、pbdesignPlackett-Burman 设计、bbdesignBox-Behnken等函数。但这里有个实际使用中会撞上的情况pyDOE2 并没有直接生成三水平正交表L9、L27 这类的接口。它偏向两水平的筛选设计因为那类设计在工业界用得最多。如果你需要 L9(3^4) 或者 L18(2^1×3^7) 这样的三水平表最省事的做法是把标准表以数组的形式直接写进代码。这不是偷懒。正交表是有限个标准化的表格工业界用了几十年就那么几十张硬编码反而更可控——不会有版本升级改接口的风险也方便你在注释里标清楚水平对应的实际物理值。pip install numpy pandas statsmodels matplotlib5.2 用 numpy 硬编码标准正交表并拼装试验记录表import numpy as np import pandas as pd # 标准 L9(3^4) 正交表每一行是一组试验每一列是一个因素位 L9 np.array([ [1, 1, 1, 1], [1, 2, 2, 2], [1, 3, 3, 3], [2, 1, 2, 3], [2, 2, 3, 1], [2, 3, 1, 2], [3, 1, 3, 2], [3, 2, 1, 3], [3, 3, 2, 1], ], dtypeint) # 水平到实际物理量的映射表 level_map { A_温度: {1: 150, 2: 170, 3: 190}, B_时间: {1: 20, 2: 30, 3: 40}, C_助剂: {1: 1.0, 2: 1.5, 3: 2.0}, } df pd.DataFrame(L9, columns[A, B, C, D]) df.insert(0, run, np.arange(1, 10)) # 生成给实验员看的操作版记录表直接给出物理量 for col, mapping in level_map.items(): key col.split(_)[0] df[f{col}(实际)] df[key].map(mapping) # 指标列先留空实验做完再回填 df[y] np.nan print(df.to_string(indexFalse))这段代码有两个设计上的用意。一是我保留了原始的水平编号列A、B、C、D因为后面的极差分析和方差分析都要基于水平编号做分组而不是基于物理量——物理量不等间距的时候用编号分组更不容易出错。二是额外生成了一个实际值列给实验员用避免他们拿着 1/2/3 去猜对应多少度。回填数据的时候直接赋值df[y] [62, 71, 68, 75, 82, 70, 80, 73, 79]5.3 极差分析的向量化写法def range_analysis(df, factor_cols, respy): 极差分析返回各因素的 k 值、R 值和优水平 records [] for col in factor_cols: grouped df.groupby(col)[resp] k grouped.mean() # 各水平指标均值 records.append({ 因素: col, k1: k.get(1, np.nan), k2: k.get(2, np.nan), k3: k.get(3, np.nan), R: k.max() - k.min(), 优水平: int(k.idxmax()), # 望大特性取最大 }) out (pd.DataFrame(records) .sort_values(R, ascendingFalse) .reset_index(dropTrue)) out.insert(0, 主次, np.arange(1, len(out) 1)) return out print(range_analysis(df, [A, B, C, D]))跑出来的结果和手算一致A 的 R10.33 排第一C 的 R8.33 排第二B 的 R3.00空列 D 的 R2.33。如果指标是望小特性把idxmax改成idxmin就行。如果用的是信噪比就先把 y 换算成 S/N 值再传进这个函数逻辑完全一样。画趋势图也很直接import matplotlib.pyplot as plt fig, axes plt.subplots(1, 3, figsize(12, 3.5)) for ax, col in zip(axes, [A, B, C]): k df.groupby(col)[y].mean() ax.plot(k.index, k.values, markero) ax.set_title(f因素 {col}) ax.set_xlabel(水平) ax.set_ylabel(指标均值) ax.set_xticks([1, 2, 3]) ax.grid(alpha0.3) plt.tight_layout() plt.savefig(trend.png, dpi150)5.4 手算方差分析 statsmodels 交叉验证我习惯把方差分析写两遍一遍手算一遍用 statsmodels两边对得上才放心用。def anova_manual(df, factor_cols, respy): y df[resp].values N len(y) ct y.sum() ** 2 / N # 校正项 ss_total ((y - y.mean()) ** 2).sum() rows [] for col in factor_cols: g df.groupby(col)[resp] n_i g.size().iloc[0] # 每个水平的重复次数 ss (g.sum() ** 2).sum() / n_i - ct rows.append({来源: col, SS: ss, df: g.ngroups - 1}) ss_err ss_total - sum(r[SS] for r in rows) df_err N - 1 - sum(r[df] for r in rows) ms_err ss_err / df_err for r in rows: r[MS] r[SS] / r[df] r[F] r[MS] / ms_err rows.append({来源: 误差e, SS: ss_err, df: df_err, MS: ms_err, F: np.nan}) return pd.DataFrame(rows)[[来源, SS, df, MS, F]] print(anova_manual(df, [A, C, B]).round(3))注意这里我把空列 D 排除在因素之外误差项就来自总平方和减去 A、B、C 三者的残差结果正好等于 D 列的平方和 8.67。再用 statsmodels 交叉验证一遍from statsmodels.formula.api import ols from statsmodels.stats.anova import anova_lm model ols(y ~ C(A) C(B) C(C), datadf).fit() print(anova_lm(model, typ2))用C()把数值列声明为分类变量很关键。如果不加statsmodels 会把 A 的 1/2/3 当成连续变量做线性回归自由度变成 1 而不是 2出来的 F 值完全不是一回事。三水平因素在模型里会被编码成 2 个哑变量所以 ANOVA 表里的 df2和手算一致。有一点要提醒L9 排 3 个因素的情况下残差自由度只有 2。这时候 statsmodels 算出来的 F 值和手算完全一样但不要因为跑出来了就认为结论可靠。第 4.3 节讲过的问题在代码里一样存在。6. Excel 与 MATLAB 两条替代路线各自的主场在哪6.1 ExcelAVERAGEIF 与数据透视表的高效组合先说清楚 Excel 的主场实验数据量在几十到几百行、需要频繁给人看、团队里其他人不会写代码的场合。这种场景下Excel 的沟通效率比 Python 高一个数量级——你把表发过去对方点开就能看懂。极差分析在 Excel 里基本不用手算。k 值用一个AVERAGEIF就够了AVERAGEIF($B$2:$B$10, 1, $F$2:$F$10)这里 B 列是因素 A 的水平号F 列是指标。把条件从 1 改成 2、3就得到 k1、k2、k3。R 值直接MAX(...)-MIN(...)。更省事的办法是用数据透视表把因素列拖到行指标列拖到值并设为平均值一秒出结果。四列因素就是四张透视表或者用多重合并计算区域拼在一起。有个小技巧能省掉大量重复劳动把列标题做成变量。用MATCH函数找到因素列的位置然后配合INDEX动态生成 AVERAGEIF 的列引用这样换一张正交表只需要改一处公式不用重写。我做过一次 L27(3^13) 的分析13 列因素如果每个都手写公式光是调试引用就要半小时。6.2 Excel 里最常见的三类公式事故第一类是引用错位。横向拖公式的时候本该锁定行的地方没加$结果第二行开始全部算错。判断方法很简单把任意一个 k 值手算一遍和公式结果对一下不一致就是引用问题。第二类是把水平号当数值参与排序。数据透视表默认对行标签按升序排1、2、3 没问题但如果你的水平号写成了低中高这种文本排序顺序就是按拼音来的不一定符合你的预期最好加一个辅助的排序列。第三类是公式覆盖范围没有跟着数据更新。新补做了一组验证实验加在第 10 行但 AVERAGEIF 的范围还是$B$2:$B$10新数据没被纳入。这个错误特别隐蔽因为结果看起来是正常的只是略微偏了一点。我的做法是把数据区转成表格CtrlT公式范围会自动扩展。6.3 MATLAB anovan一次把主效应和交互作用都算出来如果实验里有重复测量或者你需要同时考察主效应和交互作用MATLAB 的anovan是最省事的工具。% A、B、C 是三个因素的水平号向量长度 9 % y 是指标向量 [p, tbl, stats] anovan(y, {A, B, C}, ... model, linear, ... % 只算主效应不带交互 varnames, {温度, 时间, 助剂}, ... display, on);anovan返回的tbl单元格数组里直接给出了平方和、自由度、均方、F 值和 p 值不用自己算。model参数是关键——默认值是full会把三阶交互也算进去在 L9 这种小表上会直接报错说自由度为负。改成linear只算主效应或者写interaction算二阶交互并手动排除三阶[p, tbl, stats] anovan(y, {A, B, C}, ... model, [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1], ... varnames, {温度, 时间, 助剂}, ... display, on);model矩阵的每一行是一个效应项列对应因素。第一行[1 0 0]表示温度的主效应第四行[1 1 0]表示温度和时间的交互作用。这个写法比字符串方式灵活尤其是因素数量多的时候。MATLAB 还提供了fracfact生成两水平部分因子设计和fracfactgen生成对应的生成元这一点和 pyDOE2 的定位类似都是针对两水平筛选设计的。6.4 R 与 Minitab、Design-Expert 的定位如果你更习惯 RDoE.base包里的oa.design可以直接生成正交表library(DoE.base) design - oa.design(nfactors 3, nlevels 3, factor.names c(A, B, C)) summary(design)它会自动选一张最小的可行正交表并完成表头设计。FrF2包则专注于两水平部分因子和分辨度分析做筛选实验的时候很好用。商业软件方面Minitab 的正交设计模块和 Design-Expert 的响应面模块在这块做得非常成熟尤其是交互作用图和等高线图的自动生成省掉了很多画图代码。如果项目预算允许、团队里没人写代码用它们是完全合理的选择。我用过的体会是设计阶段用商业软件省时间数据处理和分析阶段的自动化用 Python 更可控。两者不是替代关系很多人把它们对立起来其实没必要。7. 复盘六个真实踩过的坑7.1 交互作用被挤进空列一个隐蔽的混杂这个坑我在第 2.5 节埋了伏笔这里展开说。L9(3^4) 上排 A、B、C 三个因素第 4 列当空列估误差。如果 A 和 B 之间实际上存在交互作用而这个交互作用的效应会分别落入第 3 列和第 4 列参照前面的交互作用表1×2 的交互占 3、4 列。后果是双重的。第 4 列本来应该只含纯误差现在混进了交互效应误差平方和被高估F 检验的灵敏度进一步降低第 3 列本来代表因素 C 的主效应现在混进了交互效应C 的显著性判断可能是假的。你在报告里写C 因素显著实际上可能是 A 和 B 的联合效应在起作用。排查的办法是对显著因素做追加验证实验把 A 和 B 分别固定在各自的两个不同水平上其他因素不变做成一个 2×2 的小实验。如果 A 的效应在 B 的不同水平下差异明显交互作用就确认存在。这时候原来那张 L9 的结论要作废重新用 L27(3^13) 或 L18(2^1×3^7) 做一轮。7.2 多指标打架时综合评分法的权重从哪来一个实验测三个指标A 因素让指标 1 变好、指标 2 变差怎么选优水平标准做法是综合评分把每个指标归一化到 0~1 区间然后按权重加权求和得到一个综合分再做极差分析。麻烦在于权重。我在项目里见过三种确定权重的方式各有适用场景| 方式 | 具体做法 | 适用场景 | 风险 | | 主观赋权 | 项目组集体讨论定权重 | 指标之间重要性差异明确 | 容易被话语权大的人主导 | | 变异系数法 | 按各指标的标准差占比赋权 | 数据驱动无需先验判断 | 波动大的指标会被过度加权 | | 熵权法 | 按指标的信息熵赋权 | 客观适合样本量较大的情况 | 需要较大样本小实验里不稳 |归一化的时候要注意方向。望大、望小、望目的归一化公式不一样如果搞混了综合分算出来是反的。而且归一化的参考范围用本批数据的最大最小值还是用工程允许的上下限会显著影响结果。我一般用工程允许范围而不是本批数据的极值因为后者会让这一批数据里最好的一个得满分掩盖了它可能根本达不到要求的事实。7.3 试验顺序随机化不是形式主义正交表的试验号顺序是固定的 1 到 9但实际执行顺序必须随机化。很多人觉得反正这 9 组都要做先后无所谓——这个想法会在两种情况下出问题。第一种是设备漂移。如果设备在实验过程中有缓慢的升温趋势或者刀具磨损按表顺序执行的话前几组和后几组之间的差异里就混进了时间效应而这个效应会被错误地归因到排在前面的因素上。L9 的标准表里第 1 列的第 1 水平恰好出现在前 3 行如果时间效应和 A 的效应方向一致A 的 R 值会被放大。第二种是操作员的学习效应。前几组做得生疏后面越做越熟练指标系统性变好。这个问题在手工操作占比较高的实验里特别明显。解决办法很直接用一个随机数生成器打乱试验号并且把实际执行顺序记录在案。分析的时候依然按正交表原编号分组顺序信息只用于事后排查异常——如果发现指标和实际执行顺序有明显相关性就该考虑加一个批次作为区组因素重新分析。7.4 重复实验做成了平行样而不是独立重复重复实验的价值取决于重复的独立性。同一个样品测三次得到的是测量重复性分别制备三个样品各测一次得到的是制备测量的综合波动。这两者反映的误差来源完全不同。如果你的实验目的是判断这个工艺组合在实际生产中能不能稳定复现那必须用后一种重复方式。前一种方式算出来的误差方差会严重偏小让所有因素都变得显著看起来结果很漂亮但放大到生产就崩了。这个坑我发现得比较晚。早年做的一批实验误差方差小得异常F 值个个都很高当时还挺高兴。后来同事提醒我看原始记录发现重复的三组是同一次配料的样品——等于把测量噪声当成了工艺噪声。重新按独立制备做了重复之后显著性结果有一半变了。7.5 把极差分析选出的优组合直接当结论极差分析给出的优组合是各个因素最优水平的拼接而这个拼接出来的组合很可能根本没在你的 9 组实验里出现过。它是在效应可加这个假设下外推出来的预测不是观测值。这就是为什么验证实验不能省。我一般会在报告里同时给出三个数字极差分析预测的优组合、你实际做过的 9 组里表现最好的那一组、验证实验的实测值。如果验证值和预测值差得远说明有交互作用或者非线性效应这时候就得回到设计阶段重新考虑。还有一种情况也要留意优组合的预测值虽然高但成本也高得离谱。比如温度取最高档能耗成本翻倍。这种情况下最优要重新定义变成一个带约束的优化问题而不是单纯看指标均值。工程上的最优从来不只是一个统计概念。7.6 数据搬来搬去导致的记录错位最后说一个看似低级但发生频率极高的问题数据在多个文件之间传递时错位。典型场景是实验员在纸质记录本上写数据助理录进 Excel分析师把 Excel 的数据复制到 Python 脚本里。每经过一次手动搬运就多一次错位的风险。L9 的 9 行数据里如果第 5 行和第 6 行搞反了极差分析的结果会有轻微变化但往往不会到一眼看出不对的程度——直到某天你重做实验验证发现复现不了。我现在坚持三个做法。第一正交表的试验号在纸质记录、Excel、脚本里必须完全一致不做任何重编号。第二原始记录只转录一次后续所有工具都从同一个源头读取不要各自复制一份。第三分析脚本跑完先打印一张试验号-指标值的对照表肉眼扫一遍确认和纸质记录对得上再往下算。这个流程听起来笨但省下来的返工时间远超过它花掉的几分钟。批量录入的活儿如果量很大用一点轻量的自动化脚本去解析固定格式的记录文件也是个办法——关键是在每个环节都留下可以对照的中间产物而不是让数据在流水线上无声地变形。正交实验设计这套方法本身不复杂难的是每一环都不打折扣地执行。因素水平的确定要认真调研正交表的选型要考虑交互作用和误差自由度极差分析要配合趋势图一起看方差分析要注意误差自由度够不够最后一定要有验证实验收尾。我个人的体会是这五个环节里任何一个打折最后都会在放大或量产阶段以某种方式还回来——可能是复现不了可能是某个因素的显著性判断错了也可能是最优参数其实一直在一个危险的边界上。数据的价值不在于它有多漂亮而在于每一步推导都能被追溯、被复现。