R统计绘图-One-Way MANOVA

最新推荐文章于 2023-06-01 08:30:00 发布

EcoEvoPhylo

最新推荐文章于 2023-06-01 08:30:00 发布

阅读量2.6k

点赞数 2

分类专栏： R统计绘图方差分析箱线图文章标签： r语言数据分析

本文链接：https://blog.csdn.net/qq_39859424/article/details/120195526

版权

R统计绘图同时被 3 个专栏收录

17 篇文章 24 订阅

订阅专栏

方差分析

4 篇文章 1 订阅

订阅专栏

箱线图

2 篇文章 1 订阅

订阅专栏

存在两个及以上连续结果(或响应)变量的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。让我们大家一起学习，互相交流，共同进步。

EcoEvoPhylo

关注

2
点赞
踩
16

收藏

觉得还不错? 一键收藏
打赏
1
评论
R统计绘图-One-Way MANOVA

存在两个及以上连续结果(或响应)变量的ANOVA被称为多变量方差分析(Multivariate Analysis of Variance,MANOVA)。例如，将小鼠分为处理A和处理B两组后，测量小鼠的长度和高度。此时小鼠的长度和高度就是结果变量，假设长度和高度都因为处理不同产生差异。此时可以使用MANOVA检测上述假设。MANOVA分析过程总结如下：1)将所有结果变量经过线性组合形成一个新的复合变量；2)比较新变量在分类变量中的均值。本文接下来描述在R中如何进行MANOVA。一、
复制链接

扫一扫