Проблема подглядывания и последовательное А-Б тестирование

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

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

Проблема подглядывания

Допустим, изучаем химические реакции. Нас интересует, какие реагенты выделяют больше газа. Чтобы сравнить два реагента проводим эксперименты поочерёдно то с одним реагентом, то с другим. Один эксперимент занимает 3 часа. Реагенты дорогие и их количество ограниченно. Сможем провести максимум по 128 измерений для каждого реагента — это займёт целый месяц.

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

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

• данные из нормального распределения;

• проверяется односторонняя гипотеза о равенстве средних тестом Стьюдента;

• допустимые вероятности ошибок первого и второго рода равны 0.05 и 0.1 соответственно;

• ожидаемый эффект подобран так, чтобы размер групп был равен 128;

• вероятности ошибок оцениваем для разного количества подглядываний от 1 до 64 с шагом по степеням двойки.

Построим график оценки вероятности ошибки первого рода в зависимости от количества подглядываний.

from collections import defaultdict
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats

ALPHA, BETA = 0.05, 0.1    # вероятности ошибок
MEAN = 0                   # исходное среднее
STD = 1                    # исходное стандартное отклонение
EFFECT = 0.366             # эффект

sum_var = 2 * STD ** 2
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 ttest_ind(a, b, **kwargs) -> float:
    """Односторонний ttest, возвращает pvalue."""
    return stats.ttest_ind(a, b, alternative='less', **kwargs).pvalue

# количество подглядываний
array_peekings = 2 ** np.arange(7)
# словарь со значениями pvalue
peekings_2_pvalues_aa = defaultdict(list)
peekings_2_pvalues_ab = defaultdict(list)

for _ in range(1000):
    a, b = np.random.normal(MEAN, STD, (2, 1000, SIZE))
    for peekings in array_peekings:
        step = SIZE // peekings
        list_pvalues_aa = []
        list_pvalues_ab = []
        for size_ in range(step, SIZE + 1, step):
            a_ = a[:, :size_]
            b_ = b[:, :size_]
            list_pvalues_aa.append(ttest_ind(a_, b_, axis=1))
            list_pvalues_ab.append(ttest_ind(a_, b_+EFFECT, axis=1))
        peekings_2_pvalues_aa[peekings] += list(zip(*list_pvalues_aa))
        peekings_2_pvalues_ab[peekings] += list(zip(*list_pvalues_ab))


# оцениваем вероятности ошибок первого рода
list_errors = []
for peekings in array_peekings:
    pvalues_aa = peekings_2_pvalues_aa[peekings]
    errors = np.mean(np.min(pvalues_aa, axis=1) < ALPHA)
    list_errors.append(errors)

# строим график
X = array_peekings
Y = list_errors
plt.plot(X, Y, '-')
plt.hlines(Y, 0, X, colors='k', linestyles='dashed', alpha=0.5)
plt.vlines(X, 0, Y, colors='k', linestyles='dashed', alpha=0.5)
plt.plot(X, Y, 'ok', markersize=5)
for x, y in zip(X, Y):
    t_ = plt.text(x*0.95+2, y*0.94-0.003, f'{y:0.2f}')
    t_.set_bbox(dict(facecolor='white', alpha=1, edgecolor='none'))
plt.ylim([0, 0.31])
plt.xlim([0, 68])
plt.grid()
plt.xlabel('Количество подглядываний')
plt.title('Вероятность совершить ошибку первого рода')
plt.show()

Из графика видно, что при увеличении количества подглядываний вероятность ошибки первого рода увеличивается. Например, если подглядываем раз в 2 дня (16 раз за эксперимент), то вероятность ошибки первого рода равна примерно 0.2 — в четыре раза выше допустимой вероятности ошибки первого рода. Подглядывание в лоб делает проверку эксперимента некорректной.

При большом количестве проверок вероятность допустить хотя бы одну ошибку растёт. Аналогичная проблема есть во множественном тестировании. Отличие подглядывания от множественного тестирования заключается в зависимости данных. Для вычисления p-value при очередном подглядывании используются данные, которые использовались для прошлого подглядывания. 

Корректируем уровень значимости

Для контроля вероятности ошибки первого рода во множественном тестировании есть метод Бонферрони. Он предлагает сравнивать p-value не с уровнем значимости, а с, где — количество экспериментов. Метод Бонферрони для последовательного тестирования работает некорректно, так как данные при подглядываниях зависимые.

Сконструируем свой критерий аналогичный критерию Бонферрони. Зафиксируем количество подглядываний и подберём значение, с которым будем сравнивать p-value так, чтобы вероятность ошибки первого рода была равна. Оценим численно для разного числа подглядываний. Используя полученные значения, оценим вероятности ошибок.

res = defaultdict(list)
for peekings in array_peekings:
    res['peekings'].append(peekings)
    min_pvalues_aa = np.min(peekings_2_pvalues_aa[peekings], axis=1)
    quantile = np.quantile(min_pvalues_aa, ALPHA)
    res['alpha*'].append(quantile)
    res['I type error'].append(np.mean(min_pvalues_aa < quantile))
    min_pvalues_ab = np.min(peekings_2_pvalues_ab[peekings], axis=1)
    res['II type error'].append(np.mean(min_pvalues_ab > quantile))
pd.DataFrame(res).set_index('peekings').round(3)

Несколько наблюдений по полученным результатам:

• Скорректированные уровни значимости отличаются от значений для метода Бонферрони. Например, для двух подглядываний у нас получилось 0.03, а по методу Бонферрони было бы 0.025.

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

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

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

Критерий Вальда

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

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

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

• Если значение меньше левой границы, то отклоняем альтернативную гипотезу (говорим, что эффекта нет) и останавливаем эксперимент.

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

• Иначе продолжаем эксперимент.

Пример реализации критерия Вальда:

def test_sequential_wald_norm(a, b, std, effect, alpha, beta):
    """Применяем критерий Вальда для нормального распределения.

    a, b - данные контрольной и экспериментальной групп;
    std - стандартное отклонение;
    effect - размер эффекта;
    alpha, beta - допустимые вероятности ошибок первого и второго рода.

    return (decision, value) - решение и значение статистики.
        Значения decision:
            - 0 - эффекта нет
            - 0.5 - границы не пересечены
            - 1 - эффект есть
    """
    coef = 2 * std ** 2 / effect
    lower_bound = coef * np.log(beta / (1 - alpha))
    upper_bound = coef * np.log((1 - beta) / alpha)
    value = len(a) * (np.mean(b) - np.mean(a) - effect / 2)
    if value < lower_bound:
        return 0, value
    elif value > upper_bound:
        return 1, value
    return 0.5, value

Построим визуализацию статистик критерия Вальда в течение эксперимента. Также отобразим границы критерия Вальда и границу теста Стьюдента. Границу теста Стьюдента получим на основе распределения Стьюдента:

Откуда выразим:

size_ = 230

def run_experiment(a, b):
    """Имитируем постепенное получение данных."""
    values = []
    for n in range(1, size_ + 1):
        decision, value = test_sequential_wald_norm(
            a[:n], b[:n], STD, EFFECT, ALPHA, BETA
        )
        values.append(value)
        if decision in [0, 1]:
            break
    return decision, values

curves_aa = []
curves_ab = []
for _ in range(10):
    a_one, a_two = np.random.normal(MEAN, STD, (2, size_))
    b = np.random.normal(MEAN + EFFECT, STD, size_)
    _, values_aa = run_experiment(a_one, a_two)
    _, values_ab = run_experiment(a_one, b)
    curves_aa.append(values_aa)
    curves_ab.append(values_ab)

plt.rcParams['figure.figsize'] = (10, 6)
# строим кривые статистики Вальда
curves = [
    [curves_aa, 'b', 'AA test'],
    [curves_ab, 'r', 'AB test'],
]
for list_curve, color, label in curves:
    for idx, curve in enumerate(list_curve):
        plt.plot(
            np.arange(len(curve)),
            curve,
            color,
            label=label if idx==0 else '',
            alpha=0.6
        )
        plt.scatter([len(curve) - 1], curve[-1], color=color)

# отображаем границы критерия Вальда
coef = 2 * STD ** 2 / EFFECT
lower_bound = coef * np.log(BETA / (1 - ALPHA))
upper_bound = coef * np.log((1 - BETA) / ALPHA)
plt.hlines(upper_bound, 0, size_, color='r', linestyle='--')
plt.hlines(lower_bound, 0, size_, color='b', linestyle='--')
plt.fill_between(
    np.arange(size_ + 1), upper_bound, upper_bound + 5,
    color='r', alpha=0.1
)
plt.fill_between(
    np.arange(size_ + 1), lower_bound, lower_bound - 5,
    color='b', alpha=0.1
)
plt.vlines(
    SIZE, lower_bound - 5, upper_bound + 5,
    linestyle='--', label='sample size'
)
plt.text(150, upper_bound+2, 'Reject $H_0$', color='r', size=16)
plt.text(150, lower_bound-3, 'Reject $H_1$', color='b', size=16)

# Строим границы критерия Стьюдента
array_n = np.arange(2, size_ + 1)
student_bounds = [
    ((2/n)**0.5 * STD * stats.t.ppf(1-ALPHA, df=2*(n-1)) - EFFECT/2) * n
    for n in array_n
]
plt.plot(array_n, student_bounds, '-k', label='student bounds', alpha=0.8)

plt.xlabel('Количество данных')
plt.title('Траектории статистики критерия Вальда')
plt.legend()
plt.show()

Рассмотрим полученные графики:

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

• Красные и синие кривые — траектории значений статистики для синтетических АА и АВ экспериментов. В двух синтетических АВ экспериментах были допущены ошибки: красные траектории пересекли синюю границу. В остальных экспериментах получены верные результаты.

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

• Чёрная линия — приближённая граница критерия Стьюдента. Если бы применяли его при подглядывании, то допустили бы ошибку первого рода: одна из синих кривых пересекает чёрную линию.

Проверим корректность работы критерия Вальда и посмотрим распределение продолжительности экспериментов. Проведёт 100 000 синтетических АА и АВ тестов. Для каждого запуска будем сохранять вердикт критерия Вальда, использованный объём данных и p-value теста Стьюдента для размера групп, равного необходимому размеру групп.

def pdf_a(x):
    """Функция плотности разницы средних при верности нулевой гипотезы."""
    return stats.norm.pdf(x, 0, np.sqrt(2) * STD)

def pdf_b(x):
    """Функция плотности разницы средних при верности альтернативной гипотезы."""
    return stats.norm.pdf(x, EFFECT, np.sqrt(2) * STD)

def test_sequential_wald(a, b, pdf_a, pdf_b, alpha, beta):
    """Проводим последовательное тестирование критерием Вальда.

    a, b - данные контрольной и экспериментальной групп
    pdf_a, pdf_b - функции плотности распределения при нулевой и альтернативной гипотезах
    alpha, beta - допустимые вероятности ошибок первого и второго рода

    return:
        - decision: 0 - незначимые отличия, 1 - значимые отличия, 0.5 - границы не пересечены
        - length - кол-во объектов при принятии решения
    """
    size = len(a)
    lower_bound = np.log(beta / (1 - alpha))
    upper_bound = np.log((1 - beta) / alpha)

    deltas = b - a
    pdf_a_values = pdf_a(deltas)
    pdf_b_values = pdf_b(deltas)
    z = np.cumsum(np.log(pdf_b_values / pdf_a_values))
    
    indexes_lower = np.arange(size)[z < lower_bound]
    indexes_upper = np.arange(size)[z > upper_bound]
    first_index_lower = indexes_lower[0] if len(indexes_lower) > 0 else size + 1
    first_index_upper = indexes_upper[0] if len(indexes_upper) > 0 else size + 1

    if first_index_lower < first_index_upper:
        return 0, first_index_lower + 1
    elif first_index_lower > first_index_upper:
        return 1, first_index_upper + 1
    else:
        return 0.5, size

size_ = 500
dict_res = defaultdict(list)
for _ in range(100000):
    a, a_, b_ = np.random.normal(MEAN, STD, (3, size_))
    b_ += EFFECT
    for key, b in [('AA', a_), ('AB', b_)]:
        wald_decision, length = test_sequential_wald(
            a, b, pdf_a, pdf_b, ALPHA, BETA
        )
        ttest_pvalue = ttest_ind(a[:SIZE], b[:SIZE])
        dict_res[key].append((wald_decision, length, ttest_pvalue))


import seaborn as sns
plt.rcParams['figure.figsize'] = (6, 4)
lengths = sum([[x[1] for x in v] for v in dict_res.values()], start=[])
part_more_size = np.mean(np.array(lengths) > SIZE)
sns.kdeplot(lengths)
plt.vlines(SIZE, 0, 0.012, linestyle='--', color='r', label='sample size')
plt.text(170, 0.0015, f'Доля экспериментов с размером группnбольше sample size равна {part_more_size:0.2f}')

plt.grid()
plt.title('Распределение размера групп остановки критерия Вальда')
plt.xlabel('Размер групп')
plt.legend()
plt.show()


agg_res = defaultdict(list)
for key, values in dict_res.items():
    agg_res['method'].append(f'part error {key}')
    wald_decisions, lengths, ttest_pvalues = np.array(values).T
    agg_res['ttest'].append(np.mean(ttest_pvalues < ALPHA))
    agg_res['wald'].append(np.mean(wald_decisions))
    wald_ttest_decisions = wald_decisions.copy()
    wald_ttest_decisions[lengths > SIZE] = ttest_pvalues[lengths > SIZE] < ALPHA
    agg_res['wald+ttest'].append(np.mean(wald_ttest_decisions))
lengths = [x[1] for x in dict_res['AA']] + [x[1] for x in dict_res['AB']]
agg_res['method'].append('mean length')
agg_res['ttest'].append(SIZE)
agg_res['wald'].append(np.mean(lengths))
agg_res['wald+ttest'].append(np.mean(np.clip(lengths, 0, SIZE)))
agg_res['method'].append('max length')
agg_res['ttest'].append(SIZE)
agg_res['wald'].append(np.max(lengths))
agg_res['wald+ttest'].append(SIZE)

df_agg = pd.DataFrame(agg_res).set_index('method').T.round(3)
df_agg['part error AB'] = 1 - df_agg['part error AB']
df_agg['mean length'] = df_agg['mean length'].round(1)
df_agg

Посмотрим на график распределения продолжительности экспериментов при использовании критерия Вальда. Порядка 88% экспериментов остановились раньше достижения sample size. Чаще всего эксперименты останавливались на размере групп, равном 34 — почти в 4 раза меньше sample size. Иногда случаются долгие эксперименты, наш максимум равен 679.

Далее рассмотрим три подхода к принятию решений:

1. Тест Стьюдента для размера групп, равного 124.

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

3. Применяем критерий Вальда. Если критерий не останавливается при размере групп, равном 124, то применяем тест Стьюдента и сравниваем его p-value с уровнем значимости 0.05. Такой способ может быть полезен, если нет возможности получить дополнительные данные.

Результаты оценок вероятностей ошибок разных методов:

• Тест Стьюдента работает корректно, имеет предсказуемую фиксированную продолжительность, но не допускает подглядывания.

• Тест Вальда работает корректно, в среднем требует меньше данных, но иногда данных нужно больше, чем тесту Стьюдента.

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

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

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

Если хотите использовать какие-либо техники подглядывания для фиксированного объёма данных, то придётся мириться с увеличением вероятностей ошибок.

Ограничения

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

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

Итоги

Основные тезисы статьи:

• Мы хотим подглядывать на данные во время эксперимента, чтобы быстрее получить результат. Если делать это бездумно, то растут вероятности ошибок — это явление называется проблемой подглядывания (peeking problem).

• Есть методы, которые позволяют подсматривать фиксированное количество раз, контролировать вероятность ошибки первого рода, но мощность теста уменьшается. Реализовали численный аналог метода Покока.

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

• Последовательное тестирование сложнее в реализации.