Хай!

Пишу эту статью для тех, кто уже знаком с 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 аб-тесты в один шаг