Play with R 第16期：多元方差分析（MANOVA）

16.1 这一章会告诉我什么？

由养宠物有利于心理健康，引出心理健康一词涵盖了广泛的概念，包括（仅举几个例子）焦虑、抑郁、一般痛苦和精神病。当我们想在几个结果变量上比较时，我们可以将方差分析扩展到MANOVA。这就是这一章的主题。

16.2 什么时候使用MANOVA？

MANOVA可以被认为是对于多个因变量的情况的方差分析。方差分析的原理可以扩展到MANOVA，因为当只有一个自变量时，我们可以使用MANOVA，或者当有几个自变量时，我们可以看到自变量之间的相互作用，甚至可以对比组间的差异。方差分析只能在有一个因变量（或结果）的情况下使用，因此被称为单变量检验（单变量很明显意味着“一个变量”）；MANOVA是可以同时考查几个因变量（结果），是一个多元检验（多元意味着“许多变量”）。

16.3 导言：与方差分析的异同

如果对于多个因变量使用方差分析，那么所测量的因变量越多，需要进行的方差分析就越多，而且发生I型错误的可能性就越大。还有一些原因使用MANOVA。一方面，MANOVA可以得到重要的附加信息，而对每个因变量进行单独的方差分析就会忽视因变量之间的关系。另一方面，方差分析只能告诉我们一个维度上的组别差异，而MANOVA可以考察一个维度的组合上的组别差异。因此，MANOVA更具有检验效力。

16.3.1 需要注意的方面

除非有很好的理论或经验基础，否则把所有的因变量集中在一个MANOVA中不是一个好主意。

16.3.2 本章的例子

考察认知行为疗法（CBT）对强迫症（OCD）的影响。我们可以将在CBT之后和行为治疗（BT）之后的OCD患者和一组仍在等待治疗的OCD患者（无治疗条件，NT）进行比较。对于强迫症的治疗，不仅要看行为结果（Actions），也要看认知结果（Thoughts），也就是有两个因变量。

16.4 MANOVA理论

16.4.1 矩阵介绍

详见p 701.

16.4.2. 一些重要的矩阵及其函数

表示系统方差的矩阵（或所有变量的平方和模型）用字母H表示，称为平方和和交叉乘积矩阵的假设（或假设SSCP）。表示非系统方差的矩阵（所有变量的残差平方和）用字母E表示，称为平方和和交叉乘积矩阵的误差和（或误差SSCP）。最后，有一个矩阵表示每个因变量的方差总量（每个因变量的平方和总和），用T表示的，是称为平方和和交叉积矩阵的总和（或总SSCP）。

之后将展示这些矩阵如何与方差分析中的简单平方和（SSM、SSR和SST）完全相同地使用，推导出一个表示模型中系统方差与非系统方差之比的检验统计量。本文中描述的矩阵都被称为平方和交叉乘积（SSCP）矩阵。

16.4.3 手动计算方差分析：一个案例

16.4.3.1 DV1单因素方差分析（Actions）

计算三种平方和，SST, SSM, SSR。计算F值。

16.4.3.2 DV2单因素方差分析（Thoughts）

同上，计算三种平方和及F值。

16.4.3.3 DVs之间的关系：交叉积

有三种不同的交叉积，这三种交叉积与我们计算的一元方差分析的三个平方和有关：也就是说，有一个总交叉积，一个由模型而产生的交叉积，和残差的交叉积。

总交叉积（CPT）、由模型而产生的交叉积（CPM）、残差交叉积（CPR）的计算见p 705-707。

16.4.3.4 SSCP总矩阵（T）

在这个例子中只有两个因变量，所以所有的SSCP矩阵都是2×2矩阵。总SSCP矩阵（T）包含每个因变量的平方和和两个因变量之间的总交叉乘积。

16.4.3.5 残差SSCP矩阵（E）

类似于SSCP总矩阵。

16.4.3.6 模型SSCP矩阵（H）

当计算单变量方差分析时，我们看到平方的总和是模型的平方和和残差的平方和之和（即SST = SSM + SSR）。MANOVA也是类似的道理，是矩阵之和相加。

16.4.4 MANOVA检验统计原理

在单因素方差分析中，我们计算了系统方差与非系统方差的比值（即SSM除以SSR）。因此同样道理，将矩阵H除以矩阵E。然而，有一个问题，即矩阵不能被其他矩阵整除。然而，有一个等价于除数的矩阵，也就是乘以所谓的矩阵的逆。因此，如果想用H除以E，我们必须用E的逆（表示为E-1）和H相乘。因此，检验统计量是基于模型SSCP与残差SSCP的逆相乘得到的矩阵。这个矩阵称为HE-1。

HE-1表示模型中的系统方差与模型中的非系统方差之比，因此得到的矩阵在概念上与单因素方差分析中的F比值相同。

16.4.4.1 判别函数变量

可以计算因变量的潜在线性维度。这些因变量的线性组合称为变量（或有时称为潜在变量或因素）。在这种情况下，我们希望使用这些线性变量来预测一个人属于哪一组（即，他们是否接受CBT、BT或不治疗），所以我们可以用以辨别不同组别的个体。因此，这些变量称为判别函数或判别函数变量。

16.4.4.2 Pillai-Bartlett trace（V）

Pillai-Bartlett trace是解释方差在判别函数上所占比例的总和。因此，它类似于SSM/SST的比率，即R2。

16.4.4.3 Hotelling's T2

Hotelling-Lawley trace（也称为Hotelling's T2），只是每个变量的特征值之和。这个检验统计量是每个变量的SSM/SSR之和，因此它直接与方差分析中的F比值进行比较。

16.4.4.4 Wilks’s lambda (Λ)

Wilks’s lambda是每个变量上无法解释的方差的乘积。Wilks’s lambda表示每个变量的误差方差与总方差（SS R/SS T）之比。

16.4.4.5. Roy’s largest root

Roy’s largest root是第一个变量的特征值。因此，在某种意义上，它与Hotelling-Lawley trace相同，但仅适用于第一个变量。因此，Roy’s largest root表示第一判别函数的解释方差与未解释方差（SSM/SSR）的比例。

16.5 Practical issues when conducting MANOVA

在运行MANOVA之前，需要考虑以下三个方面。

16.5.1 Assumptions and how to check them

首先要考虑的是有关MANOVA统计检验的假设前提是否满足。MANOVA的统计前提与先前介绍过的ANOVA类似，但略有不同:

* Independence（相关性）: 结果变量之间须存在统计意义上的相关性。

* Random Sampling（随机抽样）: 数据采样随机。

* Multivariate Normality(多变量分布正态): 不同组被试间，多个结果变量的数据分布呈现正态。

* Homogeneity of covariance matrices(协方差齐性): 不同的结果变量下的组内方差不存在显著差异，且任意两两结果变量之间的协方差也不存在显著差异。即对于不同组的被试来说，“方差-协方差矩阵”中的值相等。

其中，关于多个结果变量的分布正态检验，可使用R中的Shapiro test，也可用mvoutlier包中的aq.plot()函数画图来查看。

其中，关于协方差矩阵的齐性检验，可使用Box's test，如果满足假设前提的话，则检验结果应为不显著。但是Box's test的结果不是很稳定，易受到样本量的影响，如果样本量较大，则即便协方差相差不多，结果也有可能显著。因此，如果不同实验组之间的样本量相同的话，一般可不用参考Box's test的结果，因为一方面它不稳定，另一方面后续的Hotelling's和Pillai‘s统计检验的结果更可靠。

16.5.2 Choosing a test statistic

其次要考虑的是在MANOVA的显著性检验中，存在4种不同的统计检验值方法，具体应该如何选择统计检验值？本文从研究方法的统计效率和研究方法的可靠性（当统计假设前提不满足时，统计结果是否可靠）两个方面做了评估。

从统计效力来看, Olson(1974)研究发现，对于相对较小或中等的样本量来说，4种统计检验方式并无统计效力上的显著差异。如果预测变量的不同水平组之间的差异集中在多个结果变量中的一个上的话（通常很多社会科学研究只会关注一个结果变量指标），则Roy's statistic test是最有统计效力的，其次是Hotelling's trace, 再次是Wilk's lambda, 最后是Pillai's trace。但如果预测变量的不同水平组之间的差异集中在多个结果变量上的话，则统计效力的排序和上述情况完全相反。不同方法的统计效力还受到样本量和结果变量数目的影响。因此Stevens(1980)建议如果样本量不够大的话，则结果变量的数目不要超过10个。

从统计方法的可靠性来看，4种统计检验方法都能比较好地对抗统计假设前提不满足的情况。不过，其中Roy's root会受到platykurtic分布的影响，此外它在协方差齐性这一统计假设不满足的情况下，可靠性也会降低。Olson和Stevens的研究均表明，如果样本量相同的话，则Pillai's trace是最可靠的统计检验方法。但是如果预测变量的不同水平组之间的样本量不同的话，则需要考虑协方差矩阵的同质性和多个预测变量的正态分布，如果这两个假设满足的话，Pillai's trace方法才可靠。

 16.5.3 Follow-up analysis

第三个要考虑的因素是随后检验方法，即在运行完MANOVA的统计检验测试之后，后续应该如何分析数据。主要有两种后续分析方式。

一般情况下，建议MANOVA分析结束后，针对不同的结果变量依次进行后续的ANOVA分析。但是如果用单独分开的ANOVA分析的话，则局限在于无法整合多个结果变量来解释最终呈现出来的组间差异。

因此，还有学者建议MANOVA分析结束后，使用Discriminant Analysis方法，因为此法可以综合考虑多个结果变量，并且找到能最大化反映组间差异的结果变量权重组合。这个方法与MANOVA更为搭配，能把整组的结果变量都纳入考量。

16.6 MANOVA using R

此处将沿用16.1中的例子，用OCD.dat来进行分析。其中，预测变量（或自变量）有三个不同的水平: CBT组，BT组，NT组；结果变量有两个: 行为（Actions），认知（Thoughts）。

16.6.1 Packages for factorial ANOVA in R 

运行MANOVA需要用到以下的包，其中除了MASS已经自动安装外，其他的都要手动安装。安装指令如下:

各个包的用途介绍如下：

16.6.2 General procedure for MANOVA

运行MANOVA包含以下几个步骤：

1. Enter data (数据输入)

2. Explore your data (数据探索): 包括用图形来呈现原始数据；进行描述性统计分析；检验MANOVA的两个统计前提（多个结果变量的分布正态性，方差-协方差矩阵的齐性检验）。

3. Set contrasts for all predictor variables (设置用于显著性差异比较的预测变量水平)

4. Compute the MANOVA (运行MANOVA)

5. Run univariate ANOVAs (运行随后的分析 -- 选择项)

6. Discrimination function analysis (运行随后的分析 -- 选择项)

后续的16.6.4 - 16.6.9即分别介绍了这些步骤

16.6.3 MANOVA using R Commander

16.6.4 Enter the data

先读取wide格式的数据文件OCD.dat:

数据如下:

16.6.5 Exploring the data

通过散点图，可以知晓两个不同的结果变量之间的相关关系：

通过柱状图或箱图，可以知晓不同实验组的两个不同结果变量的均值、标准误、分布等信息:

通过by()和stat.desc()函数,可以查看不同实验组在不同结果变量上的描述性统计数据:

通过by()和cov()函数，可以对协方差进行齐性检验。如下图，分不同的实验组，来对位于第2、3列的Actions和Thoughts数据指标进行协方差的计算：

进而得到方差-协方差矩阵，如下图：

从上图中，我们可以看到不同组之间的方差和协方差是存在一定差别的（前提不满足）。但是，由于这里不同组之间的样本量相等，因此我们可以忽略方差不齐的情况，继续后续的MANOVA分析。而如果不同组之间的样本量不同，那么就要分情况讨论：1） 如果样本量越大的那个实验组的方差和协方差也越大的话，则后续的统计检验相对比较保守，因此后续的统计检验的显著性结果是相对可信的。2）如果样本量越小的那个实验组的方差和协方差却越大的话，则后续的统计检验相对比较宽松，因此后续的统计检验呈现出的显著结果需要谨慎对待。

通过mshapiro.test()函数，可以对多个结果变量的分布正态性进行检验。步骤如下，首先提取不同实验组的得分数据:

由于mshapiro.test()只能识别按行排列的数据，因此要对数据进行转置：

之后进行正态分布检验, 如果p值小于0.05的话，就代表该组数据偏离了正态分布，违背了多变量正态分布这一前提要求：

还可以用aq.plot()函数画图，进而找到极端值：

上图的数据表示的是每个数据点的行号编码（一个被试的数据即一行），在上图中，除了右上的图以外其他三幅图中用红色呈现的数据即代表该数据是极端值（上图是黑白版本，没有显示红色，但26原本是红色，说明编号26的被试是一个极端值）。而在右上图中，在97.5%的分割线右边的数据，即代表极端值（此处还是显示为26）。综上说明，26号这个被试的数据可能是个极端值，因此需要考虑是否将被试剔除，进而查看剩下来的被试数据是否符合多变量正态前提；或者考虑是否留下这个被试，看后续的统计显著性检验是否能消除这个极端被试的影响，依旧显著。

16.6.6 Setting Contrasts

在这一步骤中，需要设置用来比对的数据组，如下:

其中，contrasts(ocdData$Group)代表将在Group这个自变量的不同水平之间设置比对数据组。contr.treatment代表本次数据比对是处理组之间的比对，base=3代表将第三组（NT)设置为基线对照组。

一般建议在比对数据的时候，手动进行定义，这样命名方便，一目了然，如“CBT_vs_NT”代表CBT组和NT进行比对：

16.6.7 The MANOVA model

我们使用manova()函数进行正式的MANOVA分析，

其中：

*newModel代表输出结果，通过后续执行summary(newModel)命令，即可查看MANOVA的分析结果。

*outcome代表MANOVA测试中的被预测变量（即结果变量），在本文例子中指的是Actions和Thoughts这两个结果变量。

*predictor(s)类似于因变量，在本文例子中指的是Group变量。

*dataFrame指的是统计检验中预测变量和结果变量的数据来源。

*na.action是一个备选命令，如果数据完整无缺的话则可遵循默认值；如果数据有缺失值的话，则设置为na.action = na.exclude（剔除缺失值）。

因此，在本文的例子中，我们先对outcome赋值：

然后执行分析：

然后选择前述4种检验方法中的一种来进行统计显著性检验，默认的检验方法是Pillai's trace, 想用其它方法的话可以自定义：

结果如下：

在上述结果中，Pillai’s trace (p = .049),Wilks’s lambda (p = .050), Roy’s largest root (p = .020)都达到了验证实验组组间差异为显著的标准，而Hotelling’s trace (p = .051) 没有达到判定实验组组间差异显著的标准。但由于本文例子中不同实验处理组的样本量一致，因此我们选择相信统计效力较强的Pillai’s trace，认为这是一个显著的结果。

MANOVA的结果只能告诉我们不同实验处理组在多个结果变量指标上存在差异，但是不能告诉我们这种差异究竟是谁大谁小，也不能告诉我们这种差异究竟体现在哪一个结果变量上。因此，我们要进行进一步的单变量检验分析。

16.6.8 Follow-up analysis: univariate test statistics

如果想在MANOVA的基础上进行进一步的方差分析的话，可执行以下命令：

结果如下：

上述结果可以发现，如果按照单个的结果变量来分开做方差分析（ANOVA）的话，实验处理之间的差异是不显著的。这与MANOVA分析得出的显著结果完全不同。这是因为，MANOVA把结果变量之间的相关性纳入了考量，因此数据能够更为敏感地捕捉到组间差异。这也说明在本例中，不同的处理组之间的差异显现在综合的多个结果变量指标上，而不是单个的结果变量指标。

除了进一步做ANOVA之外，如果想要查看不同的结果变量是如何交互地体现实验处理组之间差异的话，可进行进一步的discriminant function analysis分析（在接下来的章节中会有介绍）。

16.6.9 Contrasts

注意，如果先前分不同的结果变量做的univariate ANOVAs得出显著性结果的话，则可以通过线性建模来查看不同处理组之间的具体差异模式。这一方法类似于单因素方差分析。可用lm()函数来实现：

得出的结果如下：

16.7 鲁棒（稳健）的多元方差分析

Wilcox为MANOVA提供了两种鲁棒方法，这两种方法均基于数据排序（见15章）。我们需要加载WRS package，使用其中的两个函数：

mulrank()：对排序数据执行MANOVA

cmanova()：对排序数据执行稳健性检验（robust test）

这两个函数只能在仅有一个自变量时使用。

将原始数据（长数据）整理成宽数据，像下图

1.设置路径，导入数据并赋值为OcdData：

OcdData<-read.delim(“the.dat”,header=TRUE)

2.数据框增加变量row，1:10重复三次是因为三个治疗组都是包含10行数据

ocdData$row<-rep(1:10, 3)

3.使用melt()生成新的数据框

ocdMelt<-melt(ocdData, id = c("Group", "row"), measured = c("Actions","Thoughts"))

4.使用names()命名

names(ocdMelt)<-c("Group", "row", "Outcome_Measure", "Frequency")

5.使用cast()转换数据框格式

ocdRobust<-cast(ocdMelt, row ~ Group + Outcome_Measure, value =

"Frequency")

6.清除变量row

ocdRobust$row<-NULL

7.输出最终数据

ocdRobust

mulrank()和cmanova()的使用格式：

mulrank(number of groups, number of outcome measures, data)

cmanova(number of groups, number of outcome measures, data)

根据ocdRobust，执行下列操作：

mulrank(3, 2, ocdRobust)

cmanova(3, 2, ocdRobust)

结果为：

左侧mulrank()的输出告诉我们：

F （$test.stat）= 1.64，p（$p.value ）= .168：治疗方式对强迫症预后的主效应不显著。

$q.hat可体现相对效应，将其重新标志为：

结果为：NT组中Actions和Thoughts的等级（0.57和0.54）相似，BT组中Actions（0.38）的等级低于Thoughts（0.59），CBT组中则是Thoughts（0.37）的等级低于Actions（0.55）。这说明CBT疗法对想法的影响大于行为，BT疗法影响相反。但是，总体效应不显著。

右侧cmanova()的输出告诉我们类似的结果：

H(4)=9.06，p=.060：治疗方式对强迫症预后的主效应不显著。

16.8 报告多元方差分析的结果

作者认为，报告MANOVA很像报告ANOVA（方差分析）：给出F和df。因此，可报告为：

治疗对强迫想法和行为的数量有显著影响，F(4, 54) = 2.56, p < .05。

作者指出，多元检验统计量也应该被引用，Output16.5的四个测验（实际中报告其中一个即可）依次被报告为：

我们可用常规方式报告随后的ANOVAs（见Output 16.6）：

若使用了robust MANOVA(Output 16.9)可如下报告：

16.9多元方差分析后的判别分析

显著的MANOVA后可进行单变量方差分析和判别分析（有时被称为判别函数分析或简称为DFA）。本章例子因变量之间显著相关，所以不能用单变量方差分析。判别分析是调查因变量关系性质的极佳方法。作者强烈建议，如果想充分理解数据，那么应该在MANOVA后进行单变量检验和判别分析。

在DFA（判别分析）中，我们观察如何使用多个预测器来最好地区分各组（它有点像逻辑回归logistic regression，但它有多个组而非2组）。某种角度上，这似乎是在做与MANOVA相反的事情：MANOVA是从一组变量来预测一组结果，而DFA是从一组结果来预测一组变量。然而，这些检验的基本原则是相同的：当我们研究MANOVA理论时，我们看到它是通过识别线性变量来工作的，这些线性变量是DFA中的函数。

DFA在R中非常容易实现，使用MASS package的lda()函数，函数的使用格式为：

newModel<-lda(Group ~ Predictor(s), data = dataFrame, prior = prior probabilities, na.action = "na.omit")

也可以选用其它的函数执行，但lda可获得更多的信息，这是此时我们真正需要的。在lda()函数中，Group是dataframe中变量名，该变量包含要区分的组，而Predictor(s)是要进行区分的连续变量列表，这就为线性模型创建了一个公式。如果使用单个预测器，是Group ~ Predictor；使用多个预测器，则用”+”连接，如Group ~ Predictor1+ Predictor2。根据本章例子，Group ~ Actions+ Thoughts。Data选项制定数据名称（本例为ocdData），no.action选项决定如何处理丢失的数据

根据当前数据，执行：

ocdDFA<-lda(Group ~ Actions + Thoughts, data = ocdData)

查看模型，执行：

ocdDFA

1.  Prior probabilities of groups:先验概率，每组均是0.33，这些是组样本大小除以总样本大小（即10/30=.33）。

2.  Group means:组平均值

3.  Coefficients of the linear discriminants：线性判别公式的系数，即公式中的b值，告诉我们每个变量的相对贡献。

变量1中想法和行为有相反的作用（变量1与行为正相关，变量1与想法负相关），因此变量1可被看作是一个区分想法和行为的变量（它以相反的方式影响想法和行为）同理，变量2以相同的方式影响想法和行为

4.  Proportion of trace:变量1占变异的82.2%，变量2仅占17.8%。这些比例是按比例表示的每个变量的特征值（the values of the diagonal elements of the matrix HE−1）。

判别分数是每个人在各个变量上的得分，可能代表潜在的社会或心理结构，为此执行：

predict(ocdDFA)

Output16.11显示每个参与者在LD1和LD2这两个变量上的得分

比分数本身更有用的是按组成员分列的分数图。这可以通过使用我们模型上的 plot()函数来获得，执行：

plot(ocdDFA)

得到上图顶部的图，这个图表绘制出每个人的变量得分，并根据这个人所属的实验条件进行分组。为解释这个图，作者将其分解变量1图和变量2图。

我们关注LD1时，忽略LD2，即忽略图中每个点的垂直位置，看这些组是如何沿着水平轴(LD1)分布的。一个简单粗暴的方法是将垂直轴沿中间分割(图中垂直蓝线)，然后看各组如何分布在垂直轴的两侧。用浅蓝色圈出 BT 组，用黑色圈出 CBT 组。希望这能够清楚地表明，在蓝色垂直线的左边有很多来自 BT 组的病例，但是在蓝色垂直线的右边有很多来自 CBT 组的病例。换句话说，蓝色的垂直线似乎把 BT 组和 CBT 组分开了。这告诉我们变量1区别了 BT 组和 CBT组。

我们关注LD2时，忽略LD1，即忽略图中每个点的水平位置，看这些组是如何沿着水平轴(LD1)分布的。通过垂直轴（图中中点0点处水平线），然后看各组如何分布在垂直轴的两侧。

16.10 报告判别分析的结果

呈现数据的指导原则是为读者提供足够的信息，使他们能够自己判断数据的含义。作者建议报告解释的方差百分比和线性判别系数，可如下报告：

16.11  最后的说明

综合我们所做的数据分析，进行结论的说明。

16.11.2  单变量方差分析还是判别分析？

单变量方差分析和判别分析是回答一个显著的MANOVA所引发的不同问题的方法。作者建议运行两种分析以全面了解数据。判别分析的优势在于，它能告诉你一些关于数据中潜在维度的信息(如果你采用了几个独立的测量方法，试图得到一些社会或心理结构，那么它尤其有用)。作者强烈建议即使单变量ANOVAs是重要的，判别分析提供了对数据的有用的洞察力，应该被使用。

What have I discovered about statistics?

测量多个组的多个结果的情况下，ANOVA可以扩展为MANOVA。相比于运行多个ANOVA，MANOVA可以保留对类型I错误的控制，并且可将结果变量之间的关系合并到分析中。ANOVA和MANOVA工作方式相似，只是前者用的是单值而后者用的时矩阵。

MANOVA会提供四个与相同效果相关的统计量，作者认为Pillai’s trace是最安全的选择。

跟踪MANOVA有两种方法：运行大量ANOVA或进行判别分析，判别分析会提供更多的信息，但解释很麻烦。

本章使用的R数据包

本章使用的R函数

关键术语

本期特别负责的小伙伴~

16.1-16.4 (p696-716) 代嘉幸，辽宁师范大学，研究兴趣：社会认知，人格心理学

16.5-16.6 (p717-p733)  施赛男，华东师范大学，研究兴趣：情绪，经验采样，多水平建模

16.7-16.11 (p733 –p745) 张莹莹，东北师范大学，研究兴趣：社会认知，决策，亲社会行为

本期排版：

杨伟文

往期热门：