高通量数据的多重检验问题

Bonferroni 方法

思路很简单，就是控制显著性，例如单次检验假阳性比率\(\alpha\)控制在0.05，那么n次检验假阳性比率控制为\(\frac{\alpha}{n}\)。这样实际是对整体采用了个体控制的控制思路：

\[ P(至少一个显著)=1-P(无显著差异) = 1-(1-\alpha/n)^n \]

我们来看下\(\alpha = 0.05\)随比较数增加的效果：

n <- c(1:10 %o% 10^(1:2))
p0 <- 1-(1-0.05)^n
p <- 1-(1-0.05/n)^n
# 不进行控制
plot(p0~n,ylim = c(0,1))
# Bonferroni方法控制
points(p~n,pch=19)

plot of chunk Bonferroni

其实，这样的控制得到的整体错误率是略低于0.05的，并且数目越大，整体错误率越低。这个方法十分保守，有可能什么差异你都看不到，因为都变成假阴性了。在实际应用中一般不调节p值的假阳性比率而直接调节p值，取原始p值跟整体检验数目的乘积与1的最小值作为调节p值，还可以用0.05或0.01进行判断，不过这时候控制的整体而不是单一检验了。

当然这只是最原始的Bonferroni方法，后来Holm改进了这种一步法为逐步法，此时我们需要首先对原始p值进行排序，然后每个原始p值乘上其排序作为调节p值。例如三次多重检验的p值分别是0.01、0.03与0.06，其调节后的p值为0.03，0.06，0.06。如果我们控制整体假阳性比率低于0.05，那么调解后只有第一个检验可以拒绝空假设。值得注意的是Holm的改进是全面优于原始方法的，也就是说当你一定要去用Bonferroni方法控制整体错误率，优先选Holm的改进版。

Sidak 方法

上面那种方法其实有点非参的意思，其实数学上我们是可以精确的把假阳性比率控制在某个数值的：

\[ P(至少一个显著)=1-P(无显著差异) = 1-(1-\alpha')^n = 0.05 \]

求解可得到\(\alpha' = 1-0.95^{\frac{1}{n}}\)，此时我们就可以比较精确的控制整体错误率了，但是，这个方法有个前提就是各个检验必须是独立的，这在生物学实验里几乎不可能，所以这个方法的应用远没有Bonferroni方法广。

Benjamini-Hochberg方法

这个方法跟Holm方法很像，也是先排序，但之后p值并不是简单的乘排序，而是乘检验总数后除排序：

\[ p_i \leq \frac{i}{m} \alpha \]

举例来说就是假设三次多重检验的p值分别是0.01、0.03与0.06，其调节后的p值为0.03，0.045，0.06。那么为什么说这种方法控制的是错误发现率呢？我们来看下\(\alpha\)是如何得到的：p值乘总数m得到的是在该p值下理论发现数，而除以其排序实际是该p值下实际发现数，理论发现数基于在这里的分布是均匀分布，也就是空假设的分布，这两个的比值自然就是错误发现率。下面我用仿真实验来说明一下：

pbh <- p.adjust(pvalue,method = 'BH')
ph <- p.adjust(pvalue,method = 'holm')
plot(pbh~pvalue)
points(ph~pvalue,col='red')

plot of chunk BH

从上面图我们可以看出，如果控制整体错误率（红色），那么p值很容易就到1了，过于严格。而如果用BH方法控制错误发现率，那么原始p值越大，调节后的错误发现率也逐渐递增，这就符合了区分真实差异与随机差异就要假设真实差异更可能出现更小的p值这个现象。当然至于这个方法的推演细节，可以去读原始论文。值得注意的是这个错误发现率求的是有差异存在的情况，不然零发现就出现除数为零了。

Storey方法（q值）

如果说BH方法还算是调节了p值，那么Storey提出的方法则直接去估计了错误发现率本身。刚才介绍BH算法时我提到总数m与p值的乘积是基于这里的分布是均匀分布，但实际上按照错误发现率的定义，这里应该出现的是空假设总数。直接使用所有检验数会造成一个问题，那就是对错误发现率的高估，为了保证功效，这里应该去估计空假设的总体比例。这里我们去观察混合分布会发现在p值较大的时候基本可以认为这里分布的都是空假设的p值，那么我们可以用：

\[ \hat\pi_0 = \frac{\#\{p_i>\lambda\}}{(1-\lambda)m} \]

估计这个比例\(\hat\pi_0\)，其中参数\(\lambda\)的跟\(\hat\pi_0\)的关系可以用一个三阶方程拟合，然后计算出整体假阳性比例。有了这个比例，我们再去按照BH方法计算p值，然后两个相乘就会得到q值，而q值的理论含义就是在某一概率上低于这个概率所有数里假阳性的比重。打个比方，我测到某个指标的q值是0.05，这意味着q值低于这个数所有检验中我有0.05的可能性得到的是假阳性。。但我们会发现当空假设比重较高时BH结果跟q值很接近，而比重很低的话q值会变得更小，功效会提高，基本也符合我们对错误发现率的预期。

library(qvalue)
q <- qvalue(pvalue)
# Q值
plot(q$qvalues~pvalue,col='blue')

plot of chunk qvalues

如上图所示，q值增大后会最终逼近到0.5，而我们的模拟中空假设的比例就设定就是50%。我们重新模拟一个空假设比例5%的实验：

set.seed(42)
pvalue <- NULL
for (i in 1:500){
  a <- rnorm(10,1)
  b <- a+1
  c <- t.test(a,b)
  pvalue[i] <- c$p.value
}
for (i in 1:9500){
  a <- rnorm(10,1)
  b <- rnorm(10,1)
  c <- t.test(a,b)
  pvalue[i+500] <- c$p.value
}
pbh <- p.adjust(pvalue,method = 'BH')
ph <- p.adjust(pvalue,method = 'holm')
q <- qvalue(pvalue)
plot(pbh~pvalue)
# Holm 方法
points(ph~pvalue,col='red')
# Q值
points(q$qvalues~pvalue,col='blue')

plot of chunk qbh

此时我们可以看到两者结果较为接近，q值理论上更完备，功效也更强，但算法上对\(\hat\pi_0\)的估计并不稳定，特别是比例靠近1的时候，所以BH方法可能还是更容易让人接受的保守错误发现率控制。详细的估计方法还得去啃Storey的论文。