6.2 2^{2} 设计 - 试验设计与分析中文翻译

$2^{k}$ 系列中的第一个设计只有两个因子，比如说 A 和 B，每个因子都在两个水平上运行。这种设计称为 2² 因子设计。因子的两个水平可以随意称为"低"和"高"。作为一个例子，考虑一项关于反应物浓度与催化剂用量对某个化学过程转化率（产率）影响的研究。该试验的目的是确定调整这两个因子中的任何一个是否能提高产率。设反应物浓度为因子 A，感兴趣的两个水平为 15% 和 25%。催化剂为因子 B，高水平表示使用 2 磅催化剂，低水平表示只使用 1 磅。试验重复三次，因此共有 12 次试验。试验的进行顺序是随机的，所以这是一个完全随机化试验。所得数据如下：

因子 A	因子 B	处理组合	重复 I	重复 II	重复 III	总计
−	−	A 低，B 低	28	25	27	80
+	−	A 高，B 低	36	32	32	100
−	+	A 低，B 高	18	19	23	60
+	+	A 高，B 高	31	30	29	90

该设计中的四个处理组合如图 6.1 所示。按照约定，我们用大写拉丁字母表示一个因子的效应。因此，"A"指的是因子 A 的效应，"B"指的是

图 6.1 2² 设计中的处理组合

因子 B 的效应，而"AB"指的是 $AB$ 交互作用。在 2² 设计中，A 轴和 B 轴上的"−"和"+"分别表示 A 与 B 的低水平和高水平。因此，A 轴上的 − 表示浓度的低水平（15%），而 + 表示高水平（25%）；B 轴上的 − 表示催化剂的低水平，+ 表示高水平。

设计中的四个处理组合也用相应的小写字母表示，如图 6.1 所示。从图中可以看到，处理组合中某个因子的高水平用对应的小写字母表示，而某个因子的低水平用相应字母的缺失表示。因此， $a$ 表示 A 取高水平、B 取低水平的处理组合， $b$ 表示 A 取低水平、B 取高水平的处理组合， $ab$ 表示两个因子都取高水平。按照约定，(1) 用来表示两个因子都取低水平。这一记号贯穿整个 $2^{k}$ 系列。

在两水平因子设计中，我们可以把一个因子的平均效应定义为：在其他因子的水平上取平均后，该因子水平改变所引起的响应变化。此外，符号 (1)、 $a$ 、 $b$ 和 $ab$ 现在表示该处理组合处全部 $n$ 次重复的响应观测值总和，如图 6.1 所示。于是，B 取低水平时 A 的效应为 $[a - (1)]/n$ ，B 取高水平时 A 的效应为 $[ab - b]/n$ 。把这两个量取平均就得到 A 的主效应：

\begin{array}{r l} A & = \frac{1}{2 n} \{[a b - b] + [a - (1)]\} \\ & = \frac{1}{2 n} [a b + a - b - (1)] \end{array} \tag{6.1}

B 的平均主效应由 A 取低水平时 B 的效应（即 $[b-(1)]/n$ ）和 A 取高水平时 B 的效应（即 $[ab-a]/n$ ）给出，为

\begin{array}{r l} & B = \frac{1}{2 n} \{[a b - a] + [b - (1)]\} \\ & \quad = \frac{1}{2 n} [a b + b - a - (1)] \end{array} \tag{6.2}

我们把交互效应 $AB$ 定义为 B 取高水平时 A 的效应与 B 取低水平时 A 的效应之差的平均。于是

\begin{array}{r} AB = \frac{1}{2 n} \{[a b - b] - [a - (1)]\} \\ = \frac{1}{2 n} [a b + (1) - a - b] \end{array} \tag{6.3}

另外，我们也可以把 $AB$ 定义为 A 取高水平时 B 的效应与 A 取低水平时 B 的效应之差的平均。这同样会导出式 6.3。

A、B 和 AB 效应的公式还可以用另一种方法导出。A 的效应可以由图 6.1 中方框右侧两个处理组合的平均响应（记为 $\overline{y}_{A^{+}}$ ，因为它是 A 取高水平时各处理组合的平均响应）与左侧两个处理组合的平均响应（即 $\overline{y}_{A^{-}}$ ）之差得到。即

\begin{array}{r l} A & = \overline{{y}}_{A^{+}} - \overline{{y}}_{A^{-}} \\ & = \frac{a b + a}{2 n} - \frac{b + (1)}{2 n} \\ & = \frac{1}{2 n} [a b + a - b - (1)] \end{array}

这与式 6.1 的结果完全相同。B 的效应即式 6.2，由方框顶部两个处理组合的平均值（ $\overline{y}_{B^{+}}$ ）与底部两个处理组合的平均值（ $\overline{y}_{B^{-}}$ ）之差得到，即

\begin{array}{r l} & B = \overline{{y}}_{B^{+}} - \overline{{y}}_{B^{-}} \\ & \quad = \frac{a b + b}{2 n} - \frac{a + (1)}{2 n} \\ & \quad = \frac{1}{2 n} [a b + b - a - (1)] \end{array}

最后，交互效应 AB 等于方框中从右到左对角线上两个处理组合 $[ab$ 与 (1) $]$ 的平均值减去从左到右对角线上两个处理组合（ $a$ 与 $b$ ）的平均值，即

\begin{array}{r l} AB & = \frac{a b + (1)}{2 n} - \frac{a + b}{2 n} \\ & = \frac{1}{2 n} [a b + (1) - a - b] \end{array}

这与式 6.3 完全一致。

利用图 6.1 中的试验，我们可以把平均效应估计为

A = \frac{1}{2(3)} (90 + 100 - 60 - 80) = 8.33

B = \frac{1}{2(3)} (90 + 60 - 100 - 80) = - 5.00

AB = \frac{1}{2(3)} (90 + 80 - 100 - 60) = 1.67

A（反应物浓度）的效应为正，这说明把 A 由低水平（15%）提高到高水平（25%）会提高产率。B（催化剂）的效应为负，这说明增加加入到过程中的催化剂量会降低产率。相对于两个主效应而言，交互效应看起来很小。

在涉及 $2^{k}$ 设计的试验中，考察因子效应的大小和方向以确定哪些变量可能重要，始终是很重要的。方差分析一般可以用来确认这一解释（也可以用 $t$ 检验）。效应的大小和方向应当始终与方差分析一起考虑，因为仅靠方差分析并不能传递这些信息。有几种优秀的统计软件包可用于建立和分析 $2^{k}$ 设计。此外，手工计算时还有一些节省时间的特殊方法。

下面考虑 A、B 和 AB 的平方和如何确定。注意由式 6.1 可见，估计 A 时用到了一个对照（contrast），即

\text{Contrast}_{A} = a b + a - b - (1) \tag{6.4}

我们通常把这个对照称为 A 的总效应。由式 6.2 和 6.3 可见，估计 B 和 AB 时也用到了对照。此外，这三个对照是正交的。任何对照的平方和都可以由式 3.29 计算，该式指出：任何对照的平方和等于该对照的平方除以对照中每个总和所含的观测值个数再乘以对照系数的平方和。于是有

SS_{A} = \frac{[a b + a - b - (1)]^{2}}{4 n} \tag{6.5}

SS_{B} = \frac{[a b + b - a - (1)]^{2}}{4 n} \tag{6.6}

以及

SS_{AB} = \frac{[a b + (1) - a - b]^{2}}{4 n} \tag{6.7}

分别作为 A、B 和 AB 的平方和。注意这些式子是多么简单：只需对一个数取平方就能算出平方和。

利用图 6.1 中的试验，由式 6.5 至 6.7 可得各平方和为

SS_{A} = \frac{(50)^{2}}{4(3)} = 208.33

SS_{B} = \frac{(- 30)^{2}}{4(3)} = 75.00 \tag{6.8}

以及

SS_{AB} = \frac{(10)^{2}}{4(3)} = 8.33

总平方和按通常方式求得，即

SS_{T} = \sum_{i = 1}^{2} \sum_{j = 1}^{2} \sum_{k = 1}^{n} y_{i j k}^{2} - \frac{y_{\cdots}^{2}}{4 n} \tag{6.9}

一般地， $SS_{T}$ 有 $4n - 1$ 个自由度。误差平方和通常按相减的方式计算，它有 $4(n - 1)$ 个自由度，即

SS_{E} = SS_{T} - SS_{A} - SS_{B} - SS_{AB} \tag{6.10}

对图 6.1 中的试验，我们得到

\begin{array}{r l} SS_{T} & = \sum_{i = 1}^{2} \sum_{j = 1}^{2} \sum_{k = 1}^{3} y_{i j k}^{2} - \frac{y_{\dots}^{2}}{4(3)} \\ & = 9398.00 - 9075.00 = 323.00 \end{array}

以及

\begin{array}{r l} SS_{E} & = SS_{T} - SS_{A} - SS_{B} - SS_{AB} \\ & = 323.00 - 208.33 - 75.00 - 8.33 \\ & = 31.34 \end{array}

其中 $SS_{A}$ 、 $SS_{B}$ 和 $SS_{AB}$ 取自式 6.8。完整的方差分析汇总于表 6.1。根据 $P$ 值，我们得出主效应在统计上显著、而这两个因子之间不存在交互作用的结论。这证实了我们最初基于因子效应大小对数据所作的解释。

把处理组合按 (1)、 $a$ 、 $b$ 、 $ab$ 的顺序写下来往往很方便。这称为标准顺序（standard order，也叫耶茨顺序，得名于 Frank Yates，他是 Fisher 的同事之一，对试验的设计与分析作出了许多重要贡献）。

表 6.1 图 6.1 中试验的方差分析

变异来源	平方和	自由度	均方	$F_{q}$	$P$ 值
A	208.33	1	208.33	53.15	0.0001
B	75.00	1	75.00	19.13	0.0024
AB	8.33	1	8.33	2.13	0.1826
误差	31.34	8	3.92
总计	323.00	11

使用这一标准顺序可以看到，估计效应时所用的对照系数为

效应	(1)	$a$	$b$	$ab$
A	−1	+1	−1	+1
B	−1	−1	+1	+1
AB	+1	−1	−1	+1

注意估计交互效应的对照系数恰好是相应两个主效应系数的乘积。对照系数总是 +1 或 −1，因此可以用表 6.2 这样的正负号表来确定每个处理组合的恰当符号。表 6.2 的列标题是主效应（A 和 B）、 $AB$ 交互作用，以及 I，它代表整个试验的总和或平均值。注意与 I 对应的那一列只有正号。行标记是各处理组合。要得到估计任一效应的对照，只需把表中相应列的符号与对应的处理组合相乘再相加。例如，为估计 A，对照为 $-(1) + a - b + ab$ ，这与式 6.1 一致。注意效应 A、B、AB 的对照是正交的。因此，2² 设计（以及所有 $2^{k}$ 设计）是正交设计。因子低水平与高水平的 ±1 编码通常称为正交编码（orthogonal coding）或效应编码（effects coding）。

表 6.2 2² 设计中计算效应的代数符号

处理组合	I	A	B	AB
(1)	+	−	−	+
$a$	+	+	−	−
$b$	+	−	+	−
$ab$	+	+	+	+

**回归模型。**在 $2^{k}$ 因子设计中，很容易用回归模型来表示试验结果。因为 $2^{k}$ 设计只是一个因子设计，我们也可以使用效应模型或均值模型，但回归模型的方法要自然和直观得多。对图 6.1 中的化学过程试验，回归模型为

y = \beta_{0} + \beta_{1} x_{1} + \beta_{2} x_{2} + \epsilon

其中 $x_{1}$ 是表示反应物浓度的编码变量， $x_{2}$ 是表示催化剂用量的编码变量，各 $\beta$ 是回归系数。自然变量（反应物浓度和催化剂用量）与编码变量之间的关系为

x_{1} = \frac{\operatorname{Conc} - \left(\operatorname{Conc}_{\text{low}} + \operatorname{Conc}_{\text{high}}\right) / 2}{\left(\operatorname{Conc}_{\text{high}} - \operatorname{Conc}_{\text{low}}\right) / 2}

以及

x_{2} = \frac{\text{Catalyst} - (\text{Catalyst}_{\text{low}} + \text{Catalyst}_{\text{high}}) / 2}{(\text{Catalyst}_{\text{high}} - \text{Catalyst}_{\text{low}}) / 2}

当自然变量只有两个水平时，这种编码对编码变量的水平给出我们熟悉的 $\pm 1$ 记号。为对本例加以说明，注意

\begin{array}{r l} x_{1} & = \frac{\text{Conc} - (15 + 25) / 2}{(25 - 15) / 2} \\ & = \frac{\text{Conc} - 20}{5} \end{array}

因此，如果浓度处于高水平（Conc = 25%），则 $x_{1} = +1$ ；如果浓度处于低水平（Conc = 15%），则 $x_{1} = -1$ 。此外，

\begin{array}{r l} x_{2} & = \frac{\text{Catalyst} - (1 + 2) / 2}{(2 - 1) / 2} \\ & = \frac{\text{Catalyst} - 1.5}{0.5} \end{array}

因此，如果催化剂处于高水平（Catalyst = 2 磅），则 $x_{2} = +1$ ；如果催化剂处于低水平（Catalyst = 1 磅），则 $x_{2} = -1$ 。

拟合的回归模型为

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) x_{1} + \left(\frac{-5.00}{2}\right) x_{2}

其中截距是全部 12 个观测值的总平均，回归系数 $\hat{\beta}_{1}$ 与 $\hat{\beta}_{2}$ 是对应因子效应估计值的一半。回归系数之所以是效应估计值的一半，是因为回归系数度量的是 $x$ 每变化一个单位对 $y$ 均值的影响，而效应估计是基于两个单位的变化（从 −1 到 +1）。这种估计回归系数的简便方法所得结果恰好是最小二乘参数估计。我们将在 6.7 节再次回到这个话题。另见本章的补充材料。

**需要多少次重复？**几乎每个试验都会遇到的一个标准问题是：需要多少次重复？我们在前面几章已讨论过这个问题，但其中有些方面在大量用于因子筛选的 $2^{k}$ 设计中特别有用。所谓因子筛选，就是研究一组 $k$ 个因子以确定哪些是活跃的。回顾前面的讨论可知，设计试验时样本量的恰当选择取决于感兴趣的效应有多大、统计检验的功效以及第一类错误的选择。一个重要效应的大小显然依赖于具体问题，但在许多实际情形中，试验者感兴趣的是检出至少两倍于误差标准差（ $2\sigma$ ）那么大的效应。更小的效应通常不太受关注，因为改变与这么小的效应相关联的因子，所引起的响应变化相对于系统中背景噪声往往非常小。足够的功效也依赖于具体问题，但在许多实际情形中，应以达到至少 0.80（即 80%）的功效为目标。

我们将用 2² 化学过程试验来说明如何确定样本量的恰当选择。假设我们感兴趣的是检出大小为 $2\sigma$ 的效应。如果基本的 2² 设计重复两次，共 8 次试验，则有 4 个自由度可用于估计与模型无关的误差估计（纯误差）。如果试验者使用 $\alpha = 0.05$ 的显著性水平或第一类错误率，该设计的功效为 0.572，即 57.2%。这太低了，试验者应考虑更多重复。另一种在筛选试验中可能有用的选择是采用更高的第一类错误率。在筛选试验中，第一类错误（认为某个因子活跃而实际上它不活跃）通常不像第二类错误（未能识别出一个活跃因子）那样具有同样的影响。如果一个因子被误认为活跃，这一错误会在后续工作中被发现，因此第一类错误的后果通常很小。然而，未能识别出一个活跃因子通常问题很大，因为该因子会被搁置一旁，往往再也不会被考虑。所以在筛选试验中，试验者常常愿意考虑更高的第一类错误率，比如说 0.10 或 0.20。

假设我们在化学过程试验中使用 $\alpha = 0.10$ 。这会使功效达到 75%。使用 $\alpha = 0.20$ 则把功效提高到 89%，这是一个非常合理的值。另一种选择是通过增加重复来增大样本量。如果我们使用三次重复，就有 8 个自由度用于纯误差；如果希望在 $\alpha = 0.05$ 下检出大小为 $2\sigma$ 的效应，该设计的功效为 85.7%。这是一个非常好的功效值，因此试验者决定采用 2² 设计的三次重复。

可以用软件包作出上述功效计算。下面的方框显示给出 JMP 的功效计算结果。该模型同时含两个主效应和二因子交互作用；大小为 $2\sigma$ 的效应是这样设定的：把误差均方的平方根（Anticipated RMSE）设为 1，并把每个预期模型系数的大小设为 1。

Evaluate Design

Model

Intercept

X1*X2

Power Analysis

Significance Level		0.05
Anticipated RMSE		1

Term	Anticipated Coefficient	Power
Intercept	1	0.857
X1	1	0.857
X2	1	0.857
X1*X2	1	0.857

**残差与模型适合性。**回归模型可用于得到设计中四个点处 $y$ 的预测值或拟合值。残差是 $y$ 的观测值与拟合值之差。例如，当反应物浓度处于低水平（ $x_{1} = -1$ ）且催化剂处于低水平（ $x_{2} = -1$ ）时，预测产率为

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) (- 1) + \left(\frac{-5.00}{2}\right) (- 1) = 25.835

在该处理组合处有三个观测值，其残差为

\begin{array}{l} e_{1} = 28 - 25.835 = 2.165 \\ e_{2} = 25 - 25.835 = - 0.835 \\ e_{3} = 27 - 25.835 = 1.165 \end{array}

其余预测值和残差用类似方法计算。对反应物浓度取高水平、催化剂取低水平，

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) (+ 1) + \left(\frac{-5.00}{2}\right) (- 1) = 34.165

以及

e_{4} = 36 - 34.165 = 1.835

e_{5} = 32 - 34.165 = - 2.165

e_{6} = 32 - 34.165 = - 2.165

对反应物浓度取低水平、催化剂取高水平，

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) (- 1) + \left(\frac{-5.00}{2}\right) (+ 1) = 20.835

以及

\begin{array}{l} e_{7} = 18 - 20.835 = - 2.835 \\ e_{8} = 19 - 20.835 = - 1.835 \\ e_{9} = 23 - 20.835 = 2.165 \end{array}

最后，对两个因子都取高水平，

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) (+ 1) + \left(\frac{-5.00}{2}\right) (+ 1) = 29.165

以及

\begin{array}{l} e_{10} = 31 - 29.165 = 1.835 \\ e_{11} = 30 - 29.165 = 0.835 \\ e_{12} = 29 - 29.165 = - 0.165 \end{array}

图 6.2 给出这些残差的正态概率图以及残差对预测产率的图。这些图形看起来是令人满意的，所以我们没有理由怀疑所得结论的有效性方面存在任何问题。

图 6.2 化学过程试验的残差图

**响应曲面。**回归模型

\hat{y} = 27.5 + \left(\frac{8.33}{2}\right) x_{1} + \left(\frac{-5.00}{2}\right) x_{2}

可用于生成响应曲面图。如果希望用自然因子水平来表示这些图，只需把前面给出的自然变量与编码变量之间的关系代入回归模型，得到

\begin{array}{l} \hat{y} = 27.5 + \left(\frac{8.33}{2}\right) \left(\frac{\text{Conc} - 20}{5}\right) + \left(\frac{-5.00}{2}\right) \left(\frac{\text{Catalyst} - 1.5}{0.5}\right) \\ = 18.33 + 0.8333 \text{Conc} - 5.00 \text{Catalyst} \end{array}

图 6.3 化学过程试验产率的响应曲面图与等高线图

图 6.3a 给出该模型产率的三维响应曲面图，图 6.3b 是等高线图。因为模型是一阶的（即只含主效应），拟合的响应曲面是一个平面。考察等高线图可以看到，随着反应物浓度增大、催化剂量减少，产率提高。我们常常用这样的拟合曲面来找出过程潜在改进的方向。一种正规的做法叫最速上升法（method of steepest ascent），将在第 11 章讨论系统探索响应曲面的方法时介绍。