Теорема о разбиении регрессоров: делаем CUPED аб-тесты в один шаг / Хабр

Хай!

Пишу эту статью для тех, кто уже знаком с CUPED, но ищет больше понимания этого метода и взгляда на него с другой стороны. Здесь я не буду детально объяснять базовый алгоритм CUPED аб-тестирования: про это уже достаточно материала в сети. Основное внимание уделим рассмотрению метода через призму регрессий. Цель статьи - познакомить читателя с теоремой, безумно полезной для понимания работы регрессий, а главное - продемонстрировать, как с помощью этой теоремы проводить CUPED тесты не в три последовательных шага (как в базовом алгоритме), а с помощью одной регрессии.

Содержание:

Освежаем в голове CUPED
Теорема: частый случай и интуиция
Связь теоремы с CUPED
CUPED тесты в один шаг!
Приложение: общий случай теоремы

Освежаем в голове CUPED

Для начала быстро освежим, что делает CUPED. CUPED - это Controlled-experiment Using Pre-Experiment Data. То есть метод аб-тестирования с использованием данных об участниках теста, известных о них на пред-тестовом периоде. Чаще всего в качестве таких данных используются значения тестируемой метрики, измеренные на участниках теста до начала эксперимента. Использование этой информации позволяет уменьшить дисперсию оцениваемой статистики, что повышает чувствительность теста. Оригина��ьный алгоритм CUPED выглядит так:

Построить регрессию вида $metric_i = \hat\alpha + \hat\theta * preMetric_i + \epsilon_i,$ где — значение метрики i‑го пользователя, полученное в период теста, а — значение метрики этого же пользователя, измеренное на пред‑тестовом периоде;
Используя оцененное значение $\hat\theta$ , посчитать следующий остаток: $metric^{(cuped)}_i = metric_i - \hat\theta * preMetric_i$ ;
Применить обычный стат. тест для сравнения посчитанных остатков в тестовой и контрольной группах.
На выходе получаем несмещенную оценку разницы между группами, но с меньшей дисперсией.

Я сам некоторое время не мог разобраться в "физике" этого алгоритма. Я понимаю математическое доказательство, что дисперсия оценки падает (почитать можно тут). Но какой смысл несет величина $metric^{(cuped)}_i$ ? Почему к ней можно применить стат-тест и получить адекватный и интерпретируемый результат сравнения двух групп? не математически, а так... чисто по-человечески? Ответ на эти вопросы мне дала теорема о разделении регрессоров, или как в английской литературе Frisch–Waugh–Lovell theorem. Заодно, она показала, как анализировать CUPED-тесты проще, чем это описано в оригинальной статье.

Теорема: частный случай и интуиция

При построении регрессии мы объясняем часть дисперсии зависимой переменной с помощью дисперсий объясняющих регрессоров. Простой случай - регрессия с двумя переменными:

$y_i = \hat\alpha + \hat{\beta} x_i + \hat{\gamma} z_i +\epsilon_i , \quad cov(x,y) \neq 0$

где, - объясняющие переменные, с некоторой ненулевой корреляцией; $\hat\alpha, \hat\beta, \hat\gamma$ - оценки коэффициентов регрессии; $\epsilon_i$ - остаток. Строя подобную регрессию, мы говорим, что вариация объясняется с помощью двух слагаемых: независимой (от ) частью дисперсии и независимой (от ) частью дисперсии . Плюс ошибка. Оценка коэффициентов регрессии с помощью МНК четко разделит, какая часть вариации игрик связана только с икс, а какая - только с зет. Коэффициенты бета и гамма "возьмут на себя" только нужную часть дисперсии зависимой переменной, даже если регрессоры немного скореллированы (то есть имеют общую часть дисперсии). Можно убедиться в этом на практике, используя теорему о разбиении регрессоров (для простого случая регрессии с двумя переменными).

Утверждение: оценка коэффициента $\hat{\beta}$ , полученная во множественной регрессии с двумя переменными вида
$y_i = \alpha + \hat{\beta} x_i + \hat{\gamma} z_i +\epsilon_i , \quad cov(x,y) \neq 0$
полностью идентична оценке $\hat{\beta}'$ , полученной с помощью следующей последовательности регрессий с одной переменной:
построить регрессию на : $y_i = \hat{\alpha}' + \hat{\gamma}' z_i + \delta^{yz}_i$ , посчитать остатки регрессии $\delta^{yz}_i = y_i - \hat{y_i}$ ;
построить регрессию на : $x_i = \hat{\alpha}'' + \hat{\gamma}'' z_i + \delta^{xz}_i$ , посчитать остатки $\delta^{xz}_i = x_i - \hat{x_i}$ ;
построить регрессию остатков из первого шага, на остатках из второго шага, то есть $\delta^{yz}_i = \hat{\beta}' \delta^{xz}_i + \epsilon_i$

Оценка $\hat\beta'$ в точности равна $\hat{\beta}$ по величине и значимости!

Теперь разберем, что это все значит. В регрессии на в первом шаге мы связываем вариацию игрик с дисперсией зет. Вычисляя остатки от такой регрессии, мы получаем часть дисперсии игрик, которая не связана с зет. На втором шаге мы аналогично вычисляем часть вариации икс, независимую от зет, то есть остаток регрессии на . Последний шаг - самый интересный. Берем часть , несвязанную с , и берем часть , нескореллированную с . Регрессируем первое на второе, то есть $\delta^{yz}$ на $\delta^{xz}$ , и оп, получаем ровно тот же коэффициент, что был при икс в изначальной регрессии с двумя переменными.

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

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

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
np.random.seed(321)

# setting coefficients 
alpha = 1
beta = 1.2
gamma = 0.8

# number of observations and noise
n_obs = 100
std_noise = 5
noise = np.random.normal(0, std_noise, n_obs)

# generating correlated x and z
mean_x, mean_z = 0, 0
std_x, std_z = 1, 1
corr_xz = 0.1
cov_m = [[std_x**2        , std_x*std_z*corr_xz], 
        [std_x*std_z*corr_xz,          std_z**2]] 
x, z = np.random.multivariate_normal([mean_x, mean_z], cov_m, n_obs).T

# constructing y
y = alpha + beta * x + gamma * z + noise

# transforming to dataframe
data = pd.DataFrame({'x': x, 'z': z, 'y': y})

# regression with two variables 
reg_0 = smf.ols('y ~ 1 + x + z', data=data).fit()

# algorithm from theorem
reg_yz = smf.ols('y ~ 1 + z', data=data).fit()
data['delta_yz'] = data['y'] - reg_yz.predict()

reg_xz = smf.ols('x ~ 1 + z', data=data).fit()
data['delta_xz'] = data['x'] - reg_xz.predict()

reg_deltas = smf.ols('delta_yz ~ 1 + delta_xz', data=data).fit()

# comparing results 
print("beta =",     round(reg_0.params['x'], 5), 
      '; pvalue =', round(reg_0.pvalues['x'], 3))
print("beta' =",    round(reg_deltas.params['delta_xz'], 5), 
      "; pvalue' =", round(reg_deltas.pvalues['delta_xz'], 3))

# output:
# beta = 1.21274 ; pvalue = 0.014
# beta' = 1.21274 ; pvalue' = 0.014

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

Алгоритм для случая, когда и нескореллированы:
построить регрессию на : $y_i = \hat{\alpha}' + \hat{\gamma}' z_i + \delta^{yz}_i$ , посчитать остатки регрессии $\delta^{yz}_i = y_i - \hat{y_i}$ ;
построить регрессию этих остатков на : $\delta^{yz}_i = \hat{\beta}' x_i + \epsilon_i$ .
Всё, полученная оценка $\hat\beta'$ идентичная оценке $\hat\beta$ из регрессии: $y_i = \hat\alpha + \hat\beta x_i + \hat\gamma z_i + \epsilon_i$ .

И именно этот случай имеет непосредственное отношение к CUPED аб-тестам.

Связь теоремы с CUPED

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

$metric_i = \hat\alpha + \hat\beta * \mathbb I [user_i \in testSample] + \epsilon_i ,$

где metirc_i - метрика i-го пользователя (две группы объединены в один вектор), $\mathbb I$ - индикатор, принимающий значение 1, если элемент относится к тестовой выборке, и 0 - если к контрольной. Коэффициент $\beta$ будет равен разнице средних между группами, а его значимость будет в точности совпадать с p-value обычного t-теста.

Окей, значит алгоритм CUPED на самом деле полностью совпадает с алгоритмом из теоремы (для случая когда регрессоры нескореллированы). Значит к CUPED можно применить те же рассуждения об интуиции. В тесте мы измеряем разницу в выборочных средних, а дисперсия выборочного среднего прямо пропорциональна дисперсии самой метрики. Дисперсию метрики можно в свою очередь представить в виде двух слагаемых: некоторой постоянной части $\mathbb Var(Sample)$ , связанной с тем, что выборка состоит из разных людей с разным уровнем активности, и некоторой компоненты $\mathbb Var(PeriodNoise)$ , вызванной тем, что активность людей из выборки не постоянна, а меняется от периода к периоду. То есть:

$\mathbb Var(\hat\beta) \propto \mathbb Var(metric) = \mathbb Var (Sample) + \mathbb Var(PeriodNoise),$

где за $\hat\beta$ обозначена оцениваемая разница в средних между группами. При фиксированной выборке первое слагаемое постоянно от периода к периоду, то есть одинаково как на тестовом, так и на предтестовом периодах. Когда при CUPED тестировании мы строим регрессию вида $metric_i = \hat\alpha + \hat\theta * preMetric_i + \epsilon_i$ , компонента $\hat\theta * preMetric_i$ описывает именно часть дисперсии выборки, связанную с неоднородностью участников. Тогда получается, что смысл величины $metric_i^{(cuped)} = metric_i - \hat\theta * preMetric_i$ - это часть вариации метрики, очищенная от эффекта неоднородности выборки, и включающая в себя только вариацию связанную с изменениями во время эксперимента.

$\mathbb Var(\hat\beta') \propto \mathbb Var(metric^{(cuped)}) = \mathbb Var(metric) - \mathbb Var(Sample) = \mathbb Var(PeriodNoise)$

С интерпретацией разобрались. А значит наконец можно перейти практической части - альтернативному алгоритму CUPED тестов, который состоит из одного шага.

CUPED тесты в один шаг!

Соберем вместе три вывода, которые мы получили выше:

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

После прочтения этих фактов уже должен напрашиваться вопрос: а можем ли мы развернуть теорему в другую сторону и представить алгоритм CUPED в виде множественной регрессии? Конечно можем! Ради этого и затевалась вся статья!

Утверждение: оригинальный алгоритм CUPED полностью идентичен линейной регрессии следующего вида:
$metric_i = \hat\alpha + \hat\beta * \mathbb I [i \in testSample] + \hat\gamma * preMetric_i + \epsilon_i$ ,
где - значение метрики i-го пользователя во время эксперимента, - значение метрики этого же пользователя на предтестовом периоде, $\mathbb I [i \in testSample]$ - индикатор принадлежности пользователя тестовой группе (1 если пользователь принадлежит тестовой группе, 0 - если контрольной).

Оценка коэффициента $\hat\beta$ будет отражать величину и значимость разницы средних между группами, а ее значимость и дисперсия будет аналогична дисперсии оценки из стандартного алгоритма CUPED. Вот пруф:

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
from scipy.stats import ttest_ind
np.random.seed(321)

n_obs = 1000
uplift = 0.12   # +12%

# generate pre-experiment metrics of clients
pre_c = np.random.lognormal(mean=0, sigma=1, size=n_obs)
pre_t = np.random.lognormal(mean=0, sigma=1, size=n_obs)

# generate noise 
noise_c = np.random.lognormal(mean=0, sigma=0.8, size=n_obs)
noise_t = np.random.lognormal(mean=0, sigma=0.8, size=n_obs)

# construct experiment metrics, using pre-experiment, uplift and noise
c = pre_c * (1 + noise_c)
t = pre_t * (1 + uplift + noise_t)

# pack into dataframe
df = pd.DataFrame({'preMetric': np.append(pre_t, pre_c),
                   'metric': np.append(t, c), 
                   'testDummy': np.append([1] * len(t), [0] * len(c))})


# CUPED basic algorithm
cuped_reg = smf.ols('metric ~ 1 + preMetric', data=df).fit()
delta = (df.metric - cuped_reg.params['preMetric'] * df.preMetric).values
delta_t, delta_c = delta[0:len(t)], delta[len(t):]
t_stat_cuped, pvalue_cuped = ttest_ind(delta_t, delta_c, equal_var=False)

# CUPED via single regression using theorem
theorem_reg = smf.ols('metric ~ 1 + testDummy + preMetric', data=df).fit()
pvalue_theorem = theorem_reg.pvalues['testDummy']

# compare results
print('p-value cuped =', round(pvalue_cuped, 5))
print("p-value regression =", round(pvalue_theorem, 5))

# Output: 
# p-value cuped = 0.013
# p-value regression = 0.013

Вот и всё! Теперь вы умеете анализировать CUPED тесты в одну регрессию. Весь анализ умещается в строчке statsmodels.ols('metric ~ 1 + testDummy + preMetric', data=df).fit(). И никаких трехшаговых алгоритмов с вычислением остатков. Понятно, что это не прям безумно полезный лайфхак - базовый алгоритм CUPED умещается в 3-4 строчки кода. Основная ценность этого умения, да и всей этой статьи в целом, заключается в другом: узнать о теореме, демонстрирующей принципы работы регрессии, узнать о ее связи с CUPED аб-тестами, взглянуть на этот метод под другим углом и как следствие углубиться в его понимании. Надеюсь, было полезно!

Приложение: общий случай теоремы

Привести формулировку теоремы о разбиении регрессоров только для частного случая было бы непростительно. Закроем этот гештальт тут. Пусть есть множественная регрессия с двумя (условными) группами регрессоров:

$y_i = \alpha + \beta_1 x_{1,i} + ...+ \beta_n x_{n,i} +\gamma_1 z_{1,i} + ... + \gamma_m z_{m,i} + \epsilon_i$

Перепишем уравнение в векторном виде, соединив все и в две группы (константу $\alpha$ также для удобства можно отнести к одной и групп). Получим следующий вид уравнения:

$\mathbb Y = \mathbb X \beta + \mathbb Z \gamma + \epsilon,$

где $\beta$ и $\gamma$ - это векторы коэффициентов регрессии при матрицах регрессоров $\mathbb X$ и $\mathbb Z$ .

Утверждение: оценка вектора коэффициентов $\beta$ , полученная в данной регрессии, идентична по величине и значимости оценке коэффициента $\beta'$ в регрессии следующего вида:
$\mathbb M_Z \mathbb Y = \mathbb M_Z \mathbb X \beta' + \epsilon'$
где за $\mathbb M_Z \mathbb Y$ обозначено ортогональное дополнение к образу $\hat{ \mathbb Y} = \mathbb Z \hat{\gamma}.$ Проще говоря, остатки от регрессии $\mathbb Y$ на группе регрессоров $\mathbb Z$ . Аналогичная интерпретация у $\mathbb M_Z \mathbb X$ - остатки от регрессии $\mathbb Y$ на $\mathbb Z$ .

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

Полезные материалы и источники

Вывод теоремы о разбиении регрессоров (Frisch-Waugh-Lovell theorem).
Детальный разбор CUPED тестирования с упоминанием с теоремой о разбиении регрессоров

Теорема о разбиении регрессоров: делаем CUPED аб-тесты в один шаг