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。让我们大家一起学习,互相交流,共同进步。

  • 2
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

EcoEvoPhylo

值得点赞吗?

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值