如何优雅地校正多重检验的p值？

      组学派的大家庭越来越热闹，派派也是干劲满满！作为一个公益服务性的平台，为组学战友们排忧解难始终是派派小组至高无上的责任和追求！更加要庆贺的是加入派派的小伙伴越来多，我们有了总群，还有了总群2号（二维码见文末），还有各个分群。我们有了新的LOGO，相信omicsPie会越来越好！

       上周，群里酷爱科研的大佬们对差异性分析结果p value的校正问题掀起了热情洋溢的民主大讨论。既然问题出现了，那么这宝贵的科普机会派派又怎么能够放过，今天就让苯宝宝系统地为大家梳理一下有关False Positive（假阳性）的前生今世，为后入的小伙伴们点亮这盏指坑。。。啊，不，脱坑明灯。

一、什么是False positive（假阳性）？

提及假阳性，不得不从统计假设检验的两类错误讲起。第一类错误（Type I error），就是原假设H0成立（也就是检验的两组之间没有差别），但是统计结果却显示拒绝原假设（可以理解为结果p<0.05），这样就得到了“假”的“阳性”结果。相应的第二类错误（Type II error），则是原假设H0不成立（也就是检验的两组之间有差别，是不同的），但是统计结果却显示接受原假设（可以理解为结果p>0.05），这样就得到了“假”的“阴性”结果。

今天，我们主要谈一谈“假阳性”。关于“假阴性”和“把握度”的问题，下次咱们再接着聊。

    

 

二、假阳性是如何产生的？

科学研究中，研究者往往很难得到全部的研究对象的实验数据，所以随机抽样、统计推断的思想贯穿于实验设计的始终。既然样本是从总体中抽出来的，即使是随机化做得非常好，依然难免会存在抽样误差，导致抽样数据的随机波动，不能完全准确地反应总体的数据水平。这种情况下，很显然，即使是从同一个总体中反复抽样，也有机会得到样本之间的差异经统计学检验p<0.05的结果（如图1）。因此，不难得到这样的结论：假阳性结果会随着抽样次数、比较次数的增多而大量出现。这里派派做了个简单的模拟实验，从同一组正态分布总体中随机抽取两组进行t-test，比较1000次累积的阳性结果（p<0.05）次数如图2所示。

另外，如果同时分析的指标增加，假阳性率又会有怎样的变化呢？若只分析一个指标，假阳性率为0.05，那么不犯I类错误的概率为95%；若同时分析两个指标，假阳性率则变为1-（0.95）²=0.0975；若同时分析三个指标，假阳性率则变为1-（0.95）³=0.0143......若同时分析十个指标，假阳性率则变为1-（0.95）¹⁰=0.4013，假阳性率已经超过了40%！

在组学数据的研究过程中，由于统计分析的变量维度较高，往往会进行指数级的假设检验比较，假阳性的结果则更加常见。

图1

图2

三、假阳性的校正方法

（一）传统统计上的多重检验结果校正方法，主要用于与多组间方差分析（ANOVA）相结合的两两比较。简单总结为以下几种：

1.LSD（least significant difference）最小显著差异t检验校正。这种校正方法一般用于事前比较，也就是实验设计时已经确定进行某些组之间的比较（如A和B，C和D），而其他组间不必进行比较。假如所有组间都要比较用LSD的话会增加假阳性的概率。

2. Dunnet-t校正。这种方法适用于多个实验组均数与对照组均数间的比较，也就是指定其中一组如A组为对照组，B、C、D组均与A组比较，但B、C、D组之间不进行比较。

3.SNK（Student-Newman-Keuls）检验。此方法适用于ANOVA之后的多组间两两比较（事后的全局比较，不指限定比较分组），常用于探索性研究。但结果只告诉有无差异，不提供精确p值。常见的结果形式为

分组	1	2
A	mean(A)
B	mean(B)	mean(B)
C		mean(C)

A组和B组在类别上都被定义为1类，说明两组mean值之间的差异经校正后没有统计学意义；同理B、C两组间的差异也没有统计学意义。但A、C两组分别被定义为1类和2类，说明A、C两组间mean值的差异具有统计学意义。

以上三种校正方法，使用的前提是各组均满足正态性和方差齐性，即仅适用于与ANOVA、多重t检验结合使用。而另外两种适用性更强的万能校正方法，也备受推崇：

4.Bonferroni 校正。这种校正方法比较简单：校正检验水准α'=α/m（这里α通常为0.05，m为检验次数），也相当于如果保持α=0.05不变，校正p值为p'=m*p。这种校正方式最为保守严格，也就是得到的阳性结果把握度会很高，但灵敏度比较低，往往不适用于潜在研究目标的筛选。

5.Sidak法校正。既可用于事前比较，也可用于事后比较。校正检验水准α'=1-（1-α）^1/m，m为比较次数，更适用于比较次数较多时使用。

     当然，还有例Tukey法、Duncan法等诸多方法，小编研究得不多就暂不在此枚举了。

（二）更加适用于高维数据的校正方法Benjamini-Hochberg false discovery rate (FDR)，是基于对假阳性发现率的控制来决定p值的阈值。

相对Bonferroni来说，FDR校正更加温和。其目标是在假阳性和假阴性间达到平衡，将假/真阳性的比例控制在一定范围之内。例如，如果检验100次，我们设定的FDR阈值为0.05（5%），那么无论我们得到多少个差异特征变量，这些差异特征变量中出现假阳性的概率将保持在5%之内，这就是控制FDR＜5%。这种校正方法意义明确，更加适用于新的差异物质的发现，例如用FDR<5%的标准筛选到100个差异基因、蛋白或代谢物，那么可以理解为这100个物质里面大约有95个是在两组间存在真实差异。

FDR的计算公式也十分简单：q=p *m/rank，这里的m为p值的总个数，rank为p值从大到小的排序。简单来说，若进行了100次差异分析得到了100个p值（m=100），则只有最小的p值rank=100，校正之后保持不变，而其他较大的p值排序靠前<100，都会由于检验次数的增加而受到惩罚。

当然，FDR还有很多其他的估计方法和计算公式如：SAM法、经验贝叶斯法等，感兴趣的童鞋们可以自行查阅相关文献，如哈医大卫生统计教研室李康教授的一篇中文综述《多重假设检验中FDR的控制与估计方法》中就有更多细节的介绍。

四、FDR校正的R语言实现

能做FDR校正的R packages还是有很多的，派派给大家介绍一个自己常用的“fdrtool”：

1.计算p value，合并为一个向量。

如果有100个p value,则可记录为data = c(p1:p100)

2. 载入程序包

library('fdrtool')

3.计算FDR校正后的p value

    FDR <- fdrtool(data,statistic='pvalue',plot=F)$qval

4.计算结果中FDR的顺序与data中p value的顺序一一对应。

五、p value校正小建议

千万不要迷信p value！千万不要迷信p value！千万不要迷信p value！

无论是p value还是校正后的p value，都只是代表了统计学上的概率，而概率这个东西在一定程度上并不意味着最终结果的真实与否。所以从统计学的角度，派派更倾向于在研究设计之前就确认好研究的主要目标，围绕主要目标进行的合理严谨的科研设计才是提高结果可靠度的最佳途径，而在研究过程中意外发现的阳性结果，只有经过更多的重复试验，或者更有针对性的实验设计重复研究，才可以认为是靠谱的结果。

     同样的，即使发现的差异物质p value不够显著，但呈现了预期的变化趋势，在考虑到样本量可能不足，实验设计的混杂因素控制得不够理想的情况下，建议不要轻易放弃，最好进行下实验研究的把握度评价。若本次实验的把握度不够理想，而从专业角度上该物质发挥效用的可能性甚高，那么也值得重新设计或扩大实验样本量，进行深入的验证和探索。