我正在使用 R（还不是 4 版本啊）我被建议在我的线性模型上使用 FDR 校正。我有超过 200 名参与者，140 个连续结果变量，每个结果变量都在相同的 4 个预测变量上进行测试。所以所有模型都是：Y ~ x1 + x2 + x3 + x4，对于所有 140 个变量，其中 x1 是我感兴趣的预测变量，而其他 (x2,x3,x4) 我只是用来控制它们对 Y 的影响。我如何应用 FDR？我要纠正什么？我必须纠正所有 140 个结果变量吗？我只需要纠正 4 个预测变量吗？如果您可以解释该过程以及如何决定在 fdr 中纠正什么，那将非常好，因为我正在努力理解它。非常感谢您的帮助，最好的
问问题
193 次
1 回答
0
因此，您需要控制预测变量和结果之间的 140 次测试，并对每个预测变量进行 FDR。我们可以尝试一个示例，其中 x1 对响应 y 1 到 30 有影响，而对其他响应没有影响，而 x2,x3,x4 则没有，首先是数据：
set.seed(111)
X = matrix(runif(200*4),ncol=4)
colnames(X) = paste0("x",1:4)
Y = matrix(rnorm(140*200),ncol=140)
colnames(Y) = paste0("y",1:140)
Y[,1:30] = 1.5*X[,1]+Y[,1:30]
很好用broom，整理一下，我们可以拟合一个多响应线性模型，但是每个 Y 都是单独回归的，输出是这样的：
library(broom)
library(dplyr)
model = lm(Y ~ X)
tidy(model)
# A tibble: 700 x 6
response term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 y1 (Intercept) 0.288 0.268 1.07 0.285
2 y1 Xx1 1.22 0.248 4.93 0.00000178
3 y1 Xx2 -0.356 0.251 -1.42 0.158
现在我们清理一些术语，按响应分组，我们可以使用 p.adjust 应用 FDR，“BH”代表Benjamini-Hochberg：
adjusted = tidy(model) %>%
mutate(term=gsub("X","",term)) %>%
filter(term!="(Intercept)") %>%
group_by(term) %>%
mutate(padj = p.adjust(p.value,"BH")) %>%
ungroup()
因此，在我们查看 FDR 结果之前，我们可以考虑这样的多重测试。如果预测变量对任何响应都没有影响，并且您进行了 140 次检验，那么您预计大约 0.05*140 = 7 次的检验会为您提供 0.05 的 p 值。我们可以检查每个预测变量，其中有多少 p < 0.05：
adjusted %>% group_by(term) %>% summarize(sig=sum(p.value<0.05))
# A tibble: 4 x 2
term sig
<chr> <int>
1 x1 36
2 x2 7
3 x3 6
4 x4 7
p 值分布如何？所以你可以在上面看到 x1 逆势而上，我们可以通过绘制 pvalue 分布来可视化这一点：
library(ggplot2)
adjusted %>%
ggplot(aes(x=p.value)) + geom_histogram() +
facet_wrap(~term) + theme_bw()
对于 x2、x3 和 x4，我们在 null 下模拟它们，对任何响应没有影响，您可以看到 p 值遵循均匀分布。
如果我们简单地使用 0.05 的截止值，我们将在其他预测变量 x1-x4 中得到所有 7 个误报，而其中一些在 x1 中是正确的。FDR 基本上纠正了这种预期的 p 值分布，我们可以检查其中有多少在 5% FDR 时显着：
adjusted %>% group_by(term) %>% summarize(sig=sum(padj<0.05))
# A tibble: 4 x 2
term sig
<chr> <int>
1 x1 31
2 x2 0
3 x3 0
4 x4 0
因此，我们不会再使用 x2、x3、x4 获得任何命中，而 x1，我们在 30 个真实效果下模拟得到 31 个命中。您还可以查看此视频，该视频更详细地解释了上述工作原理
于 2020-04-27T14:48:16.880 回答