Бонферрони против Бенджамини-Хохберга: выбор поправки на p-значение.

Время α = 10⁻⁹⁹ было слишком большим: сверхтяжелые элементы и обман.

Делиться

Значения p могут быть деликатной темой. Возможно, лучше избегать их при первом знакомстве со статистиком. Такое отношение к теме привело к негласному согласию, что α = 0,05 является золотым стандартом — по сути, «удобной условностью», эмпирическим правилом, установленным самим Рональдом Фишером.

Кто это?? Не знаете его? Не волнуйтесь.

Он первым ввел метод максимального правдоподобия (MLE), дисперсионный анализ (ANOVA) и информацию Фишера (последнее, как вы, возможно, догадались). Фишер был не просто значимой фигурой в научном сообществе, отцом статистики. Он проявлял глубокий интерес к менделевской генетике и эволюционной биологии, в которые внес ряд ключевых вкладов. К сожалению, у Фишера было и непростое прошлое. Он был связан с Евгеническим обществом и его политикой добровольной стерилизации «умственно отсталых».

Да, не бывает известных статистиков.

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

Существует один ключевой момент, когда вы не просто вынуждены, а обязаны изменить этот альфа-уровень, и все это сводится к проверке множества гипотез .

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

H₁: Солнце голубое

Просто повторяйте наш эксперимент до тех пор, пока не повезет. Но как работают эти поправки? И какую из них следует использовать? Они не взаимозаменяемы!

Значения p и проблема

Чтобы понять почему, нам нужно взглянуть на то, что именно говорит нам наше p-значение. Чтобы понять это глубже, чем просто «маленькое — хорошо, а большое — плохо». Но для этого нам понадобится эксперимент, и ничто не сравнится по захватывающему — и спорному — характеру с открытием сверхтяжелых элементов.

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

Однако их нестабильность становится преимуществом при обнаружении, поскольку новый сверхтяжелый элемент будет демонстрировать отчетливый радиоактивный распад. Последовательность распада, зафиксированная датчиками в реакторе, может показать, присутствует ли новый элемент.

В качестве нулевой гипотезы мы формулируем:

H₀ = Последовательность представляет собой затухание фонового шума. (Нового элемента нет)

Теперь нам нужно собрать доказательства того, что H₀ неверна, если мы хотим доказать, что создали новый элемент. Это делается с помощью нашей тестовой статистики T(X) . В общих чертах, она отражает разницу между тем, что наблюдают датчики, и тем, что ожидается от фонового излучения. Все тестовые статистики являются мерой «неожиданности» между тем, что мы ожидаем наблюдать, если H₀ верна, и тем, что фактически показывают наши выборочные данные. Чем больше T(X) , тем больше у нас доказательств того, что H₀ ложна.

Именно это и делает статистический критерий Шмидта при анализе последовательности времен радиоактивного распада.

[
sigma_{obs} = sqrt{frac{1}{n-1} sum_{i=1}^{n} (ln t_i – overline{ln t})^2}
[

Критерий Шмидта использовался при открытии: хассия (108), мейтнерия (109) в 1984 году, дармштадтия (110), элемента 111 рентгения (111), коперника (112) с 1994 по 1996 год, московского (115), теннессина (117) с 2003 по 2016 год.

Для того чтобы рассчитать вероятность того, что тестовая статистика будет столь же экстремальной, как и тестовая статистика наблюдаемых данных, крайне важно указать распределение для H₀.

Мы предполагаем, что затухание шума подчиняется экспоненциальному распределению . Существует миллион причин, почему это хорошее предположение, но давайте не будем углубляться в это. Если бы у нас не было распределения для H₀, вычисление значения нашей вероятности было бы невозможно!

[
H_0^{(Шмидт)}​:t_1​,…,t_n​ iid ∼ Exp(λ)
[

Значение p представляет собой вероятность получения тестовой статистики, по меньшей мере столь же экстремальной, как и вычисленная на основе выборочных данных, в рамках нулевой модели. Чем менее вероятен наш тестовый показатель, тем вероятнее, что нулевая гипотеза H₀ ложна.

[
p ;=; Pr_{H_0}!big( T(X) ge T(x_{mathrm{obs}}) big).
[

Конечно, это поднимает интересный вопрос. Что, если мы наблюдаем редкую фоновую скорость распада, скорость распада, которая просто напоминает скорость распада неизвестной частицы? Что, если наши датчики обнаруживают маловероятную, хотя и возможную, последовательность распада, которая дает большую статистику теста? Каждый раз, когда мы проводим тест, существует небольшая вероятность получить выброс просто случайно. Этот выброс даст большую статистику теста, поскольку он будет сильно отличаться от того, что мы ожидаем увидеть, когда H₀ истинно. Большое значение T(x) будет находиться в хвостах нашего ожидаемого распределения H₀ и даст малое значение p. Малая вероятность наблюдать что-либо более экстремальное, чем этот выброс. Но никакого нового элемента не существует! Мы просто получили 31 красный, играя в рулетку миллион раз.

Это кажется маловероятным, но если учесть, что протоны направляются на частицы-мишени в течение нескольких месяцев подряд , такая возможность существует. Так как же это объяснить?

Существует два способа: консервативный и менее консервативный. Ваш выбор зависит от эксперимента. Мы можем использовать:

• Показатель ошибки первого рода для всей совокупности критериев (FWER) и поправка Бонферрони.
• Частота ложных обнаружений (FDR) и процедура Бенджамини-Хохберга

Эти варианты не взаимозаменяемы! Вам необходимо тщательно обдумать свой выбор и сделать правильный выбор.

Если вас интересует физика этого явления:

Новые элементы создаются путем ускорения более легких ионов до 10% скорости света. Эти ионные пучки бомбардируют более тяжелые атомы-мишени. Невероятные скорости и кинетическая энергия необходимы для преодоления кулоновского барьера (огромной силы отталкивания между двумя положительно заряженными частицами).

Показатель ошибки первого рода (Бонферрони)

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

[
Pr(text{по крайней мере одна ошибка первого типа в семействе}) leq alpha
[

Это также более простая поправка. Просто разделите уровень значимости альфа на количество проведенных экспериментов. Таким образом, для каждого теста нулевая гипотеза отклоняется тогда и только тогда, когда:

[
p_i leq frac{alpha}{m}
[

Аналогично, вы можете скорректировать свои p-значения. Если вы проводите m тестов, возьмите:

[
p_i^{text{adj}} = min(1, m p_i)
[

И отклонить нулевую гипотезу, если:

[
p_i^{(text{Bonf})} le alpha
[

Здесь мы всего лишь умножили обе стороны неравенства на m .

Доказательство этого также состоит из одной простой строки. Если обозначим событие Aᵢ как ложноположительный результат в тесте i, то вероятность наличия хотя бы одного ложноположительного результата будет равна вероятности объединения всех этих событий.

[
text{Pr}(text{по крайней мере один ложноположительный результат}) = text{Pr}left(bigcup_{i=1}^{m} A_iright) le sum_{i=1}^{m} text{Pr}(A_i) le m cdot frac{alpha}{m} = alpha
[

Здесь мы используем границу объединения. Это фундаментальное понятие в теории вероятностей, которое гласит, что вероятность наступления событий A₁, A₂ или Aₖ должна быть меньше или равна сумме вероятностей наступления каждого из этих событий.

[
text{Pr}(A_1 cup A_2 cup cdots cup A_k) le sum_{i=1}^{k} text{Pr}(A_i)
[

Частота ложных срабатываний (Бенджамини-Хохберга)

Процедура Бенджамини-Хохберга тоже не слишком сложна. Всё просто:

• Отсортируйте значения p: p₁ ≤ … ≤ pₘ.
• Примите первое k, где pₖ ​> α/(m−k+1)

В этом подходе цель состоит в контроле частоты ложных срабатываний (FDR).

[
text{FDR} = Eleft[ frac{V}{max(R, 1)} right]
[

Где R — количество случаев отклонения нулевой гипотезы, а V — количество отклонений, которые (к сожалению) являются ложноположительными (ошибками первого типа). Цель состоит в том, чтобы поддерживать этот показатель ниже определенного порогового значения q = 0,05 .

Пороговые значения BH следующие:

[
frac{1}{m}q, frac{2}{m}q, dots, frac{m}{m}q = q
[

И мы отвергаем первые наименьшие значения p, где:

[
P_{(k)} leq frac{k}{m}q
[

Используйте этот метод, когда вас устраивает некоторое количество ложноположительных результатов. Когда ваша главная задача — минимизировать частоту ошибок второго типа , то есть вы хотите убедиться, что количество ложноотрицательных результатов минимально, и нет случаев, когда мы принимаем H₀, когда H₀ на самом деле ложна.

Представьте это как геномное исследование, цель которого — выявить всех, у кого есть определенный ген, делающий их более восприимчивыми к конкретному виду рака. Было бы менее вредно лечить людей, у которых этого гена нет, чем рисковать тем, что человек, у которого он есть, останется без лечения.

Быстрое сравнение бок о бок

Бонферрони :

• Контролирует уровень ошибок первого рода (FWER).
• Гарантирует вероятность единичного ложного обнаружения ≤ α
• Более высокий процент ложноотрицательных результатов ⇒ Более низкая статистическая мощность
• Нулевая терпимость к риску

Бенджамини-Хохберг

• Контролирует частоту ложных срабатываний (FDR).
• гарантирует, что среди всех обнаруженных случаев количество ложных срабатываний ≤ q
• Меньше ложноотрицательных результатов ⇒ Более высокая статистическая мощность
• Некоторая терпимость к риску

Сверхмаленькая p-структура для сверхтяжелого атома

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

Физики, как правило, используют ожидаемое количество случайных цепочек, полученных в результате всего поиска по всему набору данных:

[
Pr(ge 1 text{ случайная цепочка}) approx 1 – e^{-n_b}
[

[
1 – e^{-n_b} leq alpha_{text{family}} Rightarrow n_b approx alpha_{text{family}} quad (text{приблизительно, для редких событий})
[

Количество случайных цепочек определяется путем наблюдения за фоновыми данными в отсутствие эксперимента. На основе этих данных мы можем построить нулевое распределение H₀ с помощью моделирования методом Монте-Карло.

Мы оцениваем количество случайных цепочек, моделируя фоновые частоты событий и перевыбирая наблюдаемые фоновые события. В условиях H₀ (отсутствие цепочек распада тяжелых элементов) мы используем метод Монте-Карло для моделирования множества нулевых реализаций и вычисляем, как часто алгоритм поиска выдает цепочку, столь же экстремальную, как наблюдаемая.

Точнее:

H₀: фоновые события поступают в виде пуассоновского процесса со скоростью λ ⇒ интервалы между поступлениями имеют экспоненциальное распределение.

Тогда случайная цепочка представляет собой k последовательных попаданий за время τ. Мы сканируем данные, используя нашу тестовую статистику, чтобы определить, существует ли экстремальный кластер.

lambda_rate = 0.2 # событий в секунду T_total = 2_000.0 # секунд сбора данных (среднее количество событий ~ 400) k = 4 # длина цепочки tau_obs = 0.20 # «наблюдаемый экстремум»: 4 события в течение 0.10 сек Nmc = 20_000 rng = np.random.default_rng(0) def dmin_and_count(times, k, tau): if times.size < k: return np.inf, 0 spans = times[k-1:] - times[:-(k-1)] return float(np.min(spans)), int(np.sum(spans <= tau)) ...

Моделирование методом Монте-Карло на GitHub

Если вас интересуют цифры, то при открытии элемента 117 теннессина (Ts) использовалось значение p, равное 5×10−16 . Я предполагаю, что если бы никаких поправок никогда не применялось, наша периодическая таблица, к сожалению, не достигла бы размеров плаката, а химия была бы в полном беспорядке.

Заключение

Вся эта концепция поиска чего-либо во многих местах, а затем обработки особенно значимого всплеска как результата одного наблюдения, обычно называется эффектом «поиска в другом месте» . И есть два основных способа, как мы можем это учесть:

• Исправление Бонферрони
• Процедура Бенджамини-Хохберга

Наш выбор полностью зависит от того, насколько консервативными мы хотим быть.

Но даже при значении p, равном 5×10−16 , вы можете задаться вопросом, когда же значение p, равное 10^-99, все же следует отбросить. И все это связано с Виктором Ниновым, физиком из Национальной лаборатории имени Лоуренса в Беркли. Именно он – на короткое время – открыл элемент 118.

Однако внутреннее расследование показало, что он сфабриковал цепочку альфа-распада. В данном случае, с учетом нарушений в проведении исследований и фальсификации данных, даже значение p, равное 10⁻⁹⁹, не оправдывает отклонение нулевой гипотезы.

Ссылки

Бодмер, В., Бейли, Р.А., Чарльзворт, Б., Эйр-Уокер, А., Фэруэлл, В., Мид, А., и Сенн, С. (2021). Выдающийся ученый Р.А. Фишер: его взгляды на евгенику и расу. Наследственность, 126(4), 565-576.

Хуягбаатар, Й., Якушев, А., Дюльманн, С. Э., Аккерманн, Д., Андерссон, Л. Л., Асаи, М., … и Якушева, В. (2014). Реакция слияния Ca 48+ Bk 249, приводящая к элементу Z= 117: долгоживущий α-распадный Db 270 и открытие Lr 266. Physical review letters, 112(17), 172501.

Положительные результаты, множественные сравнения HMF: поправки Бонферрони и частота ложных срабатываний.

Источник: towardsdatascience.com