Сегодня мы поговорим о корректности статистических критериев в контексте A/B-тестирования. Мы научимся проверять, является ли критерий правильным или нет. Мы рассмотрим пример, в котором критерий Стьюдента не работает.

Правильный статистический критерий

В A/B-тестировании при проверке гипотез по статистическим критериям может быть допущена одна из двух ошибок:

  • Ошибка I рода — отклонение нулевой гипотезы, когда она действительно верна. То есть сказать, что есть эффект, когда на самом деле его нет.
  • Ошибка II рода — невозможность отвергнуть нулевую гипотезу, когда она на самом деле ложна. То есть сказать, что эффекта нет, когда на самом деле он есть.

Невозможно избежать всех ошибок. Для получения результатов, которые на 100% надежны, требуется бесконечное количество данных. На практике получить такой объем данных сложно. Если невозможно полностью избежать ошибок, желательно делать ошибки как можно реже и контролировать вероятность ошибок.

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

Предположим, мы решили, что допустимые вероятности ошибок первого и второго рода равны 0,1 и 0,2 соответственно. Будем называть статистический критерий корректным, если его вероятности ошибок I и II рода равны допустимым вероятностям ошибок I и II рода соответственно.

Как нам создать критерий, в котором вероятности ошибок равны допустимым вероятностям ошибок?

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

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

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

где α и β — допустимые вероятности ошибок I и II рода,
ε — ожидаемый эффект (насколько среднее значение изменится), а σ — стандартное отклонение случайной величины в контрольной и экспериментальной группах соответственно.

Проверка на правильность

Предположим, мы работаем в интернет-магазине с доставкой. Мы хотим исследовать, как новый алгоритм ранжирования товаров на сайте влияет на средний доход на одного покупателя в неделю. Продолжительность эксперимента — одна неделя. Ожидаемый эффект +100 руб. Допустимая вероятность ошибки первого рода равна 0,1, а ошибки второго рода — 0,2.

Рассчитаем необходимый размер группы по формуле:

import numpy as np
from scipy import stats

alpha = 0.1                     # acceptable probability of a type I error
beta = 0.2                      # acceptable probability of a type II error
mu_control = 2500               # average revenue per user in the control group.
effect = 100                    # expected effect size
mu_pilot = mu_control + effect  # average revenue per user in the experimental group
std = 800                       # standard deviation

# historical data for 10000 of clients
values = np.random.normal(mu_control, std, 10000)

def estimate_sample_size(effect, std, alpha, beta):
    """Estimation of the required group size"""
    t_alpha = stats.norm.ppf(1 - alpha / 2, loc=0, scale=1)
    t_beta = stats.norm.ppf(1 - beta, loc=0, scale=1)
    var = 2 * std ** 2
    sample_size = int((t_alpha + t_beta) ** 2 * var / (effect ** 2))
    return sample_size

estimated_std = np.std(values)
sample_size = estimate_sample_size(effect, estimated_std, alpha, beta)
print(f'Estimation of the required group size = {sample_size}')

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

Поскольку это было в прошлом, мы знаем, какие покупки совершали пользователи, можем рассчитать метрики и оценить значимость различий. Кроме того, мы знаем, что фактически никакого эффекта в то время не было, так как эксперимент фактически не был запущен. Если обнаруживались существенные различия, мы совершали ошибку I рода. В противном случае мы получили правильный результат.

Далее нам нужно повторить эту процедуру с гипотетическим запуском эксперимента в прошлом на разных группах и временных интервалах много раз, например, 1000.

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

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

Давайте посмотрим, как оценить вероятность ошибок в коде. Используя численные синтетические эксперименты A/A и A/B, мы оценим вероятности ошибок и построим доверительные интервалы.

def run_synthetic_experiments(values, sample_size, effect=0, n_iter=10000):
    """We conduct synthetic experiments and return a list of p-values"""
    pvalues = []
    for _ in range(n_iter):
        a, b = np.random.choice(values, size=(2, sample_size,), replace=False)
        b += effect
        pvalue = stats.ttest_ind(a, b).pvalue
        pvalues.append(pvalue)
    return np.array(pvalues)

def print_estimated_errors(pvalues_aa, pvalues_ab, alpha):
    """Estimating error probabilities"""
    estimated_first_type_error = np.mean(pvalues_aa < alpha)
    estimated_second_type_error = np.mean(pvalues_ab >= alpha)
    ci_first = estimate_ci_bernoulli(estimated_first_type_error, len(pvalues_aa))
    ci_second = estimate_ci_bernoulli(estimated_second_type_error, len(pvalues_ab))
    print(f'estimate of Type I error probability = {estimated_first_type_error:0.4f}')
    print(f'confidence interval = [{ci_first[0]:0.4f}, {ci_first[1]:0.4f}]')
    print(f'estimate of Type II error probability = {estimated_second_type_error:0.4f}')
    print(f'confidence interval = [{ci_second[0]:0.4f}, {ci_second[1]:0.4f}]')

def estimate_ci_bernoulli(p, n, alpha=0.05):
    """Confidence interval for a Bernoulli random variable"""
    t = stats.norm.ppf(1 - alpha / 2, loc=0, scale=1)
    std_n = np.sqrt(p * (1 - p) / n)
    return p - t * std_n, p + t * std_n

pvalues_aa = run_synthetic_experiments(values, sample_size, effect=0)
pvalues_ab = run_synthetic_experiments(values, sample_size, effect=effect)
print_estimated_errors(pvalues_aa, pvalues_ab, alpha)

Оценки вероятностей ошибок примерно равны 0,1 и 0,2, как и должно быть. Все верно, тест Стьюдента на этих данных работает корректно.

Распределение p-значений

Выше мы рассмотрели случай, когда тест контролирует вероятность ошибки I рода на фиксированном уровне значимости. Если мы решим изменить уровень значимости с 0,1 на 0,01, будет ли тест контролировать вероятность ошибки I рода? Было бы хорошо, если бы тест контролировал вероятность ошибки первого рода на любом заданном уровне значимости. Формально это можно записать так:

Для любого α∈[0,1] верно, что P(pvalue‹α∣H0​)=α.

Обратите внимание, что левая часть уравнения представляет собой выражение для функции распределения p-значения. Из уравнения следует, что функция распределения p-значения в точке X равна X для любого X от 0 до 1. Эта функция распределения является функцией распределения равномерного распределения от 0 до 1. Мы только что показали, что статистическая Критерий контролирует вероятность ошибки первого рода на заданном уровне для любого уровня значимости тогда и только тогда, когда p-значение равномерно распределено от 0 до 1 при нулевой гипотезе.

При нулевой гипотезе p-значение должно быть распределено равномерно. И как должно быть распределено p-значение согласно альтернативной гипотезе? Из условия вероятности ошибки II рода P(pvalue≥α∣H1)=β следует, что P(pvalue‹α∣H1)=1-β.

Таким образом, график функции распределения p-значений при альтернативной гипотезе должен проходить через точку [α, 1-β], где α и β — допустимые вероятности ошибки для конкретного эксперимента.

Давайте проверим, как распределяется p-значение в численном эксперименте. Мы построим эмпирические функции распределения p-значения:

import matplotlib.pyplot as plt

def plot_pvalue_distribution(pvalues_aa, pvalues_ab, alpha, beta):
    """p-value distribution chart"""
    estimated_first_type_error = np.mean(pvalues_aa < alpha)
    estimated_second_type_error = np.mean(pvalues_ab >= alpha)
    y_one = estimated_first_type_error
    y_two = 1 - estimated_second_type_error
    X = np.linspace(0, 1, 1000)
    Y_aa = [np.mean(pvalues_aa < x) for x in X]
    Y_ab = [np.mean(pvalues_ab < x) for x in X]

    plt.plot(X, Y_aa, label='A/A')
    plt.plot(X, Y_ab, label='A/B')
    plt.plot([alpha, alpha], [0, 1], '--k', alpha=0.8)
    plt.plot([0, alpha], [y_one, y_one], '--k', alpha=0.8)
    plt.plot([0, alpha], [y_two, y_two], '--k', alpha=0.8)
    plt.plot([0, 1], [0, 1], '--k', alpha=0.8)

    plt.xlabel('p-value', size=12)
    plt.legend(fontsize=12)
    plt.grid()
    plt.show()

plot_pvalue_distribution(pvalues_aa, pvalues_ab, alpha, beta)

Значение Р для синтетических А/А-тестов оказалось равномерно распределенным от 0 до 1, а для синтетических А/В-тестов оно проходит через точку [α,1-β].

Кроме оценок распределения на график нанесены еще четыре пунктирные линии:

  • диагональная линия от точки [0, 0] до точки [1, 1] представляет кумулятивную функцию распределения равномерного распределения на интервале от 0 до 1, что позволяет визуально оценить равномерность распределения p-значения;
  • вертикальная линия с x=α представляет собой пороговое значение p для принятия решения, отклонять нулевую гипотезу или нет. Проекция на ось Y точки пересечения вертикальной линии с функцией распределения p-значения для тестов A/A представляет вероятность ошибки I рода. Проекция точки пересечения вертикальной линии с функцией распределения p-значения для тестов A/B представляет мощность теста (мощность = 1 — β).
  • две горизонтальные линии представляют собой проекции на ось Y точек пересечения вертикальной линии с функцией распределения p-значения для тестов A/A и A/B.

График с оценками распределения p-значений для синтетических тестов A/A и A/B позволяет проверить корректность теста для любого уровня значимости.

Неверный критерий

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

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

def run_synthetic_experiments_two(values, sample_size, effect=0, n_iter=10000):
    """2-weeks synthetic experiment"""
    pvalues = []
    for _ in range(n_iter):
        a, b = np.random.choice(values, size=(2, sample_size,), replace=False)
        b += effect
        # дублируем данные
        a = np.hstack((a, a,))
        b = np.hstack((b, b,))
        pvalue = stats.ttest_ind(a, b).pvalue
        pvalues.append(pvalue)
    return np.array(pvalues)

pvalues_aa = run_synthetic_experiments_two(values, sample_size)
pvalues_ab = run_synthetic_experiments_two(values, sample_size, effect=effect)
print_estimated_errors(pvalues_aa, pvalues_ab, alpha)
plot_pvalue_distribution(pvalues_aa, pvalues_ab, alpha, beta)

Мы получили оценку вероятности ошибки I рода около 0,25, что намного выше уровня значимости 0,1. На графике видно, что распределение p-значений для синтетических A/A-тестов неравномерно и отклоняется от диагонали. В этом примере t-критерий Стьюдента некорректен, поскольку данные зависимы (зависимы затраты на покупки одного человека). Если бы мы сразу не осознали зависимость данных, оценка вероятностей ошибок помогла бы нам понять, что такой тест некорректен.

В итоге:

  • Правильный статистический тест — это тест, в котором вероятности ошибок первого и второго рода равны допустимым вероятностям ошибок первого и второго рода соответственно.
  • Для контроля вероятности ошибки I рода для любого уровня значимости необходимо и достаточно, чтобы p-значение было равномерно распределено от 0 до 1, когда нулевая гипотеза верна.