Применить коррекцию FDR к большому количеству исходных переменных

Question

Я использую R (не 4-версию, пока ахах). Мне посоветовали использовать коррекцию FDR на моих линейных моделях. У меня> 200 участников, 140 непрерывных переменных результата, и каждая переменная результата тестируется на тех же 4 предикторах. Итак, все модели: Y ~ x1 + x2 + x3 + x4, для всех 140 переменных, где x1 - интересующий меня предиктор, а остальные (x2, x3, x4) я просто использую для управления их влияние на Y. Как я могу применить FDR? Для чего я должен исправить? Должен ли я исправить все 140 переменных результата? Нужно ли исправлять только 4 предиктора? Если бы вы могли объяснить процесс и как решить, что нужно исправить в fdr, это было бы действительно хорошо, так как я изо всех сил пытаюсь понять это. Большое спасибо за помощь, Best

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 

Теперь мы очистим некоторые термины, сгруппировав по ответ, и мы можем применить FDR, используя p.adjust, «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 <0,05: </p>

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 противостоит тренду, описанному выше, и мы можем визуализировать это, построив график распределения значений:

library(ggplot2)
adjusted %>%
ggplot(aes(x=p.value)) + geom_histogram() +
facet_wrap(~term) + theme_bw()

Для x2, x3 и x4, мы смоделировали их под нулем, не влияя ни на какие ответы, и вы можете видеть, что значение p следует за равномерным распределением.

Если мы просто используем отсечение 0,05, мы получим все 7 ложных срабатываний в других предикторах x1-x4, в то время как некоторые из них в 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 попадание. Вы также можете посмотреть это видео , в котором более подробно объясняется, как это работает