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

博文

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

已有 7104 次阅读 2021-8-15 19:45 |个人分类:R|系统分类:科研笔记

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

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

非参数检验的应用

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

测试数据及代码链接:https://pan.baidu.com/s/1s0IX2eVe34jL2DViWCvZfQ
提取码:dunw

链接失败,请回复”多组比较”,获得更新链接或永久链接。

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

参数检验是在已知总体分布的条件下(一般要求总体服从正态分布)对一些主要的参数(如均值、百分数、方差、相关系数等)进行的检验,有时还要求某些总体参数满足一定的条件。如独立样本的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编辑器中调)

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

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

image

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

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

image

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

image

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

image

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

image

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

image

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

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

颜色设置

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

image

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

image

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

image

#再挑选喜欢的色号
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

image

分组的图形艺术

数据用的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

image

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、旭日阳光

责编:马腾飞 南京农业大学

审核:刘永鑫 中科院遗传发育所



https://blog.sciencenet.cn/blog-3334560-1299908.html

上一篇:Microbiome:生态中心张丽梅组-植物发育时期驱动玉米微生物组生态角色的分化
下一篇:Nature:分离到一种位于原核生物-真核生物“交界”的古菌
收藏 IP: 210.75.224.*| 热度|

1 刘强

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

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

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

GMT+8, 2024-3-28 19:07

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部