10  方差分析

10.1 引言

目前市面上有 4 款软件来教授用户学习一门新的变成语言,某大公司正在考虑大规模购买这 4 款软件中的一款。公司内部一些资深人士声称,这几款软件没有太大的区别,它们对用户的最终学习效果影响很小。为了检验这一假设,公司决定选择 160 名工程师,并将他们分成 4 组,每组 40 人。第 \(i\) 组的工程师使用第 \(i\)\(i = 1, 2, 3, 4\))款软件来学习新的语言。当所有工程师完成学习后,参与学习的所有工程师将进行一次全面考试。公司希望使用这次考试的结果来确定不同的教学软件是否真的可以互换。该公司需要怎么做才能实现这一点?

在回答这个问题之前,我们需要明确,显然我们希望得出这样的结论:当所有分组的平均成绩相似时,可以认为这几款软件确实是可以互换的;而当不同分组的平均成绩之间有很大差异时,可以认为这几款软件在本质上是存在区别的。要得出这样的结论,我们必须注意到:将这 160 名工程师分成 4 组的方法至关重要。例如,假设第一组的成员成绩显著高于其他组,我们能从中得出什么结论?具体来说,这个结果是由于 第 1 款软件的教学效果更好,还是第一组的工程师的学习能力更好?为了能够得出这几款软件是否可以互换的结论,需要一种合理的分组方式以保证不同分组的工程师不会存在显著的差异,这一点至关重要。以完全随机的方式将工程师分成 4 组是久经考验的一种方法。也就是说,完全随机的分组方式可以使得所有可能的分组组合都有相同的概率,从而保障任何一组明显优于其他组的可能性非常小。假设我们确实使用“随机”的方式对工程师分组,虽然目前还不清楚如何才能做到。一种有效的方法是首先对这 160 名工程师进行 1~160 的编号,然后生成 1~160 的随机排列,并将排列中的前 40 名工程师放分到 第 1 组,排列中的第 41 到 80 名工程师分到 第 2 组,以此类推。

现在,我们可以合理的假设:每一位工程师的考试成绩应该大致服从正态分布,该分布的参数取决于他所使用的教学软件。同样,也可以合理的假设:尽管工程师的平均考试成绩取决于他所使用的教学软件,但考试成绩的变异性(variability)将源自 160 名工程师的内在变异(inherent variation),而不是来自所使用的软件。因此,如果我们令 \(X_{ij},\ i=1,...4, \ j=1,...,40\) 表示第 \(i\) 组中第 \(j\) 名工程师的考试成绩。一个合理的模型可能是假设 \(X_{ij}\) 是独立的随机变量,其中 \(X_{ij}\) 服从均值为 \(\mu_i\) 且方差为 \(\sigma^2\) 的正态分布。不同的教学软件可以互换的假设相当于对原假设 \(H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\) 进行检验。

变异量

正如在 小节 9.5 中对 变异 这个术语的解释,此处再次做如下的强调。

变异量the amount of variation)是指一组数据中的变异性,即数据的波动或分散程度。这种波动可以是由多种因素引起的,包括个体差异、测量误差、随机误差等。在后续的翻译章节中,我们会根据具体的场景交叉使用 变异量数据波动,如无特殊说明,这些指的都是一个概念。

同时,我们会避免使用 数据差异difference) 来指代 变异量 的概念,因为在数学中,数据差异 更多指代两个数据之间的差值,因此和我们所讨论的数据的波动和分布程度并不是一个概念。

在本章中,我们将介绍一种技术可以用来检验如上案例中的假设。这种技术非常通用,可以用于推断与总体均值相关的大量参数,这种技术称之为 方差分析ANOVA, analysis of variance)。

10.2 综述

本章内容还需要继续修改。

小节 8.4小节 8.5.1小节 8.6.1小节 8.7.1 章节中,我们研究了关于两个总体均值的假设检验。在本章中,我们将研究关于多个总体均值的假设检验。

小节 10.3 中,假设我们有来自 \(m\) 个不同总体的样本,每个样本的大小为 \(n\),我们希望使用这些数据来检验 \(m\) 个总体均值相等的假设。对于这种场景,由于随机变量的均值仅取决于一个因素——即随机变量来自哪个样本——因此称这种情况为 单因素方差分析one-way analysis of variance)。在 小节 10.3.2 中,我们展示了如何获得 \(\binom{m}{2}\) 对总体均值之间的 多重比较multiple comparisons)。在第 小节 10.3.3 中展示了当 \(m\) 个样本之间的样本大小不同时,如何检验总体均值相等的假设。

小节 10.4小节 10.5 中,我们将介绍由两个因素决定因变量均值的模型。在这些模型中,我们可以把自变量排列在一个矩阵中,在这个矩阵中,每一行和每一列都代表着一个因素或条件,矩阵中的元素代表着自变量的观测数据。此时,某个因变量的均值同时依赖于矩阵中的行和列。这种模型称为 双因素方差分析two-way analysis of variance)。我们假设因变量的均值受到自变量的 加性方式additive fashion)影响,具体地,第 \(i\) 行、第 \(j\) 列的变量对应的响应均值可以表示为 \(\mu + \alpha_i + \beta_j\)。在 小节 10.4 中,我们将介绍如何估计 \(\alpha_i\)\(\beta_j\)。在 小节 10.5 中将介绍如何进行检验假设,即某个特定因素不会影响结果的均值。

小节 10.6 中,我们将介绍因变量的均值允许以非线性方式取决于自变量所在的行和列的情况,从而允许两个因素之间存在 交互作用interaction)。我们还将介绍如何检验不存在交互作用的假设,以及关于不存在行效应和列效应对因变量均值影响的假设。

在本章中,我们假设所有的数据都服从正态分布,并且都具有相同(虽然未知)的方差 \(\sigma^2\)。使用 ANOVA 来检验和总体均值有关的多个参数的原假设 \(H_0\),需要基于从共同的总体方差 $ ^2 $ 中推导出的两个估计量。第一个估计量是 \(\sigma^2\) 的估计量,并且无论原假设是否为真,该估计量都是 \(\sigma^2\) 的有效估计量。第二个估计量也是 \(\sigma^2\) 的估计量,只是仅在 \(H_0\) 为真时,该估计量才是 \(\sigma^2\) 的有效估计量。此外,当 \(H_0\) 不成立时,第二个估计量往往会超过 \(\sigma^2\)。检验方法就是比较这两个估计量的值,并在第二个估计量与第一个估计量的比率足够大时拒绝 \(H_0\)。换句话说,当 \(H_0\) 为真时,这两个估计量应彼此接近(因为它们在这种情况下都是 \(\sigma^2\) 的估计量),而当 \(H_0\) 不为真时,第二个估计量通常应比第一个大,因此当第二个估计量显著大于第一个时,拒绝 \(H_0\) 是自然的选择。

我们将使用和卡方随机变量有关的事实来获得方差 \(\sigma^2\) 的估计量。假设 \(X_1, \dots, X_N\) 是独立的正态分布随机变量,它们可能具有不同的均值但却具有相同的方差 \(\sigma^2\),并且设 \(\mu_i = E[X_i]\),其中 \(i = 1, \dots, N\)

由于随机变量

\(Z_i = (X_i - \mu_i)/\sigma, \quad i = 1, \dots, N\)

服从标准正态分布,因此根据卡方分布随机变量的定义可知:

\[ \sum_{i=1}^{N} Z_i^2 = \sum_{i=1}^{N} \frac{(X_i - \mu_i)^2}{\sigma^2} \tag{10.1}\]

是一个自由度为 \(N\) 的卡方分布随机变量。假设每一个 \(\mu_i, \ i = 1, \dots, N\) 都可以表示为 \(k\) 个未知参数的线性函数。同时,假设我们可以确定这 \(k\) 个参数的估计量,从而给出 \(\mu_i\) 的估计量。如果我们用 \(\hat{\mu}_i\) 表示 \(\mu_i, \ i=1,...,N\) 的估计量,那么可以证明:

\(\sum_{i=1}^{N} \frac{(X_i - \hat{\mu}_i)^2}{\sigma^2}\)

服从自由度为 \(N - k\) 的卡方分布。

换句话说,

\(\sum_{i=1}^{N} \frac{(X_i - E[X_i])^2}{\sigma^2}\)

是一个自由度为 \(N\) 的卡方分布随机变量。

如果我们用 \(k\) 个参数的线性函数来表示每个 \(E[X_i]\),并用这 \(k\) 个参数的估计量来替换这 \(k\) 个参数,那么我们得到的表达式仍然是卡方分布,但是其自由度将因为每个参数由其估计量来代替而减少 1。

为了说明如上的结论,我们考虑所有均值已知并且相等的情况,即:

\(E[X_i] = \mu, \quad i = 1, \dots, N\)

因为此时只有一个需要估计的参数,因此 \(k = 1\)。把公共均值 \(\mu\) 的估计量 \(\overline{X}\) 代入 式 10.1 得到如下的量:

\[ \sum_{i=1}^{N} \frac{(X_i - \overline{X})^2}{\sigma^2} \tag{10.2}\]

式 10.2 是一个自由度为 \(N - 1\) 的卡方分布随机变量。

在所有均值都相等的情况下,数据 \(X_1, \dots, X_N\) 均来自同一个正态分布总体,因此 式 10.2 等价与 \((N - 1)S^2/\sigma^2\),其中 \(S^2\) 是样本方差。换句话说,在这种场景下的结论就是:\((N - 1)S^2/\sigma^2\) 是自由度为 \(N - 1\) 的卡方分布随机变量(参见 定理 6.2)。

10.3 单因素方差分析

\(m\) 个独立的样本,每个样本的大小为 \(n\),其中第 \(i\) 个样本——\(X_{i1}, X_{i2}, \dots, X_{in}\)——是均值 \(\mu_i\)、方差 \(\sigma^2\) 均未知的正态分布随机变量。也就是说,

\(X_{ij} \sim \mathcal{N}(\mu_i, \sigma^2), \quad i = 1, \dots, m, \quad j = 1, \dots, n\)

我们将对以下假设进行检验:

\(H_0: \mu_1 = \mu_2 = \dots = \mu_m \quad vs \quad H_1: \text{并非所有的均值都相等}\)

也就是说,我们检验的原假设为所有总体的均值都是相等的,对应的备择假设为:至少存在两个总体的均值不同。如何理解我们所要进行的假设检验?想象一下,对于同一个物品而言,我们有 \(m\) 种不同的处理方法,第 \(i\) 种处理方法的结果是具有均值为 \(\mu_i\)、方差为 \(\sigma^2\) 的正态分布随机变量。然后我们感兴趣的是:把每种处理方法应用于 \(n\) 个物品并分析结果,以检验 所有处理方法的的效果相同 的假设。

总共有 \(nm\) 个独立的正态分布随机变量 \(X_{ij}\),因此 \(X_{ij}\) 标准化之后的随机变量的平方和将是 自由度为 \(nm\) 的卡方分布随机变量。

\[ \sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - E[X_{ij}])^2}{\sigma^2} = \sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - \mu_i)^2}{\sigma^2} \sim \chi^2_{nm} \tag{10.3}\]

为了获得 \(m\) 个未知参数 \(\mu_1, \dots, \mu_m\) 的估计量,设 \(X_i.\) 表示第 \(i\) 个样本中所有元素的平均值,即:

\(X_i. = \frac{1}{n} \sum_{j=1}^{n} X_{ij}\)

\(X_i.\) 是第 \(i, \ i = 1, \dots, m\) 个总体的样本均值,因此是总体均值 \(\mu_i\) 的估计量。因此,用 \(\mu_i\) 的估计量 \(X_i.\) 替代 式 10.3 中的 \(\mu_i\),得到:

\[ \sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - X_i.)^2}{\sigma^2} \tag{10.4}\]

式 10.4 是自由度为 \(nm - m\) 的卡方分布。(每估计一个参数,自由度就会减少 1。)

\(SS_W = \sum_{i=1}^{m} \sum_{j=1}^{n} (X_{ij} - X_i.)^2\)

因此,式 10.4 变为 \(SS_W/\sigma^2\)

由于卡方分布随机变量的期望值等于其自由度,因此:\(E[SS_W/\sigma^2] = nm - m\),即 \(E\left[\frac{SS_W}{nm - m}\right] = \sigma^2\)

我们因此得到了方差 \(\sigma^2\) 的第一个估计量,即 \(\frac{SS_W}{nm - m}\)。注意,对于该估计量而言,我们并没有对原假设的真假作任何的假设。也就是说,无论原假设是否为真,我们都将有该估计量。

定义 10.1 组内平方和within samples sum of squares

\(SS_W = \sum_{i=1}^{m} \sum_{j=1}^{n} (X_{ij} - X_{i.})^2\)

因为 \(SS_W\) 是利用样本均值代替 式 10.3 中的总体均值而得到的统计量,因此,我们称 \(SS_W\)组内平方和within samples sum of squares)。

\(\frac{SS_W}{nm - m}\)\(\sigma^2\) 的一个估计量。\(\blacksquare\)

方差 \(\sigma^2\) 的第二个估计量是只有在原假设为真时的估计量。假设 \(H_0\) 为真,因此所有总体均值 \(\mu_i\) 相等,即 \(\mu_i = \mu\)。此时,所有 \(m\) 个样本均值 \(X_1., X_2., \dots, X_m.\) 都将具有相同的均值 \(\mu\) 和相同的方差 \(\sigma^2/n\)。因此,对 \(X_i.\) 标准化后的随机变量 \(\frac{X_i. - \mu}{\sqrt{\sigma^2 / n}} = \frac{\sqrt{n}(X_i. - \mu)}{\sigma}\) 的平方和是一个自由度为 \(m\) 的卡方分布随机变量。

即当 \(H_0\) 为真时,

\[ \frac{n \sum_{i=1}^{m} {(X_i. - \mu)^2}}{\sigma^2} \sim \chi^2_m \tag{10.5}\]

当所有的总体均值都等于 \(\mu\) 时,\(\mu\) 的估计量是所有 \(nm\) 个数据值的平均值。也就是说,\(\mu\) 的估计量是 \(X_{..}\),即:

\(X_{..} = \frac{1}{nm} \sum_{i=1}^{m} \sum_{j=1}^{n} X_{ij} = \frac{1}{m} \sum_{i=1}^{m} X_i.\)

\(X_{..}\) 代替 式 10.5 中的位置参数 \(\mu\),则当 \(H_0\) 为真时,有:

\(\frac{n \sum_{i=1}^{m} {(X_i. - X_{..})^2}}{\sigma^2} \sim \chi_{m-1}^2\)

如果我们定义 \(SS_B\)

\(SS_B = n \sum_{i=1}^{m} (X_i. - X_{..})^2\)

那么,当 \(H_0\) 为真时,\(SS_B/\sigma^2\) 是自由度为 \(m - 1\) 的卡方分布随机变量。

所以,当 \(H_0\) 为真时,有:

\(E\left[\frac{SS_B}{\sigma^2}\right] = m - 1\)

即,

\[ E\left[\frac{SS_B}{(m - 1)}\right] = \sigma^2 \tag{10.6}\]

因此,当 \(H_0\) 为真时,\(\frac{SS_B}{m - 1}\) 也是 \(\sigma^2\) 的估计量。

定义 10.2 组间平方和between samples sum of squares

\(SS_B = n \sum_{i=1}^{m} (X_i. - X_{..})^2\)

称为 组间平方和between samples sum of squares)。当 \(H_0\) 为真时,\(\frac{SS_B}{m - 1}\)\(\sigma^2\) 的估计量。\(\blacksquare\)

因此,我们得到:

\(\begin{align} \frac{SS_W}{nm - m}, \quad &\text{总是} \sigma^2 \text{的估计量} \\ \frac{SS_B}{m - 1}, \quad &\text{仅当} H_0 \text{为真时} \sigma^2 \text{的估计量} \end{align}\)

可以证明(在本节的末尾会对其进行证明),当 \(H_0\) 不为真时,\(\frac{SS_B}{m - 1}\) 往往会超过 \(\sigma^2\),因此合理的做法是将检验统计量定义为:

\(TS = \frac{SS_B/(m - 1)}{SS_W/(nm - m)}\)

并且当 \(TS\) 足够大时,拒绝 \(H_0\)

接下来的问题就是:当 \(TS\) 需要多大时才能拒绝 \(H_0\)?可以证明:如果 \(H_0\) 为真,那么 \(SS_B\)\(SS_W\) 是独立的。这意味着,当 \(H_0\) 为真时,\(TS\) 服从分子自由度为 \(m - 1\)、分母自由度为 \(nm - m\)\(F-\text{分布}\)。我们使用符号 \(F_{r,s}\) 表示分子自由度为 \(r\)、分母自由度为 \(s\)\(F-\text{分布}\) 随机变量,令 \(F_{\alpha, m-1, nm-m}\) 表示这个 \(F-\text{分布}\)\(100(1 - \alpha)\%\) 分位数,即:

\(P\left \{F_{m-1, nm-m} > F_{\alpha, m-1, nm-m}\right\} = \alpha\)

在显著性水平 \(\alpha\) 时,\(H_0\) 的检验为:

\[ \begin{align} \text{拒绝} H_0&, \quad \text{如果} \frac{SS_B/(m - 1)}{SS_W/(nm - m)} \gt F_{\alpha, m-1, nm-m} \\ \text{不拒绝} H_0&, \quad \text{其他} \end{align} \tag{10.7}\]

检验 所有总体均值是否相等 的假设的另一种计算方法是计算 \(p-\text{value}\)。如果检验统计量 \(TS = v\),则 \(p-\text{value}\) 由以下公式给出:

\(p-\text{value} = P\{F_{m-1,nm-m} \geq v\}\)

在 R 中可以使用 pf() 计算 \(p-\text{value}\)

# 方式 1
1 - pf(v, m - 1, nm - m)

# 方式 2
pf(v, m - 1, nm - m, lower.tail = FALSE)
平方和恒等式the sum of squares identity

如下的 代数恒等式algebraic identity),称为 平方和恒等式sum of squares identity),在手动计算时非常有用。

\(\sum_{i=1}^{m} \sum_{j=1}^{n} X_{ij}^2 = nm X^2_{..} + SS_B + SS_W\)

根据 定义 10.2,我们定义 \(SS_B\) 为:

\(SS_B = n \sum_{i=1}^{m} (X_i. - X_{..})^2\)

在手动计算时,应首先计算 \(SS_B\) 的值。一旦得到 \(SS_B\) 后,可以根据 平方和恒等式 确定 \(SS_W\)。在计算 \(SS_W\) 的同时,我们也得到了 \(\sum_{i=1}^{m} \sum_{j=1}^{n} X_{ij}^2\)\(\overline{X}^2_{..}\)。然后根据如下公式确定 \(SS_W\)

\(SS_W = \sum_{i=1}^{m} \sum_{j=1}^{n} X_{ij}^2 - nm X^2_{..} - SS_B\)

习题 10.1 一家汽车租赁公司使用 15 台相同的发动机在固定速度下测试 3 种不同的汽油品牌。每种汽油品牌都分配了其中的 5 台发动机,并且每台发动机分配了10 加仑汽油。以下的数据表示了不同发动机在汽油耗尽时所获得的总里程数:

汽油品牌 里程(英里)
汽油 1 220, 251, 226, 246, 260
汽油 2 244, 235, 232, 242, 225
汽油 3 252, 272, 250, 238, 256

在 5% 的显著性水平下,检验假设:发动机的平均里程不受汽油品牌的影响。

解 10.1. 我们手动计算 习题 10.1 的结果。

首先需要注意的是,每个数据值减去一个常数不会影响检验统计量的值。因此,我们从每个数据值中减去 220,得到以下信息:

汽油品牌 里程(英里) \(\sum_j X_{ij}\) \(\sum_j X_{ij}^2\)
1 0, 31, 6, 26, 40 103 3273
2 24, 15, 12, 22, 5 78 1454
3 32, 52, 30, 18, 36 168 6248

因为 \(m = 3\)\(n = 5\),并且

\(\begin{align} X_1. &= \frac{103}{5} = 20.6 \\ X_2. &= \frac{78}{5} = 15.6 \\ X_3. &= \frac{168}{5} = 33.6 \\ X_{..} &= \frac{(103 + 78 + 168)}{15} = 23.2667 \\ X_{..}^2 &= 541.3393 \end{align}\)

因此,

\(SS_B = 5[(20.6 - 23.2667)^2 + (15.6 - 23.2667)^2 + (33.6 - 23.2667)^2] = 863.3335\)

此外,

\(\sum \sum X_{ij}^2 = 3273 + 1454 + 6248 = 10975\)

根据平方和恒等式,

\(SS_W = 10975 - 15(541.3393) - 863.3335 = 1991.5785\)

因此,检验统计量的值为:

\(TS = \frac{863.3335 / 2}{1991.5785/12} = 2.60\)

使用 R 语言计算得到 \(p-\text{value}\)

1 - pf(2.60, 2, 12)
[1] 0.1153232

我们同样可以使用 R 来计算整个过程:

gas1 <- c(0, 31, 6, 26, 40)
gas2 <- c(24, 15, 12, 22, 5)
gas3 <- c(32, 52, 30, 18, 36)
data <- matrix(c(gas1, gas2, gas3), ncol = 5, byrow = TRUE)
x1. <- mean(data[1,])
x2. <- mean(data[2,])
x3. <- mean(data[3,])
x.. <- mean(data)
SS_B <- 5 * sum(apply(data, 1, function(x) (mean(x) - mean(data))^2))
x_ij_ss <- sum(data^2)
SS_W <- x_ij_ss - 3 * 5 * x..^2 - SS_B
TS = (SS_B / 2) / (SS_W / 12)
1 - pf(TS, 2, 12)
[1] 0.1152489

由于 \(p-\text{value}\) 大于 0.05,因此在 5% 的显著性水平下,我们不能拒绝汽品牌不会影响发动机里程的原假设。\(\blacksquare\)

现在让我们证明:

\(E\left[\frac{SS_B}{m-1}\right] \geq \sigma^2\)

当且仅当 \(H_0\) 为真时,等号成立。

因此,我们必须证明:

\(E\left[\frac{\sum_{i=1}^{m} (X_{i.} - X_{..})^2}{m-1}\right] \geq \frac{\sigma^2}{n}\)

当且仅当 \(H_0\) 为真时,等号成立。

\(\mu_. = \frac{\sum_{i=1}^{m} \mu_i}{m}\)\(m\) 个总体均值的平均值。此外,对于 \(i = 1, \dots, m\),令:

\(Y_i = X_{i.} - \mu_i + \mu_.\)

由于 \(X_{i.}\) 服从均值为 \(\mu_i\)、方差为 \(\frac{\sigma^2}{n}\) 的正态分布,因此 \(Y_i\) 服从均值为 \(\mu_.\)、方差为 \(\frac{\sigma^2}{n}\) 的正态分布。因此,\(Y_1, \dots, Y_m\) 构成了一个方差为 \(\frac{\sigma^2}{n}\) 的正态分布总体的样本。

\(Y_{.}\)\(Y_i\) 的均值:

\(Y_{.} = \frac{\sum_{i=1}^{m} Y_i}{m} = X_{..} - \mu_. + \mu_. = X_{..}\)

则:

\(X_{i.} - X_{..} = Y_i + \mu_i - \mu_. - Y_{.}\)

因此,

\(\begin{align} E\left[\sum_{i=1}^{m} (X_{i.} - X_{..})^2\right] &= E\left[\sum_{i=1}^{m} (Y_i - Y_{.} + \mu_i - \mu_.)^2\right] \\ &= E\left[\sum_{i=1}^{m} \left[(Y_i - Y_{.})^2 + (\mu_i - \mu_.)^2 + 2(\mu_i - \mu_.)(Y_i - Y_{.})\right]\right] \\ &= E\left[\sum_{i=1}^{m} (Y_i - Y_{.})^2\right] + \sum_{i=1}^{m} (\mu_i - \mu_.)^2 + 2\sum_{i=1}^{m} (\mu_i - \mu_.) E[Y_i - Y_{.}] \\ &= (m-1)\frac{\sigma^2}{n} + \sum_{i=1}^{m} (\mu_i - \mu_.)^2 + 2\sum_{i=1}^{m} (\mu_i - \mu_.) E[Y_i - Y_{.}] \\ &= (m-1)\frac{\sigma^2}{n} + \sum_{i=1}^{m} (\mu_i - \mu_.)^2\end{align}\)

在倒数第二个等式中,\(\sum_{i=1}^{m} \frac{(Y_i - Y_{.})^2}{(m-1)}\) 是总体方差 \(\frac{\sigma^2}{n}\) 的无偏估计量,并且最后一个等式是因为 \(E[Y_i - Y_{.}] = E[Y_i] - E[Y_{.}] = \mu_. - \mu_. = 0\)。最后,在结果的等式两边都除以 \(m-1\) 得到:

\(E\left[\frac{\sum_{i=1}^{m} (X_{i.} - X_{..})^2}{m-1}\right] = \frac{\sigma^2}{n} + \frac{\sum_{i=1}^{m} (\mu_i - \mu_.)^2}{(m-1)}\)

因为 \(\sum_{i=1}^{m} (\mu_i - \mu_.)^2 \geq 0\),当且仅当所有 \(\mu_i\) 相等时等号成立,所以 \(E\left[\frac{SS_B}{m-1}\right] \geq \sigma^2\)

表 10.1 总结了本节的结果:

表 10.1: 单因素方差分析表;在显著性水平 \(\alpha\) 下的检验,当 \(TS \geq F_{\alpha, m-1, nm-m}\) 时拒绝 \(H_0\),否则不拒绝 \(H_0\);如果 \(TS = v\),则 \(p-\text{value} = P\{F_{m-1, nm-m} \geq v\}\)
变异来源 平方和 自由度 检验统计量的值
组间 \(SS_B = n\sum_{i=1}^{m} (X_{i.} - X_{..})^2\) \(m-1\)
组内 \(SS_W = \sum_{i=1}^{m} \sum_{j=1}^{n} (X_{ij} - X_{i.})^2\) \(nm-m\)
\(TS = \frac{SS_B/(m-1)}{SS_W/(nm-m)}\)

10.3.1 使用 R 计算统计量

首先,需要注意的是:在 R 中定义矩阵时,矩阵元素默认按照列的方向来填充。例如,假设我们想要用 R 构建如下的矩阵并命名为 A

\[\begin{pmatrix} 2 & 1 & 3 & 4 \\ 1 & 5 & 2 & 6 \\ 5 & 3 & 7 & 9 \end{pmatrix} \]

我们可以按如下方式进行:

A <- matrix(c(2, 1, 5, 1, 5, 3, 3, 2, 7, 4, 6, 9), nrow = 3, ncol = 4)

R 中 matrix() 的定义是:在给出矩阵元素之后,给出行数,然后给出列数。因此,我们可以省略 nrow =ncol =,直接定义 A 如下:

A <- matrix(c(2, 1, 5, 1, 5, 3, 3, 2, 7, 4, 6, 9), 3, 4)

要显示矩阵,我们会输入:

A <- matrix(c(2, 1, 5, 1, 5, 3, 3, 2, 7, 4, 6, 9), 3, 4)
A
     [,1] [,2] [,3] [,4]
[1,]    2    1    3    4
[2,]    1    5    2    6
[3,]    5    3    7    9

当然,也可以通过列向量来初始化一个矩阵,例如:

v1 <- c(2, 1, 5)
v2 <- c(1, 5, 3)
v3 <- c(3, 2, 7)
v4 <- c(4, 6, 9)
A <- matrix(c(v1, v2, v3, v4), 3, 4)
A
     [,1] [,2] [,3] [,4]
[1,]    2    1    3    4
[2,]    1    5    2    6
[3,]    5    3    7    9

如果 \(A\) 是一个矩阵,\(x\) 是一个标量(scalar),那么 \(A + x\) 是一个矩阵,其结果是将 \(x\) 加到 \(A\) 的每个元素上。矩阵的幂运算意味着矩阵中的每个元素都执行该幂次方运算。例如,如果 \(A\) 是矩阵 \([a_{i,j}]\)),那么 \(A^2\) 就是矩阵 \([a_{i,j}^2]\)。(此处的这种约定与“矩阵乘法”的概念并不一致。)

要对矩阵 \(A\) 的所有元素求和,只需输入 R 命令 sum(A);要获得矩阵 \(A\) 中所有元素的平均值,使用 mean(A)。例如,如果我们想要计算矩阵 \(A\) 中所有元素的平方求和,并且想将结果命名为 ss,我们可以输入:

A <- matrix(c(2, 1, 5, 1, 5, 3, 3, 2, 7, 4, 6, 9), 3, 4)
ss <- sum(A^2)
ss
[1] 260

如果我们想计算矩阵 \(A\) 中的每一行、或者每一列的元素的和,可以使用 rowSums()colSums()。例如:

A <- matrix(c(2, 1, 5, 1, 5, 3, 3, 2, 7, 4, 6, 9), 3, 4)
ss <- sum(A^2)
ss
[1] 260
rowSums(A)
[1] 10 14 24
colSums(A)
[1]  8  9 12 19

如果想计算矩阵 \(A\) 中的每一行、或者每一列的元素的平均值,可以使用 rowMeans()colMeans()

假设 \(A\) 是一个具有 \(m\) 行和 \(n\) 列的矩阵,并且对于 \(b_1, \dots, b_m\),我们想要得到一个矩阵 \(B\)。矩阵 \(B\) 是从 \(A\) 的第一行的每个元素中减去 \(b_1\),从 \(A\) 的第二行的每个元素中减去 \(b_2\),……可以使用 R 中的 sweep() 来完成该操作:

A <- matrix(...)
b <- c(b1,...,bm)
B <- sweep(A, 1, b)

例如:

A <- matrix(c(2, 3, 3, 5, 1, 4), nrow = 2)

# 对矩阵 A 的每一行中的所有元素,都减去该行对应的均值
b <- rowMeans(A)
B <- sweep(A, 1, b)
B
     [,1] [,2] [,3]
[1,]    0    1   -1
[2,]   -1    1    0

同样,假设 \(A\) 是一个具有 \(m\) 行和 \(n\) 列的矩阵,并且对于 \(c_1, \dots, c_n\),我们想要得到一个矩阵 \(C\)。矩阵 \(C\) 是从 \(A\) 的第一列的每个元素中减去 \(c_1\),从 \(A\) 的第二列的每个元素中减去 \(c_2\),……同样可以使用 R 中的 sweep() 来完成该操作:

A <- matrix(...)
c <- c(c1,...,cn)
C <- sweep(A, 2, c)

例如我们想让矩阵 \(A\) 的每一列的所有元素都加上该列对应的均值:

A <- matrix(c(2, 3, 3, 5, 1, 4), nrow = 2)

# 对矩阵 A 的每一列中的所有元素,都加上该列对应的均值
c <- colMeans(A)
C <- sweep(A, 2, -c)
C
     [,1] [,2] [,3]
[1,]  4.5    7  3.5
[2,]  5.5    9  6.5

接下来,我们将使用 R 来解决 习题 10.1 中所示的单因素方差分析问题。

gas1 <- c(220, 251, 226, 246, 260)
gas2 <- c(244, 235, 232, 242, 225)
gas3 <- c(252, 272, 250, 238, 256)
G <- matrix(c(gas1, gas2, gas3), nrow = 3, byrow = TRUE)

m <- nrow(G)
n <- ncol(G)

rm <- rowMeans(G)
mm <- mean(G)
d <- rm - mm
SS_B <- n * sum(d^2)
SS_W <- sum(G^2) - m * n * mm^2 - SS_B
TS = (SS_B / (m - 1)) / (SS_W / (m * n - m))
1 - pf(TS, m - 1, m * n - m)
[1] 0.1152489

10.3.2 样本均值的多重比较

当我们拒绝了总体均值相等的原假设时,我们需要对不同总体的样本均值 \(\mu_1, \dots, \mu_m\) 进行比较。实现这个目的常用方法是 \(T-\text{方法}\)T-method)。对于指定的 \(\alpha\) 值,\(T-\text{方法}\) 为所有的 \(\binom{m}{2}\) 个差异 \(\mu_i - \mu_j\)\(i \neq j, i, j = 1, \dots, m\))提供了 联合置信区间joint condifence intervals),从而使得这些所有的置信区间都将以 \(1 - \alpha\) 的概率包含其对应的 \(\mu_i - \mu_j\)\(T-\text{方法}\) 基于以下结果:

对于任何 \(i \neq j\)\(\mu_i - \mu_j\) 都将以 \(1 - \alpha\) 的概率包含于如下的区间:

\(X_i. - X_j. - W < \mu_i - \mu_j < X_i. - X_j. + W\)

其中,

\(W = \frac{1}{\sqrt{n}}C(m, nm - m, \alpha)\sqrt{\frac{SS_W}{nm - m}}\)

并且,附录中的 表 A.1 给出了 \(\alpha = 0.05\)\(\alpha = 0.01\) 时,\(C(m, nm - m, \alpha)\) 的值。

习题 10.2 某大学声称,对于该大学而言,来自三个不同高中的学生的第一学年的平均成绩没有差异。下表给出了随机选择的 12 名学生的第一学年的平均成绩,其中每个高中选取 4 名学生。假设显著性水平为 5%,这些数据是否推翻了该大学的主张?如果是,请确定来自不同高中的学生的平均成绩之间的差异的置信区间,以便我们可以在 95% 的置信水平下相信所有的置信区间都是成立的。

学校 1 学校 2 学校 3
3.2 3.4 2.8
3.4 3.0 2.6
3.3 3.7 3.0
3.5 3.3 2.7

解 10.2. 首先,共有 \(m = 3\) 个样本,每个样本的大小为 \(n = 4\)。使用 R 进行计算得到的结果如下:

school1 <- c(3.2, 3.4, 3.3, 3.5)
school2 <- c(3.4, 3.0, 3.7, 3.3)
school3 <- c(2.8, 2.6, 3.0, 2.7)

G <- matrix(c(school1, school2, school3), nrow = 3, byrow = TRUE)
m <- nrow(G)
n <- ncol(G)

rm <- rowMeans(G)
mm <- mean(G)
d <- rm - mm
SS_B <- n * sum(d^2)
SS_W <- sum(G^2) - m * n * mm^2 - SS_B
TS = (SS_B / (m - 1)) / (SS_W / (m * n - m))
1 - pf(TS, m - 1, m * n - m)
[1] 0.00480163

因此,我们会拒绝:来自三个高中的学生在第一学年的平均成绩相等的假设。

为了确定总体均值差异的置信区间,首先注意到样本均值为:

\(X_1 = 3.350, \quad X_2 = 3.350, \quad X_3 = 2.775\)

从附录 表 A.1 中,我们看到 \(C(3, 9, 0.05) = 3.95\),因此 \(W = \frac{3.95}{\sqrt{4}} \cdot \sqrt{0.0431} = 0.410\),我们得到以下置信区间:

\(\begin{align} -0.410 < &\mu_1 - \mu_2 < 0.410 \\ 0.165 < &\mu_1 - \mu_3 < 0.985 \\ 0.165 < &\mu_2 - \mu_3 < 0.985 \end{align}\)

因此,在 95% 的置信水平下,我们可以得出结论:

  • 第 3 所高中的学生平均成绩比第 1 所或第 2 所高中的学生的平均成绩低,差异在 0.165 到 0.985 之间
  • 第 1 所高中和第 2 所高中的学生的平均成绩差异小于 0.410。\(\blacksquare\)

10.3.3 样本大小不同的单因素方差分析

在上一节的模型中,我们都假设每个样本的样本大小均相同。虽然在理想情况下,样本大小都是相同的,但实际中却无法总是保证不同样本之间的样本大小都相同。假设我们有 \(m\) 个样本,样本的大小分别为 \(n_1, n_2, \dots, n_m\)。也就是说,数据由独立随机变量 \(X_{ij}\) 组成,其中 \(i = 1, \dots, m\)\(j = 1, \dots, n_i\),并且:

\(X_{ij} \sim \mathcal{N}(\mu_i, \sigma^2)\)

我们仍然希望检验假设 \(H_0\),即所有样本的均值都相等。

\(\sum_{i=1}^{m} \sum_{j=1}^{n_i} \frac{(X_{ij} - E[X_{ij}])^2}{\sigma^2} = \sum_{i=1}^{m} \sum_{j=1}^{n_i} \frac{(X_{ij} - \mu_i)^2}{\sigma^2}\)

是自由度为 \(\sum_{i=1}^{m} n_i\) 的卡方分布随机变量。因此,用 \(\mu_i\) 的估计量 \(X_{i.}\)(第 \(i\) 个样本中的元素的平均值)代替每个总体的均值 \(\mu_i\),我们得到:

\(\sum_{i=1}^{m} \sum_{j=1}^{n_i} \frac{(X_{ij} - X_{i.})^2}{\sigma^2}\)

是自由度为 \(\sum_{i=1}^{m} n_i - m\) 的卡方分布。因此,定义:

\(SS_W = \sum_{i=1}^{m} \sum_{j=1}^{n_i} (X_{ij} - X_{i.})^2\)

\(SS_W / \left( \sum_{i=1}^{m} n_i - m \right)\)\(\sigma^2\) 的无偏估计量。

此外,如果 \(H_0\) 为真,并且 \(\mu\) 是公同均值,则随机变量 \(X_{i.}\) 将是独立的正态分布随机变量,且满足:

\(E[X_{i.}] = \mu, \quad \text{Var}(X_{i.}) = \frac{\sigma^2}{n_i}\)

因此,当 \(H_0\) 为真时:

\(\sum_{i=1}^{m} \frac{(X_{i.} - \mu)^2}{\sigma^2 / n_i} = \sum_{i=1}^{m} \frac{n_i (X_{i.} - \mu)^2}{\sigma^2}\)

是自由度为 \(m\) 的卡方分布随机变量。因此,用估计量 \(X_{..}\)(所有 \(X_{ij}\) 的平均值)代替 \(\mu\),得到统计量:

\(\sum_{i=1}^{m} \frac{n_i (X_{i.} - X_{..})^2}{\sigma^2}\)

是自由度为 \(m - 1\) 的卡方分布。因此,定义:

\(SS_B = \sum_{i=1}^{m} n_i (X_{i.} - X_{..})^2\)

可以得出,当 \(H_0\) 为真时,\(SS_B / (m - 1)\) 也是 \(\sigma^2\) 的无偏估计量。

可以证明,当 \(H_0\) 为真时,\(SS_B\)\(SS_W\) 是独立的,所以统计量:

\(\frac{SS_B / (m - 1)}{SS_W / \left( \sum_{i=1}^{m} n_i - m \right)}\)

是一个分子自由度为 \(m - 1\)、分母自由度为 \(\sum_{i=1}^{m} n_i - m\)\(F-\text{分布}\) 随机变量。因此,我们可以得出显著性水平 为 \(\alpha\) 下的零假设的检验结论:

\[ \begin{align} \text{拒绝} H_0,& \quad \text{如果} \frac{SS_B / (m - 1)}{SS_W / \left( \sum_{i=1}^{m} n_i - m \right)} > F_{m-1, N, \alpha} \\ \text{接受} H_0,& \quad \text{否则} \end{align} \]

其中,\(N = \sum_{i=1}^{m} n_i - m\)

备注

当样本大小不同时,我们称之为 不均衡场景unbalanced case)。因为在样本均衡设计中,检验统计量对于总体方差不相同的场景相对不那么敏感。因此,样本均衡设计通常更稳健。所以,与不均衡场景相比,在设计检验时,我们需要尽可能的选择均衡的样本。

10.4 双因素方差分析:引言和参数估计

虽然我们可以用 小节 10.3 中的模型研究数据集中的单个因素的影响,但我们还可以研究多个因素的影响。在本节中,我们假设每个数据值都受到两个因素的影响。

例 10.1 对 5 名学生进行了 4 种不同标准的阅读测试,测试的结果如下表所示。表中的 20 个数据受到两个因素的影响:考试因素和学生因素。考试因素有 4 个可能的取值或 水平levels),学生因素有 5 个可能的取值或 水平levels)。

考试标准 学生 1 学生 2 学生 3 学生 4 学生 5
1 75 73 60 70 86
2 78 71 64 72 90
3 80 69 62 70 85
4 73 67 63 80 92

通常情况下,我们假设第一个因素有 \(m\) 个可能的 水平,第二个因素有 \(n\) 个可能的 水平。令 \(X_{ij}\) 表示当第一个因素处于第 \(i\) 个水平,第二个因素处于第 \(j\) 个水平时获得的值。我们经常将数据集表示为如下的数组:

\[ \begin{matrix} X_{11} & X_{12} & \dots & X_{1j} & \dots & X_{1n} \\ X_{21} & X_{22} & \dots & X_{2j} & \dots & X_{2n} \\ \vdots & \vdots & \dots & \vdots & \dots & \vdots \\ X_{i1} & X_{i2} & \dots & X_{ij} & \dots & X_{in} \\ \vdots & \vdots & \dots & \vdots & \dots & \vdots \\ X_{m1} & X_{m2} & \dots & X_{mj} & \dots & X_{mn} \end{matrix} \]

因此,我们称第一个因素为 “行” 因素,称第二个因素为 “列” 因素。

小节 10.3 中,我们假设数据 \(X_{ij}, i = 1, \dots, m; j = 1, \dots, n\) 是具有相同方差 \(\sigma^2\) 的独立正态分布随机变量。然而,在第 小节 10.3 中,我们假设,只有一个因素(即它所属样本的均值)会影响数据点的均值。而在本节中,我们则假设,数据点的均值以 加性方式additive manner)取决于其所在的行和列。

小节 10.3 的模型中,我们令 \(X_{ij}\) 表示第 \(i\) 个样本中第 \(j\) 个成员的值,那么模型可以符号化表示为:

\(E[X_{ij}] = \mu_i\)

如果我们令 \(\mu\) 表示 \(\mu_i\) 的平均值,即

\(\mu = \frac{1}{m} \sum_{i=1}^{m} \mu_i\)

那么我们可以将模型重写为

\(E[X_{ij}] = \mu + \alpha_i\)

其中 \(\alpha_i = \mu_i - \mu\)。根据 \(\alpha_i\) 的定义:\(\mu_i\) 相对于平均值 \(\mu\) 的偏差,我们很容易得到:

\(\sum_{i=1}^{m} \alpha_i = 0\)

一个 双因素加性模型 也可以用行偏差和列偏差来表示。如果我们令 \(\mu_{ij} = E[X_{ij}]\),则加性模型假设对于某些常数 \(a_i\)\(i = 1, \dots, m\))和 \(b_j\)\(j = 1, \dots, n\))有:

\(\mu_{ij} = a_i + b_j\)

令:

\(\mu_{i.} = \frac{1}{n} \sum_{j=1}^{n} \mu_{ij}, \quad \mu_{.j} = \frac{1}{m} \sum_{i=1}^{m} \mu_{ij}, \quad \mu_{..} = \frac{1}{mn} \sum_{i=1}^{m} \sum_{j=1}^{n} \mu_{ij}\)

同样地,我们定义

\(a_{.} = \frac{1}{m} \sum_{i=1}^{m} a_i, \quad b_{.} = \frac{1}{n} \sum_{j=1}^{n} b_j\)

因为:

\(\mu_{i.} = \frac{1}{n} \sum_{j=1}^{n} (a_i + b_j)/n = a_i + b_{.}\)

同理,

\(\mu_{.j} = a_{.} + b_j, \quad \mu_{..} = a_{.} + b_{.}\)

如果我们令:

\(\mu = \mu_{..} = a_{.} + b_{.}\)

\(\alpha_i = \mu_{i.} - \mu = a_i - a_{.}\)

\(\beta_j = \mu_{.j} - \mu = b_j - b_.\)

则我们有:

\(\mu_{ij} = E[X_{ij}] = \mu + \alpha_i + \beta_j\)

其中,\(\sum_{i=1}^{m}{\alpha_i} = \sum_{j=1}^{n}{\beta_j} = 0\)

我们称 \(\mu\) 被为 总均值the grand mean),\(\alpha_i\) 为由第 \(i\) 行导致的相对于 总均值 的偏差,\(\beta_j\) 为由第 \(j\) 列导致的相对于 总均值 的偏差。

接下来,我们确定参数 \(\mu, \alpha_i, \beta_j\)\(i = 1, \dots, m\), \(j = 1, \dots, n\))的估计量。

\(X_{i.} = \frac{1}{n} \sum_{j=1}^{n} X_{ij} \quad \text{表示第} i \text{行中的值的平均值}\)

\(X_{.j} = \frac{1}{m} \sum_{i=1}^{m} X_{ij} \quad \text{表示第} j \text{列中的值的平均值}\)

\(X_{..} = \frac{1}{nm} \sum_{i=1}^{m} \sum_{j=1}^{n} X_{ij} \quad \text{表示所有数据值的平均值}\)

则:

\[ \begin{align} E[X_{i.}] &= \frac{1}{n} \sum_{j=1}^{n} E[X_{ij}] \\ &= \mu + \frac{1}{n} \sum_{j=1}^{n} \alpha_i + \frac{1}{n} \sum_{j=1}^{n} \beta_j \\ &= \mu + \alpha_i \quad \because \sum_{j=1}^{n} \beta_j = 0 \end{align} \]

同理有:

\(E[X_{.j}] = \mu + \beta_j\)

\(E[X_{..}] = \mu\)

因为上述等式等价于:

\[ \begin{align} E[X_{..}] &= \mu \\ E[X_{i.} - X_{..}] &= \alpha_i \\ E[X_{.j} - X_{..}] &= \beta_j \end{align} \]

我们发现,参数 \(\mu, \alpha_i, \beta_j\) 的无偏估计量,记作 \(\hat{\mu}, \hat{\alpha_i}, \hat{\beta_j}\) 为:

\[ \begin{align} \hat{\mu} &= X_{..} \\ \hat{\alpha_i} &= X_{i.} - X_{..} \\ \hat{\beta_j} &= X_{.j} - X_{..} \end{align} \]

\(\blacksquare\)

习题 10.3 下表的数据来自 例 10.1,请使用这些数据估计模型的参数。

考试 学生 1 学生 2 学生 3 学生 4 学生 5 行总计 \(X_{i.}\)
1 75 73 60 70 86 364 72.8
2 78 71 64 72 90 375 75
3 80 69 62 70 85 366 73.2
4 73 67 63 80 92 375 75
列总计 306 280 249 292 353 1480 (grand total)
\(X_{.j}\) 76.5 70 62.25 73 88.25 \(X_{..} = \frac{1480}{20} = 74\)

解 10.3. 估计量为:

\(\hat{\mu} = 74\)

\(\hat{\alpha_1} = 72.8 - 74 = -1.2 \quad \hat{\beta_1} = 76.5 - 74 = 2.5\)

\(\hat{\alpha_2} = 75 - 74 = 1 \quad \hat{\beta_2} = 70 - 74 = -4\)

\(\hat{\alpha_3} = 73.2 - 74 = -0.8 \quad \hat{\beta_3} = 62.25 - 74 = -11.75\)

\(\hat{\alpha_4} = 75 - 74 = 1 \quad \hat{\beta_4} = 73 - 74 = -1\)

\(\hat{\beta_5} = 88.25 - 74 = 14.25\)

因此,如果随机选择一名学生并随机指定一场考试,那么我们对其平均分的估计是 \(\hat{\mu} = 74\)。如果我们知道该学生选择了考试 \(i\),则这会使我们对其平均分的估计增加 \(\hat{\alpha_i}\);如果我们知道学生的编号是 \(j\),则这会使我们对其平均分的估计增加 \(\hat{\beta_j}\)

例如,我们会估计第 1 场考试中第 2 名学生的分数是一个随机变量,其均值为:

\(\hat{\mu} + \hat{\alpha_1} + \hat{\beta_2} = 74 - 1.2 - 4 = 68.8\)

\(\blacksquare\)

10.5 双因素方差分析:检验假设

考虑一个双因素模型,数据为 \(X_{ij}\)\(i = 1, \dots, m\)\(j = 1, \dots, n\)。假设 \(X_{ij}\) 是具有共同方差 \(\sigma^2\) 的独立正态分布随机变量,且其均值满足以下条件:

\(E[X_{ij}] = \mu + \alpha_i + \beta_j\)

其中,

\(\sum_{i=1}^{m} \alpha_i = \sum_{j=1}^{n} \beta_j = 0\)

在本节中,我们将关注检验如下的假设:

\(H_0: \text{所有的} \ \alpha_i = 0 \quad vs \quad H_1: \text{并非所有的} \ \alpha_i \ \text{都等于} 0\)

该原假设表明 \(X_{ij}\) 不存在行效应,即某个数据的值不受其所在行的因素水平的影响。

类似的,我们还将关注有关列的假设,即

\(H_0: \text{所有的} \ \beta_j = 0 \quad vs \quad H_1: \text{并非所有的} \ \beta_j \ \text{都等于} 0\)

我们将使用 方差分析analysis of variance) 的方法来检验上述原假设。在 方差分析 中,我们会得到方差 \(\sigma^2\) 的两个不同的估计量:

  • 第一个估计量总是有效估计量。
  • 第二个估计量只有在原假设为真时才是有效估计量。
  • 此外,当原假设不为真时,第二个估计量往往大于 \(\sigma^2\)

为了获得 \(\sigma^2\) 的第一个估计量,我们需要基于如下的事实:

\(\sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - E[X_{ij}])^2}{\sigma^2} = \sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - \mu - \alpha_i - \beta_j)^2}{\sigma^2}\)

是自由度为 \(nm\) 的卡方随机变量。如果用未知参数的估计量替换对应的未知参数 \(\mu, \alpha_1, \alpha_2, \dots, \alpha_m, \beta_1, \beta_2, \dots, \beta_n\),则结果仍然是一个卡方分布随机变量,但其自由度将随着估计参数的增加而对应的减少(每估计一个参数,自由度减 1)。由于 \(\sum_{i=1}^{m} \alpha_i = \sum_{j=1}^{n} \beta_j = 0\),这意味着当我们估计了 \(m-1\)\(\alpha_i\) 后,最后一个 \(\alpha_i\) 也随之确定;同理,我们仅需要估计 \(n-1\)\(\beta_j\) 就可以确定所有的 \(\beta_j\)。最后,因为 \(\mu\) 也需要估计,所以需要估计的参数总数为:

\(1 + m - 1 + n - 1 = n + m - 1\)

因此,

\(\sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - \hat{\mu} - \hat{\alpha_i} - \hat{\beta_j})^2}{\sigma^2}\)

是一个自由度为 \(nm - (n + m - 1) = (n-1)(m-1)\) 的卡方分布随机变量。

因为 \(\hat{\mu} = X_{..}, \hat{\alpha_i} = X_{i.} - X_{..}, \hat{\beta_j} = X_{.j} - X_{..}\),所以有:

\[ \sum_{i=1}^{m} \sum_{j=1}^{n} \frac{(X_{ij} - X_{i.} - X_{.j} + X_{..})^2}{\sigma^2} \tag{10.8}\]

是一个自由度为 \((n-1)(m-1)\) 的卡方分布随机变量。

定义 10.3 误差平方和error sum of squares

我们称统计量 \(SS_e\)误差平方和

\(SS_e = \sum_{i=1}^{m} \sum_{j=1}^{n} (X_{ij} - X_{i.} - X_{.j} + X_{..})^2\) \(\blacksquare\)

如果我们认为某个值与我们所估计的均值之间的差异是“误差”(error),那么 \(SS_e\) 就等于这些误差的平方的总和。由于 \(SS_e / \sigma^2\) 就是 式 10.8,所以我们知道 \(SS_e / \sigma^2\) 是一个自由度为 \((n-1)(m-1)\) 的卡方分布随机变量。因为卡方分布随机变量的期望值等于其自由度,所以:

\(E[SS_e/\sigma^2] = (n-1)(m-1)\)

即:

\(E\left[\frac{SS_e}{(n-1)(m-1)}\right] = \sigma^2\)

所以:

\(\frac{SS_e}{(n-1)(m-1)}\)\(\sigma^2\) 的无偏估计量。

假设我们现在想要检验 不存在行效应 的原假设,即我们想要检验

\(H_0: \text{所有的} \alpha_i = 0 \quad vs \quad H_1: \text{并不是所有的} \alpha_i \text{都为 0}\)

为了获得 \(\sigma^2\) 的第二个估计量,考虑行均值 \(X_{i.}\)\(i = 1, \ldots, m\)。当 \(H_0\) 为真时,每个 \(\alpha_i = 0\),因此

\(E[X_{i.}] = \mu + \alpha_i = \mu\)

因为每个 \(X_{i.}\)\(n\) 个随机变量的均值,而每个随机变量的方差为 \(\sigma^2\),因此有

\(\text{Var}(X_{i.}) = \frac{\sigma^2}{n}\)

因此当 \(H_0\) 为真时:

\(\sum_{i=1}^{m}(X_{i.} - E[X_{i.}])^2 / \text{Var}(X_{i.}) = n \sum_{i=1}^{m}(X_{i.} - \mu)^2 / \sigma^2\)

将服从自由度为 \(m\) 的卡方分布。如果我们现在用 \(X_{..}\)\(\mu\) 的估计量)代替 \(\mu\),那么得到的表达式仍将是卡方分布,但自由度会减少 1。因此,当 \(H_0\) 为真时,有以下结论:

\(n \sum_{i=1}^{m}(X_{i.} - X_{..})^2 / \sigma^2\)

服从自由度为 \(m-1\) 的卡方分布。

定义 10.4 行平方和row sum of squares

我们称统计量 \(SS_r\)行平方和

\(SS_r = n \sum_{i=1}^{m}(X_{i.} - X_{..})^2\)

\(H_0\) 为真时,\(\frac{SS_r}{\sigma^2}\) 服从自由度为 \(m-1\) 的卡方分布。因此,当 \(H_0\) 为真时,有

\(E[SS_r / \sigma^2] = m-1\)

即:

\(E[SS_r / (m-1)] = \sigma^2\)

此外,可以证明,当 \(H_0\) 不为真时,\(\frac{SS_r}{(m-1)}\) 的值将大于 \(\sigma^2\)

因此,我们再次获得了 \(\sigma^2\) 的两个估计量:

  • 无论原假设是否为真,第一个估计量 \(\frac{SS_e}{(n-1)(m-1)}\) 都是有效估计量。
  • 而仅在 \(H_0\) 为真时,第二个估计量 \(\frac{SS_r}{(m-1)}\) 才是有效估计量;当 \(H_0\) 不为真时,\(\frac{SS_r}{(m-1)}\) 将大于 \(\sigma^2\)

基于 \(\sigma^2\) 的两个估计量(\(SS_e\)\(SS_r\))的比值对 \(H_0\) 进行检验。具体来说,我们使用如下的检验统计量:

\(TS = \frac{SS_r/(m-1)}{SS_e/[(n-1)(m-1)]}\)

可以证明,当 \(H_0\) 为真时,估计量是独立的,所以在显著性水平为 \(\alpha\) 下的检验为:

\(\begin{align} \text{拒绝} \quad &H_0, \quad \text{如果} \ TS \geq F_{\alpha, m-1, (n-1)(m-1)} \\ \text{不拒绝} \quad &H_0, \quad \text{否则}\end{align}\)

或者,可以通过计算检验统计量 \(v\)\(p-\text{value}\) 来进行检验:

\(p-\text{value} = P(F_{m-1,(n-1)(m-1)} \geq v)\)

类似的,我们同样可以使用如上的方法检验列效应的原假设,即检验所有的 \(\beta_j = 0\) 的原假设。表 10.2表 10.3 对本节介绍的 行效应列效应 的假设检验方法进行了总结。同时,可以使用 R 来计算相关的 \(p-\text{value}\)

表 10.2: 行平方和,列平方和以及误差平方和
平方和 自由度
\(SS_r = n \sum_{i=1}^{m} (X_{i.} - X_{..})^2\) \(m - 1\)
\(SS_c = m \sum_{j=1}^{n} (X_{.j} - X_{..})^2\) \(n - 1\)
误差 \(SS_e = \sum_{i=1}^{m} \sum_{j=1}^{n} (X_{ij} - X_{i.} - X_{.j} + X_{..})^2\) \((n - 1)(m - 1)\)
表 10.3: 双因素方差分析(Two-Factor ANOVA),\(N = (n-1)(m-1)\)
原假设 检验统计量 显著性水平 \(\alpha\) 下的检验 \(TS=v\) 时的 \(p-\text{value}\)
所有 \(\alpha_i = 0\) \(\frac{SS_r/(m-1)}{SS_e/N}\) \(TS \geq F_{\alpha,m-1,N}\),则拒绝 \(H_0\) \(P\{F_{m-1,N} \geq v\}\)
所有 \(\beta_j = 0\) \(\frac{SS_c/(n-1)}{SS_e/N}\) \(TS \geq F_{\alpha,n-1,N}\),则拒绝 \(H_0\) \(P\{F_{n-1,N} \geq v\}\)

例 10.2 以下数据表示在热排放附近的 6 个站点收集的、大型无脊椎动物的不同种类的数量,数据收集的时间为 1970 年至 1977 年。

注记

热排放通常是指来自工业设施(如发电厂、工厂或其他热力系统)的废热,通常以温水或蒸汽的形式排放到附近的水体或空气中。因为水温的上升可能改变栖息地的条件,影响物种的多样性和数量,所以这种排放可能会对周围的环境、生态系统以及水生生物产生影响。

Year Station 1 Station 2 Station 3 Station 4 Station 5 Station 6
1970 53 35 31 37 40 43
1971 36 34 17 21 30 18
1972 47 37 17 31 45 26
1973 55 31 17 23 43 37
1974 40 32 19 26 45 37
1975 52 42 20 27 26 32
1976 39 28 21 21 36 28
1977 40 32 21 21 36 35

我们使用 R 来检验以下假设:

    1. 数据是否随时间变化
    1. 数据是否在各站点之间变化
v1 <- c(53, 36, 47, 55, 40, 52, 39, 40)
v2 <- c(35, 34, 37, 31, 32, 42, 28, 32)
v3 <- c(31, 17, 17, 17, 19, 20, 21, 21)
v4 <- c(37, 21, 31, 23, 26, 27, 21, 21)
v5 <- c(40, 32, 45, 43, 45, 26, 36, 36)
v6 <- c(43, 18, 26, 37, 37, 32, 28, 35)
X <- matrix(c(v1, v2, v3, v4, v5, v6), 8, 6)
H <- sweep(X, 1, rowMeans(X)) + mean(X)
C <- sweep(H, 2, colMeans(X))
SSe <- sum(C^2)
D <- rowMeans(X) - mean(X)
SSr <- ncol(X) * sum(D^2)
E <- colMeans(X) - mean(X)
SSc <- nrow(X) * sum(E^2)
N <- (nrow(X) - 1) * (ncol(X) - 1)
TSrow <- (SSr / (nrow(X) - 1)) / (SSe / N)
TSrow
[1] 3.602264
1 - pf(TSrow, nrow(X) - 1, N)
[1] 0.005037051
TScol = (SSc / (ncol(X) - 1)) / (SSe / N)
TScol
[1] 22.67062
1 - pf(TScol, ncol(X) - 1, N)
[1] 4.413647e-10

因此,数据分布是否不依赖于年份的 F 统计量的值为 3.602,其对应的 \(p-\text{value}\) 为 0.005;数据分布是否不依赖于站点的 F 统计量的值为 22.67,其对应的 \(p-\text{value}\)\(4 \times 10^{-10}\)。因此,在几乎所有的显著性水平下,这两个原假设我们都将拒绝。 \(\blacksquare\)

10.6 有交互效应的双因素方差分析

小节 10.4小节 10.5 中,我们讨论了观察数据的分布依赖于两个因素——称之为 行因素列因素——的场景。具体来说,我们假设 \(X_{ij}\) 的均值可以表示为两部分的和——一部分取决于数据所在的行,另一部分取决于数据所在的列。即我们假设

\(X_{ij} \sim \mathcal{N}(\mu + \alpha_i + \beta_j, \sigma^2), \quad i = 1, \ldots, m, \quad j = 1, \ldots, n\)

然而,上述模型的问题是:上述模型总是假设 行因素列因素 的效应是可加的,它不允许 行因素列因素 之间存在 交互效应

例如,考虑如下的实验,在该实验中,我们需要比较四个不同的工人在使用三台不同机器时,生产的次品数量的均值。在分析结果时,我们可能假设,对于不同的工人而言,使用某台机器产生的次品数量的增量是相同的。但是,也存在如下的可能:同一机器对于不同工人而言所产生的次品数量的增量并不一致,也就是说,可能存在 \(\text{工人-机器}\)交互效应interation)。而在加性模型中,不允许存在这种 交互效应

交互效应

交互效应或者交互作用指的是:在某个自变量的不同水平上,另一个自变量对因变量的效应不同。例如上面提到的,在工人这个自变量的不同水平上,同一台机器对于产品次品数量的效应并不相同。

为了考虑行和列之间可能的 交互效应,令

\(\mu_{ij} = E[X_{ij}]\)

并定义 \(\mu, \alpha_i, \beta_j, \gamma_{ij}, i = 1, \ldots, m; j = 1, \ldots, n\)

\(\begin{align} \mu &= \mu_{..} \\ \alpha_i &= \mu_{i.} - \mu_{..} \\ \beta_j &= \mu_{.j} - \mu_{..} \\ \gamma_{ij} &= \mu_{ij} - \mu_{i.} - \mu_{.j} + \mu_{..} \end{align}\)

显然有:

\(\mu_{ij} = \mu + \alpha_i + \beta_j + \gamma_{ij}\)

并且,也可以很容易的证明:

\(\sum_{i=1}^{m} \alpha_i = \sum_{j=1}^{n} \beta_j = \sum_{i=1}^{m} \gamma_{ij} = \sum_{j=1}^{n} \gamma_{ij} = 0\)

  • 参数 \(\mu\)\(nm\) 个均值的平均值,我们称之为总体均值grand mean)。
  • 参数 \(\alpha_i\) 表示第 \(i\) 行的均值与 总体均值 之间的差值,我们被之称为 行效应effect of row)。
  • 参数 \(\beta_j\) 表示第 \(j\) 列的均值与 总体均值 之间的差值,我们称之为 列效应effect of column)。
  • 参数 \(\gamma_{ij} = \mu_{ij} - (\mu + \alpha_i + \beta_j)\) 表示 \(\mu_{ij}\)总体均值 和 第 \(i\) 行效应、第 \(j\) 列效应 的和的差值,是对第 \(i\) 行和第 \(j\) 列的加性偏离的度量,我们称之为第 \(i\) 行和第 \(j\) 列的 交互效应

为了能够检验假设:行和列之间没有交互效应——也就是说,所有 \(\gamma_{ij} = 0\)——每对因素需要有多个观察值。假设每行和每列都有 \(l\) 个观察值。也就是说,假设数据是 \(\{X_{ijk}, i=1, \dots, m; j=1, \dots, n; k=1, \dots, l\}\),其中 \(X_{ijk}\) 是第 \(i\) 行和第 \(j\) 列中的第 \(k\) 个观察值。由于假设所有观察值都是相互独立的正态分布随机变量且具有相同的方差 \(\sigma^2\),所以:

\(X_{ijk} \sim \mathcal{N}(\mu + \alpha_i + \beta_j + \gamma_{ij}, \sigma^2)\)

其中:

\[ \sum_{i=1}^{m} \alpha_i = \sum_{j=1}^{n} \beta_j = \sum_{i=1}^{m} \gamma_{ij} = \sum_{j=1}^{n} \gamma_{ij} = 0 \tag{10.9}\]

我们将对上述参数进行估计,并检验以下原假设:

\(\begin{align} H_0^r: \alpha_i &= 0, \quad \text{for all} \ i \\ H_0^c: \beta_j &= 0, \quad \text{for all} \ j \\ H_0^{int}: \gamma_{ij} &= 0, \quad \text{for all} \ i, j \end{align}\)

也就是说,\(H_0^r\) 是没有行效应的假设,\(H_0^c\) 是没有列效应的假设,\(H_0^{int}\) 是没有行和列之间交互效应的假设。

根据 式 10.9\(E[X_{ijk}] = \mu_{ij} = \mu + \alpha_i + \beta_j + \gamma_{ij}\) 可以得到:

\(\begin{align} E[X_{ij.}] &= \mu_{ij} = \mu + \alpha_i + \beta_j + \gamma_{ij} \\ E[X_{i..}] &= \mu + \alpha_i \\ E[X_{.j.}] &= \mu + \beta_j \\ E[X_{...}] &= \mu \end{align}\)

因此,各参数的无偏估计量为:

\(\begin{align} \hat{\mu} &= X_{...} \\ \hat{\beta}_j &= X_{.j.} - X_{...} \\ \hat{\alpha}_i &= X_{i..} - X_{...} \\ \hat{\gamma}_{ij} &= X_{ij.} - \hat{\mu} - \hat{\beta}_j - \hat{\alpha}_i = X_{ij.} - X_{i..} - X_{.j.} + X_{...} \end{align}\)

为了为原假设 \(H_0^{int}\)\(H_0^r\)\(H_0^c\) 构建检验,我们首先注意到:

\(\sum_{k=1}^{l} \sum_{j=1}^{n} \sum_{i=1}^{m} \frac{(X_{ijk} - \mu - \alpha_i - \beta_j - \gamma_{ij})^2}{\sigma^2}\)

是一个自由度为 \(nml\) 的卡方分布随机变量。因此:

\(\sum_{k=1}^{l} \sum_{j=1}^{n} \sum_{i=1}^{m} \frac{(X_{ijk} - \hat{\mu} - \hat{\alpha}_i - \hat{\beta}_j - \hat{\gamma}_{ij})^2}{\sigma^2}\)

也将是卡方分布,但是其自由度将因为每估计一个参数而减少 1。因为 \(\sum_{i} \alpha_i = 0\),所以只需要估计 \(m-1\)\(\alpha_i\);同样地,由于 \(\sum_{j} \beta_j = 0\),因此只需要估计 \(n-1\)\(\beta_j\)。另外,由于 \(\sum_{i} \gamma_{ij} = \sum_{j} \gamma_{ij} = 0\),我们只需要估计 \((m-1)(n-1)\)\(\gamma_{ij}\)。因为还需要估计 \(\mu\),因此我们总共需要估计的参数个数为:

\(n - 1 + m - 1 + (n-1)(m-1) + 1 = nm\)

因为:

\(\hat{\mu} + \hat{\alpha}_i + \hat{\beta}_j + \hat{\gamma}_{ij} = X_{ij\cdot}\)

因此,如果我们令:

\(SS_e = \sum_{k=1}^{l} \sum_{j=1}^{n} \sum_{i=1}^{m} (X_{ijk} - X_{ij \cdot})^2\)

\(\frac{SS_e}{\sigma^2}\) 是自由度为 \(nm(l-1)\) 的卡方分布。因此:\(\frac{SS_e}{nm(l-1)}\)\(\sigma^2\) 的无偏估计量。

假设现在我们想要检验行和列之间没有 交互效应 的假设,即检验:

\(H_0^{int}: \gamma_{ij} = 0, \quad i=1,\dots,m, \quad j=1,\dots,n\)

如果 \(H_0^{int}\) 为真,则随机变量 \(X_{ij.}\) 是均值为 \(E[X_{ij.}] = \mu + \alpha_i + \beta_j\) 的正态分布随机变量。

此外,因为 \(X_{ij.}\)\(l\) 个方差为 \(\sigma^2\) 的正态分布随机变量的平均值,因此:

\(\text{Var}(X_{ij.}) = \frac{\sigma^2}{l}\)

因此,如果原假设 \(H_0^{int}\) 为真:

\(\sum_{j=1}^{n} \sum_{i=1}^{m} \frac{l(X_{ij.} - \mu - \alpha_i - \beta_j)^2}{\sigma^2}\)

是一个自由度为 \(nm\) 的卡方分布。因为总共 \(1 + n - 1 + m - 1 = n + m - 1\) 个参数 \(\mu, \alpha_i, \beta_j\) 必须估计,所以:

\(SS_{int} = \sum_{j=1}^{n} \sum_{i=1}^{m} l(X_{ij.} - \hat{\mu} - \hat{\alpha}_i - \hat{\beta}_j)^2 = \sum_{j=1}^{n} \sum_{i=1}^{m} l(X_{ij.} - X_{i..} - X_{.j.} - X_{...})^2\)

所以,在 \(H_0^{int}\) 为真的情况下:

\(\frac{SS_{int}}{\sigma^2}\)

是自由度为 \((n-1)(m-1)\) 的卡方分布。

因此,在没有行和列的 交互效应 的假设下:

\(\frac{SS_{int}}{(n-1)(m-1)}\)

\(\sigma^2\) 的无偏估计量。

可以证明,在没有 交互效应 的假设下,\(SS_e\)\(SS_{int}\) 是相互独立的,故当 \(H_0^{int}\) 为真时,

\(F_{int} = \frac{SS_{int}/[(n-1)(m-1)]}{SS_e/[nm(l-1)]}\)

是分子自由度为 \((n-1)(m-1)\)、分母自由度为 \(nm(l-1)\)\(F-\text{分布}\) 随机变量。在显著性水平为 \(\alpha\) 下,原假设 \(H_0^{int} : 所有 \gamma_{ij} = 0\) 的检验为:

\(\begin{align} \text{拒绝} \quad & H_0^{int} \quad \text{如果} \ \frac{SS_{int}/[(n-1)(m-1)]}{SS_e/[nm(l-1)]} > F_{\alpha, (n-1)(m-1), nm(l-1)} \\ \text{不拒绝} \quad & H_0^{int} \quad \text{否则} \end{align}\)

或者,我们可以计算 \(p-\text{value}\)。如果 \(F_{int} = v\),那么检验所有交互项等于 0 的原假设的 \(p-\text{value}\) 是:

\(p-\text{value} = P\{F_{(n-1)(m-1), nm(l-1)} > v\}\)

\(H_0^r\) 为真时,\(X_{i..}\)\(nl\) 个正态分布随机变量的均值,每个随机变量均值为 \(\mu\)、方差为 \(\sigma^2\)。因此,在 \(H_0^r\) 为真时:

\(E[X_{i..}] = \mu, \quad \text{Var}(X_{i..}) = \frac{\sigma^2}{nl}\)

因此,\(\sum_{i=1}^{m} \frac{nl(X_{i..} - \mu)^2}{\sigma^2}\) 是自由度为 \(m\) 的卡方分布。因此,如果我们令:

\(SS_r = \sum_{i=1}^{m} nl(X_{i..} - \hat{\mu})^2 = \sum_{i=1}^{m} nl(X_{i..} - X_{...})^2\)

则当 \(H_0^r\) 为真时,\(\frac{SS_r}{\sigma^2}\) 是自由度为 \(m-1\) 的卡方分布。因此:\(\frac{SS_r}{m-1}\) 是当 \(H_0^r\) 为真时,\(\sigma^2\) 的无偏估计量。

可以证明,在 \(H_0^r\) 为真时,\(SS_e\)\(SS_r\) 是独立的,故当 \(H_0^r\) 为真时:

\(\frac{SS_r/(m-1)}{SS_e/[nm(l-1)]}\)

是一个 \(F_{m-1, nm(l-1)}\) 分布随机变量。

因此,在显著性水平为 \(\alpha\) 时,

\(H_0^r : 所有 \alpha_i = 0 \quad vs \quad H_1^r : 至少一个 \alpha_i \neq 0\)

的检验为:

\(\begin{align} \text{拒绝} \quad & H_0^r \quad \text{如果} \quad \frac{SS_r/(m-1)}{SS_e/[nm(l-1)]} > F_{\alpha, m-1, nm(l-1)} \\ \text{不拒绝} \quad & H_0^r \quad \text{否则} \end{align}\)

如果令:\(\frac{SS_r/(m-1)}{SS_e/[nm(l-1)]} = v\),则 \(p-\text{value}\) 为:

\(p\text{-value} = P\{F_{m-1, nm(l-1)} > v\}\)

同理,也可以得到 \(H_0 : \beta_j = 0\) 的检验结果。表 10.4 对本节的 方差分析ANOVA)相关信息进行了总结。

表 10.4: 对于 \(m\)\(n\) 列的每个 \(X_{ij}\) 都有 \(l\) 个观察数据时的双因素方差分析,其中 \(N=nm(l-1)\)
类型 自由度 平方和 \(F-\text{统计量}\) 显著性水平 \(\alpha\) 下的检验 \(F = v\) 时的 \(p-\text{value}\)
行效应 \(m - 1\) \(SS_r = ln \sum_{i=1}^{m} (X_{i..} - X_{...})^2\) \(F_r = \frac{SS_r/(m-1)}{SS_e/N}\) 拒绝 \(H_0^r\) 如果 \(F_r > F_{\alpha, m-1, N}\) \(P\{F_{m-1, N} > v\}\)
列效应 \(n - 1\) \(SS_c = lm \sum_{j=1}^{n} (X_{.j.} - X_{...})^2\) \(F_c = \frac{SS_c/(n-1)}{SS_e/N}\) 拒绝 \(H_0^c\) 如果 \(F_c > F_{\alpha, n-1, N}\) \(P\{F_{n-1, N} > v\}\)
交互效应 \((n - 1)(m - 1)\) \(\begin{align} &SS_{int} = l \sum_{j=1}^{n} \\ & \times \sum_{i=1}^{m} (X_{ij.} - X_{i..} - X_{.j.} + X_{...})^2 \end{align}\) \(F_{int} = \frac{SS_{int}/[(n-1)(m-1)]}{SS_e/N}\) 拒绝 \(H_0^{int}\) 如果 \(F_{int} > F_{\alpha, (n-1)(m-1), N}\) \(P\{F_{(n-1)(m-1), N} > v\}\)
误差 \(N\) \(\begin{align} & SS_e = \sum_{k=1}^{l} \sum_{j=1}^{n} \\ & \times \sum_{i=1}^{m} (X_{ijk} - X_{ij.})^2 \end{align}\)

注意,对于前面介绍的所有的检验而言,只有在相关的 \(F-\text{统计量}\) 较大时,我们才会拒绝原假设。因为当 \(H_0\) 不成立时,\(F-\text{统计量}\) 的分子的分布往往会变大;而不论 \(H_0\) 是否成立,\(F-\text{统计量}\) 的分母的分布都是相同的。所以只有较大的(而不是较小的)\(F-\text{统计量}\) 的值才会让我们拒绝原假设。

习题 10.4 我们认为,某种特定类型的发电机的寿命会受到所使用的材料以及所使用场所的温度的影响。下表显示了由三种不同材料制造的发电机在两种不同温度下使用时的 24 台发电机的寿命数据。该数据是否可以表明,发电机的材料和工作场所的温度确实会影响发电机的寿命?是否有证据可以表明材料和工作场所的温度对发电机的寿命存在交互效应?

材料 10°C 18°C
1 135, 150
176, 85
50, 55
64, 38
2 150, 162
171, 120
76, 88
91, 57
3 138, 111
140, 106
68, 60
74, 51

解 10.4. 计算得出用于检验材料和温度之间没有交互效应的 \(F-\text{统计量}\) 的值为 0.64525,其 \(p-\text{value}\) 为:

\(p-\text{value} = P(F_{2,18} > 0.64525) > 1 - pf(0.64525, 2, 18)\)

代码
1 - pf(0.64525, 2, 18)
[1] 0.5362424

因此,我们无法拒绝材料和温度之间没有交互效应的原假设。

计算得出用于检验发电机寿命不依赖于所使用材料(即没有行效应)的 \(F-\text{统计量}\) 值为 2.47976,其 \(p-\text{value}\) 为:

\(p-\text{value} = P(F_{2,18} > 2.47976) > 1 - pf(2.47976, 2, 18)\)

代码
1 - pf(2.47976, 2, 18)
[1] 0.111889

因此,在 10% 的显著性水平下,我们无法拒绝发电机寿命不依赖于所使用材料的原假设。

计算表明,用于检验发电机寿命不依赖于温度(即没有列效应)的 \(F-\text{统计量}\) 值为 69.63223,其 \(p-\text{value}\) 为:

\(p\text{-value} = P(F_{1,18} > 69.63223) = 1 - pf(69.63223, 1, 18)\)

代码
1 - pf(69.63223, 1, 18)
[1] 1.337273e-07

因此,我们应该拒绝发电机寿命不依赖于温度的原假设。\(\blacksquare\)

Problems