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

如果进行m次假设检验， ：?=0， ：?≠0。可能出现的结果如下：

?0为H0为真的次数，R为拒绝H0的次数。 I型错误或假阳性错误：?=0时判断结果为?≠0，发生V次。 II型错误或假阴性错误：?≠0时判断结果为?=0，发生T次。

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

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

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

在统计分析时进行多次假设检验，多重检验校正可以降低假阳性结果的发生。

➢校正显著性水平?

控制FWER——Bonferroni 校正

控制FWER的典型方法是Bonferroni 校正 ，它拒绝了所有的假阳性结果发生的可能性。

基本思路

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

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

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

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

‍‍‍‍控制FDR——BH校正‍‍‍‍

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

基本思路

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

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

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

缺点：允许较Bonferroni校正更多的假阳性结果。

➢校正P值（Adjusted P-values）

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

控制 FWER

假设进行m次检验，P值为 ，…， ，校正的P值 =min ?× , 1。校正P值不大于1，?× 大于1时校正P值取值为1。

回顾Bonferroni校正， =?/?。

当 <?时，认为结果是阳性的，即把FWER控制在?水平。

控制FDR

假设进行?次检验，计算每次检验的P值，结果按由小到大进行排序 ，…， ，校正的P值 = × ，此时校正的P值又称为Q值。

回顾BH校正， =?× 。

当 <?时，认为结果是阳性的，即把 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相关（阳性TRUE）的结论，即有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)

Bonferroni校正（图左），校正P值?×P?不大于1，故图的上方可以看到一条水平为1的线。

BH校正（图右），校正P值是关于P值的递增函数，校正P值比实际P值本身稍大，但并没有特别大，因为需要得到较多阳性结果。

编辑：李雪纯 冯文清

校审：张健 罗鹏