Play with R 第11期：协方差分析
11.1 什么是协方差分析(ANCOVA)
在研究中，不属于主要实验干预(experimental manipulation)但是对因变量有影响的连续变量，被称为协变量(covariates)，它们可以被包括在方差(ANOVA)分析里。当我们测量协变量，并将其纳入方差分析中，我们就称之为协方差分析(ANCOVA)。
在方差分析中纳入协变量的两个原因：
1） 减少组内误差(within-group error variance). 在讨论方差分析和t检验时，我们习惯了这样的想法：我们通过比较“实验能解释的数据变化量”和“不能解释的变化量”来评估实验的效果如果，而我们可以用其他变量（协变量）来解释这种“不能解释的变化量”（SSR），那么我们就可以减少误差方差，从而更准确地评估自变量（SSM）的影响。
2） 消除混淆：在任何实验中，都可能有未测量的可能混淆结果的变量（即，除了影响结果变量的实验干预之外的变量）。如果已知任何变量会影响因变量，那么协方差分析非常适合去除这些变量的偏差。一旦一个潜在的混淆变量被识别出来，它就可以被测量，并作为协变量纳入分析。
还有很多其他相关的原因，解释了协方差分析的必要性，感兴趣可以参考(Stevens, 2002; Wildt & Ahtola, 1978)。
在前面的章节中，如果在方程10.2中纳入协变量，会如下所示：
11.3 协方差分析中的假设和问题
协方差分析和方差分析的假设基本相同，但协方差分析中还需要考虑两个额外因素：
(1) 协变量和处理效果之间的独立性；
(2) 回归斜率的同质性
11.3.1 协变量和处理效果之间的独立性
前面提到，协方差分析的一个重要作用，是通过协变量来解释一些误差方差，减少组内误差。但是，如果需要这一点成立的话，协方差必须和处理效果独立。
图11.2 现实了三个不同的情境。
Part A显示了基本的方差分析，与第10章中的10.4相似；它展示了实验中的处理小于可以被分解为两部分，一部分是实验或处理额效果，另一部分是误差或未被解释的方差（即，能够影响因变量但未曾测量的因素）。
Part B 显示了协方差分析的理想情境，其中协方差解释了违背解释的方差的很小一部分。换言之，它完全独立于处理效果。这是协方差分析唯一合适的情况。
Part C显示了人们经常在不该使用协方差分析的情况下使用它。协变量与实验中的处理效果存在重叠，实验干预产生的效果被协方差的效果混淆了。在这种情况下，协方差会在统计上减少实验效果，因为它解释了一些本该由实验干预来解释的变异。
那么，什么情况下人们容易误用协方差分析呢？
一是被试没有被随机分配到实验干预的条件中。例如，焦虑和抑郁是高度相关的，焦虑的人们很容易抑郁。所以如果你想要比较焦虑组和非焦虑组的被试在某项任务上的表现，我们可能会想到把抑郁当做一个协变量，以便看到焦虑对表现的纯影响。但是不能这样做，因为这就会造成Part C所示的情况：抑郁会包含一部分焦虑所产生的效果。在统计上来看，我们知道焦虑和抑郁共享一部分的变异，我们不能把共享的变异划分为“焦虑带来的变异”和“抑郁带来的变异”，它们只能是“共享的变异”。
另一种情况是如果你凑巧发现你的实验组、控制组、对照组等等在年龄上存在差异，你把年龄作为一个协变量，也不会解决问题——它仍然被实验干预混淆了。协方差分析不是解决此类问题的万灵药。
这个问题只能通过随机分配被试到各个实验的组别中来解决，或者通过按照协变量来匹配实验对象来解决。例如，你可以尝试着发现低焦虑、高抑郁的被试。在实验前，我们可以先通过t检验或方差分析，比较高焦虑和低焦虑的被试是否在抑郁程度上存在显著差异。如果没有显著差异，我们就可以使用焦虑作为协变量。
11.3.2 回归斜率的同质性
当进行协方差分析的时候，我们要看到结果变量和协方差之间的总体关系。我们用整个数据集来拟合回归线。在拟合整体模型的时候，我们假设整体关系对所有组别都是一样的。然而，如果协变量与结果变量的关系（比如说，正相关），在某一组中成立，在其他组中不成立，那么整体模型就不正确。这个假设非常重要，被称为回归斜率同质性假设
Figure 11.3展示了.伴侣的libido，协方差，被试的libido，以及结果变量，这四个变量在三种实验条件下的关系。不同的标志代表不同的实验组别，如黑色圆圈代表安慰剂组，浅蓝色三角形代表低水平自变量的组，深蓝色正方形代表高水平自变量的组。我们发现三个组的斜率不太一样：浅蓝色三角形拟合出来的线条的斜率，与黑色原型拟合出来的线条斜率很接近，它们的斜率是同质的；但是浅蓝色三角形，与深蓝色正方形，拟合出来的线条斜率大相径庭，它们的斜率是异质的。
虽然，在传统的ANCOVA中，回归斜率的异质性是一件坏事，但在某些情况下，您可能真的期望回归斜率在组内不同，这本身就是一个有趣的假设。当研究在不同的地点进行时，你可以合理地预期你在这些地点得到的效果会略有不同。例如，如果你有一个新的治疗背痛的方法，你可能会让几位理疗师在不同的医院试用。你可能认为这些医院的治疗效果不同（因为治疗师的专业知识不同，他们看到的病人会有不同的问题等等）。回归斜率的非均质性本身并不是一件坏事。如果您违反了回归斜率同质性的假设，或者回归斜率的可变性本身是一个有趣的假设，那么您可以使用多级线性模型(multilevel linear models)对这种变化进行建模（参见第19章）
11.4 ANCOVA using R
11.4.1 协方差的R包
我们需要用到的R包有car (for Levene’ s test, Type III sums of squares, 俗称方差齐性检验)，compute.es (for effect sizes, 效应量)，effects (for adjusted means, 调整平均数)，ggplot2 (for graphs, 画图)， multcomp (for post hoc tests, 事后检验)， pastecs (for descriptive statistics, 描述性统计), and WRS (for robust tests)。
你需要事先安装这些包，通过以下代码来执行安装命令：
然后加载这些包，通过一下代码来执行命令：
11.4.2 协方差分析的基本流程
1. 导入数据
2. 探索数据：从画图和计算描述性统计开始。你还需要检查分布假设（distributional assumptions）和通过Levene’s test 来进行方差齐性检验。
3. 检查协变量和任一自变量是独立的：你需要将协变量作为结果结果，自变量作为预测变量来进行方差分析，以检验协变量与自变量的各水平确实没有显著性差异。如果你的到一个显著的结果，那么就可以伤心地停止分析。
4. 做协方差分析：如果步骤2和3 都很顺利，就可以进行主要的协变量分析。基于步骤2的发现，你需要做robust test。
5. 比较差异或事后检验: 通过这些分析你可以进一步发现组间差异。
6. 检验回归斜率同质性：重新进行包含有检验自变量和协变量交互作用的协方差分析。如果他们有显著性的交互作用，那么他们就不具有回归斜率同质性。
你需要依次进行以上步骤。
11.4.3 导入数据
数据示例见Table 11.1。在文件ViagraCovariate.dat 可以看到该数据。
Table 11.1 展示了被试的力比多和他们伴侣的力比多。Table 11.2 展示了这些数据的平均值和标准差。
通过以下代码， 执行导入数据命令：
也可以自己输入数据，为了让导入的数据和Table 11.1 的格式一样，你需要创建一个编码变量Dose， 并编码为1= placebo， 2= low dose, 3= high dose。在每个条件下有不同的被试数， 因此在这一列中你需要先输入9个1（因此前九行就被赋值为1），接下来8行赋值为2， 最后13行赋值为3。此时，你需要输入有一列（DOSE）包含30行的数据。接着，创建第二个变量libido，并且输入30个对应被试力比多的数值。最后创建第三个变量partnerLibido。然后输入与其伴侣的力比多对应的数值。
我们可以直接在R中输入：
通过这一命令创建有30个力比多数值的变量libido，变量 partnerLibido和变量dose （运用rep() 功能重复数字1九次，数字2八8次和数字3十三次）。我们还需要把数值型变量dose转成一个因子（例如类别变量）。我们可以通过以下代码来执行此命令：
记住！我们已经把变量dose的水平定义为1，2和3（levels = c(1;3)），并且我们需要把这些水平标记为Placebo， Low Dose 和 High Dose（labels = c(“Placebo”, ”Low Dose” , ”High Dose”)）。最终，我们可以通过以下命令，把这些变量融合为一个数据框viagraData：
最终的数据看起来如下图：
11.4.4 使用R Commander 进行协方差分析
R Commander 没有菜单可以直接进行协方差分析。然而，协方差分析是简单的回归， 因此你可以通过菜单Statistics⇒Fit models⇒Linear regression...。然而，R Commander 不能很好地处理分类变量，而且你不能控制变量的顺序（这是非常重要的），因此不推荐使用R Commander进行协方差分析。
11.4.5 探索数据
这一小结通过箱线图对比每个药剂组力比多的量以及协变量（父母力比多的量）对因变量的作用（这一作用能揭示回归斜率的同质性）。
如图所示：
如上图所示，受试者的力比多水平似乎随着用药剂量的增加而增加，而这一效应在受试者父母上呈现出了相反的模式。此外，受试者的力比多分数的波动范围似乎比家长的更大。
接下来，使用 by() 以及 pastecs 包中的stat.desc()功能得到每个分组的描述统计结果。具体代码如下：
结果如下表：
接着计算Levene’s测验的结果（方差齐性检验）。在本例中，我们可以通过如下代码了解各组间不同药物剂量对力比多水平的影响：
结果如下：
结果显示样本数据没有违背方差齐性假设，若p值显著，则需要进行稳健方差分析（robust version of ANOVA）。
有一个简便的方法也可以用来判断数据是否违背方差齐性假设：
在本例中，可以先求三个不同药剂组的方差，然后用最大的方差除以最小的方差（3.20(placebo), 2.13(low dose)，4.49 (high dose)，得到的结果为2.11。接着可以通过下图来进行简略的判断：
由上图可知，样本数为10的情况下，比较3个组（本例为3个药剂组）不同方差变异的决断值大概为5， 而之前求得的最大变异为2.11，故可以推断没有违反方差齐性假设。
11.4.6 判断预测变量和协变量的独立性
在本例中，协变量为父母的力比多水平。我们可以通过以其为因变量，以药物剂量为自变量的方差分析来判断二者是否独立。
结果显示P值并没有达到显著性水平，故原假设不能被推翻，即药剂组之间的父母力比多水平没有明显差异，故自变量（不同药剂组）和协变量（父母力比多水平）之间相互独立。
11.4.7 构建协方差分析模型
R soul tip 11.1
Order matters
变量进入模型顺序的重要性
我们将预测变量放入模型的顺序对于整体方差分析会产生重大影响——这点非常令人费解。幸运的是，这并不会影响模型参数（即b值）。我们来看一个例子。
首先，我们来拟合一个模型covariateFirst，先放入partnerLibido，随后放入dose ，预测libido.
covariateFirst<-aov(libido ~ partnerLibido + dose, data = viagraData)
summary(covariateFirst)
运行以上代码，得到结果
> covariateFirst<-aov(libido ~ partnerLibido + dose, data = viagraData)
> summary(covariateFirst)
可以看出，partnerLibido对于被试libido的影响并不显著，但是dose对于被试libido具有显著影响。接下来，我们调整放入变量的次序，将dose先放入，其次partnerLibido,运行模型结果如下：
> doseFirst<-aov(libido ~ dose + partnerLibido, data = viagraData)
> summary(doseFirst)
我们得到完全相反的结果!partnerLibido对于被试的libido具有显著效应，而dose的影响变的不显著。原因在于，当R在进行模型拟合的时候默认使用Type I来计算sums of squares, 另外可替代方法（也是被很多统计包采用）的是用Type III 来计算总方差(sums of squares). 对于以上的两个模型可以加入一行语句：
> Anova(covariateFirst, type = "III")
Anova Table (Type III tests)
Response: libido
> Anova(doseFirst, type = "III")
Anova Table (Type III tests)
Response: libido
然后，可以看出两个模型结果完全一致。或许你对Type I III是啥感到困惑……简的超级大脑来告诉你！
JANE SUPERBRAIN 简的超级大脑
11.1 不同类型的方差和
我们可以用四种不同的方法来计算方差和，即Type I,II, III, IV. 接着前边的例子，我们用协变量伴侣的libido和自变量（dose）,以及二者的交互作用（partnerLibido × dose）来预测被试的libodo.
Type I 方差和简单来说就是类似于进行层级回归hierarchical regression: 我们将一个预测变量f放入模型，接着再放入第二个预测变量。第二个变量的预估将会在第一个变量之后，如果我们放入第三个变量，对其估算就会在第一和第二变量之后。换而言之，我们放入变量的顺序会对模型产生影响（参见上面R soul tip）。
Type III,计算方差和不同于Type I 之处在于所有效应的估算都会考虑到模型中所有其它变量（而非仅仅是先放入的）。这一过程可以看作是回归模型中的迫入法forced entry 。
Type II计算方差居于Type I 和Type III之间，既考虑了所以其它效应的同时，对于higher-order高阶变量不予估算 。
all effects are evaluated taking into consideration all other effects in the model except for higher-order effects that include the effect being evaluated.
在上面的例子中，意味着does的效应计算将会在ParternerLibo效应计算之后，（注意不同与Type III 方差和，不会考虑交互作用）；类似地，ParternerLibo效应的计算，将只会在does效应的计算之后。
Type IV,与Type III完全相同，主要是为了解决变量中有缺失值的问题。
那么，我们到底该使用哪种方差和计算方式呢？
· Type I: 除非所有的变量完全独立（现实中很难出现），不然这一方法很难准确估算每个变量的主效应。
· Type II: 如果你对主效应感兴趣，那么你应该使用Type II 方差和。不同于Type III, 这一方法能准确的计算主效应，因为其计算会忽略任何与主效应相关交互作用。
· Type III: 在很多统计包中使用，其优势在于与交互作用相关变量的主效应仍然有意义（因为其计算考虑到了交互作用）。
如果我们要计算Type I 方差和，在协方差分析中，我们首先纳入协变量，
viagraModel<-aov(libido ~ partnerLibido + dose, data = viagraData)
注意：在这一模型中我们首先是放入协变量（partnerLibido）然后是自变量dose, 这意味着dose 效应的计算将会是在partnerLibido之后。然后，我们可以选择以上四种的任意一种方差和计算方式查看结果。
Anova(modelName, type = "III")
第二个问题是，我们选择何种的对比方式。这一问题与上一章节讨论的哑变量的比较相关，R默认为使用非正交对比，因此如果我们不改变对比方式，使其正交化，Type III计算的方差和将会出错。因此，我们必须要么用Helmert 对比
contrasts(viagraData$dose)<-contr.helmert(3)
或者通过第10章节中讲述的方法：
contrasts(viagraData$dose)<-cbind(c(-2,1,1), c(0,-1,1))
本例中我们使用aov 功能构建协方差分析模型，如下图：
加入协变量之后的模型如下：
注意在使用aov功能的过程中，变量放入模型的顺序可能会对模型的结果产生影响，这是因为aov默认计算 Ⅰ型平方和。因此每个变量的效应是在之前纳入变量效应的结果之上得出的，因而顺序会有影响。我们可以通过设置III型平方和来规避这一影响。具体代码如下：
Anova (viagraModel, type = "III")
11.4.8 协方差模型的解读
11.3 显示了协方差模型，首先看显著值，显然协变量对于因变量具有显著的预测效应，因其显著性小于 .05。因此，个体的libido水平受到伴侣libido的影响。更为有趣的是，剔除伴侣libido之后，Viagra的效应变的显著。
> Anova(viagraModel, type="III")
Anova Table (Type III tests)
Response: libido
虽然我们似乎可以看出安慰剂组合两个实验组存在显著差异，实际上我们还不能解释这两个组的平均数，因为这一效应并未考虑协变量的影响。为了得到准确的调整后的组均值，我们需要使用以下函数
object<-effect("name of effect", modelName, se=TRUE)
summary(object)
object$se
> adjustedMeans<-effect("dose", viagraModel, se=TRUE)
> summary(adjustedMeans)
> adjustedMeans$se
[1] 0.5532329 0.4060868 0.3219445 0.3497028 0.4699339
11.4.9 planned contrasts in ANCOVA
协方差整体模型并不能告诉我们哪些平均数存在差异，因此我们需要将整体的效应分化，然后按照我们在运行协方差模型之前的设定进行对比。为了进行对比我们可以使用summary.lm方程
summary.lm(viagraModel)
11.4.10 解释协变量
参数估计告诉我们如何解释协变量
(1). 如果协变量的 b值为正，这意味着协变量和结果变量有正向关系(随着协变量的增加，结果变量增加)。
(2). 如果 b值为负，则意味着协变量和结果变量具有负相关关系(协变量增加， 结果变量减少)。
11.4.11 ANCOVA分析中的事后检验（Post hoc tests in ANCOVA）
在进行事后检验时，由于我们想测试调整后的平均值之间的差异，只能使用glht()函数，不能用test()函数。只能进行Tukey或Dunnett’s事后检验。有几个步骤：首先将ANCOVA模型输入R；其次，使用summary()和confint()函数查看控制台中的事后检验。
对于viagraModel，我们可以执行下面的程序：
postHocs<-glht(viagraModel, linfct = mcp(dose = "Tukey"))
summary(postHocs)
confint(postHocs)
11.4.12 ANCOVA分析中的图（Plots in ANCOVA）
之前讲到使用aov()函数可以自动出来一些图以便于我们进行假设检验，我们也可以执行下面的程序来生成一些图：
plots(viagraModel)
执行程序就会生成四个图，看前两个比较重要的图：
左边的图主要检验方差齐性，如果该图呈现漏斗形状 (分数在某些点上的差距比其他点上的差距更大)，这意味着残差可能是异方差的(比较麻烦，一件坏事)。
右边的图是正态分位数图（Q-Q图，看第五章），它告诉我们模型中残差的正态性。我们一般希望残差是正态分布的，也就是说图上的点应该在对角线上。在右图中，有些点明显偏离了对角线，这对模型来说不是好消息。
两幅图表明，进行稳健型的协方差检验比较合适。
11.4.13 总结
这个例子说明了ANCOVA如何通过考虑混杂变量来帮助我们实施更严格的实验控制，从而让我们对实验操作的效果有一个“更纯粹”的测量。
11.4.14 回归斜率的齐性假设的检验
回归斜率齐性：这意味着协变量和结果变量之间的关系在预测变量的不同水平上应该是相似的。
图11.3显示了三组中伴侣的性欲和被试者的性欲之间关系的散点图。在本例中协变量和结果变量分别为伴侣的性欲和被试者的性欲，预测变量分为三个水平，即安慰剂组，低水平自变量的组，高水平自变量的组。这个散点图显示，虽然这种关系在低剂量组和安慰剂组中具有可比性，但在高剂量组中出现了不同。
想要检验回归斜率齐性的假设，需要再次运行ANCOVA，但要包括协变量和预测变量之间的相互作用，有三种方法可以实现：
第一种：重新指定整个模型。我们可以通过用冒号连接变量名来包含交互项（比如伴侣性欲和伟哥剂量之间的相互作用可以用R表示为partnerLibido:dose，两者顺序可调换)。因此，将这种交互包含在 ANCOVA模型中，我们可以执行下面的命令。:
该命令创建一个称为hoRS(回归斜率同质性的缩写)的模型，该模型包括协变量、自变量及其相互作用。
第二种：以通过指定variable1*variable2作为预测变量，将变量及其交互作用包含在同一个模型中。这样做不仅会影响交互作用，还会影响单个变量的效果。我们可以执行命令:
第三种：更新原始的ANCOVA模型(viagraModel)，使用update()函数包含交互项(参见R ' s Souls的技巧7.2)。viagraModel已经包含了伴侣者性欲和伟哥剂量，所以我们需要做的就是通过添加‘+ dose:partnerLibido’来添加交互项，如下所示:
“.~.”简单来说就是“保持和以前一样的结果变量和预测变量”，而“+ partnerLibido: dose”意味着“增加交互项”。这种方法是最快。
执行三种方法中的一种，以创建hoRS对象，然后使用anova()函数得到类型III的平方和:
输出11.9显示了ANCOVA的主要汇总表，包括交互项。伟哥的剂量和伴侣的性欲影响仍然显著，但我们感兴趣的主要是相互作用项，，所以看看相互作用的共变量的显著值(partnerLibido:dose)，如果这种影响很显著，那么回归斜率齐性的假设就被打破了。这里的影响是显著的(p<.05)，因此这个假设是站不住脚的。尽管从图 11.3 所示的关系模式来看，这一发现并不令人惊讶，但它确实引起了对主要分析的关注。这个例子说明了为什么检验假设很重要，而不是盲目地接受分析结果。
11.5 稳健型的ANCOVA检验
与单因素方差分析一样，Wilcox(2005)描述了一组用于进行单因素方差分析的稳健过程。为了访问它们，我们再次需要加载WRS包(参见5.8.4节)。我们讨论两个函数：ancova()和ancboot()，它们可用于比较两组之间的切尾均值（切尾均值是指在一个数列中，去掉两端的极端值后所计算的算术平均数）。为了使分析摆脱回归斜率齐性的限制以及其他分布假设，这些检验比较沿协变量不同点的平均值。换句话说，与其假设协变量和结果变量之间的关系在两组中是恒定的，不如找出斜率相同的 5 个点(即两组中结果变量和协变量之间的关系大致相同的 5 个协变量值)。然后比较这五个点的平均值。然后比较这五个点的切尾均值，看看它们是否不同。这个过程让我们知道组间的差异是如何随着协变量的变化而变化的。
建议使用ancova()函数执行上述分析。虽然ancboot()函数也能执行相同的分析，但是它使用bootstrap-t方法来计算置信区间。执行这些功能的基本形式如下:
命令解释：
covGrp1指第一组包含协变量数据的变量；
dvGrp1指第一组包含结果变量数据的变量；
covGrp2指第二组包含协变量数据的变量;
dvGrp2指第二组包含结果变量数据的变量。
tr =表示切尾比例，默认为20%，可更改。
nboot选项用于控制引导示例的数量(默认值为 599)。
案例：隐形衣对人们恶作剧倾向的影响。
样本：80个参与者
实验设计：
前3周：
将80个参与者安置在一个封闭的社区里。社区里到处都是隐藏的摄像机，这样我们就可以记录下恶作剧的行为。我们记录了每个人在最初的3周内做了多少恶作剧(恶作剧1)。
3周后，将被试者分为两组：
斗篷=1（没有斗篷）：样本量为34，在这一组告知他们关闭摄像头，这样没有人能看到他们在做什么。
斗篷=2（有斗篷）：样本量为46，他们每人得到一件隐形斗篷，这些穿着斗篷的人被告知不要告诉任何人他们的斗篷，他们可以随时穿斗篷。
我们记录了在接下来的3周内恶作剧的次数。
变量斗篷记录了一个人是否有斗篷(斗篷= 2)或没有(斗篷= 1)。这些数据在文件CloakofInvisibility.dat中。通过执行下面命令，将工作目录设置为正确的文件夹并执行以下操作，将该文件加载到一个名为invisibilityData的数据流中:
1. 将数值变量cloak转换成一个factor(即类别变量)
代码：
代码含义：（levels=c(1:2)）我们将cloak的两个水平定义为1和2；(labels = c(“No Cloak”, “Cloak”))为这两个水平设置标签叫No Cloak 和 Cloak。
图11.7 invisibility 数据箱线图
2. 图11.7解释：分发cloak前后淘气行为数量的箱形图，这些淘气行为是由是否得到cloak决定。两组的恶作剧水平基本差不多，而且都有所增加(结果如此并不惊讶，因为也告诉了那些没有穿cloak的人相机已关闭)。whiskers显示得到cloak的被试分数差距更大。
3.运行稳健回归的主要困难：将数据转换成正确的格式。图11.8显示了来自invisibilityData dataframe的数据。你可以看到两个组被叠起来，且协变量和因变量（结果）安排在不同的列中。稳健的ANCOVA的函数需要我们创建四个变量：
图11.8 为稳健ANCOVA获取数据
covGrp1: 该变量包含第一组(在本例中为“No Cloak”组的cloak变量)的协变量(mischief1)的得分。见图11.8中左上角数据框。
dvGrp1: 该变量包含第一组因变量/结果(mischief2)的得分(在本例中为“No Cloak”组的cloak变量)。见图11.8中右上角的数据框。
covGrp2: 该变量包含第二组(在本例中为“Cloak”组的cloak变量)的协变量(mischief1)的得分。见图11.8中左下角数据框。
dvGrp2: 该变量包含第二组因变量/结果(mischief2)的得分(在本例中为“Cloak”组的cloak变量)。见图11.8中右下角的数据框。
4. 创建数据具体操作如下：
1) 将数据分成两个dataframe：Cloak组和No-Cloak组，
语句：
含义：创建两个新的dataframe，名叫noCloak 和 invisCloak)；subset()的功能是指定原先的dataframe (invisibilityData)，设置条件来选择行(变量cloak的值在第一个dataframe中为“Cloak”，在第二个dataframe中为“No Cloak”)。
2) 从新的dataframe选择列(即变量)创建４个变量。
语句：
含义：第一行创建一个名为covGrp1的变量，该变量包含invisCloak dataframe中的mischief１变量值；第二行创建一个名为dvGrp1的变量，该变量包含invisCloak dataframe中的mischief２变量值；第三行和第四行命令同上，不过使用的是noCloak dataframe。
3) 将四个变量输入稳健ANCOVA命令中（注意将bootstrap的采样改为2000）并运行：
语句：
5.结果解释：
Output 11.10显示ancova()函数的结果，Output 11.11显示ancboot()的结果。
X列：协变量的5个值(在本例中为2、4、5、6、7)，在这5个值中，两组的baseline mischief和post-cloak mischief之间差不多；
n1和n2：数据中两个组(n1和n2)的case数量，它们的协变量值接近于x(不完全是x，但接近)；
DIF列：基于两个样本，计算trimmed均值(默认为20%)，并检验样本间差异；
se列：两样本估计标准误。
TEST列：检验统计差异（差异除以标准误）；
ci.low ci.hi ：Trimmed均值间的置信区间。
注意，置信区间是两个outputs中第一个不同的地方：这是因为在ancboot()的输出中，置信区间基于bootstrapping；
p值：用于检验trimmed均值间的差异，值小于.05就表明协变量进行调整时，trimmed均值间存在显著差异。
6 .结果的具体含义：
Output 11.10和11.11显示5个design points中有4个trimmed均值间差异显著。即多数情况下，干预后各组的平均恶作剧水平(根据恶作剧的基线水平进行调整)存在显著差异。5附近没有得到协变量值的显著差异 (五个中间design points检验)，似乎表明拥有隐形cloak让那些通常不是很调皮的人(基线分数约2和4)或通常非常调皮的人(基线分数约6和7)恶作剧行为增加，而一般/平均调皮的人并未如此。
7. 用生成的与结果变量对应的协变量图(见图11.9)帮助解释结果：
拟合了两条回归线(每组一条)，但注意不是直线(即斜率非线性)。看两个组在稳健分析中每个design point的数据点分布(即， x = 2,4,5,6,7)。注意，⚪通常比×高。唯一例外是当X = 5时，在最高点有一个×，在最低点有一个⚪。这可能解释了为什么我们在稳健分析的这个design point上没有发现组间的显著差异(在这一点上，⚪比×高得不明显)。
图11.9 ANCOVA函数中基线恶作剧数量（X）相对post-cloak恶作剧（Y）
8. SAM‘s Tips: ANCOVA
协方差分析(ANCOVA)比较几种均值，但根据一个或多个其他变量的效应(称为协变量)进行调整；如，如果你有几个实验条件，想要根据被试的年龄进行调整。
分析前，检查自变量和协变量是否独立。可以使用方差分析或t检验来检查协变量的水平在不同组之间没有显著差异。
要决定是使用Type I 还是 Type III的平方和。如果使用Type III，必须做一个正交的对比而不是非正交。
如果在实验前已经产生了特定假设，请使用planned comparisons。可以使用contrast()函数。
在ANCOVA的输出结果中，查看协变量和自变量的p值。如果该值小于.05，那么对于协变量来说，它意味着该变量与结果变量有显著的关系；对于自变量，意味着在排除协变量对结果的影响后，各实验条件下的均值有显著差异。
如果没有特定的假设，可以使用glht()函数进行事后检验。
对于对比和事后检验，看p值确定比较是否显著(显著性值小于.05表明显著)。
检验与方差分析相同的假设，除此之外，还应检验回归斜率同质性（homogeneity of regression slopes）假设。通过自定义ANCOVA模型检查自变量×协变量交互作用。
11.6. 计算效应量
在前面的章节中我们知道了可以使用η2作为方差分析的效应量测量。这个效应量只是r2的另一个名字，用感兴趣效应(SSM)除以数据中的总方差(SST)计算。因此，它是总变异被一个效应所解释的比例。在ANCOVA(以及我们在之后章节中遇到的一些更复杂的ANOVAs)中，会有不止一种效应；因此，可以计算每种效应的η2。我们也可以用η2（partial eta squared），它与平方的不同之处在于，它不是考察总变异被一个效应所解释的比例，而是考察变量所解释的在分析中其他变量不能解释的变异的比例。举个例子：比如说我们想知道伟哥的药效大小。偏η2是指，伟哥的剂量共享的、不属于伴侣的力比多(协变量) 在力比多中所占变异的比例。协方差不能解释的变异来源有两种：一是无法解释由伟哥剂量引起的变异，即SSViagr；二是无法解释误差变异性SSR。因此，我们在计算中使用这两个变异来源，而不是总变异SST。η2和偏η2的差异如下：
如果用伟哥的例子计算，用Output 11.3的平方和来计算剂量的效应(25.19)，协变量(15.08)和误差(79.05)：
以10.7节中的数据来举例，我们编写了一个名为 rcon-trast() 的函数，得到了t值和df的值。
语法操作如下：
首先，创建一个变量t，包含df ：
将变量 t 和 df 放入 rcontrast() 函数并执行:
得到输出结果：
该结果输出的协变量（.400）、联合剂量组和安慰剂（0.479）之间的效果大小都是中到大的效果，具有显著的统计学差异，而高剂量组和低剂量组之间的差异(.106)是一个相当小的影响。
另一种方法是计算所有组合之间的效应大小，就像方差分析中所做的那样。这需要再次使用 calculate.es包中的 mes()函数:
如何求得调整后的标准差：我们知道调整后的平均值(输出 11.4)和样本大小，但我们不知道调整后的标准差——可以使用未调整的标准差作为近似值，或者以调整后平均值的标准误进行估计(输出 11.4)。我们在第 2 章中发现，标准误差是标准差/平方除以平均值所依据的样本大小的平方根。换句话说，标准差是标准误差乘以样本大小的平方根。
调整后平均值的标准误在语法 adjustedMeans$se 中(参见第11.4.8节)。
首先，创建一个变量n，包含三个组的样本大小：
通过将这个变量(sqrt(n) in R-speak)的平方根乘以相应的标准错误(存储在 adjustedMeans$se 中)来估计标准偏差。执行：
得到：
1.788613 1.755879 1.812267
类似地，如果我们想将低剂量组与安慰剂组进行比较，执行:
我们进入了低剂量组的平均值(5.988117)，安慰剂组的平均值(4.151886)，相应的标准差(1.755879和 1.788613)，以及样本大小(8和 9)。
类似地，我们可以通过执行以下命令来获得高剂量组和安慰剂组之间的效应大小差异:
最后，高剂量组和低剂量组之间的差异可以通过执行以下措施加以量化:
得到输出结果：
11.7 ANCOVA如何报告结果
报告ANCOVA与报告ANOVA非常相似，只是须报告协变量的影响。对于协变量和实验效应，我们给出了F-ratio和计算F-ratio的自由度。在这两种情况下，F-ratio是从均值平方除以均值平方的影响，为残差。因此，用于评估F-ratio的自由度是模型效应的自由度（协变量为dfM=1，实验效应为2）和模型残差的自由度（协变量和实验效应均为dfR=26）——见输出11.3。
正确报告方法：
协变量，伴侣的性欲与参与者的性欲显著相关，F（1，26）=4.96，p＜.05，r=.40。在控制了性欲的影响后，伟哥的剂量对性欲水平也有显著的影响F（2，26）=4.14，p<.05，偏η2=.24。
计划对比显示，与服用安慰剂相比，服用高剂量或低剂量伟哥的性欲明显增强，t(26) = 2.79, p < .01, r = .48；高、低剂量伟哥组间无统计学差异，t(26) = 0.54, p = .59, r = .11。
事后测试可以报告如下(参见输出 11.6):
事后检验显示，高剂量组的协变量调整平均值显著高于安慰剂组(difference= 2.22,t = 2.77, p < .05, d = 1.13)。低剂量组与安慰剂组无显著统计学差异(difference = 1.79,t = 2.10, p = .11, d = 1.04)，低剂量组与高剂量组无显著统计学差异(difference = 0.44,t = 0.54, p = .85, d = 0.11)。尽管低剂量组和安慰剂组之间没有显著性差异，但效果大小值很高。
本章中使用的R数据包：
本章中使用的R函数：
我学到的关键术语：
杨伟文 中国人民大学（北京） 人力资源
刘方圆 乌特勒支大学 发展心理学
李培凯 乌特勒支大学(Utrecht) 工业组织心理学
侯中卫 河北工业大学（天津市） 组织行为学，人力资源管理
赵瑞瑛 阿姆斯特丹自由大学（荷兰阿姆斯特丹） 临床心理学
刘雨桐 哈尔滨师范大学 应用心理学（人力资源方向）
本期排版：孟雪
往期热门：
Play with R 第九期：Comparing two means