11 ноября 2013

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



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


Комментариев нет :

Отправить комментарий