ベイジアンABテストのはなし
UPDATE: 2023-08-26 00:51:32
はじめに
ここでは、いまさらではあるがベイジアンABテストの使い方を簡単にまとめておく。頻度主義とか、ベイズ主義とか、そのあたりの詳しい議論や数理的な側面は私のような一般人では立ち行かないトピックなので、統計学を専門とされている方にお譲りするとして、ここではベイジアンフレームワークのもとでのABテストを扱う方法をまとめておく。
一般的なABテストな問題点からはじめ、bayesABパッケージの基本的な使い方の紹介、ベイジアンフレームワークのもとでのABテストの利点と問題点をまとめておく。ベイズだから問題点がないわけではない。
5年前ほど前に今はなきブログに書いた内容をもとに修正、加筆しておく。bayesABパッケージの詳細は下記の通り。
一般的なABテスト
一般的な検定を用いてABテストを行う際、サンプルサイズに関して、「個人的」には1つの困難が毎回伴う。それは、サンプルサイズ設計において、効果量(例えばコンバージョンレートなどの)を決めないといけない点。
効果量については別のノートで扱っているので、ここでは詳しく扱わないが、例えば、\(\alpha=0.05\)、\(1-\beta=0.80\)として、パターンAが5％で、新しいパターンBが5％くらい高くなるだろうと見込んで10％とすると、およそ各グループで434人くらいをサンプリングすればよいということになる。
alpha <- 0.05
power <- 0.80
power.prop.test(n = NULL,
p1 = 0.05,
p2 = 0.10,
sig.level = alpha,
power = power,
alternative = 'two.sided'
)##
## Two-sample comparison of proportions power calculation
##
## n = 434.432
## p1 = 0.05
## p2 = 0.1
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* groupアカデミックな研究であれば、研究論文の過去の実験結果などから参考にできるかもしれないが、ビジネスだとそう簡単に参照にできる数字がなかったりする(ある場合も勿論ある)。それよりも何よりも、効果量という概念自体がビジネスサイドの万人に受け入れられない問題もある。広告であればサンプルサイズが足らない事例は少ないかもしれないが、設計したサンプルサイズが集まるまで、検定はできない。
bayesABパッケージ
ベイジアンフレームワークのもとで、ABテストを実行できるパッケージとしてbayesABがある。ベイズ統計の数理的な側面はここでは扱わないので、良き参考書に巡りあってください。今だとこちらの標準
ベイズ統計学でしょうか。
さておき、パッケージの使い方は簡単で、bayesTest()関数に、事前分布の設定とデータを渡せばよい。とくに事前情報がないので、ここでの事前分布の設定は無情報事前分布を利用。事前情報を利用することのの利点は下記の分析を読むとわかるので、下記を読んでください。女性の生理を題材に、妊娠しているかどうかを、事前分布を利用して、確率を推定するベイジアンモデルの話です。
library(tidyverse)
library(scales)
library(bayesAB)
set.seed(1)
# aのほうがレートが高い設定
a <- rbinom(100, 1, .20)
b <- rbinom(100, 1, .15)
ab <- bayesTest(a,
b,
priors = c('alpha' = 1, 'beta' = 1),
n_samples = 1e5,
distribution = 'bernoulli')print()を利用することで、ABテストの要約値を得られる。
print(ab)## --------------------------------------------
## Distribution used: bernoulli
## --------------------------------------------
## Using data with the following properties:
## A B
## Min. 0.00 0.00
## 1st Qu. 0.00 0.00
## Median 0.00 0.00
## Mean 0.17 0.16
## 3rd Qu. 0.00 0.00
## Max. 1.00 1.00
## --------------------------------------------
## Conjugate Prior Distribution: Beta
## Conjugate Prior Parameters:
## $alpha
## [1] 1
##
## $beta
## [1] 1
##
## --------------------------------------------
## Calculated posteriors for the following parameters:
## Probability
## --------------------------------------------
## Monte Carlo samples generated per posterior:
## [1] 1e+05summary()を利用することで、事後予想損失(Posterior
Expected
Loss)の情報が得られる。この結果をもとに、ABテストを早期に停止するかどうかを決定できる。P(A > B) by (0)%:の部分にある通り、aをbと比較した際に、aが57%の確率で良い効果をもたらすと判断できる。
summary(ab)## Quantiles of posteriors for A and B:
##
## $Probability
## $Probability$A
## 0% 25% 50% 75% 100%
## 0.0517386 0.1499328 0.1745308 0.2006559 0.3572681
##
## $Probability$B
## 0% 25% 50% 75% 100%
## 0.04751512 0.14065877 0.16452633 0.19031030 0.36070783
##
##
## --------------------------------------------
##
## P(A > B) by (0)%:
##
## $Probability
## [1] 0.57535
##
## --------------------------------------------
##
## Credible Interval on (A - B) / B for interval length(s) (0.9) :
##
## $Probability
## 5% 95%
## -0.3645588 0.7782955
##
## --------------------------------------------
##
## Posterior Expected Loss for choosing A over B:
##
## $Probability
## [1] 0.118263257%がどこから計算されたかというと、今回はモンテカルロサンプリングを100000回行っており、aとbの各試行の事後分布の値の比を計算して、bよりもaが優れている確率を計算していると思われる。
sim_len <- length(ab$posteriors$Probability$A)
sum(ab$posteriors$Probability$A / ab$posteriors$Probability$B > 1)/sim_len## [1] 0.57535この結果をわかりやすくしたのが下記の図。赤がaで、青がb。aがbよりも優れている場合、(a-b)/bはプラス、反対のaがbよりも劣っている場合、場合はマイナスになるため、このような可視化が可能になる。
plot(ab)[3]## $samples
## $samples$Probability
数字でイメージするとわかりよい。aがbよりも優れている場合はプラスになっている。
aa <- seq(0, 0.10, 0.01)
bb <- sort(aa, decreasing = TRUE)
res <- (aa-bb)/bb
flag <- sign(res)
data.frame(
aa, bb, res, flag
)## aa bb res flag
## 1 0.00 0.10 -1.0000000 -1
## 2 0.01 0.09 -0.8888889 -1
## 3 0.02 0.08 -0.7500000 -1
## 4 0.03 0.07 -0.5714286 -1
## 5 0.04 0.06 -0.3333333 -1
## 6 0.05 0.05 0.0000000 0
## 7 0.06 0.04 0.5000000 1
## 8 0.07 0.03 1.3333333 1
## 9 0.08 0.02 3.0000000 1
## 10 0.09 0.01 8.0000000 1
## 11 0.10 0.00 Inf 1sumamry関数から出力される信用区間のCredible Interval on (A - B) / B for interval length(s) (0.9):の部分が一番理解しにくい。
summary(ab)$interval## $Probability
## 5% 95%
## -0.3645588 0.7782955まず、5%: -0.36はaをbと比較した際に、aが64%以下の効果を出す(つまり悪化する)確率が5%で、95%: 0.78は、aをbと比較した際に、aが178%以上の効果を出す確率が5%。つまり、aをbと比較した際に、aが90%の確率で64%から178%の効果を出すだろうと解釈できる…ぱっと見だと混乱するが、100%を超えているのは相対値だから…という解釈で問題ないはず。
このようにベイジアンABテストであれば、P(A > B)とCredible Interval on (A-B)/Bの2つを利用して、サンプルサイズを設計せずとも、早期から分析を行って実験を停止するかどうかを判断できる。判断できるとはいえ、60%で効果ありなのか、70%で効果ありなのか、80%で効果ありなのかは分析者が決める必要がある。つまり、停止基準が必要になる。停止基準に関して検索してみると、いくつか文献がでてくるものの、効果量と同じく簡単に受け入れられるわけでもなさそうである…。
ベイジアンABテスト
ここでは、ベイジアンABテストの良さでもある停止基準を動かしながら意思決定を行う実践的(?)な使い方をまとめておく。サンプルデータとして、exploratory社のデータをお借りする。exploratory社の公開資料にベイジアンABテストの資料があり、データカタログに使用しているデータがあったので、それをお借りする。ただ、今回必要なデータは集計済みのデータではなく、集計前のログデータなので、集計値を利用して元の状態に擬似的に戻してから利用する。
集計データを集計値から擬似的に再現するスクリプトは下記のとおり。日付は14日に限定し、aパターンの方がコンバージョンレートが高くなるようにしている。
# exploratory社のデータカタログからお借りする
df_ab <- read_csv('https://exploratory.io/public/api/GMq1Qom5tS/A-B-IJp6BcB2/data')
# uniquePageViewとconversion_rateから集計前を再現するための関数
vec_gen <- function(x, y){
map2(
.x = x,
.y = y,
.f = function(x, y){rbinom(n = x, size = 1, prob = y)}
) %>% unlist()
}
df_a <- df_ab %>%
dplyr::filter(
landingPagePath == '/post?id=11' &
is_signup == TRUE &
date >= '2017-06-01' &
'2017-06-15' > date
)
df_b <- df_ab %>%
dplyr::filter(
landingPagePath == '/post?id=12' &
is_signup == TRUE &
date >= '2017-06-01' &
'2017-06-15' > date
)
dt <- seq(as.Date('2023-08-01'), as.Date('2023-08-14'), by = "day")
dt_a <- rep(dt, times = df_a$uniquePageView)
dt_b <- rep(dt, times = df_b$uniquePageView)
set.seed(1989)
cv_a <- vec_gen(x = df_a$uniquePageView, y = df_a$conversion_rate+0.015)
cv_b <- vec_gen(x = df_b$uniquePageView, y = df_b$conversion_rate)
df <- union_all(
tibble(dt = dt_a, cv = cv_a, flag = 'a'),
tibble(dt = dt_b, cv = cv_b, flag = 'b')
)
df %>%
group_by(flag) %>%
summarise(
avg_cv_rate = mean(cv)
)## # A tibble: 2 × 2
## flag avg_cv_rate
## <chr> <dbl>
## 1 a 0.112
## 2 b 0.104天下り的に、最終的なコンバージョンレートを時系列で可視化しておく。最終的な結果をABテスト前に確認することは本来は不可能ではある点は注意。ただ、一般的な統計検定を利用する場合は、仮に14日目まで待たないと必要なサンプルサイズが確保できないのであれば、14日目に判断を下すことになる。ベイジアンフレームワークであれば、3日目、5日目、7日目など自由な停止基準で分析して、意思決定に利用できる。ここでは、早期のテストの終了によって、不要な広告費を掛けずにすみ、CV数の増加も見込むことができる。
df %>%
group_by(dt, flag) %>%
summarise(
cnt = n(),
sum_cv = sum(cv),
rate = sum(cv)/n()
) %>%
ggplot(., aes(dt, rate, col = flag)) +
geom_line(size = 1) +
scale_x_date(labels = date_format("%Y/%m/%d"), breaks = date_breaks("1 day")) +
theme_bw() +
theme(axis.text.x = element_text(angle = 75, vjust = 0.2, hjust=0.2)) 
仮に3日目でABテストを分析してみたとする。結果をみると、この段階で91%の確率でAパターンのほうがよいと判断できる。停止基準がないことによる難しさは、実務ではわかりようがないが、ここでは本来はAパターンのほうが平均的に優れていないかもしれないが、この段階でAの方が良いと考えてしまい、意思決定を誤ってしまう可能性がある。
set.seed(1989)
ab3 <- bayesTest(df %>% filter('2023-08-03' >= dt & flag == 'a') %>% pull(cv),
df %>% filter('2023-08-03' >= dt & flag == 'b') %>% pull(cv),
priors = c('alpha' = 1, 'beta' = 1),
n_samples = 1e5,
distribution = 'bernoulli')
plot(ab3)[3]## $samples
## $samples$Probability
テストを継続して5日目で分析してみたとする。結果をみると、この段階で69%の確率でAパターンのほうがよいと判断できる。
set.seed(1989)
ab5 <- bayesTest(df %>% filter('2023-08-05' >= dt & flag == 'a') %>% pull(cv),
df %>% filter('2023-08-05' >= dt & flag == 'b') %>% pull(cv),
priors = c('alpha' = 1, 'beta' = 1),
n_samples = 1e5,
distribution = 'bernoulli')
plot(ab5)[3]## $samples
## $samples$Probability
テストを継続して7日目で分析してみたとする。結果をみると、この段階で95%の確率でAパターンのほうがよいと判断できる。そのため、ここでAパターンのほうが優れていると判断し、Aパターンに広告予算を割り振るという意思決定が可能である。
set.seed(1989)
ab7 <- bayesTest(df %>% filter('2023-08-07' >= dt & flag == 'a') %>% pull(cv),
df %>% filter('2023-08-07' >= dt & flag == 'b') %>% pull(cv),
priors = c('alpha' = 1, 'beta' = 1),
n_samples = 1e5,
distribution = 'bernoulli')
plot(ab7)[3]## $samples
## $samples$Probability
ベイジアンフレームワークのABテストは、サンプルサイズ設計が必要ではない一方で、停止基準を定める必要がある。また、停止基準を早めすぎてしまうと誤った意思決定にも繋がってしまう。ただ、あくまでも個人的な感覚ではあるが、ビジネスでは正解がないので、意思決定して施策を回し、改善を繰り返すことになるので、ベイジアンフレームワークのABテストの方が、ビジネスでは使いやすいと思う。また、今回はベイズの利点である事前分布をうまく使えていないので、事前分布を利用できればよりよい意思決定ができると思う。