nnazarov 11 апр в 15:23

Мифы о байесовском А/Б тестировании

11 мин

5.1K

Блог компании X5 TechPython*Математика*Статистика в IT

Хабр, привет! Сегодня сравним два подхода к А/Б тестированию: байесовский и частотный. Обсудим сложности в интерпретации p-value. Посмотрим, как можно учитывать дополнительную информацию через априорное распределение. Остановим тест раньше времени и решим проблему подглядывания.

Меня зовут Коля, я работаю аналитиком данных в X5 Tech. Мы с Сашей продолжаем писать серию статей по А/Б тестированию. Предыдущие статьи можно найти в описании профиля.

Постановка задачи

Определим задачу и требования к решению, чтобы можно было сравнивать разные подходы. Допустим, проверяем гипотезу о равенстве средних. Хотим определить, есть ли статистически значимые отличия. Если есть, то внедряем тестируемое изменение, иначе нет. В этом случае можно совершить одну из двух ошибок:

ошибку первого рода — сказать, что эффект есть, когда на самом деле его нет;
ошибку второго рода — сказать, что эффекта нет, когда на самом деле он есть.

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

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

Частотный подход

В частотном подходе, зная распределение статистики теста при верности нулевой гипотезы, по данным эксперимента вычисляется p-value. Если p-value меньше уровня значимости, то отклоняем нулевую гипотезу, иначе нет.

Для удобства сравнения с байесовским подходом рассмотрим проверку односторонней гипотезы о равенстве средних для данных из распределения Бернулли. Зафиксируем параметры эксперимента:

Уровень значимости $\alpha=0.05$
Допустимая вероятность ошибки второго рода $\beta=0.1$
Ожидаемый эффект $\epsilon=0.1$
Доля единиц в данных контрольной группы

Оценим размер групп, необходимый для получения заданной допустимой вероятности ошибки второго рода по формуле:

$n > \frac{\left[ \Phi^{-1} \left( 1-\alpha \right) + \Phi^{-1} \left( 1-\beta \right) \right]^2 (\sigma_X^2 + \sigma_Y^2)}{\varepsilon^2}$

from collections import defaultdict
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.proportion import proportions_ztest

alpha, beta = 0.05, 0.1    # вероятности ошибок
p1 = 0.3                   # исходная конверсия
effect = 0.1               # эффект
p2 = p1 + effect           # конверсия с эффектом

sum_var = p1 * (1 - p1) + p2 * (1 - p2)
t_alpha = stats.norm.ppf(1 - alpha, loc=0, scale=1)
t_beta = stats.norm.ppf(1 - beta, loc=0, scale=1)
size = int(np.ceil((t_alpha + t_beta) ** 2 * sum_var / effect ** 2))
print(size)

Напишем функцию, которая проверяет гипотезу частотным подходом. Приведём пример её применения.

def frequentist_test(sum_a, len_a, sum_b, len_b, alpha=0.05):
    pvalue = proportions_ztest(
        [sum_a, sum_b], [len_a, len_b], alternative="smaller"
    )[1]
    decision = int(pvalue < alpha)
    return pvalue, decision

sum_a = np.random.binomial(size, p1)
sum_b = np.random.binomial(size, p2)

pvalue, decision = frequentist_test(sum_a, size, sum_b, size)
print(f'sum_a={sum_a}, sum_b={sum_b}, pvalue={pvalue:0.3f}, decision={decision}')

sum_a=126, sum_b=154, pvalue=0.018, decision=1

Байесовский подход

В байесовском подходе оцениваемые параметры считают случайными величинами. Делается предположение, что они из некоторого априорного распределения. Далее численно или аналитически, используя данные эксперимента, получают оценку совместного апостериорного распределения. Чтобы принять решение, например, можно по совместному распределению оценить вероятность того, что среднее одной группы больше среднего другой группы. Если эта вероятность больше $1-\alpha$ , то говорим, что отличия значимы.

Для гипотезы о равенстве средних оцениваемые параметры — это средние значения контрольной и экспериментальной групп ( $\theta_A$ и $\theta_B$ ). Для данных из распределения Бернулли в качестве априорного распределения часто используют бета-распределение, так как с ним легко аналитически получить апостериорное распределение. Параметры бета-распределения (u,v) для априорной вероятности обычно берут равными (1, 1) или (0.5, 0.5) – их называют uninformative prior.

Если априорное распределение $P(\theta)\sim Beta(u, v)$ и за время эксперимента было получено наблюдений, из которых было единиц, то апостериорное распределение имеет вид $P(\theta |n, k)\sim Beta(u+k, v+(n-k))$ . Совместное распределение средних двух групп имеет вид $P(\theta_A, \theta_B |X) = P(\theta_A |X) * P(\theta_B |X)$ , где — данные эксперимента. На графике ниже приведены примеры плотностей бета-распределений для разных вариантов параметров распределения.

Напишем функцию, которая проверяет гипотезу байесовским подходом. Эта функция дополнительно принимает параметры априорных распределений и количество данных, которые будут сгенерированы для численной оценки вероятности. Приведём пример её применения.

def bayesian_test(
        sum_a, len_a, sum_b, len_b, prior_a, prior_b,
        n=10000, alpha=0.05
):
    posterior_a = [prior_a[0] + sum_a, prior_a[1] + len_a - sum_a]
    posterior_b = [prior_b[0] + sum_b, prior_b[1] + len_b - sum_b]
    values_a = np.random.beta(*posterior_a, size=n)
    values_b = np.random.beta(*posterior_b, size=n)
    prob_b_more_a = (values_b > values_a).mean()
    decision = int(prob_b_more_a > 1 - alpha)
    return prob_b_more_a, decision

prior_a = [1, 1]
prior_b = [1, 1]

prob, decision = bayesian_test(sum_a, size, sum_b, size, prior_a, prior_b)
print(f'sum_a={sum_a}, sum_b={sum_b}, prob_b_more_a={prob:0.3f}, decision={decision}')

sum_a=126, sum_b=154, prob_b_more_a=0.983, decision=1

Корректность критериев

Проверим, что критерии работают корректно. Проведём $10\ 000$ синтетических АА и AB экспериментов, оценим вероятности ошибок и построим доверительные интервалы. Убедимся, что они приблизительно равны значениям $\alpha$ и $\beta$ . Подробнее про проверку корректности можно прочитать в этой статье.

def show_bernulli_ci(dict_decision):
    for key, values in dict_decision.items():
        pe, std = np.mean(values), np.std(values)
        mean_std = std / len(values) ** 0.5
        left, right = pe - 1.96 * mean_std, pe + 1.96 * mean_std
        print(f'{key}: pe={pe:0.3f}, ci=[{left:0.3f}, {right:0.3f}]')

res = defaultdict(list)
for _ in range(10000):
    sum_a, sum_b = np.random.binomial(size, p1, 2)
    sum_b_effect = np.random.binomial(size, p2)
    res['freq_aa'].append(frequentist_test(sum_a, size, sum_b, size)[1])
    res['freq_ab'].append(frequentist_test(sum_a, size, sum_b_effect, size)[1])
    res['bayes_aa'].append(bayesian_test(sum_a, size, sum_b, size, prior_a, prior_b)[1])
    res['bayes_ab'].append(bayesian_test(sum_a, size, sum_b_effect, size, prior_a, prior_b)[1])
show_bernulli_ci(res)

freq_aa: pe=0.051, ci=[0.047, 0.055]
freq_ab: pe=0.899, ci=[0.893, 0.905]
bayes_aa: pe=0.051, ci=[0.046, 0.055]
bayes_ab: pe=0.899, ci=[0.893, 0.905]

Методы работают корректно.

Байесовский подход лучше

Периодически доходят слухи о том, что байесовский подход лучше частотного. Вот некоторые аргументы:

P-value в частотном подходе сложно интерпретировать.
В байесовском А/Б тестировании можно использовать дополнительную информацию через априорное распределение.
Можно остановить эксперимент раньше.
Можно проверять значимость отличий в любой момент.
Удобнее проводить множественное тестирование, можно сравнивать результаты без коррекций.

Рассмотрим каждый из аргументов подробнее.

«P-value сложно интерпретировать»

Нам нужно определить, есть ли значимые отличия между группами, чтобы принять бинарное решение: внедряем или не внедряем изменение. Для этого p-value очень удобен, достаточно сравнить его с уровнем значимости. Если p-value меньше уровня значимости, то отличия между группами статистически значимые, отвергаем нулевую гипотезу. Если p-value больше уровня значимости, то статистически значимых отличий не обнаружено. Такое решающее правило легко донести любому.

При желании погрузиться глубже в понимании p-value нужно прочитать и осознать его определение. P-value — вероятность получить такое же или более экстремальное значение статистики при верности нулевой гипотезы. Разберём определение по пунктам:

Проверяем нулевую гипотезу, например, о равенстве средних.
Для проверки гипотезы используется некоторый статистический критерий, в котором определена статистика критерия. Например, используем тест Стьюдента, в котором используется t-статистика.
По данным эксперимента вычисляем значение t-статистики, получаем конкретное число.
Для теста Стьюдента знаем функцию распределения t-статистики при верности нулевой гипотезы. По функции распределения посчитаем вероятность получить значение как в проведённом эксперименте или более экстремальное, то есть более удалённое от нуля.

Как видите, в p-value нет ничего сложно. Ниже приведена картинка, иллюстрирующая вычисление p-value.

Иногда говорят, что 1 - pvalue — это не вероятность того, что среднее в одной группе больше, чем среднее в другой группе, а хотелось бы вероятность $\mathbb{P}(m_B > m_A)$ . Действительно, — это не вероятность того, что среднее в одной группе больше, чем среднее в другой группе. Что такое p-value, объяснили несколькими абзацами выше.

В желании работать с вероятностями типа $\mathbb{P}(m_B > m_A)$ нет ничего плохого. Заметим лишь, что мнение о том, что люди хорошо оценивают вероятности в уме может быть преувеличено. Один может сказать, что вероятности $\mathbb{P}(m_B > m_A)=0.6$ достаточно, чтобы принять решение о значимости отличий, ведь это в полтора раза больше, чем $\mathbb{P}(m_B < m_A)=0.4$ . А другой, что вероятности $\mathbb{P}(m_B > m_A)=0.97$ недостаточно, ведь вероятность ошибки – целых 3%. В итоге при грамотном подходе всё сведётся к контролю вероятностей ошибок и составления решающего правила аналогично критерию для p-value.

«Дополнительные знания через prior»

Сторонники байесовского А/Б тестирования говорят, что байесовский подход позволяет учитывать дополнительные знания через априорное распределение. Рассмотрим такой пример. Допустим, конверсия до эксперимента была 0.3. Используем это знание и возьмём в качестве априорного распределения для обеих групп Beta(30, 70) . Оценим вероятности ошибок.

p1 = 0.3
effect = 0.1
p2 = p1 + effect
prior_a = [30, 70]
prior_b = [30, 70]
size = 385

res = defaultdict(list)
for _ in range(10000):
    sum_a, sum_b = np.random.binomial(size, p1, 2)
    sum_b_effect = np.random.binomial(size, p2)
    params_aa = [sum_a, size, sum_b, size, prior_a, prior_b]
    res['bayes_aa'].append(bayesian_test(*params_aa)[1])
    params_ab = [sum_a, size, sum_b_effect, size, prior_a, prior_b]
    res['bayes_ab'].append(bayesian_test(*params_ab)[1])
show_bernulli_ci(res)

bayes_aa: pe=0.028, ci=[0.025, 0.032]
bayes_ab: pe=0.859, ci=[0.852, 0.866]

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

Рассмотрим другой случай. Допустим, ожидаем, что конверсия в экспериментальной группе увеличится до 0.4. Учтём это ожидание и возьмём в качестве априорного распределения экспериментальной группы Beta(40, 60) (prior_b = [30, 70]). Проведя аналогичный эксперимент, получим следующие оценки вероятностей ошибок:

bayes_aa: pe=0.141, ci=[0.134, 0.148]
bayes_ab: pe=0.966, ci=[0.962, 0.969]

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

Из экспериментов видно, что добавление информации через априорное распределение меняет баланс вероятностей ошибок, увеличивает одну и уменьшает другую. Посмотрим на суть добавления информативного априорного распределения через формулы. Напомним, если — количество наблюдений, — количество единиц, то для априорного распределения Beta(u, v) апостериорное распределение имеет вид Beta(u+k, v+(n-k)) . Преобразуем параметры апостериорных распределений для разных априорных распределений:

$\begin{aligned} \text{При } & u=1,\ v=1: P(𝜃|n,k) \sim Beta(1+k, 1+(n-k)) \\ \text{При } & u=30,\ v=70: P(𝜃|n,k) \sim Beta (30+k, 70+(n-k)) \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = Beta (1+(k+29), 1+(n-k+69)) \end{aligned}$

Замена априорного распределения привела к переходу $k \to k+29$ и $n \to n+69$ . Получается, что использовать информативное априорное распределение вида Beta(30, 70) – то же, что использовать неинформативное априорное распределение Beta(1, 1) , добавив к данным 69 нулей и 29 единиц. По сути, мы просто добавили 98 выдуманных наблюдений к экспериментальным данным.

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

size = 385
size_ = size + 98

res = defaultdict(list)
for _ in range(10000):
    sum_a, sum_b = np.random.binomial(size, p1, 2) + 29
    sum_b_effect = np.random.binomial(size, p2) + 29
    res['freq_aa'].append(frequentist_test(sum_a, size_, sum_b, size_)[1])
    res['freq_ab'].append(frequentist_test(sum_a, size_, sum_b_effect, size_)[1])    
show_bernulli_ci(res)

freq_aa: pe=0.031, ci=[0.027, 0.034]
freq_ab: pe=0.860, ci=[0.853, 0.867]

Вероятности ошибок примерно равны вероятностям ошибок, полученным в байесовском подходе с априорными распределениями из Beta(30, 70) .

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

alpha = 0.031
p1 = 0.3
effect = 0.1
p2 = p1 + effect
size = 385

res = defaultdict(list)
for _ in range(10000):
    sum_a, sum_b = np.random.binomial(size, p1, 2)
    sum_b_effect = np.random.binomial(size, p2)
    res['freq_aa'].append(frequentist_test(sum_a, size, sum_b, size, alpha)[1])
    res['freq_ab'].append(frequentist_test(sum_a, size, sum_b_effect, size, alpha)[1])
show_bernulli_ci(res)

freq_aa: pe=0.030, ci=[0.027, 0.034]
freq_ab: pe=0.846, ci=[0.839, 0.853]

«Можно остановиться раньше»

Говорят, байесовский подход позволяет остановиться раньше запланированного. Проверим, как изменятся ошибки в байесовском и частотном подходах, если остановиться раньше. Уменьшим размер групп с 385 до 50.

p1 = 0.3                   # исходная конверсия
effect = 0.1               # эффект
p2 = p1 + effect           # конверсия с эффектом
prior_a = [0.5, 0.5]
prior_b = [0.5, 0.5]
size = 50

res = defaultdict(list)
for _ in range(10000):
    sum_a, sum_b = np.random.binomial(size, p1, 2)
    sum_b_effect = np.random.binomial(size, p2)
    res['freq_aa'].append(frequentist_test(sum_a, size, sum_b, size)[1])
    res['freq_ab'].append(frequentist_test(sum_a, size, sum_b_effect, size)[1])
    params = sum_a, size, sum_b, size, prior_a, prior_b
    res['bayes_aa'].append(bayesian_test(*params)[1])
    params = sum_a, size, sum_b_effect, size, prior_a, prior_b
    res['bayes_ab'].append(bayesian_test(*params)[1])
show_bernulli_ci(res)

freq_aa: pe=0.053, ci=[0.049, 0.058]
freq_ab: pe=0.291, ci=[0.282, 0.300]
bayes_aa: pe=0.052, ci=[0.048, 0.057]
bayes_ab: pe=0.289, ci=[0.280, 0.298]

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

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

«Можем проверять в любой момент»

Говорят, байесовский подход позволяет проверять гипотезу во время эксперимента. То есть для него нет проблемы подглядывания (peeking problem). Проверим это. Будем оценивать значимость отличий для каждой новой пары значений, начиная с размера групп, равного 30. Если значимые отличия будут найдены, то останавливаемся и говорим, что значимые отличия есть. Проведём много синтетических АА экспериментов и оценим вероятность ошибки первого рода. Будет ли такой подход на основе байесовского тестирования контролировать вероятность ошибки первого рода? Ответ под спойлером.

Скрытый текст

def bayesian_seq_test(a, b, prior_a, prior_b):
    sum_a, len_a = prior_a
    sum_b, len_b = prior_b
    for a_value, b_value in zip(a, b):
        len_a += 1
        len_b += 1
        sum_a += a_value
        sum_b += b_value
        if len_a < 30:
            continue
        _, decision = bayesian_test(
            sum_a, len_a, sum_b, len_b, prior_a, prior_b
        )
        if decision:
            return 1
    return 0

def frequentist_seq_test(a, b):
    sum_a, len_a = 0, 0
    sum_b, len_b = 0, 0
    for a_value, b_value in zip(a, b):
        len_a += 1
        len_b += 1
        sum_a += a_value
        sum_b += b_value
        if len_a < 30:
            continue
        _, decision = frequentist_test(
            sum_a, len_a, sum_b, len_b
        )
        if decision:
            return 1
    return 0

p = 0.3
sample_size = 100
prior_a = [0.5, 0.5]
prior_b = [0.5, 0.5]

res = defaultdict(list)
for _ in range(100):
    a, b = np.random.binomial(1, p, (2, sample_size))
    res['freq_aa'].append(frequentist_seq_test(a, b))
    res['bayes_aa'].append(bayesian_seq_test(a, b, prior_a, prior_b))
show_bernulli_ci(res)

freq_aa: pe=0.160, ci=[0.088, 0.232]
bayes_aa: pe=0.160, ci=[0.088, 0.232]

Вероятность ошибки первого рода не контролируется. Байесовский подход не помогает справиться с peeking problem.

Множественное тестирование

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

Итоги

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

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

Как выбрать априорное распределение и вычислить апостериорное распределение для данных из смеси различных распределений?
Как определить размер групп?
Как применить байесовский подход для метрик отношения?
Как применять способы повышения чувствительности типа CUPED и стратификации?

Хабы: