YEAST: Yet Another Sequential Test

Часть 1. Построение базового критерия

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

Определение метрики

В рамках эксперимента мы может на некоторый момент времени рассчитать разницу сумм интересующией нас метрики (прибыль, РТО и пр.) по A и B, у нас она будет обозначена как S_i, где i - индекс времени.

Но обычно данные от групп идут не одновременно, а друг за другом:

Значит расчет следует сделать инкрементальным: значения метрики из контроля мы будем складывать с предыдущим результатом суммы X, а из теста - вычитать.

Рассчитаем:

• в момент времени = 1 мы получаем данные от пользователя контрольной группы, например, 175;
• в момент времени = 2 уже от пользователя тестовой группы = 35,5

Отобразим в наших данных значение метрики от группы на некоторый момент времени как Y_i:

Следом рассчитаем результат инкремента, S_i:

• результат на момент времени = 2 следующий: 175 - 35,5 = 139,5

Чтобы не вычитать, а складывать, давайте для тестовой группы суммы умножать на (-1), тогда: 175 + (-35,5), что мы и видим в X_i. Сразу имеет смысл определить значение метрики от каждой группы с соответствующим знаком в зависимости от группы:

X_i - это метрика пользователя на некоторый момент времени

W_i - это переменная, которая показывает принадлежность к группе (0 - контроль; 1 - тестовая)

Y_i - значение метрики от этого пользователя с соответствующим знаком в зависимости от группы

Давайте рассмотрим пользователей двух групп на базе примера выше:

1. Из контрольной группы (W=0):

. Из тестовой (W=1):

Запишем результаты:

В тексте статьи таблица полностью выглядит так:

S_i=3 получается сложением 139.5+(-20), а S_i=4 сложением S_i=3 c 100: 119.5+100

Определим более формально S:

Две предварительные оценки

1) Определение количество инкрементов - N

Здесь все несложно: обычно, когда дизайним эксперимент, мы знаем кол-во дней, объем аудитории. Несложно оценить по историческим данным, сколько в среднем эта аудитория может дать значений метрики, строк, это и будет нашим N.

2) Дисперсия S при достижении N

Если знаем N, то примерно представляем себе, что при достижении этого N, наша метрика S могла бы как-то распределяться, тут сразу скажу, что нормально распределяться. Естественно, знать дисперсию мы не можем, но у нас есть из теории оценки а-ля стандартная ошибка. То есть хорошо бы иметь оценку ее дисперсии:

И... нет, Var(S_n) на момент окончания теста (=собрали N строк) у нас нет. Но у нас есть исторические данные! Проведем симуляции AA так, чтобы кол-во данных по метрике по строкам было = N, далее рассчитаем в каждой симуляции S, после - дисперсию распределения S, и именно она у нас будет Var(Sn). Да, это тоже оценка, которую мы и будем использовать сразу вместо V_n!

Дальше это V_n появится вновь, но оно нужно формально.

Критерий без последовательного тестирования (двусторонний)

b* - это у нас критическое значение на распределении S в абсолютных величинах; и да, для двустороненного теста мы будем смотреть абсолютное значение |S|

Если мы подставим V_n, то у нас будет просто дисперсия Var(Sn) на базе исторических данных. То есть вполне обычная граница значимости, разве что в абсолютной, нестандартизированной величине.

Часть 2. YEAST: Последовательное тестирование

Критерий (двусторонний)

У нас меняется граница, она становится меньше по альфе

Почему альфу в рамках двустороннего теста мы делим на 4, а не на 2? Здесь объяснение даст следующее неравенство:

Неравенство Леви

Смысл его в следующем: у нас есть кумулятивная сумма S от каждой новой строки данных, среди всех кумулятивов найдется максимальный. Например, [175, 139.5, 119.5, 219.5] -> максимальный = 219.5

Вероятность того, что максимальное значение будет больше критического порога, может быть в 2 раза больше вероятности, что конечный/последний кумулятив будет больше порогового, что определению = альфе! То есть:

Значит значение вероятности, что максимальный инкремент будет больше критического значения, не позволят нам контролировать ошибку 1-го рода согласно стандартной альфе; отсюда альфа должна быть снижена в 2 раза, чтобы вероятность относительно максимального вернулась к желаемой альфе;

2*a/2 = a

Поэтому для двустороненнего теста для левого и правого конца распределения вместо 0.025 остается 0.025 / 2 = 0.0125 (=0.05/4)

Похоже на одну и ту же границу Pocock, про которую писал в заметке про Group Sequential Testing

Алгоритм для последовательного тестирования:

1. Определяем ожидаемое кол-во N строк данных по целевой метрике, E(N)
2. Оцениваем дисперсию V распределения S через A/A при кол-ве строк = E(N)
3. Запускаем боевой тест и от первой до N строки проверяем, больше ли инкременет, чем заданный порог b*

В таком виде данное последовательное тестирование называется просто YEAST, можем делать столько проверок, сколько хотим.

Часть 2. Адаптация порогов значимости

(Вот тут будет уже немного сложно, но я. подстелю соломок по максимум!)

Наш b* рассчитан на конец сбора данных, i=N. Но что если мы хотим менять b в зависимости от того, что это за i?

Тогда нам нужно убедится, что False Discovery Rate, FDR, = alpha в рамках всех проверок с изменяемыми границами. FDR определяется так:

M_i - можно считать суммой S_i

M_k - в статье приняли за максимальную из S_i

Прочитаем вместе:

• P(M_1 > b_1) - это вероятность, что кумулятив будет больше порога b_1;
• P(M_1 <= b_1, M_2 > b_2) - вероятность, что следующий кумулятив (на момент 2) будет больше b_2, при условии, что в предыдущей проверке он был меньше или равен порогу b_1.

Зависимые вероятности

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

Пускай alpha = 0.05, а порог b_1 на базе такой альфы = 1.64 (односторонний критерий). Верна H0, поэтому мы будем прибавлять ко всем M случайную величину X из z-нормального распределения N(0, 1)

С первым членом все довольно просто, в рамках самого обычного критерия вероятность равна = alpha, то есть 0.05.

А вот со вторым все уже не так просто по двум причинам:

1.Это вероятность после НЕ наступления M_1 > b_1:

• P(M_1 <= b_1, M_2 > b_2) = P(M_2 > b_2 | M_1 <= b_1) * (1 - P(M_1 > b_1) = P(M_2 > b_2 | M1 <= b_1) * 0.95

2.Это условная вероятность: P(M_2 > b_2 | M1 <= b_1)

• Пускай порог b_2 = b_1 = 1.64. Проблема в том, что предыдущее значение, M1, это будет в наших условиях ненулевым значением. Например, M1=0.2. То есть мы уже стоит ближе к границе, а тест M2 еще даже не начался!
• Отсюда нам достаточно получить X > b_1 - M_1 > 1.64 - 0.2 > 1.46,
• Значит есть зависимость от предыдущего результата и:

P(M_2 > b_2 | M1 <= b_1) = 0.074

• Итого: P(M_2 > b_2 | M1 <= b_1) * (1 - P(M_1 > b_1) = 0.079*0.95 = 0.075

Замечу, при чем больше будет далее тестирований: M_3, M_4, ... тем меньше будет условная вероятность для каждого последующего M за счет предыдущих неудач, что логично при H0: каждая новая проверка при верности нулевой гипотезы тут не только риск ложноположительного результата, но и подтверждение самой этой гипотезы.

Двигаемся дальше: проблема - сумма всех этих вероятностей все-таки больше, чем > alpha:

P(M_1 > b1) + P(M_2 > b_2 | M1 <= b_1) * (1 - P(M_1 > b_1) = (0.05 + 0.075) > alpha

Именно поэтому мы должны двигать пороги b_1 до b_k, чтобы FDR был не более альфы.

Безусловная вероятность

Ок, пускай b_1 = 2, а b_2 = 2.5. Случайная величина X теперь принимает всего несколько значения (тут я сильно утрирую для лучшего понимания примера):

• P(M_1 = x_1 = -1.0) = 0.3
• P(M_1 = x_2 =0.5) = 0.4
• P(... x_3 = 1.8) = 0.2
• P(... x_4 = 2.2) = 0.1

Рассмотрим условные вероятности P(M_2 > b2 | M_1 = x_i):

Если M1​ = −1.0, то чтобы M2​ превысил 2.5, нужно набрать >2.5−(−1.0)=3.5. Чёт многовато, поэтому P(M_2​>2.5 | M_1​=−1.0) = 0.

Если M1​ = 0.5, то нужно набрать > 2.5−0.5 = 2.0. Это все еще много.

P(M2​>2.5 | M1​=0.5) = 0.05

Если M1​ = 1.8: >2.5−1.8=0.7. Вот это уже кажется более чем возможным:

P(M2​>2.5 | M1​=1.8) = 0.20

Если M1​ = 2.0: >2.5−2.0=0.5. Вроде изян:

P(M2​>2.5 | M1​=2.0) = 0.50

Далее рассчитаем безусловную вероятность P(M_2 > b2 | M_1 = x), когда могла наступить любая из реализаций x — это сумма всех (вероятность M1​ быть xi​) × (условная вероятность сработать на шаге 2, зная M1​=xi​).

1. Вклад от x1​ = −1.0: 0.30*0=0
2. от x2​ = 0.5: 0.40*0.05=0.020
3. x3​: 0.20*0.20=0.040
4. x4​: 0.10*0.50=0.050

P(M_2 > b2 | M_1 = x)=0+0.020+0.040+0.050=0.11, что можно выразить как:

Если у нас было бы много х был непрерывной величиной и значений было бесконечно много, тогда:

Замечу, что теперь мы умножаем на функцию плотности вероятности f_M1(x), что можно представить как степень "кучкования" значений на отрезке dx в момент M1.

Теперь можно увереннее вернутся к исходной записи:

Первый член нам уже известен:

где U_1 это var(S_N1) - просто перезапись.

Тогда:

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

Это не совсем та запись, которую ожидалось увидеть, поэтому раскрою:

• k - текущая проверка (из многих)
• k-1 - предыдущая проверка

Скажем, в k-1 было собрано N_k-1 наблюдений = 5, тогда полная сумма наблюдений:

N_1 + ... + N_k-1 = 3 + ... + 5 = 100

• X - конкретное наблюдение, оно же значение от пользователя, мы с этого начинали
• Пускай в k мы получили N_k = 2 наблюдения, то есть X_101, X_102

max{+0.5, +0.3} = +0.5, это максимальное значение, которого достигла кумулятивная сумма новых данных в момент k.

Более того, до k у нас была уже накоплена некоторая сумма (x), если

0.5+x > b_k, тогда стат. значимость

Или:

0.5 > b_k - x

Вновь используем неравенство Леви:

Отсюда мы возвращаем это в интеграл, который обозначим как I:

Вообще, то, что стоит в функции плотности f, это z-score:

Классическое отображение функции плотности для Z-значений, распределененных нормально выглядит так:

Тогда:

Возвращаем в интеграл:

Далее разложим интеграл на разницу между "1" и оставшимся:

Рассмотрим левую часть вновь через пример:

• P(x=0) = 0.1
• P(x=1) = 0.3
• P(x=4) = 0.6

Тогда:

Обобщаем:

В непрерывном случае было бы так:

Напомню, что такое Ф - это значение кумулятивной функции, то есть вероятности, от X (=b_k-1, только в стандартизированном виде), где каждая b отсортирована от большей к меньшей:

Проговорим: вероятность, что значение меньше -1 (можете представить тут b_k-1) ~ 0.18

Так как у нас:

Тогда:

Идем далее, обозначим значение Ф за некоторое Z (это не z-значение!):

Вернем в интеграл:

Вынесем Z_k за скобки, соблюдая все преобразования (см. красное):

Cделаем еще одну замену переменной:

Ну и почти готово!

Возвращаем в оценку FDR:

На само деле для всех промежуточных звеньев будет точно такое преобразование, где, выходит у нас сумма по K от 2 до k, что в итоге дает нам желанную оценку FDR!

То есть задав заранее кол-во подглядываний, оценку дисперсии кумулятива на каждое подглядывание и пр., мы можем вывести границы b, такие, что контролируют FDR, используя формулу для оценки!

Но:

1. Задать кол-во подглядываний K заранее, чтобы контролировать FDR в течение всего времени
2. Надо будет всего лишь (или - "всего лишь") рассчитать интеграл J, авторы предлагают применить методы, показанные здесь.

Алгоритм для последовательного тестирования с гибкими b:

1. Определяем кол-во подглядываний K
2. Определяем ожидаемое кол-во N строк данных на каждое K, от 1 до k
3. Оцениваем дисперсию V распределения S через A/A на каждое k
4. И некоторое маленькое значение эпислон = 0.0001, например
5. Считаем в первой итеграции b по нашим оценкам:

6. Рассчитаем FDR по формуле.

• Если больше заданной альфы, то к каждому b прибавляем малое эспилон*. Повторяем расчет FDR
• Если на уровне или меньше альфы, наши границы b подобраны!

*обратите внимание, делая b больше, мы по сути делаем жестче альфу на каждое подглядывание, что является аналогом Group Sequential Tests с alpha-spending подходом

В таком виде такое тестирование называется pYEAST.

Результаты симуляций

Ребята из Zalando решили сравнить свой метод с прочими, для каждого было по 100к симуляций, с отсутствием эффекта 0.0 и его наличием - 0.1; 0.2; 0.3. В колонке показана доля стат. значимых тестов.

• pYEAST 7  и 14 означают, что в одном случае K=7, в другом K=14

Обычный YEAST оказался самым мощным при том, что контролировал ошибку 1-го рода как положено.

А вот с точки зрения уменьшения кол-ва времени на тест, в основном выиграл pYEAST7:

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

Ниже представлены результаты "robust" и "non-robust", это значит, что для оценки дисперсии кумулятива использовались два разных способа рассчитать дисперсию: дело в том, что транзакционные данные могут зависимы от одного и того же пользователя, робастный это как раз учитывает.

Также рассмотрели случае мощность, когда кумулятив распределен ненормально, снова победитель простой YEAST

Ниже результат двух боевых тестов, в первом из которых они остановили тест раньше (на 4-ый день), что позволило им "сохранить деньги" (сохранить - не приумножить).

1. Единственное, что вызвало вопрос, так это то, что в начале статьи они определили тестовую группу как дающий X с множителем (-1). Получается, у них тестовая версия была лучше, а они говорят как будто наоборот.
2. Обратите внимание на ступенчатые границы: каждый сдвиг это момент подглядывания и, конечно, каждое следующее имеет более жесткий порог, тогда как в Group Sequential - наоборот.

А ведь у нас и по оценке FDR каждое следующая вероятность была ниже самой первой: дело не столько в порогах b, сколько в оценках дисперсии кумулятивной метрики на каждый момент.

Пускай revenue от каждой транзакции пользователей распределена N(0, sigma). Тогда сумма двух транзакций пользователей будет распределена согласно N(0, 2*sigma). То есть проблема сумм случайных величин в том, что при их сложении увеличивается дисперсия, значит, наши оценки U_k с увеличением k будут все больше и больше и именно устанавливает стартовые границы для b (далее корректируем через эпсилон)

Замечания по критерию:

1. Требует исторических данных

2. Если у вас сплит не 50/50, то значений по одной из групп в каждом моменте будет больше, что повлечет за собой ошибку 1-го рода больше ожидаемого. Авторы предлагают "down-sampling", - это в самом простом виде случайное удаление строк в большей выборке до размера меньшей. Делать это надо по-юзерно, а не по заказам: если удаляете user-111, то тогда удаляете и все его записи.

(авторы также предлагают перевзвешение результатов по сплиту (ratio) вместо down-sampling, но пока сами не придумали как)

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

4. От нашего воздействия может меняться и дисперсия поюзерной метрики X_i при том, что это не изменит среднего: это повлечет за собой больше ошибок 1-го рода

5. Воздействие может менять процесс порождения данных: например, пользователи будут реже платить, но при этом сумма оплаты станет больше, тогда конечная сумма-то не изменится, но кумулятив будет меняться куда более резко в большую/меньшую сторону с каждой строкой данных, что тоже будет давать больше ошибок 1-го рода.

Но при этом если ваша метрика счетная (count), то есть, например, просто считает, сколько у вас было конверсий, то п.3 и п.4 этой метрики не касаются.

6. Если у ваc A/B/C-тест, тогда вы можете смотреть результаты по отдельности: A/B, A/C; авторы также подчеркивают, что если есть существенные риски, что и B и C могут сделать результаты хуже, тогда стоит осуществить поглядывание А с B+C, предварительно проведя down-sampling.

Опорный материал

https://arxiv.org/pdf/2406.16523v1