Привет, Хабр! Меня зовут Игорь Пантелеев, я Applied Data Scientist в компании Garage Eight. Сейчас моя команда занимается развитием одного из разделов сайта разрабатываемого нами продукта. 

В прошлом квартале мы задались вопросом: как оценить эффект от времени, которое пользователь проводит в нашем разделе, на Retention Rate (RR)? Казалось бы, решение очевидное: провести A/B-тест, но на поверку всё оказалось не так просто. В статье разберем, как у нас получилось определить эффект, с какими сложностями столкнулись в процессе и как нам помог метод Generalized Propensity Score.

Первым делом мы подумали о том, чтобы провести A/B-тест. Но нам он не подходил по нескольким причинам: 

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

2. Если бы мы отключили его для новой аудитории, это было бы долго, дорого, а главное — рискованно. 

3. Кроме того, A/B-тест не показал бы, как увеличение времени в разделе влияет на увеличение вероятности RR. Например, если сделаем так, чтобы пользователь обязательно проводил в разделе на 10 секунд дольше, как это повлияет на вероятность RR? 

Оценив эти риски, мы поняли, что не можем провести рандомизированное исследование на текущий момент. При этом мы можем воспользоваться историческими данными, имитируя исследование, и разделить уже посещавших раздел пользователей на тестовую и контрольную группу. На основе этого мы могли бы посчитать Outcome — метрику для оценки эффекта влияния времени в разделе на Retention Rate (далее — RR). Таким образом, вместо реализации эксперимента мы обратились к его имитации и причинно-следственной аналитике. 

Propensity Score Matching 

Для начала мы решили использовать метод Propensity Score Matching. Суть метода в том, чтобы на основе конфаундеров пользователя построить модель Propensity Score, смэтчить пользователей по этому показателю и сравнить значение Outcome между группами, дополнительно проверив баланс ковариат между ними.

Этот метод очень подробно объяснен в статье The Central Role of the Propensity Score in Observational Studies for Causal Effects. Rosenbaum, P. R., & Rubin, D. B. (1983), а также представлен в публикациях на Хабре. Поэтому не буду углубляться в объяснения теории, посоветую изучить для понимания публикации от Алексея Терентьева и Вячеслава Назарова.

Важно учесть, что все методы причинно-следственной аналитики требуют включать в модель максимально доступное число конфаундеров — свойств объектов, которые лучше всего объясняют причину события. К примеру, температура на улице может быть конфаундером при построении модели продажи мороженого. 

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

Скрытый текст

Отдельная боль про потерю информации

Когда непрерывное воздействие T ∈ ℝ (например, количество часов физической активности в неделю) переводится в бинарную переменную,

мы теряем существенную часть информации о дозировке воздействия.


В терминах теории информации энтропия бинарной переменной равна

где p = P(T > c), тогда как для непрерывной величины используется дифференциальная энтропия

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

На уровне оценок это выражается в том, что вместо функции отклика

мы вынуждены анализировать

то есть усредненный эффект «выше/ниже порога». 

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

Сопоставив все преимущества и недостатки, мы поняли, что PSM нам не подходит, и обратились к GPS.

Generalized Propensity Score

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

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

Разберемся на простом примере, почему важно не терять другие факторы. Давайте представим, что у нас есть популяция людей, в которой поровну жаворонков и сов. Мы знаем, что жаворонки встают в 5–6 утра и ложатся в 8–9 вечера. Для сов типичен иной образ жизни, они встают в 12–14 и ложатся в 3–4 часа ночи. Кроме данных о режиме дня, у этой популяции есть дополнительный набор свойств, которые нам известны: пол, возраст, работа и другие. 

Если мы знаем эти характеристики и тип человека (сова или жаворонок), то можем построить модель, которая будет определять, в привычное ли для себя время встал человек или ��то нетипичное поведение. Если время стандартное, мы можем измерить дополнительную метрику, например, покупал ли человек в этот день кофе. Таким образом, сова, которая встала в 6 утра, не будет вносить никакого импакта в конечную метрику, потому что для нее такое поведение нетипично. А вот жаворонок, проснувшийся в 6 и купивший кофе, добавит к целевой метрике некое число.

Поэтому для учета и других факторов мы обратились к методу Generalized Propensity Score, который описан в статье The Propensity Score with Continuous Treatments Keisuke Hirano and Guido W. Imbens (2004). 

Чтобы применить метод GPS, нам нужно было пройти несколько этапов: 

1. Построить модель регрессии, которая бы для каждого человека по его характеристикам определяла матожидание времени, проведенного в разделе. 

2. Посчитать ошибку и дисперсию для каждого пользователя. 

3. Построить вторую модель, которая будет предсказывать дисперсию, характерную для человека с конкретными фичами. 

Посмотрим на конкретном примере. У нас есть клиент из Санкт-Петербурга, мужчина в возрасте 35 лет. С помощью первой модели мы определяем, что матожидание времени в разделе составляет 20 секунд с отклонением в 5 секунд. После этого строим нормальное распределение времени в разделе для него. Третьим шагом мы подставляем реальное значение времени в разделе в нормальное распределение, проверяем, насколько они различаются, и считаем RR спустя 4 недели (RR4W). С этими данными мы можем рассчитать GPS. 

Допустим, мы получим RR4W = 1 и GPS = 0,004. Это значит, что для этого пользователя влияние времени в разделе на RR минимально. Соответственно, далее мы сможем сравнить эффект выборочно для тех пользователей, у которых время в разделе минимально отклоняется от нормального распределения. Так мы сможем имитировать рандомизированное исследование, балансируя людей по их характеристикам. 

Таким образом, мы решили использовать GPS. Для этого взяли набор данных: 

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

• Характеристики человека — время на сайте в целом, время в других разделах, которые ведут на наш, количество сессий, возраст. 

Чтобы не создав��ть конкурирующие события в характеристиках, связанных со временем, взяли данные за предыдущую запись. 

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

Шаг 1. Модель математического ожидания Treatment по ковариатам

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

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

Шаг 2. Модель предсказания дисперсии Treatment по ковариатам

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

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

Шаг 3. GPS

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

log_gps_obs = norm.logpdf(T, loc=mu_hat, scale=sigma_hat)

gps_obs = np.exp(log_gps_obs)

df['gps'] = gps_obs

Шаг 4. Полиномиальная модель отклика

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

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

Чтобы построить модель, возьмем линейную регрессию из statsmodels, в качестве фич используем время в разделе, GPS как балансирующую фичу и добавим константу. Она нужна, чтобы модель могла правильно учитывать общий уровень (средний сдвиг) целевой переменной и не была привязана к нулю. В качестве таргета используем влияние на RR. Кроме того, так как данных немного, добавляем более консервативный метод оценки стандартных ошибок heteroskedasticity-consistent (HC3).

Шаг 5. Dose-response 

Мы получили обученную модель, которая объясняет влияние времени в разделе на RR. Посмотрим на примере: возьмем человека со матожиданием времени в разделе 30 секунд и отклонением в 10 секунд. Начинаем подставлять для него разные значения:

 Расчет приблизительный и будет зависеть от модели

Как можно заметить, с приближением к наиболее типичному времени GPS возрастает, а значит, оценка эффекта на RR для такого пользователя будет выше.

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

t_grid = np.linspace(df['T'].min(), df['T'].max(), 60)

mu_hat_t = []

n = df.shape[0]

for t_star in t_grid:

    log_gps_tstar = norm.logpdf(np.repeat(t_star, n), loc=mu_hat, scale=sigma_hat)

    gps_tstar = np.exp(log_gps_tstar)

    TP_star = poly.transform(np.column_stack([np.repeat(t_star, n), gps_tstar]))

    X_star = sm.add_constant(TP_star, has_constant='add')

    y_pred_i = X_star.dot(params)

    mu_hat_t.append(np.mean(y_pred_i))

mu_hat_t = np.array(mu_hat_t)

Шаг 6. Bootstrap

Бизнес есть бизнес, и никто нам не поверит, если мы просто скажем, что эффект есть. Поэтому поверх dose-response-кривой мы будем делать бутстрап, чтобы построить 95% доверительный интервал.

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

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

Шаг 7. График dose-response с доверительным интервалом

У нас получился такой график кривой dose-response. Из того, что мы видим, можно предположить, что первое время, проведенное в разделе, вызывает негатив, так называемый эффект новизны. Но как только пользователь проходит отметку в 100–110 секунд, привыкает и понимает преимущество, эффект начинает становиться положительным. 

Шаг 8. Балансировка

Ключевой этап во всех PS-производных методах — убедиться, что мы сравниваем между собой одинаковых по свойствам пользователей. Желательно проводить балансировку до и после добавления GPS, но в нашем примере этап «до» убран. 

Самое главное — проследить за тем, чтобы после добавления GPS у нас был баланс по ковариатам. Для этого мы вводим функцию нормализованной разницы средних. Если SMD < 0,1, значит, у нас достигнут баланс ковариат между группами. SMD > 0,25 указывает на явные проблемы с балансом. 

n_bins = 5

df['t_bin'] = pd.qcut(df['T'], q=n_bins, labels=False, duplicates="drop")

def smd(a, b):

    ma, mb = a.mean(), b.mean()

    sa, sb = a.std(ddof=1), b.std(ddof=1)

    pooled = np.sqrt((sa**2 + sb**2)/2.0)

    if pooled < 1e-8: # защита от деления на ноль

        return np.nan

    return (ma - mb) / pooled

smd_rows = []

problematic = []

for col in covariate_cols:

    cov = df[col].values

    GM = LinearRegression().fit(df[['gps']], cov.reshape(-1,1))

    cov_pred = GM.predict(df[['gps']]).ravel()

    cov_resid = cov - cov_pred

    

    raw_list, resid_list = [], []

    for b in range(n_bins):

        mask_bin = df['t_bin']==b

        mask_rest = ~mask_bin

        smd_raw = smd(cov[mask_bin], cov[mask_rest])

        smd_res = smd(cov_resid[mask_bin], cov_resid[mask_rest])

        raw_list.append(smd_raw)

        resid_list.append(smd_res)

    

    if np.any(np.isnan(raw_list)) or np.any(np.isnan(resid_list)):

        problematic.append(col)

    

    smd_rows.append({

        "covariate": col,

        "smd_raw_mean": np.nanmean(raw_list),

        "smd_resid_mean": np.nanmean(resid_list),

        "smd_raw_max": np.nanmax(raw_list),

        "smd_resid_max": np.nanmax(resid_list)

    })

smd_df = pd.DataFrame(smd_rows).sort_values('smd_raw_max', ascending=False)

print("\nBalance diagnostics (mean and max |SMD| across T bins):")

print(smd_df)

plt.figure(figsize=(7,5))

plt.scatter(smd_df['smd_raw_mean'], smd_df['smd_resid_mean'])

for , r in smddf.iterrows():

    plt.text(r['smd_raw_mean'], r['smd_resid_mean'], r['covariate'])

plt.plot([0, max(smd_df['smd_raw_mean'].max(), smd_df['smd_resid_mean'].max())],

         [0, max(smd_df['smd_raw_mean'].max(), smd_df['smd_resid_mean'].max())], linestyle='--')

plt.xlabel("Mean |SMD| before residualizing on GPS")

plt.ylabel("Mean |SMD| after residualizing on GPS")

plt.title("Does GPS reduce X–T association? (lower is better)")

plt.show()

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

После применения GPS по всем ковариатам удалось достичь баланса, а значит, сравнение пользователей внутри бинов выполнено корректно. 

Шаг 9. Второй этап оценки модели — проверках на случайных значениях

На этом этапе мы проверяем, не «подстраивается» ли на�� подход под любые данные. Для этого заменяем реальные значения Treatment (время, проведённое пользователем в SPOC) на случайные, взятые из того же распределения. Ожидаемый результат — нулевой эффект с небольшими отклонениями: при случайном воздействии влияние на модельный Outcome не должно проявляться систематически.

Это именно тот график, который мы и увидели:

• значение эффекта изменяется случайно по мере роста Treatment;

• широкий доверительный интервал;

• отклонения кривой лежат в пределах от 0.41 до 0.47.

Мы также проверяем веса Outcome модели:

Коэффициенты при T стали очень маленькими, а их p-value — больше 0.05 (0.939 и 0.144). Это говорит о том, что Treatment в таком варианте никак не объясняет Outcome, а сами коэффициенты не являются статистически значимыми.

Шаг 10. Оценка эффекта

Если нам нужно оценить эффект последних изменений, достаточно просто подставить их на нашу кривую. В нашем случае за два месяца время, проведенное в разделе, увеличилось с 90 до 120 секунд.

Интерполяция дала оценку в 0,4% прироста RR4W за два месяца, что для нашего раздела является вполне реалистичным результатом.

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

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

Посмотреть детальнее на весь процесс применения метода и покопаться в коде можно на Jupyter.

Это была моя первая статья, спасибо за внимание. Буду рад любой критике, объективной и не очень — делитесь в комментариях! =)