python怎么做q检验_FDR多重检验校正原理及Python实现⽅
歌曲曲谱
法
数据分析中常碰见多重检验问题(multiple testing).Benjamini于1995年提出⼀种⽅法，通过控制FDR(Fal Discovery Rate)来决定P值的域值。 假设你挑选了R个差异表达的基因，其中有S个是真正有差异表达的，另外有V个其实是没有差异表达的，是假阳性的.实践中希望错误⽐例Q＝V/R平均⽽⾔不能超过某个预先设定的值(⽐如0.05)，在统计学上，这也就等价于控制FDR不能超过5％。
根据Benjamini在他的⽂章中所证明的定理，控制fdr的步骤实际上⾮常简单。
设总共有m个候选基因，每个基因对应的p值从⼩到⼤排列分别是p(1),p(2),...,p(m)，则若想控制fdr不能超过q，则只需找到最⼤的正整数i，使得p(i)<= (i*q)/m.然后，挑选对应p(1),p(2),...,p(i)的基因做为差异表达基因，这样就能从统计学上保证fdr不超过q。
计算⽅法 请参考：
> p
> p
[1] 3e-04 1e-04 2e-02
> p.adjust(p,method="fdr",length(p))
[1] 0.00045 0.00030 0.02000
腾开头的成语> p*length(p)/rank(p)
[1] 0.00045 0.00030 0.02000
> length(p)
[1] 3
> rank(p)
[1] 2 1 3
sort(p)
[1] 1e-04 3e-04 2e-02
1) P-value是 (在H0 = true的情况下)得到和试验数据⼀样极端(或更极端)的统计量的概率.它不是H1发⽣的概率.假定吃苹果的⼀组和不吃苹果的⼀组的差异为D, P-value=0.2的意思是, pure randomly (即H0=true)的情况下，观察到和D⼀样或⽐D更⼤的差异的概率是20%.
2) p-value的本质是控制PFR (fal positive rate), hypothesis test的⽬的是make decision.传统上把⼩概率事件的概率定义为0.05或0.01，但不总是这样.主要根据研究⽬的.在⼀次试验中(注意：是⼀次试验，即single test), 0.05或0.01的cutoff⾜够严格了(想象⼀下,⼀个⼝袋有100个球, 95个⽩的, 5个红的，只让你摸⼀次，你能摸到红的可能性是多⼤?).我刚才强调的是single test，在multiple test中，通常不⽤p-value,⽽采⽤更加严格的q-value.与p-value不同, q-value控制的是FDR (fal discovery rate).
3)举个例⼦.假如有⼀种诊断艾滋病的试剂，试验验证其准确性为99%(每100次诊断就有⼀次fal positive).对于⼀个被检测的⼈(single test)来说，这种准确性够了.但对于医院(multiple test)来说，这种准确性远远不够，因为每诊断10 000个个体，就会有100个⼈被误诊为艾滋病.
咖啡口味
4)总之，如果你很care fal positive, p-value cutoff就要很低.如果你很care fal negative (就是"宁可错杀⼀千，也不能漏掉⼀个"情况), p-value可以适当放松到0.1, 0.2都是可以的.
珍惜这一刻q-value 是随着multipel test⽽产⽣的.在multiple test (⽐如10000次)，如果⽤p-value=0.05去cut.如果有1000次是显著的，那么在这1000中，有10000*0.05=500次是Fal positive.这显然不能接受.太宽松
了.
Bonferroni提出FWER，在上⾯的例⼦中，就是把cutoff设为: 0.05/10000 = 0.000005，这虽然能控制Fal Positive，但这只在极少数情况下有⽤.因为太严格了,⼤量的true alternatives被miss掉了
q-vlaue 实际上是上述两种⽅法的折衷.既能控制FP，有不会miss掉太多的true alternatives.
*************************学会忘记
赞同：)不过我的问题并不是关于FWER，⽽是关于FDR的控制。Benjamini and Hochberg在1995年第⼀次提出了FDR的概念，其出发点就是基于Bonferroni的保守性，并给出了控制FDR的⽅法(这算是FDR控制⽅法的祖师爷了)。不过他们的⽅法也有其保守性。所以随后⼈们开始研究更加powerful的⽅法，现有的⽅法有Storey的, Broberg的，Dalmasso的，Guan的，Strimmer的等等等等。Benjamini的⽅法是将FDR控制在⼀个level以下，⽽之后所有的⽅法都在试图精确地估计FDR。所以后来的这些⽅法都要powerful⼀些。不过他们所付出的代价就是robustness。
现有FDR控制⽅法最⼤的弊端在于，他们假设p-value's under the null hypothesis是(1)independent(2)following uniform (0,1) distribution.这两点假设从实际观察到的数据来看经常是不合理的，尤其是第⼆点。(顺便提⼀句，Storey和Leek在07年的PLOS Genetics发表了⼀篇⽂章专门解决第⼆个假设的合理性问题，很⽜，有兴趣可以看⼀下)
*******************
接下来我们以T检验结合FDR校正为例，提供Python语⾔下实现的代码
⽣成结果⽂件如表所⽰：
gene
p value
FDR
PMM2
0.006586
0.014778
ZC3H13
0.006281
0.013361
SP1
0.006837
0.016128
左肺下叶纤维灶ABHD14B
0.006727
招商推介会0.015471
CPEB3
皮蛋功效
0.003265
0.00473
BAK1
0.008951
0.033957
第⼆列为T检验后的P值，第三列为BH校正后的P值。
在R语⾔中有专门做FDR校正的函数，因此实现起来相对⽅便很多，实现代码如下data(exp.arr,design)
output
