Bonferroni与Benjamin Hochberg校正无效问题求助（附R代码）

问题背景 #

现有1020个样本、156个特征，通过Wilcoxon检验得到p值，需对总计1020次检验结果做Bonferroni和Benjamin Hochberg校正，但当前代码输出的校正后p值与原始值完全一致。

数据示例 #

原始p值样本结构：

structure(list(p_values = c(0.00551261839474566, 0.00909340979590469, 
0.42610555368556, 0.711610700326496, 0.00439218856215691, 0.859681237958105)), 
class = "data.frame", row.names = c("a","b","c","d","e","f"))

问题代码 #

当前使用的R代码：

rf0->d. # 基础数据集，含1020个样本、156个待检验特征
out <- lapply(3:158, function(x) pairwise.wilcox.test(d[[x]], d$LesionResponse,p.adjust.method="bonferroni"))
names(out) <- names(d)[3:158]
pvalue<-sapply(out, function(x) {
    p <- x$p.value
    n <- outer(rownames(p), colnames(p), paste, sep='v')
    p <- as.vector(p)
    names(p) <- n
    p
})
pvalue<-as.data.frame(pvalue)

问题原因 #

pairwise.wilcox.test中的p.adjust.method参数仅针对单个特征内部的两两比较做校正，而非对所有156个特征的全部1020次检验做全局校正。因此该参数无法实现你需要的全局多重检验校正，导致输出p值与原始值一致。

解决方法 #

需先提取所有未校正的原始p值，再用p.adjust函数对所有p值做全局校正：

修正代码 #

# 1. 提取所有特征的未校正原始p值
out <- lapply(3:158, function(x) {
  pairwise.wilcox.test(d[[x]], d$LesionResponse, p.adjust.method = "none")
})
names(out) <- names(d)[3:158]

# 2. 将原始p值整理为数据框
p_raw <- sapply(out, function(x) {
  p <- x$p.value
  n <- outer(rownames(p), colnames(p), paste, sep='v')
  p_vec <- as.vector(p)
  names(p_vec) <- n
  p_vec
})
p_raw_df <- as.data.frame(p_raw)

# 3. 对所有p值做全局校正
all_raw_p <- unlist(p_raw_df)

# Bonferroni校正
p_bonferroni <- p.adjust(all_raw_p, method = "bonferroni")
# Benjamin Hochberg校正
p_bh <- p.adjust(all_raw_p, method = "BH")

# 4. 将校正后p值还原为原数据框结构
p_bonferroni_df <- as.data.frame(
  matrix(p_bonferroni, nrow = nrow(p_raw_df), 
         dimnames = dimnames(p_raw_df))
)
p_bh_df <- as.data.frame(
  matrix(p_bh, nrow = nrow(p_raw_df), 
         dimnames = dimnames(p_raw_df))
)

说明 #

• 第一步设置p.adjust.method="none"，确保拿到的是未经过任何校正的原始p值；
• 第二步将所有特征的原始p值整理为统一的数据框；
• 第三步用p.adjust对所有p值做全局校正，这一步才是针对1020次检验的多重比较校正；
• 最后将校正后的p值还原为与原始数据框一致的结构，方便后续分析。

内容的提问来源于stack exchange，提问作者NDe