墙裂推荐！统计方法如何选以及全代码作图实现。

墙裂推荐！统计方法如何选以及全代码作图实现。 by 果子学生信

果子推荐

这一篇关于统计的帖子，我要强烈推荐，
他来自一位线下学员的投稿，这位朋友自我迭代的速度快到惊人。

文中讲到了科研绘图时统计方法的选择，如何在图上增加p值，以及全方位的代码实现。
本文的灵感来自于那张统计神图
我觉得数据挖掘的小伙伴们，可以都来参考借鉴一下，至少可以把论文中的方法选对。

到目前为止，我还是没能攻破统计大关，但是只要不下牌桌，就有赢的可能。
在生信数据挖掘的淘金路上，总得有人卖水吧，而我们就是其中一员，但我们不生产水。

以下是正文

看完了MC上那篇m6A胃癌的文章，大家是否对于里面两组或多组之间比较的统计学知识一脸懵，什么时候选择Kruskal-Wallis test、Wilcoxon test、ANOVA test、T-test等等，我们先看看下面几张图吧：

图1

图2

图3

理论讲解

下面开始正式捋一下，首先我们来看一张非常出名的神图

同样地，我们在R中常用的比较方法主要有下面几种：

由此，
当我们拿到一个要进行比较的数据时，应该先进行：正态分布检验和方差齐性检验：

library(multcomp)
#加载数据
data("cholesterol") 
head(cholesterol)

#正态分布检验：用shapiro.test()
shapiro.test(cholesterol$response)    

#方差齐性检验：用bartlett.test()或者leveneTest()
bartlett.test(response~trt,data = cholesterol)   #巴雷特检验

library(car)  #leveneTest()属于car包
leveneTest(cholesterol$response~cholesterol$trt)  #列文检验

当数据服从正态分布，且方差齐时：
两组比较用T检验，R语言中是t.test()函数。
多组比较用单因素方差分析，R中是aov()或anova()函数。

当数据非正太分布或者方式不齐时：
两组的比较，用Mann-Whitney检验（也称Wilcoxon test），在R中是wilcox.test()函数。
多组的比较用Kruskal-Wallis检验，在R中是kruskal.test()函数。

现在，我们再回头看开头三张图:

1.图1，是一个基因的表达量在正常和肿瘤两组的比较，只是作者把各个基因绘制在同一张图罢了，其本质还是基因在正常和肿瘤两组的比较。由于基因的表达量（即是y轴代表的数据）多为非正太分布或者方式不齐的（可做正态分布检验和方差齐性检验），而且compare_means()默认的检验方法是Wilcoxon test，故文章中作者虽然没有说图1用的是什么检验，但是我们可以推知用的是Wilcoxon test。

2.图2，是一种细胞的浸润量（Immune Infiltration，我们做过Cibersort，知道这是一种怎样的数据）在A、B、C三组的比较，属于多组，只是作者把各种细胞绘制在同一张图罢了，其本质还是细胞的浸润量在三组的比较。虽然文章中作者没说，但是Immune Infiltration（即是y轴代表的数据）是非正太分布或者方式不齐的数据（可做正态分布检验和方差齐性检验），多组比较，用Kruskal-Wallis检验。

3.图3B，一种细胞的Enrichment score在在A、B、C三组的比较，属于多组，只是作者把各种细胞绘制在同一张图罢了，其本质还是细胞的Enrichment score在三组的比较。文章中作者说用的是the one-way ANOVA test，而the one-way ANOVA test是包括anova()和kruskal.test()的，这样 score多为非参数的（可做正态分布检验和方差齐性检验），所以推测用的统计方法是Kruskal-Wallis。

至此，大家应该理解了什么情况选择四种中的哪种了，那么问题来了，啥是那张神图最后说的事后两两比较（事后检验）呢。我们知道，如果分组变量中包含两个以上（即是多组）的水平，那么R会自动进行pairwise test（两两比较）,stat_compare_means()函数默认方法为"wilcox.test"，但这个自动两两比较和事后检验的两两比较并非一回事。

果子老师帖子说过：

我们在科研中用到统计学，80%的情况都是比较差异。而上面这张神图中给出了多组数据求差异的流程图，我看过后，那些以前的只言片语，零星记忆都串联在了一起。
1.我们手上的数据，要首先判断是否是正态分布，方差齐不齐
2.假如是正态分布，而且方差齐，就用参数检验，此时要看究竟有几组数据
3.如果是两组，这种情况很常见，就用t检验
4.如果多组，比如基因在癌症中四个亚型中表达，要用方差分析。
5.如果没有差异就算了，如果有差异，还要进行事后检验( post-hoc test )：即是事后两两比较。
6.事后检验，如果是4组：
方法一，各组互相之间进行事后两两比较，理论上事后两两比较需要6次（有些类似下面的箱线图，会发现有6道线）。
方法二：以其中一个组作为对照，其他各组与其比较，只要3次，这两种情况用的统计方法也不太一样。
7.假如一开始数据不是正态分布，或者方差不齐，那么可以用数据变换的方式把数据变成正态方差齐性数据，也可以用非参数检验。
8.非参数检验，如果数据只有两组，用Mann-Whitney U test（Wilcoxon test），这个在文章中比较常见。
9.如果是多组，要用到Kruskal-walls检验，如果有差异，也需要进行事后检验( post-hoc test )，而且也有一系列方法可供选择。

（全局比较即是4组比较用的是Kruskal-walls检验，而这里的两两比较实际上用的是Wilcoxon test，并不就是事后检验两两比较）

代码实战

讲了那么多，下面我们一起来实战看看怎么在R中操作这些检验。

如何添加`p-value`

主要利用ggpubr包中的两个函数：

compare_means():可以进行一组或多组间的比较
stat_compare_mean():自动添加p-value、显著性标记到ggplot图中
①compare_means()函数
该函数主要用用法如下：

compare_means(formula, data, method = "wilcox.test", paired = FALSE,
  group.by = NULL, ref.group = NULL, ...)

注释：

formula:形如x~group，其中x是数值型变量，group是因子，可以是一个或者多个组
data：数据集
method:比较的方法，默认为"wilcox.test", 其他可选方法为："t.test"、"anova"、"kruskal.test"
paired:是否要进行paired test(TRUE or FALSE)，即是配对检验
group_by: 比较时是否要进行分组，下文有演示
ref.group: 是否需要指定参考组

②stat_compare_means()函数

主要用法：

stat_compare_means(mapping = NULL, comparisons = NULL hide.ns = FALSE,
                   label = NULL,  label.x = NULL, label.y = NULL,  ...)

注释：

mapping:由aes()创建的一套美学映射
comparisons:指定需要进行比较以及添加p-value、显著性标记的组
hide.ns:是否要显示显著性标记ns
label:显著性标记的类型，可选项为：p.signif(显著性标记)、p.format(显示p-value值)
label.x、label.y:显著性标签的位置调整
…:其他参数

比较独立的两组

先加载测试数据

#先加载包
library(ggpubr)
#加载数据集ToothGrowth
data("ToothGrowth")
head(ToothGrowth)

统计分析

shapiro.test(ToothGrowth$len) #正态性检验   

library(car)  #leveneTest()属于car包
leveneTest(ToothGrowth$len~ToothGrowth$supp)  #方差齐性检验
#len符合正态分布，且方差齐,两组的比较应用t.test
compare_means(len~supp, data=ToothGrowth)  #R默认是method = "wilcox.test"，两组

结果解释：

.y:测试中使用的y变量
p:p-value
p.adj:调整后的p-value。默认为p.adjust.method="holm"
p.format:四舍五入后的p-value
p.signif:显著性水平，即是p.signif(显著性标记)
method:用于统计检验的方法

绘制箱线图

p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp", 
               palette ="jco", #用jco配图，也可以改为“lancet”，柳叶刀
               add = "jitter")
#添加p-value,运行这行代码将得到下图
p+stat_compare_means()

实际上上面的两组比较应该用t.test，

compare_means(len~supp, data=ToothGrowth,method = "t.test")

也可以用下面的代码改为T检验，并且直接出图

p+stat_compare_means(method = "t.test")

上述p-value值也可以用显著性标记表示，且显著性标记可通过label.x、label.y、hjust及vjust来调整
显著性标记还可以通过aes()映射来更改：

aes(label=..p.format..)或aes(lebel=paste0("p=",..p.format..)):只显示p-value，不显示统计检验方法
aes(label=..p.signif..):仅显示显著性水平
aes(label=paste0(..method..,"\n", "p=",..p.format..)):p-value值与显著性水平分行显示

用代码来展示

#设置method为t.test
a1 <- p+stat_compare_means(aes(label=..p.signif..),
                         method = "t.test",
                         label.x = 1.5, label.y = 40)
#默认是method = "wilcox.test"
a2 <- p+stat_compare_means(aes(label=..p.signif..),
                           label.x = 1.5, label.y = 40)
library(patchwork) #加载拼图R包
a1+a2   #拼图

也可以将上面的标签指定为字符向量，不要映射，只需将p.signif两端的..去掉即可

p+stat_compare_means(label = "p.signif", label.x = 1.5, label.y = 40)

配对分析

比较两个paired sample，参数paired = TRUE

compare_means(len~supp, data=ToothGrowth, paired = TRUE) #默认是method = "wilcox.test"

#compare_means(len~supp, data=ToothGrowth,method = "t.test",paired = TRUE) 实际上这个len数据应该用T检验

#利用ggpaired()进行可视化
ggpaired(ToothGrowth, x="supp", y="len", color = "supp", line.color = "gray", 
         line.size = 0.4, palette = "nat")+ stat_compare_means(paired = TRUE)

多组比较

Global test 全局的

shapiro.test(ToothGrowth$len) #正态性检验   
bartlett.test(len~dose,data = ToothGrowth)   #巴雷特检验(方差齐性检验)
##服从正态分布，且方差齐，该数据的多组比较用anova
compare_means(len~dose, data=ToothGrowth, method = "anova")

###可视化，同理的对于多组的比较，画图是默认的方法也是Kruskal-Wallis
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", 
          palette = "lancet")+
  stat_compare_means()

可以手动改为anova

ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", 
          palette = "lancet")+
  stat_compare_means(method = "anova") #设置method为anova

Pairwise comparisons:
如果分组变量中包含两个以上的水平，那么R会自动进行pairwise test，
compare_means()的帮助文档中有这么一句：

If the grouping variable contains more than two levels, then a pairwise comparison is performed。

stat_compare_means()和compare_means()默认方法为"wilcox.test"。

compare_means(len~dose, data=ToothGrowth) #Pairwise comparisons的默认方法为wilcox.test

stat_compare_means()中可以指定比较哪些组

#先用非参方法试试
my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2"))

ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "lancet")+
  # Add pairwise comparisons p-value 即是各组两两比较的P值，默认wilcox.test
  stat_compare_means(comparisons=my_comparisons)+   
  # Add global p-value 即是全局比较，这里是3组比较的p值，默认Kruskal-Wallis
  stat_compare_means(label.y = 50)

可以通过修改参数label.y来更改标签的位置

ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "lancet")+

  # Add pairwise comparisons p-value 即是各组两两比较的P值，默认wilcox.test
  stat_compare_means(comparisons=my_comparisons,label.y = c(29, 35, 40))+ 

  # Add global p-value 即是全局比较，这里是3组比较的p值，默认Kruskal-Wallis
  stat_compare_means(label.y = 45)

至于通过添加线条来连接比较的两组，这一功能可由包ggsignif实现，
如下，感兴趣的小伙伴可以学习

还可以指定一个组作为参考组，其他各组与该组比较。用参数ref.group

compare_means(len~dose, data=ToothGrowth, ref.group = "0.5",  #以dose=0.5组为参考组 
              method = "t.test" ) #设置比较方法为t.test,这个数据应用t.test

#可视化
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "lancet")+ 

  # Add global p-value，全局比较，设定比较方法为anova,这个数据应用anova
  stat_compare_means(method = "anova", label.y = 40)+

  # Pairwise comparison against reference,增加两两比较的，这个数据应用t.test
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = "0.5")

参考组也可以设置为.all.，即所有的平均值

## 参考组也可以设置为.all.即所有的平均值
compare_means(len~dose, data=ToothGrowth, ref.group = ".all.", method = "t.test")
#可视化
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose", palette = "lancet")+

  # Add global p-value，全局比较，设定比较方法为anova,这个数据应用anova
  stat_compare_means(method = "anova", label.y = 40)+

  #Pairwise comparison against all,这个数据应用t.test
  stat_compare_means(label = "p.signif", method = "t.test", 
                     ref.group = ".all.")

为什么有时候我们需要将ref.group设置为.all.，因为如果存在多组的时候，这样子更方便

library(survminer)#使用survminer包中的数据集myeloma，先加载包
data("myeloma")
head(myeloma)[1:3,1:11]

根据患者的分组(molecular_group)来绘制TP53基因的表达谱，看其在不同组之间的表达量是否存在显著性的差异，可以在7组之间进行比较，但是这样的话组间比较的组合就太多了，因此我们可以将7组中每一组与全部平均值进行比较，看看TP53基因在不同的组中是否高表达还是低表达。

#正态分布检验：用shapiro.test()
shapiro.test(myeloma$TP53)    
#方差齐性检验：用bartlett.test()或者leveneTest()
bartlett.test(TP53~molecular_group,data = myeloma)   #巴雷特检验

可知数据不服从正态分布，方差不齐，故应该用非参方法比较，
wilcox.test（两两）和Kruskal-Wallis（全局）

先用参数方法试试吧

compare_means(TP53~molecular_group, data = myeloma, ref.group = ".all.", method = "t.test")
#可视化TP53基因表达谱
ggboxplot(myeloma, x="molecular_group", y="TP53", 
          color = "molecular_group", add = "jitter", legend="none")+ 
  rotate_x_text(angle = 45)+ 
  geom_hline(yintercept = mean(myeloma$TP53), linetype=2)+# Add horizontal line at base mean 
  stat_compare_means(method = "anova", label.y = 6000)+ # Add global annova p-value 
  stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.")# Pairwise comparison against all

实际上这个数据，应用非参方法比较wilcox.test（两两）和Kruskal-Wallis（全局）

compare_means(TP53~molecular_group, data = myeloma, ref.group = ".all.", method = "wilcox.test")
ggboxplot(myeloma, x="molecular_group", y="TP53", 
          color = "molecular_group", add = "jitter", legend="none")+ 
  rotate_x_text(angle = 45)+ 
  geom_hline(yintercept = mean(myeloma$TP53), linetype=2)+# Add horizontal line at base mean 
  stat_compare_means(method = "kruskal.test", label.y = 6000)+ # Add global kruskal.test p-value 
  stat_compare_means(label = "p.signif", method = "wilcox.test", ref.group = ".all.")# Pairwise comparison against all

从图中可以看出，TP53基因在Hyperdiploid组中显著性地过表达，而在MMSET组中显著性地低表达.

我们也可以将非显著性标记ns去掉，只需要将参数hide.ns=TRUE

ggboxplot(myeloma, x="molecular_group", y="DEPDC1", 
          color = "molecular_group", add = "jitter", legend="none")+
  rotate_x_text(angle = 45)+ 
  geom_hline(yintercept = mean(myeloma$DEPDC1), linetype=2)+# Add horizontal line at base mean 
  stat_compare_means(method = "kruskal.test", label.y = 1600)+ # Add global kruskal.test p-value 
  stat_compare_means(label = "p.signif", method = "wilcox.test", ref.group = ".all.", hide.ns = TRUE)# Pairwise comparison against all

事后检验：

Post-hoc pairwise comparisons（事后两两比较）

那么问题来了，多组比较：对于全局比较，我们知道用的是kruskal.test（非参数）或者anova（参数），同时，stat_compare_means()和compare_means()函数对于含两组以上的水平时，R会自动进行两两比较Pairwise comparisons（不管全局比较有没有发现显著性差异）。compare_means（）的帮助文档中有这么一句话：If the grouping variable contains more than two levels, then a pairwise comparison is performed。stat_compare_means()和compare_means()默认两两比较方法为"wilcox.test"，我们也可以根据数据（为参数），手动改为t.test。

问题：多组比较时，stat_compare_means()和compare_means()函数的两两比较，与那张神图和果子老师所说的事后检验两两比较是否是一回事呢？很明显，用的统计检验方法是不一样的：stat_compare_means()和compare_means()函数用的是wilcox.test(非参数)或者t.test(参数)，不管全局比较有没有发现显著性差异，R自动进行的；而事后检验，是基于多组全局比较发现显著性差异 后才需要进行的，用的统计方法如那张神图所示：

①对于参数：全局检验anova()发现有显著性差异时，事后检验可用：SNK法，Duncan法，Tukey HSD法，Scheffe法等
②对于非参数：全局检验kruskal.test()发现有显著性差异时，事后检验可用：Steel-Dwass检验，Games-Howell检验，Tamhane's T检验，Dunnett‘s T3检验等

下面来在R中体会一下吧：

参数事后检验

(1)Tukey HSD法进行事后检验

data("ToothGrowth")
head(ToothGrowth)
str(ToothGrowth)
ToothGrowth$dose<-factor(ToothGrowth$dose)

#正态性检验 
shapiro.test(ToothGrowth$len) 
#巴雷特检验(方差齐性检验)
bartlett.test(len~dose,data = ToothGrowth)

数据服从正态分布，且方差齐

table(ToothGrowth$dose)

0.5   1   2 
 20  20  20

可知各组样本数相等，都是20，对于各组样本相等的参数，首选Tukey HSD法进行事后检验

#Global test，aa下文反复要用到
aa<-aov(len~dose,data = ToothGrowth)  
summary(aa)

结果是

            Df Sum Sq Mean Sq F value   Pr(>F)    
dose         2   2426    1213   67.42 9.53e-16 ***
Residuals   57   1026      18                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

事后检验，并显示95%置信区间

TukeyHSD(aa,conf.level = 0.95)

结果是,可见两两比较还是有意义的。

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = len ~ dose, data = ToothGrowth)

$dose
        diff       lwr       upr    p adj
1-0.5  9.130  5.901805 12.358195 0.00e+00
2-0.5 15.495 12.266805 18.723195 0.00e+00
2-1    6.365  3.136805  9.593195 4.25e-05

这里的p值是多重假设检验之后的，跟我们直接来比较的不一样，比如之前图中的

library(ggpubr)
my_comparisons <- list(c("0.5", "1"), c("1", "2"), c("0.5", "2"))
ggboxplot(ToothGrowth, x="dose", y="len", color = "dose",palette = "lancet")+
  # Add pairwise comparisons p-value 即是各组两两比较的P值，默认wilcox.test
  stat_compare_means(comparisons=my_comparisons,label.y = c(29, 35, 40))+ 
  # Add global p-value 即是全局比较，这里是3组比较的p值，默认Kruskal-Wallis
  stat_compare_means(label.y = 45)

(2)LSD法进行事后检验
LSD法在R语言中可利用agricolae包中的LSD.test函数实现，其调用格式为：

LSD.test(y, trt, DFerror, MSerror, alpha = 0.05, p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr"),...)

其中y为方差分析对象，
trt为要进行多重比较的分组变量，
p.adj可以选定P值矫正方法。
p.adj="bonferroni"时为Bonferroni法。

把console参数设置为TRUE可以直接展示分析的结果

aa<-aov(len~dose,data = ToothGrowth)
library(agricolae)
aa_output <- LSD.test(aa,"dose", p.adj="none",console=TRUE)

结果如下：

Study: aa ~ "dose"
LSD t Test for len 
Mean Square Error:  17.99605 
dose,  means and individual ( 95 %) CI
       len      std  r       LCL     UCL  Min  Max
0.5 10.605 4.499763 20  8.705503 12.5045  4.2 21.5
1   19.735 4.415436 20 17.835503 21.6345 13.6 27.3
2   26.100 3.774150 20 24.200503 27.9995 18.5 33.9
Alpha: 0.05 ; DF Error: 57
Critical Value of t: 2.002465 
least Significant Difference: 2.686295 
Treatments with the same letter are not significantly different.
       len groups
2   26.100      a
1   19.735      b
0.5 10.605      c

注意：运行结果并不提供P值，只提供归类结果，看最后的a,b,c
字母不同，表示跟别人比有差异。
如果是英文字母者表示没有统计学差异。比如下面的结果

      yield groups
oo 36.90000      a
ff 36.33333      a
cc 24.40000      b
fc 12.86667      c

那么a的两组，oo，和ff就是没有差异的。
但是，ab，ac，bc之间是有差异的。
从图上隐约可以看出来。

(3)SNK法（Student-Newman-Keuls)进行事后检验
SNK法可用agricolae包中的SNK.test()函数实现，其调用格式为：

SNK.test(y, trt, alpha = 0.05,...)

其中y为方差分析对象，trt为要进行多重比较的分组变量

aa<-aov(len~dose,data = ToothGrowth)
library(agricolae)
output <- SNK.test(aa,"dose",console=TRUE)

结果是：

Study: aa ~ "dose"
Student Newman Keuls Test
for len 
Mean Square Error:  17.99605 
dose,  means
       len      std  r  Min  Max
0.5 10.605 4.499763 20  4.2 21.5
1   19.735 4.415436 20 13.6 27.3
2   26.100 3.774150 20 18.5 33.9

Alpha: 0.05 ; DF Error: 57 
Critical Range
       2        3 
2.686295 3.228195 
Means with the same letter are not significantly different.
       len groups
2   26.100      a
1   19.735      b
0.5 10.605      c

输出结果与LSD.test类似。不赘述。

(4)Duncan法进行事后检验
Duncan法可用agricolae包中的duncan.test()函数实现，其调用格式为：

duncan.test(y, trt,...)

其中y为方差分析对象，trt为要进行多重比较的分组变量

aa<-aov(len~dose,data = ToothGrowth)
library(agricolae)
output <- duncan.test(aa,"dose",console=TRUE)

输出结果与LSD.test类似。不赘述。

(5)Scheffe检验进行事后比较
各组样本数相等或不等均可以，但是以各组样本数不相等使用较多。
Scheffe法可用agricolae包中的scheffe.test()函数实现，其调用格式为：

scheffe.test(y, trt, ...)

其中y为方差分析对象，trt为要进行多重比较的分组变量

aa<-aov(len~dose,data = ToothGrowth)
library(agricolae)
output <- scheffe.test(aa,"dose",console=TRUE)

心累，输出结果与LSD.test类似。不赘述。

非参数事后检验

使用survminer包中的数据集myeloma，先加载包

library(survminer)
data("myeloma")
## 正态检验
shapiro.test(myeloma$TP53)    
#方差齐性检验：用bartlett.test()或者leveneTest()
bartlett.test(TP53~molecular_group,data = myeloma)

可知数据不服从正态分布，方差不齐，故应该用非参方法比较
用非参方法进行全局比较检验，多组，用kruskal.test()

kruskal.test(TP53 ~ molecular_group, data=myeloma)

有差异

    Kruskal-Wallis rank sum test

data:  TP53 by molecular_group
Kruskal-Wallis chi-squared = 35.133, df = 6, p-value = 4.062e-06

需要事后检验
举一个例子。
posthoc.kruskal.nemenyi.test:
用法是这样的

library(PMCMR)
posthoc.kruskal.nemenyi.test(x=myeloma$TP53, g=myeloma$molecular_group,dist="Tukey")

然后就开始输出两两比较的结果

    Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
                   with Tukey-Dist approximation for independent samples 

data:  myeloma$TP53 and myeloma$molecular_group 

                 Cyclin D-1 Cyclin D-2 Hyperdiploid Low bone disease MAF    MMSET 
Cyclin D-2       0.9928     -          -            -                -      -     
Hyperdiploid     0.1912     0.0014     -            -                -      -     
Low bone disease 1.0000     0.9971     0.0523       -                -      -     
MAF              0.9009     0.9943     0.0026       0.9215           -      -     
MMSET            0.8121     0.9820     1.9e-05      0.8270           1.0000 -     
Proliferation    0.9312     0.3869     0.8581       0.8221           0.2272 0.0762

P value adjustment method: none

可以看到两两有差异的并不多。
(果子说明:该部分还有很多风流操作，目前已经超出了我的能力范围，在后面我们会专门写一个帖子来记录学习成果)

事后检验的方法非常多，但功能比较一致，仅在个别点或使用场景上有小区别，
当前SPSSAU共提供LSD，Scheffe，Tukey，Bonferroni校正，Tamhane T2总共最为常见使用的五种方法，该五种方法如何选择说明如下表格所示：

呼~~，事后检验终于结束，我们回归正传！！

多个分组变量

按另一个变量进行分组之后进行统计检验，比如按变量dose进行分组：

data("ToothGrowth")
library(ggpubr)
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

可见默认的是Wilcoxon
上面我们已经知道len这个数据应该用t.test,两组比较

compare_means(len~supp, data=ToothGrowth,method = "t.test",,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.00636 0.013  0.0064   **       T-test
2 1     len   OJ     VC     0.00104 0.0031 0.0010   **       T-test
3 2     len   OJ     VC     0.964   0.96   0.9639   ns       T-test

画图展示

p <- ggboxplot(ToothGrowth, x="supp", y="len", color = "supp", 
               palette = "jco", 
               add = "jitter", facet.by = "dose", short.panel.labs = FALSE)#按dose进行分面
#label只绘制p-value值,不在图上展示统计方法
b1=p+stat_compare_means(label = "p.format",method = "t.test")
#label绘制显著性水平，不在图上展示统计方法
b2=p+stat_compare_means(label = "p.signif", label.x = 1.5,method = "t.test")
library(patchwork)
b1+b2  #拼图

将所有箱线图绘制在一个panel中,这个图就是我们常见的,例如m6A胃癌那里的箱线图

p <- ggboxplot(ToothGrowth, x="dose", y="len", color = "supp", palette = "lancet", add = "jitter")
m1=p+stat_compare_means(aes(group=supp)) #默认是Wilcoxon

## 手动改为t.test
p <- ggboxplot(ToothGrowth, x="dose", y="len", color = "supp", palette = "lancet", add = "jitter")
m2=p+stat_compare_means(aes(group=supp),method = "t.test")

m1+m2

在图中只显示p-value,不显示检验检验方法,加label = "p.format"即可

m1=p+stat_compare_means(aes(group=supp),method ="wilcox.test",label = "p.format") #其实默认就是method ="wilcox.test"
m2=p+stat_compare_means(aes(group=supp),method = "t.test",label = "p.format")
m1+m2

显示显著性水平,加label = "p.signif"

m1=p+stat_compare_means(aes(group=supp), label = "p.signif") #默认就是method ="wilcox.test"
m2=p+stat_compare_means(aes(group=supp),method = "t.test", label = "p.signif") #改为method = "t.test"
m1+m2

可进行配对样本检验：paired sample检验

compare_means(len~supp, data=ToothGrowth,group.by = "dose", paired = TRUE)#默认是wilcox.test
#compare_means(len~supp, data=ToothGrowth,method = "t.test",group.by = "dose", paired = TRUE) #这个数据应该用t.test

#可视化
p <- ggpaired(ToothGrowth, x="supp", y="len", color = "supp", 
              palette = "lancet", line.color="gray", 
              line.size=0.4, facet.by = "dose", 
              short.panel.labs = FALSE)#按dose分面
#只显示p-value
p+stat_compare_means(label = "p.format", paired = TRUE) #默认是wilcox.test
#p+stat_compare_means(label = "p.format",method = "t.test", paired = TRUE)

大感谢：

至此，四种比较方法理完了，我们可在熟知这几种统计方法的基础上，绘制各种SCI文章中的箱线图了！
感谢前辈们的学习经验！
感谢生信路上遇到果子老师，带我入门，前行路上给了我很多亦师亦兄的帮助和鼓励！！
谢谢！！

参考资料：

(注：资料来源与果子老师发的资料以及网上查找学习及实践总结，感谢前辈)
https://stats.idre.ucla.edu/r/faq/how-can-i-do-post-hoc-pairwise-comparisons-in-r/
https://cran.r-project.org/web/packages/PMCMR/vignettes/PMCMR.pdf
R 语言教程：方差分析与多重比较
https://mp.weixin.qq.com/s/2AHFAWn9wE_B-aYKuDVGdQ
https://www.cnblogs.com/xiaojikuaipao/p/7862990.html
单因素方差分析及事后检验在R中的实现
https://mp.weixin.qq.com/s/LGX-xnquExOJAJOalnr1ig
http://blog.sciencenet.cn/blog-3406804-1190969.html
https://blog.csdn.net/linkequa/article/details/84248536
https://blog.csdn.net/enyayang/article/details/98175441

ixxmu / mp_duty