13 октября 2013

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



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



Групповая вероятность ошибки (первого рода)

Прежде чем начать рассмотрение процедур множественной проверки статистических гипотез, стоит разобраться с тем, для чего именно эти процедуры служат с математической точки зрения. Предположим, что мы проверяем истинность \(m\) нулевых гипотез, которые можно обозначить как \(H_1, H_2, H_3, \dots, H_m\). В отношении каждой из этих гипотез мы применяем определенный статистический критерий (например, t-критерий Стьюдента) и делаем заключение о том, верна ли она (т.е. выполняем \(m\) одинаковых тестов на одних и тех же данных и либо принимаем, либо отвергаем каждую гипотезу). Результаты такого анализа можно свести в таблицу следующего вида:

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

В этой таблице
  • \(m_0\) - число верных нулевых гипотез (англ. true null hypotheses); 
  • \(m - m_0\) - число истинных альтернативных гипотез (true alternative hypotheses); 
  • \(U\) - число безошибочно принятых гипотез (true negatives); 
  • \(V\) - число ошибочно отвергнутых гипотез (false positives) (ошибка первого рода); 
  • \(T\) - число ошибочно принятых гипотез (false negatives) (ошибка второго рода); 
  • \(S\) - число безошибочно отвергнутых гипотез; 
  • \(W\) - общее число принятых гипотез; 
  • \(R\) - общее число отвергнутых гипотез.
Представленную таблицу следует читать следующим образом. В первой строке мы видим, что всего имеется \(m_0\) верных нулевых гипотез. В ходе выполнения анализа определенная их часть ошибочно отвергается (\(V\)), тогда как остальные \(U\) гипотез классифицируются правильно. Аналогично, во второй строке мы видим, что всего имеется \(m - m_0\) альтернативных гипотез, из которых в ходе анализа \(S\) гипотез безошибочно отвергаются, а \(T\) гипотез - ошибочно принимаются. Обратите внимание на то, что общие количества отвергнутых (\(R\)) и принятых (\(W\)) гипотез нам известны, тогда как \(m_0\), \(U\), \(V\), \(T\) и \(S\) - ненаблюдаемые случайные величины.

Таким образом, при одновременной проверке группы (или "семейства", англ. family) статистических гипотез задача заключается в том, чтобы минимизировать число ложных отклонений \(V\) и ложных принятий \(T\). Традиционно исследователи пытаются минимизировать величину \(V\). Если \(V \geq 1 \), мы совершаем как минимум одну ошибку первого рода. Вероятность допущения такой ошибки при множественной проверке гипотез называют "групповой вероятностью ошибки" (англ. "familywise error rate", FWER; реже используется термин "experiment-wise error rate"). По определению, \(FWER = P(V \geq 1)\). Соответственно, когда мы говорим, что хотим контролировать групповую вероятность ошибки на определенном уровне значимости \( \alpha \), мы подразумеваем, что должно выполняться неравенство \(FWER \leq \alpha\). Методы множественной проверки гипотез как раз и позволяют обеспечивать этот контроль.


Поправка Бонферрони

Метод Бонферрони (назван так в честь предложившего его итальянского математика Карло Эмилио Бонферрони; Carlo Emilio Boferroni) является одним из наиболее простых и известных способов контроля над групповой вероятностью ошибки.

Предположим, что мы применили определенный статистический критерий 3 раза (например, сравнили при помощи критерия Стьюдента средние значения групп А и В, А и С, и В и С) и получили следующие три Р-значения: 0.01, 0.02 и 0.005. Если мы хотим, чтобы групповая вероятность ошибки при этом не превышала определенный уровень значимости \(\alpha\) (например, 0.05), то, согласно методу Бонферрони, мы должны сравнить каждое из полученных Р-значений не с \(\alpha\), а с \(\alpha/m\), где \(m\) - число проверяемых гипотез. Деление исходного уровня значимости \(\alpha\) на \(m\) - это и есть поправка Бонферрони. В рассматриваемом примере каждое из полученных Р-значений необходимо было бы сравнить с 0.05/3 = 0.017. В результате мы выяснили бы, что Р-значение для второй гипотезы (0.02) превышает 0.017 и, соответственно, у нас не было бы оснований отвергнуть эту гипотезу.

Вместо деления изначально принятого уровня значимости на число проверяемых гипотез, мы могли бы умножить каждое из исходных Р-значений на это число. Сравнив такие скорректированные Р-значения (англ. adjusted P-values; обычно обозначаются буквой q) с \(\alpha\), мы пришли бы к точно тем же выводам, что и при использовании поправки Бонферрони:
  • 0.01 * 3 = 0.03 < 0.05: гипотеза отклоняется;
  • 0.02 * 3 = 0.06 > 0.05: гипотеза принимается;
  • 0.005 * 3 = 0.015 < 0.05: гипотеза отклоняется.
В ряде случаев при умножении исходных Р-значений на уровень значимости результат может превысить 1. По определению, вероятность не может быть больше 1, и если это происходит, то получаемое значение просто приравнивают к 1.

В базовой версии R имеется функция p.adjust(), которая позволяет легко применять как метод Бонферрони, так и ряд других методов, обеспечивающих контроль групповой вероятности ошибки на определенном уровне (задаются при помощи аргумента method). При реализации метода Бонферрони с помощью этой функции применяется второй из описанных выше подходов, т.е. исходные Р-значения умножаются на число проверяемых гипотез. Функция p.adjust() принимает вектор с исходными Р-значениями и возвращает скорректированные значения:

# Скорректированные Р-значения:
p.adjust(c(0.01, 0.02, 0.005), method = "bonferroni")
[1] 0.030 0.060 0.015
 
# Какие из проверяемых гипотез следует отвергнуть?
alpha <- 0.05
p.adjust(c(0.01, 0.02, 0.005), method = "bonferroni") < alpha
[1]  TRUE FALSE  TRUE

Хотя метод Бонферрони очень прост в реализации, он обладает одним существенным недостатком: при возрастании числа проверяемых гипотез мощность этого метода резко снижается. Другими словами, при возрастании числа гипотез нам будет все сложнее и сложнее отвернуть многие из них, даже если они неверны и должны быть отвергнуты. Например, при проверке 10 гипотез, применение поправки Бонферрони привело бы к снижению исходного уровня значимости до 0.05/10 = 0.005. Соответственно, для отклонения той или иной гипотезы, соответствующие Р-значения должны были бы оказаться меньше 0.005, что случалось бы нечасто. В связи с этим метод Бонферрони не рекомендуется использовать, если число проверяемых гипотез превышает 7-8.


Метод Холма

Для преодоления проблем, связанных с низкой мощностью метода Бонферрони, в 1978 г. Стур Холм (Holm 1978) предложил гораздо более мощную его модификацию (часто этот метод называют еще методом Холма-Бонферрони). Этот модифицированный метод основан на алгоритме, который включает следующие шаги:
  • Исходные Р-значения упорядочиваются по возрастанию: \(p_{(1)} \leq p_{(2)} \leq \dots \leq p_{(m)} \). Эти Р-значения соответствуют проверяемым гипотезам \(H_{(1)}, H_{(2)}, \dots H_{(m)} \).
  • Если \(p_{(1)} \geq \alpha/m\), все нулевые гипотезы \(H_{(1)}, H_{(2)}, \dots H_{(m)}\) принимаются и процедура останавливается. Иначе следует отвергнуть гипотезу \(H_{(1)}\) и продолжить.
  • Если \(p_{(2)} \geq \alpha/(m - 1)\), нулевые гипотезы \(H_{(2)}, H_{(3)}, \dots H_{(m)}\) принимаются и процедура останавливается. Иначе гипотеза \(H_{(2)}\) отвергается и процедура продолжается.
  • ...
  • Если \(p_{(m)} \geq \alpha\), нулевая гипотеза \(H_{(m)}\) принимается и процедура останавливается.
Описанную процедуру называют нисходящей (англ. step-down): она начинается с наименьшего P-значения в упорядоченном ряду и последовательно "спускается" вниз к более высоким значениям. На каждом шаге соответствующее значение \(p_{(i)}\) сравнивается со скорректированным уровнем значимости \( \alpha/(m + i -1) \). Как и в случае с поправкой Бонферрони, вместо корректировки уровня значимости мы можем скорректировать непосредственно Р-значения - конечный результат (в смысле принятия или отклонения той или иной гипотезы) окажется идентичным. Соответствующая поправка выполняется в виде \( q_i = p_{(i)} (m + i -1) \). Так, для рассмотренного выше примера с тремя Р-значениями получаем:
  • \(q_1 = p_{(1)}(m - 1 + 1) = 0.005*3 =  0.015\)
  • \(q_2 = p_{(2)}(m - 2 + 1) = 0.01*2 =  0.02\)
  • \(q_3 = p_{(3)}(m - 3 + 1) = 0.02*1 =  0.02\)
Именно последний подход реализован в R-функции p.adjust():

# Скорректированные Р-значения:
p.adjust(c(0.01, 0.02, 0.005), method = "holm")
[1] 0.020 0.020 0.015
# (обратите внимание: функция автоматически упорядочивает
# итоговые Р-значения по убыванию)
 
# Какие из проверяемых гипотез следует отвергнуть?
alpha <- 0.05
p.adjust(c(0.01, 0.02, 0.005), method = "holm") < alpha
[1] TRUE TRUE TRUE

Сравните результаты, полученные при помощи методов Бонферрони и Холма: в последнем случае мы отвергаем все три гипотезы, тогда как при использовании поправки Бонферрони отклонены были бы только две из трех проверяемых гипотез. Этот простой пример подтверждает, что метод Холма обладает большей мощностью, чем метод Бонферрони.

В следующем сообщении я опишу методы, используемые для сравнения групповых средних в ходе выполнении дисперсионного анализа (в частности, критерий Тьюки).

21 комментарий :

Анонимный комментирует...

Здравствуйте, Сергей. А если предположить ситуацию, когда мы сравниваем не 3 группы A, B и C по некоему 1 параметру, как в вашем примере, а по 3 параметрам, но при этом групп всего 2. Такая ситуация, вероятно, также предполагает применение поправки?

Sergey Mastitsky комментирует...

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

Анонимный комментирует...

Спасибо, Сергей!

Анонимный комментирует...

Сергей, спасибо Вам за очередную полезную статью.
Возможно я ошибаюсь, но по моему у Вас несоответствие между обозначениями гипотез в таблице и обозначениями в тексте.

Вот один из кусков текста: "Таким образом, при одновременной проверке группы (или "семейства", англ. family) статистических гипотез задача заключается в том, чтобы минимизировать число ложных отклонений V и ложных принятий T.".
Неверно отклоненные гипотезы в таблице обозначены буквой S, а в тексте они же обозначаются буквой V (что на самом деле, исходя из таблицы - верно отвергнутые гипотезы). В других частях текста также есть эти несоответствия.

Sergey Mastitsky комментирует...

Здравствуйте! Спасибо за коментарий.
Понимаю, что приведенная таблица может быть сложна для понимания, но нет - все обозначения в тексте и в таблице верны. Перечитайте еще раз список обозначений и сразу следующий за ним абзац. Ключ к пониманию состоит в том, что таблицу нужно "читать" по строкам. В первой строке речь идет о верных нулевых гипотезах (например, "mu_1 = mu_2"). Если эти гипотезы отвергаются, совершается ошибка 1-го рода (с частотой V). Во второй строке информация касается альтернативных гипотез ("mu_1 не равно mu_2"). Если мы их сохраняем (=не отвергаем, как должны были бы сделать), то совершаем ошибку 2-го рода (с частотой Т). Конечно, на практике мы хотим минимизировать обе ошибки, но обычно фокус направлен на ошибки первого рода.

Анонимный комментирует...

А, теперь понял! Благодарю за пояснение.
В таком случае, прошу прощения за путаницу.

Evgen Antonov комментирует...

Здравствуйте! Подскажите пожалуйста, а как применять метод Холма, если имеются два одинаковых Р-значения, например (0,0016, 0,0016, 0,006). Спасибо!

Sergey Mastitsky комментирует...

Evgen, не очень понимаю, в чем заключается Ваш вопрос. Процедура применения этой поправки полностью описана в сообщении - пожалуйста, уточните, что именно вызывает затруднение.

Анонимный комментирует...

Здравствуйте!
>> (например, сравнили при помощи критерия Стьюдента средние значения групп А и В, А и С, и В и С)

Подскажите, пожалуйста, в этом случае применяется one- или two-sample t-test?

Sergey Mastitsky комментирует...

Вид теста не имеет значения (т.е. one- или two-sided) - логика остается той же. Хотя в примере речь идет о двухстороннем тесте.

Анонимный комментирует...

А в применении к языку R: средние не считаются отдельно, а берется t.test для каждой пары АВ, ВС, АС?

Sergey Mastitsky комментирует...

Как работает t-test в R подробно описано здесь: http://r-analytics.blogspot.de/2012/03/t.html

Анонимный комментирует...

Если для попарного сравнения групп данных использовались разные двухвыборочные критерии (параметрический и непараметрический), то можно ли объединить p-value для всех этих гипотез и проводить поправку (Бонферрони, FDR или другие)? Дело в том, что одна из групп данных имеет распределение отличное от нормального. Как Вы считаете, как можно поступить в этом случае для получения корректных выводов по результатам анализа данного исследования? Спасибо. Марина

Sergey Mastitsky комментирует...

Марина, чтобы сравнивать "яблоки с яблоками", а не "яблоки с апельсинами", объединять р-значения, полученные при помощи разных методов, не стоит. Если условия применения параметрического теста не выполняются и с этим ничего поделать нельзя (например, путем трансформации данных), я бы посоветовал применять непараметрический тест во всех случаях и затем применять выбранную поправку к р-значениям. При прочих равных условиях, мощность непараметрических критериев несколько ниже, чем у параметрических, но в большинстве случаев (особенно при наличии хорошего "сигнала" в данных) это не будет большой проблемой.

Unknown комментирует...

Сергей, спасибо за статью и за сайт в целом.
У меня вопрос немного в сторону.

Предположим, мы провели серию из нескольких экспериментов. Принимаем порог статистической значимости 0.05. В каждом из экспериментов мы получили p = 0.04.

Если мы используем поправку Холма-Бонферрони (или Бонферрони), то каждый из экспериментов дал НЕЗНАЧИМЫЙ результат.

Но все-таки интуитивно-то понятно, что если у нас серия из 100 экспериментов, в каждом из них p = 0.04, то только в 5 экспериментах такая низкая p -- это результат случая, а в 95 экспериментах p = 0.04 -- это значимый результат.

Нет?

Спасибо.
Алексей.

Sergey Mastitsky комментирует...

Алексей, в том-то все и дело, что допустимая вероятность ошибки в КАЖДОМ эксперименте будет 5% (а не в 5% случаев из числа выполненных экспериментов) и, как следствие, общая ошибка окажется гораздо выше 5%. Методы Бонферрони или Бонферрони-Холма в таких ситуациях будут неэффективными (т.е. проявят очень низкую мощность). Вам нужно это: http://r-analytics.blogspot.de/2013/11/blog-post.html
или это: http://r-analytics.blogspot.de/2013/11/blog-post_11.html

Unknown комментирует...

Сергей,
Давайте рассмотрим не такую уж и гипотетическую ситуацию. Экспериментатор А получил данные, стистическая обработка которых показала, что введение препарата Х повлияло на измеряемый параметер (скажем, температуру тела), а последующий статистический тест показал, что p = 0.04. Другой экспериментатор Б, совершенно независимо получил такие же точно результаты. Оба ученых опубликовали свои результаты в престижных научных изданиях. Третий ученый В нашел эти публикации и вставил эти результаты в свой обзор, сделав вывод, что препарат Х статистически достоверно (p = 0.04) влияет на температуру тела. Однако, опубликовать свой обзор ученый В не смог, поскольку не прошел критику референта: для множественного теста (два статистических теста, по одному из каждой оригинальной статьи) должена быть использована попрака типа Бонферони, а если это сделать, то эффект оказывается нелостоверным. Что Вы об этом думаете, Сергей? Спасибо, Владимир.

Vladimir Tsibulsky комментирует...

Cергей,
Вот еще одно соображение.
Предположим, что мы сравнили 3 группы А1:А2, А1:А3, А2:А3 и получили, соответственно, следующие значения для р: 0.016, 0.024 и 0.049. Установив альфа = 0.05 и используя метод Холма-Бонферрони мы получаем:
1. альфа/3 = 0.017, следовательно первый тест выявил достоверное различие А1:А2
2. альфа/2 = 0.025, следовательно второй тест выявил достоверное различие А1:А3
3. альфа/1 = 0.050, следовательно третий тест выявил достоверное различие А2:А3
А теперь предположим, что, вместо указанных выше, получились следующие значения для р: 0.017, 0.024 и 0.049. Используя тот же метод Холма-Бонферрони мы получаем:
1. альфа/3 = 0.0167, следовательно первый тест не выявил достоверное различие в сравнении А1:А2, что также приводит к недостовености различий между группами А1:А3 и А2:А3.
Сравнение этих двух случаев показывает, что весьма скромное (на 5%) уменьшение статистической достоверности различий между группами с наибольшим из всех эфектов приводит к заключению о полном отсутствии достоверных различий во всех остальных сравнениях (а их может быть много больше двух). Мне кажется, что-то в этой картине не правильно. А вы как думаете, Сергей? Заранее спасибо за ответ. Владимир.
P.S. Прошу прощения за пару орфографических опечаток в моем первом комментарии.

Sergey Mastitsky комментирует...

Уважаемый Владимир,

Вопросов Вы задали много и для полного ответа на них пришлось бы написать целую статью. Тем не менее, вот несколько общих соображений:

1) По поводу ученых опубликовавших свои статьи в "престижных журналах":
Как таковое p-значение совершенно ни о чем не говорит, в частности оно ничего не говорит о величине эффекта исследуемого экспериментального фактора, а именно она представляет наибольший интерес. Поэтому третий ученый в своей статье и корреспонденции с референтом должен описывать и использовать в качестве аргументации не только p-значения, но и такие вещи, как объем задействованных выборок, размер достигнутого эффекта и мощность теста. Кроме того, все методы поправок исходных p-значений предполагают, что от эксперимента к эксперименту выборки имеют фиксированный объем и происходят из одних и тех же генеральных совокупностей. Так ли это во всех трех случаях с учеными А, В и С?
Ну и, конечно, кто сказал, что в качестве критического значения P нужно использовать 5%, а не, скажем, 10%? Любое пороговое р-значение является произвольным и используется нами просто для удобства. В этой связи давно нужно переходить с р-значений на оперирование доверительными интервалами, а ученому C для обобщения разных экспериментов стоит ознакомиться с методами мета-анализа. В своей корреспондеции с референтом ученому С неплохо было еще задать вопрос о том, почему "должна" использоваться именно поправка Бонферрони, а не какой-то другой метод? Никаких строгих правил для выбора метода корректировки р-значений просто не существует.
Тем не менее, при соблюдении определенных условий (см. выше), р-значения могут быть полезными в качестве инструмента, помогающего находить "сигнал" в данных. Поэтому я описал наиболее распространенные методы корректировки р-значений. Как и в случае с любым другим инструментом, пользователь должен хорошо понимать, для чего этот интструмент предназначен и в чем состоят его ограничения.

2) По поводу метода Холма:
Отвечу вопросом на вопрос, уж не обессудьте: получив р-значение, очень близкое к пороговому, станете ли вы доверять своим результатам и строить целую теорию на их основе? Или все же попробуете набрать дополнительный материал, чтобы удостовериться в воспроизводимости результата?
См. также выше по поводу принятого уровня значимости (который всегда являтся произвольным), использованного объема выборок и желаемых мощности теста и величине эффекта.

Sergei комментирует...

> Ну и, конечно, кто сказал, что в качестве критического значения P нужно использовать 5%, а не, скажем, 10%?

Это уровень вероятности который мы принимаем за уровень практически недостоверного события. Настолько можно доверять доказываемым нами суждениям. Если 19 из 20 гипотез проверенных при таком подходе окажутся истинными, то это считается достаточным для биологии с её высокой изменчивостью и трудностью учета-ограничения влияющих на результат факторов.

Артем комментирует...

В том то и дело, что уровень 5% совершенно не достаточен, для того, чтобы доверять полученным результатам, особенно в биологии. Если у одного врача умирает каждый 20 пациент, является ли это случайностью? Невозможно ответить на этот вопрос без дополнительных данных. Но в науке вдруг кем-то было заведено, совсем не автором P, что p < 0.05 это критерий истины.

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