Актуальность проблематики

Ложный положительный результат теоретически легко понять и получить

Когда проводится один статистический тест на значимость различий, всегда есть шанс (ошибка первого рода = 5%, на уровне значимости p=0.05) получить ложный положительный результат случайно. Эта ошибка означает, что мы можем ложно утверждать, что значимое различие существует, при том что в реальности этой значимости нет.

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

Предположим, что делается 20 однотипных тестов. Вероятность получить ложный положительный результат равна 1 - (1 - 0.05)^20 ≈ 64%.

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

Ложный положительный результат легко возникает в реальных аналитических задачах.

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

За десять лет наблюдений такой отчет включает в себя по 120-ть однотипных тестов значимости различий индекса NPS (количество промотеров - количество детракторов) как для каждой компании - конкурента, так и ещё 120-ть однотипных тестов для собственной компании. В реальных аналитических отчетах по показателю NPS, как для каждого конкурента, так и для собственной компании заказчика, обычно не применяются никакие методы коррекции ошибок. В таких аналитических отчетах и графиках хотя бы один тест значимых различий для индекса NPS гарантированно имеет для каждой анализируемой компании один ложный положительный результат.

Отчеты очень часто содержат и несколько сотен других показателей, измеряемых ежемесячно, так что очень легко посчитать, что образуются десятки, иногда сотни различий ежемесячно, ложно обозначенных как значимые. Чтобы самостоятельно посмотреть, как быстро вырастет ошибка от количества множественных сравнений, подставьте формулу 1 - (1 - 0.05)^120 в калькулятор. 120 однотипных сравнений — это почти гарантия одной ошибки.

Решения

Конечно, самый точный вариант — это собрать более полный или совсем другой случайный набор данных и убедиться, что никакой значимости различий там нет, только шум. А что делать, если хочется убрать ложные положительные результаты прямо сейчас, на тех данных, что уже доступны, желательно за гарантированное время обработки запроса к озеру данных? Решение, конечно, есть и я предлагаю постепенно его осознать, двигаясь от прошлого к современности. Постараюсь изложить максимально понятно, без формул, а в конце будет рабочий код с примерами решения.

Слабое и простое

Чтобы контролировать общий уровень ошибок (не превышая 5% для всей серии тестов), можно корректировать значение альфа. Одним из первых методов такой поправки является коррекция Бонферрони¹ , при которой уровень значимости p делится на количество тестов.
Для 20 тестов это будет p = 0.05 / 20 = 0.0025 для каждого отдельного теста, чтобы общий уровень ошибок всегда оставался на уровне альфа = 5% = 0.05. Вероятность найти ложный положительный результат (мощность поправки) для этой коррекции падает с увеличением количества гипотез.

На нашем примере, для индекса NPS, это значит, что все значимые различия, найденные в двадцати тестах, для которых уровень значимости теста p меньше 0.05, но больше 0.0025 считаются ложным положительным результатом, а все тесты, для которых p меньше 0.0025 истинно положительными. Методу уже пятдесят+ лет, но он все ещё актуален, если в серии единицы или пара десятков тестов.

Мощнее, сложнее, не всегда применимо

Способ посвежее называется метод Холма²: немного дополненная поправка Бонферрони, где уровни значимости сначала упорядочиваются по возрастанию значения p, p(1) <= p(2) <= ... <= p(i). Корректировки делаются пошагово. Принимается только та гипотеза H(i), чей уровень значимости p(i) первым превысит значение показателя = альфа / (кол-во тестов - i + 1). Дополнительный плюс в том, что мощность этого показателя не падает с увеличением количества гипотез.

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

Например, обманутые не очень добросовестными консультантами секты Дата Дривен Сайентизма менеджеры с энтузиазмом просят найти значимые различия изменения индекса NPS внутри групп по полу, возрасту, доходам и прочих подмножеств, отдельно. Радость клиента от того, что тест без поправки показывает значимое различие в индексе NPS у женщин, 23-35 лет, с высоким доходом, в апреле 2023 года, легко охладить этим методом, когда все p(i) оказываются меньше показателя с поправкой Холма.

Сенсационные данные о том, что именно зеленые конфеты вызывают рак, или, что знак зодиака влияет на спортивные результаты баскетболистов, при том, что на общем массиве данных никакого влияния не видно, пропущенные через поправку Холма, позволяют поймать за руку недобросовестных Дата Дривен Сайентистов или отозвать их научные работы. Подобные манипуляции с данными даже получили отдельное название "p-hacking" или взлом значимости и, конечно, все совпадения с реальностью в данных со взломанной значимостью совершенно случайны, не повторяются при повторном изучении и крайне опасны, особенно в медицинских исследованиях.

Чтобы не теряло мощность и работало всегда, используем теорию вероятности

Относительно недавно появилось понимание, что доля ложных положительных оценок в тестах, при количестве тестов, стремящимся к тысячам, это тоже случайная величина, для которой можно рассчитать математическое ожидание. Появился другой подход - FDR³Бенджамини-Хохберга или, другими словами, вместо метода оценки "групповой вероятности ошибки" был показан метод "контроля ложных открытий" на общем уровне альфа.

И опять уровни значимости всех проведенных тестов упорядочиваются по возрастанию, но теперь считаем, что контролируя во множестве тестов уровень ложных открытий FDR на уровне, не превышающем q* = 0.05 = альфа, ищем такое i, начиная с максимального значения элемента p(i), т.е.с конца, для которого "критическое значение" = q* умножить на (i / кол-во тестов) станет меньше, чем значимость теста i. После этого все гипотезы, для которых значимость различий p больше критического значения - отвергаем как ложные. Считается, что такие вычисления FDR корректно отображают отношение числа ложноположительных результатов к общему числу отвергнутых нулевых гипотез. (FDR = "Число ложноположительных результатов" / "Общее число отвергнутых нулевых гипотез")

Применяя FDR, общее количество ошибок в любом количестве тестов можно удержать в приемлемом объеме - менее 5%.

Очень хорошо работает, подходит для аналитических отчетов, качественно отсекает ложный положительный результат, можно легко удалить из отчета почти все такие значимости. Но не всё так положительно. Применяя FDR необходимо соблюдать условия, что Беты должны быть строго независимы. Что вообще есть какая-то бета или ошибка второго рода вспоминают не очень часто, но про них надо писать отдельную статью.

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

Пермутационное тестирование

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

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

Для проведения тестов на язык R нам понадобится встроенные в R наборы данных iris и mtcars. Подключите библиотеку Datasets, если они недоступны.

Чтобы работал filter() как в примере, установите и подключите пакеты

library(dplyr)
library(tidyr)

 install.packages("dplyr")
 install.packages("tidyr")

 library(dplyr)
 library(tidyr)

 iris_filtered <- filter(iris,iris$Species %in% c("setosa","virginica"))
 observed_t_test  <- t.test(Sepal.Length ~ Species, data = iris_filtered, alternative = "t")
 observed_t_test

Welch Two Sample t-test

data: Sepal.Length by Species
t = -15.386, df = 76.516, p-value < 0.00000000000000022
alternative hypothesis: true difference in means between group setosa and group virginica is not equal to 0
95 percent confidence interval:
-1.78676 -1.37724
sample estimates:
mean in group setosa mean in group virginica
5.006 6.588

Имеем очень неплохую значимость теста на различие средних. 0.00000000000000022 значительно меньше 0.05 похоже, что эти средние совершенно точно разные. Давайте проведем пермутационный тест и посмотрим так ли это. Я буду писать так, чтобы было максимально понятно, что происходит, потом можем порассуждать про оптимизацию и т.д. в комментариях.

Считаем наблюдаемое различие среднего между растениями setosa и virginica по параметру Sepal.Length

 obs_stat <- mean(filter(iris_filtered,iris_filtered$Species == "setosa")$Sepal.Length) - 
             mean(filter(iris_filtered,iris_filtered$Species == "virginica")$Sepal.Length)
 obs_stat

[1] -1.582

Перемешиваем данные и повторяем тест 10000 раз

 n_perm <- 10000
 perm_stat <- numeric(n_perm)
 
 for (i in 1:n_perm) {
      perm_data <- iris_filtered
      perm_data$Species <- sample(perm_data$Species)
      perm_data$Sepal.Length <- sample(perm_data$Sepal.Length)
      perm_stat[i] <- mean(filter(perm_data,perm_data$Species == "setosa")$Sepal.Length) 
                    - mean(filter(perm_data,perm_data$Species == "virginica")$Sepal.Length)
 }

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

hist(perm_stat, main="Гистограмма распределения всех перестановочных значений средних",
      xlab="Разница средних", breaks=100)
 abline(v = obs_stat, col = 'red', lwd = 2, lty = 2)

Получаем вот такое распределение, на котором отчетливо видно, что все значения расположены справа от -1.582, а красная линия так далеко, что даже не появилась на графике.

Рис.1 Нулевое распределение разницы средних на данных iris между растениями setosa и virginica по параметру Sepal.Length

Проверяем это формально

 p_value_adjusted <- sum(perm_stat >= obs_stat) / n_perm
 p_value_adjusted

[1] 1

100% значений находится справа и получается, что высокая значимость теста на различие средних совершенно случайна? Sepal.Length у setosa и virginica в среднем не отличаются?

На первый взгляд это не очень справедливая оценка, ведь разница отрицательная. Давайте лучше проверим по модулю.

p_val <- mean(abs(perm_stat) >= abs(obs_stat))
p_val

[1] 0

Это даёт более корректную оценку - 100% вероятности, что отвергаем нулевую гипотезу о равенстве средних и различие есть.

Вот весь код одним куском. R version 4.4.0. Кол-во пермутаций - n_perm лучше установить больше 100000 для надежности результата.

install.packages("dplyr")
install.packages("tidyr")

library(dplyr)
library(tidyr)

iris_filtered <- filter(iris,iris$Species %in% c("setosa","virginica"))
observed_t_test  <- t.test(Sepal.Length ~ Species, data = iris_filtered, alternative = "t")
observed_t_test
obs_stat <- mean(filter(iris_filtered,iris_filtered$Species == "setosa")$Sepal.Length) - mean(filter(iris_filtered,iris_filtered$Species == "virginica")$Sepal.Length)
obs_stat
n_perm <- 10000
perm_stat <- numeric(n_perm)

for (i in 1:n_perm) {
     perm_data <- iris_filtered
     perm_data$Species <- sample(perm_data$Species)
     perm_data$Sepal.Length <- sample(perm_data$Sepal.Length)
     perm_stat[i] <- mean(filter(perm_data,perm_data$Species == "setosa")$Sepal.Length) - mean(filter(perm_data,perm_data$Species == "virginica")$Sepal.Length)
}

hist(perm_stat, main="Гистограмма распределения всех перестановочных значений средних",
     xlab="Разница средних", breaks=100)
abline(v = obs_stat, col = 'red', lwd = 2, lty = 2)

p_value_adjusted <- sum(perm_stat >= obs_stat) / n_perm
p_value_adjusted

p_val <- mean(abs(perm_stat) >= abs(obs_stat))
p_val

Теперь посмотрим аналогичный случай для mtcars, на данных, для которых значение p побольше.

observed_t_test  <- t.test(mpg ~ am, data = mtcars, alternative = "t")
observed_t_test

Welch Two Sample t-test

data: mpg by am
t = -3.7671, df = 18.332, p-value = 0.001374
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-11.280194 -3.209684
sample estimates:
mean in group 0 mean in group 1
17.14737 24.39231

Считаем наблюдаемое различие среднего по параметру mpg у машин с am == 1 и am == 0

obs_stat <- mean(filter(mtcars,mtcars$am == 1)$mpg) - 
            mean(filter(mtcars,mtcars$am == 0)$mpg)
obs_stat

[1] 7.244939

Разница положительная. Выполняем те-же пермутации

n_perm <- 10000
perm_stat <- numeric(n_perm)
 
for (i in 1:n_perm) {
      perm_data <- mtcars
      perm_data$am <- sample(perm_data$am)
      perm_data$mpg <- sample(perm_data$mpg)
      perm_stat[i] <- mean(filter(perm_data,perm_data$am == 1)$mpg) - mean(filter(perm_data,perm_data$am == 0)$mpg)
}
 
hist(perm_stat, main="Гистограмма распределения всех перестановочных значений средних",
      xlab="Разница средних", breaks=100)
abline(v = obs_stat, col = 'red', lwd = 2, lty = 2)

Рис.1 Нулевое распределение разницы средних на данных mtcars по параметру mpg у машин с am == 1 и am == 0

Уже более понятно, что происходит. Сравнение по модулю и без модуля дает одинаковый ответ, отвергаем нулевую гипотезу о равенстве средних.

p_value_adjusted <- sum(perm_stat = obs_stat) / n_perm
 p_value_ajdusted

[1] 0.0003

p_val <- mean(abs(perm_stat) >= abs(obs_stat))
p_val

[1] 0.0003

Код одним куском

observed_t_test  <- t.test(mpg ~ am, data = mtcars, alternative = "t")
observed_t_test
obs_stat <- mean(filter(mtcars,mtcars$am == 1)$mpg) - mean(filter(mtcars,mtcars$am == 0)$mpg)
obs_stat
n_perm <- 10000
perm_stat <- numeric(n_perm)

for (i in 1:n_perm) {
     perm_data <- mtcars
     perm_data$am <- sample(perm_data$am)
     perm_data$mpg <- sample(perm_data$mpg)
     perm_stat[i] <- mean(filter(perm_data,perm_data$am == 1)$mpg) - mean(filter(perm_data,perm_data$am == 0)$mpg)
}

hist(perm_stat, main="Гистограмма распределения всех перестановочных значений средних",
     xlab="Разница средних", breaks=100)
abline(v = obs_stat, col = 'red', lwd = 2, lty = 2)

p_value_adjusted <- sum(perm_stat >= obs_stat) / n_perm
p_value_adjusted

p_val <- mean(abs(perm_stat) >= abs(obs_stat))
p_val

Теперь посмотрим, что будет, если в данных действительно нет различий. Создадим две группы примерно с одинаковым нормальным распределением значений и протестируем их при помощи пермутаций.

groups<-rep(c(0,1), c(100,200))
null_y_is_true<-rnorm(300)
data <- data.frame (
     groups <- groups,
     null_y_is_true <- null_y_is_true
)

obs_stat <- mean(filter(data,data$groups == 0)$null_y_is_true) - mean(filter(data,data$groups == 1)$null_y_is_true)
obs_stat
n_perm <- 10000
perm_stat <- numeric(n_perm)
for (i in 1:n_perm) {
     perm_data <- data
     perm_data$groups <- sample(perm_data$groups)
     perm_data$null_y_is_true <- sample(perm_data$null_y_is_true)
     perm_stat[i] <- mean(filter(perm_data,perm_data$groups == 0)$null_y_is_true) - mean(filter(perm_data,perm_data$groups == 1)$null_y_is_true)
}

hist(perm_stat, main="Гистограмма распределения всех перестановочных значений средних",
     xlab="Разница средних", breaks=100)
abline(v = obs_stat, col = 'red', lwd = 2, lty = 2)

p_value_adjusted <- sum(perm_stat >= obs_stat) / n_perm
p_value_adjusted

p_val <- mean(abs(perm_stat) >= abs(obs_stat))
p_val

Увидим примерно такую картинку. Примерно потому, что rnorm() будет выдавать каждый раз разные наборы данных.

Рис.3 Нулевое распределение при отсутствии разницы средних.

p_value_adjusted

[1] 0.0739

p_val

[1] 0.1465

Мне остается сделать только несколько финальных замечаний:

Чем больше n_perm, тем больше можно быть уверенным в точности теста. На 100000 пермутаций точность определения для p = 0.05 достигает ±0.1%. А 100000 однотипных тестов можно параллельно запустить в 10 различных потоках, по 10000 тестов в потоке, для гарантированного времени исполнения.

Для того, чтобы в R было легко вносить поправки на множественные сравнения, уже есть встроенная функция p.adjust(), в которую легко передать вектор x, содержащий значения p для нескольких тестов, выбрать один из методов поправки method = "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", # "fdr", "none" часть из которых рассмотрена в этой статье, и получить значения p, уже с поправкой на количество одновременных сравнений. Можете попробовать, для тренировки, применить поправку следующим образом:

iris_filtered <- filter(iris,iris$Species %in% c("setosa","virginica"))
observed_t_test  <- t.test(Petal.Length  ~ Species, data = iris_filtered, alternative = "t")
x  <- list(observed_t_test$p.value)

iris_filtered <- filter(iris,iris$Species %in% c("setosa","versicolor"))
observed_t_test  <- t.test(Petal.Length ~ Species, data = iris_filtered, alternative = "t")
x  <- c(x,observed_t_test$p.value)

iris_filtered <- filter(iris,iris$Species %in% c("versicolor","virginica"))
observed_t_test  <- t.test(Petal.Length ~ Species, data = iris_filtered, alternative = "t")
x  <- c(x,observed_t_test$p.value)
as.numeric(x)

[1] 9.269628e-50 9.934433e-46 4.900288e-22

adjusted_p_values <- p.adjust(x,method = "bonferroni")
adjusted_p_values

[1] 2.780888e-49 2.980330e-45 1.470086e-21

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

library(perm) library(multtest)

Про обязательность поправки на множественные сравнения, которая часто игнорируется адептами Data Driven методов