查看: 3535|回复: 29

[R语言] R统计绘图-One-Way MANOVA

  [复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
29
注册时间
2017.4.3
在线时间
107 小时

发表于 2021.9.9 09:55:36 | 显示全部楼层 |阅读模式
存在两个及以上连续结果(或响应)变量的ANOVA被称为多变量方差分析(Multivariate Analysis of Variance,MANOVA)。例如,将小鼠分为处理A和处理B两组后,测量小鼠的长度和高度。此时小鼠的长度和高度就是结果变量,假设长度和高度都因为处理不同产生差异。此时可以使用MANOVA检测上述假设。
MANOVA分析过程总结如下:
1)将所有结果变量经过线性组合形成一个新的复合变量;2)比较新变量在分类变量中的均值。
本文接下来描述在R中如何进行MANOVA。


图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]。
  1. # 3.1.1 检查样本量是否足够
  2. ENV %>%
  3.   group_by(tillage) %>%
  4.   summarise(N = n())
  5. ## 上面统计结果显示tillage每个组中只有9个样本,而环境变量有11个,不满足假设。
  6. ## 因此对环境变量进行筛选,只保留氮相关的环境因子
  7. env = ENV[c(2,4,7:9)]
  8. set.seed(123)
  9. env %>% sample_n_by(tillage,size = 2) # 基于tillage分组显示数据,每组随机显示两个样本。


  10. # 3.1.2 检测离群值
  11. ## 检测单变量TN的各组离群值
  12. env %>% group_by(tillage) %>%
  13.   identify_outliers(TN) # 没有


  14. ## 检测多变量离群值:多变量离群值指在结果变量上具有异常结合值得数据点。
  15. ## Mahalanobis distance经常用于在MANOVA过程中检测多变量离群值。这个检验会检测观察值距离均值有多远,还会考虑到协方差。
  16. ## rstatix包的mahalanobis_distance()用于计算Mahalanobis distance并标注多变量离群值。
  17. env %>% group_by(tillage) %>%
  18.   mahalanobis_distance(names(env[2:5])) %>%
  19.   filter(is.outlier == TRUE) %>%
  20.   as.data.frame() # 根据Mahalanobis距离(p>0.001),数据中没有多变量异常值。


  21. # 3.1.3 正态分布假设检验
  22. ## 单变量正态分布假设检验
  23. env %>%
  24.   group_by(tillage) %>%
  25.   shapiro_test(TN,Ammonia,Nitrate,AHN) %>%
  26.   arrange(variable) # Shapiro检验,p<0.05,该分组数据不满足正态分布假设。
  27. ?ggqqplot
  28. ggqqplot(env,names(env[2:5]),
  29.          facet.by = "tillage",
  30.          ggtheme = theme_bw()) # 四个变量会产生四幅组合QQ图。样本量大于50时推荐使用QQ图检测正态分布假设。
  31. ## 多变量正态分布假设
  32. env %>%
  33.   select(names(env[2:5])) %>%
  34.   mshapiro_test() # p<0.05,多变量正态分布假设不满足。可以选择其他非参数检验方法。


  35. # 3.1.4 检测多重共线性
  36. ## 结果变量间的相关性应该是居中的,如果相关性大于0.9表明有多重共线性,不适合使用MANOVA,可以考虑删除一个高度相关的结果变量。如果相关性值太低,则需要考虑对每个结果变量进行单因素ANOVA分析。
  37. #env %>% cor_test(names(env[2:5]),method = "pearson") # 适合只有两个结果变量时使用
  38. env %>% cor_mat(names(env[2:5]),method = "pearson") # 没有两个结果变量间相关性系数大于0.9。


  39. # 3.1.5 检测线性假设,每组的结果变量应该是线性关系。
  40. ?ggpairs
  41. results <- env %>% select(names(env[c(1:5)])) %>%
  42.   group_by(tillage) %>%
  43.   doo(~ggpairs(.) + theme_bw(),result = "plots")
  44. results$plots[1] # 各结果变量在各分组中没有检测到线性关系,可以转换或删除有问题的结果变量。


  45. # 3.1.6 检测协方差同质性假设
  46. box_m(env[2:5],env$tillage) #p<0.001,数据违反了方差-协方差同质性假设。
  47. ## 如果是均衡实验设计,每个分组具有相同的样本大小,违背方差-协方差同质性假设,仍然可以继续进行分析。
  48. ## 如果实验设计不均衡,则可以1)转换变量数据以满足此假设或2)使用Pillai多变量统计值代替Wilks统计值。


  49. # 3.1.7 检测方差同质性假设,检测每个变量在各分组中是否具有方差齐性。
  50. env %>% gather(key = "variable",value = "value",names(env[2:5])) %>%
  51.   group_by(variable) %>%
  52.   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统计。
  1. model <- lm(cbind(TN,Ammonia,Nitrate,AHN) ~ tillage,env)
  2. model
  3. 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 单变量单因素方差分析
  1. # 将数据由wide format转换为long format
  2. data <- env %>%
  3.   gather(key = "variables",value = "value",names(env[2:5])) %>%
  4.   group_by(variables) # 此处变量名一定要命名为variables,否则后续添加显著性标签位置时,会报错。
  5. data


  6. # anova_test
  7. data %>% anova_test(value ~ tillage) %>%
  8.   adjust_pvalue(method = "bonferroni")


  9. # welch anova test
  10. data %>% welch_anova_test(value ~ tillage) %>%
  11.   adjust_pvalue(method = "bonferroni")


  12. # Kruskal test
  13. data %>% kruskal_test(value ~ tillage) %>%
  14.   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)也可用于对不满足方差齐性假设的数据进行两两比较。
  1. res <- data %>%
  2.   games_howell_test(value ~ tillage) %>%
  3.   select(-estimate,-conf.low,-conf.high) #移除一些细节
  4. res
复制代码


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

[size=0.83em]微信店铺.jpg (87.74 KB, 下载次数: 0)
下载附件
[color=rgb(153, 153, 153) !important]4 天前 上传



扫描上方二维码,关注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。让我们大家一起学习,互相交流,共同进步。              

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
努力奋斗
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
79 小时

发表于 2021.9.9 11:01:14 | 显示全部楼层
当工人也要加油
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2021.1.28
在线时间
10 小时

发表于 2021.9.9 14:22:46 | 显示全部楼层
顶!!!!!!!!!!!!!!!!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.9 16:22:22 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 17:48:57 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 18:23:17 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 18:24:33 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 18:25:13 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 18:25:51 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
137 小时

灌水之王


发表于 2021.9.9 18:26:57 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.10 08:41:10 | 显示全部楼层
加油,加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.2.25
在线时间
2 小时

发表于 2021.9.15 22:30:20 | 显示全部楼层
加油加油
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2020.11.23
在线时间
13 小时

发表于 2021.9.20 21:56:58 | 显示全部楼层
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.23 08:13:59 | 显示全部楼层
加油,加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2021.1.28
在线时间
10 小时

发表于 2021.9.26 10:46:49 | 显示全部楼层
顶!!!!!!!!!!!!!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.27 11:22:09 | 显示全部楼层
加油,加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.28 08:12:19 | 显示全部楼层
加油,加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2021.1.28
在线时间
10 小时

发表于 2021.9.29 18:23:43 | 显示全部楼层
顶!!!!!!!!!!!!!!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.9.30 08:32:25 | 显示全部楼层
加油,加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.10.4 21:35:03 | 显示全部楼层
加油,加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表