Приветствуем всех читателей! Сегодня команда Ad-Hoc аналитики X5 Tech приоткроет дверь в увлекательный мир A/B-тестирования Causal Inference. С момента написания предыдущей статьи прошло уже 4 года. За это время наш подход к оценке инициатив значительно эволюционировал. Мы собирали бизнес-кейсы, изучали научную литературу, экспериментировали с реальными данными и в итоге пришли не только к другой модели для оценки эффекта, но и изменили методологию в целом. 

Кстати, методология позволила перейти не просто к стандартизации, но и к автоматизации оценки и дизайна экспериментов. На эту тему наши коллеги – команда A/B-платформы – выпустили свою статью, которую мы обязательно рекомендуем к прочтению. 

Но зачем нам пришлось столько всего поменять? Что не так с привычными A/B-тестами? На самом деле, с ними всё в порядке, просто они плохо дружат с жестокой действительностью тестирования в оффлайн-ритейле. 

Возьмём двухвыборочный t-тест и вспомним одну из его ключевых предпосылок – выборки взяты случайным образом из генеральной совокупности. И эта предпосылка, к сожалению, в нашем случае работает довольно редко.

Есть заказчик Вася, сотрудник X5, который обратился к нам с просьбой оценить его инициативу по повышению выручки магазинов с помощью умной системы скидок. Идея, конечно, здравая, и мы можем помочь с пилотом (так мы называем любую инициативу). Тут же в голову приходит мысль о том, что нужно засэмплировать репрезентативную выборку, внедрить на ней эту систему скидок, а потом сравнить среднюю выручку в тестовых магазинах (мы называем их пилотной группой) со средней выручкой в таким же образом случайно подобранных контрольных магазинах. Но в процессе деловой встречи выясняется, что Вася уже провёл эксперимент… в Екатеринбурге… на 20% магазинов с самой высокой выручкой за предыдущий год.

Предлагаем взглянуть, как себя поведёт t-тест в подобных ситуациях на синтетических данных. Представим, что у нас магазин в среднем генерирует 200 тысяч рублей выручки за какой-то период времени со стандартным отклонением в 30 тысяч. Будем сэмплировать значения выручки для 100 магазинов из нормального распределения. Для пилотной группы установим дополнительное ограничение при сэмплировании: в неё будут попадать только магазины с выручкой выше 150 тысяч. На графике распределения пилотной и контрольной групп выглядят схоже, но, как будет показано далее, эта схожесть обманчива.

При проведении 10 000 симуляций выясняется, что статистически значимая разница средних проявляется в ~12% случаев, хотя мы установили уровень значимости в 5%. Таким образом, ошибка I рода не контролируется на заданном уровне значимости. При более экстремальных отклонениях тестовой выборки от генеральной совокупности эта ошибка легко может стремиться к 100%.

import numpy as np
from scipy.stats import ttest_ind
 
MEAN = 200000.0
STD = 30000.0
CUTOFF = 150000.0
N_SAMPLES = 100
N_SIM = 10000
ALPHA = 0.05
 
np.random.seed(1)
 
def get_biased_sample(
    n_samples: int = N_SAMPLES,
    mean: float = MEAN,
    std: float = STD,
    cutoff: float = CUTOFF,
) -> np.ndarray:
    """Сэмплируем смещенную пилотную группу."""
    test_group = []
    while len(test_group) < n_samples:
        random_sample = np.random.normal(loc=mean, scale=std)
        if random_sample > cutoff:
            test_group.append(random_sample)
            
    return np.array(test_group)
 
def show_true_error(
    n_sim: int = N_SIM,
    n_samples: int = N_SAMPLES,
    mean: float = MEAN,
    std: float = STD,
    cutoff: float = CUTOFF,
    alpha: float = ALPHA,
) -> float:
    """Замеряем эмпирическую ошибку I рода."""
    result = []
    for _ in range(n_sim):
        test_group = get_biased_sample(
            n_samples=n_samples,
            mean=mean,
            std=std,
            cutoff=cutoff,
        )
        control_group = np.random.normal(
            size=n_samples,
            loc=mean,
            scale=std,
        )
        pval = ttest_ind(test_group, control_group).pvalue
        result.append(int(pval < alpha))
 
    return np.mean(result)
 
print('Эмпирическая ошибка 1 рода:', show_true_error())

Что следует сделать? Отказаться от оценки эффекта? Можно. Отказаться от оценки всех подобных экспериментов? Фатальная ошибка. На практике пилотов, обладающих неслучайной пилотной группой, не так уж и мало. Специфика бизнеса не всегда позволяет провести рандомизированный эксперимент. Вот лишь небольшой список примеров, когда пул пилотных магазинов заведомо отличается от генеральной совокупности:

1. Пилот вносит изменения в распределительный центр, к которому привязаны определённые магазины.

2. Пилот влияет на промо-продажи, управление которыми зависит от географии магазинов.

3. Пилот физически можно реализовать только в рамках целого региона из-за ограничений, наложенных контрагентом.

4. Пилот затрагивает наружную рекламу, в рамках которой можно делать раскатку исключительно на уровне города.

Распределения групп в реальном кейсе

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

Конечно, к нам обращаются не только за оценкой постфактум, но и за дизайном эксперимента. Однако зачастую ограничения на дизайн довольно жёсткие. Даже если у заказчика есть возможность случайным образом выбирать магазины, то пул этих магазинов может быть ограничен географией (только магазины одного города/региона) и особенностями самих магазинов. Поэтому ключевой задачей на этапе оценки пилота становится подбор сопоставимой контрольной группы при условии, что генеральная совокупность уже была отфильтрована по известным неслучайным факторам. Именно здесь заканчивается A/B-тестирование и начинается Causal Inference, основной задачей которого является выявление причинно-следственных связей.

Значит ли это, что предыдущий подход не работал, и все ранее оценённые пилоты необходимо пересчитывать? Нет, просто раньше использовалась техника мэтчинга, в рамках которой к каждому пилотному магазину подбирался уникальный контрольный магазин на основе значений целевой метрики на истории (подробнее об этом в предыдущей статье). Это позволяло получить частичную сопоставимость групп, но не позволяло учесть качество сопоставимости каждого магазина при оценке эффекта. Поэтому текущие нововведения – это путь самурая к идеальным контрольной и пилотной (да-да, с ней тоже можно проводить манипуляции) группам.

1. Модель оценки эффекта

Постановка задачи

На входе мы имеем множество магазинов, на которых был проведён эксперимент, мы будем называть их пилотной группой (ПГ). Кроме этого, мы имеем множество всех потенциально контрольных магазинов (КГ), которые не участвовали в эксперименте. Это могут быть как все остальные магазины торговой сети, так и только некоторое из них, определённые заданными ограничениями, например, географией. Как мы говорили ранее, КГ и ПГ не обязательно сопоставимы. 

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

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

Общий пайплайн модели:

1. Извлекаем признаки (фичи) из данных.

2. Понижаем размерность пространства фичей и отбираем среди них признаки, коррелирующие с целевой метрикой на периоде пилота.

3. На основе полученных ранее фичей строим модель оценки вероятности попадания магазина в ПГ или КГ – propensity score. С помощью результатов работы этой модели можно будет побороть возможную несопоставимость групп.

4. Производим процедуру тримминга propensity score для лучшего сопоставления ПГ и КГ.

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

Построенная таким образом итоговая модель оценки эффекта является композицией нескольких моделей и обладает свойством Doubly Robust: оценка эффекта является корректной, если хотя бы одна из составляющих её моделей корректно специфицирована: модель оценки propensity score или модель линейной регрессии. Иначе говоря, имея корректную модель оценки propensity score, мы будем корректно оценивать эффект в нашем эксперименте, даже если вторая модель выдаёт случайный шум. И наоборот, имея корректную модель линейной регрессии, мы будем корректно оценивать эффект в нашем эксперименте, даже если модель оценки propensity score неверна.

1.1 Признаки

Прежде, чем перейти к описанию наших моделей, необходимо разобраться с тем, какие мы используем признаки. А зачем нам вообще какие-либо признаки? Все же знают, что для простого t-test нужны только две выборки, значения целевой метрики для каждого наблюдения, да в общем-то, и всё.  

В рамках же применения более сложных моделей фичи используются для двух целей:

1. уменьшение смещения КГ относительно ПГ (то, насколько наши группы отличаются в целом, даже в случае отсутствия эффекта);

2. повышение чувствительности модели оценки эффекта (также это можно назвать снижением дисперсии).

В различных источниках обычно встречаются только методы повышения чувствительности, такие как CUPED. Но, как мы помним, в большинстве наших задач ПГ может иметь сильное смещение относительно КГ. Отметим, что задачи повышения чувствительности и уменьшения смещения должны решаться одновременно, так как иначе мы можем повышать чувствительность, не уменьшая смещение. В этом случае доверительные интервалы будут строиться с центром в ложно оценённом эффекте, что чревато высокой ошибкой I рода. Соответственно, нужны фичи, которые влияют:

1. на целевую метрику;

2. на статус присвоения к ПГ и КГ.

Целевую переменную мы будем называет таргетом (), а статус присвоения ПГ или КГ – тритментом (). В литературе фичи, которые удовлетворяют только первому свойству, называются ковариатами, только второму – инструментами, двум свойствам сразу – конфаундерами. 

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

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

Агрегация и разности

Наши данные представимы в виде набора временных рядов. Каждый временной ряд описывает динамику изменения целевой метрики (например, выручки) для одного конкретного магазина в течение последних 2-3 лет. Всего временных рядов столько же, сколько и рассматривается магазинов.

Для того, чтобы далее работать с выборкой независимых объектов, мы агрегируем данные за пилотный период до одной точки для каждого магазина. Например, на картинке ниже пилотный период имеет длину 2. Таким образом, исходный ряд преобразуется в ряд , где – значение целевой метрики за весь период пилота, – её значение за период (-1, 0) и т. д. 

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

Как учесть период раскатки

Если имеется период раскатки, то нужно его как-то учесть. Мы хотим, чтобы из значения на периоде пилота вычиталось первое (ближайшее) значение препилотного периода, то есть . На картинке выше, например, при взятии таких разностей мы сделали 4 шага назад (один период пилота + длина раскатки, равная трём длинам пилота), поэтому и т. д. 

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

Аналогия с подходом diff-in-diff

В этом процессе взятия разностей можно заметить некоторый аналог известного подхода diff-in-diff. Если мы сравним группы на основе рядов разностей целевой метрики по пилотному периоду "как есть", то мы получим разность разностей . 

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

Другие варианты фичей:

1. Фичи, определяющие тренд ряда разностей. Для их определения по ряду разностей на препилотном периоде для каждого магазина отдельно обучается полином методом МНК для степеней 0, 1 и 2. Все полученные коэффициенты из каждого полинома берутся в качестве признаков.

2. Фичи на основе повторных разностей разных степеней. По сути, функция применяется два раза для различных . Например: .

3. Можно также использовать специальные ковариаты, заданные пользователем. Они могут быть как вещественными, так и категориальными, причём, во втором случае можно применять различные алгоритмы энкодинга. Пример вещественной ковариаты – количество завезённого в магазин товара, пример категориальной – концепт магазина, уровень ассортимента.

4. Фичи с особым значением для каждого магазина:

   • среднее разностной метрики по всему препилоту;

   • последнее и предпоследнее значения разностной метрики на препилоте;

   • значение разностной метрики в те же даты пилота, но в предыдущем году; для учёта сезонных факторов.

1.2 Уменьшение размерности и отбор признаков

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

Для преобразования фичей сперва мы используем метод главных компонент (PCA). Данный метод позволяет сократить количество фичей, выполнив линейный переход из одного пространства фичей в другое, в котором новые фичи будут некоррелированными (на самом деле, гарантируется только равенство нулю оценки коэффициента корреляции). Однако, некоторые фичи, получаемые после применения PCA, могут не оказывать практически никакого влияния на и/или на . Такие фичи могут увеличить дисперсию оценок в линейных моделях, поэтому подобные фичи необходимо удалить.

После PCA мы можем выделить фичи, которые коррелируют только с или одновременно с и . Напомним, фичи, влияющие на (ковариаты), помогают повысить чувствительность метода (аналог CUPED), а фичи, влияющие ещё и на (конфаундеры), помогают снизить смещение ПГ относительно КГ. 

Почему бы тогда не оставить все фичи? Дело в том, что лишние фичи, которые не оказывают никакого влияния на или на , наоборот, могут повысить дисперсию оценки эффекта. В процессе отбора фичей мы оставляем только те фичи, которые имеют статистически значимую корреляцию только с или одновременно с и с .  Кроме того, из всех таких фичей со статистически значимой корреляцией мы берём только топ-N фичей по абсолютной величине корреляции. Так как фичи посчитаны до пилота, то причина корреляции может быть только однонаправленной: фичи влияют на , а не наоборот. Аналогичным образом фичи влияют на . 

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

В конечном итоге остаются фичи следующего вида:

• Фичи, преобразованные в рамках PCA и в дальнейшем отобранные по значениям статистически значимых коэффициентов корреляции.

• Пользовательские ковариаты, если они не участвуют в PCA.

• Все фичи с особым значением, не участвующие в PCA и не проходящие процедуру отбора.

1.3 Propensity Score

Введём понятие propensity score. Говоря неформально, propensity score – это некоторая величина, с помощью которой можно перевзвешивать наблюдения в выборке или же разделять выборку на блоки с целью избавления от изначального смещения ПГ относительно КГ, вызванного неслучайным назначением . Если же говорить формально, propensity score – это условная вероятность назначения конкретному объекту при условии его заранее отобранных фичей : .

Можно понимать propensity score как функцию, которая преобразует вектор фичей в единственное число, в котором сосредоточена вся информация о том, как признаки влияют на статус присвоения , все же остальные факторы в совокупности – непрогнозируемый шум. Было бы логично сравнивать между собой только те магазины, которые имеют примерно одинаковое значение propensity score, тем самым считая эффект среди сопоставимых магазинов. Иначе говоря, это те магазины, которые с одной и той же вероятностью могут быть отнесены к ПГ, причём, фактическая реализация статуса присвоения – чистая случайность.

Построение с помощью логистической регрессии

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

где – оценка коэффициентов в модели логистической регрессии.

Основные идеи использования

Propensity Score Matching, PSM. Для каждого объекта из ПГ найти наиболее похожий по propensity score на него объект из потенциальной контрольной группы и включить в КГ. В полученной КГ объекты могут повторяться. Далее напрямую сравниваем ПГ с подобранной КГ.

Propensity Score Blocking, PSB. Разделить объекты в выборках на основе значений оценок propensity score на блоки так, чтобы внутри каждого блока были бы достаточно сопоставимыми. Далее оценить эффект в каждом блоке независимо с помощью некоторой статистической процедуры, а затем определённым образом объединить результаты.

Propensity Score Weighting, PSW. Взвесить наблюдения так, чтобы получить выборки сопоставимых магазинов. По результатам наших экспериментов, пайплайн которых приведён далее, наилучшие результаты показывает данный подход со взвешиванием объектов, поэтому его мы и рассмотрим подробнее.

Калибровка

При оценке propensity score может возникнуть соблазн использовать различные сложные техники машинного обучения. Однако, чаще всего это неверный путь, поскольку в данном случае мы стремимся не к высокой точности прогнозирования , а к несмещённым оценкам вероятностей для получения более точной сопоставимости выборок. Как мы сказали ранее, именно это требование обеспечивает выполнение свойства Doubly Robust нашей модели. В частности, очень уверенная модель, верно классифицирующая выборку, и предсказывающая вероятности в районе нуля или единицы, не позволяет использовать такие оценки вероятностей для сопоставления выборок. Учитывая эти особенности, для получения несмещённых прогнозов вероятностей мы используем калибровку (Calibrated Classifier CV) предсказаний вероятностей, о которой можно почитать в этой статье. Кроме того, не стоит использовать сильную регуляризацию в модели логистической регрессии при оценке propensity score.

1.4 Тримминг

Для получения несмещённых оценок эффекта необходимо выполнение условия на propensity score, которое носит название Common Overlap (Common Support, Positivity Check). Это условие означает, что область пересечения значений propensity score должна быть существенной для ПГ и КГ. Например, на представленном ниже графике слева для выделенных магазинов КГ почти нет сопоставимых магазинов из ПГ:

Для выполнения условия Common Overlap мы можем очистить наблюдения-выбросы, тем самым убив сразу двух зайцев:

1. будут выброшены магазины с сильной дисперсией, следовательно, оценка эффекта будет иметь меньшую дисперсию;

2. будут выброшены несопоставимые магазины, поэтому смещение в оценке будет меньше.

Выбросы можно удалять по минимальным и максимальным значениям каждой из групп и по квантилям. Отметим, что для достижения сопоставимости групп важно применять одинаковые правила очистки к обеим группам. Например, можно посчитать квантили только по пилотной группе, но удалить по ним выбросы необходимо из обеих групп сразу. В зависимости от размера групп можно производить несколько раундов тримминга с разными правилами, каждый раз заново оценивая propensity score.

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

1.5 Propensity Score Weighting

Регрессия

В случае обычного t-test мы сравниваем средние, посчитанные по ПГ и КГ

Если эти средние отклоняются друг от друга больше, чем на некоторую величину, зависящую от их дисперсий, то мы говорим, что наш t-test детектирует статистически значимое отличие средних друг от друга:

где – уровень значимости, – квантиль нормального распределения, – оценки дисперсий в ПГ и КГ соответственно.

На самом деле t-test легко получить из обычной линейной регрессии, которая в данном случае будет выглядеть, как

где – бинарный признак (равен нулю или единице), означающий статус присвоения к ПГ и КГ.

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

Можно усовершенствовать полученную регрессию, добавив в неё ковариаты и конфаундеры

где – вектор ковариат и конфаундеров. Такая модель является аналогом CUPED и позволяет объяснить некоторую часть случайного шума, изначально заложенного в , тем самым уменьшив дисперсию оценки коэффициента . Об этом можно подробнее почитать в этой статье. Кроме того, в рамках этой модели можно полностью избавиться от возможного смещение ПГ относительно КГ, но только при включении в модель всех ковариат и конфаундеров. Однако, это довольно проблематично реализовать на практике, поэтому для получения более устойчивых к смещению оценок мы дополнительно используем модели на основе propensity score, о чём будет более подробно рассказано далее.

Определение весов

Следующий этап усовершенствования нашей модели – введение весов для каждого объекта, которые мы будем использовать в регрессии. Рассмотрим магазин и пусть – оценка его propensity score, которую мы уже получили ранее. Напомним, что – оценка вероятности того, что магазин будет отнесён к ПГ. 

Пусть магазин находится в ПГ, и достаточно мала. Это может свидетельствовать о том, что подобные магазины в большинстве своём попадают в КГ, и магазин выглядит довольно интересным для сравнения с ними, следовательно, его вес необходимо увеличить. Сопоставим ему вес . То есть, чем меньше , тем больше вес. И наоборот, чем больше , тем больше подобных магазинов имеется в ПГ, и тем их вес меньше. Если же магазин находится в КГ, то ситуация обратная, мы определяем его вес как . Отметим также, что ранее мы потребовали несмещённость оценок propensity score, тем самым действительно можно представлять как долю "примерно одинаковых магазинов", которые попали в ПГ.

Обобщая два случая, можно получить формулу весов

где – равен 1 для ПГ и 0 для КГ.

Рассмотрим простой пример, поясняющий принцип работы весов. Представим, что был проведён пилот по внедрению новой умной системы скидок, в котором ПГ состоит из магазинов, различающихся своим ассортиментом: например, 8 магазинов, где можно купить розовых единорогов, и 1 магазин, где можно купить упитанных бегемотов. При этом для КГ осталось 2 магазина с единорогами и 5 магазинов с бегемотами. Для простоты примера будем считать, что в остальном магазины идентичны, соответственно, в качестве признаков будем использовать только ассортимент магазинов. В таблице ниже приведён расчёт весов для данного пилота:

Из последней колонки таблицы следует, что после взвешивания мы получили сопоставимые ПГ и КГ. Ниже приведена соответствующая визуализация. Логистическая регрессия, по сути, обобщает данный подход на случай произвольных типов признаков, в частности, для вещественных признаков.

Оценка эффекта и стандартных ошибок 

В случае прямого аналога t-test мы должны сравнить между собой взвешенные средние целевой переменной по ПГ и КГ, то есть величины . Отметим, что если ПГ и КГ выбраны случайно, то есть каждый магазин попадает в ПГ с вероятностью 1/2, мы получаем обычный t-test. 

Иногда PSW в простом случае записывают в виде одной суммы по всем магазинам , где веса рассматриваются со знаком: .

Для учёта использования ковариат и конфаундеров будем использовать линейную регрессию вида

где – вектор ковариат и конфаундеров. Предварительно необходимо все признаки , константный признак, а также целевую величину (то есть все, кроме ) домножить на положительный вес , при этом веса необходимо предварительно отнормировать. 

Признак домножать на вес не нужно, поскольку мы хотим, чтобы разница в предсказании при изменении с 0 на 1 и фиксировании остальных признаков была бы равна эффекту, который мы хотим оценить. Иначе говоря, в такой модели равен эффекту от нашего пилота, если для КГ и для ПГ.

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

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

Построить доверительный интервал можно стандартным способом

где – уровень значимости, – квантиль распределения Стьюдента. Полученная оценка как раз является той оценкой эффекта, которую мы искали. Её и соответствующий доверительный интервал мы выдаём заказчику. Результат называем статистически значимым, если этот интервал не содержит 0, и в таком случае мы можем говорить, что эффект от инициативы присутствует. Иначе, если интервал содержит 0, мы можем заключить, что на основе имеющихся данных нет оснований полагать наличие эффекта.

2. Валидация модели на этапе оценки пилота

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

Что мы должны cделать? Во-первых, получить от заказчика все необходимые вводные, например, метрики, величину ожидаемого эффекта, ограничения на географию эксперимента, даты проведения пилота, наличие периода раскатки и т. д. Здесь полезно сразу учесть потенциальные источники смещения, не всё отдавая на откуп модели. Например, если заказчик проводит эксперимент только на магазинах с пекарнями, то имеет смысл подбирать КГ только среди всех остальных магазинов с пекарнями. Во-вторых, нам необходимо выбрать наиболее оптимальную модель, о чём мы и будем говорить в данном разделе. Наконец, мы отдаём заказчику отчёт, в котором зафиксирован полученный эффект (или же его отсутствие) — точечная оценка и доверительный интервал.

Для выбора оптимальной модели мы используем традиционные статистические показатели – ошибку I рода для статистических критериев, а также минимальный детектируемый эффект для заданной вероятности ошибки II рода. Ошибка I рода – случай детектирования эффекта нашей моделью, когда его на самом деле нет. Ошибка II рода – случай отсутствия детектирования эффекта, когда он на самом деле есть. Соответственно, можно определить вероятности этих ошибок, причём, для ошибки II рода соответствующая вероятность определяется для конкретной величины эффекта.

Зависимость вероятности ошибки I рода от распределения

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

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

Соответственно, мы ограничиваем вероятность ошибки I рода на уровне 20%, иначе говоря, мы допускаем возможность получить ложноположительный результат в одном из пяти случаев. В литературе чаще встречается ограничение на уровне 5%, но в силу специфики данных и характера бизнес-гипотез мы используем менее строгий порог. Далее мы ещё подробнее рассмотрим процесс оценки вероятностей ошибок I и II рода.

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

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

Гиперпараметры

Единственный гиперпараметр, который мы перебираем в процессе валидации методов –  это теоретический уровень значимости критерия. Здесь у некоторых читателей может возникнуть непонимание: как же его можно перебирать, его ведь нужно фиксировать до эксперимента! Вообще говоря, это верно. Особенность нашей конкретной ситуации заключается в том, что статистические свойства нашей модели справедливы при определённых теоретических предположениях, таких как достаточно большие размеры выборок и их репрезентативность (которой чаще всего добиваются рандомизированным выбором ПГ и КГ). Поэтому мы и называем соответствующий уровень значимости теоретическим. Однако, на практике данные могут не соответствовать этим предположениям, и реальная вероятность ошибки I рода оказывается выше ожидаемой. Например, если мы фиксируем теоретический уровень значимости на уровне 5%, можем получить оценку вероятности ошибки I рода на уровне 30%.

На самом деле, доверительный интервал (соответственно, и статистический критерий) построен так, что при уменьшении теоретического уровня значимости (увеличении уровня доверия интервала) доверительный интервал будет расширяться, следовательно, отвержение будет происходить реже, и вероятность ошибки I рода уменьшится. Тем самым, перебрав некоторое множество значений теоретических уровней значимости, мы можем найти то значение, для которого оценка вероятности ошибки I рода приближенно равна желаемому значению 20%.

Оценка вероятностей ошибок (AA-тестирование)

Для более точных оценок вероятностей ошибок, соответствующих распределению наших данных, мы должны приблизить данное распределение. Для этого мы пользуемся сэмплированием по методу Монте-Карло из предполагаемого распределения наших данных. Но где же нам найти достаточно данных для оценки этого распределения, подобных тем, на которых мы применяем нашу модель? На самом деле – нигде. 

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

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

1. выбираем псевдопилотный период на исторических данных;

2. для данного псевдопилотного периода применяем нашу модель;

3. замеряем эффект на псевдопилотном периоде (он либо есть, либо его нет).

В конечном итоге у нас есть бинарных флагов наличия/отсутствия эффекта. В качестве оценки вероятности ошибки I рода возьмём долю случаев, когда эффект был найден. Соответственно, оценка вероятности ошибки II рода строится похожим образом, но до применения модели нужно добавить искусственный эффект к ПГ на псевдопилотном периоде. Далее необходимо взять долю тех случаев, в которых эффект не обнаружен:

3. Дизайн пилота

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

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

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

Рандомизированная пилотная группа

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

Например, Вася решил ограничить пилотную группу магазинами Москвы, но при этом он может себе позволить взять в пилот столько магазинов, сколько мы сочтём нужным, а ожидаемого эффекта у него нет. Соответственно, подбор контрольной группы также осуществляется из Москвы. В процессе такого дизайна мы варьируем число магазинов от некоторого заданного снизу порога (например, 30 магазинов) до порядка половины от всех магазинов в Москве. 

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

Приведём пошаговое описание процедуры дизайна:

1. Выбираем размер ПГ (перебираем из заранее заданного набора).

2. Для оценки вероятности ошибки I рода и MDE повторяем раз следующие шаги: 

   • выбор псевдопилотного периода на исторических данных;

   • выбор случайной пилотной группы, оставшиеся магазины определяются в потенциальную контрольную группу;

   • для данного псевдопилотного периода применяем нашу модель;

   • замеряем эффект на псевдопилотном периоде (он либо есть, либо его нет).

3. На основе полученных расчётов оцениваем вероятность ошибки I рода и MDE.

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

Внимательный читатель заметит, что подобной процедурой мы фактически оцениваем минимальный детектируемый эффект (MDE), для которого существует точная формула 

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

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

Дизайн при отсутствии возможности рандомизации

В очередной раз Вася пришёл к нам с проблемой следующего вида: он хочет получить от нас дизайн пилота, но у него есть существенное ограничение – невозможна рандомизация пилотной группы. Например, из-за ограничений различных бизнес-процессов раскатка возможна только на все магазины одного или нескольких распределительных центров (РЦ). Это может быть тестирование нового ассортимента товаров, или же нового подхода к доставке товаров из РЦ в соответствующие ему магазины. 

В данном случае вместо варьирования размера пилотной группы нам необходимо перебрать различные пилотные группы, определяемые конкретными РЦ. Фактически, мы выбираем, какой РЦ будет выступать в качестве пилотного, работая с магазинами данного РЦ:

В случае, при котором пилот можно раскатить только на один из всех возможных РЦ, дизайн происходит следующим образом:

1. Выбираем конкретный РЦ, все соответствующие ему магазины относятся к ПГ, все остальные определяются в потенциальную КГ.

2. Повторяем раз следующие шаги:

   • выбор псевдопилотного периода по историческим данным;

   • для данного псевдопилотного периода применяем нашу модель;

   • замеряем эффект на псевдопилотном периоде (он либо есть, либо его нет).

3. Оцениваем вероятности ошибки I рода и MDE на основе полученных результатов.

4. Оценка минимального детектируемого эффекта

Одной из интересных проблем является оценка минимального детектируемого эффекта (MDE), то есть эффекта, который наша модель способна обнаружить при заданных вероятностях ошибок I и II рода и размере пилотной группы. Из предыдущего раздела стало понятно, почему аналитическая формула не работает. Однако без MDE дизайн эксперимента не сделаешь, а волшебного алгоритма для быстрого подсчёта MDE эмпирическим путем с неба не свалится. 

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

Размеры этих эффектов мы выявили эмпирически для каждой метрики, проведя множество пилотов. Таким образом, при удовлетворительных значениях оценок вероятностей ошибок I рода и фиксированном размере пилотной группы мы ищем тот эффект, который соответствует принятому порогу вероятности ошибки II рода. На выходе имеем MDE для каждого размера пилотной группы.

Без чекушки графика тут не разобраться, поэтому далее приведём пример. Представим, что мы зафиксировали размер пилотной группы, которую будем сэмплировать из некоторого подмножества магазинов. А также зафиксируем порог вероятности ошибки II рода в 20% (соответствует стандартному значению мощности в 80%). 

Далее на наши реальные данные набрасываем синтетические эффекты: 0.1%, 1%, 3%, 5%, 7%, 10%. Добросовестно оцениваем вероятности ошибки II рода для каждого из этих значений. После этого между точками проводим прямые линии и ищем такое значение на оси X, которое бы соответствовало выбранному порогу в 0.2 на оси Y. Найденное значение – 4.67% – и будет нашим MDE для заданного размера пилотной группы. Далее проделаем аналогичные действия для каждого размера пилотной группы и получим ряд значений MDE:

На этом инженерная мысль себя не исчерпывает. Ведь после получения MDE для ограниченного числа размеров пилотной группы хотелось бы получить более точную оценку количества необходимых магазинов. Например, мы оценили MDE для 50 магазинов и 100 магазинов, как 5% и 2% соответственно. А заказчик говорит, что рассчитывает увидеть эффект в 3%. 

Неужели придётся перебирать все возможные размеры пилотной группы? Не придётся, если вспомнить тот факт, что теоретическая зависимость между размером выборки и MDE описывается формулой . Остаётся только подобрать такое значение , при котором минимизируется MSE, то есть решить задачу оптимизации. Таким образом, мы избегаем перебора размеров пилотных групп и получаем вполне точные значения MDE для любого размера пилотной группы. График ниже демонстрирует, как это работает на практике:

5. Выбор наилучшей модели

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

Кроме того, при валидации на данных, соответствующих конкретному эксперименту, у нас зачастую нет информации о том, насколько ПГ смещена относительно КГ. В рамках отдельного пайплайна исследований мы можем контролировать смещение между группами.

Постановка задачи

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

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

Валидация на синтетических и реальных дополненных данных

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

Валидация происходит в два этапа: на синтетических данных и на модификациях реальных (фактических) данных. На каждом этапе производится итераций следующего вида:

1. подготавливается датасет только с помощью генерации случайных чисел или с помощью модификаций реальных данных (этот процесс отдельно для каждого типа расписан далее);

2. для каждого исследуемого метода оценивается эффект и строится доверительный интервал.

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

• размер эффекта (его мы задаём сами);

• размер выборки;

• доля пилотной группы среди всего набора данных;

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

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

В отличие от синтетических данных, в реальных дополненных данных есть меняющаяся временная структура. Например, в какие-то непоследовательные периоды времени динамика среднего в ПГ и КГ может расходиться, дисперсия может изменяться во времени и т. д. При этом варьируются следующие параметры конфигурации датасета:

• пилотный период, который выбирается из временного интервала, представленного в данных;

• количество недель пилота;

• размер эффекта (его мы задаем сами);

• размер выборки;

• доля пилотной группы среди всего набора данных;

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

Определение свойств алгоритма и выбор наилучшего из алгоритмов

Результаты анализируются по такой схеме:

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

2. Аналогичным образом исключаются методы, для которых оценка вероятности ошибки I рода по реальным дополненным данным превосходит некоторый заранее фиксированный порог. 

3. Для оставшихся моделей вычисляется оценка смещения на основе сравнения оценок эффекта по проведённым итерациям с истинными размерами эффектов, которые мы искусственно добавляли к данным. Из этих моделей выбирается множество моделей, которые не отличаются статистически значимо от наилучшей по оценке смещения.

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

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

6. Куда можно двигаться дальше?

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

Например, мы можем построить модель, которая будет обучаться на значениях целевой метрики пилотной группы на периоде до пилота, а после предсказывать значения этой метрики уже на периоде пилота. Так мы получим альтернативную реальность, где нашей инициативы не было, и сможем с ней сравниться. Конечно, такой подход наивен и не учитывает влияния каких-то глобальных трендов. Исправить это можно добавлением ковариат, помогающих объяснить поведение целевой метрики. Именно такой подход реализован в библиотеке для R и Python от Google. Он позволяет оценить эффект в тех случаях, когда кандидаты в контрольную группу отсутствуют совсем.

Если кандидаты в контроль всё же имеются, то альтернативную реальность можно создать путём их взвешенной суммы, как это происходит в алгоритме под названием Synthetic Difference In Differences. 

Придание весов отдалённо напоминает оценку propensity score, так как аналогично позволяет учесть «важность» магазина при оценке эффекта. Однако тут накладывается ограничение: сумма весов должна быть равна единице. На примере наших экспериментов это означает, что синтетический магазин может состоять на 20% из одного реального магазина и на 80% из другого реального магазина, но сумма не должна превышать 100%. 

Однако, не только контрольные магазины получают свои веса, но и временные периоды, так как каждый из них вносит неравнозначный вклад в оценку эффекта. За готовой имплементацией и простым объяснением того, что происходит под капотом данной модели, стоит обратиться к одной из глав учебника по Causal Inference. Этот учебник мы рекомендуем в принципе всем, кто хочет за короткий срок погрузиться в тему (в учебнике также есть много примеров на Python).

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

Мир Causal Inference не стоит на месте и старается идти в ногу с хайпом машинного обучения и искусственного интеллекта. Вот наглядный пример – нейросеть DragonNet, которая используется для оценки Propensity Score. И снова авторы утверждают о SOTA-результатах по сравнению с предыдущими подходами. 

В итоге практикующим специалистам остаётся только экспериментировать со всеми этими наработками в поисках наиболее оптимального подхода.

7. Заключение

Подведём итоги:

1. Ввиду различных бизнес-ограничений, зачастую невозможно провести стандартный A/B-тест, и в таком случае можно использовать методы Causal Inference.

2. Как показывает практика, для оценки нестандартных экспериментов неплохо подходит модель Propensity Score Weighting, обладающая свойством двойной устойчивости (Doubly Robust).

3. Также важным является выстраивание корректных пайплайнов дизайна и оценки экспериментов с учётом специфики поставленных задач – прежде всего, необходимо обеспечить валидацию модели на исторических данных.

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

Благодарим всех читателей за внимание! Если у вас есть вопросы, будем рады на них ответить в комментариях или личных сообщениях.

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

Наша маленькая команда по разработке методологии и инструментов для A/B-тестирования:

• Валл Федор

• Волков Никита [tg: @ nikitavolkov21]

• Лавров Денис

• Морозов Станислав