p.value和FDR
2019-05-08 本文已影响121人
Juan_NF
最近用limma作差异分析，接触p.value和adj.p.value比较多，今天就重点解释下p.value和FDR吧：
在topTable函数的结果里，我们都会看到p.value和adj.p.value：
不得不提假设检验:
其实前两天的简书里提到了这个概念https://www.jianshu.com/p/eede4ea05f59
这里不再赘述，只说下对H0和H1的一个较好的解释：
- H0:差别是由抽样误差所致；
- H1:差别不是由抽样误差所致，即并不来自同一总体；
最终，基于p.value进行H0和H1的抉择，那么它真的很重要，这意味着，得p值者得天下吗？
不不不
topTable里有adjust.method参数，而这个参数就是针对p.value的：
如此重要的p.value也是需要校正的，我们看下，究竟是如何校正的：
p.adjust(p, method = p.adjust.methods, n = length(p))
> head(tT)
logFC AveExpr t P.Value adj.P.Val B
ENSG00000129159.6 5.781884 0.2374425 13.54654 7.193427e-25 6.484393e-21 46.00225
ENSG00000139445.13 6.651408 -1.2961424 13.51254 8.523124e-25 6.484393e-21 45.74851
ENSG00000175175.4 4.901581 1.1389013 13.09707 6.834505e-24 3.466461e-20 43.80622
ENSG00000072832.10 4.746913 4.7398309 12.76099 3.726586e-23 1.098270e-19 42.10851
ENSG00000224621.1 -4.624916 -2.9418765 -12.79910 3.072963e-23 1.098270e-19 41.98058
ENSG00000267746.1 5.580776 -4.9319911 12.73131 4.330718e-23 1.098270e-19 41.74255
> dim(tT)
[1] 15216 6
> head(p.adjust(tT$P.Value,method='BH'))
[1] 6.484393e-21 6.484393e-21 3.466461e-20 1.098270e-19 1.098270e-19 1.098270e-19
> p.adjust(tT$P.Value[1],method='BH')
[1] 7.193427e-25
代码解释：
- tT是limma中topTable函数的结果；
- 将所有p.value直接用p.adjust中的‘BH’方法进行校正，head展示前六个结果，可以看出得到的结果与topTable一致；
- 仅将第一个p.value用p.adjust中的‘BH’方法进行校正，得到的结果其实与p.value一致；
综上：
- 在多重检验的时候，需要对p值进行校正；
- FDR（Benjamini and Hochberg（BH））是p值的校正方法之一；（所以，统计学家对统计方法抱有的态度比较客观，认为其本身并不完美，需要进行校正）
FDR计算方法：
- 总计m次检验的结果按由小到大进行排序，k为其中一次检验结果的P值所对应的排名;
- 找到符合原始阈值α的最大的k值，满足P(k)<=α*k/m，认为排名从1到k的所有检验存在显著差异，并计算对应的q值公式为q = p*(m/k)，q值即FDR;
- 若后一个FDR小于前一个FDR，则以前一位的计算结果为FDR，所以会在结果中看到相同的数值；
【参考内容】
1.https://scientistseessquirrel.wordpress.com/2015/02/09/in-defence-of-the-p-value/
2.https://zhuanlan.zhihu.com/p/51546651
以上内容有什么写得不对的，烦请指正，在此谢过。
课程分享
生信技能树全球公益巡讲
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小时生信工程师教学视频合辑
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招学徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）
