R-ggplot2-箱圖系列（1） basic - 簡書 (jianshu.com)
R-ggplot2-箱圖系列（2） 注釋P值與組間比較 - 簡書 (jianshu.com)
組間比較分析時可能會涉及到以下的分析情況：
1、兩組間比較：（1）選擇有參法還是無參法；（2）能否進行配對比較
2、多組間比較：（1）多組間兩兩比較/多組整體比較(方差分析)
ggpubr包提供了組間比較的分析函數(shù)與可視化函數(shù)，主要參考自http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/76-add-p-values-and-significance-levels-to-ggplots/
0、加載包與示例數(shù)據(jù)
library(ggpubr)
library(patchwork)
#組別名最好是字符型；如果是數(shù)值類型，最好轉(zhuǎn)為因子化
ToothGrowth$dose = factor(ToothGrowth$dose)
summary(ToothGrowth)
# len supp dose
# Min. : 4.20 OJ:30 0.5:20
# 1st Qu.:13.07 VC:30 1 :20
# Median :19.25 2 :20
# Mean :18.81
# 3rd Qu.:25.27
# Max. :33.90
head(ToothGrowth)
# len supp dose
# 1 4.2 VC 0.5
# 2 11.5 VC 0.5
# 3 7.3 VC 0.5
# 4 5.8 VC 0.5
# 5 6.4 VC 0.5
# 6 10.0 VC 0.5
1、組間分析函數(shù)
ggpubr::compare_means()
- 兩組的情況
# （1）Wilcoxon
compare_means(len ~ supp, data = ToothGrowth) #default
# # A tibble: 1 x 8
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 len OJ VC 0.0645 0.064 0.064 ns Wilcoxon
# （2）t檢驗
compare_means(len ~ supp, data = ToothGrowth,
method = "t.test")
# 1 len OJ VC 0.0606 0.061 0.061 ns T-test
# （3） 修改adjP的計算方法
compare_means(len ~ supp, data = ToothGrowth,
p.adjust.method = "BH") #default "holm"
# 1 len OJ VC 0.0645 0.064 0.064 ns Wilcoxon
# （4）考慮其它變量的影響
compare_means(len ~ supp, data = ToothGrowth,
group.by = "dose")
# # A tibble: 3 x 9
# dose .y. group1 group2 p p.adj p.format p.signif method
# <fct> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 0.5 len OJ VC 0.0232 0.046 0.023 * Wilcoxon
# 2 1 len OJ VC 0.00403 0.012 0.004 ** Wilcoxon
# 3 2 len OJ VC 1 1 1.000 ns Wilcoxon
# （5）如果進行配對分析
#那么需要保持每組的樣本排列順序是一致的
compare_means(len ~ supp, data = ToothGrowth,
paired = T)
# 1 len OJ VC 0.00431 0.0043 0.0043 ** Wilcoxon
# （6）修改標(biāo)簽閾值
compare_means(len ~ supp, data = ToothGrowth,
symnum.args = list(cutpoints = c(0, 0.01, 0.05, 1),
symbols = c("***", "*", "not")))
# 1 len OJ VC 0.0645 0.064 0.064 not Wilcoxon
- 多組的情況
#（1）所有兩兩間比較
compare_means(len ~ dose, data = ToothGrowth)
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 len 0.5 1 0.00000702 0.000014 7.0e-06 **** Wilcoxon
# 2 len 0.5 2 0.0000000841 0.00000025 8.4e-08 **** Wilcoxon
# 3 len 1 2 0.000177 0.00018 0.00018 *** Wilcoxon
# （2）都和0.5的組進行比較
compare_means(len ~ dose, data = ToothGrowth,
ref.group = "0.5")
# .y. group1 group2 p p.adj p.format p.signif method
# <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 len 0.5 1 0.00000702 0.000007 7.0e-06 **** Wilcoxon
# 2 len 0.5 2 0.0000000841 0.00000017 8.4e-08 **** Wilcoxon
# （3）方差分析-有參
compare_means(len ~ dose, data = ToothGrowth,
method = "anova") #有參
# .y. p p.adj p.format p.signif method
# <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 len 9.53e-16 9.5e-16 9.5e-16 **** Anova
# （4）方差分析-無參
compare_means(len ~ dose, data = ToothGrowth,
method = "kruskal.test") #無參
# .y. p p.adj p.format p.signif method
# <chr> <dbl> <dbl> <chr> <chr> <chr>
# 1 len 0.00000000148 0.0000000015 1.5e-09 **** Kruskal-Wallis
2、箱圖可視化
2.1 兩組比較
- (1) 不同比較方法
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
# 配色方案 ?ggboxplot
color = "supp", palette = "aaas",
add = "jitter")
# Add p-value
p1 = p + stat_compare_means() #default Wilcoxon
p2 = p + stat_compare_means(method = "t.test")
p1 + p2
- （2）標(biāo)簽顯示格式
#標(biāo)簽位置
p1 = p + stat_compare_means(label.x.npc = "left",
# label.x = 1.5, label.y = 40
label.y.npc = "top")
#標(biāo)簽內(nèi)容
p2 = p + stat_compare_means(aes(label = ..p.signif..))
#自定義閾值
p3 = p + stat_compare_means(aes(label = ..p.signif..),
symnum.args = list(cutpoints = c(0, 0.01, 0.05, 1),
symbols = c("***", "*", "notsig")),
label.x = 1.5, label.y = 40)
p1 | p2 | p3
- （3）配對分析
# 要確保相同樣本在不同組的排列順序相同
ggpaired(ToothGrowth, x = "supp", y = "len",
color = "supp", palette = "jco",
line.color = "gray", line.size = 0.4) +
stat_compare_means(paired = TRUE)
- （4）考慮其它分組變量的影響
p <- ggboxplot(ToothGrowth, x = "supp", y = "len",
color = "supp", palette = "jco",
#add = "jitter",
facet.by = "dose",
short.panel.labs = F)
p1 = p + stat_compare_means(label = "p.format")
# p + stat_compare_means(label = "p.signif", label.x = 1.5)
p <- ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "supp", palette = "jco")
p2 = p + stat_compare_means(aes(group = supp))
# p + stat_compare_means(aes(group = supp), label = "p.signif")
p1 / p2
2.2 多組比較
- （1）多組間比較可視化時，默認(rèn)是 default 方差分析
p1 = ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means()
p2 = ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(method = "anova")
p1 + p2
- （2）多組間兩兩比較
p1 = ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(comparisons = list( c("0.5", "1"),
c("1", "2"),
c("0.5", "2") ))
p2 = ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(comparisons = list( c("0.5", "1"),
c("1", "2"),
c("0.5", "2") ),
label.y = c(29, 35, 40))+ #指定標(biāo)簽的高度
stat_compare_means(label.y = 45) #添加方差分析結(jié)果
p1 | p2
- （3）直接指定一個參考組
ggboxplot(ToothGrowth, x = "dose", y = "len",
color = "dose", palette = "jco")+
stat_compare_means(method = "anova", label.y = 40)+ # Add global p-value
stat_compare_means(label = "p.signif", method = "t.test",
ref.group = "0.5") # Pairwise comparison against reference
