对大量结果变量应用FDR校正

Question

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

Answer 1

所以你需要在预测值和结果之间控制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个点击。您还可以查看this video，它更详细地解释了上面的工作原理