统计方法的选择(3)--非参数检验
前面探索了如果对数据进行了数据分布检验,发现符合正态分布,并且数据方法齐性,那么就可以采用参数检验的方法进行统计检验。如果数据的分布不符合正态分布或者方差齐性,那么就要采用非参数检验。
学了好久,参数检验和非参数检验当中,参数是什么意思,参数(parameters)指什么,而非参数检验又是怎么把参数去掉。
参数检验主要就是指通过对样本数据的统计参数,比如均值、方差、中位数等这些数值的分析,来估计总体参数,然后比较总体之间是否有差异。而非参数检验就是不去管这些参数,只是通过了解数据的分布状况,来判断总体是否具有差异性。
3.非参数检验
还是贯彻以实例促进认知的个人学习理论,直接用例子来开始。
这次使用使用survminer包中的数据集myeloma先加载包
library(survminer)
data("myeloma")
head(myeloma)[1:3,1:11]
查看这个数据结果
colnames(myeloma)
[1] "molecular_group" "chr1q21_status" "treatment" "event" "time"
[6] "CCND1" "CRIM1" "DEPDC1" "IRF4" "TP53"
[11] "WHSC1"
可以看到行是样本,11列中6到11列是基因名称,前面5列是其他信息,我们想要看看DEPDC1基因在几个分组中的表达情况。
3.1 正态性和方差齐性检验
在分析数据之前,尤其在比较差异之前,先对数据的分布进行检查。在统计方法选择(1)中写到的正态分布采用shapiro.test()函数,方差齐性检验使用 bartlett.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
> library(car)
> leveneTest(DEPDC1~molecular_group,data = myeloma)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 6 6.9894 6.954e-07 ***
249
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
正态性检验和方差齐性检验,只要数据正态分布不齐或者方差不齐那么就应该选择非参数检验的方法。
3.2 两组之间非参数检验
两组之间非参数检验通常选择wilcox.test方法,在统计方法的选择(2)--参数检验使用stat_compare_means()和compare_means()两个函数时发现其默认的方法就是非参数检验的方法,两者比较使用wilcox.test,多组全局比较使用kruskal.test 。
针对上面的数据,我们先进性两两比较
table(myeloma$molecular_group)
compare_means(DEPDC1~molecular_group, data = myeloma, method = "wilcox.test")
结果如下
Cyclin D-1 Cyclin D-2 Hyperdiploid Low bone disease MAF MMSET
22 43 66 31 21 44
Proliferation
29
# 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
因为总共分了7组,所以两两组合,就会很多,这个时候,其实可以将这几个组和DEPDC1的平均值进行比较,函数只需要加入一个参数,在统计方法的选择(2)--参数检验最后也谈及到。
方法如下:
compare_means(DEPDC1~molecular_group, data = myeloma,
ref.group = ".all.",
method = "wilcox.test")
结果:
# A tibble: 7 x 8
.y. group1 group2 p p.adj p.format p.signif method
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
1 DEPDC1 .all. Cyclin D-1 0.234 0.94 0.23417 ns Wilcoxon
2 DEPDC1 .all. Cyclin D-2 0.790 1 0.78957 ns Wilcoxon
3 DEPDC1 .all. Hyperdiploid 0.000639 0.0038 0.00064 *** Wilcoxon
4 DEPDC1 .all. Low bone disease 0.00457 0.023 0.00457 ** Wilcoxon
5 DEPDC1 .all. MAF 0.399 1 0.39923 ns Wilcoxon
6 DEPDC1 .all. MMSET 0.383 1 0.38321 ns Wilcoxon
7 DEPDC1 .all. Proliferation 0.000000155 0.0000011 1.5e-07 **** Wilcoxon
可以看到两两比较的差异的结果及各组与平均值之间的比较。
3.2 多组之间的差异比较
其实按照上面的思路,结合上一篇,多组之间只是在函数中选择相应的方法。将上面的进行可视化,并选用kruskal.test方法对多组进行比较,代码如下
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,na.rm = T), linetype=2)+# Add horizontal line at base mean
stat_compare_means(method = "kruskal.test",
label.y = 2000
)+ # Add global kruskal.test p-value
stat_compare_means(label = "p.signif", method = "wilcox.test", ref.group = ".all.")# Pairwise comparison against all
根据统计方法选择的那幅图,只剩下一个问题,就是多组比较之后,如果有差异,那么两两比较如何进行。还是上一篇中提到的,compare_means()和stat_compare_means()这两个函数默认的会对两两之间进行差异比较,更严格的来讲,这样其实会增加犯Ⅰ类错误的概率,所以,统计方面的专家又进行了事后两两比较方法的研究。
写在后面的话:这一篇其实就是上一篇很简单的延续,之所以独立出来,就是为了区分参数检验和非参数检验,而这一篇当中用到的函数以及作图所用的函数都是compare_means()和stat_compare_means()这两个函数。其实针对非参数检验,有专门的两个函数kruskal.test()和wilcox.test(),只不过整合到这个stat_compare_means()函数中,在函数参数选择上选择相应的函数就好。
参考文章及书籍
R语言统计分析与机器学习-薛震
R语言统计分析与应用-汪海波
白话统计-冯国双
一张图说明统计方法的选择
R语言中方差齐性检验丨数析学院
方差分析与R实现
统计方法如何选以及全代码作图实现