数模系列(6)：方差分析（ANOVA）
1. 引子：我们不一样
1.1. 先来瞎扯扯
上期的k近邻算法中说到不同的特征对于结果的影响程度，如何取舍这些特征，可以考虑下方差分析。这一期我们就来看看方差分析（Analysis of Variance，简称ANOVA）。
我是学计算数学的，本科大二上概率论课程（教材是茆诗松版）的时候，最后一章就接触到了方差分析，还有其他的一些课也包含方差分析，虽然当时没怎么学，但后来可以自己参照些教材学。陈希孺先生的《概率论与数理统计》这本书很多人推荐，我在此也推荐一番，很多地方想要打破砂锅问到底，知其然也要知其所以然的，可以去看看这本统计书，下文也参考了这本书。
方差分析是英国大统计学家R. A. Fisher在1919年的英国的一个农业试验站工作时提出来的，需要进行很多的田间试验，于是他发明了方差分析法。其实Fisher 这个名字在统计学中经常见到，最近看的周志华的西瓜书里的Fisher 判别分析（其实Fisher判别分析里用到的数学知识和方差分析很相像）。对，猜得没错，这就是那个Fisher，不然别人怎么是大统计学家呢。有兴趣的可以查下wiki，有人对他的评价是，“a genius who almost single-handedly created the foundations for modern statistical science”，“一个几乎独自打造了现代统计科学基石的天才”，“the single most important figure in 20th century statistics”,“20 世纪统计学中最重要的人物”。既然这么重要的人物，怎么能错过照片。
1.2. 多态是好事
最近有首歌比较火，堪称神曲，《我们不一样》，虽然我对此歌并不感冒，但也被动地听了几遍。其中一句很火的歌词是：“我们不一样，每个人有不同的境遇”。
我从王小波的《思维的乐趣》里面看到，“我最赞成罗素先生的一句话：‘须知层次多态，乃是幸福的本源’”。我也是赞成的，每个人都一样的世界太无趣。
同样的，方差分析其实也和我想的一样，它也喜欢多看到一些不一样的东西，并同时分析出造成这种不一样的原因。
我们看一个实例，孙杨等8人一起游400米自由泳，当然最终的结果有快有慢，孙杨毫无悬念地夺得了第1 名。我们想知道是什么造就了他们8个人的成绩差异。于是我们找来了身高，体重，臂长，游泳馆里训练的时间，每天的饭量，喜欢的音乐等等特征变量，那么怎么知道这些指标对于衡量游泳成绩是否有效果呢？
利用方差分析可以得到这些特征变量（在方差分析中称特征变量为因素或者因子）对于最终的游泳成绩的影响程度是否是显著的。当然，我们做方差分析的时候，肯定喜欢看到对于某个因素，他们8个人之间的因素值差异巨大，且对成绩的影响程度是高度显著的。
这个实例是我自己编的，由于原始数据对最终结果的影响太大，瞎编数据的话，结果不管怎样，都可信度太差，所以下面的例子还是用书中的农业种植的实例。
2. 方差分析详解（结合农业种植实例）
在现实的生产和经营管理中，经常要分析各种因素对研究对象的某些试验结果（如游泳成绩，产品质量，数量，销量，成本等等）的影响。此时就要用到方差分析，其分为单因素方差分析，无交互作用的双因素方差分析，有交互作用的双因素方差分析，多因素正交表设计及方差分析。
我们来看一个农业种植的实例。某地之前没种过小麦，现在打算种植小麦，其他地方的一些优良品种，施肥方法，种植温度和种植压力（4个因素，每个因素下有不同的水平）可供此地种植小麦的参考备选。为此，进行田间试验，得到最好的种植因素的组合，希望能有个好收成。
2.1. 假设检验(Hypothesis Tests)
先岔开方差分析不谈，说下假设检验，因为接下来会用到。
2.1.1. 假设检验的本质
简单来说，假设检验就是一句话：小概率事件原理。
再多解释下，就是“在你的假设下，这事件发生的概率太低了，我根本就不相信！”或者是“在你的假设下，这事件发生的概率不算低，我没理由拒绝你的假设，只好姑且相信。”
2.1.2. 假设检验的基本原理
对总体参数进行假设检验时，首先要给定1个原假设 H_0 ， H_0 是关于总体参数的描述，与此同时存在1个与 H_0 相对立的备择假设 H_1 ， H_0 和 H_1 有且仅有1个成立。经过1次抽样后，若发生了小概率事件（一般取p<0.05），可依据“小概率事件在一次实验中几乎不可能发生”的理由，推翻原假设，即拒绝原假设 H_0 ，从而接受备择假设 H_1 ；反之，小概率事件并没有发生，就没有理由拒绝 H_0 。
2.1.3. 假设检验的一般步骤
（1）：根据问题确立原假设 H_0 和备择假设 H_1 ；
（2）：确定一个显著性水平 \alpha ，它是衡量小概率事件的标准，常取为0.05；
（3）：选定合适的检验用统计量W（通常原假设成立时，W的分布已知），由样本观测值计算出W的观测值 W_0 和衡量结果极端性的p值（原假设成立时得到样本观测值和更极端结果的概率 p=P\{|W|\ge|W_0|\} ）。
（4）：比较p值和显著性水平 \alpha ，若p值小于显著性水平，拒绝原假设 H_0 ，若不小于，则不拒绝原假设 H_0 。
三大抽样分布，卡方分布，F分布，和t分布均对应着各自的假设检验，在实际应用中广泛出现。
2.2. 单因素方差分析
2.2.1. 单因素方差分析理论
回到方差分析的话题上，我们先来看最基本的：单因素方差分析。顾名思义，只考虑一个因素A对于所关心的试验结果的影响。A 取k 个水平，在每个水平上作 n_i 次试验（其中 i=1,2,\cdots,k 。）试验过程中除A 外其它影响试验结果的因素都保持不变（只有随机因素存在）。我们的任务是从试验结果推断，因素A 对试验结果有无显著影响，即当A 取不同水平时试验结果有无显著差别。
以 Y_{ij} 记第i个水平的第j个观测值,希望由此对不同水平下总体的均值进行比较。列表如下
常用以下的统计模型表示： Y_{ij}=a_{0i}+e_{ij} ，其中 i=1,2,\cdots,k;j=1,2,\cdots,n_i ， a_{0i} 表示水平 A_i 下的试验结果的均值， e_{ij} 表示为随机误差，在方差分析中为了得到有效的检验，还常假定 e_{ij} 满足以下的独立同分布和方差齐性的条件：
（1）： e_{ij} 为相互独立的；
（2）： e_{ij} 都服从正态分布，且 e_{ij} 的均值都为0，方差都相同。
为了和后面的多因素的情况进行统一，这里我们把统计模型改写为 \\Y_{ij}=a_{0i}+e_{ij}=u+a_{i}+e_{ij}
其中 u=\frac{1}{k}*\sum_{i=1}^{k}a_{0i} ,为总均值， a_{i}=a_{0i}-u 表示因素A的第i个水平的附加效应（和平均值比的差距）。
因素A的各个水平的优劣取决于 a_i 的大小，如果 a_1=a_1=\cdots=a_k=0 ，则表示因素A对于试验结果Y毫无影响。这是我们就说因素A的效应不显著。在实际应用中，所谓“显著”,是指 a_i 之间的差异要大到“一定的程度”（其实是指与随机误差相比）。
我们把所要检验的假设写为： \\H_0=a_1=a_2=\cdots=a_k=0
为了检验上述假设，我们做出下面的分析，为什么实际上各个 Y_{ij} 的值会有差异？从统计模型上看，不外乎两个原因：一是各个水平可能有差异，例如 A_1 水平可能远远优于 A_2 ，这就使得 Y_{1j} 很可能大于 Y_{2j} ；二是随机误差的存在。
先找一个衡量全部 Y_{ij} 变异的量，它自然而然地取为了： \\SS=\sum_{i=1}^k\sum_{j=1}^{n_i}(Y_{ij}-\overline{Y})^2 ，其中 \overline{Y}=\frac{1}{n}\sum_{i=1}^k\sum_{j=1}^{n_i}Y_{ij} ， n=n_1+\cdots+n_k 。
SS又被称为总偏差平方和，SS越大，表示 Y_{ij} 之间的差异越大。根据上面的分析，设法把把SS分解为两部分，一部分表示随机误差的影响，记为 SS_e ，又被称为组内偏差；另一部分表示因素A的各水平之间的差异带来的影响，记为 SS_A ，又被称为组间偏差。
先考虑 SS_e 的部分，第i个水平的 n_i 次试验结果依次为 Y_{i1},Y_{i2},\cdots,Y_{in_{i}} 。它们之间的差异和因素A的各个水平完全无关，只和随机误差相关。
衡量第i个水平的 n_i 次试验结果的差异程度的量是 \sum_{j=1}^{n_i}(Y_{ij}-\overline{Y_i})^2 ，其中 \overline{Y_i}=\frac{1}{n}(Y_{i1}+Y_{i2}+\cdots+Y_{in_{i}}) ， i=1,2,\cdots,k 。 \overline{Y_i} 是第i个水平的试验结果的算术平均，将所有的k个水平的平方和相加，可得： \\SS_e=\sum_{i=1}^{k}\sum_{j=1}^{n_i}(Y_{ij}-\overline{Y_i})^2
我们将SS和 SS_e 相减看看会得到什么，应该就是 SS_A 吧。可以得到 \\SS_A=SS-SS_e\\=\sum_{i=1}^{k}n_i(\overline{Y_i}-\overline{Y})^2
从 SS_A 的表达式可以看出，它是衡量因素A的不同水平间的差异的平方和。
为了将 SS=SS_A+SS_e 看得更清晰，这里可以利用一个小技巧。 令 (Y_{ij}-\overline{Y})=(Y_{ij}-\overline{Y_i})+(\overline{Y_i}-\overline{Y}) 。 两边平方，先固定i对j进行求和，注意到： \sum_{j=1}^{n_i}(Y_{ij}-\overline{Y_i})(\overline{Y_i}-\overline{Y})\\=(\overline{Y_i}-\overline{Y})\sum_{j=1}^{n_i}(Y_{ij}-\overline{Y_i})=0 再对i进行求和即可得到 SS=SS_A+SS_e 。
引出了总偏差平方和SS，因素A的平方和 SS_A 和误差平方和 SS_e ， 我们就可以得知方差分析的名称由来。是因为这三个表达式都是属于样本方差一类的形状。
若记 \\MS_A=\frac{SS_A}{(k-1)},MS_e=\frac{SS_e}{(n-k)} ，其分别表示因素A 和随机误差的平均平方和，分母k-1和n-k分别为这两个平方和的自由度。基于前面的独立同正态分布和方差齐性的假设，当原假设 H_0 成立时，可以证明： \\\frac{MS_A}{MS_e}\sim F_{k-1,n-k}
此即对应于F分布的F检验，当 \frac{MS_A}{MS_e}\le F_{k-1,n-k}(\alpha) 时，其中 \alpha 为显著性水平，接受原假设 H_0 ，反之，拒绝原假设 H_0 。
一般来说，若得到的p值小于0.05,则表示因素A的效应显著，若小于0.01，则表示因素A的效应高度显著。
在统计应用上常把上述计算列出表格，记为方差分析表。
若最终得到的结果不显著，一般就不用考虑后续了。但若得到的结果显著，后续还需要进行多重比较，进行区间估计，这里便不展开了。
2.2.2. 单因素方差分析实例
在我们的例子中，因素A即为小麦种子的品种，总共3个品种（即水平）。分别重复4 次，5次和3次，试验结果即为亩产量，列出表格如下：
全部的算数平均为377.33，总平方和SS为 \\SS=(390-377.33)^2+\cdots+(408-377.33)^2=5274.67 它的自由度为12-1=11。
3个品种各自数据的算数平均分别为389.25，360.60和401.33，因此算出误差平方和为 SS_e=(390-389.25)^2+\cdots+(385-389.25)^2\\+(370-360.60)^2+\cdots+(362-360.60)^2\\+(413-401.3)^2+\cdots+(408-401.3)^2=1686.62
它的自由度为12-3=9。
品种平方和 SS_A 可由 SS_A=SS-SS_e 算出，为了验算，单独计算下：
SS_A=4*(389.25-377.33)^2+5*(360.60-377.33)^2\\+3*(401.33-377.33)^2=3588.05
，它的自由度为3-1=2。可以看见三者之间的关系完全没有错误，当然这是正常的。
由此可得到 \\MS_A=\frac{3588.05}{2}=1794.03,MS_e=\frac{1686.62}{9}=187.40 因素A 的F比为 \frac{MS_A}{MS_e}=9.57 ，计算 F_{2,9} 的9.57对应的p值为0.0059，不仅小于0.05，还小于0.01，因此可以看出种植品种这个因素是高度显著的。
以上这些均可通过matlab的anova1计算（具体使用查看help文档），得到的表格如下。
2.3. 双因素方差分析
如果要考虑两个因素 A,B对指标的影响，A,B各划分l，m个水平。对每一个水平的组合，作n次试验，（总共进行了lmn次试验）。 对所得数据进行方差分析，检验两个因素是否分别对指标有显著影响，或者还要进一步检验两因素是否对指标有显著的交互影响。
2.3.1. 无交互作用的双因素方差分析
我们把统计模型写为 \\Y_{ijk}=u+a_{i}+b_{j}+e_{ijk} 其中 i=1,2,\cdots,l;j=1,2,\cdots,m;k=1,2,\cdots,n ，u 为总均值， a_{i} 表示因素A的第i个水平的附加效应， b_{j} 表示因素B的第j个水平的附加效应。 e_{ijk} 表示在因素A的第i个水平和因素B 的第j个水平下的第k次试验的随机误差，这里仍需使它满足独立同均值为0的正态分布和方差齐性假设。
要说明因素A有无显著影响，就是要检验如下假设： \\H_{0A}=a_1=a_2=\cdots=a_l=0
要说明因素B有无显著影响，就是要检验如下假设： \\H_{0B}=b_1=b_2=\cdots=b_m=0
而模型无显著效果是指以上两个假设的原假设同时成立。
由于原理和单因素方差分析一致，自由度、均方和和F统计量已经在方差分析表中出现，这里我们只给出方差分析表如下。（根据单因素自己想想就可以推导，如果需要深究可详细参看参考文献）。
2.3.2. 有交互作用的双因素方差分析
若要考虑交互作用，则每个组合下的重复次数n必须大于1，不然没法计算交互作用的偏差平方和。
我们把统计模型写为 \\Y_{ijk}=u+a_{i}+b_{j}+ab_{ij}+e_{ijk} ，其中 i=1,2,\cdots,l;j=1,2,\cdots,m;k=1,2,\cdots,n ，u 为总均值， a_{i} 表示因素A的第i个水平的附加效应， b_{j} 表示因素B的第j个水平的附加效应, ab_{ij} 表示在因素A的第i个水平和因素B 的第j个水平下的交互作用的附加效应。 e_{ijk} 表示在因素A的第i个水平和因素B 的第j 个水平下的第k 次试验的随机误差，这里仍需使它满足独立同均值为0的正态分布和方差齐性假设。
要说明因素A有无显著影响，就是要检验如下假设： \\H_{0A}=a_1=a_2=\cdots=a_l=0
要说明因素B有无显著影响，就是要检验如下假设： \\H_{0B}=b_1=b_2=\cdots=b_m=0
要说明因素A和因素B的交互作用有无显著影响，就是要检验如下假设： \\H_{0AB}=ab_{ij}=0 其中 i=1,2,\cdots,l;j=1,2,\cdots,m 。
而模型无显著效果是指以上三个假设的原假设同时成立。
由于原理和单因素方差分析一致，自由度、均方和和F统计量已经在方差分析表中出现，这里我们只给出方差分析表如下。（根据单因素自己想想就可以推导，需要深究可详细参看参考文献）
2.4. 多因素正交表设计及方差分析
前面介绍了一个或两个因素的试验，由于因素较少，我们可以对不同因素的所有可能的水平组合做试验，这叫做全面试验。当因素较多时，虽然理论上仍可采用前面的方法进行全面试验后再做相应的方差分析，但是在实际中有时会遇到试验次数太多的问题。如4因素3水平的问题，所有不同水平的组合有81种，在每一种组合下只进行1次试验，也需做81次。如果考虑更多的因素及水平，则全面试验的次数可能会大得惊人。因此，在实际应用中，对于多因素做全面试验是不现实的。于是我们考虑是否可以选择其中一部分组合进行试验，这就要用到试验设计方法选择合理的试验方案，使得试验次数不多，但也能得到比较满意的结果。（其实就是想办法更有效率的偷懒，还不被发现！）
在这里，介绍一种叫做“正交表”的工具，它简单易用，在实际应用中广为流传。正交表的构造较为麻烦，但是在实用上，只需将已造出的统计著作中的正交表拿来直接用即可。
正交表的特点：
（1）每列中数字出现的次数相同;
（2）任取两列数字的搭配是均衡的，即任一列中同一数字的那些位置，在其他列中被该列所有不同的数字占据，且个数相同。
说了这么久正交表，然后现在一个表格都没看到。我的错，下面就祭出正交表 L_{8}(4*2^4) 。
（这个表的斜线表头LaTex折腾很久，use package diagbox 一直报错，暂时将就着看。）
这个正交表 L_{8}(4*2^4) 一共有8行，5列。 4*2^4 把它写开更好理解些 4*2*2*2*2 ，即表示：表中第1列代表的因素有4个水平，后面4列代表的因素各有2个水平。最后一列表示对应的各个因素的水平的试验结果。
检查是否满足正交表的两个特点：
（1）每列中数字出现的次数相同。例如第1列含不同数字1，2，3，4，每个数字出现2 次，其他4列也是。
（2）任取两列数字的搭配是均衡的，即任一列中同一数字的那些位置，在其他列中被该列所有不同的数字占据，且个数相同。例如，第3列中数字1 占据1，3，6，8行的位置,而在第1列中，这4个位置恰被1，2，3，4各占据了1 次。而在第5列，这4个位置则被该列不同数字1，2各占据2次，同样的也是第2 列，第3 列，达到均衡。
对各个因素同样的构造方差分析表，查看各个因素的效应是否达到显著。
这里给出的单因素方差分析，无交互作用的双因素方差分析，有交互作用的双因素方差分析，多因素正交表设计及方差分析，讲得内容并不深。想要知道更多的详细，建议参看下面的参考文献，或者直接搜寻对应的书籍和文献钻研。
3. 参考文献
1.陈希孺，《概率论与数理统计》，第6章第5节：方差分析。
2.司守奎，《数学建模算法与程序》，第11章：方差分析。
4. 无关内容
这期的最后，把王小波的两句话贴上来分享下，这两句话我很喜欢。
“我对自己的要求很低：我活在世上，无非想要明白些道理，遇见些有趣的事。倘能如我所愿，我的一生就算成功。”
“一个人只拥有此生此世是不够的，他还应该有一个诗意的世界。”
开头说到判别分析和方差分析里用到的数学知识很相像，下期更新的主角就是判别分析吧！
声明：转载需联系本人。
