R统计绘图-One-Way MANOVA

存在两个及以上连续结果(或响应)变量的ANOVA被称为多变量方差分析(Multivariate Analysis of Variance,MANOVA)。例如，将小鼠分为处理A和处理B两组后，测量小鼠的长度和高度。此时小鼠的长度和高度就是结果变量，假设长度和高度都因为处理不同产生差异。此时可以使用MANOVA检测上述假设。

MANOVA分析过程总结如下：

1)将所有结果变量经过线性组合形成一个新的复合变量；

2)比较新变量在分类变量中的均值。

本文接下来描述在R中如何进行MANOVA。

一、设置工作路径并调用R包

# 设置工作路径
#knitr::opts_knit$set(root.dir="D:\\EnvStat")# 使用Rmarkdown进行程序运行
Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置# 调用R包
library(tidyverse)
library(ggpubr)
library(rstatix)
library(RColorBrewer)
library(ggsci)
library(car)
library(GGally)
install.packages("plotrix")
library(plotrix) # 只能绘制截断图
library(ggbreak)
library(ggsignif)

二、数据准备

options(stringsAsFactors=F)# R中环境变量设置，防止字符型变量转换为因子
# 读入环境因子数据表
ENV=read.csv("env.csv",header = T,row.names = 1,sep = ",",comment.char = "",stringsAsFactors = F,colClasses = c(rep("character",4),rep("numeric",11)))
head(ENV) # 查看数据前几行
#dim(ENV) # 查看数据行、列数
#str(ENV) # 查看数据表每列的数据形式# 描述统计数据
ENV %>% group_by(condition) %>% get_summary_stats(names(ENV[4:14]),type="mean_se")# 可视化数据
## 使用ggplot2绘制，使用ggsci中的jco配色
pdf("Visualization environmental factor.pdf",width =20,height = 15,family = "Times" )
ENV %>% gather(key="variables",value = "value",names(ENV[4:14])) %>%ggplot(aes(x=depth,y=value,color=tillage))+stat_boxplot(geom = "errorbar",position = position_dodge(1))+geom_boxplot(position = position_dodge(1))+geom_point(position = position_dodge(1))+theme_bw(base_size = 16)+scale_color_jco()+facet_wrap(~variables,scales = "free")
dev.off() # 可以设置Y轴标签随数据变化## 使用ggboxplot绘制箱线图
pdf("Visualization environmental factor.pdf",width =17,height = 13,family = "Times" )
ggboxplot(ENV,x = "depth",y =names(ENV[4:14]),combine = TRUE,fill = "tillage",palette = "jco",add = "point",error.plot = "upper_errorbar",font.label = list(size = 16, face = "plain", color ="black")) # palette参数用于设置颜色板
dev.off() # Y轴标签都是固定的，还没找到怎么调整。

注：使用ggboxplot()绘制多Y箱线图时，要注意combine和merge两个参数。combine = TRUE，多Y将绘制在不同的绘图区域，merge=TRUE,多Y将绘制在同一绘图区域。还是推荐使用ggplot2绘图。

图1|环境因子及分组信息表，env.csv。行为样品名称，列为环境因子名称和分组信息，共有11个环境变量，3个分组信息。

图2|环境因子数据可视化，分面箱线图。以土壤深度为横坐标，环境因子数值为纵坐标，不同tillage进行着色，以环境因子变量分分面。

三、单因素MANVOA

进行MANVOA分析需要满足一下假设：

1)足够的样本大小，一般每个分组的n大于结果变量的数目。

2)独立观察：每个观察值应该仅只属于一个组。各组的观察结果之间没有关系，不允许同一样本进行重复测量。样本应该是随机选择的。

3)不要出现单变量或多变量离群值。

4)多变量满足正态分布：rstatix包中的mshapiro_test()可用与进行多变量正态性分布Shapiro-Wilk检验。

5)结果变量之间不能出现多重共线性：结果变量之间不能太相关，相关系数不应超过0.9。

6)各组结果变量之间呈线性关系。

7)满足方差齐性(同质性)：Levene检验用于检测各组间方差的同质性，Levene检验结果无显著性则表明各组间方差具有同质性。

8)满足方差-协方差矩阵同质性：Box's M检验用与检测各组间协方差的同质性，相当于多元方差的同质性检验。Box's 检验高度敏感，alpha=0.001时确定检验的统计学意义。

3.1 单因素MANOVA假设检验

ENV数据中存在两个分类变量tillage和depth，这里是进行单因素MANOVA，所以下面的分析都忽略depth变量。在前面推文中介绍过的假设检验这里就不多缀述了。可以查看前面的推文：R统计绘图-协方差分析[Translation]。

# 3.1.1 检查样本量是否足够
ENV %>% group_by(tillage) %>% summarise(N = n())
## 上面统计结果显示tillage每个组中只有9个样本，而环境变量有11个，不满足假设。
## 因此对环境变量进行筛选，只保留氮相关的环境因子
env = ENV[c(2,4,7:9)]
set.seed(123)
env %>% sample_n_by(tillage,size = 2) # 基于tillage分组显示数据，每组随机显示两个样本。# 3.1.2 检测离群值
## 检测单变量TN的各组离群值
env %>% group_by(tillage) %>%identify_outliers(TN) # 没有## 检测多变量离群值：多变量离群值指在结果变量上具有异常结合值得数据点。
## Mahalanobis distance经常用于在MANOVA过程中检测多变量离群值。这个检验会检测观察值距离均值有多远，还会考虑到协方差。
## rstatix包的mahalanobis_distance()用于计算Mahalanobis distance并标注多变量离群值。
env %>% group_by(tillage) %>% mahalanobis_distance(names(env[2:5])) %>%filter(is.outlier == TRUE) %>%as.data.frame() # 根据Mahalanobis距离(p>0.001)，数据中没有多变量异常值。# 3.1.3 正态分布假设检验
## 单变量正态分布假设检验
env %>% group_by(tillage) %>%shapiro_test(TN,Ammonia,Nitrate,AHN) %>%arrange(variable) # Shapiro检验，p<0.05,该分组数据不满足正态分布假设。
?ggqqplot
ggqqplot(env,names(env[2:5]),facet.by = "tillage",ggtheme = theme_bw()) # 四个变量会产生四幅组合QQ图。样本量大于50时推荐使用QQ图检测正态分布假设。
## 多变量正态分布假设
env %>% select(names(env[2:5])) %>%mshapiro_test() # p<0.05,多变量正态分布假设不满足。可以选择其他非参数检验方法。# 3.1.4 检测多重共线性
## 结果变量间的相关性应该是居中的，如果相关性大于0.9表明有多重共线性，不适合使用MANOVA，可以考虑删除一个高度相关的结果变量。如果相关性值太低，则需要考虑对每个结果变量进行单因素ANOVA分析。
#env %>% cor_test(names(env[2:5]),method = "pearson") # 适合只有两个结果变量时使用
env %>% cor_mat(names(env[2:5]),method = "pearson") # 没有两个结果变量间相关性系数大于0.9。# 3.1.5 检测线性假设，每组的结果变量应该是线性关系。
?ggpairs
results <- env %>% select(names(env[c(1:5)])) %>%group_by(tillage) %>% doo(~ggpairs(.) + theme_bw(),result = "plots")
results$plots[1] # 各结果变量在各分组中没有检测到线性关系，可以转换或删除有问题的结果变量。# 3.1.6 检测协方差同质性假设
box_m(env[2:5],env$tillage) #p<0.001,数据违反了方差-协方差同质性假设。
## 如果是均衡实验设计，每个分组具有相同的样本大小，违背方差-协方差同质性假设，仍然可以继续进行分析。
## 如果实验设计不均衡，则可以1)转换变量数据以满足此假设或2)使用Pillai多变量统计值代替Wilks统计值。# 3.1.7 检测方差同质性假设，检测每个变量在各分组中是否具有方差齐性。
env %>% gather(key = "variable",value = "value",names(env[2:5])) %>%group_by(variable) %>% levene_test(value ~ tillage) # p<0.05，不满足方差齐性。

注：马氏距离(Mahalanobis distance)计算时要求总体样本数大于样本的维数，否则得到的总体样本协方差矩阵逆矩阵不存在，这种情况下使用欧式距离计算即可。样本点在二维空间中共线，比如(3,4)、(5,6)和(7,8)，这种情况协方差矩阵的逆方差仍然不存在，同样采用欧式距离进行计算。

图3|环境因子变量多重共线性检验。结果变量间的相关性应该是居中的，如果相关性大于0.9表明有多重共线性，不适合使用MANOVA，可以考虑删除一个高度相关的结果变量。如果相关性值太低，则需要考虑对每个结果变量进行单因素ANOVA分析。

图4|检测线性假设。每组的结果变量应该是线性关系。如果各结果变量在各分组中没有检测到线性关系，可以转换或删除有问题的结果变量。

3.2 单因素MANOVA计算

这里介绍4种MANOVA分析统计值：Pillai, Wilks, Hotelling-Lawley和Roy.Wilks’ Lambda是比较常用的方法。Pillai’s Trace结果更加稳健，推荐实验设计不均衡或不满足方差-协方差同质性假设时使用。car包中的Manova()默认使用pillai统计。

model <- lm(cbind(TN,Ammonia,Nitrate,AHN) ~ tillage,env)
model
Manova(model,test.statistic="Pillai") # p>0.05,环境因子结合变量处理间差异不具统计学意义。

图5|单因素MANOVA结果。 p>0.05,环境因子结合变量处理间差异不具统计学意义。

3.3 事后检验

为了检测每个变量对显著性结果的贡献，MANOVA分析完成后，接着进行单变量单因素方差分析。

rstatix包的anova_test()用于对满足正态分布和方差齐性假设的变量进行方差分析。变量不满足方差齐性检测则使用welch_anova_test()。非参数kruskal_test()用于对不满足正态分布或方差齐性假设的变量进行方差分析

3.3.1 单变量单因素方差分析

# 将数据由wide format转换为long format
data <- env %>% gather(key = "variables",value = "value",names(env[2:5])) %>%group_by(variables) # 此处变量名一定要命名为variables，否则后续添加显著性标签位置时，会报错。
data# anova_test
data %>% anova_test(value ~ tillage) %>%adjust_pvalue(method = "bonferroni")# welch anova test
data %>% welch_anova_test(value ~ tillage) %>%adjust_pvalue(method = "bonferroni")# Kruskal test
data %>% kruskal_test(value ~ tillage) %>%adjust_pvalue(method = "bonferroni")

图6|单因素非参数 Kruskal test结果。经bonferroni校正后的P值>0.05,接受零假设，表示不同处理间环境因子的含量差异，不具统计学意义。

3.3.2 两两比较

单因素方差分析完成后，对变量进行两两比较，检测存在显著组间差异的两组。当数据满足方差齐性假设时，rstatix包的tukey_hsd()用于计算Tukey post-hoc检验。数据不满足方差齐性假设时，则使用games_howell_test()进行Games-Howell post-hoc检验。 pairwise_t_test(pool.sd=FALSE, var.equal=FALSE)也可用于对不满足方差齐性假设的数据进行两两比较。


res <- data %>% games_howell_test(value ~ tillage) %>%select(-estimate,-conf.low,-conf.high) #移除一些细节
res

图7|games_howell_test两两比较结果。只有Nitrate在NT vs. WL时p.adj<0.05,差异具有统计学意义。

扫描上方二维码，关注EcoEvoPhylo数据分析小店，购买数据分析书籍和服务。

后台发送“One-Way MANOVA”，获取数据与代码文件。

参考文献

马氏距离及其几何解释

换个姿势看马氏距离和主成分分析

https://www.datanovia.com/en/lessons/one-way-manova-in-r/

推荐阅读

R中进行单因素方差分析并绘图

R统计-多变量单因素参数、非参数检验及多重比较

R绘图-相关性分析及绘图

R绘图-相关性系数图

R统计绘图-环境因子相关性热图

R统计绘图-corrplot绘制热图及颜色、字体等细节修改

R统计绘图-corrplot热图绘制细节调整2(更改变量可视化顺序、非相关性热图绘制、添加矩形框等)

R绘图-RDA排序分析

R统计-VPA分析(RDA/CCA)

R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程

R数据可视化之美-节点链接图

R统计绘图-rgbif包下载GBIF数据及绘制分布图

R统计-正态性分布检验[Translation]

R统计-数据正态分布转换[Translation]

R统计-方差齐性检验[Translation]

R统计-Mauchly球形检验[Translation]

R统计绘图-单、双、三因素重复测量方差分析[Translation]

R统计绘图-混合方差分析[Translation]

R统计绘图-协方差分析[Translation]

开放转载

公众号原创文章开放转载，在文末留言告知即可

EcoEvoPhylo ：主要分享微生物生态和phylogenomics的数据分析教程。扫描上方二维码，即可关注EcoEvoPhylo。让我们大家一起学习，互相交流，共同进步。

R统计绘图-One-Way MANOVA相关推荐

R统计绘图 | 物种组成冲积图(绝对/相对丰度，ggalluvial)
一.数据准备数据使用的不同处理土壤样品的微生物组成数据,包含物种丰度,分类单元和样本分组数据.此数据为虚构,可用于练习,请不要作他用. # 1.1 设置工作路径 #knitr::opts_knit$ ...
R统计绘图-PCA详解1(princomp/principal/prcomp/rda等)
此文为<精通机器学习:基于R>的学习笔记,书中第九章详细介绍了无监督学习-主成分分析(PCA)的分析过程和结果解读. PCA可以对相关变量进行归类,从而降低数据维度,提高对数据的理解.分析 ...
R统计绘图-VPA(变差分解分析)
变差分解分析(Variance Partitioning Analysis)可用于确定指定环境因子对微生物(原生生物/植物/动物等等)群落结构变化的解释比例.要计算指定环境因子与群落结构的相关性,就需 ...
R统计绘图 | 物种组成堆叠柱形图(绝对/相对丰度)
一.数据准备数据使用的不同处理土壤样品的微生物组成数据,包含物种丰度,分类单元和样本分组数据.此数据为虚构,可用于练习,请不要作他用. # 1.1 设置工作路径 #knitr::opts_knit$ ...
R统计绘图 | 物种组成堆叠面积图(绝对/相对丰度，ggalluvial)
一.数据准备数据使用的不同处理土壤样品的微生物组成数据,包含物种丰度,分类单元和样本分组数据.此数据为虚构,可用于练习,请不要作他用. # 1.1 设置工作路径 #knitr::opts_knit$ ...
R统计绘图-多元线性回归(最优子集法特征筛选及模型构建，leaps)
此文为<精通机器学习:基于R>的学习笔记,书中第二章详细介绍了线性回归分析过程和结果解读. 回归分析的一般步骤: 1. 确定回归方程中的自变量与因变量. 2. 确定回归模型,建立回归方程. ...
R统计绘图-多元线性回归(平均加权模型/最优子集筛选，MuMIn)
此文介绍如何使用MuMIn包使用最优子集法进行多重线性回归的模型筛选以及模型平均.多重线性回归需要进行的数据检验过程都写在R统计绘图-多重线性回归(最优子集法特征筛选,leaps)中了.大家可以自行查 ...
R统计绘图-PCA分析及绘制双坐标轴双序图
zhe 点击名片关注我们有师妹来咨询,怎样画类似于上图的双坐标轴PCA双序图.正好之前虽然PCA和RDA分析及绘图都写过教程,但是变量分析结果没有在图中显示,所以使用R统计绘图-环境因子相关性 ...
R统计绘图-corrplot热图绘制细节调整2(更改变量可视化顺序、非相关性热图绘制、添加矩形框等)
上一篇文章推送的是怎样调整corrplot热图的可视化参数,以修改字符和图例位置,数据可视化形式和字符小大和颜色等这篇是一个补充部分,记录怎样修改参数以变量排序方式和突出部分数据.本流程还是使用R统计 ...

R统计绘图-One-Way MANOVA

二、数据准备

R统计绘图-One-Way MANOVA相关推荐

最新文章

热门文章