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

Рассмотренные ранее процедуры множественных проверок гипотез обеспечивают контроль над групповой вероятностью ошибки первого рода (в частности, методы Бонферрони, Холма и Тьюки). Такой контроль означает, что вероятность совершить хотя бы одну ошибку первого рода удерживается на уровне \(\leq \alpha \), где \(\alpha\) - принятый в исследовании уровень значимости (например, 0.05). Хотя методы Холма и Тьюки обладают более высокой статистической мощностью, чем поправка Бонферрони, при проверке очень большого числа гипотез их мощность может оказаться недостаточной. Под "недостаточной мощностью" имеется в виду сохранение многих нулевых гипотез, которые потенциально могут представлять исследовательский интерес и которые, соответственно, следовало бы отклонить. Сегодня проверка действительно большого числа гипотез (десятков тысяч и даже миллионов) стала рутинной операцией в самых разных областях, таких как генетика (анализ данных, получаемых при помощи технологии микрочипов), протеомика (данные масс-спектрометрии), нейробиология (анализ изображений мозга), экология, астрофизика, и др. Недостаточная мощность традиционных процедур множественной проверки гипотез привела к разработке новых методов, одному из которых - методу Беньямини-Хохберга - посвящено это сообщение.

Контроль над ожидаемой долей ложных отклонений

Результаты проверки большого числа статистических гипотез принято сводить в виде следующей таблицы (содержание этой таблицы подробно обсуждалось здесь):

	Число принятых гипотез	Число отклоненных гипотез	Всего
Число верных гипотез	\(U\)	\(V\)	\(m_0\)
Число неверных гипотез	\(T\)	\(S\)	\(m - m_0\)
Всего	\(W\)	\(R\)	\(m\)

Как было отмечено ранее, групповая вероятность ошибки первого рода представляет собой вероятность отклонить как минимум одну верную нулевую гипотезу при проверке \(m\) гипотез: \(FWER = P(V \geq 1) \). Контроль над этой ошибкой на уровне значимости \(\alpha\) означает, что должно выполняться неравенство \(FWER \leq \alpha \). Проблема с контролем над FWER заключается в его консерватизме, который проявляется в следующем: при больших значениях \(m\) многие нулевые гипотезы, которые должны были бы быть отклонены, сохраняются. В итоге исследователь может пропустить интересные "открытия" (англ. discoveries), достойные более подробного изучения. Например, в ходе сравнения уровней экспрессии генов у больных и здоровых испытуемых, результаты для некоторых потенциально важных генов могли бы оказаться ложно-отрицательными.

В 1995 г. израильские исследователи Иоав Беньямини (Yoav Benjamini) и Йозеф Хохберг  (Yosef Hochberg) опубликовали статью, в которой был предложен принципиально иной подход к проблеме множественных проверок статистических гипотез (эта работа входит в список 25 наиболее цитируемых статьей по статике). Суть предложенного подхода заключается в том, что вместо контроля над групповой вероятностью ошибки первого рода выполняется контроль над ожидаемой долей ложных отклонений (англ. false discovery rate, FDR) среди всех отклоненных гипотез. Используя обозначения из приведенной выше таблицы, математически эта ожидаемая доля может быть записана следующим образом:

\[ FDR = \mathbf{E}\left ( \frac{V}{R} \right ) \]

В отличие от уровня значимости \(\alpha\), каких-либо "общепринятых" значений FDR не существует. Многие исследователи по аналогии контролируют FDR на уровне 5%. В генетических исследованиях часто встречается также уровень 10%. Интерпретация порогового значения FDR очень проста: например, если в ходе анализа данных отклонено 1000 гипотез, то при \(q = 0.10\) ожидаемая доля ложно отклоненных гипотез не превысит 100.

Процедура Беньямини-Хохберга

В оригинальной статье Беньямини и Хохберга (Benjamini and Hochberg 1995) описание процедуры контроля над FDR выглядит так:

• Исходные Р-значения упорядочивают по возрастанию: \(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}q\)
• Отклоняют все гипотезы \(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

Если бы мы осуществили контроль над групповой вероятностью ошибки, применив, например, поправку Бонферрони для уровня значимости 0.05 (т.е. 0.05/15 = 0.0033), то отклоненными оказались бы три гипотезы, которым соответствуют первые три Р-значения.

При контроле над ожидаемой долей ложных отклонений на уровне 5% мы сравниваем каждое значение \(p_{(i)}\) с \( 0.05i/15\), начиная с самого высокого - \(p_{(15)}\). В итоге мы увидим, что первое Р-значение, соответствующее указанному ограничению (5%) - это \(p_{(4)}\):

\[p_{(4)} = 0.0095 \leq (4/15)0.05 = 0.013 \]

Теперь мы отклоняем четыре гипотезы, которым соответствуют первые четыре Р-значения в приведенном выше ряду (поскольку все эти значения не превышают 0.013).


Реализация в R

Подобно тому, как это было с поправкой Бонферрони и методом Холма, контроль над ожидаемой долей ложных отклонений можно реализовать в R при помощи функции p.adjust(), подав на нее вектор с исходными Р-значениями и присвоив аргументу method значение "BH" (что значит "Benjamini-Hochberg") или значение-синоним "fdr" ("false discovery rate"). Опять-таки, подобно поправке Бонферрони и методу Холма, контроль над FDR при помощи этой функции выполняется несколько отличным от оригинально описанной процедуры образом. В частности, вместо нахождения максимального индекса \(k\) (см. выше), исходные Р-значения корректируются следующим образом:

\[q_{(i)} = \frac{p_{(i)}m}{i}\]

Например, для первых двух Р-значений из приведенного выше примера мы получили бы

• \((0.0001\times15)/1 = 0.0015\)
• \((0.0004\times15/2) = 0.003\)

Воспользовавшись функцией 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, method = "BH")
 [1] 0.00150000 0.00300000 0.00950000 0.03562500 0.06030000
 [6] 0.06385714 0.06385714 0.06450000 0.07650000 0.48600000
[11] 0.58118182 0.71487500 0.75323077 0.81321429 1.00000000

Интерпретация этих Р-значений с поправкой (в большинстве литературных источников их называют q-значениями) такова:

• Допустим, что мы хотим контролировать долю ложно отклоненных гипотез на уровне \(FDR = 0.05\).
• Все гипотезы, q-значения которых \(\leq 0.05\), отклоняются.
• Среди всех этих отклоненных гипотез доля отклоненных по ошибке не превышает 5%.

Как видим, в рассмотренном примере конечный результат после коррекции исходных Р-значений при помощи функции p.adjust() идентичен тому, что был получен при использовании оригинальной процедуры Беньямини-Хохберга: на уровне 5% отклоняются первые четыре гипотезы.


О применимости процедуры Беньямини-Хохберга

Коррекция Р-значений по методу Беньямини-Хохберга работает особенно хорошо в ситуациях, когда необходимо принять общее решение по какому-либо вопросу при наличии информации (=проверенных гипотез) по многим параметрам. Типичным примером будет одновременный анализ многих биологических параметров (вес и температура тела, клеточные показатели крови, и т.п.) в группе пациентов, которых лечили новым препаратом, и в группе, которой давали плацебо. Средние значения каждого параметра в этих группах можно было бы сравнить, например, при помощи t-теста или какого-либо из его непараметрических аналогов. В итоге в распоряжении исследователя оказалось бы большое число соответствующих Р-значений. Общий вывод, который исследователь хотел бы сделать, состоит в том, что новый препарат оказывает положительное влияние на исход лечения. Конечно, в такой ситуации, исследователь был бы заинтересован в обнаружении максимально большого числа параметров, по которым экспериментальные группы различаются (контролируя при этом долю ложных заключений на определенном уровне). Контроль над групповой вероятностью ошибки в этом случае оказался бы слишком строгим (напомню: контролировалась бы вероятность допустить как минимум одну ошибку), тогда как более мощный метод Беньямини-Хохберга допускает наличие определенной доли ложных отклонений среди всех отклоненных гипотез и тем самым способствует общему положительному выводу по поводу эффективности нового препарата.

Следует помнить, что описанный здесь метод контроля над ожидаемой долей ложных отклонений предполагает, что все тесты, при помощи которых получают Р-значения, независимы (Benjamini and Hochberg 1995). На практике в большинстве случаев это условие выполняться не будет (см. пример выше: многие биологические параметры были измерены у одних и тех же испытуемых, что вносит определенный уровень корреляции между соответствующими тестами). Учитывая эту проблему, в 2001 г. И. Беньямини в соавторстве с Д. Йекутили (Benjamini and Yekutieli 2001) предложили усовершенствованный метод, о которым пойдет речь в следующем сообщении.