lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

单因素、双因素、协方差、多元方差分析的非参数检验替代方法在R中实现

已有 6074 次阅读 2019-7-24 20:38 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _差异分析, _非参数检验

单因素、双因素、协方差、多元方差分析的非参数检验替代方法在R中实现

 

对于多组数据间比较差异,我们常会想到使用方差分析(ANOVA来实现。不过由于ANOVA的前提假设条件比较严格,要求数据必须满足正态性、方差齐性等,而很多情况下我们的数据并不符合方差分析的条件。通常情况下,我们可以考虑转换数据,例如,使用log转换等或许可以使非正态分布的原始数据转变为正态分布类型(当然,我们需要确保转换后的数据能够被合理解释,否则将无意义);另一种合适的选择是使用非参数的检验方法,替代ANOVA

本文就几种常见的ANOVA为例,介绍与其相关的几种非参数检验替代方法在R中的实现过程。包括:

单因素方差分析(单因素ANOVA)的非参数检验替代方法,如Kruskal-Wallis检验(kruskal.test())、Friedman检验(friedman.test());

单因素协方差分析(ANCOVA)的非参数检验替代方法,smsm.ancova()

双因素方差分析(双因素ANOVA)的非参数检验替代方法,Scheirer-Ray-Hare检验,rcompanionscheirerRayHare()

多元方差分析(MANOVA)的非参数检验替代方法,置换多元方差分析(PERMANOVA),veganadonis()

备注:当数据满足ANOVA的条件时,尽可能使用ANOVA分析,相较于非参数方法,能够有效地鉴别出非参数检验鉴别不到的差异。非参数的检验方法普遍保守。为了更好地和ANOVA有所比较,本文使用的测试数据和前文方差分析中的测试数据保持一致(因此测试数据即可使用ANOVA也可使用非参数方法)。可以比较查看ANOVA和非参数检验的结果,例如在下文Scheirer-Ray-Hare检验中,treattimes的交互作用不显著;而在前文双因素ANOVA中,treattimes的交互则是显著的。当然,若我们无法使用ANOVA时,就另当别论了。

 

测试数据及整理好的R代码可在以下百度盘链接中获取(提取码fjuq):

https://pan.baidu.com/s/1I2dX58-q5XJgexxfYTvpmg

 

我们首先将示例数据读到R中。对于数据中每列变量的含义,详见下文。

#读入文件
soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group  <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')

 


Kruskal-Wallis检验(单因素ANOVA的非参数替代方法)


单因素方差分析要求数据满足正态性和方差齐性,若二者之一不符合,如果各组独立,则Kruskal-Wallis检验将会是一种实用的方法。

 

示例数据说明


从总示例数据中挑选部分数据作为测试(测试数据同前文中的单因素ANOVA,以方便二者比较)。

#以 chao1 指数为例,同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'site', 'chao1')]
chao1$site <- factor(chao1$site)

1.png

2.png

假设存在这么一个研究:

我们在3个地域(ABC)分别采集了土壤样本,即获得了3种类型的土壤,并通过16S测序,获得了每种类型土壤中细菌群落的Alpha多样性指数。我们想要得知,三种土壤环境下的细菌群落的Alpha多样性指数是否存在显著不同。

对应于上述挑选出的测试数据“chao1”:sample,采集的土壤样本名称;site,土壤样本来源的环境(ABC),这列作为分组列,需要转换为因子变量类型;chao1Alpha多样性指数中的Chao1指数,数值变量。

在前文方差分析中,我们已经使用了单因素方差分析(单因素ANOVA)对3种土壤环境下的细菌群落的Chao1指数进行比较。在这里,我们更换为Kruskal-Wallis非参数检验比较差异。

 

Kruskal-Wallis检验


Kruskal-Wallis检验的调用格式为:kruskal.test(y~A, data)。其中y是数值型结果变量(因变量),A是因子类型的分组变量(因子变量)。

我们调用kruskal.test()完成检验,对于该函数详情可实用?kruskal.test查看帮助。

#当数据不满足单因素方差分析的条件(正态性、方差齐性)时,尝试非参数的方法,例如
#Kruskal-Wallis Test
kruskal.test(chao1~site, data = chao1)

Kruskal-Wallis检验结果表明,3种土壤环境下的细菌群落的Chao1指数具有显著差异,p值远低于0.05水平。

3.png

 

非正态性的数据,不太适合以均值及标准差的形式作为展示,可以使用中位数等代替。

#若想查看各组中位数、极大(小)值、(四)分位数,可使用 aggregate()
aggregate(chao1$chao1, by = list(chao1$site), FUN = max)
aggregate(chao1$chao1, by = list(chao1$site), FUN = function(x) quantile(x, 0.75))
aggregate(chao1$chao1, by = list(chao1$site), FUN = median)
aggregate(chao1$chao1, by = list(chao1$site), FUN = function(x) quantile(x, 0.25))
aggregate(chao1$chao1, by = list(chao1$site), FUN = min)

4.png

 

事后两两比较


上述Kruskal-Wallis检验告诉我们3种土壤环境下的细菌群落的Chao1指数具有显著差异,这种差异是在整体水平而言的,并没有告诉我们究竟谁和谁存在差异。在ANOVA中,我们常使用Tukey HSD检验基于ANOVA的结果继续执行组间两两比较;而在非参数检验中,则可以考虑使用Wilcoxon秩和检验。

以下使用Wilcoxon秩和检验继续探究两两差异。当分组水平较多时,多次的Wilcoxon秩和检验容易提升I类错误的风险,因此推荐引入p值校正过程。

##此时若想继续探寻两两分组间的差异,可使用 Wilcoxon 秩和检验,如果分组较多,可以使用循环来完成
#继续在上述 Kruskal-Wallis 检验的基础上,探究两两分组间差异,一个示例如下
group <- levels(chao1$site)
group1 <- NULL
group2 <- NULL
median1 <- NULL
median2 <- NULL
p <- NULL
 
for (i in 1:(length(group) - 1)) {
        for (j in (i + 1):length(group)) {
                group1 <- c(group1, group[i])
                group2 <- c(group2, group[j])
                group_ij <- subset(chao1, site %in% c(group[i], group[j]))
                group_ij$site <- factor(group_ij$site, levels = c(group[i], group[j]))
                
                wilcox_test <- wilcox.test(chao1~site, data = group_ij, alternative = 'two.sided', conf.level = 0.95)
                p <- c(p, wilcox_test$p.value)
                median1 <- c(median1, median(subset(group_ij, site == group[i])$chao1))
                median2 <- c(median2, median(subset(group_ij, site == group[j])$chao1))
        }
}
 
result <- data.frame(group1, group2, median1, median2, p)
result$padj <- p.adjust(result$p, method = 'holm') #推荐加上 p 值校正,这里使用 Holm 方法校正 p 值
result
#write.table(result, 'Wilcoxon.txt', sep = '\t', row.names = FALSE, quote = FALSE)

示例结果如下,展示两组比较的组名、中位数、Wilcoxon检验p值以及Holm校正后的p值。

5.png

 

ggpubr包箱线图可视化示例


ggpubr包是ggplot2的拓展,作图风格和ggplot2一致,并且内置了多种统计检验方法。我们可以很方便地借助该包,一次完成统计检验及可视化。

#使用 ggpubr(ggplot2 的扩展包),作图展示 wilcox 比较的结果
library(ggpubr)
 
ggboxplot(data = chao1, x = 'site', y = 'chao1', color = 'site') +
stat_compare_means(method = 'wilcox.test', comparisons = list(c('A', 'B'), c('A', 'C'), c('B', 'C')))

例如这里我们绘制了箱线图展示数据分布,并同时执行了Wilcoxon秩和检验,将显著性p值标注在图中。Wilcoxon秩和检验p值和上文一致。

备注:stat_compare_means()支持T检验(t.test)、ANOVAanova)、Wilcoxon秩和检验(wilcox.test)以及Kruskal-Wallis检验(kruskal.test)。若使用T检验或ANOVA,一定要在数据满足前提假设时使用。

6.png

 


Friedman检验


Kruskal-Wallis检验相比,如果各组不独立(如重复测量设计或随机区组设计),那么Friedman检验可能会更合适。

 

Friedman检验的调用格式为:friedman.test(y~A|B, data)。其中y是数值型结果变量(因变量),A是因子类型的分组变量(因子变量),而B是一个用以认定匹配观测的区组变量(blocking variable)。

我们直接以R自带的数据集做演示吧。使用?friedman.test查看帮助文档,最下方提供了示例。

#如果各组不独立(如重复测量设计或随机区组设计),那么 Friedman Test 可能会更合适
#这里以 ?friedman.test 中的示例为例
wb <- aggregate(warpbreaks$breaks, by = list(w = warpbreaks$wool, t = warpbreaks$tension), FUN = mean)
friedman.test(x ~ w | t, data = wb)

这个示例的p值不显著。

7.png

 


单因素协方差分析(ANCOVA)的非参数替代方法


ANCOVA除了要求数据服从正态分布,以及各组方差相等外,还假定回归斜率相同(回归斜率的同质性)。如不满足,可采用非参数的方法。

 

示例数据说明


从总示例数据中挑选部分数据作为测试(测试数据同前文中的ANCOVA,以方便二者比较)。

#以 shannon 指数为例,同时将分组列转换为因子变量
shannon <- soil[ ,c('sample', 'treat', 'shannon', 'days')]
shannon$treat <- factor(shannon$treat)

8.png9.png


假设存在这么一个研究:

我们使用来源于同一环境中的土壤进行盆栽试验(土壤类型一致),并种植了同种植株(植物类型一致)。我们将盆栽(1植株/1盆栽)分为了3组,分别在土中添加了三种化学物质(abc);然后等待植物到达花期后,收集每个植株的根际土,并通过16S测序,获得了植物根际细菌群落的Alpha多样性指数,意在探究不同类型的化学物质是如何影响植物根际细菌群落的。但由于各盆栽中植株开花期时间并非一致,不同个体间在开花时间上可能相隔数天(这个我们无法控制),而该植物开花时间又很短,可能早开花的植株还未等晚开花的植株开花就已经凋谢,因此无法保证所有植株均在同一天取样观察,于是我们就先开花的先取样后开花的后取样。尽管期间并未相隔很长时间(天),理论上单因素ANOVA就可以满足需求,但我们仍然想要将植物生长时间这个因素考虑在内(如果觉得植物生长时间相隔几天的差异可能会导致其根际菌群产生较大的变异时),即将它作为协变量处理(协变量是独立变量,实验者不能操纵,但仍影响实验结果)。

对应于上述挑选出的测试数据“shannon”:sample,试验样本名称,每一个样本即对应一个盆栽,各盆栽中土壤类型、植物类型完全相同,而且均是1植株/1盆栽;treat,在土壤中添加的三种化学物质(abc),这列作为分组列,需要转换为因子变量类型;shannonAlpha多样性指数中的Shannon指数,数值变量;days,各盆栽中植株的开花时间(即生长天数),这里为数值变量。

在前文方差分析中,我们以植物根际菌群的Shannon指数为因变量,植物生长天数(days)为协变量,使用了单因素协方差分析(ANCOVA),探究不同类型的化学物质处理下的植物根际细菌群落的Shannon指数是否存在显著不同。在这里,我们更换为非参数的方法。

 

非参数检验(不好意思,我实在没查到下面用的这个方法叫什么名字……


smsm.ancova()函数提供了单因素协方差分析的非参数替代方法,详情使用?sm.ancova查看详情。

接下来我们使用这个函数完成上述统计。

#当数据不满足单因素协方差分析的条件(正态性、方差齐性、回归斜率的同质性)时,尝试非参数的方法,例如
#sm 包 sm.ancova()
library(sm)
sm.ancova(x = shannon$days, y = shannon$shannon, group = shannon$treat, model = 'equal')

式中,x为协变量,y为因变量,group为分组因子,method定义参考模型(详见帮助)。

10.png

11.png

 


Scheirer-Ray-Hare检验(双因素ANOVA的非参数替代方法)


接下来介绍双因素。同样地,双因素方差分析要求数据满足正态性和方差齐性,若二者之一不符合,Scheirer-Ray-Hare检验将会是一种实用的方法。

写到这儿,先感谢一下大神师妹(不对,严格来说现在得喊人师姐了……)。回想当初想找个合适的方法替代双因素ANOVA(没办法,数据没法用ANOVA分析),翻书没见到有介绍相关方法的,网上也没搜到,就用Wilcoxon秩和检验逐一比较,然后把所有可能的组合都试了……后来还是师妹告知的这个检验,问师妹怎么找到的,师妹说很简单就找见了,嗯嗯,一定是我找的方式不对……

 

示例数据说明


从总示例数据中挑选部分数据作为测试(测试数据同前文中的双因素ANOVA,以方便二者比较)。

#同样以 chao1 指数为例,同时将分组列转换为因子变量
chao1 <- soil[ ,c('sample', 'treat', 'times', 'chao1')]
chao1$treat <- factor(chao1$treat)
chao1$times <- factor(chao1$times)

12.png

13.png

假设存在这么一个研究:

我们采集了来源于同一环境中的土壤(土壤类型一致),分为了3组,分别添加了三种类型的化学物质(abc),并将土壤孵育了不同的时间(10152025 天)。在不同时间段收集样本后,通过16S测序,获得了土壤细菌群落的Alpha多样性指数,期望关注化学物质类型以及处理时间对土壤细菌群落的影响(关注交互作用)。

对应于上述挑选出的测试数据“chao1”:sample,试验样本名称,土壤类型完全相同;treat,在土壤中添加的三种化学物质(abc),这列作为分组列,需要转换为因子变量类型;times,土壤孵育时间(天数),这列作为分组列,需要转换为因子变量类型(尽管它本来是数值类型的,但必需转化为因子类型后,才可作为分组变量用于方差分析);chao1Alpha多样性指数中的Chao1指数,数值变量。

此处存在“化学物质类型”以及“土壤孵育时间”两组分组变量,对应于双因素,同时还需要关注二者的交互作用。在前文方差分析中,我们已经使用双因素方差分析(双因素ANOVA)探究了化学物质类型以及处理时间是否对土壤细菌群落产生了显著影响。在这里,我们更换为Scheirer-Ray-Hare检验来实现。

 

Scheirer-Ray-Hare检验


Scheirer-Ray-Hare检验可使用rcompanionscheirerRayHare()函数来实现,详情使用?scheirerRayHare查看帮助说明。

接下来我们使用scheirerRayHare()完成上述统计。

#当数据不满足双因素方差分析的条件(正态性、方差齐性)时,尝试非参数的方法,例如
#Scheirer-Ray-Hare Test,rcompanion 包 scheirerRayHare()
library(rcompanion)
scheirerRayHare(chao1~treat*times, data = chao1)

统计结果屏幕输出如下,包含两组因素及其交互作用的检验结果。

14.png

 


置换多元方差分析(PERMANOVA)(多元方差分析(MANOVA)的非参数替代方法)


当因变量不止一个时,即一个或多个因子变量对应了多个因变量时,可使用多元方差分析(MANOVA)。但是MANOVA的条件非常苛刻,以最常见的单因素MANOVA为例,相较于上述提到的常规单因素、双因素ANOVA,要求数据满足多元正态性、方差-协方差矩阵同质性,大部分案例中都是直接拒绝的。

替代方法可以使用稳健MANOVA,这个在前文方差分析中已经提到。对于非参数的方法,常用置换多元方差分析(PERMANOVA)。

 

示例数据说明


从总示例数据中挑选部分数据作为测试(测试数据同前文中的MANOVA)。

#选择数据,并将分组列转换为因子变量
muti <- soil[ ,c('sample', 'treat', 'chao1', 'pH', 'NR')]
muti$treat <- factor(muti$treat)

15.png

16.png

假设存在这么一个研究:

我们采集了来源于同一环境中的土壤(土壤类型一致),分为了3组,分别添加了三种类型的化学物质(abc),并将土壤孵育了一定的时间(时间相同)。最后取样后,通过16S测序,获得了土壤细菌群落的Alpha多样性指数;通过土壤理化测定,获得了土壤中多种理化指标;通过土壤酶活性测定,获得了主要的几种土壤酶活性数据。通过这些数据,我们想要得知土壤细菌群落的Alpha指数、土壤理化、以及土壤酶活性是否因所添加化学物质类型不同而显著改变。

对应于上述挑选出的测试数据“muti”:sample,试验样本名称,土壤类型完全相同,土壤孵育时间完全一致;treat,在土壤中添加的三种化学物质(abc),这列作为分组列,需要转换为因子变量类型;chao1Alpha多样性指数中的Chao1指数,数值变量;pH,土壤理化数据中的pH值数据,数值变量;NR,土壤硝酸还原酶活性,数值变量。

这里存在3组因变量:“土壤菌群Chao1指数”、“土壤pH值”、“土壤硝酸还原酶(NR)活性”,对应于1组因子变量“化学物质类型”,因此我们期望使用单因素多元方差分析(单因素MANOVA),探究土壤细菌群落Chao1指数、土壤pH、土壤硝酸还原酶(NR)活性是否因所添加化学物质类型不同而发生显著改变。但是在评估检验的假设条件时,发现方差-协方差矩阵同质性并不满足,于是我们考虑更换为非参数的方法作为替代。

 

置换多元方差分析(PERMANOVA


置换多元方差分析(PERMANOVA)可以通过veganadonis()函数来实现,详细使用?adonis查看帮助。(我也曾经在这篇博文中简介了adonis()的使用,可点击查看)

以下为本文的示例参考。

#当数据不满足单因素 MANOVA 的条件(多元正态性、方差-协方差矩阵同质性)时,尝试非参数的方法,例如
#置换多元方差分析(PERMANOVA),vegan 包 adonis()
library(vegan)
 
rownames(muti) <- muti$sample
muti_dat <- muti[ ,c('chao1', 'pH', 'NR')]
#muti_dat <- decostand(muti_dat, 'standardize')   #可选 z-score 标准化,消除量纲差异
dis <- vegdist(muti_dat, method = 'euclidean')      #获取欧式距离矩阵(这里以欧式距离为例,实际分析中视情况选择合适的距离测度)
adonis(dis~treat, muti, permutations = 999)  #999 次置换获得 p 值(结合了置换检验)

结果是显著的。

17.png

 

以上示例中只有一组分组因子。当存在多组分组时(多因素时),adonis()仍然适用。

#多因素下,在 ~ 后面添加分组因子即可,例如存在 A、B 双因素时
#adonis(y~A+B, data, permutations = 999)

 



http://blog.sciencenet.cn/blog-3406804-1190972.html

上一篇:多组间差异分析之方差分析(ANOVA)在R中实现
下一篇:线粒体在线注释网站MITOS使用简介

0

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 17:47

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部