nnazarov 11 сен в 12:52

А/Б тестирование: множественная проверка гипотез

10 мин

4.3K

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

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

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

Много гипотез

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

$\mathbb{P}(FP>0) = 1 - (1 - \alpha)^m$

где — количество ложноположительных результатов (ошибок первого рода), а — количество гипотез.

Ниже приведён график этой функции для $\alpha= 0.05$ . При проверке 100 гипотез вероятность получить ошибку первого рода близка к единице.

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

Если среди 100 гипотез будет только одна гипотеза с эффектом, то с вероятностью 0.24 среди 99 гипотез без эффекта найдётся гипотеза, у которой точечная оценка эффекта будет больше точечной оценки эффекта гипотезы, у которой эффект на самом деле был.

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

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

Метрики

Перед тем как перейти к поправкам на множественное тестирование, определим метрики, которые хотим контролировать. Введём некоторые обозначения:

m — общее число гипотез;
TN — число истинно отрицательных результатов;
FN — число ложноотрицательных результатов;
FP — число ложноположительных результатов;
TP — число истинно положительных результатов.

	верна	неверна	Всего
Не отвергаем	TN	FN	TN+FN
Отвергаем	FP	TP	FP+TP
Всего	TN+FP	FN+TP	m

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

FWER (Family-wise error rate) — вероятность того, что хотя бы один из тестов даёт ложноположительный результат при истинности нулевых гипотез.

$FWER = \mathbb{P}(FP>0) = 1 - (1 - \alpha)^m$

FDR (False discovery rate) – средняя доля ложноположительных результатов среди положительных результатов.

$FDR = \mathbb{E} \left( \frac{FP}{\max(FP+TP, ~ 1)} \right)$

При прочих равных алгоритмы принятия решений, которые контролируют FWER на заданном уровне значимости, допускают меньше ложноположительных результатов, чем алгоритмы, контролирующие FDR. Допустим, в 20 гипотезах из 100 эффект на самом деле есть. Для контроля FDR на уровне значимости 0.05 достаточно выбирать 20 гипотез, допуская один ложноположительный результат. При этом FWER будет равен 1. Алгоритм, контролирующий FWER, будет выбирать меньше гипотез, что приведёт к уменьшению ложноположительных результатов, но, в то же время, и к увеличению ложноотрицательных результатов.

Также нам понадобится аналог мощности для множественного тестирования. Будем использовать среднюю долю истинно положительных результатов среди гипотез с эффектом.

$POWER = \mathbb{E} \left( \frac{TP}{\max(TP+FN, ~ 1)} \right)$

Поправки

Существует много вариантов поправок на множественную проверку гипотез. Рассмотрим три из них: Бонферрони, Холма и Бенджамини-Хохбрега.

Метод Бонферрони предлагает использовать для проверки каждой гипотезы уровень значимости, равный $\alpha / m$ , где — количество гипотез. Это позволяет контролировать FWER на уровне, не превосходящем $\alpha$ . Метод Бонферрони прост в реализации, но уступает в мощности другим поправкам.

Метод Холма также контролирует FWER на уровне, не превосходящем $\alpha$ . Опишем алгоритм его работы:

Упорядочим значения p-value по возрастанию: $p_{(1)}, \ldots, p_{(m)}$ . Обозначим нулевую гипотезу, соответствующую значению $p_{(i)}$ , как $H_{(i)}$ .
Определим уровни значимости $\alpha_i = \alpha / (m - i + 1)$ .
Последовательно перебираем значения p-value из ряда $p_{(1)}, \ldots, p_{(m)}$ . Если $p_{(i)} < \alpha_i$ , то отклоняем $H_{(i)}$ и переходим к следующему значению. Иначе говорим, что данные гипотезам $H_{(i)}, \ldots, H_{(m)}$ не противоречат.

Метод Бенджамини-Хохберга контролирует FDR на уровне, не превосходящем $\alpha$ . Данный метод работает так же, как и метод Холма, только в качестве уровней значимости используются значения $\alpha_i = \alpha * i / m$ .

Ниже приведены реализации поправок. Можно использовать готовые реализации, например, из функции statsmodels.stats.multitest.multipletests.

Код

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


def method_without_correct(pvalues, alpha=0.05):
    """Проверяет значимость отличий без поправок.

    :param pvalues (np.array): массив pvalue
    :param alpha (float): уровень значимости

    :return res (np.array): массив из нулей и единиц:
        0 - значимых отличий нет, 1 - значимые отличия есть.
    """
    res = (pvalues <= alpha).astype(int)
    return res

def method_bonferroni(pvalues, alpha=0.05):
    """Применяет метод Бонферрони для проверки значимости изменений.

    :param pvalues (np.array): массив pvalue
    :param alpha (float): уровень значимости

    :return res (np.array): массив из нулей и единиц:
        0 - значимых отличий нет, 1 - значимые отличия есть.
    """
    m = len(pvalues)
    res = (pvalues <= alpha / m).astype(int)
    return res

def method_holm(pvalues, alpha=0.05):
    """Применяет метод Холма для проверки значимости изменений.
    
    :param pvalues (np.array): массив pvalue
    :param alpha (float): уровень значимости

    :return res (np.array): массив из нулей и единиц:
        0 - значимых отличий нет, 1 - значимые отличия есть.
    """
    m = len(pvalues)
    alphas = alpha / np.arange(m, 0, -1)
    sorted_pvalue_indexes = np.argsort(pvalues)
    res = np.zeros(m, dtype=int)
    for alpha_, pvalue_index in zip(alphas, sorted_pvalue_indexes):
        pvalue = pvalues[pvalue_index]
        if pvalue < alpha_:
            res[pvalue_index] = 1
        else:
            break
    return res

def method_benjamini_hochberg(pvalues, alpha=0.05):
    """Применяет метод Бенджамини-Хохберга для проверки значимости изменений.
    
    :param pvalues (np.array): массив pvalue
    :param alpha (float): уровень значимости

    :return res (np.array): массив из нулей и единиц:
        0 - значимых отличий нет, 1 - значимые отличия есть.
    """
    m = len(pvalues)
    alphas = alpha * np.arange(1, m+1) / m
    sorted_pvalue_indexes = np.argsort(pvalues)
    res = np.zeros(m, dtype=int)
    for alpha_, pvalue_index in zip(alphas, sorted_pvalue_indexes):
        pvalue = pvalues[pvalue_index]
        if pvalue <= alpha_:
            res[pvalue_index] = 1
        else:
            break
    return res

Численные эксперименты

Проведём численные эксперименты, чтобы посмотреть, как ведут себя поправки в разных условиях. Будем проверять 100 гипотез на равенство средних тестом Стьюдента. Уровень значимости равен 0.05, допустимая вероятность ошибки второго рода — 0.1, ожидаемый эффект — 275. Размер групп подобран по формуле для оценки необходимого размера групп при проверке одной гипотезы и равен 100. В экспериментах рассматривались ситуации, когда среди 100 гипотез было от 0 до 100 гипотез с эффектом. Результаты с оценками метрик для разных методов отображены на графиках ниже.

Код

def estimate_sample_size(effect, std, alpha, beta):
    """Оценка необходимого размер групп."""
    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(np.ceil((t_alpha + t_beta) ** 2 * var / (effect ** 2)))
    return sample_size

def run(dict_methods, mean, std, sample_size, n_exp, n_iter=1000):
    """Проводит симуляцию экспериментов.
    
    :param dict_methods: словарь с методами множественного тестирования
    :param mean: среднее значение генерируемых данных
    :param std: стандартное отклонение генерируемых данных
    :param sample_size: размер групп
    :param n_exp: количество проверяемых гипотез
    :param n_iter: количество итераций симуляции

    :return method_2_result: результаты применения методов к данным симуляции
    """
    list_n_exp_with_effect = [0, 1] + list(range(10, n_exp, 10)) + [n_exp-1, n_exp]
    method_2_result = {
        method_name: defaultdict(list)
        for method_name in dict_methods
    }
    for _ in range(n_iter):
        a_values, b_values = np.random.normal(mean, std, (2, n_exp, sample_size))
        prev_n_exp_with_effect = 0
        for n_exp_with_effect in list_n_exp_with_effect:
            b_values[prev_n_exp_with_effect: n_exp_with_effect] += effect
            prev_n_exp_with_effect = n_exp_with_effect
            pvalues = stats.ttest_ind(a_values, b_values, axis=1).pvalue
            for method_name, method in dict_methods.items():
                res = method(pvalues)
                method_2_result[method_name][n_exp_with_effect].append({
                    'FN': sum(res[:n_exp_with_effect] == 0),
                    'TP': sum(res[:n_exp_with_effect] == 1),
                    'FP': sum(res[n_exp_with_effect:] == 1),
                    'TN': sum(res[n_exp_with_effect:] == 0),
                })
    return method_2_result

def plot_metrics(method_2_result, n_exp, alpha=0.05, beta=0.1):
    """Рисует графики."""
    method_2_metrics = {method_name: defaultdict(list) for method_name in method_2_result}
    markers = ['o', 'v', '<', '>', '^', 's']
    for method_name in method_2_result:
        list_n_exp_with_effect = sorted(method_2_result[method_name].keys())
        for n_exp_with_effect in list_n_exp_with_effect:
            results = method_2_result[method_name][n_exp_with_effect]
            fwer = np.mean([x['FP'] > 0 for x in results])
            fdr = np.mean([x['FP'] / max(1, x['FP'] + x['TP']) for x in results])
            power = np.mean([x['TP'] / max(1, x['FN'] + x['TP']) for x in results])
            method_2_metrics[method_name]['FWER'].append(fwer)
            method_2_metrics[method_name]['FDR'].append(fdr)
            method_2_metrics[method_name]['POWER'].append(power)
    for metric_name in ['FWER', 'FDR', 'POWER']:
        min_value = 1
        for idx, (method_name, metrics) in enumerate(method_2_metrics.items()):
            marker = markers[idx % len(markers)]
            if metric_name in ['FWER', 'FDR']:
                X = list_n_exp_with_effect[:1] + list_n_exp_with_effect[2: -1]
                Y = metrics[metric_name][:1] + metrics[metric_name][2: -1]
            elif metric_name == 'POWER':
                X = list_n_exp_with_effect[1: -2] + list_n_exp_with_effect[-1:]
                Y = metrics[metric_name][1: -2] + metrics[metric_name][-1:]
            min_value = min(min_value, min(Y))
            plt.plot(X, Y, f'-{marker}', label=method_name)
        if metric_name in ['FWER', 'FDR']:
            plt.hlines(alpha, 0, n_exp, 'k', linestyles='--', label=f'alpha={alpha}')
        else:
            plt.hlines(1-beta, 0, n_exp, 'k', linestyles='--', label=f'power={1-beta}')
        plt.title(metric_name)
        plt.xlabel('Количество гипотез с эффектом')
        plt.legend()
        plt.grid()
        plt.xlim([0, n_exp])
        plt.ylim([0.8 if min_value > 0.8 else 0, 1])
        plt.show()

alpha = 0.05              # уровень значимости
beta = 0.1                # допустимая вероятность ошибки II рода
mean = 1000               # средняя выручка с пользователя в контрольной группе
std = 600                 # стандартное отклонение
effect = 276              # ожидаемый размер эффекта
sample_size = estimate_sample_size(effect, std, alpha, beta)
n_exp = 100               # количество проверяемых гипотез
dict_methods = {
    'without_correct': method_without_correct,
    'bonferroni': method_bonferroni,
    'holm': method_holm,
    'benjamini_hochberg': method_benjamini_hochberg,
}
method_2_result = run(dict_methods, mean, std, sample_size, n_exp)
plot_metrics(method_2_result, n_exp)

Графики FWER и FDR построены для точек от 0 до 99. В точке 100 нет гипотез без эффекта, поэтому значения FWER и FDR там равны нулю для всех методов.

Из графиков для FWER видно, что методы Бонферрони и Холма контролируют FWER на уровне, не превосходящем $\alpha$ , а метод Бенджамини-Хохберга — нет. График для FDR показывает, что метод Бенджамини-Хохберга контролируют FDR на уровне, не превосходящем $\alpha$ .

Графики мощности построены для точек от 1 до 100. В точке 0 нет гипотез с эффектом, значения мощности там равны нулю для всех методов.

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

Мощность

Чтобы контролировать мощность на заданном уровне, попробуем использовать уровень значимости, равный $\alpha / m$ в формуле для оценки необходимого размера групп.

Код

sample_size = estimate_sample_size(effect, std, alpha/n_exp, beta)
method_2_result = run(dict_methods, mean, std, sample_size, n_exp)
plot_metrics(method_2_result, n_exp)

Мощность увеличилась, минимальные значения мощности находятся в районе 0.9, как мы и хотели.

Мощность метода Бонферрони чуть меньше 0.9. Это связано с тем, что при оценке необходимого размера групп используется приближение нормальным распределением вместо распределения Стьюдента. Для $\alpha=0.05$ отличия могли быть незаметны, а при $\alpha=0.0005$ оценка размера групп отличается: 215 вместо 217.

Общий контроль

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

Какой выбрать размер групп, чтобы мощность сохранилась на желаемом уровне? Вспомним формулу, связывающую параметры критерия:

$\varepsilon^2 = \left[ \Phi^{-1} \left( 1- \alpha / 2 \right) + \Phi^{-1} \left( 1-\beta \right) \right]^2 \left(\frac{\sigma_A^2}{n_A} + \frac{\sigma_B^2}{n_B}\right)$

Проверяем гипотез c объединённой контрольной группой. Положим n_A=n m , n_B=n и $\sigma=\sigma_A=\sigma_B$ , тогда

$n = \frac{m+1}{m} \frac{\sigma^2}{\varepsilon^2} \left[ \Phi^{-1} \left( 1- \alpha / 2 \right) + \Phi^{-1} \left( 1-\beta \right) \right]^2$

При m=1 получаем (m+1)/m=2 , а при m=100 получаем . Ранее размер групп был равен 215, теперь будем использовать размер групп, равный $215/2*1.01 \approx 109$ . Проведём численный эксперимент:

Код

sample_size = 109
list_n_exp_with_effect = [0, 1] + list(range(10, n_exp, 10)) + [n_exp-1, n_exp]
method_2_result = {
    method_name: defaultdict(list)
    for method_name in dict_methods
}
for _ in range(1000):
    a_values, b_values = np.random.normal(mean, std, (2, n_exp, sample_size))
    a_values = [a_values.flatten()]
    prev_n_exp_with_effect = 0
    for n_exp_with_effect in list_n_exp_with_effect:
        b_values[prev_n_exp_with_effect: n_exp_with_effect] += effect
        prev_n_exp_with_effect = n_exp_with_effect
        pvalues = stats.ttest_ind(a_values, b_values, axis=1).pvalue
        for method_name, method in dict_methods.items():
            res = method(pvalues)
            method_2_result[method_name][n_exp_with_effect].append({
                'FN': sum(res[:n_exp_with_effect] == 0),
                'TP': sum(res[:n_exp_with_effect] == 1),
                'FP': sum(res[n_exp_with_effect:] == 1),
                'TN': sum(res[n_exp_with_effect:] == 0),
            })
plot_metrics(method_2_result, n_exp)

Качество работы методов не изменилось, а необходимый объём данных уменьшился почти в 2 раза. На самом деле, можно получить ещё больший выигрыш при оптимальном разбиении пользователей на группы. Как это сделать – читайте в нашей статье Определяем оптимальный размер групп при множественном А/Б тестировании.

Итоги

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

Теги:

Хабы: