[R绘图研究 21.1]差异*和abc的显著性（参数检验）

先前有几篇帖子我们利用R也做过几个基本的显著性检验，尤其是柱状图的学习中，如T检验、Wilcoxon检验、方差分析等。你说起这些复杂的统计学问题，其实我也很头疼，需要理解很久。所以我们也陆陆续续去学习几个常见的统计学检验，尤其是针对特定生信的数据。

对于显著性水平的标注，常见的主要有两种样式，一种是“*”，另一种是“abc”。

*很好理解，一般都是根据P值，*（0.01<=p<0.05）,**（0.001<=p<0.01）,***（p<0.001）。

但在多重比较中，对于“abc”的标注，需要同时结合显著性水平以及均值等信息。一般做法如下：

（1）首先根据均值大小，将各组由高往低排排序，均值最高的组标注为“a”；

（2）将均值最高的组与第二高的组相比，若差异显著，则第二组标注为“b”；若不显著，继续比较其与均值第三高的组的差异；

（3）若均值最高的组与第二高的组不显著，均值第二高的组与第三高的组显著，则第二高的组就同样标注为“a”，第三高的组标注为“b”；若均值最高的组与第二高的组不显著、均值第二高的组与第三高的组不显著，但均值最高的组与第三高的组显著，则第二高的组就标注为“ab”，第三高的组标注为“b”；

（4）然后以标注为“b”的组的均值为标准，以此类推，继续循环往后比较，直到最小均值的组被标记，且比较完毕为止。

在这种模式下，只要两组间达到显著性差异水平，即为一个层级；对于显著性差异到底多大，并不是很侧重。

举个例子，下表是A-F6组数据的统计检验结果。

所以标记的显著性水平和abc如下图所示。

注：如果从小到大也是成立的，并不是一定要从大到小。

但是，这种方法只适用于参数检验，啥是参数检验呢？就是假定数据服从某分布（一般为正态分布），通过样本参数的估计量（x±s）对总体参数（μ）进行检验，比如t检验、u检验、方差分析。所以基本只关心均值差异。

某些非参数检验，例如Wilcoxon检验这些，它们的比较方法可以认为是关注的分位数（中位数）差异。有些非参数方法可能不便通过单一的指标（如分位数）评判“高低水平”，如置换检验这种，除非数据分布的高低水平非常明显，否则将很难通过“abc”表示出。

如果采用“*”表示法，将两两分组的显著性全部呈现出，分组少的时候还好说，标注两两差异也不算麻烦；但是当分组多的时候，就不再建议以图形的方式表示了，建议附张统计表，以表格的形式记录两两分组间的差异，如下表所示：

下面，我们就用测试数据来系统的计算一下上面的统计检验方法，以及标注显著性和abc。

如果数据满足方差分析的前提假设，正态性、方差齐性等，那么方差分析肯定就是首选了。我们先进行参数检验，执行单因素 ANOVA 比较整体差异，再执行多重比较（Tukey HSD）查看两两差异。

library(reshape2)

library(multcomp)

library(ggplot2)

library(tidyverse)

data <- read.table("test_data.txt",header=T,sep="\t")

stat_anova <- NULL

abc_list <- NULL

#单因素 ANOVA，整体差异

dat <- data[c("group1","value1")]   #取了group1和value1列

names(dat) <- c('class', 'var')

dat$class <- factor(dat$class)

fit <- aov(var~class, dat)   #group1里面不同类别做anov test

p_value <- summary(fit)[[1]][1,5]

#单因素 ANOVA 结果整理

if (p_value < 0.001) sig <- '***' else if (p_value >= 0.001 & p_value < 0.01) sig <- '**' else if (p_value >= 0.01 & p_value < 0.05) sig <- '*' else sig <- ''

stat_anova <- rbind(stat_anova, c("group1", "value1", p_value, sig))

#Tukey HSD 检验（multcomp 包），多重比较

tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(class = 'Tukey')), sig = p, decreasing = TRUE)

#cld() 自动得到了显著性“abc”水平，提取显著性标记“abc”

sig <- data.frame(tuk$mcletters$Letters, stringsAsFactors = FALSE)

names(sig) <- 'sig'

sig$class <- rownames(sig)

sig$var <- "value1"

#计算均值和标准差

dat_new <- dat  %>% group_by(class) %>% mutate(mean = mean(var), sd =sd(var))

dat_new2 <- unique(dat_new[,c(1,3,4)])

dat_new2$class <- as.character(dat_new2$class)

dat_new2$group <- "group1"

#合并均值，标准差和统计检验结果

dat_new2 <- merge(dat_new2, sig, by = 'class')

dat_new2 <-dat_new2[c(4, 6, 1, 2, 3, 5)]

#画出选的group1和value1的统计检验abc结果

ggplot(data = dat_new2, aes(x = class, y = mean)) +

geom_col(aes(fill = class), color = NA, show.legend = FALSE) +

geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +

geom_text(aes(label = sig, y = mean + sd + 0.3*(mean+sd)),size=10) +

labs(x = '', y = '')+

theme(text = element_text(size = 20))

根据上面的演示，我们可以写2个for循环，外圈来循环group1和group2，内圈循环value，计算分别计算2个group和每个value的统计检验。

groups<-c('group1', 'group2')

values<-c('value1', 'value2', 'value3', 'value4')

stat_anova <- NULL

data_list <- NULL

for(i in groups)

{

  for (j in values)

  { 

    dat <- data[c(i,j)]

    names(dat) <- c('class', 'var')

    dat$class <- factor(dat$class)

    fit <- aov(var~class, dat)

    p_value <- summary(fit)[[1]][1,5]

    #单因素 ANOVA 结果整理

    if (p_value < 0.001) {

      sig <- '***' } else if(p_value >= 0.001 & p_value < 0.01) {

      sig <- '**'  } else if(p_value >= 0.01 & p_value < 0.05)  {

      sig <- '*'  } else {

      sig <- ''    }

    stat_anova <- rbind(stat_anova, c(i, j, p_value, sig))

    #Tukey HSD 检验（multcomp 包），多重比较

  tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(class = 'Tukey')), sig = p, decreasing = TRUE)

  #cld() 自动得到了显著性“abc”水平，提取显著性标记“abc”

  sig <- data.frame(tuk$mcletters$Letters, stringsAsFactors = FALSE)

  names(sig) <- 'sig'

  sig$class <- rownames(sig)

  sig$var <- j

  dat_new <- dat  %>% group_by(class) %>% mutate(mean = mean(var), sd =sd(var))

  dat_new2 <- unique(dat_new[,c(1,3,4)])

  dat_new2$class <- as.character(dat_new2$class)

  dat_new2$group <- i

  dat_new2 <- merge(dat_new2, sig, by = 'class')

  dat_new2 <-dat_new2[c(4, 6, 1, 2, 3, 5)]

  data_list <- rbind(data_list, dat_new2)

  }

}

每组整体pvalue结果。

每组abc结果：

#我们来分面展示每个group的结果：

ggplot(data = data_list, aes(x = class, y = mean)) +

geom_col(aes(fill = class), color = NA, show.legend = FALSE) +

geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +

facet_grid(var~group, scale = 'free' , space = 'free_x') +

geom_text(aes(label = sig, y = mean + sd + 0.3*(mean+sd)),size=5) +

labs(x = '', y = '')+

theme(text = element_text(size = 20))

我们还可以用box来展示最终的结果。

dat_new3 <- melt(data[c(groups, values)], id = groups)

dat_new3 <- melt(dat_new3, id = c('variable', 'value'))

dat_new3 <- dat_new3[c(3, 1, 4, 2)]

names(dat_new3) <- c('group', 'var', 'class', 'value')

value_max <- aggregate(dat_new3$value, by = list(dat_new3$group, dat_new3$var), FUN = max)

names(value_max) <- c('group', 'var', 'value')

value_max <- merge(value_max, data_list[c(-4, -5)], by = c('group', 'var'))

ggplot(data = dat_new3, aes(x = class, y = value)) +

geom_boxplot(aes(fill = class), show.legend = FALSE) +

facet_grid(var~group, scale = 'free' , space = 'free_x') +

geom_text(data = value_max, aes(label = sig, y = value + 0.3*value),size=5) +

labs(x = '', y = '')+

theme(text = element_text(size = 20))