前面几篇通过参数检验和非参数检验对多组数据进行检验后，发现有差异，那么究竟是哪几个之间有差异，这就涉及到本篇所讲的事后检验或者事后两两检验。真如前面几篇中写的，事后检验和compare_means()和stat_compare_means()这两个函数默认的对两两之间进行差异比较是有区别的。
本篇就事后检验进行R语言的实践和操作。
4. 事后检验
还是把最开始的那幅图贴上，最后一部分就是事后检验的方法选择
统计方法选择
从图中可以看出，参数检验和非参数检验在多组比较有差异后进行两两比较方法的选择。
4.1 参数检验事后检验
数据依然使用统计方法的选择（2）--参数检验中三组小鼠的实验的数据，
先进性正态性和方差齐性检验
#正态性检验
shapiro.test(Example8_1$atp)
#方差齐性检验
bartlett.test(Example8_1$atp ~Example8_1$group)
结果如下
> #正态性检验
> shapiro.test(Example8_1$atp)
Shapiro-Wilk normality test
data: Example8_1$atp
W = 0.94636, p-value = 0.135
> #方差齐性检验
> bartlett.test(Example8_1$atp ~Example8_1$group)
Bartlett test of homogeneity of variances
data: Example8_1$atp by Example8_1$group
Bartlett's K-squared = 5.5965, df = 2, p-value = 0.06092
p值>0.05，符合正态分布和方差齐性
接下来进行方差分析
attach(Example8_1)
fit <- aov(atp~group)
summary(fit)
结果如下
> fit <- aov(atp~group)
> summary(fit)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 119.8 59.92 14.32 5.76e-05 ***
Residuals 27 113.0 4.18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
经过ANOVA方差分析,发现p<0.05，说明全局角度有差异，那么接下来对两两进行比较。
table(group)
group
1 2 3
10 10 10
每个分组的数目都是相等的，事后检验首选Tukey HSD法，并添加conf.level = 0.95置信区间。
4.1.1 Tukey HSD法
> TukeyHSD(fit,conf.level = 0.95)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = atp ~ group)
$group
diff lwr upr p adj
2-1 4.712 2.443877 6.980123 0.0000588
3-1 1.206 -1.062123 3.474123 0.3975104
3-2 -3.506 -5.774123 -1.237877 0.0019270
这里的p值和之前直接比较得来的p值是有区别的，之前直接比较的结果如下
> compare_means(atp~group, data=Example8_1,method = 't.test')
# A tibble: 3 x 8
.y. group1 group2 p p.adj p.format p.signif method
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
1 atp 1 2 0.000411 0.00120 0.00041 *** T-test
2 atp 1 3 0.100 0.1 0.10010 ns T-test
3 atp 2 3 0.00325 0.0065 0.00325 ** T-test
可以发现通过事后比较的方法得到的p值，显著差异的是明显小于直接比较的，这也就可以解释通过直接比较会增加犯Ⅰ类错误的概率。
4.1.2 LSD法进行事后检验
还可以使用agricolae包中的LSD.test函数实现参数事后检验
依然是上面的数据，操作如下
library(agricolae)
LSD.test(fit,"group", p.adj="none",console=T)
console=T可以直接展示结果
Study: fit ~ "group"
LSD t Test for atp
Mean Square Error: 4.184117
group, means and individual ( 95 %) CI
atp std r LCL UCL Min Max
1 8.043 1.808075 10 6.715779 9.370221 5.78 11.73
2 12.755 2.790047 10 11.427779 14.082221 6.94 17.72
3 9.249 1.224277 10 7.921779 10.576221 7.19 11.26
Alpha: 0.05 ; DF Error: 27
Critical Value of t: 2.051831
least Significant Difference: 1.876975
Treatments with the same letter are not significantly different.
atp groups
2 12.755 a
3 9.249 b
1 8.043 b
最后的结果只是将各组分组，不同字母表示不同组别，具有差异性，结果可以看出3组和1组不具有差异性，和之前做的图也是符合的。
4.1.3 其他方法
参数事后检验还有其他方法，比如SNK法（Student-Newman-Keuls)事后检验，Duncan事后检验，Scheffe事后比较，这几个方法的检验都在agricolae包中，展现形式不外乎以上两种方式。
4.2 非参数事后检验
非参数事后检验的数据还是使用统计方法的选择（3）--非参数检验中的数据。
先加载数据并查看数据
library(survminer)#使用survminer包中的数据集myeloma，先加载包
data("myeloma")
head(myeloma)[1:3,1:11]
colnames(myeloma)
主要探索DEPDC1基因的表达情况
先进性正态性和方差齐性检验
> #正态分布检验：用shapiro.test()
> shapiro.test(myeloma$DEPDC1)
Shapiro-Wilk normality test
data: myeloma$DEPDC1
W = 0.84077, p-value = 1.606e-15
> #方差齐性检验：用bartlett.test()或者leveneTest()
> bartlett.test(DEPDC1~molecular_group,data = myeloma) #巴雷特检验
Bartlett test of homogeneity of variances
data: DEPDC1 by molecular_group
Bartlett's K-squared = 89.032, df = 6, p-value < 2.2e-16
结果如之前所写，正态性检验和方差齐性检验，p值都小于0.05，所以不符合正态性，方差不齐，采用非参数检验。对于多组非参数比较，采用用用kruskal.test()
kruskal.test(DEPDC1 ~ molecular_group, data=myeloma)
结果如下
Kruskal-Wallis rank sum test
data: DEPDC1 by molecular_group
Kruskal-Wallis chi-squared = 57.611, df = 6, p-value = 1.374e-10
p值小于0.05，多组之间有差异，那么进行组间比较，先看看之前直接比较的结果
# A tibble: 21 x 8
.y. group1 group2 p p.adj p.format p.signif method
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
1 DEPDC1 Cyclin D-1 Cyclin D-2 0.267 1 0.2667 ns Wilcoxon
2 DEPDC1 Cyclin D-1 Hyperdiploid 0.00356 0.053 0.0036 ** Wilcoxon
3 DEPDC1 Cyclin D-1 Low bone disease 0.00521 0.065 0.0052 ** Wilcoxon
4 DEPDC1 Cyclin D-1 MAF 0.819 1 0.8193 ns Wilcoxon
5 DEPDC1 Cyclin D-1 MMSET 0.681 1 0.6805 ns Wilcoxon
6 DEPDC1 Cyclin D-1 Proliferation 0.00465 0.065 0.0046 ** Wilcoxon
7 DEPDC1 Cyclin D-2 Hyperdiploid 0.00219 0.035 0.0022 ** Wilcoxon
8 DEPDC1 Cyclin D-2 Low bone disease 0.00474 0.065 0.0047 ** Wilcoxon
9 DEPDC1 Cyclin D-2 MAF 0.563 1 0.5625 ns Wilcoxon
10 DEPDC1 Cyclin D-2 MMSET 0.502 1 0.5016 ns Wilcoxon
# ... with 11 more rows
对于非参数检验事后比较，这里采用posthoc.kruskal.nemenyi.test方法，操作如下
library(PMCMRplus)
posthoc.kruskal.nemenyi.test(x=myeloma$DEPDC1, g=myeloma$molecular_group,dist="Tukey")
结果如下
Pairwise comparisons using Tukey and Kramer (Nemenyi) test
with Tukey-Dist approximation for independent samples
data: myeloma$DEPDC1 and myeloma$molecular_group
Cyclin D-1 Cyclin D-2 Hyperdiploid Low bone disease MAF MMSET
Cyclin D-2 0.98063 - - - - -
Hyperdiploid 0.04443 0.11703 - - - -
Low bone disease 0.06004 0.16783 0.99993 - - -
MAF 0.99998 0.99797 0.11248 0.13002 - -
MMSET 0.99922 0.99931 0.02700 0.05619 1.00000 -
Proliferation 0.09810 0.00083 3.6e-10 2.7e-08 0.05438 0.00401
P value adjustment method: none
Warning message:
In posthoc.kruskal.nemenyi.test.default(x = myeloma$DEPDC1, g = myeloma$molecular_group, :
Ties are present, p-values are not corrected.
结果还是有点小复杂，好像没有之前直接比较或者与平均值比较更清晰，并且差异性把分组也加入了，这个还要继续再探索一下。
好了，基本把两组比较，多组比较的整个流程都进行了R实际操作，应该是能解决绝大部分平时遇到的问题，至于还有其他的非参数检验，比如卡方检验，列联表的检验方法可以再写。
天下大事，必作于细。天下难事，必作于易。平时看似都会的东西，不经过实践，永远都是眼睛会了、脑子可能会了，就是手不会。以上的操作，原理可能以前原理都有所了解，但操作的时候，依然会发现有好多不能执行。
无他，唯手熟而。
网友评论