4.4 平衡不完全区组设计

在某些使用随机区组设计的试验中，我们可能无法在每个区组中实施所有处理组合。出现这种情况通常是由于试验仪器或设施不足，或者区组的物理尺寸有限。例如，在血管移植物试验（例 4.1）中，假设每批材料只够测试三个挤压压力。因此每个压力无法在每批材料中测试。对这类问题，可以使用并非每个处理都出现在每个区组中的随机区组设计。这类设计称为随机不完全区组设计（randomized incomplete block design）。

当所有处理比较都同等重要时，每个区组中所用的处理组合应当以平衡的方式选取，使得任意一对处理共同出现的次数与其他任意一对相同。因此，平衡不完全区组设计（balanced incomplete block design，BIBD）是一种任意两个处理共同出现次数相同的不完全区组设计。假设有 $a$ 个处理，每个区组恰好能容纳 $k$（$k < a$）个处理。取 $\binom{a}{k}$ 个区组并为每个区组分配一种不同的处理组合，就可以构造一个平衡不完全区组设计。不过，常常用少于 $\binom{a}{k}$ 个区组就能获得平衡。Fisher and Yates (1953)、Davies (1956) 以及 Cochran and Cox (1957) 中给出了 BIBD 的表。

作为例子，假设一位化学工程师认为某个化学过程的反应时间是所用催化剂类型的函数。目前正在考察四种催化剂。试验过程是：选取一批原材料，装入中试装置，在装置上分别用每种催化剂各独立运行一次，并观测反应时间。由于原材料批次间的变异可能影响催化剂的表现，该工程师决定以原材料批次作为区组。然而每批材料只够运行三种催化剂。因此必须使用随机不完全区组设计。该试验的平衡不完全区组设计连同记录到的观测值见表 4.22。每个区组中催化剂的运行顺序是随机的。

4.4.1 BIBD 的统计分析¶

如通常一样，假定有 $a$ 个处理和 $b$ 个区组。此外假定每个区组含有 $k$ 个处理，每个处理在设计中出现 $r$ 次（或被重复 $r$ 次），观测值总数为 $N = ar = bk$。进一步地，每一对处理出现在同一区组中的次数为

若 $a = b$，则称该设计是对称的。

参数 $\lambda$ 必须是整数。为导出 $\lambda$ 的关系式，考虑任一个处理，比如说处理 1。由于处理 1 出现在 $r$ 个区组中，而这些区组中每个区组还有 $k-1$ 个其他处理，所以包含处理 1 的区组中共有 $r(k-1)$ 个观测值。这 $r(k-1)$ 个观测值还必须把其余 $a-1$ 个处理各表示 $\lambda$ 次。因此 $\lambda(a-1)=r(k-1)$。

BIBD 的统计模型为

其中 $y_{ij}$ 是第 $j$ 个区组中的第 $i$ 个观测值，$\mu$ 是总均值，$\tau_{i}$ 是第 $i$ 个处理的效应，$\beta_{j}$ 是第 $j$ 个区组的效应，$\varepsilon_{ij}$ 是 $\mathrm{NID}(0,\sigma^{2})$ 随机误差成分。数据中的总变异用总校正平方和表示：

表 4.22 催化剂试验的平衡不完全区组设计

总变异可以分解为

其中处理平方和经过调整，以区分处理效应与区组效应。这一调整是必要的，因为每种处理出现在不同的 $r$ 个区组中。因此未调整的处理总和 $y_{1.}, y_{2.}, \ldots, y_{a.}$ 之间的差异也受到区组之间差异的影响。

区组平方和为

其中 $y_{.j}$ 是第 $j$ 个区组中的总和。$SS_{\text{区组}}$ 有 $b-1$ 个自由度。经调整的处理平方和为

其中 $Q_{i}$ 是第 $i$ 个处理的经调整总和，按下式计算：

当处理 $i$ 出现在区组 $j$ 中时 $n_{ij}=1$，否则 $n_{ij}=0$。经调整的处理总和总是等于零。$SS_{\text{处理（经调整）}}$ 有 $a-1$ 个自由度。误差平方和由相减得到：

它有 $N - a - b + 1$ 个自由度。

检验各处理效应是否相等的适当统计量为

方差分析汇总于表 4.23。

表 4.23 平衡不完全区组设计的方差分析

例 4.4¶

考虑表 4.22 中催化剂试验的数据。这是一个 BIBD，其中 $a = 4$，$b = 4$，$k = 3$，$r = 3$，$\lambda = 2$，$N = 12$。该数据的分析如下。总平方和为

区组平方和由式 4.32 得到：

为计算剔除区组之后经调整的处理平方和，首先用式 4.34 确定经调整的处理总和：

经调整的处理平方和由式 4.33 计算：

误差平方和由相减得到：

方差分析见表 4.24。由于 $P$ 值很小，我们得出结论：所用催化剂对反应时间有显著影响。

表 4.24 例 4.4 的方差分析

若所研究的因子是固定的，则可能关心对各处理均值作检验。如果使用正交对照，对照必须建立在经调整的处理总和 $\{Q_{i}\}$ 上，而不是 $\{y_{i.}\}$ 上。对照平方和为

其中 $\{c_{i}\}$ 是对照系数。其他多重比较方法可以用来比较所有成对的经调整处理效应，我们将在 4.4.2 节看到这些效应由 $\hat{\tau}_{i}=kQ_{i}/(\lambda a)$ 估计。经调整处理效应的标准误为

在我们所描述的分析中，总平方和被分解为经调整的处理平方和、未调整的区组平方和以及误差平方和。有时我们想考察区组效应。为此需要对 $SS_{T}$ 作另一种分解，即

这里 $SS_{\text{处理}}$ 是未调整的。如果设计是对称的，即 $a = b$，则可以得到 $SS_{\text{区组（经调整）}}$ 的一个简单公式。经调整的区组总和为^[3]

以及

例 4.4 中的 BIBD 是对称的，因为 $a = b = 4$。因此

以及

此外，

对称 BIBD 的方差分析汇总见表 4.25。注意表 4.25 中各均方对应的平方和之和不等于总平方和，即

这是处理与区组不正交的必然结果。

表 4.25 例 4.4 的方差分析（同时含处理与区组）

**计算机输出。**有若干计算机软件包可以对平衡不完全区组设计进行分析。SAS 的一般线性模型（General Linear Models）过程是其中之一，Minitab 和 JMP 也是。表 4.26 的上半部分是例 4.4 的 Minitab 一般线性模型输出。比较表 4.26 和表 4.25 可以看出，Minitab 计算了经调整的处理平方和与经调整的区组平方和（在 Minitab 输出中称为 “Adj SS”）。^[1]

表 4.26 的下半部分是使用 Tukey 方法的多重比较分析。其中给出了所有均值对之差的置信区间以及 Tukey 检验。注意 Tukey 方法会引导我们得出催化剂 4 与其他三种不同的结论。

表 4.26 例 4.4 的 Minitab（一般线性模型）分析

General Linear Model
Factor Type Levels Values
Catalyst fixed 4 1 2 3 4
Block fixed 4 1 2 3 4

Analysis of Variance for Time, using Adjusted SS for Tests
Source DF Seq SS Adj SS Adj MS F P
Catalyst 3 11.667 22.750 7.583 11.67 0.011
Block 3 66.083 66.083 22.028 33.89 0.001
Error 5 3.250 3.250 0.650
Total 11 81.000

Tukey 95.0% Simultaneous Confidence Intervals
Response Variable Time
All Pairwise Comparisons among Levels of Catalyst

Catalyst = 1 subtracted from:
Catalyst Lower Center Upper ----+----+----+----
2 -2.327 0.2500 2.827 (----*----)
3 -1.952 0.6250 3.202 (----*----)
4 1.048 3.6250 6.202 (----*----)
0.0 2.5 5.0

Catalyst = 2 subtracted from:
Catalyst Lower Center Upper ----+----+----+----
3 -2.202 0.3750 2.952 (----*----)
4 0.798 3.3750 5.952 (----*----)
0.0 2.5 5.0

Catalyst = 3 subtracted from:
Catalyst Lower Center Upper ----+----+----+----
4 0.4228 3.000 5.577 (----*----)
0.0 2.5 5.0

Tukey Simultaneous Tests
Response Variable Time
All Pairwise Comparisons among Levels of Catalyst

Catalyst = 1 subtracted from:
Level Difference SE of Adjusted
Catalyst of Means Difference T-Value P-Value
2 0.2500 0.6982 0.3581 0.9825
3 0.6250 0.6982 0.8951 0.8085
4 3.6250 0.6982 5.1918 0.0130

Catalyst = 2 subtracted from:
Level Difference SE of Adjusted
Catalyst of Means Difference T-Value P-Value
3 0.3750 0.6982 0.5371 0.9462
4 3.3750 0.6982 4.8338 0.0175

Catalyst = 3 subtracted from:
Level Difference SE of Adjusted
Catalyst of Means Difference T-Value P-Value
4 3.000 0.6982 4.297 0.0281

4.4.2 参数的最小二乘估计¶

考虑估计 BIBD 模型中的处理效应。最小二乘正规方程为

施加约束 $\sum \hat{\tau}_{i} = \sum \hat{\beta}_{j} = 0$，可得 $\hat{\mu} = \overline{y}_{..}$。进一步地，用 $\{\beta_{j}\}$ 的方程从 $\{\tau_{i}\}$ 的方程中消去区组效应，得到

注意式 4.40 的右端就是 $kQ_{i}$，其中 $Q_{i}$ 是第 $i$ 个经调整处理总和（见式 4.34）。^[2] 现在，因为当 $p \neq i$ 时 $\sum_{j=1}^{b} n_{ij} n_{pj} = \lambda$，且 $n_{pj}^{2} = n_{pj}$（因为 $n_{pj} = 0$ 或 1），我们可以把式 4.40 改写为

最后，注意约束 $\sum_{i=1}^{a} \hat{\tau}_{i} = 0$ 意味着 $\sum_{\substack{p=1 \\ p \neq i}}^{a} \hat{\tau}_{p} = -\hat{\tau}_{i}$，并回想 $r(k-1) = \lambda(a-1)$，可得

因此，平衡不完全区组模型中处理效应的最小二乘估计为

作为例示，考虑例 4.4 中的 BIBD。因为 $Q_{1} = -9/3$，$Q_{2} = -7/3$，$Q_{3} = -4/3$，$Q_{4} = 20/3$，我们得到

与我们在 4.4.1 节得到的结果一致。

4.4.3 BIBD 中区组间信息的恢复¶

4.4.1 节给出的 BIBD 分析通常称为区组内分析（intrablock analysis），因为区组差异被消除，而处理效应的所有对照都可以表示为同一区组内观测值之间的比较。无论区组是固定还是随机，这一分析都是恰当的。Yates (1940) 指出，如果区组效应是均值为零、方差为 $\sigma_{\beta}^{2}$ 的不相关随机变量，则可以获得关于处理效应 $\tau_{i}$ 的额外信息。Yates 把获得这种额外信息的方法称为区组间分析（interblock analysis）。

把区组总和 $y_{.j}$ 看作 $b$ 个观测值的集合，这些观测值的模型［取自 John (1971)］为

其中括号内的项可以视为误差。$\mu$ 和 $\tau_{i}$ 的区组间估计通过最小化最小二乘函数

得到。这给出如下最小二乘正规方程：

其中 $\tilde{\mu}$ 和 $\tilde{\tau}_{i}$ 表示区组间估计。施加约束 $\sum_{i=1}^{a} \tilde{\tau}_{i} = 0$，得到式 4.45 的解为

可以证明区组间估计 $\{\tilde{\tau}_{i}\}$ 与区组内估计 $\{\hat{\tau}_{i}\}$ 不相关。

区组间估计 $\{\tilde{\tau}_{i}\}$ 可能与区组内估计 $\{\hat{\tau}_{i}\}$ 不同。例如，例 4.4 中 BIBD 的区组间估计计算如下：

注意 $\sum_{j=1}^{b} n_{ij} y_{.j}$ 的这些值前面（原书第 164 页）在计算区组内分析中的经调整处理总和时已经用过。

现在假设我们希望把区组间估计与区组内估计结合起来，得到每个 $\tau_{i}$ 的单个无偏最小方差估计。可以证明 $\hat{\tau}_{i}$ 和 $\tilde{\tau}_{i}$ 都是无偏的，并且^[4]

以及

我们用两个估计的线性组合，比如说

来估计 $\tau_{i}$。对这一估计方法，最小方差无偏组合估计 $\tau_{i}^{*}$ 的权重应为 $\alpha_{1}=u_{1}/(u_{1}+u_{2})$ 和 $\alpha_{2}=u_{2}/(u_{1}+u_{2})$，其中 $u_{1}=1/V(\hat{\tau}_{i})$，$u_{2}=1/V(\tilde{\tau}_{i})$。因此最优权重与 $\hat{\tau}_{i}$ 和 $\tilde{\tau}_{i}$ 的方差成反比。这意味着最佳组合估计为

它可以简化为

遗憾的是，式 4.49 不能用来估计 $\tau_{i}$，因为方差 $\sigma^{2}$ 和 $\sigma_{\beta}^{2}$ 是未知的。通常的做法是由数据估计 $\sigma^{2}$ 和 $\sigma_{\beta}^{2}$，并用这些估计代替式 4.49 中的参数。$\sigma^{2}$ 通常取区组内方差分析中的误差均方，即区组内误差。因此

$\sigma_{\beta}^{2}$ 的估计则由对处理调整后的区组均方得到。一般地，对平衡不完全区组设计，这一均方为

其期望值［由 Graybill (1961) 导出］为

因此，若 $MS_{\text{区组（经调整）}} > MS_{E}$，$\hat{\sigma}_{\beta}^{2}$ 的估计为

而若 $MS_{\text{区组（经调整）}} \leq MS_{E}$，则取 $\hat{\sigma}_{\beta}^{2} = 0$。这样就得到组合估计

现在计算例 4.4 数据的组合估计。由表 4.25 得到 $\hat{\sigma}^{2} = MS_{E} = 0.65$ 和 $MS_{\text{区组（经调整）}} = 22.03$。（注意在计算 $MS_{\text{区组（经调整）}}$ 时我们用到了该设计是对称的这一事实。）一般地必须使用式 4.50。因为 $MS_{\text{区组（经调整）}} > MS_{E}$，我们用式 4.51 估计 $\sigma_{\beta}^{2}$：

因此，我们可以把 $\hat{\sigma}^{2}=0.65$ 和 $\hat{\sigma}_{\beta}^{2}=8.02$ 代入式 4.52a，得到下面列出的组合估计。为方便起见，同时也给出区组内估计和区组间估计。在本例中，组合估计与区组内估计接近，因为区组间估计的方差相对较大。