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.4 模型适合性检验

通过方差分析恒等式(式 3.6)对观测值中的变异进行分解,纯粹是一个代数关系。然而,要用这一分解来形式化地检验处理均值无差异,就必须满足某些假定。具体地说,这些假定是:观测值能被模型

yij=μ+τi+ϵijy_{ij} = \mu+\tau_{i}+\epsilon_{ij}

充分描述,且误差是均值为零、方差为常数但未知的正态独立随机变量。如果这些假定成立,方差分析程序就是处理均值无差异这一假设的精确检验。

然而在实践中,这些假定通常不会严格成立。因此,在检验这些假定的有效性之前就依赖方差分析,通常是不明智的。违背基本假定以及模型不适合的情形,都可以通过考察残差(residual)来方便地研究。我们把处理 ii 中第 jj 个观测值的残差定义为

eij=yij−y^ij(3.16)e_{ij} = y_{ij}-\hat{y}_{ij} \tag{3.16}

其中 y^ij\hat{y}_{ij} 是相应观测值 yijy_{ij} 的估计,按下式得到:

y^ij=μ^+τ^i=y‾..+(y‾i.−y‾..)=y‾i.(3.17)\begin{aligned} \hat{y}_{ij} &= \hat{\mu}+\hat{\tau}_{i} \\ &= \overline{y}_{..}+(\overline{y}_{i.}-\overline{y}_{..}) \\ &= \overline{y}_{i.} \end{aligned} \tag{3.17}

式 3.17 给出了一个直观上很合理的结果:第 ii 个处理中任一观测值的估计,就是相应的处理平均。

考察残差应当是任何方差分析自动包含的一部分。如果模型适合,残差应当是无结构的,也就是说,它们不应当含有明显的模式。通过残差分析,可以发现许多类型的模型不适合以及对基本假定的违背。本节说明如何通过对残差作图形分析来方便地进行模型诊断检验,以及如何处理若干常见的异常情形。

3.4.1 正态性假定

检验正态性假定的一种做法是画出残差的直方图。如果误差满足 NID(0,σ2)\mathrm{NID}(0,\sigma^{2}) 假定,该图应当看起来像取自以零为中心的正态分布的样本。遗憾的是,在小样本下直方图的形状往往波动很大,因此出现中度偏离正态的外观并不一定意味着严重违背假定。对正态性的严重偏离则可能很严重,需要作进一步分析。

一个极其有用的做法是构造残差的正态概率图。回忆第 2 章,我们在使用 tt 检验时曾用原始数据的正态概率图来检验正态性假定。在方差分析中,用残差来做这件事通常更有效(也更直接)。如果潜在的误差分布是正态的,该图就会近似一条直线。在目测这条直线时,应更看重图中间部分的点,而不是两端的极点。

表 3.6 给出了例 3.1 刻蚀速率数据的原始数据与残差。正态概率图见图 3.4。考察这一图形得到的总体印象是:误差分布近似正态。正态概率图在左侧略微下弯、在右侧略微上翘的趋势意味着,误差分布的尾部比正态分布所预期的略细一些;也就是说,最大的残差(绝对值)不像预期的那样大。不过,该图并没有严重偏离正态。

一般来说,在固定效应方差分析中,中度偏离正态性关系不大(回忆 3.3.2 节关于随机化检验的讨论)。误差分布的尾部比正态分布厚得多或薄得多,比偏斜分布更值得关注。由于 FF 检验只受到轻微影响,我们说方差分析(以及多重比较等相关程序)对正态性假定是稳健的(robust)。偏离正态性通常会使真实显著性水平和功效都与标称值略有差异,功效一般偏低。我们将在 3.9 节和第 13 章讨论的随机效应模型,受非正态性的影响要严重得多。

表 3.6 例 3.1 的刻蚀速率数据与残差 a

功率 (W)观测值 (j) 12345y^ij=y‾i.\hat{y}_{ij}=\overline{y}_{i.}
16023.8-9.2-21.2-12.218.8
575 (13)542 (14)530 (8)539 (5)570 (4)551.2
180-22.45.62.6-8.422.6
565 (18)593 (9)590 (6)579 (16)610 (17)587.4
200-25.425.6-15.411.63.6
600 (7)651 (19)610 (10)637 (20)629 (1)625.4
22018.0-7.08.0-22.03.0
725 (2)700 (3)715 (15)685 (11)710 (12)707.0

a 每个单元格方框内是残差,括号中的数字表示该次试验运行的实施顺序。

图 3.4 例 3.1 残差的正态概率图

正态概率图上一种很常见的缺陷是:某一个残差远大于其他所有残差。这样的残差通常称为离群值(outlier)。一个或多个离群值的存在会严重扭曲方差分析,因此一旦发现潜在的离群值,就需要仔细查证。离群值的原因常常是计算错误或数据编码、抄录错误。如果不是这些原因,就必须仔细研究该次运行周围的试验情形。如果这个异常的响应是特别理想的值(强度高、成本低等),那么该离群值可能比其余数据更有信息量。除非有合理的非统计依据,我们应当谨慎,不要轻易剔除或舍弃离群观测值。最坏的情况下,你可以保留两个分析结果:一个包含该离群值,一个不包含。

可以用若干形式化的统计程序来检测离群值[例如见 Stefansky(1972)、John 和 Prescott(1975)、Barnett 和 Lewis(1994)]。有些统计软件包会在残差正态概率图上报告正态性统计检验(如 Anderson–Darling 检验)的结果。对此应持谨慎态度,因为这些检验通常假定所应用的数据是独立的,而残差并不独立。

离群值的粗略检查可以通过考察标准化残差(standardized residual)来做:

dij=eijMSE(3.18)d_{ij} = \frac{e_{ij}}{\sqrt{MS_{E}}} \tag{3.18}

如果误差 ϵij\epsilon_{ij} 服从 N(0,σ2)N(0,\sigma^{2}),标准化残差应当近似服从均值为零、方差为 1 的正态分布。因此,约 68% 的标准化残差应落在 ±1\pm1 之内,约 95% 应落在 ±2\pm2 之内,几乎全部应落在 ±3\pm3 之内。偏离零超过 3 或 4 个标准差的残差就是潜在的离群值。

对例 3.1 的刻蚀速率数据,概率图没有显示出任何离群值的迹象。此外,最大的标准化残差为

d1=e1MSE=25.6333.70=25.618.27=1.40d_{1} = \frac{e_{1}}{\sqrt{MS_{E}}} = \frac{25.6}{\sqrt{333.70}} = \frac{25.6}{18.27} = 1.40

这不应引起任何担心。

3.4.2 残差随时间顺序的图形

按数据采集的时间顺序画出残差,有助于发现残差之间的强相关性。出现正残差与负残差成串交替的倾向,说明存在正相关。这意味着误差的独立性假定被违背了。这是一个潜在的严重问题,而且难以纠正,因此在采集数据时尽可能防止这一问题就很重要。对试验进行恰当的随机化是获得独立性的重要一步。

有时试验者(或受试者)的技术会随着试验的推进而变化,或者所研究的过程会“漂移”或变得更加不稳定。这往往导致误差方差随时间发生变化。这种情形常常使残差对时间的图形呈现出一端的散布比另一端更大。非常数方差是一个潜在的严重问题。我们将在 3.4.3 节和 3.4.4 节对此作更多讨论。

表 3.6 给出了刻蚀速率数据的残差及其采集的时间顺序。图 3.5 是这些残差对运行顺序(时间)的图形。没有理由怀疑独立性假定或常数方差假定被违背。

3.4.3 残差对拟合值的图形

如果模型正确、假定满足,残差应当是无结构的;特别地,它们应当与任何其他变量(包括预测响应)都无关。一个简单的检验方法是把残差对拟合值 y^ij\hat{y}_{ij} 作图。(对单因子试验模型,记住 y^ij=y‾i.\hat{y}_{ij} = \overline{y}_{i.},即第 ii 个处理平均。)这一图形不应当显示出任何明显的模式。图 3.6 把例 3.1 刻蚀速率数据的残差对拟合值作图,没有发现异常结构。

这一图形上偶尔会出现的一种缺陷是非常数方差。有时观测值的方差会随着观测值大小的增大而增大。如果试验中的误差或本底噪声是观测值大小的一个恒定百分比,就会出现这种情形。(许多测量仪器都是这样——误差是量程读数的某一百分比。)如果是这样,残差就会随着 yijy_{ij} 增大而变大,残差对 y^ij\hat{y}_{ij} 的图形就会像一个向外张开的漏斗或喇叭。当数据服从非正态的偏斜分布时也会出现非常数方差,因为在偏斜分布中方差往往是均值的函数。

图 3.5 残差对运行顺序(时间)的图形

图 3.6 残差对拟合值的图形

如果违背了方差齐性假定,在平衡的(各处理样本量相等)固定效应模型中,FF 检验只受到轻微影响。然而在非平衡设计中,或者当某个方差远大于其他方差时,问题就更严重。具体地说,如果方差较大的那些因子水平同时也具有较小的样本量,实际第 I 类错误率就会大于预期(或者说置信区间的实际置信水平低于标称值)。反之,如果方差较大的因子水平同时也具有较大的样本量,显著性水平就会小于预期(置信水平更高)。这是尽可能选择等样本量的一个很好理由。对随机效应模型,即使采用平衡设计,误差方差不等也会显著干扰关于方差分量的推断。

方差不等偶尔也会出现在残差对运行顺序的图形上。向外张开的漏斗形说明变异随时间而增大。这可能源于操作者或受试者疲劳、设备上累积的应力、材料性质的变化(例如催化剂劣化)、刀具磨损,或许多其他原因。

当非常数方差源于上述原因时,通常的处理办法是施加一个方差稳定化变换(variance-stabilizing transformation),然后对变换后的数据作方差分析。采用这种做法时应当注意:方差分析的结论适用于变换后的总体来说。

关于如何选择合适的变换已有大量研究。如果试验者知道观测值的理论分布,就可以利用这一信息来选择变换。例如,若观测值服从泊松分布,应使用平方根变换 yij∗=yijy_{ij}^{*} = \sqrt{y_{ij}} 或 yij∗=1+yijy_{ij}^{*} = \sqrt{1+y_{ij}}。若数据服从对数正态分布,对数变换 yij∗=log⁡yijy_{ij}^{*} = \log y_{ij} 是合适的。对表示为分数的二项数据,反正弦变换 yij∗=arcsin⁡yijy_{ij}^{*} = \arcsin\sqrt{y_{ij}} 很有用。当没有明显的变换时,试验者通常凭经验寻找一个能使方差与均值取值无关地变得相等的变换。本节末尾会给出一些指导。在第 5 章介绍的因子试验中,另一种做法是选择使交互作用均方最小的变换,从而得到更容易解释的试验。第 15 章将更详细地讨论解析地选择变换形式的方法。为处理方差不等而作的变换也会影响误差分布的形式,在大多数情况下,变换会使误差分布更接近正态。关于变换的更多讨论,参见 Bartlett(1947)、Dolby(1963)、Box 和 Cox(1964)以及 Draper 和 Hunter(1969)。

方差是否相等的统计检验。 尽管残差图常用来诊断方差不等,人们也提出了若干统计检验。这些检验可以看作对假设

H0:σ12=σ22=⋯=σa2H1:至少有一个 σi2 不满足上述关系\begin{array}{l} H_{0}: \sigma_{1}^{2} = \sigma_{2}^{2} = \cdots = \sigma_{a}^{2} \\ H_{1}: \text{至少有一个 } \sigma_{i}^{2} \text{ 不满足上述关系} \end{array}

的形式化检验。

应用广泛的一种程序是 Bartlett 检验。该程序要计算一个统计量,当 aa 个随机样本来自独立正态总体来说时,该统计量的抽样分布能被自由度为 a−1a-1 的卡方分布很好地近似。检验统计量为

χ02=2.3026qc(3.19)\chi_{0}^{2} = 2.3026\frac{q}{c} \tag{3.19}

其中

q=(N−a)log⁡10Sp2−∑i=1a(ni−1)log⁡10Si2q = (N-a)\log_{10} S_{p}^{2}-\sum_{i=1}^{a} (n_{i}-1)\log_{10} S_{i}^{2}
c=1+13(a−1)(∑i=1a(ni−1)−1−(N−a)−1)c = 1+\frac{1}{3(a-1)}\left(\sum_{i=1}^{a} (n_{i}-1)^{-1}-(N-a)^{-1}\right)
Sp2=∑i=1a(ni−1)Si2N−aS_{p}^{2} = \frac{\sum_{i=1}^{a} (n_{i}-1)S_{i}^{2}}{N-a}

而 Si2S_{i}^{2} 是第 ii 个总体的样本方差。

当样本方差 Si2S_{i}^{2} 相差很大时,量 qq 很大;当所有 Si2S_{i}^{2} 都相等时,qq 等于零。因此,我们应当在 χ02\chi_{0}^{2} 取值过大时拒绝 H0H_{0};也就是说,只有当

χ02>χα,a−12\chi_{0}^{2} > \chi_{\alpha, a-1}^{2}

时才拒绝 H0H_{0},其中 χα,a−12\chi_{\alpha,a-1}^{2} 是自由度为 a−1a-1 的卡方分布的上侧 α\alpha 分位点。也可以采用 PP 值方法作决策。

Bartlett 检验对正态性假定非常敏感。因此,当这一假定的有效性可疑时,就不应使用 Bartlett 检验。

例 3.4

在等离子体刻蚀试验中,正态性假定不成问题,因此我们可以对刻蚀速率数据应用 Bartlett 检验。先计算各处理的样本方差,得 S12=400.7S_{1}^{2}=400.7、S22=280.3S_{2}^{2}=280.3、S32=421.3S_{3}^{2}=421.3、S42=232.5S_{4}^{2}=232.5。于是

Sp2=4(400.7)+4(280.3)+4(421.3)+4(232.5)16=333.7S_{p}^{2} = \frac{4(400.7)+4(280.3)+4(421.3)+4(232.5)}{16} = 333.7
q=16log⁡10(333.7)−4[log⁡10400.7+log⁡10280.3+log⁡10421.3+log⁡10232.5]=0.21\begin{aligned} q ={}& 16\log_{10}(333.7)-4[\log_{10}400.7+\log_{10}280.3 \\ &+\log_{10}421.3+\log_{10}232.5] = 0.21 \end{aligned}
c=1+13(3)(44−116)=1.10c = 1+\frac{1}{3(3)}\left(\frac{4}{4}-\frac{1}{16}\right) = 1.10

检验统计量为

χ02=2.3026(0.21)(1.10)=0.43\chi_{0}^{2} = 2.3026\frac{(0.21)}{(1.10)} = 0.43

由附录表 III 查得 χ0.05,32=7.81\chi_{0.05,3}^{2}=7.81(PP 值为 P=0.934P=0.934),因此我们不能拒绝原假设。没有证据可以反驳全部四个方差相同这一说法[1]。这与分析残差对拟合值的图形所得结论一致。

由于 Bartlett 检验对正态性假定敏感,在某些情形下另一种程序会更有用。Anderson 和 McLean(1974)对相等方差的统计检验作了有益的讨论。修正 Levene 检验[见 Levene(1960)以及 Conover、Johnson 和 Johnson(1981)]是一个很好的程序,它对偏离正态性具有稳健性。为检验所有处理的方差是否相等,修正 Levene 检验使用每个处理中观测值 yijy_{ij} 相对该处理中位数的绝对偏差,记该中位数为 y~i\tilde{y}_{i}。把这些偏差记为

dij=∣yij−y~i∣{i=1,2,…,aj=1,2,…,nid_{ij} = |y_{ij}-\tilde{y}_{i}| \quad \left\{ \begin{array}{ll} i = 1, 2, \ldots, a \\ j = 1, 2, \ldots, n_{i} \end{array} \right.

修正 Levene 检验随后考察这些偏差的均值对所有处理是否相等。可以证明,如果偏差均值相等,则所有处理中观测值的方差都相同。Levene 检验的检验统计量就是把通常用于检验均值相等的方差分析 FF 统计量应用于这些绝对偏差。

例 3.5

一位土木工程师想确定,四种不同的洪水流量频率估计方法应用于同一流域时,是否给出等效的洪峰流量估计。每种方法在该流域上使用六次,所得的流量数据(单位:立方英尺每秒)见表 3.7 的上半部分。表 3.8 总结了该数据的方差分析,结果表明四种方法给出的平均洪峰流量估计存在差异。图 3.7 给出的残差对拟合值的图形令人不安,因为向外张开的漏斗形说明常数方差假定不满足。

我们将对洪峰流量数据应用修正 Levene 检验。表 3.7 的上半部分给出了各处理的中位数 y~i\tilde{y}_{i},下半部分给出了围绕中位数的偏差 dijd_{ij}。Levene 检验就是对 dijd_{ij} 进行标准的方差分析。

由此得到的 FF 检验统计量为 F0=4.55F_{0}=4.55,相应的 PP 值为 P=0.0137P=0.0137。因此,Levene 检验拒绝了方差相等的原假设,基本证实了我们从图 3.7 目视检查得出的判断。洪峰流量数据很适合做数据变换。

表 3.7 洪峰流量数据

估计方法观测值 123456yˉi.\bar{y}_{i.}y~i\tilde{y}_{i}SiS_{i}
10.340.121.230.701.750.120.710.5200.66
20.912.942.142.362.864.552.632.6101.09
36.318.379.756.099.827.247.937.8051.66
417.1511.8210.9517.2014.3516.8214.7215.592.77
估计方法修正 Levene 检验的偏差 dijd_{ij}
10.180.400.710.181.230.40
21.700.330.470.250.251.94
31.4950.5651.9451.7152.0150.565
41.563.774.641.611.241.23

表 3.8 洪峰流量数据的方差分析

变异来源平方和自由度均方F0F_{0}PP 值
方法708.34713236.115776.07< 0.001
误差62.0811203.1041
总计770.428223

图 3.7 例 3.5 中残差对 y^ij\hat{y}_{ij} 的图形

凭经验选择变换。 上面我们看到,如果试验者知道观测值方差与均值之间的关系,就可以利用这一信息来指导选择变换的形式。现在我们对这一点加以展开,并说明一种从数据凭经验选择所需变换形式的方法。

设 E(y)=μE(y)=\mu 为 yy 的均值,并假设 yy 的标准差正比于 yy 的均值的某个幂,即

σy∝μα\sigma_{y} \propto \mu^{\alpha}

我们希望找到一个对 yy 的变换,使方差为常数。假设该变换是原始数据的幂,即

y∗=yλ(3.20)y^{*} = y^{\lambda} \tag{3.20}

则可以证明

σy∗∝μλ+α−1(3.21)\sigma_{y^{*}} \propto \mu^{\lambda+\alpha-1} \tag{3.21}

显然,若取 λ=1−α\lambda = 1-\alpha,变换后数据 y∗y^{*} 的方差就是常数。

前面讨论过的若干常用变换汇总在表 3.9 中。注意 λ=0\lambda=0 对应于对数变换。这些变换按强度递增的顺序排列。所谓变换的强度,指的是它所带来的弯曲程度。温和的变换施于跨度较窄的数据对分析影响很小,而强变换施于较大跨度则可能产生显著的效果。除非 ymax⁡/ymin⁡y_{\max}/y_{\min} 大于 2 或 3,变换往往影响不大。

在许多有重复的试验设计情形中,我们可以从数据凭经验估计 α\alpha。因为在第 ii 个处理组合中 σyi∝μiα=θμiα\sigma_{y_{i}} \propto \mu_{i}^{\alpha} = \theta\mu_{i}^{\alpha},其中 θ\theta 是比例常数,两边取对数得

log⁡σyi=log⁡θ+αlog⁡μi(3.22)\log \sigma_{y_{i}} = \log\theta+\alpha\log\mu_{i} \tag{3.22}

因此,log⁡σyi\log \sigma_{y_{i}} 对 log⁡μi\log \mu_{i} 的图形将是一条斜率为 α\alpha 的直线。由于 σyi\sigma_{y_{i}} 和 μi\mu_{i} 未知,我们可以在式 3.22 中代入它们的合理估计,并把所得直线拟合的斜率作为 α\alpha 的估计。通常我们用第 ii 个处理(更一般地,第 ii 个处理组合或试验条件组)的标准差 SiS_{i} 和平均 y‾i.\overline{y}_{i.} 来估计 σyi\sigma_{y_{i}} 和 μi\mu_{i}。

为考察对例 3.5 的洪峰流量数据使用方差稳定化变换的可能性,我们在图 3.8 中画出 log⁡Si\log S_{i} 对 log⁡y‾i.\log \overline{y}_{i.} 的图形。穿过这四个点的直线斜率接近 1/2,由表 3.9 可知这意味着平方根变换可能是合适的。变换后数据 y∗=yy^{*}=\sqrt{y} 的方差分析见表 3.10,残差对预测响应的图形见图 3.9。与图 3.7 相比,这一残差图有了很大改善,因此我们得出结论:平方根变换是有帮助的。注意表 3.10 中我们把误差和总和的自由度各减少了一个,以计入利用数据估计变换参数 α\alpha 这一事实。

表 3.9 方差稳定化变换

σy\sigma_{y} 与 μ\mu 的关系α\alphaλ=1−α\lambda = 1-\alpha变换说明
σy∝\sigma_{y} \propto 常数01不作变换
σy∝μ1/2\sigma_{y} \propto \mu^{1/2}1/21/2平方根泊松(计数)数据
σy∝μ\sigma_{y} \propto \mu10对数
σy∝μ3/2\sigma_{y} \propto \mu^{3/2}3/2−1/2-1/2倒数平方根
σy∝μ2\sigma_{y} \propto \mu^{2}2-1倒数

图 3.8 例 3.5 洪峰流量数据的 log⁡Si\log S_{i} 对 log⁡yˉi.\log\bar{y}_{i.} 的图形

图 3.9 变换后数据的残差对 y^ij∗\hat{y}_{ij}^{*} 的图形(例 3.5 洪峰流量数据)

表 3.10 变换后洪峰流量数据 y∗=yy^{*}=\sqrt{y} 的方差分析

变异来源平方和自由度均方F0F_{0}PP 值
方法32.6842310.894776.99< 0.001
误差2.6884190.1415
总计35.372622

在实践中,许多试验者通过简单试用若干备选变换、并观察每种变换对残差对预测响应图形的影响,来选择变换的形式,然后选用产生最满意残差图的那个变换。作为替代,还有一种称为 Box-Cox 方法的形式化方法,用来选择方差稳定化变换。第 15 章将讨论并说明这一程序。它被广泛使用,并在许多软件包中实现。

3.4.4 残差对其他变量的图形

如果还收集了任何可能影响响应的其他变量的数据,就应当把残差对这些变量作图。例如,在例 3.1 的刻蚀速率试验中,刻蚀速率可能显著受气压影响,因此应当画出残差与气压的关系图。如果数据是用不同的刻蚀机采集的,就应当画出残差对这些机器的关系图。这类残差图中出现模式,意味着该变量影响响应。这提示我们:在未来的试验中应当更仔细地控制该变量,或者把它纳入分析。

Footnotes
  1. 原文此处写作“all five variances”(全部五个方差),但试验只有四个处理,应为四个方差。——译者注