先前有幾篇帖子我們利用R也做過幾個基本的顯著性檢驗(yàn)，尤其是柱狀圖的學(xué)習(xí)中，如T檢驗(yàn)、Wilcoxon檢驗(yàn)、方差分析等。你說起這些復(fù)雜的統(tǒng)計學(xué)問題，其實(shí)我也很頭疼，需要理解很久。所以我們也陸陸續(xù)續(xù)去學(xué)習(xí)幾個常見的統(tǒng)計學(xué)檢驗(yàn)，尤其是針對特定生信的數(shù)據(jù)。
對于顯著性水平的標(biāo)注，常見的主要有兩種樣式，一種是“*”，另一種是“abc”。
*很好理解，一般都是根據(jù)P值，*（0.01<=p<0.05）,**（0.001<=p<0.01）,***（p<0.001）。
但在多重比較中，對于“abc”的標(biāo)注，需要同時結(jié)合顯著性水平以及均值等信息。一般做法如下：
（1）首先根據(jù)均值大小，將各組由高往低排排序，均值最高的組標(biāo)注為“a”；
（2）將均值最高的組與第二高的組相比，若差異顯著，則第二組標(biāo)注為“b”；若不顯著，繼續(xù)比較其與均值第三高的組的差異；
（3）若均值最高的組與第二高的組不顯著，均值第二高的組與第三高的組顯著，則第二高的組就同樣標(biāo)注為“a”，第三高的組標(biāo)注為“b”；若均值最高的組與第二高的組不顯著、均值第二高的組與第三高的組不顯著，但均值最高的組與第三高的組顯著，則第二高的組就標(biāo)注為“ab”，第三高的組標(biāo)注為“b”；
（4）然后以標(biāo)注為“b”的組的均值為標(biāo)準(zhǔn)，以此類推，繼續(xù)循環(huán)往后比較，直到最小均值的組被標(biāo)記，且比較完畢為止。
在這種模式下，只要兩組間達(dá)到顯著性差異水平，即為一個層級；對于顯著性差異到底多大，并不是很側(cè)重。
舉個例子，下表是A-F6組數(shù)據(jù)的統(tǒng)計檢驗(yàn)結(jié)果。
所以標(biāo)記的顯著性水平和abc如下圖所示。
注：如果從小到大也是成立的，并不是一定要從大到小。
但是，這種方法只適用于參數(shù)檢驗(yàn)，啥是參數(shù)檢驗(yàn)?zāi)?？就是假定?shù)據(jù)服從某分布（一般為正態(tài)分布），通過樣本參數(shù)的估計量（x±s）對總體參數(shù)（μ）進(jìn)行檢驗(yàn)，比如t檢驗(yàn)、u檢驗(yàn)、方差分析。所以基本只關(guān)心均值差異。
某些非參數(shù)檢驗(yàn)，例如Wilcoxon檢驗(yàn)這些，它們的比較方法可以認(rèn)為是關(guān)注的分位數(shù)（中位數(shù)）差異。有些非參數(shù)方法可能不便通過單一的指標(biāo)（如分位數(shù)）評判“高低水平”，如置換檢驗(yàn)這種，除非數(shù)據(jù)分布的高低水平非常明顯，否則將很難通過“abc”表示出。
如果采用“*”表示法，將兩兩分組的顯著性全部呈現(xiàn)出，分組少的時候還好說，標(biāo)注兩兩差異也不算麻煩；但是當(dāng)分組多的時候，就不再建議以圖形的方式表示了，建議附張統(tǒng)計表，以表格的形式記錄兩兩分組間的差異，如下表所示：
下面，我們就用測試數(shù)據(jù)來系統(tǒng)的計算一下上面的統(tǒng)計檢驗(yàn)方法，以及標(biāo)注顯著性和abc。
如果數(shù)據(jù)滿足方差分析的前提假設(shè)，正態(tài)性、方差齊性等，那么方差分析肯定就是首選了。我們先進(jìn)行參數(shù)檢驗(yàn)，執(zhí)行單因素 ANOVA 比較整體差異，再執(zhí)行多重比較（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 結(jié)果整理
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 檢驗(yàn)（multcomp 包），多重比較
tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(class = 'Tukey')), sig = p, decreasing = TRUE)
#cld() 自動得到了顯著性“abc”水平，提取顯著性標(biāo)記“abc”
sig <- data.frame(tuk$mcletters$Letters, stringsAsFactors = FALSE)
names(sig) <- 'sig'
sig$class <- rownames(sig)
sig$var <- "value1"
#計算均值和標(biāo)準(zhǔn)差
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"
#合并均值，標(biāo)準(zhǔn)差和統(tǒng)計檢驗(yàn)結(jié)果
dat_new2 <- merge(dat_new2, sig, by = 'class')
dat_new2 <-dat_new2[c(4, 6, 1, 2, 3, 5)]
#畫出選的group1和value1的統(tǒng)計檢驗(yàn)abc結(jié)果
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))
根據(jù)上面的演示，我們可以寫2個for循環(huán)，外圈來循環(huán)group1和group2，內(nèi)圈循環(huán)value，計算分別計算2個group和每個value的統(tǒng)計檢驗(yàn)。
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 結(jié)果整理
? ? 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 檢驗(yàn)（multcomp 包），多重比較
? tuk <- cld(glht(fit, alternative = 'two.sided', linfct = mcp(class = 'Tukey')), sig = p, decreasing = TRUE)
? #cld() 自動得到了顯著性“abc”水平，提取顯著性標(biāo)記“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結(jié)果。
每組abc結(jié)果：
#我們來分面展示每個group的結(jié)果：
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來展示最終的結(jié)果。
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))
