Процедуры множественных проверок гипотез: метод Беньямини-Йекутили

В предыдущем сообщении был описан метод Беньямини-Хохберга, широко используемый для контроля над ожидаемой долей ложных отклонений при проверке большого числа статистических гипотез. Одно из условий применимости этого метода заключается в том, что все проверяемые гипотезы должны быть независимы. На практике это условие выполняется редко, поскольку в большинстве случаев гипотезы проверяются на одном и том же наборе данных. Понимая важность этого ограничения, И. Беньямини и Д. Йекутили (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() для преобразования всех Р-значений из рассмотренного примера:

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-пакета, реализующего методы множественных сравнений на основе теории линейных моделей (см. введение здесь и здесь).