浅析R语言非参数检验的多组比较及分面与分组的图形艺术

浅析R语言多组定量资料非参数检验的多组比较及簇状柱形图显著性字母标记之分面与分组的图形艺术

R语言多组定量资料非参数检验的多组比较

非参数检验的应用

本流程是在刘永鑫老师提供的代码资料指导下完成

先简单说明一下非参数检验

参数检验是在已知总体分布的条件下(一般要求总体服从正态分布)对一些主要的参数(如均值、百分数、方差、相关系数等)进行的检验，有时还要求某些总体参数满足一定的条件。如独立样本的t检验和方差分析不仅要求总体符合正态分布，还要求各总体方差齐性。

而当我们的数据不满足参数检验时我们可以选择非参数检验-秩和检验(Wilcoxon/Kruskal-Wallis H)

非参数检验方法简便，不依赖于总体分布的具体形式，因而适用性强，但灵敏度和精确度不如参数检验。一般而言，非参数检验适用于以下3种情况

(1)顺序类型的数据资料，这类数据的分布形态一般是未知的。

(2)虽然是连续数据，但总体分布形态未知或者非正态，这和卡方检验一样，称为自由分布检验。

(3)总体分布虽然正态，数据也是连续类型，但样本容量极小，如10以下(虽然T检验被
称为小样本统计方法，但样本容量太小时，代表性毕竟很差，最好不要用要求较严格的参数检验法)。

非参数检验的缺点与不足

因为有这些特点，加上非参数检验法一般的原理和计算比较简单，因此常用于一些为正式研究进行探路的预备性研究的数据统计中。当然，由于非参数检验许多牵涉不到参数计算，对数据中的信息利用不够，因而其统计检验力相对参数检验则差得多。

关于使用非参数检验的纠结

在很多人看来，参数检验似乎是统计方法的“主流”，而非参数检验往往被当成“非主流”。大家似乎更喜欢用t检验、方差分析这样的参数检验。即使在数据不满足正态分布的条件下，仍然有使用参数检验的执念；总觉得非参数检验不够正式，或者就是不喜欢用。

而现实中的很多数据都不服从正态分布，这时候非参数检验的统计学效能要高于参数检验。并且，非参数检验更加稳健。参数统计方法是建立在严格假设条件基础上，一旦假设条件不符合，其推断的正确性就不存在了。非参数检验带有最弱的假设，对模型的限制很少，因而天然的具有稳健性。

多组定量资料非参数检验及多组比较的实现

kruskal.test()函数可以确定总体差异是否有统计学意义，但是不知道具体那些组之间存在差异。可以使用wilcox.test()函数每次进行两组数据比较。一种更为有效的方法是在控制犯第一类错误的概率前提下,执行可以同步进行的多组比较,这样可以直接完成所有组之间的成对比较。nparcomp包提供了所需要的非参数多组比较程序。此包中的nparcomp()函数接受的输入为一个两列的数据框,其中一列为因变量,另一列为分组变量。==本推文的第三种方法采用此方法==。

本文分别采用三个统计R包的三种不同函数方法进行非参数多组比较：

  • “pgirmess”统计包的kruskalmc()函数

  • “PMCMRplus”统计包的posthoc.kruskal.nemenyi.test()函数

  • “nparcomp”统计包的nparcomp()函数

自己挑选喜欢的用

library(tidyverse)#数据整理与数据转换用了一些更好用更易懂的函数
library(ggprism)
dat <- read.delim('./vegan.txt')#读入α多样性数据
head(dat, n = 3)
design <- read.delim('./metadata.txt')#读入试验设计文件
head(design, n = 3)
dat <- design %>%
dplyr::select(SampleID,Group) %>%
inner_join(dat,by='SampleID')#按照SampleID内连接文件
head(dat, n = 3)
dat <- gather(dat,alpha,v,-(SampleID:Group))#宽数据变长数据
head(dat, n = 3)

函数和参数简介

函数参数设置：

  • data就是上面整好的数据,

  • group是你的分组信息列，比如α多样性的种类(或不同的基因),

  • compare是每个α多样性要比较的不同处理(或每个gene要比较的不同处理),

  • value 值就是要比较的α多样性/gene拷贝数的数值。

整体思想如下(例如本数据)：
首先给输入数据dat，根据alpha列分成不同的小子集，每个小子集比较不同Group下v值的差异情况，最后汇总输出。

“pgirmess”统计包的kruskalmc()函数

# 1 -----------------------------------------------------------------------
kruskalmc_compare1 <- function(data,group,compare,value){
library(multcompView)
library(pgirmess)##多组两两比较函数用到的R包
library(multcomp)
a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
g1=compare
sub_dat <- data[data[,group]==i,]
names(sub_dat)[names(sub_dat)==compare] <- 'g1'
names(sub_dat)[names(sub_dat)==value] <- 'value'
sub_dat$g1 <- factor(sub_dat$g1)
options(warn = -1)

k <- kruskalmc(value ~ g1, data=sub_dat, probs=0.05)
dif <- k$dif.com[['difference']]
names(dif) <- rownames(k$dif.com)
difL <- multcompLetters(dif)
label <- data.frame(difL['Letters'], stringsAsFactors = FALSE)
label$compare = rownames(label)
label$type <- i

mean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),
aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1"
)
names(mean_sd) <- c('compare','std','mean')
a <- rbind(a,merge(mean_sd,label,by='compare'))
}
names(a) <- c(compare,'std','mean','Letters',group)
return(a)
}

“PMCMRplus”统计包的posthoc.kruskal.nemenyi.test()函数

# 2 -----------------------------------------------------------------------
PMCMR_compare1 <- function(data,group,compare,value){
library(multcompView)
library(multcomp)
library(PMCMRplus)
library(PMCMR)##多组两两比较函数用到的R包
a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
g1=compare
sub_dat <- data[data[,group]==i,]
names(sub_dat)[names(sub_dat)==compare] <- 'g1'
names(sub_dat)[names(sub_dat)==value] <- 'value'
sub_dat$g1 <- factor(sub_dat$g1)
options(warn = -1)

k <- posthoc.kruskal.nemenyi.test(value ~ g1,data=sub_dat)
n <- as.data.frame(k$p.value)
h <- n %>%
mutate(compare=rownames(n)) %>%
gather(group,p,-compare,na.rm = TRUE) %>%
unite(compare,group,col="G",sep="-")
dif <- h$p
names(dif) <- h$G
difL <- multcompLetters(dif)
K.labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)
K.labels$compare = rownames(K.labels)
K.labels$type <- i

mean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),
aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1"
)
names(mean_sd) <- c('compare','std','mean')
a <- rbind(a,merge(mean_sd,K.labels,by='compare'))
}
names(a) <- c(compare,'std','mean','Letters',group)
return(a)
}

“nparcomp”统计包的nparcomp()函数

# 3 -----------------------------------------------------------------------
nparcomp_compare1 <- function(data,group,compare,value){
library(nparcomp)##多组两两比较函数用到的R包
library(multcompView)
library(do)
a <- data.frame(stringsAsFactors = F)
type <- unique(data[,group])
for (i in type)
{
g1=compare
sub_dat <- data[data[,group]==i,]
names(sub_dat)[names(sub_dat)==compare] <- 'g1'
names(sub_dat)[names(sub_dat)==value] <- 'value'
sub_dat$g1 <- factor(sub_dat$g1)

k <- nparcomp::nparcomp(value ~ g1,data=sub_dat, alternative = "two.sided")
k$Analysis
b <- k$Analysis
dif <- b$p.Value
names(dif) <- b$Comparison
difname <- names(dif) %>%
Replace(from = ' , ',to='-') %>%
Replace(from=c('p\\(','\\)'),to='')#正则表达式替换
temp_name <- data.frame(tn=difname) %>%
separate(col = 'tn',sep = '-',into = c('n1','n2'))
difname <- paste(temp_name$n2,temp_name$n1,sep = '-') %>%
Replace(from = ' - ',to='-')#正则表达式替换
names(dif) <- difname
difL <- multcompLetters(dif)
labels <- data.frame(difL['Letters'], stringsAsFactors = FALSE)
labels$compare = rownames(labels)
labels$type <- i

mean_sd <- merge(aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=sd),
aggregate(sub_dat[['value']],by=list(sub_dat[,'g1']),FUN=mean),by="Group.1"
)
names(mean_sd) <- c('compare','std','mean')

a <- rbind(a,merge(mean_sd,labels,by='compare'))
}
names(a) <- c(compare,'std','mean','Letters',group)
return(a)
}

alpha多样性在不同处理下的差别

运行，这里拿alpha多样性测试，看不同alpha多样性在不同处理下的差别。

#df1 <- kruskalmc_compare1(data=dat,group='alpha',compare='Group',value='v')
df1 <- kruskalmc_compare1(dat,'alpha','Group','v')
df1######字母是反着标的(a<b<c)(测试时的字母顺序，下同，可能在跑你自己的数据时会出现反着标的可能，不过不影响使用，可在AI或PDF编辑器中调)

在此可以查看各个α多样性下不同处理间的多组比较字母标注结果，这也是本脚本的亮点之一

数据量很大的情况下，可以直接查看差异情况，不用一个个的做出图再点开看，很是方便。

p1 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
geom_text(data=df1,aes(x=Group,y=mean+1.3*std,label=Letters))+
facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p1

本图一张即可包含所有数据情况，方便查看

#df2 <- PMCMR_compare1(data=dat,group='alpha',compare='Group',value='v')
df2 <- PMCMR_compare1(dat,'alpha','Group','v')
df2########字母标正着标(a>b>c)

p2 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
geom_text(data=df2,aes(x=Group,y=mean+1.3*std,label=Letters))+
facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p2

#df3 <- nparcomp_compare1(data=dat,group='alpha',compare='Group',value='v')
df3 <- nparcomp_compare1 (dat,'alpha','Group','v')
df3########字母标正着标(a>b>c)

p3 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
geom_text(data=df3,aes(x=Group,y=mean+1.3*std,label=Letters))+
facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p3

簇状柱形图显著性字母标记之分面与分组的图形艺术

只是展示方法(多样性数据可能不适合做这个图，但是纵轴标度一致的数据且有需求要做，适合做这个图)。

颜色设置

library(RColorBrewer)
library("scales")
#自己选颜色
display.brewer.all()

#选取颜色
display.brewer.pal(9,"Set1")

#可以直接选取
yanse=brewer.pal(9,"Set1")
#也可以显示出色号后
show_col(brewer.pal(9,"Set1"))

#再挑选喜欢的色号
yanse=c("#E41A1C","#377EB8","#4DAF4A",
"#984EA3","#FF7F00","#FFFF33",
"#A65628","#F781BF","#999999")

分面的图形艺术

数据用的df2

df1$alpha <- factor(df2$alpha,levels = c("richness","chao1","ACE","shannon","simpson","invsimpson"))
H = max(df2$mean)*1.2
#可以根据需求调整分面及分组顺序
p4 = ggplot(df2,aes(x =Group,y = mean,fill = Group))+
facet_wrap(~alpha,nrow=1)+
geom_bar(stat = "identity" ,width=1, position = position_dodge(0.9)) +
geom_text(aes(label = Letters),size=7,vjust=-1.2,position = position_dodge(0.9))+
geom_errorbar(aes(ymin = mean - std/3^0.5,ymax = mean+std/3^0.5),
position = position_dodge(0.9), width = .1,color="red",lwd=1)+
labs(x='Group',y='AlphaDiv')+
scale_y_continuous(expand = c(0,0),limits = c(0,H))+
scale_fill_manual(values = yanse, guide = guide_legend(title = NULL))+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p4

分组的图形艺术

数据用的df2

p5 = ggplot(df2,aes(x =alpha,y = mean,fill = Group))+
geom_bar(stat = "identity", width=0.9, position = position_dodge(0.9)) +
geom_text(aes(label = Letters),size=7,vjust=-1.2,position = position_dodge(0.9))+
geom_errorbar(aes(ymin = mean - std/3^0.5,ymax = mean+std/3^0.5),
position = position_dodge(0.9), width = .1,color="red",lwd=1)+
labs(x='Group',y='AlphaDiv')+
scale_y_continuous(expand = c(0,0),limits = c(0,H))+
scale_fill_manual(values = yanse, guide = guide_legend(title = NULL))+
ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p5

Output figure width and height
Letter纸图片尺寸为单栏89 mm，双栏183 mm，页面最宽为247 mm
推荐比例16：10，即半版89 mm x 56 mm; 183 mm x 114 mm

# ggsave("./alpha1.pdf", p1, width = 350, height = 200, units = "mm")
# ggsave("./alpha2.pdf", p2, width = 350, height = 200, units = "mm")
# ggsave("./alpha3.pdf", p3, width = 350, height = 200, units = "mm")
# ggsave("./alpha4.pdf", p4, width = 600, height = 200, units = "mm")
# ggsave("./alpha5.pdf", p5, width = 600, height = 200, units = "mm")

猜你想学

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

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

参考资料

《R语言统计分析与应用》

浅析R语言单因素方差分析中的多重比较

在R语言上做Kruskal–Wallis秩和检验

Wilcoxon, Friedman…扒一扒非参数检验名称中的牛人

作者：flyfly、旭日阳光

赞 (0)

相关推荐

  • 大样本量多分组表达量矩阵分析你难道没想到单细胞吗

    前面我们演示了:泛癌分析时候的大样本量多分组建议选择tSNE而不是PCA,整合全部的33种癌症的仅仅是蛋白质编码基因的表达量矩阵,进行降维聚类分群可以看到并不是严格的各个癌症泾渭分明. 其中很明显乳腺 ...

  • date_sub函数

    作用:从某日期减去指定的时间间隔后的日期 用法:date_sub(date,  interval expr type); date 参数是合法的日期表达式.expr 参数是希望添加的时间间隔. sel ...

  • 【R分享|实战】科白君教你相关性分析及其绘图

    "R实战"专题·第7篇   编辑 | 科白维尼   1898字 | 4分钟阅读 本期推送内容 相关性分析通常用来定量描述两个变量之间的联系,正相关?负相关?不存在相关?等.常见相关 ...

  • 怎么样才能正确的学习生信分析呢？—从学徒做起

    昨天的我可以为你做些什么好像阅读量不高,不过效果还是蛮显著的,跟部分粉丝聊了聊,希望对他们有帮助吧! 我们的生信故事会一直在征稿,详见:生信故事会栏目征稿启事 下面分享一下最新投稿,讲述跨专业转行学习 ...

  • 泛癌水平的批量生存分析

    肿瘤免疫微环境我们讲了很多了,目录是: estimate的两个打分值本质上就是两个基因集的ssGSEA分析 针对TCGA数据库全部的癌症的表达量矩阵批量运行estimate 不同癌症内部按照estim ...

  • HNSCC数据分析-GSE2379-GPL830-GPL91

    五月份的学徒专注于GEO数据库里面的表达量芯片数据处理,主要的难点是表达量矩阵获取和探针的基因名字转换,合理的分组后就是标准的差异分析,富集分析.主要是参考我八年前的笔记: 解读GEO数据存放规律及下 ...

  • 关于批次效应矫正后出现负值

    学徒已经陆续出师,是时候把生信技能树的舞台交给后辈了! 下面是YuanSH的分享 首先要了解一下什么叫批次效应 那么如何解决批次效应呢? limma 包中 removeBatchEffect 函数中出 ...

  • 三阴性乳腺癌表达矩阵探索笔记之差异性分析

    以第一个基因为例,根据group_list来检验在分组之间是否存在差异 load(file='step1-output.RData') dat[1:4,1:4] table(group_list) # ...

  • Seurat包的findmarkers函数只能根据划分好的亚群进行差异分析吗

    结业考核20题:https://mp.weixin.qq.com/s/lpoHhZqi-_ASUaIfpnX96w 课程配套资料文档在:https://docs.qq.com/doc/DT2NwV0F ...

  • R语言分层线性模型案例

    原文 http://tecdat.cn/?p=3740 有许多分层数据的例子.例如,地理数据通常按层次分组,可能是全球数据,然后按国家和地区分组 .一个生物学的例子是按物种分组的动物或植物的属性,或者 ...