统计学基础之单因子方差分析

方差分析

  前面讨论的都是一个总体或两个总体的统计分析问题，在实际工作中我们还会经常碰到多个总体均值的比较问题，处理这类问题通常采用所谓的方差分析方法。

问题的提出

  在饲料养鸡增肥的研究中,某研究所提出三种饲料配方：$A_1$是以鱼粉为主的饲料，$A_2$是以槐米粉为主的饲料，$A_3$是以苜蓿粉为主的饲料。为比较三种饲料的效果,特选24只相似的雏鸡随机均分为三组，每组各喂一种饲料，60天后观察它们的重量。试验结果如下所示:

饲料A	鸡重(单位:g)
$A_1$	1073 1009 1060 1001 1002 1012 1009 1028
$A_2$	1107 1092 990 1109 1090 1074 1122 1001
$A_3$	1093 1029 1080 1021 1022 1032 1029 1048

  本例中，我们要比较的是三种饲料对鸡的增肥作用是否相同。为此，把饲料称为因子，记为A，三种不同的配方称为因子A的三个水平，记为$A_1,A_2,A_3$，使用配方$A_i$下第j只鸡60天后的重量用$y_{ij}$表示，$i=1，2，3，j=1，2...8$。我们的目的是比较三种饲料配方下鸡的平均重量是否相等，为此，需要做一些基本假定，把所研究的问题归结为一个统计问题，然后用方差分析的方法进行分析。若相等，可任选一种饲料，特别可选廉价饲料。若不等，应选增肥效果好的饲料。

单因子方差的统计模型

  上述试验只考察了一个因子，称其为单因子试验。通常，在单因子试验中，记因子为A，设其有$r$个水平，记为$A_1，A_2，...，A_r$在每一水平下考察的指标可以看成一个总体，现有r个水平，故有r个总体，假定：

• 每一总体均为正态总体，记为$N(\mu_i，\sigma^2)，i=1，2，...，r$。
• 各总体的方差相同，记为$\sigma_1^2=\sigma_2^2=...=\sigma_r^2=\sigma^2$
• 从每一总体中抽取的样本是相互独立的，即所有的试验结果y，都相互独立。

  我们的目的是验证不同水平的均值是否有显著性差异，即进行假设检验：$$H_0:\mu_1=\mu_2=...=\mu_r=\mu$$$$H_1:均值不全相等$$

  为对假设进行验证需要进行抽样，现从每个水平下抽取$m$个样本，在每个水平$A_i$下的个体$y_{ij}$与水平均值$\mu_i$通常是有差距的，这是由随机误差引起的，记$$\epsilon_{ij}=y_{ij}-\mu_i$$

  则单因子方差的统计模型可表示为：$$y_{ij}=\epsilon_{ij}+\mu_i$$$$\epsilon_{ij}独立同分布于N(0,\sigma^2)$$

• 总均值$$\mu = \frac 1r\sum_{i=1}^{i=r}\mu_i$$
• 水平效应：又称第$i$水平的主效应$$a_i=\mu_i-\mu$$

  则统计模型可改写为：$$y_{ij} = \mu+a_i+\epsilon_{ij}$$$$\sum_{i=1}^{r}a_i=0$$$$\epsilon_{ij}独立同分布于N(0,\sigma^2)$$

模型假设可改写为：$$H_0:a_1=a_2=...=a_r$$

平方和分解

  在完成对总体的建模后，需要对样本数据进行处理。

组内偏差与组间偏差

• 组内偏差：每个组内的单个数据$y_{ij}$与缩在组的均值$\bar y_{i·}$的差$y_{ij}-\bar y_{i·}$
• 组间偏差：每个组的平均数据$\bar y_{i·}$与总平均值$\bar y$的差$\bar y_{i·}-\bar y$

偏差平方和及自由度

• 偏差平方和：一组数据中每个数据与均值之差的平方和$$Q=\sum_{i=1}^k(y_i-\bar y)^2$$
• 自由度：独立数据个数，由于偏差平方和由$k$个数据组成，其中有约束$$\sum(y_i-\bar y)=0$$

因此，独立偏差只有$k-1$个，故自由度为$$f_Q=k-1$$

总平方和分解公式

• 总偏差平方和$S_T$：$$S_T=\sum_{i=1}^r\sum_{j=1}^m(y_{ij}-\bar y)^2,f_T=n-1$$
• 组内偏差平方和：也称误差偏差平方和$$S_e=\sum_{i=1}^r\sum_{j=1}^m(y_{ij}-\bar y_{i·})^2，f_e=n-r$$
• 组间偏差平方和：也称因子A的偏差平方和$$S_A=\sum_{i=1}^r\sum_{j=1}^m(\bar y_{i·}-\bar y)^2，f_A=r-1$$
• 定理：在上述符号下，总偏差平方和可分解为组内和组间偏差平方和的的和，即$$S_T=S_A+S_e，f_T=f_A+f_e$$

证明：$$\begin{split}S_T&=\sum_{i=1}^r\sum_{j=1}^m(y_{ij}-\bar y)^2\\&=\sum_{i=1}^r\sum_{j=1}^m(y_{ij}-\bar y_{i·}+\bar y_{i·}-\bar y)^2\\&=S_A+S_e+\sum_{i=1}^r\sum_{j=1}^m2(y_{ij}-\bar y_{i·})(\bar y_{i·}-\bar y)\end{split}$$

由于$$\sum_{i=1}^r\sum_{j=1}^m(y_{ij}-\bar y_{i·})(\bar y_{i·}-\bar y)=\sum_{i=1}^r[(\bar y_{i·}-\bar y)\sum_{j=1}^m(y_{ij}-\bar y_{i·})]=0$$故原命题得证。

检验方法

• 均方：$$MS = \frac Q {f_Q}$$

• 定理：在单因子方差模型及前述符号下，有

  1. $\frac {S_e}{\sigma^2}\sim \chi^2(n-r)，E(S_e)=(n-r)\sigma^2$
  2. $E(S_A)=(r-1)\sigma^2+m\sum a_i^2$，若$H_0$成立，则$\frac{S_A}{\sigma^2}\sim\chi^2(r-1)$
  3. $S_A$与$S_e$独立

证明：

由正态性假设可知，$$\frac{\sum_{j=1}^{m}(y_{ij}-\bar y_{i·})^2}{\sigma^2}\sim \chi^2(m-1)$$

由卡方分布的可加性可知$$\frac{\sum_{i=1}^r\sum_{j=1}^{m}(y_{ij}-\bar y_{i·})^2}{\sigma^2}\sim \chi^2(rm-r)$$

即$$\frac{S_e}{\sigma^2}\sim\chi^2(n-r)$$定理1得证。

$$S_A=m\sum_{i=1}^r(\bar y_{i·}-\bar y)^2=m\sum_{i=1}^r(\mu_{i}+\bar\epsilon_{i·}-\mu-\bar \epsilon)^2=m\sum_{i=1}^r(a_{i}+\bar\epsilon_{i·}-\bar \epsilon)^2$$

  由于正态总体的样本均值和方差相互独立，因此，$S_e$与$\bar\epsilon_{i·}$独立，而$S_A$仅是$\bar\epsilon_{i·}$的函数，故定理3得证。

$$E(S_A)=m\sum_{i=1}^ra_i^2+E[m\sum_{i=1}^r(\bar\epsilon_{i·}-\bar \epsilon)^2]$$

由于$$\frac{\sum_{i=1}^rm(\bar\epsilon_{i·}-\bar \epsilon)^2}{\sigma^2}\sim \chi^2(r-1)$$

故$$E[m\sum_{i=1}^r(\bar\epsilon_{i·}-\bar \epsilon)^2]=(r-1)\sigma^2$$

定理2得证。

• 构造检验统计量

  由上述定理可知，$$F=\frac{MS_A}{MS_e}\sim F(r-1,n-r)$$

  F统计量表示的是组间差异与组内差异的比值，当原假设不满足时，组间差异较大，因此，当F统计量较大时拒绝原假设，即拒绝域$$\{F>F_{1-\alpha}(r-1,n-r)\}$$

使用python完成检验

from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
import pandas as pd


#转换数据形式
df = pd.DataFrame({'A1':(1073,1009,1060,1001,1002,1012,1009,1028),
                   'A2':(1107,1092,990,1109,1090,1074,1122,1001),
                   'A3':(1093,1029,1080,1021,1022,1032,1029,1048)})
df = pd.DataFrame([[i[0] for i in df.unstack().index],
                   list(df.unstack())],index=['group','values']).T
df['values']=pd.to_numeric(df['values'])
df.head()

	group	values
0	A1	1073
1	A1	1009
2	A1	1060
3	A1	1001
4	A1	1002

#输出结果
anova_lm(ols('values~C(group)',df).fit())

	df	sum_sq	mean_sq	F	PR(>F)
C(group)	2.0	9660.083333	4830.041667	3.594816	0.045432
Residual	21.0	28215.875000	1343.613095	NaN	NaN

  可以看出，p-value为0.045，因此，在0.05的显著性水平下，拒绝原假设，即认为三种鸡饲料对鸡的增肥有明显差异。

参数估计

  在完成显著性检验后，即可对参数进行估计，由于各水平为正态总体，因此，有如下估计：

• 各水平均值$$\hat\mu_i=\bar y_{i·}$$
• 总均值$$\hat\mu=\bar y$$
• 水平效应$$\hat a_i=\bar y_{i·}-\bar y$$
• 误差方差$$\hat{\sigma^2}=MS_e$$
• 置信区间：各水平均值$\mu_i$的置信区间$$\bar y_{i·}\sim N(\mu_i,\frac{\sigma^2}{m})$$$$\frac {\bar y_{i·}-\mu_i}{\sqrt{MS_e/m}}\sim t(f_e)$$

故置信区间为$$\bar y_{i·}\pm\sqrt{MS_e/m}t_{1-\alpha/2}(f_e)$$