Play with R 第10期：方差分析
第十章
10.1这章会讲什么？
10.2 ANOVA背后理论
10.2.1 膨胀的错误率：为什么不能直接用t检验比较三组间的差异？
10.2.2 关于F值的解释
10.2.3 作为回归的ANOVA
10.2.4 F比率（F-ratio）的逻辑
10.2.5 总平方和SST
10.2.6 模型平方和SSM
10.2.7 残差平方和SSR
10.2.8 均方根
10.2.9 F比率
10.3 ANOVA的可靠性
10.4 ANOVA的计划对比（其实是线性回归的思路做成组比较）
10.4.1 如何选择比较
10.4.2使用加权方法定义对比
10.4.3非正交的比较
10.4.4 标准比较
10.4.5 多项式比较：趋势分析
10.5 事后分析（post-hoc）
10.6. One-way ANOVA using R 使用R进行单因素方差分析
10.6.1. Packages for one-way ANOVA in R 单因素方差分析包
10.6.2. General procedure for one-way ANOVA 单因素方差分析的一般程序
10.6.3. Entering data 输入数据
10.6.4. One-way ANOVA using R Commander 使用R指令进行单因素方差分析
10.6.5. Exploring the data 探索数据
10.6.6. The main analysis 主要分析
10.6.6.1. When the test assumptions aremet 当满足测试假设时
10.6.6.2. When variances are not equal acrossgroups 当组间的方差不相等时
10.6.6.3. Robust ANOVA – it’s not for the weakof heart 方差分析鲁棒性—不适合脆弱的小心脏
10.6.7. Planned contrasts using R 使用R进行有计划的对比
10.6.7.1.Your own contrasts
10.6.7.2.Trend analysis
10.6.8.Post hoc tests using R 使用R进行事后检验
10.6.8.1. Bonferroni and related methods Bonferroni和相关方法
10.6.8.2. Tukey and Dunnett
10.7 计算效应量
10.8报告单因素方差分析的结果
10.1 这章会讲什么？
在人生的关键节点上，我们会考虑作出某个决定对未来的影响，当只有两种选择的时候（第九章），比如清华还是北大，我们可以通过t检验来比较两种决定的结果，然而，如果这个时候多了种选择呢？比如MIT也给了你offer，这时就需要用方差分析(analysis of variance,ANOVA)来比较这三种选择。
10.2 ANOVA背后理论
10.2.1膨胀的错误率：为什么不能直接用t检验比较三组间的差异？
如果用t检验比较1，2，3三组的差异那我们就需要检验：1和2的组间差异；2和3的组间差异；1和3的组间差异。假设我们每次t检验都采用0.05的显著性水平，那么每次检验错误地拒绝零假设（I类错误）的可能性是5%，那么不存在I类错误的可能性是95%，假设每次检验都是独立事件，那么三次检验均不存在I类错误的可能性是0.953=0.857，那么至少出现一次I类错误的可能性为1-0.857=0.143，即14.3%，而这已经高于只做一次t检验的I类错误率（5%），这一误差称之为族系（familywise）或者实验误差（experimentwise error）。将做3次t检验的族系误差推广到做n次t检验的族系误差，可得到如下计算公式（显著性水平为0.05）：
族系误差=1-(0.95)n
那么，随着检验次数的增多，族系误差越来愈大。
10.2.2关于F值的解释
如果我们做t检验，会得到t值，t值代表两组的均值是否存在显著性差异。相应地，我们做ANOVA，则会得到F值。而F值则代表多组的均值之间是否存在显著性差异，但F值代表了整体的一种效应，以三组比较为例，当F值显示应该拒绝零假设（即每组的均值是相等的）时，我们只知道至少这三组中的两组均值之间存在差异，但并不代表我们知道这种差异具体是什么样的。
10.2.3作为回归的ANOVA
第七章中提到了，F比率是评估回归方程在多大程度上能预测实际结果的指标，而当回归方程中的预测变量只有分类变量（如性别，实验操作）时，我们就可以把t检验看成是线性回归方程（详见9.4.2），ANOVA可以被看成是多重回归方程，且预测变量比自变量的类别少1（即自由度）。
举例如下，假设我们想要检验Viagra这种药是否能增强力比多（libido），我们在三组被试中进行研究：安慰剂组（吃糖片），低剂量组（吃低剂量的Viagra），高剂量组（吃高剂量的Viagra），并测量被试在一周时间内的力比多水平（因变量），结果如表10.1。
想要预测力比多的水平，那么我们应当采用如下方程：
如果我们想要采用线性模型，那么在只有两组的情况下，我们可以采用一个虚无变量（dummy variable）来描述两组，这样一组编码为0，另一组编码为1（详见9.4.2）。而在有三组的情况下，我们需要采用两个虚无变量b2（称之为高）和b1（称之为低）来进行编码，此时回归方程如下：
在知道了截距（b0）和某个被试的分组，我们可以通过将被试的分组编码代入这个回归方程，进而预测其力比多水平。
三组实验条件的虚无编码如表10.2所示，其中，安慰剂组为对照组（reference）或者基线（base）。
在代入虚无编码后，安慰剂组的回归模型如下，安慰剂组的力比多均值为b0。
在代入虚无编码后，高剂量组的回归模型如下：
因为安慰剂组的力比多均值为b0，所以我们可以得到：
所以，b2为高剂量组和安慰剂组力比多均值之差。
在代入虚无编码后，低剂量组的回归模型如下：
同样的，我们可以得到：
即b1为低剂量组和安慰剂组力比多均值之差。
结果如下：
结合表10.1里的数据来看，第一个虚无编码（b2）的回归系数为高剂量组与安慰剂组的力比多平均值差（5.0−2.2 = 2.8），第二个虚无编码的回归系数是第一个虚无编码（b1）的回归系数为高剂量组与安慰剂组的力比多平均值差（3.2 − 2.2 = 1），而根据t检验结果，高剂量组和安慰剂的差异（第一个虚无编码b2）达到了显著性水平（p < 0.05），低剂量组和安慰剂的差异（第二个虚无编码b2）未达到显著性水平（p = 0.282）.
将此种虚无编码方法推广到四个实验组中时，编码方法如表10.3所示。
10.2.4 F比率（F-ratio）的逻辑
在对于Viagra的效果研究中，我们可以通过图10.2来直观展示组间差异，其中横轴代表被试，纵轴代表力比多水平，每个数据点代表一个被试的力比多水平（圆点为安慰剂组，三角代表低剂量组，方块代表高剂量组），中间的黑线代表总体平均值，三组内的三条横线分别代表三组的均值，低剂量组部分的竖向双箭头b1代表低剂量组和安慰剂组力比多均值之差，高剂量组部分的竖向双箭头b2为高剂量组和安慰剂组力比多均值之差（注意：图中标识有错误，正方形点区域的‘b1’应该为‘b2’）。
当我们用ANOVA检验回归模型时，我们会将总平均值（grand mean）作为模型来预测实际数据，但是在我们的实际数据和我们的预测模型（总平均值）之间会存在差异，即模型无法解释的变异（方差），而F比率即为模型可以解释的变异（方差）除以模型不能解释的变异（方差）而得到的比率。
10.2.5 总平方和SST
平方和（SST）：首先计算每个观测到的数据点与均值之差，然后对这些差异求平方，并将它们相加，得出平方和。
方差与平方和是息息相关的，方差的计算公式是s2= SS /（N-1），其中N是观测值的数量。根据此公式，平方和的计算为SS = s2（N-1）。表10.1给出了Viagra数据，它的SST计算如下：
图10.3以图形的方式描述了不同平方和的表示方式，其中SST指的是观测数据与平均值Y之间的差异；SSR指的是观测数据和组均值之间的差异；SSM指的是平均值Y与组均值之间的差异。
图10.3ANOVA中不同平方和的图形表示
10.2.6 模型平方和SSM
SSM计算的步骤如下：
1、计算组均值和总均值之间的差
2、将这些差异平方
3、将每个结果乘以该组内的参与者数（nk）
4、将每个组的值加在一起。
此过程的数学表达式为：
因此，Viagra数据的SSM的计算如下：
SSR指的是观测数据和组均值之间的差异。因此，它的计算公式是：
每组的平方和表示该组中每个参与者的得分与组均值之间的平方差之和。因此，我们可以将SSR表示为SSR = SSgroup1+ SSgroup2 + SSgroup3 +…。依据上面方差和平方和之间的关系，我们就可以使用Viagra数据中每个组的方差来创建一个方程。因此，SSR可以表示为：
Viagra数据SSR的计算如下：
10.2.8 均方根
SSM和SSR都是求和值，因此它们将受到求和的数量的影响，比如，SSR和SST分别使用了12和15个数量，而SSM仅使用了3个数量（组均值）。因此，为了减少这种由于数值不同带来的偏差，我们可以计算平方和的平均值（即均方根，MS），计算方法是用平方和除以自由度。Viagra数据的均方根为：
10.2.9 F比率
F比率指的是模型解释的变异与非系统因素解释的变异之比。换句话说，它是模型的好与模型的坏（有多少错误）的比率。它的表达式为：
Viagra数据的F比率为：
10.3 ANOVA的可靠性
F统计检验是基于正态分布。也就是说，需要保证每种实验条件下的方差都必须非常相似（方差齐性）。我们可以采用Levene检验来看方差是否齐性，如果Levene检验显著（即p小于0.05），则说明方差不齐。
关于方差分析的可靠性。Glass等（1972）发现F检验在峰度和偏度不太好、或非正态条件下依然可以很好地控制I型错误率。另外，当各组样本量相同时，ANOVA特别可靠。但是，当各组样本量不等同时，ANOVA的可靠性就会降低，因此建议大家尽量在各组中采集相同的数据样本。
10.4 ANOVA的计划对比分析（其实是线性回归的思路做成组比较）
F比率只告诉我们，与无关的因素相比，适合数据的模型是否有更多的变化，但它不能告诉我们组间的差异在哪里。因此，如果F-ratio足够大，具有统计学意义，那么我们只知道均值之间的一个或多个差异具有统计学意义(例如，b2或b1具有统计学意义)。但我们不知具体的差异体现在哪些对比中。在多元回归中，每个b系数分别用t检验进行检验，我们也可以对方差分析进行同样的检验。然而，我们需要进行两次t检验，这将提高family wise的错误率(参见10.2节)。因此，我们需要一种方法来对比不同的组，而不增加类型I的错误率。有两种方法来实现这一目标:第一种方法是将模型所考虑的方差分解为组成部分:第二个是比较每组(如果进行多个t),但是使用一个更严格的矫正标准,familywise错误率不超过0.5。第一个选项可以使用计划的比较(也称为计划的对比)来完成，而第二个选项是使用事后比较来完成的(参见下一节)。计划比较和事后测试之间的区别可以比作单尾测试和双尾测试之间的区别，因为计划比较是在您想要测试特定假设时进行的，而事后测试是在您没有特定假设时进行的。
10.4.1 如何选择比较
在伟哥的例子中，我们可以有非常具体的假设。首先，与安慰剂组相比，我们认为任何剂量的伟哥都会改变性欲。作为第二个假设，我们可能认为高剂量比低剂量更能增加性欲。为了进行有计划的比较，必须在收集数据之前推导出这些假设。将实验条件与对照条件作为第一个对比，然后观察实验组之间的差异，这在科学上是相当标准的。方差分析的基础是将总变异分为两部分:由实验操作引起的变异(SSM)和由非系统因素引起的变异(SSR)，如下图所示。
计划中的比较通过将实验引起的变化分解为组成部分(参见图10.5)，进一步采用了这种逻辑。准确的比较取决于你想要测试的假设。图10.5显示了实验方差分解的情况看多少变异是由两个药物条件相比安慰剂条件(对比1)。然后服用伟哥的变异解释分解多少是由服用高剂量与低剂量(对比2)
图10.5 进行计划比较所假设的contrast
通常情况下，学生们会纠结于计划中的比较，但有三条规则可以帮助你想出该怎么做:
1 .如果我们有一个对照组，这通常是因为我们想把它与其他组进行比较。
2.每个对比只能比较两个变化的模块。
3.一旦一个群体在一个对比中被挑选出来，它就不能被用在另一个对比中
在10.6部分对如何根据这三个原则进行选择做了具体的解释，我们可以根据10.6部分的内容来进行计划比较的具体操作。
10.4.2使用加权方法定义对比
在这部分，我们需要学会如何根据回归中的思想来为想要进行的比较进行虚拟变量的权重赋值。这个过程看起来非常混乱，但是有一些基本的规则可以为虚拟变量赋值，以获得您想要的比较。作者解释了这些简单的规则，然后展示这个过程实际上是如何工作的。当你通读这些规则时，记住前面的部分，并提醒自己作者所说的变化是什么意思。
规则1:选择合理的比较。记住，你只想比较两个组的差异，如果一个组在一个比较中被挑选出来，那么这个组应该被排除在任何后续的对比之外。
规则2:编码为正权值的组将与编码为负权值的组进行比较。所以，分配一组正权值，另一组负权值。
规则3:比较的权值之和应该为零。如果将给定对比度的权重相加，结果应该是零。
规则4:如果一个组没有参与比较，则自动为它分配0的权重。如果我们给一个组的权值为0那么这个组就会被所有比较中剔除。
规则5:对于给定的对比，在一组变量中分配给组的权重应该等于在另一组变量中分配给组的权重。
下面我们具体来看：
使用规则5，我们分配给组块1中的组的权重应该等于组块2中的组的数量。在组块2中只有一个组，所以我们给组块1中的每个组分配一个1的权重。同样，我们给组块2中的组分配一个与组块1中的组数相同的权重。第一组有两组，所以我们给安慰剂组的权重是2。然后我们将权重的符号与大小结合起来，得出-2(安慰剂)、1(低剂量)和1(高剂量)的权重。如下图所示：
规则3规定，对于给定的对比，权重相加应该为0，通过遵循规则2和规则5，这个规则将始终被遵守(如果您没有适当地遵循这些规则，那么当您添加权重时，这一点将变得清晰)。那么，让我们通过添加权值来检查:权值的和= 1 + 1 -2 = 0。
第二个对比是比较两个实验组，所以我们要忽略安慰剂组。规则4告诉我们，我们应该自动为该组分配一个0的权重(因为这将从任何计算中消除该组)。我们只剩下两组变量:组块1包含低剂量组，组块2包含高剂量组。按照规则2和规则5，很明显一个组分配一个值：1和-1，如下图所示：
最后，我们用一个表格来总结赋值分配：
在这些比较中，重要的是将比较和的权值设置为0，因为这样可以确保比较的是两个不同的变量块。因此，我们可以进行t检验。一个更重要的考虑是，当您为一个特定的组乘以权重时，这些产品加起来应该是零(参见表10.4的最后一列)。如果乘积的和为零，那么我们可以确定对比是独立的或正交的。这非常重要。计算公式的内容其实是用数学方法展示了我们这一过程，这一部分在10.6代码实现过程中可以详细见到。重要的是理解设置contrast的权重的原理和过程。
在这些比较中，重要的是将比较和的权值设置为0，因为这样可以确保比较的是两个不同的变量块。因此，我们可以进行t检验。一个更重要的考虑是，当您为一个特定的组乘以权重时，这些产品加起来应该是零(参见表10.4的最后一列)。如果乘积的和为零，那么我们可以确定对比是独立的或正交的。这非常重要。计算公式的内容其实是用数学方法展示了我们这一过程，这一部分在10.6代码实现过程中可以详细见到。重要的是理解设置contrast的权重的原理和过程。
10.4.3非正交的比较
我花了很多时间来努力设计适当的正交比较，却没有提到非正交对比提供的可能性。非正交对比是在某种程度上相关的比较，获得它们的最佳方法是违反前一节中的规则1。用蛋糕类比来看，非正交比较是你把蛋糕切成片然后再把它们粘在一起。因此，对于伟哥的数据，一组非正交对比可能具有相同的初始对比(将实验组与安慰剂进行比较)，但随后将高剂量组与安慰剂进行比较。这违反了规则1，因为安慰剂组在第一个对比中被挑选出来，但在第二个对比中再次使用。表10.5显示了这组对比的编码，通过查看最后一列可以清楚地看到，当您将这两个对比的编码相乘并相加时，总和不是零。这告诉我们，对比不是正交的。
表10.5 非正交的对比矩阵
作者认为，执行非正交对比并没有本质上的错误。然而，如果你选择执行这种类型的对比，你必须非常小心如何解释结果。对于非正交对比，你所做的比较是相关的，因此结果的测试统计数据和p值将在某种程度上相关。因此，你应该使用一个更保守的概率水平来接受一个给定的对比是有统计意义的。
10.4.4 标准比较
虽然在大多数情况下，你会设计自己的对比，但有些特殊的对比是为了比较某些情况而设计的。这些对比中的一些是正交的，而另一些是非正交的。表10.6显示了在R中使用contrast()函数可获得的对比度。这个函数用于编码任何分类变量，得到的编码几乎可以用于任何线性模型(方差分析、回归、逻辑回归等)。虽然表10.6中没有提供准确的编码，但给出了在三组和四组情况下进行比较的例子(分别将这些组标记为1、2、3和1、2、3、4)。当您对变量进行编码时，R将把值最低的代码视为组1，将值次高的代码视为组2，依此类推。因此，根据您希望进行的比较，您应该对分组变量进行适当的编码(然后使用表10.6作为指南，说明R将执行哪些比较)。
表10.6 标准的比较所设置contrast的例子（三组和四组情况时）
10.4.5 多项式比较：趋势分析
表10.6中故意省略的一种对比是多项式对比，可以使用contr.poly()获得。这种对比测试了数据中的趋势，以其最基本的形式，它寻找的是线性趋势（例如，这个组意味着按比例增加)。然而，也有其他的趋势，如二次、三次和四次趋势可以检查。图10.8显示了数据集中可能存在的趋势类型的示例。线性趋势现在应该对你们所有人都很熟悉了，它表示因变量在有序类别中的值的简单比例变化(图中显示了一个正的线性趋势，但当然也可能是负的)。
二次趋势是指直线的方向有一个变化(例如，直线在一个地方被弯曲)。这方面的一个例子可能是一种情况，一种药物起初提高了一项任务的表现，但随着剂量的增加，表现再次下降。要找到一个二次趋势，你至少需要三组(因为在两组的情况下，自变量的类别数量不够，因变量的均值无法变化)。三次趋势是指趋势的方向有两个变化。举个例子，因变量的均值一开始在自变量的前几类中上升，然后在之后的几类中下降，然后在最后几类中上升。要使均值方向有两次改变至少要有四类自变量。最后一个可能遇到的趋势是四次趋势，这个趋势有三个方向的变化(所以至少需要五类自变量)。
多项式趋势应该在对自变量的类别排序有意义的数据集中进行检查(例如，如果你给了一种药物五次剂量，按大小顺序检查五次剂量是有意义的)。对于伟哥的数据，只有三组，所以我们只能期望找到一个线性或二次趋势(测试任何高阶趋势是没有意义的)。这些趋势中的每一个都有一组用于回归模型中的虚拟变量的代码，所以我们所做的事情与我们对计划的对比所做的事情是一样的，只是我们已经设计了代码来表示感兴趣的趋势类型。
图10.8 不同的趋势图（线性、二次、三次、四次）
10.5 事后分析（post-hoc）
我们知道当方差分析F检验显著时需要应用多重比较的方法进行检验。这里介绍三个方法Bonferronicorrection，Holm，Benjamini andHochberg。在介绍这三种方法之前先讲一下假设检验中的两类错误Type Ierror和Type II error ，当虚无假设为真却被拒绝时所犯的错误为Type I error，这类错误的概率为alpha，当备择假设为真却被拒绝时所犯的错误为Type II error，这类错误的概率为beta。1-beta为统计检验力。当alpha增大时1-beta减少，当alha减少时1-beta增大。也就是说对Type I error越保守，越容易错过发现数据间真正有差异的机会。
Bonferroni correction：具体为pcrit k（k为比较的次数）。以下这个图是解释了三种比较方法，其中有四组要进行比较Spiderman,Superman, the Hulk, and Teenage Mutant Ninja Turtles (NT)，一共需要比6次，/k.05/6.0083。Holm的方法是先将每个实际得出的p值从小到大排列，然后分别赋予j（6至1）的数值，对应的矫正p值界限为pcrit=j,而且从上到下开始比较，如果显著将继续比较，如果不显著将停止继续比较（比如到第三行发现不显著之后，下面的就认为也不显著不需要继续比较了）。前两个关注familywise error rate，Benjamini and Hochberg关注的是false discovery rate (FDR)，FDR= numberof falsely rejected null hypotheses / total number of rejected null hypotheses。这个方法具体的算法和Holm的差不多，只不过将实际得出来的p值从小到大排列后，分别从1到6进行赋值，pcrit=j/k。还有一个不同的是，它是从下往上的顺序和矫正后的P值进行比较，如果发现显著性后往上就不比较了（即从6比到4发现第4这显著，321这些就不比较了，认为他们都是显著的）。
Bonferroni’s and Tukey’s HSD tests在控制一类错误中相对严格，但是会减少统计检验力。在统计检验力方面FDR>Holm> Bonferroni。
当数据不是正态分布（偏离的不严重）时，大多数多重比较方法表现的还不错，但是当群体大小或者方差不同时，他们的表现就不如人意。这个时候就需要要bootstrap, trimmed means, or M-estimators这三种方法了。
10.6.1.Packages for one-way ANOVA in R 单因素方差分析包
本章中将用到以下的包：
car (for Levene’s test),
compute.es (for effect sizes)
ggplot2 (for graphs),
multcomp (for post hoc tests事后检验),
pastecs (for descriptive statistics描述性统计),
WRS (for robust tests鲁棒性检验).
第一步，通过执行以下命令来安装这些包：
install.packages("compute.es")install.packages("car")install.packages ("ggplot2")install.packages("multcomp")install.packages("pastecs")install.packages("WRS", repos="http://R-Forge.R-project.org")第二步，通过执行这些命令来加载这些包:
library(compute.es)library(car)library(ggplot2)library(multcomp)library(pastecs)library(WRS)10.6.2.General procedure for one-way ANOVA 单因素方差分析的一般程序
1、输入数据；
2、探索数据:数据作图和描述性统计；检查数据分布和方差齐性（使用Levene）；
3、基本方差分析:主要的方差分析，根据结果，判断是否需要进行鲁棒性检验；
4、计算对比或事后测试:在进行了主要的方差分析之后，可进一步进行对比或事后检验。具体选择哪一种方法取决于步骤2中的结果。
10.6.3.Entering data输入数据
样本量很小，采用如下方法输入数据：
libido<-c(3,2,1,1,4,5,2,4,2,3,7,4,5,3,6)#创建了一个libido的变量，并为其赋值
dose<-gl (3,5, labels = c("Placebo", "Low Dose", "High Dose"))#使用gl函数创建了一个dose的变量，将其分为三个组，每组有5个样本
viagraData<-data.frame(dose, libido)#以上两个变量和其样本编入viagraData中
10.6.4.One-way ANOVA using R Commander 使用R指令进行单因素方差分析
1、从Viagra.dat加载数据。（Data⇒Import data⇒from text file, clipboard, or URL… menu）
2、描述性统计和验证假设
（1）Levene检验：即方差齐性检验。（Statistics⇒Variances⇒Levene’s test… menu）。
labelled Groups中填写因素（即自变量，本例中为dose）；labelled Response Variable中填写结果变量（即因变量，本例中为libido）。Levene检验默认基于偏离median的离差程度，因此默认median选项，但是可以更改。结果输出见10.6.5.
（2）ANOVA的初步检验：首先按照如下菜单栏进行点击,（Statistics⇒Means⇒One-way ANOVA… menu）；然后，在Enter name for model这个选项框中，为模型起一个名字（如：viagraModel），并将libido和dose选入合适的位置；最后，你不能使用R commander进行详细的事后比较，但是，如果想要基础的事后比较，请点击pairwise选项。
10.6.5.Exploring the data 探索数据
1、可以使用by()函数获得每组的部分描述性信息。这个函数的一般形式如下：by(variable, group, output)。通常来讲，Variable是因变量（在本例中为libido），group是自变量（在本利中为dose），output是想要得到的描述性统计结果（比如说均值）。
2、by(viagraData$libido, viagraData$dose, stat.desc)
#结合使用by()函数和pastecs包的stat.desc()函数可以得到每组的描述性统计。
3、输出结果：
nbr.val（有效以样本个数）；min（最小值）；max（最大值）；range（范围）；median（中位数）；mean（平均值）；var（方差）；std.dev（标准差）；SE.mean（标准误） and CI（置信区间）。
4、正式的ANOVA之前需要用car包进行Levene’s test（leveneTest(outcome variable, group, center = median/mean)）。我们期待的Levene’s test结果是不显著的，这说明没有明显的组间差异。但是如果是显著的，那么我们应该进行Welch’s F检验或者ANOVA稳健性检验。
10.6.6.1. When the test assumptions are met 当满足测试假设时
1、由于方差分析是一般线性模型的特例，因此可以使用一般线性模型函数lm()进行分析：
viagraModel<-lm(libido~dose,data = viagraData)除此之外，还可以使用专用的aov()进行分析。一般形式：
newModel<-aov(outcome ~ predictor(s), data = dataFrame, na.action = an action))lm()和aov()输出的结果是一样的，实际上，工作是lm()做的，但是，aov()输出的结果更加符合我们的阅读习惯。
2、aov()函数中的na.action是一个可选命令，当数据集有缺失时，使用：
na.action = na.exclude可以将有缺失值的样本剔除。
3、得出模型之后，执行：
summary(viagraModel)得出模型信息，具体包括自由度、SS、MS、F值和p值。其中，SS（residual sum of squares）表征非系统性变异；MS（mean squares）表征均方差；F值表示模型是否显著；p值表征模型出现错误的可能性。
4、作图：plot(viagraModel)。第一张图用于检验方差齐性，显示的是三个组的平均分布，如果是漏斗形的，则结果不佳，如图所示，则结果很好；第二张图是Q-Q图，用于检验模型中残差是否为正态分布，残差越靠近对角线则越接近正态分布，由于这种方法（是否靠近对角线）具有主观判断特征，因此，有时会有一些争议，因此，或许应该有一个更加具有鲁棒性的方法进行残差正态性检验。
另外：单尾检验和双尾检验的问题。不能按照双尾的结果（不显著）得出p之后除以2得出单尾的结果（显著），从而得出显著的结论，这是因为，如果是单尾的，则是有方向性的（如：猫的尾巴比狗的长），但是很明显，方差分析不具有方向性，因此，将显著性值减半的做法是无效的。
10.6.6.2. When variances are not equal acrossgroups 当组间的方差不相等时
1、如果Levene’s test没有通过，那么使用：
oneway.test(outcome ~predictor, data = dataframe)进行检验，得出Welch’s F值，并在报告中说明，自由度已经被调整。最后依据运行结果中的p值，确定因变量在各个组之间是否有显著差异。
10.6.6.3.Robust ANOVA – it’s not for the weak of heart 方差分析鲁棒性—不适合脆弱的小心脏
注：本节实际上是用各种均值的替代统计量进行差异性检验。
1、加载检验鲁棒性的命令参见5.8.4
2、Wilcox’s 函数（需要将常用的长数据转换为宽数据）
（1）转化为宽数据的命令一般格式：
newDataFrame<-unstack(oldDataFrame,scores ~ columns)通过这一转变，可以将不同自变量（组别）按列分开。
（2）鲁棒性检验的第一个命令的一般格式：
t1way(dataFrame, tr = .2,grp = c(x, y, …, z))这个命令基于被修剪的均值。其中，tr是需要修剪的比例，grp用于指定想要分析的列（组别）。检验结果用于说明各组之间被修剪的均值是否有显著性差异。
（3）鲁棒性检验的第二个命令的一般格式：
med1way(dataFrame, grp =c(x, y, …, z))这个命令基于中位数而不是均值。检验结果用于说明各组之间中位数是否有显著性差异。
（4）鲁棒性检验的第二个命令的一般格式：
t1waybt(dataFrame, tr =.2, alpha = .05, grp = c(x, y, …, z), nboot = 599)这个命令是在修剪平均值中加入bootstrap（可理解为有放回抽样）。Alpha用于设置第I类错误错误率，nboot用于设置bootstrap次数。检验结果用于说明各组之间在bootstrap下，均值是否有显著性差异。
10.6.7.Plannedcontrasts using R 使用R进行计划比较
1、在有viagraModel的情况下，执行summary.lm(viagraModel)，这样能够对线性模型的各个参数进行对比。（结果怎么看，以Output 10.8为例：low dose效应是指对应于安慰剂的低剂量的效应，这里是不显著的）。
2、进行趋势分析或者有计划的分析的一般形式：
contrasts(predictor variable)<-contrast instructions具体的函数形式如下：
contr.helmert(n)contr.poly(n)contr.treatment(n, base = x)contr.SAS(n)n是预测变量的个数；contr.treatment和contr.SAS允许设置基线组
10.6.7.1. Your owncontrasts
1、我们需要告诉R给每个组分配什么权重。赋权重，如：contrast1<-c(-2,1,1)
10.6.7.2. Trendanalysis（趋势分析）
为了进行趋势分析，我们可以使用con .poly()。重要的是，我们已经按照有意义的顺序对预测变量组进行了编码。我们希望安慰剂组的性欲最小，低剂量组性欲增加，高剂量组性欲再次增加。为了发现有意义的趋势，我们需要按升序对这些组进行编码。我们将安慰剂组编码为最低值1，低剂量组编码为中值2，高剂量组编码为最高值3。如果我们以不同的方式对这些组进行编码，这将影响是否检测到一种趋势，以及如果检测到一种趋势，它是否具有统计意义。
首先第一步，我们需要使用：
contrasts(viagraData$dose)<-contr.poly(3)来设定对比中有多少组，这里我们是3组，然后是：
viagraTrend<-aov(libido ~ dose,data = viagraData)来设置模型参数。最后，使用：
summary.lm(viagraTrend)来获得结果，输出的结果10.10（下图）分解了实验效果，我们看看它是否可以用数据中的线性(剂量. l)关系或二次(剂量. q)关系来解释。我们可以看到，这个比较测试了组间均值是否以线性方式增加。需要注意的最重要的事情是t的值和相应的显著性值。对于线性趋势t = 3.16，且该值在p = .008时显著。因此，我们可以说，当伟哥的剂量从零增加到低剂量到高剂量时，性欲也相应增加。而二次关系是不显著的。
10.6.8 执行post-hoc
您如何在R中执行事后测试取决于您想要执行哪个测试。Bonferroni和相关的方法(如Holm和Benjamini Hochberg)是使用pairwise.t.test()函数完成的，该函数是R base系统的一部分。但是，可以使用multcomp()包中的glht()函数来完成Tukey和Dunnetts测试(以及其他一些我们不打算讨论的测试)。最后，Wilcox(2005)在他的函数lincon()和mcpp20()中实现了一些鲁棒的方法。本节按照这些不同的方法进行划分。
10.6.8.1. Bonferroni and related methods
使用这句：pairwise.t.test(outcome, predictor, paired = FALSE, p.adjust.method ="method")来完成Bonferroni矫正，其中，结果是结果变量的名称(在本例中是libido (viagraData$libido))，预测变量是分组变量的名称(在本例中是dose (viagraData$dose))，paired是一个逻辑语句，默认情况下为FALSE，但可以设置为TRUE(大写字母matter)。这指定了是否需要配对t检验。对于这些数据，我们有独立的组，所以我们不希望配对t检验，默认为FALSE也可以，但是我们将在第13章中重新讨论这个选项。p.adjust.method指的是用那种方法来进行p值的矫正，因此，我们写的代码就是：
pairwise.t.test(viagraData$libido,viagraData$dose, p.adjust.method ="bonferroni")或者是：
pairwise.t.test(viagraData$libido,viagraData$dose, p.adjust.method = "BH")我们得到的结果会如下图10.11所示：首先，让我们来看看Bonferroni的修正值:将安慰剂组与低剂量组进行比较，并发现一个不显著的差异(但与高剂量组相比，有显著性差异0.025小于0.05)。同时可以看到，BH的结果其实和Bonferroni的结果是一样的。
10.6.8.2. Tukey and Dunnett
Tukey和Dunnett可以使用multcomp包中的glht()函数来实现(请记住安装和加载它)。这个函数的一般形式如下：
newModel<-glht(aov.Model, linfct = mcp(predictor = "method"), base = x)其中，newModel是一个包含来自特定测试的信息的对象。要查看这些信息，我们可以使用summary(newModel)来查看基本的事后测试，使用confint(newModel)来查看置信区间。Aov.Model是已经用aov()函数创建的模型的名称(在本例中是viagraModel)。预测变量是分组变量的名称(在本例中是dose (viagraData$dose))，linfct = mcp(predictor
= “method”)指定要对p值应用哪种校正。您可以将上面命令中的方法替换为Dunnett、Tukey、Sequen、AVE、Changepoint、Williams、Marcus、McDermott、lawilliams和GrandMean。
只有在指定Dunnett时才使用base。此选项允许您使用组号指定基线组。在这种情况下，如果我们想要安慰剂作为基线，我们将使用base = 1，但如果我们想要高剂量组，我们可以指定base = 3。
基于以上的信息，我们的代码将会如下面所示：
postHocs<-glht(viagraModel, linfct = mcp(dose = "Tukey"))summary(postHocs)confint(postHocs)结果将如下图所示：
结果图显示了三个比较(低剂量和安慰剂,高剂量和安慰剂,高剂量与低剂量), 估计值(组间均值之差)、与均值之差相关的标准误差、t检验(t检验就是均值之差除以标准误差，所以第一次对比是1/0.8869 = 1.127)及其相关的p值。与测试前一节中一致,该输出确认高剂量和安慰剂组之间的显著差异,t = 3.16, p < . 05,但不是低剂量组和安慰剂之间,t = 1.13, p =0.52,和高剂量,t = 2.03, p =0.15。置信区间(输出10.13)也证实了这一点，因为对于高剂量组和安慰剂组的比较，它们没有跨越零，这意味着组间均值的真实差异可能不是零(即,没有区别);相反地，对于另一种情况，置信区间跨越0，意味着均值之间的真实差可以是0。
输出10.13
10.6.8.3.Run for cover – it’s robust post hoc tests（鲁棒的post-hoc分析方法）
与稳健的方差分析一样，为了运行稳健的事后测试，我们需要：
(1)来源Rand Wilcox的函数(有关如何做到这一点，请参阅第10.6.6.3节);
(2)输入宽格式数据因此，
因此，我们将使用在10.6.6.3节中创建的viagraWide对象。我们将使用两个函数:lincon()，它基于裁剪方法;和mcppb20()，它使用百分比引导来计算p值，并调整组平均值。尤其是后一种方法，似乎很擅长控制I类型的错误率。这些函数的一般形式类似于t1way()和t1waybt()，我们在本章前面已经介绍过了。
这是这两个函数的一般使用方法：
lincon(dataframe, tr= .2, grp = c(x, y, …, z))mcppb20(dataframe, tr= .2, nboot = 2000, grp = c(x, y, …, z))在我们这里如下设置即可，其他保持默认即可：
lincon(viagraWide)mcppb20(viagraWide)lincon结果如下，可以看到结果和前面的一样，1组和3组的差异显著：
mcppb20的结果如下：可以看到结果和前面的一样，也是1组和3组的差异显著：
10.7效应量的计算
在R中比较两组别之间不同的效应量用如下语句
计算正交对比的时候可以用以下语句
10.8报告单因素方差分析的结果
当我们报告一个方差分析时，我们必须给出F比率的细节以及计算它的自由度。对于这些数据的实验效果，通过用效果的均方除以残差的均方得到F比值。因此，用来评估F-ratio的自由度是模型效果的自由度(dfM = 2)和模型残差的自由度(dfR= 12)。因此，正确的报告主要发现的方法是：
There was a significant effect of Viagra on levels of libido, F(2, 12) = 5.12, p < .05,ω = .60。
请注意，F-ratio的值之前是该效果的自由度的值。此外，我们很少声明F比率的确切显著性值:相反，我们报告显著性值p小于标准值.05，并包含一个效应大小度量：ω。线性对比也可以用同样的方法来描述:
There was a significant linear trend, F(1, 12) = 9.97, p < .01, ω = .62, indicating that as the dose of Viagra increased, libido increased proportionately.
要注意的是，自由度已经发生了变化，以反映F比率是如何计算的。我还包括了一个效果大小度量(试着计算一下这个，就像我们计算主要的F比率一样，看看是否得到相同的值)。此外，我们现在已经报告了F值在小于标准值0.01时是显著的。我们还可以报告我们计划的对比或组比较：
首先是计划的成组比较：
Planned contrasts revealed that taking any dose of Viagra significantly increased libido compared to having a placebo, t(12) = 2.47, p < .05 (one-tailed), and that taking a high dose significantly increased libido compared to taking a low dose, t(12) =2.03, p < .05 (one-tailed)
其次是事后分析比较：
Despite fairly large effect sizes, Bonferroni tests revealed non-significant differences between the low-dose group and both the placebo, p = .845, d = -0.77, and highdose, p = .196, d = -1.24, groups. The high-dose group, however, had a mean almost 2 standard deviations bigger than the placebo group, p = .025, d = -1.93.
~~~~~~~~~~~~本期分享人~~~~~~~~~~
苏双 阿姆斯特丹大学荷兰 发展心理学
王翠翠 北京师范大学 心理学
袁树歌 西南大学（重庆） 心理学
郑美慧 哈尔滨工业大学 神经经济学
吴国伟 江苏师范大学(江苏省徐州市) 神经语言学
全湘 武汉大学（湖北省武汉市） 心理学
~~~~~~~~~~~~~排版：杨伟文~~~~~~~~~~
往期热门：
Play with R 第九期：Comparing two means