数据科学23 | 统计推断-多重检验

当我们进行数据分析时，有时候需要反复进行假设检验，使用多重检验校正可以避免假阳性的发生，主要包括误差测量和校正。

错误类型

假设检验H0：?=0，H1：?≠0。可能出现的结果如下：

I型错误或假阳性错误(V为发生的次数)：参数?=0时主张?≠0。II型错误或假阴性错误(T为发生的次数)：参数?≠0时主张?=0。

• 假阳性率(False positive rate, FPR)：误报率，将“阴性”的错误结果(?=0)称为“阳性”的概率，E[V/?0]。
• FWER (Family-Wise Error Rate) ：出现至少一次假阳性错误的概率，Pr（V≥1）。
• FDR (False Discovery Rate) ：所有主张“阳性”的次数中，错误主张（假阳性）所占的比例E[V/R]。

目的：控制假阳性率FPR

如果正确计算了P值，所有P值小于?时被称为阳性，假阳性率FPR将为?。

假设进行10000次假设检验?=0，其中所有P<0.05的认为阳性，?=0.05，预期的误报数为10000×0.05=500。所以如果进行了10000次假设检验并获得500个阳性结果，其中很有可能有大部分的结果是假阳性。

用多重检验来进行校正，减低假阳性结果出现的次数。

校正?值（Adjusted P-values）

• 控制FWER (Family-Wise Error Rate)

控制 FWER 比较典型的方法就是 Bonferroni 校正 ，它拒绝了所有的假阳性结果发生的可能性，通过对p值的阈值进行校正来实现消除假阳性结果。

基本思路

假设进行?次检验，希望控制FWER使Pr（V≥1）<?，将每次检验的 I 型错误率控制在?/?之内。经过Bonferroni校正，?fwer=?/?，每次检验的P值小于?fwer时认为阳性。

进行10000检验，?fwer=0.05/10000=0.000005，当P<0.000005时认为是阳性的，假阳性结果出现的次数为10000×0.000005=0.5，不到1次。

优点：计算简单，犯I型错误的概率很小；

缺点：过于保守，犯II型错误的概率（假阴性率）增大，很多阳性结果无法被检测出来。

• 控制FDR (False Discovery Rate)

FDR（False Discovery Rate）校正没有 FWER 法那么保守严苛 ，但又不像 t 检验那样有很高的假阳性率。常用的方法为BH（Benjaminiand Hochberg）校正。

基本思路

假设进行?次检验，希望控制FDR使E[V/R]<?。计算每次检验的P值，结果按由小到大进行排序P(1)，…，P(?)，找到第?个P值，当P(?)≤?×?/?时，认为是阳性的，此时第1到第?个P值对应的检验都认为是阳性的。

进行10000次检验，控制E[V/R]<0.05，如果得到500个阳性结果，其中假阳性结果小于500×0.05=25个。

优点：计算简单，不会过于保守；

缺点：允许更多的假阳性结果

校正P值（Adjusted P-values）

除了校正?水平，另一种多重检验校正方式是校正P值。

• 控制 FWER ：

假设进行m次检验，P值为P1，…，P?，校正的P值=max ?×P?, 1。校正P值不大于1，?×P?大于1时校正P值取值为1。回顾Bonferroni校正，?fwer=?/?。当<?时，认为是阳性的，即把 FWER控制在?水平。

• 控制FDR：

假设进行?次检验，计算每次检验的P值，结果按由小到大进行排序P(1)，…，P(?)，校正的P值=P(?)×(?/?)，此时校正的P值又称为Q值。回顾BH校正，?fdr=?×?/?。当Q?<?时，认为是阳性的，即把 FDR控制在?水平。

例：模拟生成1000个没有阳性结果的数据集

 set.seed(1010093)
 pValues <- rep(NA, 1000)
 for (i in 1:1000) {
   y <- rnorm(20)
   x <- rnorm(20)
   pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
 }

生成1000个数据集，每个数据集中生成互不相关的正态随机数y和x。建立变量x和y之间的线性相关模型，并得到它们的相关系数矩阵，矩阵的第二行第四列的元素即为P值。

没有校正，查看小于0.05的P值的数量：

 sum(pValues < 0.05)
 [1] 51

实际上所有数据集中变量x和y是不相关的，但仍有51个数据集得到x与y相关的结论，即有51个假阳性结果。

Bonferroni校正控制FWER：

 sum(p.adjust(pValues, method = "bonferroni") < 0.05)
 [1] 0

BH校正控制FDR：

 sum(p.adjust(pValues, method = "BH") < 0.05)
 [1] 0

经过两种方法校正之后，减少了假阳性结果。

例：模拟生成1000个有50%阳性结果的数据集

 set.seed(1010093)
 pValues <- rep(NA, 1000)
 for (i in 1:1000) {
   x <- rnorm(20)
   # First 500 beta=0, last 500 beta=2
   if(i<=500){
     y <- rnorm(20) }
   else{
     y<-rnorm(20,mean=2*x)
  }
   pValues[i] <- summary(lm(y ~ x))$coeff[2, 4]
 }
 trueStatus <- rep(c("zero", "not zero"), each = 500)

生成1000个数据集，前500个仍旧包含不相关的随机变量x和y，后500个生成y的均值为x的2倍，y和x之间存在关系。前500个?等于零，而后500个?等于2，定义trueStatus声明实际?前500个为"zero"，后500个"not zero"。

没有校正，查看小于0.05的P值的数量：

 table(pValues < 0.05, trueStatus)
        trueStatus
         not zero zero
   FALSE        0  476
   TRUE       500   24

500个阳性结果全部被检测到；但实际x与y不相关时，有24个数据集得到x与y相关的结论，即有24个假阳性结果。

Bonferroni校正控制FWER：

 table(p.adjust(pValues, method = "bonferroni") < 0.05, trueStatus)
        trueStatus
         not zero zero
   FALSE       23  500
   TRUE       477    0

经过Bonferroni校正之后，没有假阳性结果，但出现了23个假阴性结果。

BH校正控制FDR：

 table(p.adjust(pValues, method = "BH") < 0.05, trueStatus)
        trueStatus
         not zero zero
   FALSE        0  487
   TRUE       500   13

经过BH校正之后，500阳性结果全部被检测到，假阳性结果较校正前减少到13个。

 par(mfrow = c(1, 2))
 plot(pValues, p.adjust(pValues, method = "bonferroni"), pch = 19)
 plot(pValues, p.adjust(pValues, method = "BH"), pch = 19)

图1.两种校正方法的校正P值

对于Bonferroni校正，校正后的P值?×P?不大于1，故图的上方可以看到一条水平为1的线；对于BH校正，校正P值是关于P值的递增函数，校正后的P值比实际P值本身稍大，但并没有特别大，因为需要得到许多阳性结果。

编辑：冯文清 李雪纯

校审：张健 罗鹏