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

1. Предыстория

Начнем с того, что в любом магазине товары делятся по категориям разного уровня. Например, на верхнем уровне могут быть такие категории, как «Продукты», «Свежее» и «Непродовольственные товары». В свою очередь, категория «Продукты» может разделяться на подкатегории, такие как «Молоко», «Мясо» и «Рыба». Категория «Мясо» может быть еще более детализирована, включая подкатегории «Охлажденное мясо», «Замороженное мясо» и т. д.

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

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

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

• Можно отловить статистически значимый эффект на тотал-категории, даже если его не получилось найти для отдельных категорий.

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

А вот и нет…

Точнее, в реальном мире так и будет. Однако в статистике так будет не всегда. И в этой статье мы выясним, где деньги, как такое возможно и что с этим делать.

2. Согласованность и напоминание о методах A/B-тестирования

Начнем с самого любимого — дадим определение.

Определение. Пусть – оценка эффекта на категорию, – оценка эффекта на тотал-категорию. Будем называть оценки согласованными, если выполняется равенство

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

Чтобы понять, в каких случаях оценки эффектов перестают быть согласованными, вспомним методику А/B-тестирования на магазинах.

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

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

В таких случаях простой t-тест не может контролировать вероятность ошибки первого рода на заданном уровне значимости, поскольку выборки являются зависимыми из-за неслучайного выбора. В упомянутой выше статье мы описывали метод, основанный на композиции propensity score weighting (PSW) и линейной регрессии (многомерный CUPED), обладающий свойством Doubly Robust (двойная устойчивость) – устойчивости модели к неверным предположениям о любой из этих двух моделей.

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

Подробнее об этом можно прочитать в нашей статье о методике А/B-тестирования на магазинах.

2.1 Но где здесь проблема с согласованностью?

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

2.1.1 T-test

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

где и – размеры групп. Доверительный интервал для этого эффекта имеет вид:

где — выборочные дисперсии КГ и ПГ соответственно, а — квантиль распределения Стьюдента.

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

Таким образом, применение t-критерия Стьюдента не нарушает согласованность оценок эффектов.

2.1.2 CUPED

Метод CUPED использует ковариаты (признаки, которые могут влиять на целевую метрику) для снижения дисперсии, что повышает чувствительность теста. Напомним, что при наличии ковариат можно обучить линейную регрессию следующего вида:

где — данные по обеим выборкам, — метка группы (0 для КГ, 1 для ПГ), а — вектор признаков-ковариат. Тогда коэффициент будет равен эффекту в нашем А/B-тесте.

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

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

Давайте рассмотрим пример кода, чтобы проиллюстрировать эту проблему. Этот код выполняет оценку эффектов четырьмя разными методами (t-test, CUPED, PSW и Doubly Robust):

from typing import List, Tuple

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.calibration import CalibratedClassifierCV
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler


class TreatmentEffectModel:
    def __init__(
        self,
        df: pd.DataFrame,
        treatment_name: str,
        target_name: str,
        covariate_list: List[str] = None,
        prop_score_features: List[str] = None,
    ) -> None:
        """Инициализация модели оценки эффекта тритмента.

        Параметры
        ----------
        df : pandas DataFrame
            DataFrame с данными, содержащий метку тритмента и целевую переменную.
        treatment_name : str
            Название столбца, содержащего метку тритмента.
        target_name : str
            Название столбца, содержащего целевую переменную.
        covariate_list : List[str], по умолчанию=None
            Список дополнительных ковариат.
        prop_score_features : List[str], по умолчанию=None
            Список признаков, используемых для оценки propensity score.
        """
        self.df = df
        self.treatment_name = treatment_name
        self.target_name = target_name
        self.covariate_list = covariate_list
        self.prop_score_features = prop_score_features
        self.result_model = None

    def predict_propensity_score(self) -> np.ndarray:
        """Оценивает propensity score для заданного тритмента.

        Возвращает
        -------
        np.ndarray
            Массив, содержащий оцененные значения propensity score для каждого наблюдения.
        """
        features_matrix = self.df[self.prop_score_features]
        treatment_col = self.df[self.treatment_name]
        scaler = StandardScaler()
        scaled_features_matrix = scaler.fit_transform(features_matrix)

        ps_model = CalibratedClassifierCV(LogisticRegression(C=1e6), method="sigmoid", cv=3)
        ps_model.fit(scaled_features_matrix, treatment_col)
        propensity_score = ps_model.predict_proba(scaled_features_matrix)[:, 1]

        return propensity_score

    def calculate_weights(self) -> np.ndarray:
        """Вычисляет веса для анализа.

        Возвращает
        -------
        np.ndarray
            Массив весов для каждого наблюдения.
        """
        if self.prop_score_features is not None:
            prop_score = self.predict_propensity_score()
            weights = (self.df[self.treatment_name] - prop_score) / (prop_score * (1 - prop_score))
            weights = np.sqrt(np.abs(weights.values))
            weights /= weights.mean()
        else:
            weights = np.ones(len(self.df))

        return weights

    def predict_effect(self) -> None:
        """Предсказывает эффект с помощью линейной регрессии."""
        if self.covariate_list is not None:
            features = self.df[[self.treatment_name] + self.covariate_list]
        else:
            features = self.df[[self.treatment_name]]
        weights = self.calculate_weights()

        weighted_target = self.df[self.target_name] * weights
        weighted_features = features * weights.reshape((-1, 1))

        self.result_model = sm.OLS(weighted_target, weighted_features).fit(cov_type="HC0")

    def get_summary(self) -> Tuple[float, Tuple[float, float], float]:
        """Возвращает сводку результатов модели.

        Возвращает
        -------
        Tuple[float, Tuple[float, float], float]
            Средний эффект тритмента (ATE), доверительный интервал и p-значение.
        """
        stat_table = self.result_model.summary().tables[1]

        ate = float(stat_table.data[1][1])
        confidence_interval = (float(stat_table.data[1][5]), float(stat_table.data[1][6]))
        pvalue = float(stat_table.data[1][4])

        return {"ate": ate, "confidence interval": confidence_interval, "pvalue": pvalue}


def calculate_categories_effect(
    df: pd.DataFrame,
    categories: List[str],
    mode: str = "one",
    treatment_name: str = "is_pilot",
    covariate_list: List[str] = None,
    prop_score_weighting: bool = True,
) -> List[Tuple[float, Tuple[float, float], float]]:
    """Выполняет анализ эффектов по категориям.

    Параметры
    ----------
    df : pandas DataFrame
        DataFrame с данными, содержащий метку тритмента и целевую переменную.
    categories : List[str]
        Список категорий для анализа.
    mode : str, по умолчанию='one'
        Режим выбора признаков для оценки propensity score.
    treatment_name : str, по умолчанию='is_pilot'
        Название столбца, содержащего метку тритмента.
    covariate_list : List[str], по умолчанию=None
        Список ковариат, которые могут быть включены в модель.
    prop_score_weighting: bool, по умолчанию=True
        Указывает, следует ли использовать взвешивание на основе propensity score.

    Возвращает
    -------
    List[Tuple[float, Tuple[float, float], float]]
        Список эффектов по каждой категории, включая доверительные интервалы и p-значения.
    """
    summary_list = []
    summary_table = {}

    for category in categories:
        target_name = f"rto_{category}"

        if prop_score_weighting:
            prop_score_features = (
                [f"prepilot_rto_{category}"] if mode == "one" else ["prepilot_rto_" + cat for cat in categories[:-1]]
            )
        else:
            prop_score_features = None

        model = TreatmentEffectModel(
            df,
            treatment_name,
            target_name,
            covariate_list,
            prop_score_features,
        )
        model.predict_effect()
        summary = model.get_summary()
        summary_list.append({f"ATE_{category}": summary})
        summary_table[
            f"Category_{category}"
        ] = f"{summary['ate']} ({summary['confidence interval'][0]}, {summary['confidence interval'][1]})"

    print(pd.DataFrame([summary_table]))

    return summary_list

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

В этой таблице мы видим несогласованность эффектов, например, в последней строчке 18.9 + 218.3 + 694.0 = 931.2 ≠ 1161.5, значения суммарного эффекта различаются!

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

2.1.3 Propensity Score Weighting

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

Формально propensity score — это условная вероятность назначения объекту при условии его признаков :

Эту вероятность мы не знаем, но можем оценить ее с помощью логистической регрессии, используя влияние признаков на назначение :

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

Далее можно применить t-тест с взвешенными средними. Более продвинутый метод — doubly robust, который использует взвешенную линейную регрессию:

а в качестве весов используются .

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

Результаты предсказания эффектов в методе PSW следующие:

Видно, что, как и в случае с CUPED, сумма эффектов по категориям не равна эффекту для тотал-категории.

3. Решение проблемы

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

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

3.1 Случай CUPED

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

где – тотал-метрика, – метрика на категории где , а — общий вектор признаков-ковариат. При этом известно, что

Тогда мы имеем следующую оценку коэффициента для тотал-категории

и для -ой категории

где матрица и векторы , имеют вид

Проверим согласованность оценок

Следовательно, равенство верно и для второй компоненты вектора

Таким образом, при использовании одинакового набора признаков в методе CUPED мы достигаем согласованности оценок эффектов.

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

3.2 Случай PSW

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

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

Важно отметить, что в модели doubly robust согласованность гарантируется только в том случае, если обе модели (propensity score и линейная регрессия) обучаются на наборах признаков, которые не различаются между экспериментами. Если это условие выполняется только для одной из моделей, то согласованность оценок эффектов не достигается, что будет продемонстрировано далее на примерах.

4. Эксперименты

Мы провели эксперименты на синтетических и реальных данных. Синтетические данные для значений метрик во время эксперимента и до него генерировались в двух форматах: для сопоставимых (data=unbiased) и несопоставимых выборок (data=biased). Кратко опишем процесс генерации синтетических данных:

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

• Для каждой категории товаров вычисляются параметры для гамма-распределения на основе вычисленных ранее средних значений и стандартных отклонений.

• Генерируются данные продаж согласно гамма-распределению с вычисленными параметрами и добавляется нормальный шум.

• Для пилотного периода может добавляться случайный эффект от проведенного эксперимента.

• Вычисляются значения тотал-метрики путем суммирования значений метрик по категориям.

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

• Простой t-test

• CUPED

• PSW

• Doubly Robust

Для каждого метода, кроме t-test, рассматривались два режима работы:

• В качестве ковариаты модель использует данные за предыдущей период только по той же метрике, эффект для которой мы и оцениваем (mode=one).

• Модель использует данные по всем метрикам (mode=all).

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

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

Кроме того, мы провели эксперименты по оценке вероятности ошибки первого рода и мощности для наших тестов, сравнивая следующие случаи для всех метрик (как по отдельным категориям, так и для тотал-метрики):

• В качестве признаков используются только значения рассматриваемой метрики на исторических данных.

• В качестве признаков используются значения всех метрик на исторических данных.

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

5. Итог

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

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

Отдельно хотелось бы выразить благодарность нашим коллегам, которые на разных этапах оказывали помощь в исследовании и написании статьи, в особенности Семену Дипнеру, Анастасии Майстер, Станиславу Морозову.

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

• Никита Волков

• Мария Милованова

• Диана Наниева