R语言 MANOVA测试

多变量方差分析（MANOVA） 是简单的带有多个因变量的方差分析（Analysis of Variance）。它是方差分析的一个延续。在方差分析中，我们通过一个独立分组变量测试一个连续因变量的统计差异。MANOVA是这种分析的延续，它采用多个连续因变量，并将它们共同捆绑成一个加权的线性综合变量。MANOVA比较了新创建的组合是否因自变量的不同水平或组别而不同。人们可以在R编程中非常容易地进行这种MANOVA测试。

例如，让我们进行一个实验，我们给两组大鼠进行两种处理，我们采取的是大鼠的体重和身高。在这种情况下，大鼠的体重和身高是两个因变量，假设两者都受到处理方法的影响。可以用多元方差分析来检验这一假设。

MANOVA的解释

如果全局多变量检验是重要的，那么假设相应的影响是重要的。在这种情况下，随后的问题是要决定治疗是否只影响身高、只影响体重或两者都影响。换句话说，我们要区分哪些特定的因变量促成了显著的全局效应。而为了澄清这个问题，使用单因素方差分析来分别检验每个因变量。

MANOVA的假设

MANOVA可以在特定条件下使用。

因变量在组内应该是正态分布。
预测因素范围内变异的均匀性。
所有的协变量对、所有的因变量对、以及每个单元中所有的因变量-协变量对之间的线性关系。

R语言的实现

R提供了一个方法 manova() 来执行MANOVA检验。类 “manova “与类 “aov “的不同之处在于选择了不同的总结方法。函数 manova() 调用aov，然后将类 “manova “添加到每个分层的结果对象中。

语法：

manova(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, …)

参数：

formula： 一个指定模型的公式。
data ： 一个数据框，在其中可以找到公式中指定的变量。如果缺失，将以标准方式搜索变量。

projections： 逻辑标志

qr： 逻辑标志

contrasts： 用于公式中某些因素的对比列表。

…：传递给lm的参数，如subset或na.action

例子：

为了在R中进行MANOVA测试，让我们采取虹膜数据集。

# R program to illustrate
# MANOVA test
  
# Import required library
library(dplyr)
  
# Taking iris data set
myData = iris
  
# Show a random sample
set.seed(1234)
dplyr::sample_n(myData, 10)

输出

 Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1           5.5         2.5          4.0         1.3 versicolor
2           5.6         2.5          3.9         1.1 versicolor
3           6.0         2.9          4.5         1.5 versicolor
4           6.4         3.2          5.3         2.3  virginica
5           4.3         3.0          1.1         0.1     setosa
6           7.2         3.2          6.0         1.8  virginica
7           5.9         3.0          4.2         1.5 versicolor
8           4.6         3.1          1.5         0.2     setosa
9           7.9         3.8          6.4         2.0  virginica
10          5.1         3.4          1.5         0.2     setosa

为了了解不同物种之间在萼片和花瓣的长度上是否存在任何重要的差异，可以进行MANOVA测试。因此，可以使用函数 manova() ，如下所示。

# Taking two dependent variable
sepal = irisSepal.Length
petal = irisPetal.Length
  
# MANOVA test
result = manova(cbind(Sepal.Length, Petal.Length) ~ Species, 
                                                data = iris)
summary(result)

输出

            Df Pillai approx F num Df den Df    Pr(>F)    
Species     2 0.9885   71.829      4    294 < 2.2e-16 ***
Residuals 147                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1