Линейная регрессия на стероидах: Double Machine Learning для устранения смещений в данных

Любой аналитик знает, что самым надёжным способом проверки гипотез являются рандомизированные контролируемые эксперименты (RCT), или, как их называют в народе — A/B-тесты. На практике часто возникают ситуации, когда провести A/B-тест невозможно — в основном это происходит по этическим или техническим причинам. Однако бывают кейсы, когда рандомизация невозможна потому, что treatment-ом является определённое действие пользователя. Например, treatment-ом может быть оформление платной подписки или отмена бронирования на сервисе. Давайте назовём такой вид воздействия добровольным.

В русскоязычном пространстве, и в частности на Хабре, достаточно много статей, посвящённых таким методам Causal Inference, как DiD, PSM и Causal Impact. Тем не менее, к моему удивлению, практически нет статей, посвящённых методам на основе ортогонализации и regression adjustment, хотя, на мой взгляд, именно эти методы являются самыми удобными для оценки эффекта от добровольного treatment-а. Пришло время исправить это недоразумение и разобрать метод Double/Debiased Machine Learning (DML) и Partial Linear Regression для задач Causal Inference!

Всем привет! Меня зовут Максим Кан, и я продуктовый аналитик в Островке. В этой статье мы начинаем серию постов, посвящённых Causal Inference и Causal ML как альтернативным методам проверки гипотез и оценки эффекта от treatment-а при отсутствии возможности провести A/B-тест.

Вымышленный кейс

Представим, что вы работаете аналитиком в неком онлайн-сервисе для бронирования отелей. После успешного оформления бронирования вашим клиентам предлагаются дополнительные ancillary-услуги (например, экскурсии, трансфер, аренда машины и т.д.). Вашего продакт-менеджера мучает вопрос: «Влияет ли покупка ancillary на дальнейшую активность пользователя на сервисе?» Возможно, опыт, связанный с вашими ancillary-услугами, настолько положительный, что пользователи начинают чаще использовать ваш сервис в дальнейшем. С другой стороны, покупка дополнительных услуг после бронирования жилья подразумевает дополнительные траты со стороны клиента, что может, в свою очередь, привести к меньшим тратам в ближайшем будущем. Давайте введём несколько обозначений:

•  — значение переменной воздействия (treatment) для i-го объекта. Переменная принимает значение 1, если объект подвергся воздействию, и 0 — в противном случае. В контексте нашего кейса: для пользователей, которые приобрели ancillary​ = 1, а для остальных = 0.

• ​ — фактическое значение переменной результата (outcome) для i-го объекта. В нашем случае это будет количество забронированных ночей в течение 6 месяцев после изначального бронирования. Требуется определить, существует ли причинно-следственная связь между и.

Главная причина, по которой трудно проверить гипотезу о наличии причинно-следственной связи между и с помощью A/B-теста, — это невозможность рандомизации. Воздействием в данном случае является добровольное действие пользователя, и мы не можем заставить случайную часть пользователей оформлять ancillary после бронирования жилья. Более того, даже если мы проведём A/B-тест, где будем предлагать ancillary случайной половине пользователей и не предлагать другой, — наблюдаемая оценка эффекта хоть и будет несмещённой, но будет являться оценкой эффекта от предложения ancillary, а не от покупки ancillary. Поэтому мы прибегнем к альтернативным методам Causal Inference!

Confounding Bias

Известная причина, по которой Causal Inference-анализ без A/B-тестирования не так тривиален, — это наличие различного рода смещений в данных. Самый часто встречающийся вид смещения — это Confounding Bias. Данный вид смещения возникает, когда существуют переменные (конфаундеры), одновременно влияющие как на treatment , так и на результат . Давайте обозначим все переменные-конфаундеры вектором . Например, в нашем случае одним из конфаундеров может являться тип поездки (командировка, отдых), ведь он может влиять на решение о покупке ancillary, т.к. на отдыхе люди больше склонны бронировать экскурсии, брать в аренду автомобиль, оформлять страховку и т.д. В то же время тип поездки также может влиять на количество бронирований в будущем, т.к. люди, которые ездят на отдых, как правило, делают это с большей регулярностью, чем те, кто ездит в командировки.

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

Популярный класс методик, который используется в Causal Inference для устранения Confounding Bias, — это Regression Adjustment. Идея данного подхода заключается в неком контроле (или фиксировании) всех переменных-конфаундеров . Мы уже обозначили, что одним из конфаундеров в нашем случае является тип поездки пользователя: если не учитывать эту переменную при оценке причинно-следственной связи, мы рискуем получить смещённую оценку эффекта. Например, можем ошибочно приписать высокую частотность последующих бронирований влиянию ancillary, тогда как на самом деле она объясняется тем, что отдыхающие пользователи просто по природе своей деятельности бронируют жильё чаще.

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

Обычная OLS и её проблемы

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

• Отсутствие ненаблюдаемых конфаундеров (unconfoundedness). Данное предположение говорит о том, что оценка эффекта, полученная через Regression Adjustment, будет несмещённой только в том случае, если в модель включены все переменные-конфаундеры.

• Stable Unit Treatment Value Assumption (SUTVA). Если вы хорошо знакомы с A/B-тестированием, то наверняка слышали и про SUTVA. Данное предположение говорит о том, что значение пользователя не должно зависеть от того, в какую группу попали другие пользователи. Например, факт того, что какой-то пользователь купил ancillary, не должен повлиять на поведение тех, кто ancillary не покупал.

• Overlap Assumption. Предположение о позитивности (или перекрытии) говорит о том, что для каждого значения должны существовать как , так и наблюдения. Другими словами, среди похожих наблюдений (одинаковых ) должны быть как treated-, так и control-объекты. Давайте вспомним идею контроля конфаундеров: она заключается в том, что, сгруппировав наблюдения по одинаковому уровню и делая сравнения только в рамках одной группы, мы устраняем смещение, так как сравниваем сопоставимых пользователей. Звучит логично: ведь если условный тип поездки является конфаундером, то разумно сравнивать тех, кто ездит отдыхать, с себе подобными, а тех, кто ездит в командировки, — с себе подобными. Но что, если среди пользователей, которые ездят в командировку, нет тех, кто покупал ancillary? Или, наоборот, если среди пользователей, которые ездили в отпуск, нет тех, кто ancillary не покупал? В таком случае в этих подгруппах просто не будет наблюдений для сравнения, а значит, эффект будет проблематично идентифицировать.

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

Представим, что мы определили список переменных-конфаундеров и собрали все необходимые данные в pd.DataFrame() 👇

import pandas as pd
import numpy as np

df.sample(5)

• ancillary_purchased — индикатор покупки ancillary при изначальном бронировании. Переменная воздействия D.

• nights_booked_6m — количество забронированных ночей на сервисе в течение 6 месяцев после изначального бронирования. Переменная результата Y.

• device_type — устройство с которого сделано бронирование (desktop / mobile / app). Является конфаундером, т.к. Клиенты через приложение более вовлечены в платформу → с большей вероятностью купят ancillary (D=1) → и с большей вероятностью вернутся снова (высокий Y).

• trip_type — цель поездки (leisure / business / mixed). Конфаундер, потому что leisure-пользователи системно отличаются по обоим направлениям — они чаще покупают ancillary и имеют более высокую частотность бронирований.

• has_loyalty — индикатор участия в программе лояльности платформы (0/1). Участники программы также более склонны покупать ancillary и также имеют более высокую частотность бронирований.

• bookings_prev_365 — количество бронирований за предыдущие 365 дней. Является конфаундером, потому что пользователи с большим количеством бронирований больше знакомы с сервисом -> и более высокую вероятность покупки ancillary.

• avg_order_value — средний чек за предыдущие 365 дней. Конфаундер аппроксимирующий доход и платежеспособность. Клиент с высоким средним чеком предположительно больше готов платить за дополнительные услуги → чаще берёт ancillary.

Данные выше являются сгенерированными, а истинный эффект от на заложенный на DGP равен +4.63 ночей.

Давайте дадим оценку эффекта используя длинную регрессию OLS 👇

import statsmodels.formula.api as smf

model = smf.ols('nights_booked_6m ~ ancillary_purchased + C(trip_type) + C(device_type) + has_loyalty + bookings_prev_365 + avg_order_value', data=df).fit()
print(model.summary())

OLS оценила эффект от покупки ancillary в +12.74 ночи, что значительно выше истинного эффекта, который равен +4.63 ночи. Более того, 95%-й доверительный интервал оценки также не накрывает истинное значение эффекта.

Почему же так произошло? Мы включили все переменные-конфаундеры в модель, но всё же получили смещённую оценку эффекта. Дело в том, что на этапе генерации данных была заложена нелинейная связь конфаундеров с переменными и . В свою очередь, функциональная форма OLS предполагает линейные и аддитивные эффекты регрессоров на — поэтому, если зависимость конфаундеров с и нелинейная, регрессия даст смещённую оценку эффекта. Один из способов решить данную проблему — Double/Debiased Machine Learning, предложенный Виктором Черножуковым в оригинальной статье.

Идея Double/Debiased Machine Learning

Прежде чем мы приступим к решению проблемы выше, — интересный факт: такую же оценку эффекта (+12.74), которую мы ранее получили длинной регрессией, можно получить иначе, применив разложение Фриша — Во — Ловелла:

1. Регрессия D на X и получение остатков  D.

2. Регрессия Y на X и получение остатков Y.

3. Регрессия остатков Y на остатки D.

# получение остатков D
debiase_model = smf.ols('ancillary_purchased ~ C(trip_type) + C(device_type) + has_loyalty + bookings_prev_365 + avg_order_value', data=df).fit()
df['d_resid'] = debiase_model.resid

# получение остатков Y
denoise_model = smf.ols('nights_booked_6m ~ C(trip_type) + C(device_type) + has_loyalty + bookings_prev_365 + avg_order_value', data=df).fit()
df['y_resid'] = denoise_model.resid

# регрессия остатков y на остатки d
final_model = smf.ols('y_resid ~ d_resid', data=df).fit()
beta_d = final_model.params['d_resid']
print(beta_d.round(4))

# 12.7388

Модели debiase_model и denoise_model называются nuisance-моделями и контролируют конфаундеры. На первом этапе мы строим регрессию, которая предсказывает значение воздействия с помощью конфаундеров, а затем берём остатки этой модели. Полученные остатки d_resid — это та часть переменной воздействия, которую мы не смогли предсказать с помощью конфаундеров, а значит, остатки d_resid по определению ортогональны (некоррелированы) конфаундерам. Можно сказать, что d_resid — это ортогонализованный вариант воздействия.

На втором этапе мы строим регрессию nights_booked_6m на trip_type, device_type, has_loyalty, bookings_prev_365 и avg_order_value, а затем получаем остатки этой модели. Логика такая же: полученные остатки y_resid ортогональны конфаундерам, так как всю изменчивость, обусловленную конфаундерами, мы уже объяснили, а остатки — это то, что невозможно предсказать с помощью .

Идея Double Machine Learning (DML) состоит в том, чтобы использовать для ортогонализации ML модели (например, решающие деревья) для учёта нелинейных взаимосвязей конфаундеров с и . Для этого достаточно использовать ML-алгоритмы в качестве nuisance-моделей, а в качестве final_model использовать OLS — такая комбинация имеет название Partial Linear Regression (PLR). Формально PLR записывается как два уравнения:

где и — произвольные нелинейные функции конфаундеров, а — эффект от воздействия , который нас как раз и интересует. Главная проблема обычной OLS-регрессии состоит в том, что OLS ищет линейную аппроксимацию и — то есть предполагает, что каждый конфаундер влияет на и одинаково при любом своём значении. В нашем кейсе это предположение нарушено: на этапе DGP мы задали более сложные зависимости между переменными. Линейная модель аппроксимирует эти зависимости прямой, оставляя в остатках систематическую компоненту — часть конфаундинга, которую она не смогла объяснить. Эта необъяснённая часть коррелирует с , что нарушает условие ортогональности остатков и приводит к смещению оценки . DML решает эту проблему, заменяя линейные модели на ML-алгоритмы в первом и втором шаге разложения Фриша — Во — Ловелла: в нашем случае LightGBM захватывает нелинейные зависимости, т.к. является непараметрической моделью, и при достаточном объёме данных вычищает остатки от влияния конфаундеров вне зависимости от формы их связи с и .

Известное свойство ML-моделей — они чаще склонны к переобучению (overfitting), что, в свою очередь, приводит к смещённой оценке эффекта. В связи с этим при работе с DML применяются out-of-fold-предсказания, которые можно получить с помощью класса cross_val_predict из библиотеки sklearn.

Выбор модели для DML

Давайте попробуем LightGBM в качестве nuisance-моделей. Работа с DML отличается от классического ML отсутствием этапов feature-инжиниринга и feature-selection. Согласно теории Causal Inference, вам необходимо проконтролировать в модели все переменные-конфаундеры для получения несмещённой оценки эффекта воздействия D, так что список фичей в обеих nuisance-моделях фиксируется на этапе построения причинно-следственного графа.

Тем не менее нам ничего не мешает заняться тюнингом модели. На практике чаще всего оказывается, что чем выше скорость сходимости обеих nuisance-функций и чем меньше ошибка предсказания nuisance-моделей, тем меньше смещение финальной оценки эффекта. Поэтому, подбирая гиперпараметры в DML, используют RMSE в качестве параметра оптимизации для вещественных переменных и LogLoss — для бинарных. Воспользуемся optuna с KFold для подбора 👇

import lightgbm as lgb
from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import KFold
from sklearn.metrics import root_mean_squared_error, log_loss
import optuna

kf = KFold(n_splits=5, shuffle=True, random_state=42)

# создадим dummy переменные для категориальных фичей
for varib in ['trip_type', 'device_type']:
    df_dummy = pd.get_dummies(df[[varib]], drop_first=True)
    df = pd.concat([df, df_dummy], axis=1)


# датафрейм с конфаундерами
X = df[['trip_type_leisure', 'trip_type_mixed', 'device_type_desktop', 'device_type_mobile', 'has_loyalty', 'bookings_prev_365', 'avg_order_value']].copy()

def objective_y(trial, X, y):
    params = {
        'n_estimators':      trial.suggest_categorical('n_estimators', [50, 100, 200, 300]),
        'max_depth':         trial.suggest_categorical('max_depth', [3, 4, 6]),
        'num_leaves':        trial.suggest_categorical('num_leaves', [15, 31, 63]),
        'learning_rate':     trial.suggest_categorical('learning_rate', [0.05, 0.1]),
        'min_child_samples': trial.suggest_categorical('min_child_samples', [20, 50]),
        'random_state': 42,
        'verbose': -1,
    }
    scores = []
    for train_idx, val_idx in kf.split(X):
        X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
        y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

        model = lgb.LGBMRegressor(**params)
        model.fit(X_train, y_train)
        preds = model.predict(X_val)
        scores.append(root_mean_squared_error(y_val, preds))

    return np.mean(scores)


def objective_d(trial, X, d):
    params = {
        'n_estimators':      trial.suggest_categorical('n_estimators', [50, 100, 200, 300]),
        'max_depth':         trial.suggest_categorical('max_depth', [3, 4, 6]),
        'num_leaves':        trial.suggest_categorical('num_leaves', [15, 31, 63]),
        'learning_rate':     trial.suggest_categorical('learning_rate', [0.05, 0.1]),
        'min_child_samples': trial.suggest_categorical('min_child_samples', [20, 50]),
        'random_state': 42,
        'verbose': -1,
    }
    scores = []
    for train_idx, val_idx in kf.split(X):
        X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
        d_train, d_val = d.iloc[train_idx], d.iloc[val_idx]

        model = lgb.LGBMClassifier(**params)
        model.fit(X_train, d_train)
        proba = model.predict_proba(X_val)[:, 1]
        scores.append(log_loss(d_val, proba))

    return np.mean(scores)


# тюнинг для модели предсказания D
study_debiase = optuna.create_study(direction='minimize')
study_debiase.optimize(
    lambda trial: objective_d(trial, X, df['ancillary_purchased']),
    n_trials=50,
    show_progress_bar=True,
)

# тюнинг для модели предсказания Y
study_denoise = optuna.create_study(direction='minimize')
study_denoise.optimize(
    lambda trial: objective_y(trial, X, df['nights_booked_6m']),
    n_trials=50,
    show_progress_bar=True,
)

Обратите внимание, что в случае бинарной переменной я использую классификатор, а не регрессор, и в качестве предсказания возвращаю вероятность положительного класса методом .predict_proba().

После тюнинга nuisance-моделей мы можем получить oof-предсказания переменных D и Y, а затем вычислить их остатки 👇

denoise_model = lgb.LGBMRegressor( **study_denoise.best_params, random_state=42, verbose=-1)
debiase_model = lgb.LGBMClassifier(**study_debiase.best_params, random_state=42, verbose=-1)

preds_d = cross_val_predict(debiase_model, X, df['ancillary_purchased'], method='predict_proba', cv=5)[:, 1]
df['d_resid_lgbm'] = df['ancillary_purchased'] - preds_d

preds_y = cross_val_predict(denoise_model, X, df['nights_booked_6m'], cv=5)
df['y_resid_lgbm'] = df['nights_booked_6m'] - preds_y

Наконец, мы можем регрессировать остатки y_resid_lgbm на остатки d_resid_lgbm и посмотреть на оценку эффекта воздействия D. Обратите внимание, что в качестве финальной модели для регрессии остатков мы используем именно OLS.

model = smf.ols('y_resid_lgbm ~ d_resid_lgbm', data=df).fit()
print(model.summary())

Полученная оценка эффекта очень близка к истинному значению +4.63! Использование LightGBM в качестве nuisance-моделей позволило учесть сложные связи конфаундеров с и , что, в свою очередь, позволило получить несмещённую оценку каузального эффекта.

Применить DoubleML и PLR можно также, воспользовавшись классом DoubleMLPLR из библиотеки doubleml. В таком случае вам не придётся самостоятельно получать oof-предсказания, вычислять остатки и фитить OLS 👇

from doubleml import DoubleMLPLR, DoubleMLData

data_dml = DoubleMLData(
    df,
    y_col='nights_booked_6m',
    d_cols='ancillary_purchased',
    x_cols=['trip_type_leisure', 'trip_type_mixed', 'device_type_desktop', 'device_type_mobile', 'has_loyalty', 'bookings_prev_365', 'avg_order_value'],
)

dml = DoubleMLPLR(
    data_dml,
    ml_l=denoise_model,
    ml_m=debiase_model,
    n_folds=5,
    score='partialling out',
)

Более того, тюнинг ML-моделей через optuna вы также можете реализовать, используя метод .tune_ml_models(). По умолчанию метод будет использовать KFold-валидацию, RMSE для вещественных переменных и LogLoss — для бинарных.

def lgbm_params(trial):
    return {
        'n_estimators':      trial.suggest_categorical('n_estimators', [50, 100, 200, 300]),
        'max_depth':         trial.suggest_categorical('max_depth', [3, 4, 6]),
        'num_leaves':        trial.suggest_categorical('num_leaves', [15, 31, 63]),
        'learning_rate':     trial.suggest_categorical('learning_rate', [0.05, 0.1]),
        'min_child_samples': trial.suggest_categorical('min_child_samples', [20, 50]),
    }
    
dml.tune_ml_models(
    ml_param_space={
        'ml_l': lgbm_params,  
        'ml_m': lgbm_params, 
    },
    cv=5,
    optuna_settings={
        'n_trials': 50,
        'show_progress_bar': True,
    }
)

dml.fit()
print(dml.summary)

DoubleMLPLR с параметром score='partialling out' проделывает практически идентичный пайплайн с нашим «ручным DoubleML» — точечная оценка эффекта также практически совпадает.

Проверим OLS и DML на А/А симуляциях

На одном конкретном DGP выше DoubleML с использованием LightGBM в качестве nuisance-моделей действительно позволил нам правильно оценить эффект тритмента, в отличие от обычной OLS. Давайте на симуляциях сравним методы с точки зрения вероятности ошибок первого рода. Для этого мы отойдем от ранее рассматриваемого кейса и снова обратимся к библиотеке doubleml, которая имеет свои методы для DGP разных видов.

Мы будем использовать метод .make_plr_CCDDHNR2018(), который возвращает pd.DataFrame() с указанным количеством наблюдений и множеством ковариат, две из которых являются конфаундерами, а также вещественные переменные воздействия и результата. В случае make_plr_CCDDHNR2018() связь конфаундеров с и нелинейна — процесс DGP описывается в документации как:

Ковариаты . А nuisance-функции задаются как:

При вызове можно регулировать «уровень нелинейности» через параметр a1 для переменной воздействия и параметр b0 для переменной результата .

Давайте проделаем 1500 симуляций, в которых будем выполнять DGP выше, при этом эффект от тритмента будем закладывать равным нулю. После этого мы будем применять обычную OLS, «ручной DoubleML», а также DoubleMLPLR из библиотеки doubleml для оценки эффекта. В итоге мы сравним распределения полученных оценок и распределение p-value от каждого метода.

# функции для оценки эффектов

# обычная ols регрессия
def get_ols_ate(df, formula, d_col):
    import statsmodels.formula.api as smf
    model = smf.ols(formula, data=df).fit()
    
    ate = model.params[d_col]
    pvalue = model.pvalues[d_col]
    ci = model.conf_int(alpha=0.05).loc[d_col]
    ci_lower = ci[0]
    ci_upper = ci[1]
    
    return ate, pvalue, ci_lower, ci_upper

# ручной DoubleML, обе nuisance модели регрессионные
def manual_doubleml_ate(df, debiase_model, denoise_model, X, y_col, d_col):
    from sklearn.model_selection import cross_val_predict
    import statsmodels.formula.api as smf
    df = df.copy()
    preds_d = cross_val_predict(debiase_model, X, df[d_col], cv=5)
    df['d_resid'] = df[d_col] - preds_d
    preds_y = cross_val_predict(denoise_model, X, df[y_col], cv=5)
    df['y_resid'] = df[y_col] - preds_y
    model = smf.ols('y_resid ~ d_resid', data=df).fit()
    ate = model.params['d_resid']
    pvalue = model.pvalues['d_resid']
    ci = model.conf_int(alpha=0.05).loc['d_resid']
    ci_lower = ci[0]
    ci_upper = ci[1]
    return ate, pvalue, ci_lower, ci_upper

# DoubleML из пакета doubleml
def package_doubleml_ate(df, debiase_model, denoise_model, X, y_col, d_col):
    from doubleml import DoubleMLPLR, DoubleMLData
    df = df.copy()
    data_dml = DoubleMLData(
        df,
        y_col=y_col,
        d_cols=d_col,
        x_cols=X.columns.to_list(),
    )
    dml = DoubleMLPLR(
        data_dml,
        ml_l=denoise_model,
        ml_m=debiase_model,
        n_folds=5,
        score='partialling out',
    )
    dml.fit()
    ate = dml.coef[0]
    pvalue = dml.pval[0]
    ci = dml.confint(level=0.95)
    ci_lower = ci.iloc[0, 0]
    ci_upper = ci.iloc[0, 1]
    return ate, pvalue, ci_lower, ci_upper

Проделаем 1500 итераций, предварительно затюнив debiase_model и denoise_model модели 👇

# фиксируем тюненые nuisance-модели (обе регрессионные)
denoise_model = lgb.LGBMRegressor(**study_denoise.best_params, random_state=42, verbose=-1)
debiase_model = lgb.LGBMRegressor(**study_debiase.best_params, random_state=42, verbose=-1)

from collections import defaultdict
from tqdm import tqdm

results_h0 = defaultdict(list)

# H0: истинный эффект равен нулю (alpha=0.0)
for sim in tqdm(range(1500)):
    data = make_plr_CCDDHNR2018(n_obs=50000, alpha=0.0, return_type='DataFrame', dim_x=3, a_1=5.0, b_0=5.0, random_state=sim)
    x_cols = [c for c in data.columns if c not in ['y', 'd']]
    X = data[x_cols]

    # прогоняем три метода на одном и том же датасете
    for method_name, func, kwargs in [
        ('ols', get_ols_ate, dict(df=data, formula='y ~ d + ' + ' + '.join(x_cols), d_col='d')),
        ('manual_dml', manual_doubleml_ate, dict(df=data, debiase_model=debiase_model, denoise_model=denoise_model, X=X, y_col='y', d_col='d')),
        ('package_dml', package_doubleml_ate, dict(df=data, debiase_model=debiase_model, denoise_model=denoise_model, X=X, y_col='y', d_col='d')),
    ]:
        ate, pval, ci_lower, ci_upper = func(**kwargs)
        results_h0[method_name].append({'ate': ate, 'pval': pval, 'ci_len': ci_upper-ci_lower})

# собираем результаты по каждому методу
df_ols = pd.DataFrame(results_h0['ols'])
df_manual_dml = pd.DataFrame(results_h0['manual_dml'])
df_package_dml = pd.DataFrame(results_h0['package_dml'])

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

Как и ожидалось, при наличии нелинейных связей обычная OLS-регрессия имеет неравномерное распределение p-value, а доля ложноположительных результатов > alpha. В то же время ручной DML и DoubleML из пакета doubleml показали равномерное распределение p-value, что говорит о том, что мы не завышаем вероятность ошибки первого рода.

Давайте также рассмотрим случай, когда в данных всё-таки имеется эффект от тритмента . Зададим параметр alpha = 0.5 при вызове метода .make_plr_CCDDHNR2018() и тем самым добавим в пилотную группу эффект, равный +0.5. В данном случае по каждому методу нас интересует bias — разница между истинным и наблюдаемым эффектом.

Среднее значение смещения (bias) находится около нуля в случае DML-оценщиков и достаточно сильно отклоняется от нуля в случае OLS-регрессии. Получается, что на DGP с нелинейными связями конфаундеров с и использование обычной OLS-регрессии, скорее всего, завысит долю ложноположительных результатов, а также сделает саму оценку эффекта систематически смещённой. В то же время применение Double Machine Learning значительно сокращает количество «ложных срабатываний» и уменьшает bias.

Заключение и рекомендации

Мы разобрали Double/Debiased Machine Learning как инструмент для оценки эффекта от добровольного воздействия — в тех ситуациях, когда рандомизация невозможна по природе самого treatment’а. Ключевая идея метода проста: заменить линейные nuisance-модели на ML-алгоритмы в разложении Фриша — Во — Ловелла и тем самым снять ограничение на функциональную форму связи конфаундеров с D и Y.

• Важно сказать, что DML — далеко не единственный инструмент для таких задач. Например, если у вас есть хороший инструмент — переменная, которая влияет на , но не влияет напрямую на — можно использовать инструментальные переменные (IV). Каждый метод имеет свои предпосылки и области применения, и DML не является универсальным ответом на все вопросы.

• Отдельно хочется сказать про выбор ML-алгоритмов для nuisance-моделей. На практике оказывается, что разные алгоритмы — LightGBM, CatBoost, Random Forest — дают похожие результаты, если данных достаточно и конфаундеры правильно определены.

• Гораздо важнее то, что происходит до написания кода: построение причинно-следственного графа (DAG) и правильное определение списка конфаундеров. DoubleML никак не поможет, если вы забыли включить сильный конфаундер или включили медиатор вместо конфаундера. В этой статье мы намеренно сконцентрировались на самом DoubleML и меньше затрагивали такие темы, как выбор переменных для контроля и базовая теория Causal Inference в целом. Если вам интересно быстро с ними ознакомиться, то рекомендую отличную статью на Хабре, которая разбирает Regression Adjustment именно с этой стороны.

• Поскольку выполнение unconfoundedness в observational данных проверить невозможно, важно оценивать, насколько устойчив полученный эффект к возможным ненаблюдаемым конфаундерам — то есть насколько сильным должен быть пропущенный конфаундер, чтобы свести ваш результат к нулю. Хорошая новость в том, что в doubleml анализ чувствительности зашит «из коробки».

На практике самое важное ограничение, которое стоит держать в голове: DML решает проблему нелинейного конфаундинга через наблюдаемые переменные, но по умолчанию никак не поможет в случае наличия сильных ненаблюдаемых конфаундеров — тех переменных, которых у вас в данных просто нет. Хорошая новость в том, что DML стакается с другими методами Causal Inference — например, с Two-Way Fixed Effects (TWFE), который контролирует ненаблюдаемые time-invariant-конфаундеры через панельные данные.

Ссылки на дополнительные материалы:

• Метод DML впервые описан в работе Chernozhukov et al. (2017)

• Тюнинг гиперпараметров для DML

• Causal Inference for brave and true