Storey FDR矫正方法

简介
零假设(null hypothesis)的P-values为均匀分布
真实多次测试的p-values将在0处出现峰值
从一个p-value阈值获得FDR
p-values映射到q-values
R语言实现

简介

John Storey 创建了一种将系列p值转换为q值的方法，不同之处在于p值用于衡量单个测试为空假设模型的累积概率，而q值用于衡量False Discovery。通过接受给定的测试以及每个测试具有较小p值。

q值需要解决的主要问题是多重假设检验问题。问题在于，长时间执行假设检验的标准p值在同时执行多个检验时并不可靠。这是个常见的生物学问题，例如当我们要问基因表达在整个基因组中是否显著（比如在酵母中检测了6,000多次，而在人类中检测了20,000多次）。所以在单次检验中，你可能因为p小于0.05而拒绝零假设（意思就是零假设产生的观测值或者极值只有小于5%的可能性）。但是在多重检验中，有可能出现任何假设错误的几率将会随着测试次数增加而极具增加，所以将p值的阈值设置为0.05将会导致错误的测试超过5%。

假设我们真正需要的是FDR值，FDR定义为通过阈值但是为假的测试比例。初次逼近时，似乎通过乘以测试次数来矫正p值可以比较好的近似所需FDR，而且实际上这正是Bonferroni矫正所作的。但是不幸的是，这个做法太苛刻而导致丢失太多好的测试结果。此后更高级的FDR矫正被开发出来，而Storey矫正被认为是最好的方法以取得真实测试和减少错误测试之间的平衡。Storey给出了很多很好的理论，到那时我们只用图形的方法对他的原理进行补充阐述，这样可是使大家直观理解并完成q值的过程。有兴趣可以参阅Storey 原文（Storey，2002）。

零假设(null hypothesis)的P-values为均匀分布

想象一下进行多个测试时，我们确认每次测试零假设为true。所有p值的直方图如上图所示。对于真实的数据，我们将会期望出现小的波动，上图没有显示，但是趋势仍然时均匀分布的。如果可以的话，推荐查看一下p值的直方图：收集一个希望零假设为真的数据集，如果p值不是均匀分布，那么就要查看有没有正确计算p值，或者检查零假设是否有误。p值正确对分析过程至关重要。

真实多次测试的p-values将在0处出现峰值

这里的假设是，在零检验为false的多次测试中，最小的p值会被富集。重要的是，p值的分布在远离p=0处应该保持均匀，当然仍然会有一些小波动，因为大部分的测试仍然为符合零检验。否则与之前同样需要警惕，要么p值计算有误，要么零假设有问题。

Storey FDR过程的关键在于：估计在p=0附近有多少false预测。想象p值的分布其实是两个基础分布的混合，一个是零假设为true，则p值分布均匀；另一个是零假设为false，则在p=0附近存在峰值。

Storesy FDR过程需要 π0\pi0π0，定义为零检验为真在所有检验中的比例。在上图中，大约为零检验分布的高度，很容易估计该值。估计时越远离p=1，那么估计π0\pi0π0所需要的数据越多，估计的方差则低，但是可能会冒零假设为false的风险，也就是获得的值会大于真实的π0\pi0π0。

如果确实觉得无法可靠的估计π0\pi0π0，那么就直接将其设置为1.这样将会减少FDR的功效，因为实际FDR将会小于估计，所以将会丢失些真实预测。但是通常π0\pi0π0非常接近于1，如果大多数检验满足零假设的话。实际上，设置为1是将Storey的过程简化为BH方法。

从一个p-value阈值获得FDR

如上图所示，假设t为p值的阈值，也就是说每个p<t的检验将被通过。FDR包含两个部分。假设总面积被归一化为1，分母时p<t的面积；或者是p<t的检验个数与总检验个数的比值。分子是p<t false tests的面积。为了重申我们前面的所说的false test总面积为π0\pi0π0，并考虑p<t面积的分数，因为是均匀的，所以就是t∗π0t*\pi0t∗π0。最后的估计false的公式显示在上图中。
Arandom(t)=tπ0Arandom(t)=t\pi0Arandom(t)=tπ0 为p<t时false的面积
KaTeX parse error: Expected 'EOF', got '#' at position 9: Aall(t)=#̲{p<t}/#{p} 为p<t时所有面积
KaTeX parse error: Expected 'EOF', got '#' at position 31: …t)/All(t)=t\pi0#̲{p}/#{p<t} 即为FDR

p-values映射到q-values

所以如果先选好了t再获得FDR就非常方便了。但是我们通常都是先想要FDR，然后再找到我们想要的t。甚至，我们是不是可以不设定FDR的情况下分析数据呢？Storey提供的定理表明我们可以很方便的完成上述任务。

首先，我们生成从 t 到 FDR(t) 的映射。本质上就是我们从t = 0步进到 t = 1，并且随之存储FDR值。如果仔细编写代码，可以非常有效地完成此操作。我们可以将 FDR§ 视为p的q值，但我们可以做得更好。通常 FDR§ 随p的增加而增加（如上图所示），但是实际数据波动很大时，情况并非总是如此。在这种情况下，t 阈值的 FDR(t) 可能小于我们查找的p！将FDR（t）用作p的q值是有意义的，因为我们将同时获得更多的预测和更低的FDR！

所以，q§= p<t 时 t 的最小值；FDR(t) 时 p值的最终定义。这样，q值时随p单调变化的。为了计算，一旦 FDR(t) 映射被计算，然后就反向沿着 t 从1到0的方向，以查找最小值。现在，如果你需要一个0.05的FDR，那么你就可以接受所有 q<0.05；而如果你需要一个0.1的FDR，那么就使用 q<0.1, 如此即可！

R语言实现

BiocManager::install("qvalue")
library(qvalue)
qvalue(p, fdr.level = NULL, pfdr = FALSE, lfdr.out = TRUE, pi0 = NULL, ...)

Arguments	—
p	A vector of p-values (only necessary input).
fdr.level	A level at which to control the FDR. Must be in (0,1]. Optional; if this is selected, a vector of TRUE and FALSE is returned that specifies whether each q-value is less than fdr.level or not.
pfdr	An indicator of whether it is desired to make the estimate more robust for small p-values and a direct finite sample estimate of pFDR – optional.
lfdr.out	If TRUE then local false discovery rates are returned. Default is TRUE.
pi0	–It is recommended to not input an estimate of pi0. Experienced users can use their own methodology to estimate the proportion of true nulls or set it equal to 1 for the BH procedure.
…	Additional arguments passed to pi0est and lfdr.

*** http://viiia.org/fdrFigs/?l=en-us