5.6 因子设计中的区组化
前面我们在完全随机化试验的背景下讨论了因子设计。有时把因子设计中的所有运行完全随机化并不可行或不实际。例如，某个干扰因子的存在可能要求试验分若干区组进行。我们在第 4 章以单因子试验为背景讨论了区组化的基本概念。现在说明如何把区组化纳入因子设计。因子设计中区组化的其他一些方面将在第 7—9 章和第 13 章介绍。
考虑一个有 n n n
y i j k = μ + τ i + β j + ( τ β ) i j + ϵ i j k { i = 1 , 2 , … , a j = 1 , 2 , … , b k = 1 , 2 , … , n (5.33) y_{ijk} = \mu + \tau_{i} + \beta_{j} + (\tau\beta)_{ij} + \epsilon_{ijk} \qquad \left\{ \begin{array}{l} i = 1, 2, \ldots, a \\ j = 1, 2, \ldots, b \\ k = 1, 2, \ldots, n \end{array} \right. \tag{5.33} y ijk = μ + τ i + β j + ( τ β ) ij + ϵ ijk ⎩ ⎨ ⎧ i = 1 , 2 , … , a j = 1 , 2 , … , b k = 1 , 2 , … , n ( 5.33 ) 其中 τ i \tau_{i} τ i β j \beta_{j} β j ( τ β ) i j (\tau\beta)_{ij} ( τ β ) ij A B AB A B a b n abn abn a b ab ab n n n
y i j k = μ + τ i + β j + ( τ β ) i j + δ k + ϵ i j k { i = 1 , 2 , … , a j = 1 , 2 , … , b k = 1 , 2 , … , n (5.34) y_{ijk} = \mu + \tau_{i} + \beta_{j} + (\tau\beta)_{ij} + \delta_{k} + \epsilon_{ijk} \quad \left\{ \begin{array}{l} i = 1, 2, \dots, a \\ j = 1, 2, \dots, b \\ k = 1, 2, \dots, n \end{array} \right. \tag{5.34} y ijk = μ + τ i + β j + ( τ β ) ij + δ k + ϵ ijk ⎩ ⎨ ⎧ i = 1 , 2 , … , a j = 1 , 2 , … , b k = 1 , 2 , … , n ( 5.34 ) 其中 δ k \delta_{k} δ k k k k
模型（式 5.34）假定区组与处理之间的交互作用可以忽略。这一点在随机区组设计的分析中已经假定过。如果这些交互作用确实存在，它们无法与误差成分分离。事实上，该模型中的误差项实际上由 ( τ δ ) i k (\tau\delta)_{ik} ( τ δ ) ik ( β δ ) j k (\beta\delta)_{jk} ( β δ ) jk ( τ β δ ) i j k (\tau\beta\delta)_{ijk} ( τ β δ ) ijk n n n { y . . k } \{y_{..k}\} { y .. k } σ δ 2 \sigma_{\delta}^{2} σ δ 2
σ ^ δ 2 = M S 区组 − M S E a b \hat{\sigma}_{\delta}^{2} = \frac{MS_{\text{区组}} - MS_{E}}{ab} σ ^ δ 2 = ab M S 区组 − M S E 在前面的例子中，随机化仅限于一批原材料之内。在实践中，各种现象都可能造成随机化限制，例如时间和操作人员。例如，如果我们无法在一天之内完成整个因子试验，那么试验者可以在第 1 天完成一次完整的重复，第 2 天完成第二次重复，依此类推。于是每一天就是一个区组。
表 5.20 随机完全区组中两因子因子设计的方差分析
变异来源 平方和 自由度 期望均方 F 0 F_{0} F 0 区组 1 a b ∑ k y . . k 2 − y . . . 2 a b n \frac{1}{ab}\sum_{k} y_{..k}^{2} - \frac{y_{...}^{2}}{abn} ab 1 ∑ k y .. k 2 − abn y ... 2 n − 1 n-1 n − 1 σ 2 + a b σ δ 2 \sigma^{2} + ab\sigma_{\delta}^{2} σ 2 + ab σ δ 2 A 1 b n ∑ i y i . . 2 − y . . . 2 a b n \frac{1}{bn}\sum_{i} y_{i..}^{2} - \frac{y_{...}^{2}}{abn} bn 1 ∑ i y i .. 2 − abn y ... 2 a − 1 a-1 a − 1 σ 2 + b n ∑ τ i 2 a − 1 \sigma^{2} + \frac{bn \sum \tau_{i}^{2}}{a-1} σ 2 + a − 1 bn ∑ τ i 2 M S A M S E \frac{MS_{A}}{MS_{E}} M S E M S A B 1 a n ∑ j y . j . 2 − y . . . 2 a b n \frac{1}{an}\sum_{j} y_{.j.}^{2} - \frac{y_{...}^{2}}{abn} an 1 ∑ j y . j . 2 − abn y ... 2 b − 1 b-1 b − 1 σ 2 + a n ∑ β j 2 b − 1 \sigma^{2} + \frac{an \sum \beta_{j}^{2}}{b-1} σ 2 + b − 1 an ∑ β j 2 M S B M S E \frac{MS_{B}}{MS_{E}} M S E M S B AB 1 n ∑ i ∑ j y i j . 2 − y . . . 2 a b n − S S A − S S B \frac{1}{n}\sum_{i}\sum_{j} y_{ij.}^{2} - \frac{y_{...}^{2}}{abn} - SS_{A} - SS_{B} n 1 ∑ i ∑ j y ij . 2 − abn y ... 2 − S S A − S S B ( a − 1 ) ( b − 1 ) (a-1)(b-1) ( a − 1 ) ( b − 1 ) σ 2 + n ∑ ∑ ( τ β ) i j 2 ( a − 1 ) ( b − 1 ) \sigma^{2} + \frac{n \sum \sum (\tau\beta)_{ij}^{2}}{(a-1)(b-1)} σ 2 + ( a − 1 ) ( b − 1 ) n ∑∑ ( τ β ) ij 2 M S A B M S E \frac{MS_{AB}}{MS_{E}} M S E M S A B 误差 相减 ( a b − 1 ) ( n − 1 ) (ab-1)(n-1) ( ab − 1 ) ( n − 1 ) σ 2 \sigma^{2} σ 2 总计 ∑ i ∑ j ∑ k y i j k 2 − y . . . 2 a b n \sum_{i}\sum_{j}\sum_{k} y_{ijk}^{2} - \frac{y_{...}^{2}}{abn} ∑ i ∑ j ∑ k y ijk 2 − abn y ... 2 a b n − 1 abn-1 abn − 1
例 5.6 ¶ 一位工程师正在研究提高在雷达屏上探测目标能力的各种方法。她认为两个重要因子是屏上背景噪声（即"地面杂波"）的强度和置于屏幕上的滤波器类型。试验设计使用三个水平的地面杂波和两种滤波器类型。我们把它们当作固定型因子。试验的做法是随机选取一个处理组合（地面杂波水平与滤波器类型），然后把代表目标的信号引入屏上。该目标的强度不断增加，直到操作人员观察到它为止。检测时的强度水平作为响应变量来测量。由于操作人员的可用性，方便的做法是选定一名操作人员，让他或她一直守在屏前，直到所有必要的运行都完成为止。此外，操作人员使用屏幕的技能和能力各不相同。因此，把操作人员作为区组看来是合理的。随机选取四名操作人员。一旦选定某名操作人员，六个处理组合的实施顺序就随机确定。这样我们就有一个在随机完全区组中实施的 3 × 2 3 \times 2 3 × 2
该试验的线性模型为
y i j k = μ + τ i + β j + ( τ β ) i j + δ k + ϵ i j k { i = 1 , 2 , 3 j = 1 , 2 k = 1 , 2 , 3 , 4 y_{ijk} = \mu + \tau_{i} + \beta_{j} + (\tau\beta)_{ij} + \delta_{k} + \epsilon_{ijk} \qquad \left\{ \begin{array}{l} i = 1, 2, 3 \\ j = 1, 2 \\ k = 1, 2, 3, 4 \end{array} \right. y ijk = μ + τ i + β j + ( τ β ) ij + δ k + ϵ ijk ⎩ ⎨ ⎧ i = 1 , 2 , 3 j = 1 , 2 k = 1 , 2 , 3 , 4 其中 τ i \tau_{i} τ i β j \beta_{j} β j ( τ β ) i j (\tau\beta)_{ij} ( τ β ) ij δ k \delta_{k} δ k ϵ i j k \epsilon_{ijk} ϵ ijk N I D ( 0 , σ 2 ) \mathrm{NID}(0,\sigma^{2}) NID ( 0 , σ 2 ) { y . . k } \{y_{..k}\} { y .. k }
S S 区组 = 1 a b ∑ k = 1 n y . . k 2 − y . . . 2 a b n = 1 ( 3 ) ( 2 ) [ ( 572 ) 2 + ( 579 ) 2 + ( 597 ) 2 + ( 530 ) 2 ] − ( 2278 ) 2 ( 3 ) ( 2 ) ( 4 ) = 402.17 \begin{array}{r l} SS_{\text{区组}} & = \frac{1}{ab} \sum_{k=1}^{n} y_{..k}^{2} - \frac{y_{...}^{2}}{abn} \\ & = \frac{1}{(3)(2)} \left[ (572)^{2} + (579)^{2} + (597)^{2} + (530)^{2} \right] - \frac{(2278)^{2}}{(3)(2)(4)} \\ & = 402.17 \end{array} S S 区组 = ab 1 ∑ k = 1 n y .. k 2 − abn y ... 2 = ( 3 ) ( 2 ) 1 [ ( 572 ) 2 + ( 579 ) 2 + ( 597 ) 2 + ( 530 ) 2 ] − ( 3 ) ( 2 ) ( 4 ) ( 2278 ) 2 = 402.17 表 5.21 目标探测时的强度水平
滤波器类型 操作人员 1 操作人员 1 操作人员 2 操作人员 2 操作人员 3 操作人员 3 操作人员 4 操作人员 4 1 2 1 2 1 2 1 2 地面杂波：低 90 86 96 84 100 92 92 81 地面杂波：中 102 87 106 90 105 97 96 80 地面杂波：高 114 93 112 91 108 95 98 83
表 5.22 例 5.6 的方差分析
变异来源 平方和 自由度 均方 F 0 F_{0} F 0 P P P 地面杂波 (G) 335.58 2 167.79 15.13 0.0003 滤波器类型 (F) 1066.67 1 1066.67 96.19 <0.0001 GF 77.08 2 38.54 3.48 0.0573 区组 402.17 3 134.06 误差 166.33 15 11.09 总计 2047.83 23
该试验的完整方差分析汇总于表 5.22。表 5.22 的呈现方式表明，所有效应都通过把其均方除以误差均方来检验。地面杂波水平和滤波器类型在 1% 水平上显著，而它们的交互作用只在 10% 水平上显著。因此我们得出结论：地面杂波水平和所用屏幕滤波器类型都影响操作人员探测目标的能力，并且有证据表明这两个因子之间存在轻微的交互作用。区组方差分量的方差分析估计为
σ ^ δ 2 = M S 区组 − M S E a b = 134.06 − 11.09 ( 3 ) ( 2 ) = 20.50 \hat{\sigma}_{\delta}^{2} = \frac{MS_{\text{区组}} - MS_{E}}{ab} = \frac{134.06 - 11.09}{(3)(2)} = 20.50 σ ^ δ 2 = ab M S 区组 − M S E = ( 3 ) ( 2 ) 134.06 − 11.09 = 20.50 该试验的 JMP 输出见表 5.23。输出中给出了区组方差分量的残差最大似然（REML）估计，由于这是一个平衡设计，REML 估计与方差分析估计一致。JMP 还给出了两个方差分量 σ 2 \sigma^{2} σ 2 σ δ 2 \sigma_{\delta}^{2} σ δ 2
表 5.23 例 5.6 的 JMP 输出
Whole Model
Actual by Predicted Plot
Mean of Response 94.91667
Summary of Fit
Root Mean Square Error 3.329998
Observations (or Sum Wgts) 24
RSquare 0.917432
RSquare Adj 0.894497
表 5.23（续）
REML Variance Component Estimates
Random Effect Var Ratio Var Component Std Error 95% Lower 95% Upper Pct of Total Operators (Blocks) 1.8481964 20.494444 18.255128 -15.28495 56.273839 64.890 Residual 11.088889 4.0490897 6.0510389 26.561749 35.110 Total 31.583333 100.000
-2 LogLikelihood = 118.73680261
Covariance Matrix of Variance Component Estimates
Random Effect Operators (Blocks) Residual Operators (Blocks) 333.24972 -2.732521 Residual -2.732521 16.395128
Fixed Effect Tests
Source Nparm DF DFDen F Ratio Prob > F Clutter 2 2 15 15.1315 0.0003* Filter Type 1 1 15 96.1924 <.0001* Clutter*Filter Type 2 2 15 3.4757 0.0575
Residual by Predicted Plot
现在假设操作人员被视为区组。再假设由于所需的准备时间，每天只能做六次运行。于是天数成为第二种随机化限制，从而得到表 5.24 所示的 6 × 6 6 \times 6 6 × 6 [1] f i f_{i} f i g j g_{j} g j i i i j j j f 1 g 2 f_{1}g_{2} f 1 g 2 3 × 2 3 \times 2 3 × 2 3 × 2 = 6 3 \times 2 = 6 3 × 2 = 6 A = f 1 g 1 A = f_{1}g_{1} A = f 1 g 1 B = f 1 g 2 B = f_{1}g_{2} B = f 1 g 2 C = f 1 g 3 C = f_{1}g_{3} C = f 1 g 3 D = f 2 g 1 D = f_{2}g_{1} D = f 2 g 1 E = f 2 g 2 E = f_{2}g_{2} E = f 2 g 2 F = f 2 g 3 F = f_{2}g_{3} F = f 2 g 3
六个拉丁字母之间的五个自由度对应于滤波器类型的主效应（一个自由度）、地面杂波的主效应（两个自由度）以及它们的交互作用（两个自由度）。该设计的线性统计模型为
y i j k l = μ + α i + τ j + β k + ( τ β ) j k + θ l + ϵ i j k l { i = 1 , 2 , … , 6 j = 1 , 2 , 3 k = 1 , 2 l = 1 , 2 , … , 6 (5.35) y_{ijkl} = \mu + \alpha_{i} + \tau_{j} + \beta_{k} + (\tau\beta)_{jk} + \theta_{l} + \epsilon_{ijkl} \qquad \left\{ \begin{array}{l} i = 1, 2, \ldots, 6 \\ j = 1, 2, 3 \\ k = 1, 2 \\ l = 1, 2, \ldots, 6 \end{array} \right. \tag{5.35} y ijk l = μ + α i + τ j + β k + ( τ β ) jk + θ l + ϵ ijk l ⎩ ⎨ ⎧ i = 1 , 2 , … , 6 j = 1 , 2 , 3 k = 1 , 2 l = 1 , 2 , … , 6 ( 5.35 ) 其中 τ j \tau_{j} τ j β k \beta_{k} β k α i \alpha_{i} α i θ l \theta_{l} θ l
地面杂波 滤波器类型 1 滤波器类型 2 y j . . y_{j..} y j .. 低 560 512 1072 中 607 528 1135 高 646 543 1189 y . . k . y_{..k.} y .. k . 1813 1583 3396 = y . . . . 3396 = y_{....} 3396 = y ....
此外，行总和与列总和为
行（y j k l y_{jkl} y jk l 563 568 568 568 565 564 列（y i j k y_{ijk} y ijk 572 579 597 530 561 557
表 5.24 在 6 × 6 6 \times 6 6 × 6
天数 操作人员 1 操作人员 2 操作人员 3 操作人员 4 操作人员 5 操作人员 6 1 A ( f 1 g 1 = 90 ) A(f_{1}g_{1}=90) A ( f 1 g 1 = 90 ) B ( f 1 g 2 = 106 ) B(f_{1}g_{2}=106) B ( f 1 g 2 = 106 ) C ( f 1 g 3 = 108 ) C(f_{1}g_{3}=108) C ( f 1 g 3 = 108 ) D ( f 2 g 1 = 81 ) D(f_{2}g_{1}=81) D ( f 2 g 1 = 81 ) F ( f 2 g 3 = 90 ) F(f_{2}g_{3}=90) F ( f 2 g 3 = 90 ) E ( f 2 g 2 = 88 ) E(f_{2}g_{2}=88) E ( f 2 g 2 = 88 ) 2 C ( f 1 g 3 = 114 ) C(f_{1}g_{3}=114) C ( f 1 g 3 = 114 ) A ( f 1 g 1 = 96 ) A(f_{1}g_{1}=96) A ( f 1 g 1 = 96 ) B ( f 1 g 2 = 105 ) B(f_{1}g_{2}=105) B ( f 1 g 2 = 105 ) F ( f 2 g 3 = 83 ) F(f_{2}g_{3}=83) F ( f 2 g 3 = 83 ) E ( f 2 g 2 = 86 ) E(f_{2}g_{2}=86) E ( f 2 g 2 = 86 ) D ( f 2 g 1 = 84 ) D(f_{2}g_{1}=84) D ( f 2 g 1 = 84 ) 3 B ( f 1 g 2 = 102 ) B(f_{1}g_{2}=102) B ( f 1 g 2 = 102 ) E ( f 2 g 2 = 90 ) E(f_{2}g_{2}=90) E ( f 2 g 2 = 90 ) F ( f 2 g 3 = 95 ) F(f_{2}g_{3}=95) F ( f 2 g 3 = 95 ) A ( f 1 g 1 = 92 ) A(f_{1}g_{1}=92) A ( f 1 g 1 = 92 ) D ( f 2 g 1 = 85 ) D(f_{2}g_{1}=85) D ( f 2 g 1 = 85 ) C ( f 1 g 3 = 104 ) C(f_{1}g_{3}=104) C ( f 1 g 3 = 104 ) 4 E ( f 2 g 2 = 87 ) E(f_{2}g_{2}=87) E ( f 2 g 2 = 87 ) D ( f 2 g 1 = 84 ) D(f_{2}g_{1}=84) D ( f 2 g 1 = 84 ) A ( f 1 g 1 = 100 ) A(f_{1}g_{1}=100) A ( f 1 g 1 = 100 ) B ( f 1 g 2 = 96 ) B(f_{1}g_{2}=96) B ( f 1 g 2 = 96 ) C ( f 1 g 3 = 110 ) C(f_{1}g_{3}=110) C ( f 1 g 3 = 110 ) F ( f 2 g 3 = 91 ) F(f_{2}g_{3}=91) F ( f 2 g 3 = 91 ) 5 F ( f 2 g 3 = 93 ) F(f_{2}g_{3}=93) F ( f 2 g 3 = 93 ) C ( f 1 g 3 = 112 ) C(f_{1}g_{3}=112) C ( f 1 g 3 = 112 ) D ( f 2 g 1 = 92 ) D(f_{2}g_{1}=92) D ( f 2 g 1 = 92 ) E ( f 2 g 2 = 80 ) E(f_{2}g_{2}=80) E ( f 2 g 2 = 80 ) A ( f 1 g 1 = 90 ) A(f_{1}g_{1}=90) A ( f 1 g 1 = 90 ) B ( f 1 g 2 = 98 ) B(f_{1}g_{2}=98) B ( f 1 g 2 = 98 ) 6 D ( f 2 g 1 = 86 ) D(f_{2}g_{1}=86) D ( f 2 g 1 = 86 ) F ( f 2 g 3 = 91 ) F(f_{2}g_{3}=91) F ( f 2 g 3 = 91 ) E ( f 2 g 2 = 97 ) E(f_{2}g_{2}=97) E ( f 2 g 2 = 97 ) C ( f 1 g 3 = 98 ) C(f_{1}g_{3}=98) C ( f 1 g 3 = 98 ) B ( f 1 g 2 = 100 ) B(f_{1}g_{2}=100) B ( f 1 g 2 = 100 ) A ( f 1 g 1 = 92 ) A(f_{1}g_{1}=92) A ( f 1 g 1 = 92 )
表 5.25 把雷达探测试验作为拉丁方中的 3 × 2 3 \times 2 3 × 2
变异来源 平方和 自由度 自由度的一般公式 均方 F 0 F_{0} F 0 P P P 地面杂波 (G) 571.50 2 a − 1 a-1 a − 1 285.75 28.86 <0.0001 滤波器类型 (F) 1469.44 1 b − 1 b-1 b − 1 1469.44 148.43 <0.0001 GF 126.73 2 ( a − 1 ) ( b − 1 ) (a-1)(b-1) ( a − 1 ) ( b − 1 ) 63.37 6.40 0.0071 天数（行） 4.33 5 a b − 1 ab-1 ab − 1 0.87 操作人员（列） 428.00 5 a b − 1 ab-1 ab − 1 85.60 误差 198.00 20 ( a b − 1 ) ( a b − 2 ) (ab-1)(ab-2) ( ab − 1 ) ( ab − 2 ) 9.90 总计 2798.00 35 ( a b ) 2 − 1 (ab)^{2}-1 ( ab ) 2 − 1
方差分析汇总于表 5.25。我们在表中加了一列，说明每个平方和的自由度数目是如何确定的。[2]
