В предыдущем сообщении был описан метод Беньямини-Хохберга, широко используемый для контроля над ожидаемой долей ложных отклонений при проверке большого числа статистических гипотез. Одно из условий применимости этого метода заключается в том, что все проверяемые гипотезы должны быть независимы. На практике это условие выполняется редко, поскольку в большинстве случаев гипотезы проверяются на одном и том же наборе данных. Понимая важность этого ограничения, И. Беньямини и Д. Йекутили (Benjamini and Yekutieli 2001) предложили усовершенствованный метод, учитывающий наличие корреляции между проверяемыми гипотезами (подробное описание метода и соответствующие доказательства см. в указанной оригинальной статье).
Процедура Беньямини-Йекутили
Процедура Беньямини-Йекутили очень похожа на процедуру Беньямини-Хохберга. Основное отличие заключается во введении поправочной константы \(c(m) = \sum_{i=1}^{m}\frac{1}{i}\):
- Исходные Р-значения упорядочивают по возрастанию: \(p_{(1)} \leq p_{(2)} \leq \dots \leq p_{(m)}\). Пусть \(H_{(i)}\) обозначает нулевую гипотезу, которой соответствует i-тое значение в этом упорядоченном ряду, т.е. \(p_{(i)}\).
- Находят максимальное значение \(k\) среди всех индексов \(i = 1, 2, \dots, m\), для которого выполняется неравенство \(p_{(i)} \leq \frac{i}{m}\times\frac{q}{c(m)}\)
- Отклоняют все гипотезы \(H_{(i)}\) с индексами \(i = 1, 2, \dots, k\)
В качестве примера рассмотрим следующий ряд из 15 упорядоченных по возрастанию Р-значений (из статьи Benjamini and Hochberg 1995):
0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000
При контроле над ожидаемой долей ложных отклонений на уровне 5% мы сравниваем каждое значение \(p_{(i)}\) с \( \frac{i}{15} \times \frac{0.05}{c(m)}\), начиная с самого высокого, т.е. \(p_{(15)}\). В данном случае константа \(c(m) = \sum_{i=1}^{15}\frac{1}{i} = 3.31823\) В итоге мы увидим, что первое Р-значение, соответствующее указанному ограничению (5%), - это \(p_{(3)}\):
\[p_{(3)} = 0.0019 \leq (3/15)\times (0.05/3.31823) = 0.003014 \]
Таким образом, следует отклонить три гипотезы, которым соответствуют первые три Р-значения в приведенном выше ряду (поскольку все эти значения не превышают 0.003014).
Реализация в R
Подобно тому, как это было с другими методами множественных проверок гипотез (см. здесь и здесь), контроль над ожидаемой долей ложных отклонений по методу Беньямини-Йекутили можно реализовать в R при помощи функции p.adjust(), подав на нее вектор с исходными Р-значениями и присвоив аргументу method значение "BY" (от "Benjamini-Yekutieli"). BY-контроль над FDR при помощи этой функции выполняется несколько отличным от описанной выше процедуры образом. В частности, вместо нахождения максимального индекса \(k\), выполняется следующее преобразование исходных Р-значений:
\[q_{(i)} = \frac{p_{(i)}\times m \times с(m)}{i}\]
Если то или иное преобразованное Р-значение превышает 1, его просто приравнивают к 1.
Например, для первых двух Р-значений из приведенного выше примера мы получили бы
Например, для первых двух Р-значений из приведенного выше примера мы получили бы
- \((0.0001 \times 15 \times 3.31823)/1 = 0.004977\)
- \((0.0004 \times 15 \times 3.31823/2) = 0.009955\)
Воспользуемся функцией p.adjust() для преобразования всех Р-значений из рассмотренного примера:
Интерпретация скорректированных Р-значений выполняется тем же образом, что и в случае с методом Беньямини-Хохберга.
Помимо функции p.adjust(), корректировку Р-значений для контроля над ожидаемой долей ложных отклонений можно также выполнять при помощи других функций, входящих в состав соответствующих пакетов. Список этих пакетов можно найти, например, здесь.
Этим сообщением я завершаю рассмотрение процедур множественных сравнений, реализованных в функции p.adjust(). Помимо рассмотренных, эта функция предлагает также методы Хохберга (не путать с методом Беньямини-Хохберга) и Хоммеля (Hommel) (подробнее см. stats.stackexchange.com). Следующие несколько сообщений будут посвящены функциональным возможностям multcomp - очень мощного R-пакета, реализующего методы множественных сравнений на основе теории линейных моделей (см. введение здесь и здесь).
pvals <- c(0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000) p.adjust(pvals, "BY") [1] 0.004977343 0.009954687 0.031523175 0.118211908 0.200089208 0.211892623 [7] 0.211892623 0.214025770 0.253844518 1.000000000 1.000000000 1.000000000 [13] 1.000000000 1.000000000 1.000000000
Интерпретация скорректированных Р-значений выполняется тем же образом, что и в случае с методом Беньямини-Хохберга.
Помимо функции p.adjust(), корректировку Р-значений для контроля над ожидаемой долей ложных отклонений можно также выполнять при помощи других функций, входящих в состав соответствующих пакетов. Список этих пакетов можно найти, например, здесь.
Этим сообщением я завершаю рассмотрение процедур множественных сравнений, реализованных в функции p.adjust(). Помимо рассмотренных, эта функция предлагает также методы Хохберга (не путать с методом Беньямини-Хохберга) и Хоммеля (Hommel) (подробнее см. stats.stackexchange.com). Следующие несколько сообщений будут посвящены функциональным возможностям multcomp - очень мощного R-пакета, реализующего методы множественных сравнений на основе теории линейных моделей (см. введение здесь и здесь).
Отправить комментарий