与“六西格玛绿带手册”DOE回归建模商榷

◆俞钟行/文

《六西格玛绿带手册》[1]第6章对一个全因子设计实例用回归建模取得最优解。本文提出不同的方法，可以获得更合理、更精确的答案。

1.原例简介

为降低加工的平面度，研究3个因子各2个水平，又设0水平，见表1。

实际做了23+4次试验，具体方案和结果见表2。

及它们的交互作用AB是显著的，继而对这3项回归建模，得到回归方程：

Y=-219.412+0.267A+109.625B-0.116250AB（使用未编码单位的数据）

最后用MINITAB的响应优化器，得因子最佳组合为A=900、B=2.8，预测最佳值y=32.8875。

原例讨论说：表1中运行序8的试验已达到32.3，为什么预测最佳值明显地比它大？这就还要看运行序4的试验结果为33.5，这两次试验的条件只差一个非显著因子C（倾角），况且（32.3+33.5）/2=32.9，很接近预测最佳值。

《六西格玛绿带手册》的242页倒数第4行起有 一段话：“要注意的是：如果一个高阶项（交互作用或二次项）是显著的，则此高阶项所包含的主因子也必须包含在模型中。例如，二阶交互作用BC项显著，则B及C这两个主效应也一定要被包含在模型中，即使表面上看这两个主效应项本身并不显著。”这段话对回归建模有影响，无以名之，姑且称其为“金科玉律”。

2.是否一定要对因子编码

3因子2水平是很普通的试验，即使加几个中心点也如此。原例用MINITAB软件作分析，其实用excel也可以。原例强调，回归建模前必须先对因子编码，即高水平=1，低水平=-1，0水平=0。试以原例最后得到的回归方程加以说明：其实编码不是必须的。若以表3所示的excel电子表格，这里因子并未编码，用excel“数据分析”中“回归”模块分析，可得图1所示结果。图中最下表的Coefficients这列，就给出该回归方程的各项系数，与前述相同。而且，《六西格玛绿带手册》第249页图中给出原例的S=1.10997、R-Sq=88.11%、R-Sq=83.66%等值，也与图1最上表相同。在excel中直接以原始数据做回归分析，显然是简捷可靠的，至少省略了从编码后获得的回归方程再返回到原始数据回归方程的麻烦。

表1 因素水平表

表2 试验方案和结果

表3 原例最佳方程用excel做回归分析的界面

3.因子C（倾角）是非显著因子吗

对表2可以用excel画出因素趋势图如图2。

各因子的趋势图有显著的交叉，提示它们之间可能有显著的交互作用。实际上以表4所示的excel电子表格界面，用“回归”分析，可以得到拟合很好的回归方程，而且因子C是显著的（图3）。

从图3的最上表看，诸如S=1.070933、R-Sq=90.3192%、R-Sq=84.7873%等值都比原例好；最下表的4项因子或交互作用的p值也全都小于0.05。

图1 原例最佳方程的回归结果（直接用原始数据）

图2 因素趋势图（自左向右为A、B和C）

表4 参考图2的excel电子表格界面

图3 由表4所得的回归结果

若看《六西格玛绿带手册》的246页所示图，当对原例3个因子及其2阶交互作用进行回归建模时，标准误差S=1.23956，故因子C、交互作用BC都不是显著的。但在图3所示情况下，标准误差S降为1.070933，C和BC成为显著就不足为奇了。

4.“金科玉律”对吗

按照“金科玉律”，图3所得的回归方程是不可取的。因为此方程里有2阶交互作用BC和AB，却没有主效应B。更何况从图2看，B是最强的因子。

有些人认为“金科玉律”来自“效应排序原则”和“效应遗传原则”。“效应排序原则”（i）为“低阶效应应比高阶效应更重要”。根据这条，似乎可依逻辑推得：2阶交互作用没有相关的主效应重要，所以既然有2阶交互作用，就必须保留相关的主效应，而这就是“金科玉律”所述。但是，这条原则接着解释：“在因子效应的数目较大而不能全部进行估计时特别有效，这是一个经验原则”[2]。应当说现在已有丰富的经验证明，像图3那样判断效应显著性是合理与可行的，如“六西格玛管理”（第二版）例7-7就是。

“效应遗传原则”为“要使一个交互作用是显著的，至少它的一个亲本因子应该是显著的”[2]。按这个原则看，回归方程里没有显著的主效应不一定不允许，这个原则是针对交互作用项的。而且用生活常识作比喻：肉、菜和萝卜熬成的汤里肉看不见了，但菜、萝卜和汤里充满了肉味，这是可以的吧。

总之，应当认为“金科玉律”的要求并不合理，它会导致回归建模受到不应有的束缚，影响改进。有的MINITAB专家说明：“MINITAB15以前不能删除交互作用重要而主效应不重要的主效应项的，从16版本开始就增加了一个选项，由使用者自行决定删不删除这样的主效应项，这一改进带来了灵活性，方便读者使用！”还有知情方说明：金科玉律“是在用MINITAB的DOE模型构建时所必须的”，如果“采用GLM（广义线性模型）做的，当然不受刚才这句话（指金科玉律）的约束了”。因此，我们应当呼吁：让原例也享受采用GLM（广义线性模型）的待遇，如“六西格玛管理”（第二版）例7-7所示，否则就有削足适履之谬。

5.对残差的正态性检验

对图3最下表获得的回归方程y=74.325+0.258058A+1.301116 BC-3.48915C-0.11281AB，用excel的“规划求解”选优，得到最佳组合为A=900、B=2.4、C=84，即表2中运行序8的试验。此时预测值y=32.12218，比实测值32.3还小些，可视作误差。

那么这个回归方程是否真的拟合得很好，除了图3表明的各项参数外（包括中间表的Significance F=0.001175、残差占比=8.028286/82.93=9.68%），下面以W检验法（又名“夏皮罗-威尔克正态检验法”）[3]对它的残差作检验，以增加信心。这个检验当8<=n<=50时可以利用，适合本例n=12。表5是残差计算表，残差的平均值=1.0445E-12。

这个检验是建立在次序观测值的基础上，将残差按升序排列，记为x(1),x(2),…x(n),然后计算：

这里k=n/2=6(n为偶数时)，系数ak从文献[3]查得。S2的计算见表6。

表5 残差计算

表6 S2的计算表