0

我正在使用 R（还不是 4 版本啊）我被建议在我的线性模型上使用 FDR 校正。我有超过 200 名参与者，140 个连续结果变量，每个结果变量都在相同的 4 个预测变量上进行测试。所以所有模型都是：Y ~ x1 + x2 + x3 + x4，对于所有 140 个变量，其中 x1 是我感兴趣的预测变量，而其他 (x2,x3,x4) 我只是用来控制它们对 Y 的影响。我如何应用 FDR？我要纠正什么？我必须纠正所有 140 个结果变量吗？我只需要纠正 4 个预测变量吗？如果您可以解释该过程以及如何决定在 fdr 中纠正什么，那将非常好，因为我正在努力理解它。非常感谢您的帮助，最好的

4

1 回答 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 回答