Я не знал с чего начать статью, поэтому попросил ChatGPT сгенерировать анекдот:

Продуктовый аналитик говорит менеджеру: «Если не рассчитаем размер выборки перед экспериментом, это как пытаться испечь торт без рецепта!»
Менеджер отвечает: «Но иногда получаются вкусные эксперименты!»
Аналитик: «Да, но чаще всего — просто кекс с сюрпризом, который никто не хочет есть!»

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

Привет, Хабр! В процессе написания предыдущей статьи меня всё время не покидал вопрос: "А обязательно ли нам увеличивать размеры групп для наших экспериментов? Ведь гораздо выгоднее использовать меньшее количество пользователей продукта (или магазинов, курьеров и т.д.), сохраняя при этом репрезентативность выборок и сопоставимый уровень мощности теста". К тому же мы заранее понимаем, что для анализа проведенного эксперимента будет использоваться метод CUPED, который опирается на исторические данные о метрике. В связи с этим возникает вопрос: можем ли мы предварительно вычислить необходимый размер групп для A/B-тестирования с учетом применения этого метода анализа? Сегодня я постараюсь максимально информативно и без громоздких математических формул ответить на этот вопрос.

Если вы долгое время работаете в аналитике, часто запускаете эксперименты (или наоборот только начинаете погружаться в это искусство), но никогда не слышали о CUPED'е, то настоятельно рекомендую ознакомиться с принципами его работы. В рамках данной статьи мы сразу приступим к его использованию, поэтому если необходимо, то лучше для начала погрузиться в основы: можно почитать у ребят из X5 Tech тут, а также в статье ребят из Авито. Некоторые фрагменты моего кода взяты из этих статей.

Применимость

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

Собственно, име��но поэтому важно иметь продукт, который существует уже достаточно длительное время. Метод CUPED особенно полезен в ситуациях, когда мы стремимся выявить небольшие эффекты. В случае новых продуктов значительный рост метрики в начале может достигать 20-30%, что делает использование CUPED менее актуальным.

Генерируем данные

Перед проведением наших экспериментов создадим данные вида:

user_id — наши пользователи продукта
metric — значения метрики для каждого пользователя в период эксперимента
covariate — ковариата (значение метрики на предпериоде)

import numpy as np
import pandas as pd

def generate_data(sample_size, correlation, mean=200, sigma=30):
    """
    Генерирует DataFrame с данными для A/B тестирования.

    Parameters:
    — sample_size: Количество образцов (магазинов).
    — correlation: Коэффициент корреляции между метрикой и ковариатой.
    — mean: Среднее значение для метрики и ковариаты.
         sigma: Стандартное отклонение для метрики и ковариаты.

    Returns:
    - DataFrame с уникальными идентификаторами магазинов, метрикой и ковариатой.
    """
    # Установка начального значения генератора случайных чисел для воспроизводимости
    np.random.seed(52) # 🇲🇽
    
    # Определение средних значений и ковариационной матрицы
    means = np.array([mean, mean])
    covariance_matrix = sigma ** 2 * np.array([[1, correlation], [correlation, 1]])
    
    # Генерация многомерного нормального распределения
    data = np.random.multivariate_normal(means, covariance_matrix, sample_size)
    
    # Генерация уникальных идентификаторов пользователей
    user_ids = np.arange(1, sample_size + 1)
    
    # Создание DataFrame
    df = pd.DataFrame({
        'user_id': user_ids,
        'metric': data[:, 0],
        'covariate': data[:, 1]
    })
    
    # Упорядочивание по возрастанию user_id
    df.sort_values(by='user_id', inplace=True)
    df.reset_index(drop=True, inplace=True)
    
    return df

# Генерация данных
mean = 200
sigma = 35
correlation = 0.75
sample_size = 5000

df = generate_data(sample_size, correlation, mean, sigma)
df.head(10)

Создаем id'шки 5000 пользователей, для каждого пользователя генерим метрику и ковариату из нормального распределения со средним значением равным 200, стандартным отклонением 35 и коэффициентом корреляции 0.75.

Вычисление размеров групп

Напишем функцию, которая поможет рассчитать 2 параметра: при заданном значении наблюдений в группе на выходе получим эффект в процентах, который сможем пронаблюдать при заданных α и β; если зададим эффект — получим размер одной группы (тоже для фиксированных α и β).

from typing import Optional, Union
from pingouin import power_ttest

def get_stat(
    array: Union[pd.Series, np.array],
    effect_percent: Optional[float] = None,
    n: Optional[int] = None,
    alpha: Optional[float] = None,
    power: Optional[float] = None,
    alternative: str = 'two-sided',
    paired: bool = False
) -> float:
  
    mean = array.mean()
    var = array.var()  
    cohen_effect = effect_percent * (mean / np.sqrt(var)) if effect_percent is not None else None
    contrast = 'paired' if paired else 'two-samples'
    
    result = power_ttest(
        d=cohen_effect,
        n=n,
        power=power,
        alpha=alpha,
        contrast=contrast,
        alternative=alternative,
    )
    return round(result * (np.sqrt(var) / mean) * 100, 2) if effect_percent is None else result

Сперва рассчитаем необходимое количество наблюдений в одной группе стандартным методом (без использования CUPED):

def get_stat_table(
    array: Union[pd.Series, np.array],
    effect_percent: Union[np.array, list, float, None] = None,
    n: Union[np.array, list, float, None] = None,
    alpha: Union[np.array, list, float, None] = None,
    power: Union[np.array, list, float, None] = None,
    alternative: str = 'two-sided',
    paired: bool = False
) -> pd.DataFrame:
  
    stats = n if n else effect_percent
    result = pd.DataFrame(index=stats)
    
    for p in power:
        for a in alpha:
            outputs = []
            for stat in stats:
                if n:
                    output = get_stat(
                        array=array,
                        effect_percent=None,
                        n=stat,
                        alpha=a,
                        power=p,
                        alternative=alternative,
                        paired=paired
                    )
                else:
                    output = get_stat(
                        array=array,
                        effect_percent=stat,
                        n=None,
                        alpha=a,
                        power=p,
                        alternative=alternative,
                        paired=paired
                    )
                outputs.append(round(output, 2))
            
            result[(f'α={round(a, 3)}', f'β={round(1-p, 3)}')] = outputs
            result.index.name = 'n' if n else 'MDE'
    
    result.columns = pd.MultiIndex.from_tuples(result.columns)
    
    return result

get_stat_table(np.array(df['metric']), 
               effect_percent = [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1], #эффект в долях
               n = None, 
               alpha = [0.01, 0.05, 0.1], 
               power = [0.8, 0.9]) #non-CUPED

Смотрим:

Размеры выборок для различных MDE, α и β (non-CUPED)

Предположим, что мы стремимся выявлять малые эффекты, при этом наш тест должен обладать высокой мощностью и минимизировать вероятность ложноположительных результатов. Из таблицы выберем MDE = 2%, допустимые уровень значимости α = 0.05 и ошибку второго рода β = 0.2: для таких параметров нам необходимо по ~ 1200 юзеров на группу. Многовато... Изначально мы всего сгенерировали 5000 пользователей, и такой тест займет почти половину имеющегося трафика, а нам хочется еще какие-то тесты в параллель запускать. Что делать будем? Увольняться))

И тут к нам приходит спасение: мы можем перейти к новой метрике $\overline{Y}_{\text{CUPED}}$ , которая представляет собой исходную метрику $\overline{Y}$ , преобразованную с помощью нашей ковариаты $\overline{X}$ .

$\overline{Y}_{\text{CUPED}} = \overline{Y} - \theta \times \overline{X}$

$\theta = \frac{\text{Cov}(X,Y)}{\text{Var}(X)}$

Дисперсия преобразованной метрики при этом будет уменьшена на квадрат коэффициента пирсоновской корреляции между исходной метрикой и ковариатой:

$\text{Var}\left(\overline{Y}_{\text{CUPED}}\right) = \text{Var}(\overline{Y})\left(1 - \rho^2\right)$

$\rho$ — коэффициент корреляции Пирсона

В математику, стоящую за этой формулой, можно погрузиться тут: раз и два.

Теперь модифицируем функцию, рассчитав дисперсию согласно формуле:

def get_stat(
    array: Union[pd.Series, np.array],
    effect_percent: Optional[float] = None,
    n: Optional[int] = None,
    alpha: Optional[float] = None,
    power: Optional[float] = None,
    alternative: str = 'two-sided',
    paired: bool = False,
    cuped: Optional[bool] = False #флаг CUPED'a
) -> float:
    mean = array.mean()
    var = array.var() * (1 - (df['metric'].corr(df['covariate']))**2) if cuped else array.var() #считаем дисперсию в зависимости от применения CUPED
    cohen_effect = effect_percent * (mean / np.sqrt(var)) if effect_percent is not None else None
    contrast = 'paired' if paired else 'two-samples'
    
    result = power_ttest(
        d=cohen_effect,
        n=n,
        power=power,
        alpha=alpha,
        contrast=contrast,
        alternative=alternative,
    )
    return round(result * (np.sqrt(var) / mean) * 100, 2) if effect_percent is None else result

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

def get_stat_table(
    array: Union[pd.Series, np.array],
    effect_percent: Union[np.array, list, float, None] = None,
    n: Union[np.array, list, float, None] = None,
    alpha: Union[np.array, list, float, None] = None,
    power: Union[np.array, list, float, None] = None,
    alternative: str = 'two-sided',
    paired: bool = False,
    cuped: bool = False #флаг CUPED'a
) -> pd.DataFrame:
  
    stats = n if n else effect_percent
    result = pd.DataFrame(index=stats)
    
    for p in power:
        for a in alpha:
            outputs = []
            for stat in stats:
                if n:
                    output = get_stat(
                        array=array,
                        effect_percent=None,
                        n=stat,
                        alpha=a,
                        power=p,
                        alternative=alternative,
                        paired=paired,
                        cuped=cuped #флаг CUPED'a
                    )
                else:
                    output = get_stat(
                        array=array,
                        effect_percent=stat,
                        n=None,
                        alpha=a,
                        power=p,
                        alternative=alternative,
                        paired=paired,
                        cuped=cuped 
                    )
                outputs.append(round(output, 2))
            
            result[(f'α={round(a, 3)}', f'β={round(1-p, 3)}')] = outputs
            result.index.name = 'n' if n else 'MDE'
    
    result.columns = pd.MultiIndex.from_tuples(result.columns)
    
    return result

Размеры выборок для различных MDE, α и β (CUPED)

А теперь нам потребуется всего лишь 525 наблюдений в группе для тех же MDE, α, β. Мы уменьшили размеры выборок на 56%! Неплохо, правда?

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

Для доказательства корректности критерия проведем достаточно большое количество
симуляций (пусть будет 10000) A/A- и A/B-тестов и посмотрим на поведение p-value. Будем отбирать в группы по 525 пользователей:

import scipy.stats as sps
from tqdm import tqdm

class ExperimentComparisonResults:
    def __init__(self, pvalue, effect, ci_length, left_bound, right_bound):
        self.pvalue = pvalue
        self.effect = effect
        self.ci_length = ci_length
        self.left_bound = left_bound
        self.right_bound = right_bound

def absolute_ttest(control, test):
    mean_control = np.mean(control)
    mean_test = np.mean(test)
    var_mean_control = np.var(control) / len(control)
    var_mean_test = np.var(test) / len(test)
    
    difference_mean = mean_test - mean_control
    difference_mean_var = var_mean_control + var_mean_test
    difference_distribution = sps.norm(loc=difference_mean, scale=np.sqrt(difference_mean_var))

    left_bound, right_bound = difference_distribution.ppf([0.025, 0.975])
    ci_length = (right_bound - left_bound)
    pvalue = 2 * min(difference_distribution.cdf(0), difference_distribution.sf(0))
    effect = difference_mean
    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)

def cuped_ttest(control, test, control_before, test_before):
    theta = (np.cov(control, control_before)[0, 1] + np.cov(test, test_before)[0, 1]) /\
             (np.var(control_before) + np.var(test_before))

    control_cup = control - theta * control_before
    test_cup = test - theta * test_before
    return absolute_ttest(control_cup, test_cup)

def run_tests(df, num_tests=10000, num_users_per_group=525): #отбираем по 525 юзеров, как и посчитали для CUPED'a
    p_values_aa = []
    p_values_ab = []
    
    unique_users = df['user_id'].unique()
    
    for _ in tqdm(range(num_tests), desc="Running tests"):
        # Случайным образом выбираем 2*num_users_per_group элементов из уникальных пользователей
        selected_users = np.random.choice(unique_users, size=2 * num_users_per_group, replace=False)
        
        # Случайным образом делим выбранных пользователей на контрольную и тестовую группы
        np.random.shuffle(selected_users)  # Перемешиваем
        control_users = selected_users[:num_users_per_group]  # Первые num_users_per_group для контрольной группы
        test_users = selected_users[num_users_per_group:]  # Остальные для тестовой

        # Извлекаем данные для контроля и теста
        control = df[df['user_id'].isin(control_users)]['metric'].values
        test = df[df['user_id'].isin(test_users)]['metric'].values
        control_before = df[df['user_id'].isin(control_users)]['covariate'].values
        test_before = df[df['user_id'].isin(test_users)]['covariate'].values

        # Считаем и сохраняем p-value для A/A и A/B
        result_aa = cuped_ttest(control, test, control_before, test_before)
        result_ab = cuped_ttest(control, test * 1.02, control_before, test_before) #примешиваем те самые 2% эффекта
        p_values_aa.append(result_aa.pvalue)
        p_values_ab.append(result_ab.pvalue)

    return p_values_aa, p_values_ab

pvalues_aa, pvalues_ab = run_tests(df)

# Running tests: 100%|██████████| 10000/10000 [00:37<00:00, 267.97it/s]

import matplotlib.pyplot as plt

def plot_pvalue_distribution(pvalues_aa, pvalues_ab, alpha, beta):
    """Рисует графики распределения p-value для A/A и A/B тестов."""
    # Преобразование списков в массивы NumPy, если это необходимо
    pvalues_aa = np.array(pvalues_aa)
    pvalues_ab = np.array(pvalues_ab)
    
    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.title('Оценка распределения p-value', size=16)
    plt.xlabel('p-value', size=12)
    plt.legend(fontsize=12)
    plt.grid(which='both', linestyle='--', linewidth=0.5)
    
    plt.show()

plot_pvalue_distribution(pvalues_aa, pvalues_ab, alpha=0.05, beta=0.2) # α и β берем те, которые установили перед тестом

Распределение p-value для проведенных тестов

Все работае�� корректно: p-value на А/A-тестах распределены равномерно, а для АБ-тестов зафиксировали в 80% случаев значимые изменения там, где они действительно были (напомню, мы руками добавили эффект, умножив метрику на 1.02).

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

Подводим итог

В современном мире аналитики и A/B-тестирования метод CUPED становится незаменимым инструментом. Этот метод не только позволяет более точно анализировать результаты экспериментов, но и может быть эффективно применен на этапе их дизайна.

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

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

А закончить хочется красивой пословицей:

"Расчет размера групп перед до эксперимента — это не просто формальность, это способ избежать ситуации, когда твои выводы основаны на мнении бабушки и ее кота!"

Да, это тоже ChatGPT придумал...

tg @dezluvv

Размер имеет значение: как исторические данные помогают на этапе дизайна A/B-теста