5.1 引言

在工业生产或者实验中经常要比较若干个因素对业务指标的影响。比如，为了比较药物A,B,C对治疗某疾病的疗效，我们将实验对象分成三组，分别记录服用三种药物的治疗效果，得到三组样本

\[X_1,\dots,X_{n_1};\ Y_1,\dots,Y_{n_2};\ Z_1,\dots,Z_{n_3}.\]

通过这些实验数据，我们希望回答：

这个例子中，药物称为 因子 ，A,B,C称为该因子的 水平 。

由于这个实验只涉及单个因子——“药物”，我们称之为 单因子实验 。此外，如果比较不同的药物和性别对疗效的影响，这就是 两因子实验 。不难推广到 多因子实验 。

本章介绍 方差分析方法 (Analysis of Variance, ANOVA)，研究因子不同水平的差异性，不同因子交互作用的显著性。这些研究有助于搭配有利于指标的不同因子的水平。虽然本章叫做方差分析，但实际上关心的是不同总体的均值比较，而不是它们的方差。

5.2 单因子方差分析

单因子实验设计(one-way layout)在一个因子的不同水平下分别进行独立的观测，是两个独立样本比较方法的推广。

模型假设 ：考虑一个因子A，有 \(r\) 个水平 \(A_1,\dots,A_r\) , \(r\ge 2\) . 设在水平 \(A_i\) 下重复进行了 \(n_i\) 次实验( \(n_i\ge2\) )，数据是 \(y_{i1},y_{i2},\dots,y_{in_i}\) . 假设这些数据之间相互独立且 \(y_{ij}\sim N(\mu_i,\sigma^2)\) ，其中 \(\sigma\) 未知。我们关心的问题是 \(\mu_i\) 是否全相等，即要检验

\[H_0:\mu_1=\dots=\mu_r\ vs.\ H_1:\mu_i\text{不全相等}.\]

令 \(n=\sum_{i=1}^r n_i\) ，

\[\bar y = \frac{1}{n}\sum_{i=1}^r\sum_{j=1}^{n_i}y_{ij},\ \bar y_{i\cdot}=\frac 1{n_i}\sum_{j=1}^{n_i}y_{ij},i=1,\dots,r.\]

这里 \(\bar y\) 表示所有观测的平均值， \(\bar y_{i\cdot}\) 表示水平 \(A_i\) 下的观测平均值。令

\[S_T^2=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y )^2,\ S_e^2=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot} )^2,\]

\[S_A^2 = \sum_{i=1}^r n_i(\bar y_{i\cdot} -\bar y)^2.\]

其中， \(S_T^2\) 刻画全部数据的波动程度， \(S_e^2\) 刻画组内数据的波动程度， \(S_A^2\) 刻画不同组均值的差异引起的波动程度。这三者满足以下三角分解：

\[S_T^2 = S_e^2+S_A^2.\]

\[\begin{align*} S_T^2&=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot}+\bar y_{i\cdot}-\bar y )^2\\ &=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot})^2+\sum_{i=1}^r\sum_{j=1}^{n_i}(\bar y_{i\cdot}-\bar y )^2+2\sum_{i=1}^r(\bar y_{i\cdot}-\bar y )\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot})\\ &=S_e^2+S_A^2. \end{align*}\]

这表明，总的平方和等于组内平方和加上组间平方和。注意到 \(\bar y_{i\cdot}\) 是 \(\mu_i\) 的无偏估计，如果 \(H_0\) 成立， \(\bar y_{i\cdot}\) 应该接近，那么 \(S_A^2\) 相对 \(S_e^2\) 小得多。也就意味着两者的比值 \(S_A^2/S_e^2\) 大到一定程度就有理由拒绝 \(H_0\) . 为给出确切的拒绝域， 我们需要知道在 \(H_0\) 成立下， \((S_A^2, S_e^2)\) 的分布。

定理 5.1 考虑上述模型假设，有 \(S_e^2/\sigma^2\sim \chi^2(n-r)\) 且 \(S_e^2\) 与 \(S_A^2\) 独立。如果 \(H_0:\mu_1=\dots=\mu_r\) 成立，则有

\[\frac{S_A^2}{\sigma^2}\sim \chi^2(r-1),\ F=\frac{S_A^2/(r-1)}{S_e^2/(n-r)}\sim F(r-1,n-r).\]

证明. 由单个正态总体的抽样分布定理有，

\[V_i:=\frac{1}{\sigma^2} \sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot} )^2\sim \chi^2(n_i-1).\]

由于 \(y_{ij}\) 之间独立， 所以 \(V_i\) 相互独立。由卡方分布的可加性，我们有 \(S_e^2/\sigma^2=\sum_{i=1}^r V_i\sim \chi^2(n-r)\) . 此外， \(\{V_1,\dots,V_r\}\) 与 \((\bar y_{1\cdot},\dots,\bar y_{r\cdot})\) 独立。注意到 \(S_A^2\) 为 \(\bar y_{i\cdot}\) 的函数，所以 \(S_e^2\) 与 \(S_A^2\) 独立。 当 \(\mu_1=\dots=\mu_r=\mu\) 成立时， \(\bar y_{i\cdot}\stackrel{iid}{\sim} N(\mu,\sigma^2/n_i)\) . 令 \(x_i = \sqrt{n_i}(\bar y_{i\cdot}-\mu)/\sigma\stackrel{iid}\sim N(0,1)\) . 注意到，

\[\begin{align*} S_A^2&=\frac{\sum_{i=1}^rn_i(\bar y_{i\cdot} -\bar y)^2}{\sigma^2}\\ &=\frac{\sum_{i=1}^r(\sqrt{n_i}\bar y_{i\cdot} -\sqrt{n_i}\sum_{j=1}^r \frac{n_j}{n}\bar y_{j\cdot} )^2}{\sigma^2}\\ &=\frac{\sum_{i=1}^r[\sqrt{n_i}(\bar y_{i\cdot}-\mu) -\sqrt{n_i}\sum_{j=1}^r \frac{n_j}{n}(\bar y_{j\cdot}-\mu) ]^2}{\sigma^2}\\ &=\sum_{i=1}^r(x_i-\sqrt{n_i}\sum_{j=1}^r \frac{\sqrt{n_j}}{n}x_j)^2\\ &=\sum_{i=1}^rx_i^2-(\sum_{i=1}^r\sqrt{n_i/n} x_i)^2\\ &=||x_{1{:}r}||^2-(\alpha^\top x_{1{:}r})^2, \end{align*}\] 其中 \(\alpha=(\sqrt{n_1/n},\dots,\sqrt{n_r/n})^\top\) . 注意到 \(||\alpha||=1\) ，参考定理1.4的证明，可以构造一个正交矩阵 \(A\) 使得 \(A\) 的第一行为 \(\alpha^\top\) 。令 \(z_{1{:}r}=Ax_{1{:}r}\sim N(0,I_r)\) ，此时， \(||z_{1{:}r}||^2=||x_{1{:}r}||^2\) , \(z_1=\alpha^\top x_{1{:}r}\) . 所以，

\[S_A^2=||z_{1{:}r}||^2-z_1^2=\sum_{i=2}^r z_i^2\sim \chi^2(r-1).\]

为此，我们采用F检验，拒绝域为 \(W=\{F>F_{1-\alpha}(r-1,n-r)\}\) . 这种方法为方差分析法，是R. A. Fisher在1923年提出来的。在实践中常用以下方差分析表格。

library(ggplot2)
mouse <- data.frame(
  X = c(2, 4, 3, 2, 4, 7, 7, 2, 2, 5, 4, 
        5, 6, 8, 5, 10, 7, 12, 12, 6, 6, 
        7, 11, 6, 6, 7, 9, 5, 5, 10, 6, 3, 10),
  A = factor(rep(1:3,c(11,10,12)))
ggplot(mouse,aes(A,X,fill=A)) + geom_boxplot()

mouse.aov <- aov(X~A, data=mouse)
summary(mouse.aov)

##             Df Sum Sq Mean Sq F value Pr(>F)   
## A            2   94.3    47.1    8.48 0.0012 **
## Residuals   30  166.7     5.6                  
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## A            2    353   176.3    8.96 0.00049 ***
## B            3     88    29.2    1.48 0.23108    
## A:B          6     72    12.0    0.61 0.72289    
## Residuals   48    944    19.7                    
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1