p值校正，FDR（BH法）的实现过程

原理：我们要看下最常用的BH法的论文

做m次无效假设作物的数量

 

那么，被错误地拒绝了的无效假设的比例Q=V/(V+S)=V/R

所谓的FDR值就是Q的期望值，E（Q）=E（V/R）

如果无效假设是正确的（s=0且v=r），FDR值就和FWER（familywise error rate）一样，假设v=0(结合s=0,就是说假设结果一次错误都没有出),Q就为0；如果v>0(出现了假阳性)，那么Q=1，P(v≥1）=E（Q）。这样的结果提示一般的FWER法控制假阳性结果的能力是很弱的；

只要有m=m0，就会出现Q=1，所以E（V/R|R>0）是不能作为控制，另外一种无效假设的控制方式是P（R>0）E（V/R|R>0）.

如果P（i）≤i/m*q(i,第i次假设检验;q是设定阈值、一般为0.05),则拒绝所有无效假设Hi(数学证明请参照原文献)

这个网站写的比较清楚https://www.statisticshowto.datasciencecentral.com/false-discovery-rate/

 

方法1：JMP软件

将p值复制到软件中，找到False Discovery Rate PValue插件

按照以下设置，点OK

校正后的p值就显示在原数据表中了

方法2 R

p.adjust(p,"BH")，其中p是原始p值，BH是选定的方法

可以发现这和JMP给出的方法是一致的

如果是多因素的分析结果，p值保存在矩阵或者数据框中，可以用一下循环一次解决（6个因素）

d<-p.adjust(p[,1],"BH") #先给d赋值
i=2     #已经有了第一个因素的结果，所以从第二列开始循环
while(i<=6){
  d<-cbind(d,d<-p.adjust(p[,i],"BH"));i<-i+1    #利用cbind()实现列追加
}
write.csv(d,"1.fdr.csv")   #输出结果