FDR计算 - 孟德尔随机化教程

📖 什么是FDR？

FDR（False Discovery Rate）即错误发现率，是多重检验校正的常用方法。

控制错误拒绝零假设的比例
比Bonferroni校正更宽松，统计效力更高
适用于大规模假设检验（如GWAS、代谢物分析）

🔬 Benjamini-Hochberg方法

BH方法步骤：

将所有p值从小到大排序：p(1) ≤ p(2) ≤ ... ≤ p(m)
对于每个p(i)，计算校正后的阈值：q × i / m
找到最大的i，使得 p(i) ≤ q × i / m
该i之前的所有检验被认为是显著的

其中：q = 目标FDR阈值（通常0.05），m = 检验总数

💻 R代码实现

方法1：使用p.adjust函数

# 原始p值向量
p_values <- c(0.001, 0.005, 0.01, 0.03, 0.05, 0.1, 0.5, 0.8)

# 使用BH方法进行FDR校正
p_adjusted <- p.adjust(p_values, method = "BH")

# 查看结果
data.frame(
  p_original = p_values,
  p_adjusted = p_adjusted,
  significant = ifelse(p_adjusted < 0.05, "Yes", "No")
)

方法2：使用qvalue包（更精确）

# 安装qvalue包
if (!require("qvalue")) {
  install.packages("qvalue")
  library(qvalue)
}

# 使用qvalue进行FDR估计
qobj <- qvalue(p = p_values)

# 查看q值
qvalues(qobj)

# 识别显著结果
sig <- sum(qvalues(qobj) < 0.05)

# 查看结果摘要
summary(qobj)

方法3：手动实现BH方法

# 手动实现Benjamini-Hochberg方法
bh_correction <- function(p_values, q = 0.05) {
  m <- length(p_values)
  ranks <- rank(p_values, ties.method = "min")
  
  # 计算每个p值的临界值
  critical_values <- q * ranks / m
  
  # 找到最大的满足p(i) <= q*i/m的i
  significant <- p_values <= critical_values
  
  # 调整后的p值（最小化FDR）
  adjusted <- pmin(cummin(rev(p_values * m / rev(ranks))), 1)
  
  return(list(
    adjusted = adjusted,
    significant = significant
  ))
}

# 使用
result <- bh_correction(p_values, q = 0.05)
data.frame(
  p_original = p_values,
  p_adjusted = result$adjusted,
  significant = result$significant
)

📊 不同校正方法比较

方法	控制	特点
Bonferroni	FWER	最保守，易产生假阴性
Benjamini-Hochberg	FDR	平衡敏感性和特异性
q-value	FDR	允许估计显著比例

⚠️ 注意事项

• FDR适用于多个独立检验的情况
• MR分析中通常对多个SNP或多个结局进行校正时使用
• 0.05的FDR阈值意味着预期5%的显著结果是假阳性
• 如果需要更严格控制错误，可使用更小的FDR阈值（如0.01）

📚 相关教程

📊

🔗