CUPED на практике: гайд на R и Python

Пошаговый разбор на синтетических данных: проверка на истории и MDE, новые пользователи, винсоризация, CUPED руками и через регрессию, ARPPU, утечка ковариаты, проверка на A/A
statistics
ab-testing
cuped
R
python
Author

Elena U

Published

October 12, 2026

Это практическое дополнение к посту в телеграм-канале про то, почему CUPED – это старая добрая ANCOVA. Там теория, а здесь пошаговый разбор с кодом на R и Python: как применить CUPED к A/B-тесту и не наступить на типичные грабли.

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

NoteПлан
  1. Данные (и как собрать такую таблицу у себя)
  2. Готовим ковариату: новые пользователи и выбросы
  3. До запуска теста: проверяем на истории, есть ли смысл в CUPED, и считаем MDE
  4. Тест закончился: CUPED руками
  5. То же самое через регрессию (и Multi-CUPED)
  6. ARPPU: то же самое на платящих
  7. Как сломать CUPED: ковариата из периода теста
  8. Проверяем на симуляции: A/A и мощность

В конце приложение с кодом симуляции из поста в канале.

1. Данные

Онлайн-библиотека, тест нового экрана оплаты, по 2 недели на тест. Данные лежат в двух файлах:

  • cuped_experiment.csv – сам эксперимент, 40 000 пользователей;
  • cuped_history.csv – те же метрики за прошлый период, до запуска теста. Понадобится в шаге 3.

Колонки:

  • group – контроль или тест;
  • is_new – пользователь пришёл недавно, истории до теста у него нет (таких 20%);
  • rev_pre – выручка за 2 недели до теста, у новых пользователей пропуск;
  • revenue – выручка за 2 недели теста, это наша метрика;
  • rev_week1 – выручка за первую неделю теста. Понадобится в шаге 7, чтобы показать, как делать не надо.

Истинный эффект в тестовой группе: +3% к выручке.

library(tidyverse)
library(estimatr)   # lm_robust(), lm_lin() с робастными ошибками и tidy()
df <- read_csv("data/cuped_experiment.csv")
head(df)
# A tibble: 6 × 6
  user_id group   is_new rev_pre revenue rev_week1
    <dbl> <chr>   <lgl>    <dbl>   <dbl>     <dbl>
1       1 control FALSE    149.     181.      75.0
2       2 control FALSE      0        0        0  
3       3 control FALSE     96.6      0        0  
4       4 control TRUE      NA      199.     107. 
5       5 test    FALSE    146.       0        0  
6       6 test    FALSE      0      168.      69.7
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.formula.api as smf

df = pd.read_csv("data/cuped_experiment.csv")
df.head()
   user_id    group  is_new     rev_pre     revenue   rev_week1
0        1  control   False  149.011673  180.641194   75.046791
1        2  control   False    0.000000    0.000000    0.000000
2        3  control   False   96.624863    0.000000    0.000000
3        4  control    True         NaN  199.284841  106.593654
4        5     test   False  146.227997    0.000000    0.000000

Связь между «до» и «во время» задаёт скрытая склонность пользователя платить (activity): кто платил раньше, скорее заплатит и сейчас. В реальных данных её не видно, у нас есть только выручка. Генератор понадобится ещё раз в шаге 8, для симуляции.

generate_data <- function(n = 40000, effect = 0.03, seed = 42) {
  set.seed(seed)
  tibble(
    user_id  = 1:n,
    group    = sample(c("control", "test"), n, replace = TRUE),
    is_new   = rbinom(n, 1, 0.2) == 1,   # нет истории до теста
    activity = rnorm(n)                  # скрытая склонность платить
  ) |>
    mutate(
      # предпериод: платил ли и сколько
      pay_pre   = rbinom(n, 1, plogis(-0.5 + 3 * activity)),
      rev_pre   = if_else(is_new, NA_real_,
                          pay_pre * exp(5 + 0.8 * activity + rnorm(n, 0, 0.5))),
      # период теста: та же склонность + эффект в тестовой группе
      pay       = rbinom(n, 1, plogis(-0.5 + 3 * activity)),
      revenue   = pay * exp(5 + 0.8 * activity + rnorm(n, 0, 0.5)) *
                  if_else(group == "test", 1 + effect, 1),
      # часть выручки теста, пришедшаяся на первую неделю
      rev_week1 = revenue * runif(n, 0.3, 0.7)
    ) |>
    select(user_id, group, is_new, rev_pre, revenue, rev_week1)
}

# так получены файлы:
# write_csv(generate_data(), "data/cuped_experiment.csv")
# write_csv(generate_data(effect = 0, seed = 2025) |> select(-group, -rev_week1),
#           "data/cuped_history.csv")
def generate_data(n=40_000, effect=0.03, seed=42):
    rng = np.random.default_rng(seed)
    activity = rng.normal(size=n)                     # скрытая склонность платить
    p_pay = 1 / (1 + np.exp(-(-0.5 + 3 * activity)))
    group = rng.choice(["control", "test"], size=n)
    is_new = rng.random(n) < 0.2                      # нет истории до теста

    # предпериод: платил ли и сколько
    pay_pre = rng.random(n) < p_pay
    rev_pre = pay_pre * np.exp(5 + 0.8 * activity + rng.normal(0, 0.5, n))
    rev_pre = np.where(is_new, np.nan, rev_pre)

    # период теста: та же склонность + эффект в тестовой группе
    pay = rng.random(n) < p_pay
    revenue = pay * np.exp(5 + 0.8 * activity + rng.normal(0, 0.5, n))
    revenue = revenue * np.where(group == "test", 1 + effect, 1)
    rev_week1 = revenue * rng.uniform(0.3, 0.7, n)    # часть выручки за первую неделю

    return pd.DataFrame({"user_id": np.arange(1, n + 1), "group": group, "is_new": is_new,
                         "rev_pre": rev_pre, "revenue": revenue, "rev_week1": rev_week1})

Нужна одна строка на пользователя эксперимента. Главное правило: предпериод отсчитывается от момента, когда этот пользователь попал в тест, а не от даты старта теста. Набросок на SQL:

with exposure as (            -- первое попадание пользователя в тест
  select user_id, group_name, min(event_time) as exposed_at
  from ab_exposures
  where test_id = 'paywall_v2'
  group by user_id, group_name
)
select
  e.user_id,
  e.group_name as "group",
  -- метрика за период теста: от попадания в тест до конца теста
  coalesce(sum(case when p.paid_at >= e.exposed_at
                     and p.paid_at < '2026-10-01' then p.amount end), 0) as revenue,
  -- ковариата: такое же окно (14 дней) ДО попадания в тест
  coalesce(sum(case when p.paid_at >= e.exposed_at - interval '14 days'
                     and p.paid_at < e.exposed_at then p.amount end), 0) as rev_pre,
  -- нет истории: зарегистрировался позже начала окна предпериода
  u.registered_at > e.exposed_at - interval '14 days' as is_new
from exposure e
join users u using (user_id)
left join payments p using (user_id)
group by e.user_id, e.group_name, u.registered_at, e.exposed_at

Для новых пользователей rev_pre здесь получится 0, а не пропуск. Это нормально, если есть флаг is_new, дальше мы делаем ровно так же.

2. Готовим ковариату

Сначала договоримся о подготовке данных: ровно её же мы применим и к истории в шаге 3, и к самому тесту.

Новые пользователи. У них rev_pre отсутствует. Выкидывать их нельзя: это часть эксперимента, и без них оценка эффекта будет про другую аудиторию. Стандартное решение, оно же описано в оригинальной статье про CUPED: заполнить пропуск константой (например, нулём) и добавить индикатор «нет истории» как отдельную ковариату. Константа может быть любой, индикатор всё равно отделит этих пользователей.

Выбросы. Пара пользователей с огромной выручкой раскачивает и среднее, и саму оценку θ. Поэтому до расчёта CUPED метрику и ковариату винсоризируют: обрезают по высокому перцентилю. Порог считаем по объединённым данным контроля и теста, а не отдельно по группам, иначе группы будут обработаны по-разному.

Warning

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

winsorize <- function(x, p = 0.995) pmin(x, quantile(x, p, na.rm = TRUE))

prepare <- function(df) {
  df |>
    mutate(
      revenue_w = winsorize(revenue),
      x_pre     = coalesce(winsorize(rev_pre), 0)   # пропуск у новых → 0, флаг is_new уже есть
    )
}

df <- prepare(df)
def winsorize(x, p=0.995):
    return np.minimum(x, np.nanquantile(x, p))

def prepare(df):
    df = df.copy()
    df["revenue_w"] = winsorize(df["revenue"])
    df["x_pre"] = winsorize(df["rev_pre"]).fillna(0)   # пропуск у новых → 0, флаг is_new уже есть
    return df

df = prepare(df)

3. До запуска теста: есть ли смысл в CUPED

CUPED умножает дисперсию на \(1 - \rho^2\), то есть снижает её на \(\rho^2\), где \(\rho\) – корреляция метрики с ковариатой. Решать, стоит ли игра свеч, нужно до запуска теста, а не по его результатам. Иначе это то же подглядывание: попробовали с CUPED и без, выбрали, где p-value меньше.

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

  • какую долю дисперсии метрики объясняет ковариата. Это и есть ожидаемое снижение дисперсии. Для одной ковариаты оно равно \(\rho^2\), для нескольких (Multi-CUPED) это \(R^2\) регрессии;
  • какой MDE получится с CUPED и без него при запланированном размере выборки.
history <- read_csv("data/cuped_history.csv") |> prepare()

rho_old  <- with(filter(history, !is_new), cor(x_pre, revenue_w))
r2_cuped <- cor(history$x_pre, history$revenue_w)^2
r2_multi <- summary(lm(revenue_w ~ x_pre + is_new, data = history))$r.squared

# MDE для двустороннего теста: alpha = 5%, мощность 80%
mde <- function(sd, n_per_group, alpha = 0.05, power = 0.8) {
  (qnorm(1 - alpha / 2) + qnorm(power)) * sd * sqrt(2 / n_per_group)
}
n_per_group <- 20000
sd_y <- sd(history$revenue_w)
mu_y <- mean(history$revenue_w)

tibble(
  method             = c("t-test", "CUPED: x_pre", "Multi-CUPED: x_pre + is_new"),
  variance_reduction = c(0, r2_cuped, r2_multi),
  mde_relative       = mde(sd_y * sqrt(1 - variance_reduction), n_per_group) / mu_y
) |>
  mutate(variance_reduction = scales::percent(variance_reduction, accuracy = 1),
         mde_relative       = scales::percent(mde_relative, accuracy = 0.1)) |>
  knitr::kable(col.names = c("Метод", "Снижение дисперсии", "MDE"))
Метод Снижение дисперсии MDE
t-test 0% 4.7%
CUPED: x_pre 33% 3.9%
Multi-CUPED: x_pre + is_new 36% 3.8% 
from scipy.stats import norm

history = prepare(pd.read_csv("data/cuped_history.csv"))

old = history[~history["is_new"]]
rho_old = old["x_pre"].corr(old["revenue_w"])
r2_cuped = history["x_pre"].corr(history["revenue_w"]) ** 2
r2_multi = smf.ols("revenue_w ~ x_pre + is_new", data=history.astype({"is_new": int})).fit().rsquared

# MDE для двустороннего теста: alpha = 5%, мощность 80%
def mde(sd, n_per_group, alpha=0.05, power=0.8):
    return (norm.ppf(1 - alpha / 2) + norm.ppf(power)) * sd * np.sqrt(2 / n_per_group)

n_per_group = 20_000
sd_y, mu_y = history["revenue_w"].std(), history["revenue_w"].mean()
pd.DataFrame({
    "Метод": ["t-test", "CUPED: x_pre", "Multi-CUPED: x_pre + is_new"],
    "Снижение дисперсии": [0, r2_cuped, r2_multi],
}).assign(MDE=lambda d: mde(sd_y * np.sqrt(1 - d["Снижение дисперсии"]), n_per_group) / mu_y).round(3)
                         Метод  Снижение дисперсии    MDE
0                       t-test               0.000  0.047
1                 CUPED: x_pre               0.332  0.039
2  Multi-CUPED: x_pre + is_new               0.356  0.038

Что видно:

  • У пользователей с историей корреляция \(\rho \approx\) 0.67, это обещало бы снижение дисперсии на 45%. Но по всей аудитории выходит меньше, 33%: у новых пользователей предпериода нет, и для них CUPED ничего не даёт. Поэтому ожидаемый выигрыш нужно считать на всех пользователях, которые попадут в тест, а не только на «удобных».
  • Индикатор is_new добавляет немного: новые пользователи платят иначе, чем старые с нулевой выручкой в предпериоде.
  • MDE с CUPED заметно меньше. На практике это значит, что тот же тест можно провести быстрее или на меньшей выборке.

Если бы на истории получилось \(\rho \approx 0.2\)–\(0.3\), снижение дисперсии было бы 4–9%, и возиться с CUPED, скорее всего, не стоило бы.

4. Тест закончился: CUPED руками

Формула одна: \(Y_{cuped} = Y - \theta (X - \bar{X})\), где \(\theta = \text{Cov}(X, Y) / \text{Var}(X)\). Два важных момента:

  • θ считаем один на обе группы, по объединённым данным. Если считать отдельно в контроле и тесте, оценка эффекта становится смещённой.
  • Центрирование \(X - \bar{X}\) нужно только для того, чтобы среднее метрики не поменялось. На оценку разницы между группами оно не влияет.

Для сравнения посчитаем и разность «после минус до» (\(\theta = 1\)), про которую была ловушка в посте.

theta <- cov(df$x_pre, df$revenue_w) / var(df$x_pre)

df <- df |>
  mutate(revenue_cuped = revenue_w - theta * (x_pre - mean(x_pre)),
         revenue_diff  = revenue_w - x_pre)

# t.test считает «контроль − тест», переворачиваем в «тест − контроль»
summarise_t <- function(tt, method) {
  tibble(method  = method,
         effect  = unname(diff(tt$estimate)),
         ci_low  = -tt$conf.int[2],
         ci_high = -tt$conf.int[1],
         p_value = tt$p.value)
}

t_plain <- t.test(revenue_w ~ group, data = df)
t_diff  <- t.test(revenue_diff ~ group, data = df)
t_cuped <- t.test(revenue_cuped ~ group, data = df)

bind_rows(summarise_t(t_plain, "t-test"),
          summarise_t(t_diff,  "Разность «после минус до»"),
          summarise_t(t_cuped, "CUPED")) |>
  mutate(variance_vs_ttest = ((ci_high - ci_low) / (ci_high[1] - ci_low[1]))^2) |>
  knitr::kable(digits = 3)
method effect ci_low ci_high p_value variance_vs_ttest
t-test 3.338 -2.039 8.715 0.224 1.000
Разность «после минус до» 3.971 -0.690 8.631 0.095 0.751
CUPED 3.749 -0.581 8.078 0.090 0.648 
is_test = df["group"] == "test"
theta = np.cov(df["x_pre"], df["revenue_w"])[0, 1] / df["x_pre"].var()
df["revenue_cuped"] = df["revenue_w"] - theta * (df["x_pre"] - df["x_pre"].mean())
df["revenue_diff"] = df["revenue_w"] - df["x_pre"]

def welch(col):
    t, c = df.loc[is_test, col], df.loc[~is_test, col]
    res = stats.ttest_ind(t, c, equal_var=False)
    ci = res.confidence_interval()
    return {"effect": t.mean() - c.mean(), "ci_low": ci.low, "ci_high": ci.high, "p_value": res.pvalue}

out = pd.DataFrame({"t-test": welch("revenue_w"),
                    "Разность «после минус до»": welch("revenue_diff"),
                    "CUPED": welch("revenue_cuped")}).T
out["variance_vs_ttest"] = ((out.ci_high - out.ci_low) / (out.ci_high.iloc[0] - out.ci_low.iloc[0])) ** 2
out.round(3)
                           effect  ci_low  ci_high  p_value  variance_vs_ttest
t-test                      3.338  -2.039    8.715    0.224              1.000
Разность «после минус до»   3.971  -0.690    8.631    0.095              0.751
CUPED                       3.749  -0.581    8.078    0.090              0.648

Эффект +3% в данных заложен, но обычный t-тест его не видит (p = 0.22): выручка слишком шумная. Оценка CUPED почти та же, а доверительный интервал заметно уже: дисперсия уменьшилась на 35%, p-value снизился с 0.22 до 0.09. Это близко к тому, что мы ожидали по истории в шаге 3 (33%), так что прогноз на истории работает.

Разность «после минус до» тоже немного помогла (дисперсия 75% от t-теста), потому что корреляция здесь больше 0.5. Но CUPED всё равно точнее. При слабой корреляции разность оказалась бы хуже обычного t-теста.

Значимым эффект не стал и с CUPED. Окажется ли он значимым на одном конкретном датасете – дело случая: CUPED не гарантирует «прокрас», а повышает его вероятность. На симуляции в шаге 8 мощность растёт примерно с 39% до 55%, то есть такой эффект CUPED будет находить заметно чаще.

TipОтносительный эффект

Продакты обычно спрашивают «на сколько процентов». Делить разницу нужно на среднее исходной (не скорректированной) метрики в контроле: здесь это 2.3%. Для строгого доверительного интервала относительного эффекта нужен дельта-метод или бутстреп.

5. То же самое через регрессию

CUPED – это регрессия revenue ~ group + x_pre, посчитанная в два шага. В одну строку это удобнее, а заодно сразу получаются робастные стандартные ошибки и можно добавить сколько угодно ковариат.

# CUPED = ANCOVA
m_cuped <- lm_robust(revenue_w ~ group + x_pre, data = df)

# Multi-CUPED: + индикатор новых пользователей
m_multi <- lm_robust(revenue_w ~ group + x_pre + is_new, data = df)

# поправка Лина: взаимодействие с центрированными ковариатами
m_lin <- lm_lin(revenue_w ~ group, covariates = ~ x_pre + is_new, data = df)

bind_rows(
  tidy(m_cuped) |> mutate(model = "CUPED: ~ group + x_pre"),
  tidy(m_multi) |> mutate(model = "Multi-CUPED: + is_new"),
  tidy(m_lin)   |> mutate(model = "Lin: group * (x_pre + is_new)")
) |>
  filter(term == "grouptest") |>
  select(model, estimate, std.error, conf.low, conf.high, p.value) |>
  knitr::kable(digits = 3)
model estimate std.error conf.low conf.high p.value
CUPED: ~ group + x_pre 3.749 2.209 -0.581 8.078 0.090
Multi-CUPED: + is_new 3.992 2.168 -0.258 8.241 0.066
Lin: group * (x_pre + is_new) 3.991 2.168 -0.259 8.242 0.066 
d = df.astype({"is_new": int})
# центрируем ковариаты для модели Лина
d["x_c"] = d["x_pre"] - d["x_pre"].mean()
d["new_c"] = d["is_new"] - d["is_new"].mean()

models = {
    "CUPED: ~ group + x_pre":        "revenue_w ~ group + x_pre",
    "Multi-CUPED: + is_new":         "revenue_w ~ group + x_pre + is_new",
    "Lin: group * (x_pre + is_new)": "revenue_w ~ group * (x_c + new_c)",
}
rows = []
for name, formula in models.items():
    fit = smf.ols(formula, data=d).fit(cov_type="HC2")
    ci = fit.conf_int().loc["group[T.test]"]
    rows.append({"model": name, "estimate": fit.params["group[T.test]"],
                 "std.error": fit.bse["group[T.test]"], "conf.low": ci[0],
                 "conf.high": ci[1], "p.value": fit.pvalues["group[T.test]"]})
pd.DataFrame(rows).round(3)
                           model  estimate  ...  conf.high  p.value
0         CUPED: ~ group + x_pre     3.749  ...      8.078    0.090
1          Multi-CUPED: + is_new     3.992  ...      8.241    0.066
2  Lin: group * (x_pre + is_new)     3.991  ...      8.241    0.066

[3 rows x 6 columns]

Оценка из первой модели совпадает с CUPED руками с точностью до третьего знака. Небольшая разница есть потому, что в CUPED θ считается по объединённым данным, а в регрессии наклон оценивается с учётом группы. Индикатор is_new добавляет ещё немного точности, как и обещала история.

Какую модель брать:

  • Сплит 50/50 – подойдёт любая из трёх, разница будет минимальной.
  • Неравный сплит (например, 10/90) и эффект может быть разным у разных пользователей – берите модель Лина. Модель без взаимодействия в таком случае может оказаться даже хуже обычного t-теста.
  • Стандартные ошибки всегда робастные (HC): lm_robust() в R, cov_type="HC2" в Python.

6. ARPPU: то же самое на платящих

До сих пор метрикой была выручка на пользователя (ARPU). С ARPPU (выручка на платящего) работает тот же пайплайн, только на подвыборке платящих: рандомизация по пользователям, у каждого платящего одно значение, так что обычный t-тест и CUPED применимы как есть. Отличий три.

1. Выигрыш от CUPED нужно пересчитать на истории, по платящим. Без нулей корреляция с предпериодом обычно ниже, чем у ARPU: часть «объяснения» у ARPU давал сам факт платил/не платил. Порог винсоризации тоже считаем по платящим, поэтому заново прогоняем prepare() на подвыборке.

2. Новые пользователи среди платящих. Их доля среди платящих может заметно отличаться от доли во всей аудитории (у нас она такая же, около 20%, но в реальных данных бывает по-разному). Индикатор is_new оставляем.

3. Подвыборка определяется после воздействия. «Платящий» считается по периоду теста. Если фича меняет конверсию в оплату, платящие в контроле и тесте – уже разные по составу люди, и сравнение ARPPU становится сравнением разных аудиторий. Это проблема и для t-теста, и для CUPED. Поэтому сначала проверяем конверсию.

# 1. На истории: сколько дисперсии объясняет предпериод у платящих
payers_hist <- history |> filter(revenue > 0) |> prepare()

tibble(
  metric    = c("ARPU: все пользователи", "ARPPU: платящие"),
  share_new = c(mean(history$is_new), mean(payers_hist$is_new)),
  variance_reduction = c(
    summary(lm(revenue_w ~ x_pre + is_new, data = history))$r.squared,
    summary(lm(revenue_w ~ x_pre + is_new, data = payers_hist))$r.squared
  )
) |>
  knitr::kable(digits = 2)
metric share_new variance_reduction
ARPU: все пользователи 0.2 0.36
ARPPU: платящие 0.2 0.24 
# 3. Тест: сначала конверсия в оплату
m_conv <- lm_robust(I(revenue > 0) ~ group, data = df)
tidy(m_conv) |> filter(term == "grouptest") |> select(term, estimate, p.value)
       term     estimate   p.value
1 grouptest -0.001604407 0.7466574
# ARPPU: t-test и Multi-CUPED на платящих
payers <- df |> filter(revenue > 0) |> prepare()
m_arppu       <- lm_robust(revenue_w ~ group, data = payers)
m_arppu_cuped <- lm_robust(revenue_w ~ group + x_pre + is_new, data = payers)

bind_rows(tidy(m_arppu)       |> mutate(model = "ARPPU: t-test"),
          tidy(m_arppu_cuped) |> mutate(model = "ARPPU: Multi-CUPED")) |>
  filter(term == "grouptest") |>
  select(model, estimate, std.error, conf.low, conf.high, p.value) |>
  knitr::kable(digits = 3)
model estimate std.error conf.low conf.high p.value
ARPPU: t-test 9.848 4.880 0.283 19.414 0.044
ARPPU: Multi-CUPED 10.768 4.184 2.567 18.968 0.010 
def r2(d):
    return smf.ols("revenue_w ~ x_pre + is_new", data=d.astype({"is_new": int})).fit().rsquared

# 1. На истории: сколько дисперсии объясняет предпериод у платящих
payers_hist = prepare(history[history["revenue"] > 0])
print(pd.DataFrame({"metric": ["ARPU: все пользователи", "ARPPU: платящие"],
                    "share_new": [history["is_new"].mean(), payers_hist["is_new"].mean()],
                    "variance_reduction": [r2(history), r2(payers_hist)]}).round(2))
                   metric  share_new  variance_reduction
0  ARPU: все пользователи        0.2                0.36
1         ARPPU: платящие        0.2                0.24
# 3. Тест: сначала конверсия в оплату
conv = smf.ols("paid ~ group", data=df.assign(paid=(df["revenue"] > 0).astype(int))).fit(cov_type="HC2")
print(f"\nКонверсия: разница {conv.params['group[T.test]']:.4f}, p = {conv.pvalues['group[T.test]']:.3f}\n")

Конверсия: разница -0.0016, p = 0.747
# ARPPU: t-test и Multi-CUPED на платящих
payers = prepare(df[df["revenue"] > 0]).astype({"is_new": int})
rows = []
for name, formula in {"ARPPU: t-test": "revenue_w ~ group",
                      "ARPPU: Multi-CUPED": "revenue_w ~ group + x_pre + is_new"}.items():
    fit = smf.ols(formula, data=payers).fit(cov_type="HC2")
    ci = fit.conf_int().loc["group[T.test]"]
    rows.append({"model": name, "estimate": fit.params["group[T.test]"], "std.error": fit.bse["group[T.test]"],
                 "conf.low": ci[0], "conf.high": ci[1], "p.value": fit.pvalues["group[T.test]"]})
pd.DataFrame(rows).round(3)
                model  estimate  std.error  conf.low  conf.high  p.value
0       ARPPU: t-test     9.848      4.880     0.283     19.413    0.044
1  ARPPU: Multi-CUPED    10.768      4.184     2.568     18.967    0.010

Конверсия от теста не изменилась (p = 0.75), так что ARPPU сравнивать можно. У платящих CUPED снизил дисперсию на 27%: меньше, чем у ARPU, как и предсказывала история. p-value снизился с 0.044 до 0.01.

lm_robust(revenue_w ~ group) без ковариат – это тот же t-тест Уэлча, только записанный как регрессия: так удобнее сравнивать с CUPED в одной таблице.

7. Как сломать CUPED: ковариата из периода теста

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

df_leak <- df |> mutate(x_leak = winsorize(rev_week1))
cat("Корреляция с метрикой:", round(cor(df_leak$x_leak, df_leak$revenue_w), 2), "\n")
Корреляция с метрикой: 0.97 
theta_leak <- cov(df_leak$x_leak, df_leak$revenue_w) / var(df_leak$x_leak)
df_leak <- df_leak |>
  mutate(revenue_leak = revenue_w - theta_leak * (x_leak - mean(x_leak)))
t_leak <- t.test(revenue_leak ~ group, data = df_leak)

bind_rows(summarise_t(t_cuped, "CUPED: выручка до теста"),
          summarise_t(t_leak,  "CUPED: выручка за 1-ю неделю теста")) |>
  knitr::kable(digits = 3)
method effect ci_low ci_high p_value
CUPED: выручка до теста 3.749 -0.581 8.078 0.09
CUPED: выручка за 1-ю неделю теста -1.086 -2.381 0.210 0.10 
df["x_leak"] = winsorize(df["rev_week1"])
print("Корреляция с метрикой:", round(df["x_leak"].corr(df["revenue_w"]), 2))
Корреляция с метрикой: 0.97
theta_leak = np.cov(df["x_leak"], df["revenue_w"])[0, 1] / df["x_leak"].var()
df["revenue_leak"] = df["revenue_w"] - theta_leak * (df["x_leak"] - df["x_leak"].mean())

pd.DataFrame({"CUPED: выручка до теста": welch("revenue_cuped"),
              "CUPED: выручка за 1-ю неделю теста": welch("revenue_leak")}).T.round(3)
                                    effect  ci_low  ci_high  p_value
CUPED: выручка до теста              3.749  -0.581    8.078     0.09
CUPED: выручка за 1-ю неделю теста  -1.086  -2.381    0.210     0.10

Корреляция почти идеальная, интервал стал втрое уже… а оценка эффекта стала -1.09 вместо 3.75: эффект не просто пропал, а сменил знак. Фича повлияла и на выручку первой недели, поэтому CUPED вычел эффект вместе с шумом (и даже с запасом). Самое неприятное, что выглядит это как очень точный тест. Авторы CUPED в оригинальной статье описывают ровно такой случай: ковариата из периода эксперимента дала оценку с неверным знаком.

Important

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

8. Проверяем на симуляции

Перед тем как внедрять CUPED в пайплайн, стоит прогнать его на A/A и A/B-симуляциях: ошибка первого рода должна остаться около 5%, а мощность должна вырасти. Генерируем 1000 датасетов без эффекта и 1000 с эффектом +3% и на каждом прогоняем тот же пайплайн.

analyze <- function(df) {
  df <- prepare(df)
  theta <- cov(df$x_pre, df$revenue_w) / var(df$x_pre)
  df <- mutate(df, revenue_cuped = revenue_w - theta * (x_pre - mean(x_pre)))
  tibble(
    method  = c("t-test", "CUPED"),
    p_value = c(t.test(revenue_w ~ group, data = df)$p.value,
                t.test(revenue_cuped ~ group, data = df)$p.value)
  )
}
n_sims <- 1000

sim <- expand_grid(effect = c(0, 0.03), sim = seq_len(n_sims)) |>
  mutate(res = map2(effect, sim, \(e, s) analyze(generate_data(effect = e, seed = s)))) |>
  unnest(res)

sim |>
  group_by(effect, method) |>
  summarise(share_significant = mean(p_value < 0.05), .groups = "drop") |>
  mutate(
    test = if_else(effect == 0, "A/A: ошибка I рода", "A/B +3%: мощность"),
    share_significant = scales::percent(share_significant, accuracy = 0.1)
  ) |>
  select(test, method, share_significant) |>
  pivot_wider(names_from = method, values_from = share_significant) |>
  knitr::kable()
test CUPED t-test
A/A: ошибка I рода 4.4% 5.2%
A/B +3%: мощность 55.6% 38.5% 
def analyze(df):
    df = prepare(df)
    y, x = df["revenue_w"], df["x_pre"]
    theta = np.cov(x, y)[0, 1] / x.var()
    y_cuped = y - theta * (x - x.mean())
    t = df["group"] == "test"
    return {
        "t-test": stats.ttest_ind(y[t], y[~t], equal_var=False).pvalue,
        "CUPED":  stats.ttest_ind(y_cuped[t], y_cuped[~t], equal_var=False).pvalue,
    }

rows = []
for effect in [0, 0.03]:
    for seed in range(1, 1001):
        p = analyze(generate_data(effect=effect, seed=seed))
        rows += [{"effect": effect, "method": m, "significant": p[m] < 0.05} for m in p]

(pd.DataFrame(rows)
   .groupby(["effect", "method"])["significant"].mean()
   .unstack())
method  CUPED  t-test
effect               
0.00    0.043   0.048
0.03    0.527   0.393

Здесь каждый прогон генерирует новые данные, а генераторы случайных чисел в R и Python разные, поэтому цифры во вкладках немного отличаются. Выводы одни и те же. В A/A оба метода держат ошибку первого рода около 5% (при 1000 прогонах это ±1.4 п.п.). А мощность у CUPED заметно выше: тот же эффект находится чаще на той же выборке. Если на ваших данных CUPED в A/A даёт заметно больше 5%, где-то ошибка: чаще всего ковариата посчитана не до теста.

Чеклист

Приложение: симуляция из поста

Код симуляции, цифры из которой приведены в посте в канале: логнормальная выручка без нулей, 5000 пользователей в группе, по 2000 экспериментов на сценарий, сравнение t-теста, разности «после минус до» и CUPED при корреляции около 0.3 и 0.7.

Показать код 
gen_revenue <- function(n, rho, effect = 0) {
  z_pre  <- rnorm(n)
  z_post <- rho * z_pre + sqrt(1 - rho^2) * rnorm(n)
  tibble(x_pre = exp(1 + 0.8 * z_pre),
         y     = exp(1 + 0.8 * z_post) * (1 + effect))
}

post_experiment <- function(n, rho, effect) {
  d <- bind_rows(gen_revenue(n, rho)         |> mutate(group = "control"),
                 gen_revenue(n, rho, effect) |> mutate(group = "test"))
  theta <- cov(d$x_pre, d$y) / var(d$x_pre)
  d <- mutate(d, y_change = y - x_pre, y_cuped = y - theta * (x_pre - mean(x_pre)))
  welch <- function(f) { r <- t.test(f, data = d); tibble(p_value = r$p.value, se = r$stderr) }
  bind_rows(welch(y ~ group)        |> mutate(method = "t-test"),
            welch(y_change ~ group) |> mutate(method = "Разность"),
            welch(y_cuped ~ group)  |> mutate(method = "CUPED")) |>
    mutate(corr = cor(d$x_pre, d$y))
}
Показать код 
set.seed(2026)
# rho задаётся на лог-шкале; наблюдаемая корреляция выручки получается ≈ 0.3 и ≈ 0.7
post_sim <- expand_grid(rho = c(0.38, 0.76), effect = c(0, 0.03), sim = 1:2000) |>
  mutate(res = map2(rho, effect, \(r, e) post_experiment(5000, r, e))) |>
  unnest(res)

post_sim |>
  group_by(rho) |>
  mutate(scenario = paste0("ρ ≈ ", sprintf("%.1f", mean(corr)))) |>
  group_by(scenario, effect, method) |>
  summarise(rate = mean(p_value < 0.05), se2 = mean(se^2), .groups = "drop") |>
  group_by(scenario, effect) |>
  mutate(var_ratio = se2 / se2[method == "t-test"]) |>
  ungroup() |>
  mutate(method = factor(method, levels = c("t-test", "Разность", "CUPED")),
         test   = if_else(effect == 0, "A/A", "A/B +3%")) |>
  arrange(scenario, test, method) |>
  transmute(scenario, test, method,
            rate      = scales::percent(rate, accuracy = 0.1),
            var_ratio = sprintf("×%.2f", var_ratio)) |>
  knitr::kable(col.names = c("Сценарий", "Тест", "Метод", "Доля p < 0.05", "Дисперсия vs t-test"))
Сценарий Тест Метод Доля p < 0.05 Дисперсия vs t-test
ρ ≈ 0.3 A/A t-test 4.8% ×1.00
ρ ≈ 0.3 A/A Разность 5.4% ×1.39
ρ ≈ 0.3 A/A CUPED 4.6% ×0.91
ρ ≈ 0.3 A/B +3% t-test 32.2% ×1.00
ρ ≈ 0.3 A/B +3% Разность 26.6% ×1.37
ρ ≈ 0.3 A/B +3% CUPED 35.4% ×0.91
ρ ≈ 0.7 A/A t-test 4.3% ×1.00
ρ ≈ 0.7 A/A Разность 5.2% ×0.60
ρ ≈ 0.7 A/A CUPED 5.3% ×0.51
ρ ≈ 0.7 A/B +3% t-test 34.5% ×1.00
ρ ≈ 0.7 A/B +3% Разность 52.8% ×0.59
ρ ≈ 0.7 A/B +3% CUPED 58.9% ×0.51
Показать код 
rng = np.random.default_rng(2026)

def gen_revenue(n, rho, effect=0.0):
    z_pre = rng.normal(size=n)
    z_post = rho * z_pre + np.sqrt(1 - rho**2) * rng.normal(size=n)
    return np.exp(1 + 0.8 * z_pre), np.exp(1 + 0.8 * z_post) * (1 + effect)

def post_experiment(n, rho, effect):
    xc, yc = gen_revenue(n, rho)
    xt, yt = gen_revenue(n, rho, effect)
    x, y = np.r_[xc, xt], np.r_[yc, yt]
    theta = np.cov(x, y)[0, 1] / x.var(ddof=1)
    out = {}
    for method, (a, b) in {"t-test":   (yc, yt),
                           "Разность": (yc - xc, yt - xt),
                           "CUPED":    (yc - theta * (xc - x.mean()), yt - theta * (xt - x.mean()))}.items():
        res = stats.ttest_ind(b, a, equal_var=False)
        out[method] = (res.pvalue, a.var(ddof=1) / n + b.var(ddof=1) / n)
    return out, np.corrcoef(x, y)[0, 1]

rows = []
for rho in [0.38, 0.76]:
    for effect in [0, 0.03]:
        for _ in range(2000):
            res, corr = post_experiment(5000, rho, effect)
            rows += [{"rho": rho, "effect": effect, "method": m, "sig": p < 0.05, "se2": v, "corr": corr}
                     for m, (p, v) in res.items()]

ps = pd.DataFrame(rows)
ps["scenario"] = "ρ ≈ " + ps.groupby("rho")["corr"].transform("mean").round(1).astype(str)
summary = ps.groupby(["scenario", "effect", "method"], sort=False).agg(rate=("sig", "mean"), se2=("se2", "mean"))
summary["var_ratio"] = summary["se2"] / summary.groupby(level=[0, 1])["se2"].transform("first")
summary[["rate", "var_ratio"]].round(2)
                          rate  var_ratio
scenario effect method                   
ρ ≈ 0.3  0.00   t-test    0.05       1.00
                Разность  0.05       1.39
                CUPED     0.05       0.91
         0.03   t-test    0.36       1.00
                Разность  0.26       1.37
                CUPED     0.39       0.91
ρ ≈ 0.7  0.00   t-test    0.05       1.00
                Разность  0.05       0.60
                CUPED     0.05       0.51
         0.03   t-test    0.36       1.00
                Разность  0.53       0.59
                CUPED     0.60       0.51
NoteЧто почитать дальше
  • Deng, Xu, Kohavi, Walker, 2013 – оригинальная статья про CUPED
  • Vickers & Altman, 2001, BMJ – почему в клинических исследованиях с измерением «до» и «после» ANCOVA лучше разностей
  • Lin, 2013 – регрессионная поправка в рандомизированных экспериментах
  • estimatr – R-пакет с lm_robust() и lm_lin()

Подписывайтесь на телеграм-канал!