先前有几篇帖子我们利用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))
网友评论