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

Я не знал с чего начать статью, поэтому попросил 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 = 2%, допустимые уровень значимости α = 0.05 и ошибку второго рода β = 0.2: для таких параметров нам необходимо по ~ 1200 юзеров на группу. Многовато… Изначально мы всего сгенерировали 5000 пользователей, и такой тест займет почти половину имеющегося трафика, а нам хочется еще какие-то тесты в параллель запускать. Что делать будем? Увольняться))

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

,

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

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

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

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

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

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 на А/A-тестах распределены равномерно, а для АБ-тестов зафиксировали в 80% случаев значимые изменения там, где они действительно были (напомню, мы руками добавили эффект, умножив метрику на 1.02).

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

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

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

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

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

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

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

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

tg @dezluvv