Кластерные эксперименты, часть 1. Юниты тестирования и анализа, симуляции зависимых данных.

Про юниты тестирования, ratio-метрики, вот это вот все.

Алоха 🖐🏾! Сегодня у нас первая статья про кластерные эксперименты, разберем что это такое и сделаем базовые симуляции.

Линеаризацию, перевзвешивание и дельта-метод часто обсуждали лет 5 назад, обычно ввиду их применения к анализу ratio-метрик и/или увеличения чувствительности. Да, уже 5 лет назад вышла отличная статья от команды VK Core ML, в которой разобрано большинство этих методов.

Но к сожалению, вокруг этих методов уже скопилось несколько мифов, которые пора развеять. Мало кто рассказывает о том, для чего вообще придумали дельта-метод (спойлер, не для ratio-метрик), а для чего линеаризацию (почему-то Бабушкин тут c ее помощью увеличивает чувствительность, но это мы тоже разберем 🤫). Так что статья от VK определенно требует продолжения.

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

Содержание:

· 1. Юниты тестирования и анализа ∘ Почему важно различать юниты · 2. Количественная оценка ∘ Intra-Class Correlation (ICC) · 3. Метрики в кластерных А/В-тестах ∘ Про ratio-метрики ∘ Симуляции ∘ Датасет ∘ A/A-тесты ∘ Кластерные тесты ∘ Дельта-метод и двойное среднее ∘ Просто глобальное среднее · 4. Скоррелированные данные · 5. Как считать деньги ∘ Числитель, зависящий от знаменателя · 6. Не-ratio метрики

1. Юниты тестирования и анализа

Начнем с азов. Немного терминологии:

• Юнит тестирования (unit of randomisation, sample unit), R — это бизнес-сущность, по которой мы раздаем тест и вешаем метки контрольной/таргетной группы. Обычно это пользователь/cookie/девайс.
• Юнит анализа (population unit), A — сущность, на чьем уровне мы считаем метрики. Например, в метрике ARPU юнит анализа — клиент (хотя тест может быть роздан на регион/девайс/…). Вообщем это то, что стоит в знаменателе при расчете метрики.

Обычно у наших бизнес сущностей есть какая-то иерархия: сессия < клиенто-день < клиент < регион.

В дальнейшем мы будем использовать 2 примера:

• E-commerce: юнит тестирования — пользователь, юнит анализа — сессия или покупка.
• EdTech: юнит тестирования — учитель, юнит анализа — студент.

Юниты тестирования и анализа могут не совпадать. Во-первых, так приходится делать чтобы не было не-консистентности клиентского опыта. Пример e-commerce: представь, что мы тестируем конверсию в интернете-магазине (юнит анализа — сессия), и пользователь при повторном заходе увидит другую цену/дизайн, или что учитель (пример EdTech) будет учить половину учеников по одной программе, а вторую половину по другой. Во-вторых, пресловутый сетевой эффект. Я уже писал об этом в прошлой статье (второй абзац), и можешь прочитать прекрасную статью от Lyft.

Почему важно различать юниты

1. В простом кейсе (R = A) имеем обычный А/В-тест. Раздаем на клиента и метрики считаем тоже по-клиентно. Тут прекрасно работают обычные стат-критерии, например t-тест.
2. Если R < A (юнит тестирования более гранулярен, чем юнит анализа) — тут все плохо, т.к. в юните анализа (его метрике) будут как наблюдения из теста, так и из контроля, что сводит на нет весь смысл A/B-теста. Так делать нельзя, т.е. нельзя рандомизировать сессии и считать поюзерные метрики! Один из примеров этой проблему в тестах на двустороннем маркетплейсе я разбирал в своем докладе на Матемаркетинге.
3. Теперь о кейсе R > A. Думаю, тебе прекрасно знакома данная ситуация: несколько юнитов анализа (сессии) принадлежат одному юниту тестирования (клиенту). Данный вариант рандомизации в большинстве источников называется кластерным. Также, иногда его называют гнездовым. Поскольку сессии одного юнита между собой зависимы, обычные стат-критерии на уровне сессии не подходят! Если ты еще не в курсе этой проблемы — можно начать со статьи от Scyscanner.

Картинка 1. Не-кластерная и Кластерная рандомизации

Если лень читать статью, на Картинке 1 интуитивное объяснение, (котики 😺 из популярной книги). Серыми квадратиками отмечены кластеры (по цвету кота), и пусть желтые коты — это юниты с более высоким значением метрики (допустим, денег), а светлые и темные — с более низким. Вхождение кота в конкретный кластер (цвет котика) — это ковариат, т. е. свойство котика, влияющее на целевую метрику. При рандомизации по юнитам анализа (левая часть), каждый кластер разделится примерно пополам, и между контролем и таргетом в среднем доля каждого класса будет совпадать, а значит совпадет и целевая метрика по группам. Все это я рассказывал в первой половине своей статьи про causal inference. Если же мы будем рандомизировать кластеры (правая часть), то все богатые попадут в одну группу, а бедные — в другую.

Мы все еще можем использовать среднеквадратичное отклонение выборки как несмещенную оценку дисперсии между юнитами анализа. Однако, ЦПТ перестает работать из-за корреляции между юнитами. А значит дисперсию статистики (среднего) мы не можем посчитать классическим способом. Именно это и проводит к тому, что считать ее на уровне котика (т.н. ванильным методом рассчета дисперсии) нельзя, и это инвалидирует использование классического t-теста ввиду слишком высокого FPR.

2. Количественная оценка

Другое объяснение может дать Картинка 2. На ней показана формула среднеквадратичного отклонения разницы двух случайных величин, например средних значений контроля и таргета. В случае сложения дисперсий там стоит “плюс”.

Картинка 2. Дисперсия зависимых выборок

В обычном А/В нет нарушения SUTVA, и ковариация (последний член под корнем) равна нулю. В кластерном анализе это не так, поскольку если один желтый котик оказался в одной из групп, все остальные желтые котики окажутся в той же. Т.е. в контроль попадут все богатые, а таргету останутся только бедные — поэтому ковариация не равна нулю, а std/дисперсия увеличиваются.

Степень увеличения дисперсии имеет свое название — **design effect (deff). Данное понятие встречается не только в кластерных А/В-тестах, но и в любых экспериментах, в которых семплирование отлично от простого рандомайзера (в котором все юниты имеет одинаковую вероятность оказаться в одной из групп). Deff — это отношение дисперсии экспериментальной статистики к дисперсии “эталонного” рандомизатора (т.н. [Simple Random Sampling](https://en.wikipedia.org/wiki/Simple_random_sample), или SRS**). Ниже покажем на сколько мы ошибаемся в ‘ванильном” методе в кластерном тесте.

А пока важно запомнить: сама проблема, которую решают дельта-метод и линеаризация — это зависимость данных в кластерном эксперименте, а ratio-метрики это лишь следствие. Подобная проблема возникает также в метриках-квантилях (например медиана) и при использовании других стат-критериев (не только t-теста). Они решаются например кластерным бутстрапом.

Intra-Class Correlation (ICC)

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

Если мы возьмем всех пользователей и объединим их в рандомные кластеры, то при достаточном кол-ве юнитов центры кластеров будут равны, а дисперсия — определяться по ЦПТ.

То есть величина нашей недооценки дисперсии зависит от того, насколько близки между собой юниты в кластере относительно их близости с юнитами из других кластеров. Логично было бы посчитать корреляцию между юнитами внутри кластера, но как это сделать? Для этого введен специальный параметр — Коэффициент Внутри-кластерной Корреляции, или в англоязычном мире **Intra-Class Correlation (ICC**).

Картинка 3. Для любой темы есть индийские парни на Youtube

По картинке 3 становиться понятен бизнесовый смысл. Например, у нас есть выборка из 100 пользователей, но в этой выборке есть 20 пар близнецов (вырожденный случай когда блинецы прям одинаковы и ICC в паре равен 1). В таком случае мы возьмем из каждой пары по одному близнецу, а второй нам не нужен т.к. не несет в себе никакой информации (полностью аналогичен тому что уже есть). Итого, когда мы по выборке будем считать дов-интервал, мы должны в размере выборке брать не 100, а уже 80 человек (т.н. Effective Sample Size, ESS = Deff Sample Size*). В таком случае у нас возникает design effect = 100/80. В других кейсах, когда 0 < ICC < 1, мы бы получили Effective Sample Size что-то среднее между 80 и 100, то есть среднее между кол-вом кластеров при ICC=1 и общему кол-ву юниту анализа при ICC=0.

Изначально для расчета ICC пытались адаптировать корреляцию Пирсона, однако она требует упорядоченных пар юнитов, а также не умеет работать с кол-вом кластеров больше 2х. Сейчас используется такая мера — доля межкластерной дисперсии от общей.

Картинка 4. Формула коэффициента внутриклассовой корреляции

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

Картинка 5. Разложение дисперсии в ANOVA

На википедии доказано, что ICC совпадает с корреляцией Пирсона юнитов внутри одного кластера (ковариация делить на произведение дисперсий). Однако, в отличии от корреляции, данный коэффициент принимает значения только от 0 до 1.

Как посчитать ICC в Python/R — посмотри мой ответ на этот вопрос на cross-validated. Можно использую пакеты [pingouin](https://pingouin-stats.org/generated/pingouin.intraclass_corr.html) (Python) и [irr](https://www.rdocumentation.org/packages/irr/versions/0.84.1/topics/icc) (R), но они работают только на кластерах одинакового размера (все это увидим в следующих частях серии). Вообще, тут можно найти расчет на чистом Python, но ниже в блоке 4 (про скоррелированные данные) я покажу как его считать.

В отличии от Корреляции Пирсона, ICC можно посчитать и в случае кластеров разного размера (клиенты имеют разное количество чеков/сессий). ICC, определяемый через ANOVA, работает всегда:

Картинка 6. Доказательство разложения в ANOVA

К сожалению, ICC — это общий показатель для всех кластеров, и в нем не учитывается похожесть юнитов отдельного кластера.

3. Метрики в кластерных А/В-тестах

Давай поймем, зачем нам вообще это знать и как применить в А/В. Я буду использовать терминологию из статьи Яндекса:

• Метрика — это любая количественная величина, измеряемая или рассчитываемая на уровне юнита анализа (важно, не кластера!)
• Evaluation Statistic (ES) — это статистика (функция от выборки) на этой метрике, например среднее, медиана, std и тд.
• OAC — метрика+ES называются Overall Evaluation Criterion (OEC); а если к ним добавить статистический критерий, который будет использовать для расчёта стат-значимости разницы в OEC, то все вместе это называется Overall Acceptance Criterion (OAC).

Ron Kohavi нас учит, что считать эффект в А/В-тесте для средней величины (например, среднего чека) можно двумя разными способами (2 разных OEC):

Картинка 7. Разные виды средних

• т.н. “наивное, или глобально среднее”, когда мы суммируем метрики от всех юнитов и делим на их общее количество юнитов во всех кластерах.
• т.н. “нормализованное, или двойное среднее”, когда мы считаем среднее значение по каждому кластеру и затем усредняем уже кластеры.

Хороший пример наивного среднего — это глобальный CTR, когда мы делим все клики на объявления на все показы. В качестве примера нормализованного среднего, также называемого “двойным”, можно привести поюзерный CTR, когда вначале считаем отношение всех кликов к показам для пользователя (т.е. кластера), а затем усредняем между пользователями. Тут стоит дать комментарий — в известном докладе Яндекса наивным наоборот назвали нормализованное среднее, поэтому я вижу что на собесах люди ошибаются, но не нужно их путать!

Как это соотносится к сказанному выше:

1. Наивное среднее обладает тем самым недостатком “наивной” дисперсии в кластерном тесте, поскольку юниты анализа зависимы, и для OAC не подходит t-критерий.
2. Нормализованное среднее фактически делает юнитом анализа кластер, изменяя метрику (теперь метрика по-кластерная — это среднее внутри кластера), но при этом юниты анализа уже независимы и можно использовать t-критерий.

Вопрос: какое среднее выбрать? Кажется, что второй вариант лишен недостатков и стоит всегда использовать его. Но на самом деле нет, ответ зависит от ситуации. Для понимания сути вопроса обратимся к статье от Ebay, где рассматривается общий случай средней метрики, взвешенной по кластерам:

Картинка 8. Взвешенное среднее метрики

Если wi равно 1/k для всех кластеров, то Rweighted = Rb, т.е. мы получаем нормализованное среднее. Если wi = n_i/N, то Rweighted=Ra и взвешенное среднее равно наивному.

Теперь становится понятным бизнес-смысл каждой из величин. Нормализованное среднее дает каждому кластеру одинаковый вес в общем среднем, в то время как в наивном вес каждого кластера равен его размеру. В ситуации, когда нам одинаково важен каждый кластер, стоит использовать нормализованное среднее: допустим, мы сделали новый дизайн и смотрим на среднее время сессии в поюзерном тесте (каждый пользователь для нас одинаково важен). Однако, если мы считаем денежные метрики (средний чек в E-commerce, или доход на студента в тесте на учителя в EdTech), правильнее будет использовать наивное среднее, поскольку эта метрика точно выражает среднюю метрику из юнит-экономики: все деньги, деленные на общее кол-во юнитов.

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

Про ratio-метрики

Раньше можно было встретить множество статей о ratio-метриках. Откуда они взялись и что это такое?

Несколько команд (Microsoft, Uber, Yandex, Facebook) рассказывают об использовании уже имеющихся или о разработке новых методов для анализа OEC в виде наивного среднего, поэтому можно предположить что это правильный выбор при оценке средних значений метрик. Чтобы перейти от наивного среднего к ratio, можно воспользоваться следующим наблюдением: количество юнитов в каждом кластере (ni) — это тоже случайная величина, а значит ее можно рассматривать как отдельную метрику. Обозначим за X с 3мя точками суммарное значение метрики по всем юнитам кластера:

Картинка 9. Переход к ratio-метрике

Получается, что наивное среднее есть отношение средних 2х других метрик, для каждой из которых как юнитом анализа, так и sample-юнитом является кластер. Для примера EdTech — этими 2мя метриками будут Доход на учителя и Размер класса (кол-во студентов).

Важно понимать, что наивное среднее в кластерном эксперименте — это только частный случай ratio-метрик.

Симуляции

Для начала в целях повторяемости введем импортнем пакеты и сделаем пару функций.

import numpy as np
import pandas as pd

import collections
import math

from scipy import stats
from scipy.stats import multivariate_normal
import statsmodels.stats.weightstats as ws
# будем так делить на КГ/ТГ, хорошо показало себя в А/А-тестах
# см. https://koch-kir.medium.com/%D0%BE%D1%87%D0%B5%D0%BD%D1%8C-%D0%BC%D0%BD%D0%BE%D0%B3%D0%BE-%D1%81%D0%BF%D0%BE%D1%81%D0%BE%D0%B1%D0%BE%D0%B2-%D1%83%D1%81%D0%BA%D0%BE%D1%80%D0%B8%D1%82%D1%8C-pandas-db74efceb086
from sklearn.model_selection import train_test_split 

import matplotlib.pyplot as plt
import seaborn as sns
from tqdm.notebook import tqdm

N = 10000    # размер выборки (кол-во пользователей/кластеров)
B = 1000     # кол-во подвыборок для бутстрэпа
alpha = 0.05 # порог p-value

def plot_pvalue(pvalues):
    # Рисуем распределение QQ
    X = np.linspace(0, 1, B)
    for name, pvalues in pvalues.items():
        Y = [np.mean(pvalues < x) for x in X]
        plt.plot(X, Y, label=name)
    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()
    plt.show()

def hist_pvalue(pvalues):
    # Гистограмма распределений
    plt.hist(pvalues)
    plt.title('Оценка распределения p-value', size=16)
    plt.xlabel('p-value', size=12)
    plt.ylabel('Частота', size=12)
    plt.show()

def check_p_values(p_values):
    # Проверяем нулевую гипотезу о том, что p-value распределен неравномерно
    error_rate = np.mean(np.array(p_values) < alpha)
    print(f'Доля ошибок первого рода: {error_rate:0.2f}')

    if error_rate < 0.05:
        print("Различия не статистически значимы!")
    else:
        print("Различия статистически значимы!")
    plot_pvalue({'A/A': p_values})
    hist_pvalue(p_values)

Датасет

Моя идея для датафрейма заключается в следующем — выбираем какое-то натуральное число No, которое равно максимальному числу строчек (сессий, заказов) для пользователя. Далее мы генерим для каждого пользователя из многомерного нормального распределения это кол-во строчек.

No = 15      # Максимальное кол-во заказов юзера/кластера
avg = 1500   # Средний чек заказов внутри юзера/кластера
sigma = 500  # Стандартное отклонение заказов внутри юзера/кластера

cv = np.identity(No)*sigma

x = np.random.multivariate_normal(mean=np.full(No, avg), cov=cv, size=1)
for i in tqdm(range(N-1)):
    x_i = np.random.multivariate_normal(mean=np.full(No, avg), cov=cv, size=1)
    x = np.vstack([x, x_i])

# Датафрейм для N пользователей с No колонками заказов/сессий
df = pd.DataFrame(
  x,
  columns = ['item' + str(i) for i in  range(No)],
  index = ['user_' + str(i) for i in  range(N)]
).reset_index(inplace=True)

# Раскладываем по строчкам
y = pd.melt(
  df,
  id_vars=['index'],
  value_vars=['item' + str(i) for i in  range(No)]
)

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

# По определению - https://ru.wikipedia.org/wiki/%D0%9B%D0%BE%D0%B3%D0%BD%D0%BE%D1%80%D0%BC%D0%B0%D0%BB%D1%8C%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5
# Если случайная величина имеет логнормальное распределение,
# то её логарифм имеет нормальное распределение.

# Значит в обратную сторону:
# Экспонента от нормального имеет логнормальное распределение.
s = np.random.normal(1.5, 0.5, N) ; n = np.exp(s).astype(np.int8)# + 10

data_dict = {
    'user_id': y['index'].unique(),
    'n_orders': n
}
df = pd.DataFrame(data_dict) # датафрейм с кол-вом заказов для каждого юзера

# Ограничиваем кол-во заказов от 1 до No
df['new_n_orders'] = df['n_orders'].apply(lambda x: x if x <= No else No)
df['new_n_orders'] = df['new_n_orders'].apply(lambda x: 1 if x <= 0 else x)

df = y.merge(
    df,
    how = 'inner',
    left_on = 'index',
    right_on = 'user_id'
)

# нумеруем все заказы
df['rn'] = (
    df.sort_values(['user_id', 'variable',], ascending=[True,True])
    .groupby('user_id')
    .cumcount() + 1
)

# оставляем только первые n заказов
data = df.loc[df.new_n_orders >= df.rn].copy()
data.drop(columns=['index', 'variable'], inplace=True)
data.columns = ['cheks', 'user_id', 'n_orders', 'new_n_orders', 'rn']

Картинка 10. Полученный дата-сет

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

A/A-тесты

Теперь перейдем непосредственно к симуляциям.

# Для проверки проводим классический АА тест 1000 раз
# При этом сохраняем p-value всех тестов для картинки распределение

pvalues1 = []
for el in tqdm(range(B)):
    # Разбиение выборки на группы
    half = int(len(data['cheks']) / 2)
    list = data['cheks'].tolist()
    np.random.shuffle(list)

    a1 = list[0:half]
    a2 = list[half:len(list)]

    # Сохранение p-value
    pvalue = stats.ttest_ind(a1, a2).pvalue
    pvalues1.append(pvalue)

check_p_values(pvalues1)

Получаем данные, как на картинке ниже:

Картинка 11. Симуляции классического t-теста

Кластерные тесты

Все бы хорошо, но тут мы просто делим строчки пополам без учета того, к какому пользователю они относиться. Для правильного подсчета воспользуемся одним из двух способов — глобальным или двойным средним. Вначале сделаем поюзерный дата-фрейм:

# Средний чек по юзерам
avg_by_user = data.pivot_table(
  index='user_id',
  values='cheks',
  aggfunc=['sum', 'count']
)
avg_by_user = avg_by_user['sum'].join(
  avg_by_user['count'],
  how='inner',
  lsuffix='_sum',
  rsuffix='_count'
).reset_index()
avg_by_user['cheks_mean'] = avg_by_user['cheks_sum'] * 1.0000 / avg_by_user['cheks_count']

Теперь у нас в дата-фрейме user_id является ключом; в колонках указано кол-во чеков, их сумма и средний чек.

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

Для каждого способа посчитать среднее используем адаптированный t-тест. Сделаем функцию для подсчета дисперсии дельта-методом.

Дельта-метод и двойное среднее

Ниже код для подсчета дельта-методом дисперсии среднего одной группы. Самое интересное — если не указать что степени свободы равны 1, то мы получим смещенные оценки и часть симуляций отрабатывали бы неправильно.

def get_ratio_var(num, denom):
    """
    Возвращает оценку дисперсии для ratio.

    num: float,int, метрика в числителе.
    denom: float,int, метрика в знаменателе.
    """
    cov = np.cov(num, denom, ddof=1)[0, 1]
    var = (
        (np.std(num,ddof=1) ** 2) / (np.mean(denom) ** 2) +
        (np.mean(num) ** 2) * (np.std(denom,ddof=1) ** 2) / (np.mean(denom) ** 4) -
        2 * np.mean(num) / (np.mean(denom) ** 3) * cov
    ) #/ len(num)

    return var

Так мы считаем только почти-дисперсию чеков одной группы (КГ или ТГ), для полной дисперсии нам нужно во-первых разделить на кол-во кластеров (как в ЦПТ), во-вторых сложить получившиеся дисперсии средних КГ и ТГ.

pvalues = []

for el in tqdm(range(B)):
    # Разбиение выборки на группы
    s1, s2 = train_test_split(avg_by_user, test_size = 0.5, random_state=el*10, shuffle=True)

    num0 = s1['cheks_sum']
    denom0 = s1['cheks_count']
    num1 = s2['cheks_sum']
    denom1 = s2['cheks_count']

    ratio0 = np.sum(num0) / np.sum(denom0)
    ratio1 = np.sum(num1) / np.sum(denom1)
    v0 = get_ratio_var(num0, denom0)
    v1 = get_ratio_var(num1, denom1)
    n0 = len(num0)
    n1 = len(num1)

    # Выше описали что тут происходит
    se = np.sqrt(v0/n0 + v1/n1)
    # так считается кол-во степеней свободы в тесте Уэлча
    # https://ru.wikipedia.org/wiki/T-%D0%BA%D1%80%D0%B8%D1%82%D0%B5%D1%80%D0%B8%D0%B9_%D0%A3%D1%8D%D0%BB%D1%87%D0%B0
    df_ = (v0/n0 + v1/n1)**2 / ( (v0/n0)**2 / (n0-1) ) + ( (v1/n1 )**2 / (n1-1) )

    delta = ratio1 - ratio0
    statistic = delta / se
    # Делаем "вручную" тест Уэлча на агрегированных данных
    pvalue = 2 * ( 1 - stats.t.cdf(np.abs(statistic), df=math.ceil(df_)) )
    pvalues.append(pvalue)

check_p_values(p_values)

Небольшой дисклеймер: не так давно появлилась статья, в которой говорятся что использовать по умолчанию тест Уэлача (он же equal_var = False в ttest_ind в scipy и statsmodels) не стоит. Так вот — успокойтесь, на больших выборках в сотни-тысяч пользователей никакой проблемы с ним нет.

Картинка 12. Симуляции дельта-метода

Проверим что линеаризация отрабатывает также как и дельта-метод:

pvalues = []

for el in tqdm(range(B)):
    # Разбиение выборки на группы
    s1, s2 = train_test_split(avg_by_user, test_size = 0.5, random_state=el*10, shuffle=True)

    # линеаризация
    alpha_coeff = s1['cheks_sum'].sum() / s1['cheks_count'].sum()
    s1['aov_lin'] = s1['cheks_count']*(s1['cheks_mean'] - alpha_coeff)
    s2['aov_lin'] = s2['cheks_count']*(s2['cheks_mean'] - alpha_coeff)

    _, pvalue = stats.ttest_ind(s1['aov_lin'], s2['aov_lin'])
    pvalues.append(pvalue)

check_p_values(pvalues)

Картинка 13. Ожидаемо, линеаразация таже сошлась.

Также, сделаем классический t-тест для двойного (нормализованного) среднего, ведь оно не подвержено проблеме скоррелированных выборок!

pvalues = []

for el in tqdm(range(B)):
    # Разбиение выборки на группы
    half = int(len(avg_by_user['cheks_mean']) / 2)
    list = avg_by_user['cheks_mean'].tolist()
    np.random.shuffle(list)

    a1 = list[0:half]
    a2 = list[half:len(list)]

    # Рассчитываем p-value
    pvalue = stats.ttest_ind(a1, a2).pvalue
    pvalues.append(p_value)

check_p_values(p_values)

Даже не стал приводить последние графики. Ожидаемо, и дельта-метод для глобального, и t-критерий для двойного дадут красивые картинки с равномерным p-value, аналогичные тем что были в самом первом примере. Все это вы уже читали в других статья. Но теперь давайте перейдем уже к интересным примерам.

Просто глобальное среднее

Итак, попробуем поделить пользователей на 2 группы по уникальным useri_id, потом сделать по группам чеков t-тест (глобальное среднее). Как нам говорят в других статьях — не должны сойтись А/А-тесты, ведь это ratio-метрика. Точно ли это так:

pvalues = []
for el in tqdm(range(B)):
    # Разбиение выборки на группы
    half = int(data['user_id'].nunique() / 2)
    list = data['user_id'].unique()
    np.random.shuffle(list)

    a1 = data.loc[data.user_id.isin( list[0:half] ), 'cheks']
    a2 = data.loc[data.user_id.isin( list[half:] ), 'cheks']
    # Сохранение p-value
    pvalue = stats.ttest_ind(a1, a2).pvalue
    pvalues.append(pvalue)

check_p_values(p_values)

И как же получилось, что у нас сошлись А/А а не завышенные ошибки первого рода?

Картинка 14. Симуляции классического t-теста

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

4. Скоррелированные данные

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

rho = 0.2    # Это и есть коэффициент корреляции заказов внутри пользователя
No = 15      # Максимальное кол-во заказов юзера/кластера
avg = 1500   # Средний чек заказов внутри юзера/кластера
sigma = 500  # Стандартное отклонение заказов внутри юзера/кластера

# это старая ковариационная матрица
# cv = np.identity(No)*sigma
# а вот новая
cv = (np.full((No, No), rho) + np.identity(No)*(1-rho))*sigma

x = np.random.multivariate_normal(mean=np.full(No, avg), cov=cv, size=1)
for i in tqdm(range(N-1)):
    x_i = np.random.multivariate_normal(mean=np.full(No, avg), cov=cv, size=1)
    x = np.vstack([x, x_i])

# Дальше весь остальной код, который был во второй ячейке.
.....

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

Картинка 15. По-юзерный t-тест на скоррелированных данных.

Ну и для понимания давайте посчитаем корелляцию на наших внутри-пользовательских данных:

def calculate_icc(data, cluster_col, value_col):
    # Разбиваем данные по кластерам
    clusters = data.groupby(cluster_col)[value_col].apply(np.array).values

    # Вычисляем ANOVA
    n_groups = len(clusters)
    n_obs = sum(len(c) for c in clusters)

    # Общее среднее
    grand_mean = np.concatenate(clusters).mean()

    # Межгрупповая дисперсия (MSB)
    ss_between = sum(len(c) * (c.mean() - grand_mean)**2 for c in clusters)
    df_between = n_groups - 1
    ms_between = ss_between / df_between

    # Внутригрупповая дисперсия (MSW)
    ss_within = sum(np.sum((c - c.mean())**2) for c in clusters)
    df_within = n_obs - n_groups
    ms_within = ss_within / df_within

    # Поправка на неравные размеры кластеров
    n_i_squared = sum(len(c)**2 for c in clusters)
    k0 = (n_obs**2 - n_i_squared) / (n_obs * (n_groups - 1))

    # ICC (формула Shrout & Fleiss, 1979)
    icc = (ms_between - ms_within) / (ms_between + (k0 - 1) * ms_within)

    return icc

icc = calculate_icc(y, 'index', 'value')

Картинка 16. Корелляция внутри пользователя.

Как видите, корреляция равна 0.2, то как мы ее и задавали. Если бы мы посчитали это для первоначального датафрейма, то получили бы ICC равным примерно 0. Можете выгрузить свои данные из БД и посчитать этой функцией ICC для них — если он около нуля, то никакой дельта-метод вам не нужен, используйте классический t-тест.

5. Как считать деньги

Итак, назовем это парадоксом среднего и рассмотрим его на примере e-commerce.

Рон К. в книге “Доверительное А/В-тестирование” раскладывает кол-во поисковых запросов в месяц как произведение MAU, среднего кол-ва сеансов на пользователя и среднего кол-ва запросов на сеанс. Подобно этому, мы можем разложить GMV на денежные метрики:

Картинка 17. Как мы раскладываем деньги через наивное среднее

Можем сделать вывод, что изменение в наивном (глобальном) среднем коррелирует с изменением в GMV в случае, если кол-во заказов на пользователя не изменяется. VK в свой статье также отмечали, что в большинстве случаев предполагали что не изменяют знаменатель в CTR.

Попробуем сделать тоже самое в случае нормализованного среднего:

Картинка 18. Как мы раскладываем деньги через двойное среднее

Если вдруг вы не поняли что происходит на картинке выше — мы попробовали разложить GMV через двойное среднее, и оказалось что там есть еще один элемент в сумме — ковариация кол-во заказов и среднего чего на пользователя. То есть — если средней чек пользователя зависит с кол-во заказов, то GMV и двойное среднее в эксперименте могут быть разнонаправлены. Например, когда богатые покупают и дороже, и чаще. Или наоборот когда у пользователя есть выбор — покупать все сразу большим чеком или покупать чаще но мелким. Сейчас покажу на симуляциях.

Числитель, зависящий от знаменателя

Давай промоделируем данные: возьмем логнормальное кол-во заказов пользователя, а AOV пользователя — через нормальное распределение.

[embed]Код 1. Данные для симуляции

Представим такую ситуацию: наше изменение дало гетерогенный эффект — кому-то мы улучшили UX, кому-то ухудшили, но в среднем на клиентские метрики мы не повлияли. В таком случае наш анализ не должен показать стат-значимых отличий. Давай смоделируем такую ситуацию. На gist’е выше добавляем к AOV каждого пользователя эффект, распределенный нормально со средним 0. Оба средних на A/A (нормализованное и наивное) дадут хороший FPR, т.е. немного меньше 5% при p-value=0.05 (см. Код 2). Можно сделать упражнение и проверить, что суммарная выручка на этих 1000 итерациях в среднем также не увеличилась.

[embed]Код 2. Симуляции нулевого эффекта

Оба графика выше дадут равномерное распределение p-value.

Но, у нас может быть неслучайный гетерогенный эффект на AOV: средние также 0, но при этом он распределен неравномерно. Одна группа пользователей всегда имеет положительный эффект, а другая отрицательный. Если эти группы скорее всего имеют разный вклад в общие деньги, наш анализ уже должен показать стат-значимые отличия. Давай смоделируем! Добавим нормальный эффект со средним 0 к AOV пользователя, но так чтобы более активные (с большим кол-вом заказов) получили больший эффект. Самый простой вариант — преобразовать кол-во заказов в нормальное распределение монотонной функцией, и тогда клиент с большим кол-во заказов получит больший эффект:

[embed]Код 3. Неравномерный нулевой эффект.

Нормализованное среднее все еще имеет FPR~5%, в том время как наивное среднее в основном дает эффект. А, эффект на деньги, как мы понимаем, тоже есть (условно у пользователя с одним заказом стало -10% каждый чек, а у пользователя с 20 заказами +10% каждый чек, очевидно сумарно деньги вырастут).

Наивное среднее сонаправлено с общим кол-во денег — мы увидим, что начнем зарабатывать больше 🤑. В то же время нормализованное среднее проверяет совсем другую гипотезу. И используя его можно отвергнуть хорошие фичи, или, что еще хуже, раскатить в прод ухудшающие изменения💩. Если тебе кажется, что тут есть математический диссонанс, то обратись к последнему разделу статьи.

Метод выше называется добавлением шума, об говорил Бабушкин в своем видео. Это правильный способ делать А/А тесты!

6. Не-ratio метрики

Сегодня мы узнали, что ratio-метрики — это частный случай метрик в кластерном эксперименте (глобальное среднее). На самом деле могут быть другие и другие типы метрик, например медиана.

Если апроксимировать медиану или другой квантиль какой-либо гладкой функцией (например, функцией включения), то ее можно будет разложить в ряд Тейлора, а значит и применить дельта-метод:

Картинка 19. Дельта-метод для медианаы

А теперь сделаем это на Python.

from statsmodels.nonparametric.kde import KDEUnivariate

def cluster_median_variance(data, clusters, median_est):
    """Оценка дисперсии медианы с учетом кластеров."""
    # Оценка плотности в медиане
    kde = KDEUnivariate(data)
    kde.fit()
    f_median = kde.evaluate(median_est)[0]

    # Функция влияния для медианы
    influence = (np.array(data) <= median_est) - 0.5
    influence /= f_median

    # Учет кластеров
    cluster_sums = pd.Series(influence).groupby(clusters).sum()
    n = len(data)
    variance = (cluster_sums ** 2).sum() / (n ** 2)
    return variance

# Генерация данных
np.random.seed(42)
cluster_ids = np.repeat([1, 2, 3, 4, 5], 10)  # 5 кластеров
data_A = np.random.normal(10, 2, 50)  # Группа A
data_B = np.random.normal(11, 2, 50)  # Группа B

# Медианы
median_A = np.median(data_A)
median_B = np.median(data_B)
delta = median_B - median_A

# Дисперсии с учетом кластеров
var_A = cluster_median_variance(data_A, cluster_ids, median_A)
var_B = cluster_median_variance(data_B, cluster_ids, median_B)
se_delta = np.sqrt(var_A + var_B)

# Доверительный интервал (95%)
ci_lower = delta - 1.96 * se_delta
ci_upper = delta + 1.96 * se_delta

print(f"Разность медиан: {delta:.3f}")
print(f"95% ДИ: [{ci_lower:.3f}, {ci_upper:.3f}]")

Итог (повторим ключевые пункты)

• Ratio-метрики это частный случай метрик в кластерных экспериментах, где отличаются между собой юниты анализа и тестирования. Главная проблема таких экспериментов — то что юниты анализа внутри кластера как правило зависимы (их метрики скоррелированы) между собой.
• Можно использовать как наивное (глобальное), так и нормализованное (двойное) среднее, у каждого есть свои плюсы, но деньги лучше считать глобальным.
• Правильнее делать не просто А/А, а добавлять шум с нулевым средним, в том числе шум который зависит от каких-либо параметров метрик. Стат-критерии должны правильно отрабатывать в этих кейсах.
• Ну и главное— чтобы использовать данные методы, нужно вначале сделать проверку на отсутствие корреляции между кол-во элементов в кластере и распределением метрики внутри кластера (статья Deng в 2011г), и вдобавок i.i.d..

P.S.: Если тебе понравилась стать, в качестве благодарности можешь купить мне кофе на https://www.donationalerts.com/r/stats_data_ninja или https://buymeacoffee.com/koch.