Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

3.9 随机效应模型

3.9.1 单个随机因子

试验者常常关注某个有大量可能水平的因子。如果试验者从该因子的水平总体中随机选取其中 aa 个水平，我们就说这个因子是随机的。由于试验中实际使用的因子水平是随机选定的，因此所作的推断是针对因子水平的整个总体的。我们假定因子水平总体要么是无限大的，要么大到可以看作无限大。因子水平总体小到需要采用有限总体方法的情形并不常见。关于有限总体情形，参见 Bennett 和 Franklin（1954）以及 Searle 和 Fawcett（1970）。

线性统计模型为

yij=μ+τi+ϵij{i=1,2,…,aj=1,2,…,n(3.44)y_{ij} = \mu+\tau_{i}+\epsilon_{ij} \quad \left\{ \begin{array}{l} i = 1, 2, \ldots, a \\ j = 1, 2, \ldots, n \end{array} \right. \tag{3.44}

其中处理效应 τi\tau_{i} 和 ϵij\epsilon_{ij} 都是随机变量。我们假定处理效应 τi\tau_{i} 是 NID(0,στ2)\mathrm{NID}(0,\sigma_{\tau}^{2}) 随机变量[1]，误差是 NID(0,σ2)\mathrm{NID}(0,\sigma^{2}) 随机变量，且 τi\tau_{i} 与 ϵij\epsilon_{ij} 独立。由于 τi\tau_{i} 与 ϵij\epsilon_{ij} 独立，任一观测值的方差为

V(yij)=στ2+σ2V(y_{ij}) = \sigma_{\tau}^{2}+\sigma^{2}

方差 στ2\sigma_{\tau}^{2} 和 σ2\sigma^{2} 称为方差分量（variance component），模型（式 3.44）称为方差分量模型或随机效应模型。随机效应模型中的观测值服从正态分布，因为它们是两个正态独立随机变量 τi\tau_{i} 与 ϵij\epsilon_{ij} 的线性组合。然而，与所有观测值 yijy_{ij} 都相互独立的固定效应情形不同，在随机模型中观测值 yijy_{ij} 只有在来自不同因子水平时才相互独立。具体地说，可以证明任意两个观测值的协方差为

Cov⁡(yij,yij′)=στ2j≠j′Cov⁡(yij,yi′j′)=0i≠i′\begin{array}{ll} \operatorname{Cov}(y_{ij},y_{ij'}) = \sigma_{\tau}^{2} & j \neq j' \\ \operatorname{Cov}(y_{ij},y_{i'j'}) = 0 & i \neq i' \end{array}

注意同一因子水平内的观测值都具有相同的协方差，因为在实施试验之前，我们预期该因子水平下的观测值彼此相似，因为它们具有相同的随机分量。一旦试验实施完毕，我们就可以假定所有观测值都相互独立，因为参数 τi\tau_{i} 已经确定，该处理中的观测值只因随机误差而不同。

单因子随机效应模型中观测值的协方差结构可以用观测值的协方差矩阵来表示。为说明这一点，假设我们有 a=3a=3 个处理、n=2n=2 次重复，共有 N=6N=6 个观测值，可以写成一个向量

y=[y11y12y21y22y31y32]\mathbf{y} = \left[ \begin{array}{c} y_{11} \\ y_{12} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32} \end{array} \right]

这些观测值的 6×66 \times 6 协方差矩阵为

Cov⁡(y)=[στ2+σ2στ20000στ2στ2+σ2000000στ2+σ2στ20000στ2στ2+σ2000000στ2+σ2στ20000στ2στ2+σ2]\operatorname{Cov}(\mathbf{y}) = \left[ \begin{array}{cccccc} \sigma_{\tau}^{2}+\sigma^{2} & \sigma_{\tau}^{2} & 0 & 0 & 0 & 0 \\ \sigma_{\tau}^{2} & \sigma_{\tau}^{2}+\sigma^{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & \sigma_{\tau}^{2}+\sigma^{2} & \sigma_{\tau}^{2} & 0 & 0 \\ 0 & 0 & \sigma_{\tau}^{2} & \sigma_{\tau}^{2}+\sigma^{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \sigma_{\tau}^{2}+\sigma^{2} & \sigma_{\tau}^{2} \\ 0 & 0 & 0 & 0 & \sigma_{\tau}^{2} & \sigma_{\tau}^{2}+\sigma^{2} \end{array} \right]

该矩阵的主对角线元素是每个观测值的方差，每个非对角线元素是一对观测值的协方差。

3.9.2 随机模型的方差分析

基本的方差分析平方和恒等式

SST=SS处理+SSE(3.45)SS_{T} = SS_{\text{处理}}+SS_{E} \tag{3.45}

仍然成立。也就是说，我们把观测值的总变异分解为一个度量处理之间变异的成分（SS处理SS_{\text{处理}}）和一个度量处理内部变异的成分（SSESS_{E}）。检验关于各个处理效应的假设意义不大，因为它们是随机选定的；我们更关心处理的总体，因此我们检验关于方差分量 στ2\sigma_{\tau}^{2} 的假设：

H0:στ2=0H1:στ2>0(3.46)\begin{array}{l} H_{0}: \sigma_{\tau}^{2} = 0 \\ H_{1}: \sigma_{\tau}^{2} > 0 \end{array} \tag{3.46}

若 στ2=0\sigma_{\tau}^{2}=0，则所有处理都相同；而若 στ2>0\sigma_{\tau}^{2}>0，则处理之间存在变异[2]。和前面一样，SSE/σ2SS_{E}/\sigma^{2} 服从自由度为 N−aN-a 的卡方分布；并且在原假设下，SS处理/σ2SS_{\text{处理}}/\sigma^{2} 服从自由度为 a−1a-1 的卡方分布。这两个随机变量相互独立。因此，在原假设 στ2=0\sigma_{\tau}^{2}=0 下，比值

F0=SS处理a−1SSEN−a=MS处理MSE(3.47)F_{0} = \frac{\frac{SS_{\text{处理}}}{a-1}}{\frac{SS_{E}}{N-a}} = \frac{MS_{\text{处理}}}{MS_{E}} \tag{3.47}

服从分子自由度为 a−1a-1、分母自由度为 N−aN-a 的 FF 分布。不过，我们需要考察期望均方才能完整地描述检验程序。

考虑

E(MS处理)=1a−1E(SS处理)=1a−1E[∑i=1ayi.2n−y..2N]=1a−1E[1n∑i=1a(∑j=1nμ+τi+ϵij)2−1N(∑i=1a∑j=1nμ+τi+ϵij)2]\begin{aligned} E(MS_{\text{处理}}) &= \frac{1}{a-1}E(SS_{\text{处理}}) = \frac{1}{a-1}E\left[\sum_{i=1}^{a}\frac{y_{i.}^{2}}{n}-\frac{y_{..}^{2}}{N}\right] \\ &= \frac{1}{a-1}E\left[\frac{1}{n}\sum_{i=1}^{a}\left(\sum_{j=1}^{n}\mu+\tau_{i}+\epsilon_{ij}\right)^{2}-\frac{1}{N}\left(\sum_{i=1}^{a}\sum_{j=1}^{n}\mu+\tau_{i}+\epsilon_{ij}\right)^{2}\right] \end{aligned}

对括号内的量作平方并取期望时，我们看到涉及 τi2\tau_{i}^{2} 的项被 στ2\sigma_{\tau}^{2} 取代，因为 E(τi)=0E(\tau_{i})=0。此外，涉及 ϵi.2\epsilon_{i.}^{2}、ϵ..2\epsilon_{..}^{2} 和 ∑i=1a∑j=1nτi2\sum_{i=1}^{a}\sum_{j=1}^{n}\tau_{i}^{2} 的项分别被 nσ2n\sigma^{2}、anσ2an\sigma^{2} 和 an2an^{2} 取代[3]。再者，所有涉及 τi\tau_{i} 与 ϵij\epsilon_{ij} 的交叉乘积项的期望都为零。由此得到

E(MS处理)=1a−1[Nμ2+Nστ2+aσ2−Nμ2−nστ2−σ2]E(MS_{\text{处理}}) = \frac{1}{a-1}[N\mu^{2}+N\sigma_{\tau}^{2}+a\sigma^{2}-N\mu^{2}-n\sigma_{\tau}^{2}-\sigma^{2}]

即

E(MS处理)=σ2+nστ2(3.48)E(MS_{\text{处理}}) = \sigma^{2}+n\sigma_{\tau}^{2} \tag{3.48}

类似地，可以证明

E(MSE)=σ2(3.49)E(MS_{E}) = \sigma^{2} \tag{3.49}

由期望均方可见，在 H0H_{0} 下检验统计量（式 3.47）的分子和分母都是 σ2\sigma^{2} 的无偏估计量；而在 H1H_{1} 下，分子的期望值大于分母的期望值。因此，我们应当在 F0F_{0} 取值过大时拒绝 H0H_{0}。这意味着使用上尾的单侧临界域，即若 F0>Fα,a−1,N−aF_{0} > F_{\alpha,a-1,N-a} 则拒绝 H0H_{0}。

随机效应模型的计算程序和方差分析与固定效应情形完全相同。然而结论差别很大，因为它们适用于处理的整个总体。

3.9.3 模型参数的估计

我们通常关注模型中方差分量（σ2\sigma^{2} 和 στ2\sigma_{\tau}^{2}）的估计。可以用来估计 σ2\sigma^{2} 和 στ2\sigma_{\tau}^{2} 的一个非常简单的程序称为方差分析法（analysis of variance method），因为它利用方差分析表中的各行。该程序把期望均方与方差分析表中的观测值相等同，并解出方差分量。在单因子随机效应模型中把观测均方与期望均方相等，我们得到

MS处理=σ2+nστ2MS_{\text{处理}} = \sigma^{2}+n\sigma_{\tau}^{2}

和

MSE=σ2MS_{E} = \sigma^{2}

因此，方差分量的估计量为

σ^2=MSE(3.50)\hat{\sigma}^{2} = MS_{E} \tag{3.50}

和

σ^τ2=MS处理−MSEn(3.51)\hat{\sigma}_{\tau}^{2} = \frac{MS_{\text{处理}}-MS_{E}}{n} \tag{3.51}

对样本量不等的情形，把式 3.51 中的 nn 换成

n0=1a−1[∑i=1ani−∑i=1ani2∑i=1ani](3.52)n_{0} = \frac{1}{a-1}\left[\sum_{i=1}^{a} n_{i}-\frac{\sum_{i=1}^{a} n_{i}^{2}}{\sum_{i=1}^{a} n_{i}}\right] \tag{3.52}

方差分量估计的方差分析法是一种矩法（method of moments）程序。它不需要正态性假定，并且给出的 σ2\sigma^{2} 和 στ2\sigma_{\tau}^{2} 的估计量是最佳二次无偏的（也就是说，在观测值的所有无偏二次函数中，这些估计量的方差最小）。还有一种基于最大似然的方法可以用来估计方差分量，稍后将会介绍。

有时方差分析法会产生方差分量的负估计。显然，按定义方差分量是非负的，因此方差分量的负估计令人有些担忧。一种做法是接受该估计，并把它当作方差分量真值为零的证据，假定负估计是抽样变异造成的。这有直观吸引力，但在理论上有些困难。例如，用零代替负估计会干扰其他估计的统计性质。另一种做法是用一种总能给出非负估计的方法重新估计该负的方差分量。还有一种做法是把负估计视为所假定的线性模型不正确的证据，并重新审视问题。关于方差分量估计的全面论述，见 Searle（1971a，1971b）、Searle、Casella 和 McCullogh（1992）以及 Burdick 和 Graybill（1992）。

例 3.10

一家纺织公司在许多台织机上织造某种织物。它希望各织机是齐同的，以便获得强度一致的织物。工艺工程师怀疑，除了同一织机上织物样本内通常存在的强度变异之外，织机之间可能还存在显著的强度变异。为考察这一点，她随机选取四台织机，并在每台织机生产的织物上作四次强度测定。该试验按随机顺序进行，所得数据见表 3.18。所作方差分析见表 3.19。由方差分析我们得出结论：厂内的织机之间存在显著差异。

方差分量估计为 σ^2=1.90\hat{\sigma}^{2}=1.90 和

σ^τ2=29.73−1.904=6.96\hat{\sigma}_{\tau}^{2} = \frac{29.73-1.90}{4} = 6.96

因此，任一强度观测值的方差估计为

σ^y=σ^2+σ^τ2=1.90+6.96=8.86\hat{\sigma}_{y} = \hat{\sigma}^{2}+\hat{\sigma}_{\tau}^{2} = 1.90+6.96 = 8.86

这一变异的大部分可归因于织机之间的差异。

表 3.18 例 3.10 的强度数据

织机观测值 1234yi.y_{i.}
198979996390
291909392366
396959795383
495969998388
y..=1527y_{..} = 1527

表 3.19 强度数据的方差分析

变异来源平方和自由度均方F0F_{0}PP 值
织机89.19329.7315.68<0.001
误差22.75121.90
总计111.9415

图 3.19 纤维强度问题中的过程输出

这个例子说明了方差分量的一个重要用途——分离影响产品或系统的不同变异来源。产品变异问题在质量保证中经常出现，而且往往难以分离出变异来源。例如，这项研究可能源于一个观察：织物强度的变异太大，如图 3.19a 所示。该图把过程输出（纤维强度）表示为方差 σ^y2=8.86\hat{\sigma}_{y}^{2}=8.86 的正态分布。（这是例 3.10 中任一强度观测值方差的估计。）图 3.19a 还显示了强度的上下规格限，可以比较容易地看出有相当大一部分过程输出落在规格限之外（图 3.19a 中带阴影的尾部区域）。工艺工程师问：为什么这么多织物不合格，必须报废、返工或降级为低质量产品？答案是：产品强度变异的大部分来自织机之间的差异。织机性能不同可能是由于安装不当、维护不良、监督不力、操作人员培训不足、输入纤维有缺陷等原因造成的。

工艺工程师现在必须设法分离出织机性能差异的具体原因。如果她能找出并消除这些织机间变异来源，过程输出的方差就可以大为降低，或许能低至 σ^y2=1.90\hat{\sigma}_{y}^{2}=1.90，即例 3.10 中织机内部（误差）方差分量的估计。图 3.19b 给出了 σ^y2=1.90\hat{\sigma}_{y}^{2}=1.90 时纤维强度的正态分布。注意输出中不合格产品的比例已被大幅降低。尽管不可能消除所有织机间变异，但显然这一方差分量的显著降低会大大提高所生产纤维的质量。

我们很容易求出方差分量 σ2\sigma^{2} 的置信区间。如果观测值是正态独立分布的，那么 (N−a)MSE/σ2(N-a)MS_{E}/\sigma^{2} 服从 χN−a2\chi_{N-a}^{2} 分布。因此

P[χ1−(α/2),N−a2≤(N−a)MSEσ2≤χα/2,N−a2]=1−αP\left[\chi_{1-(\alpha/2),N-a}^{2} \leq \frac{(N-a)MS_{E}}{\sigma^{2}} \leq \chi_{\alpha/2,N-a}^{2}\right] = 1-\alpha

σ2\sigma^{2} 的 100(1−α)100(1-\alpha) 百分置信区间为

(N−a)MSEχα/2,N−a2≤σ2≤(N−a)MSEχ1−(α/2),N−a2(3.53)\frac{(N-a)MS_{E}}{\chi_{\alpha/2,N-a}^{2}} \leq \sigma^{2} \leq \frac{(N-a)MS_{E}}{\chi_{1-(\alpha/2),N-a}^{2}} \tag{3.53}

由于 MSE=1.90MS_{E}=1.90、N=16N=16、a=4a=4、χ0.025,122=23.3367\chi_{0.025,12}^{2}=23.3367、χ0.975,122=4.4038\chi_{0.975,12}^{2}=4.4038，σ2\sigma^{2} 的 95% 置信区间为 0.9770≤σ2≤5.17750.9770 \leq \sigma^{2} \leq 5.1775[4]。

现在考虑方差分量 στ2\sigma_{\tau}^{2}。它的点估计量为

σ^τ2=MS处理−MSEn\hat{\sigma}_{\tau}^{2} = \frac{MS_{\text{处理}}-MS_{E}}{n}

随机变量 (a−1)MS处理/(σ2+nστ2)(a-1)MS_{\text{处理}}/(\sigma^{2}+n\sigma_{\tau}^{2}) 服从 χa−12\chi_{a-1}^{2} 分布，而 (N−a)MSE/σ2(N-a)MS_{E}/\sigma^{2} 服从 χN−a2\chi_{N-a}^{2} 分布。因此 σ^τ2\hat{\sigma}_{\tau}^{2} 的概率分布是两个卡方随机变量的线性组合，即

u1χa−12−u2χN−a2u_{1}\chi_{a-1}^{2}-u_{2}\chi_{N-a}^{2}

其中

u1=σ2+nστ2n(a−1)和u2=σ2n(N−a)u_{1} = \frac{\sigma^{2}+n\sigma_{\tau}^{2}}{n(a-1)} \quad \text{和} \quad u_{2} = \frac{\sigma^{2}}{n(N-a)}

遗憾的是，这一卡方随机变量线性组合的分布无法得到闭式表达式。因此，无法构造 στ2\sigma_{\tau}^{2} 的精确置信区间。Graybill（1961）和 Searle（1971a）给出了近似程序。另见第 13 章 13.6 节。

对于比值 στ2/(στ2+σ2)\sigma_{\tau}^{2}/(\sigma_{\tau}^{2}+\sigma^{2}) 的置信区间，则很容易求出精确表达式。这一比值称为组内相关系数（intraclass correlation coefficient），它反映观测值方差［回忆 V(yij)=στ2+σ2V(y_{ij})=\sigma_{\tau}^{2}+\sigma^{2}］中由处理之间差异造成的比例。为在平衡设计情形下导出这一置信区间，注意 MS处理MS_{\text{处理}} 和 MSEMS_{E} 是相互独立的随机变量，此外可以证明

MS处理/(nστ2+σ2)MSE/σ2∼Fa−1,N−a\frac{MS_{\text{处理}}/(n\sigma_{\tau}^{2}+\sigma^{2})}{MS_{E}/\sigma^{2}} \sim F_{a-1,N-a}

因此

P(F1−α/2,a−1,N−a≤MS处理MSEσ2nστ2+σ2≤Fα/2,a−1,N−a)=1−α(3.54)P\left(F_{1-\alpha/2,a-1,N-a} \leq \frac{MS_{\text{处理}}}{MS_{E}}\frac{\sigma^{2}}{n\sigma_{\tau}^{2}+\sigma^{2}} \leq F_{\alpha/2,a-1,N-a}\right) = 1-\alpha \tag{3.54}

对式 3.54 作整理，可以得到

P(L≤στ2σ2≤U)=1−α(3.55)P\left(L \leq \frac{\sigma_{\tau}^{2}}{\sigma^{2}} \leq U\right) = 1-\alpha \tag{3.55}

其中

L=1n(MS处理MSE1Fα/2,a−1,N−a−1)(3.56a)L = \frac{1}{n}\left(\frac{MS_{\text{处理}}}{MS_{E}}\frac{1}{F_{\alpha/2,a-1,N-a}}-1\right) \tag{3.56a}

和

U=1n(MS处理MSE1F1−α/2,a−1,N−a−1)(3.56b)U = \frac{1}{n}\left(\frac{MS_{\text{处理}}}{MS_{E}}\frac{1}{F_{1-\alpha/2,a-1,N-a}}-1\right) \tag{3.56b}

注意 LL 和 UU 分别是比值 στ2/σ2\sigma_{\tau}^{2}/\sigma^{2} 的 100(1−α)100(1-\alpha) 百分置信下限和置信上限。因此，στ2/(στ2+σ2)\sigma_{\tau}^{2}/(\sigma_{\tau}^{2}+\sigma^{2}) 的 100(1−α)100(1-\alpha) 百分置信区间为

L1+L≤στ2στ2+σ2≤U1+U(3.57)\frac{L}{1+L} \leq \frac{\sigma_{\tau}^{2}}{\sigma_{\tau}^{2}+\sigma^{2}} \leq \frac{U}{1+U} \tag{3.57}

为说明这一程序，我们对例 3.10 的强度数据求出 στ2/(στ2+σ2)\sigma_{\tau}^{2}/(\sigma_{\tau}^{2}+\sigma^{2}) 的 95% 置信区间。回忆 MS处理=29.73MS_{\text{处理}}=29.73、MSE=1.90MS_{E}=1.90、a=4a=4、n=4n=4、F0.025,3,12=4.47F_{0.025,3,12}=4.47、F0.975,3,12=1/F0.025,12,3=1/14.34=0.070F_{0.975,3,12}=1/F_{0.025,12,3}=1/14.34=0.070。因此由式 3.56a 和 b，

L=14[(29.731.90)(14.47)−1]=0.625L = \frac{1}{4}\left[\left(\frac{29.73}{1.90}\right)\left(\frac{1}{4.47}\right)-1\right] = 0.625
U=14[(29.731.90)(10.070)−1]=55.633U = \frac{1}{4}\left[\left(\frac{29.73}{1.90}\right)\left(\frac{1}{0.070}\right)-1\right] = 55.633

由式 3.57，στ2/(στ2+σ2)\sigma_{\tau}^{2}/(\sigma_{\tau}^{2}+\sigma^{2}) 的 95% 置信区间为

0.6251.625≤στ2στ2+σ2≤55.63356.633\frac{0.625}{1.625} \leq \frac{\sigma_{\tau}^{2}}{\sigma_{\tau}^{2}+\sigma^{2}} \leq \frac{55.633}{56.633}

即

0.38≤στ2στ2+σ2≤0.980.38 \leq \frac{\sigma_{\tau}^{2}}{\sigma_{\tau}^{2}+\sigma^{2}} \leq 0.98

我们得出结论：织机之间的变异占所观测到的织物强度变异的 38% 到 98%。由于试验中使用的织机数量很少，这一置信区间相对较宽。然而显然，织机之间的变异（στ2\sigma_{\tau}^{2}）不可忽略。

总均值 μ\mu 的估计。 在许多随机效应试验中，试验者关注总均值 μ\mu 的估计。由基本的模型假定容易看出，任一观测值的期望值就是总均值。因此，总均值的一个无偏估计量为

μ^=y‾..\hat{\mu} = \overline{y}_{..}

所以例 3.10 中总平均强度的估计为

μ^=y‾..=y..N=152716=95.44\hat{\mu} = \overline{y}_{..} = \frac{y_{..}}{N} = \frac{1527}{16} = 95.44

也可以求出总均值的 100(1−α)%100(1-\alpha)\% 置信区间。yˉ\bar{y} 的方差为

V(y‾..)=V(∑i=1a∑j=1nyijan)=nστ2+σ2anV(\overline{y}_{..}) = V\left(\frac{\sum_{i=1}^{a}\sum_{j=1}^{n} y_{ij}}{an}\right) = \frac{n\sigma_{\tau}^{2}+\sigma^{2}}{an}

这一比值的分子由处理均方来估计，因此 V(y‾)V(\overline{y}) 的一个无偏估计量为

V^(y‾..)=MS处理an\hat{V}(\overline{y}_{..}) = \frac{MS_{\text{处理}}}{an}

因此总均值的 100(1−α)%100(1-\alpha)\% 置信区间为

y‾..−tα/2,a(n−1)MS处理an≤μ≤y‾..+tα/2,a(n−1)MS处理an(3.58)\overline{y}_{..}-t_{\alpha/2,a(n-1)}\sqrt{\frac{MS_{\text{处理}}}{an}} \leq \mu \leq \overline{y}_{..}+t_{\alpha/2,a(n-1)}\sqrt{\frac{MS_{\text{处理}}}{an}} \tag{3.58}

为求例 3.10 织物强度试验中总均值的 95% 置信区间，需要 MS处理=29.73MS_{\text{处理}}=29.73 和 t0.025,12=2.18t_{0.025,12}=2.18。置信区间按式 3.58 计算如下：

y‾..−tα/2,a(n−1)MS处理an≤μ≤y‾..+tα/2,a(n−1)MS处理an\overline{y}_{..}-t_{\alpha/2,a(n-1)}\sqrt{\frac{MS_{\text{处理}}}{an}} \leq \mu \leq \overline{y}_{..}+t_{\alpha/2,a(n-1)}\sqrt{\frac{MS_{\text{处理}}}{an}}
95.44−2.1829.7316≤μ≤95.44+2.1829.731695.44-2.18\sqrt{\frac{29.73}{16}} \leq \mu \leq 95.44+2.18\sqrt{\frac{29.73}{16}}
92.47≤μ≤98.4192.47 \leq \mu \leq 98.41

所以，在 95% 置信水平下，该厂织机所生产织物的平均强度在 92.47 与 98.41 之间。这是一个相对较宽的置信区间，因为抽样选取的织机数量很少，而且织机之间存在很大的差异——这体现在织机间差异解释了总变异中很大的一部分。

方差分量的最大似然估计。 本节前面介绍了方差分量估计的方差分析法。该方法应用起来相对直接，并利用了熟悉的量——方差分析表中的均方。然而它也有一些缺点。如前所述，它是一种矩法估计，而数理统计学家一般不喜欢用矩法做参数估计，因为它往往给出不具备良好统计性质的参数估计。一个明显的问题是，它并不总能提供构造所关注方差分量置信区间的简便途径。例如在单因子随机模型中，没有构造 στ2\sigma_{\tau}^{2} 置信区间的简单方法，而 στ2\sigma_{\tau}^{2} 恰恰是试验者主要关注的参数。更受青睐的参数估计技术称为最大似然法（method of maximum likelihood）。该方法的具体实施可能有些繁复，尤其是对试验设计模型而言；但它已被纳入一些支持设计试验的现代计算机软件包中，包括 JMP。

完整介绍最大似然法超出本书范围，但其基本思想可以很容易地说明。假设 xx 是一个概率分布为 f(x,θ)f(x,\theta) 的随机变量，其中 θ\theta 是未知参数。设 x1,x2,…,xnx_{1}, x_{2}, \ldots, x_{n} 是 nn 个观测值的随机样本。样本的联合概率分布为 ∏i=1nf(xi,θ)\prod_{i=1}^{n} f(x_{i},\theta)。似然函数就是把样本观测值视为固定、参数 θ\theta 视为未知时的这一联合概率分布。注意似然函数，例如

L(x1,x2,…,xn;θ)=∏i=1nf(xi,θ)L(x_{1}, x_{2}, \ldots, x_{n};\theta) = \prod_{i=1}^{n} f(x_{i},\theta)

现在只是未知参数 θ\theta 的函数。θ\theta 的最大似然估计量就是使似然函数 L(x1,x2,…,xn;θ)L(x_{1}, x_{2}, \ldots, x_{n};\theta) 达到最大的 θ\theta 值。为说明它如何应用于含随机效应的试验设计模型，设 y\mathbf{y} 是单因子随机效应模型（aa 个处理、nn 次重复）的 an×1an \times 1 观测值向量，Σ\Sigma 是观测值的 an×anan \times an 协方差矩阵。参见 3.9.1 节，那里我们针对 a=3a=3、n=2n=2 的特殊情形导出了这一协方差矩阵。似然函数为

L(x11,x12,…,xa,n;μ,στ2,σ2)=1(2π)N/2∣Σ∣1/2exp⁡[−12(y−jNμ)′Σ−1(y−jNμ)]L(x_{11}, x_{12}, \ldots, x_{a,n};\mu,\sigma_{\tau}^{2},\sigma^{2}) = \frac{1}{(2\pi)^{N/2}|\Sigma|^{1/2}}\exp\left[-\frac{1}{2}(\mathbf{y}-\mathbf{j}_{N}\mu)'\Sigma^{-1}(\mathbf{y}-\mathbf{j}_{N}\mu)\right]

其中 N=anN=an 是观测值总数，jN\mathbf{j}_{N} 是 N×1N \times 1 的全 1 向量，μ\mu 是模型中的总均值。参数 μ\mu、στ2\sigma_{\tau}^{2} 和 σ2\sigma^{2} 的最大似然估计就是使似然函数达到最大的这些量的取值。

最大似然估计量（maximum likelihood estimator，MLE）有一些非常有用的性质。对大样本，它们是无偏的，并且服从正态分布。此外，似然函数的二阶导数矩阵的逆（乘以 -1）就是 MLE 的协方差矩阵。这使得求 MLE 的近似置信区间相对容易。

用于估计方差分量的标准最大似然估计变体称为残差最大似然法（residual maximum likelihood，REML）。它很流行，因为它给出无偏估计量，并且像所有 MLE 一样容易求置信区间。REML 的基本特点是：在估计随机效应时它把模型中的位置参数考虑在内。举一个简单的例子，假设我们希望用最大似然法估计一个正态分布的均值和方差。容易证明 MLE 为

μ^=∑i=1nyin=y‾\hat{\mu} = \frac{\sum_{i=1}^{n} y_{i}}{n} = \overline{y}
σ^2=∑i=1n(yi−y‾)2n\hat{\sigma}^{2} = \frac{\sum_{i=1}^{n} (y_{i}-\overline{y})^{2}}{n}

注意 MLE σ^2\hat{\sigma}^{2} 并不是我们熟悉的样本标准差，它没有把位置参数 μ\mu 的估计考虑在内。REML 估计量则是

S2=∑i=1n(yi−y‾)2n−1S^{2} = \frac{\sum_{i=1}^{n} (y_{i}-\overline{y})^{2}}{n-1}

REML 估计量是无偏的。

表 3.20 例 3.10 织机试验的 JMP 输出

Response Y
Summary of Fit
RSquare0.793521
RSquare Adj0.793521
Root Mean Square Error1.376893
Mean of Response95.4375
Observations (or Sum Wgts)16
Parameter Estimates
TermEstimateStd ErrorDFDent RatioProb > |t|
Intercept95.43751.363111370.01< .0001*
REML Variance Component Estimates
Random EffectVar RatioVar ComponentStd Error95% Lower95% Upper
X13.67032976.95833336.0715247-4.94163618.858303
Residual1.89583330.77397070.97486085.1660065
Total8.8541667
Covariance Matrix of Variance Component Estimates
Random EffectX1Residual
X136.863412-0.149758
Residual-0.1497580.5990307

为说明 REML 方法，表 3.20 给出了例 3.10 织机试验的 JMP 输出。输出中显示了模型参数 μ\mu、στ2\sigma_{\tau}^{2} 和 σ2\sigma^{2} 的 REML 估计。注意方差分量的 REML 估计与前面用方差分析法得到的完全相同。对平衡设计，这两种方法会给出相同结果。然而 REML 输出还包含方差分量的协方差矩阵。该矩阵主对角线元素的平方根就是方差分量的标准误。若 θ^\hat{\theta} 是 θ\theta 的 MLE，σ^(θ^)\hat{\sigma}(\hat{\theta}) 是其估计标准误，则 θ\theta 的近似 100(1−α)100(1-\alpha) 百分置信区间为

θ^−Zα/2σ^(θ^)≤θ≤θ^+Zα/2σ^(θ^)\hat{\theta}-Z_{\alpha/2}\hat{\sigma}(\hat{\theta}) \leq \theta \leq \hat{\theta}+Z_{\alpha/2}\hat{\sigma}(\hat{\theta})

JMP 就用这一途径求出输出中所示的 στ2\sigma_{\tau}^{2} 和 σ2\sigma^{2} 的近似置信区间。REML 给出的 σ2\sigma^{2} 的 95% 置信区间与 3.9 节前面基于卡方算出的区间非常相似。

Footnotes

  1. 假定 {τi}\{\tau_{i}\} 是相互独立的随机变量，这意味着固定效应模型中常用的约束 ∑i=1aτi=0\sum_{i=1}^{a}\tau_{i}=0 并不适用于随机效应模型。

    ↩

  2. 原文此处第二个条件印作 στ2=0\sigma_{\tau}^{2}=0，与前半句重复且与上下文矛盾，按文意应为 στ2>0\sigma_{\tau}^{2}>0。——译者注

    ↩

  3. 原文此处第三个量印作 an2an^{2}，与推导结果不符：按 E(τi2)=στ2E(\tau_{i}^{2})=\sigma_{\tau}^{2}，该项应为 anστ2an\sigma_{\tau}^{2}（即 Nστ2N\sigma_{\tau}^{2}），后面的 Nστ2−nστ2=n(a−1)στ2N\sigma_{\tau}^{2}-n\sigma_{\tau}^{2}=n(a-1)\sigma_{\tau}^{2} 正由此而来。译文按原文印出。——译者注

    ↩

  4. 原文此处印作 MSE=190MS_{E}=190，按上下文应为 MSE=1.90MS_{E}=1.90；同处 χ0.025,122\chi_{0.025,12}^{2} 印作 23,3367，应为 23.3367。——译者注

    ↩
Design and Analysis of Experiments（中文翻译）
3.8 单因子试验的其他例子
Design and Analysis of Experiments（中文翻译）
3.10 方差分析的回归方法