はじめに

ここでは、以前おさらいしたベイジアンフレームワークのもとでのABテストをrstanパッケージを利用して再現する。RStanはStanをRから使うためのインターフェースで、Stanは確率的プログラミング言語である。密度関数を指定して、MCMCサンプリングによるベイジアン統計推論による結果が得ることが可能。詳しい話は良き書籍やWebサイトがあるので、そちらにお任せする。

サンプルデータの準備

ABテストで使用するサンプルデータはBayesABパッケージのまとめの時に使用したサンプルデータを利用する。日付は14日に限定し、aパターンの方がコンバージョンレートが高くなるようにしている。

library(tidyverse)
library(bayesAB)
library(rstan)
library(scales)
# 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_summary <- df %>% 
  group_by(flag) %>% 
  summarise(
    cnt = n(),
    cv = sum(cv),
    not_cv = cnt - cv,
    cvr = cv / cnt
  )

df_summary

## # A tibble: 2 × 5
##   flag    cnt    cv not_cv   cvr
##   <chr> <int> <int>  <int> <dbl>
## 1 a     15575  1737  13838 0.112
## 2 b     15182  1586  13596 0.104

サンプルデータの準備が整ったので、これからベイジアンフレームワークのABテストを行う。

ベイジアンフレームワークのABテスト

事前分布として前回同様Beta(2, 10)を利用する。これは、コンバージョンレートは10%前後であると考えており、この確信度合いを観察されたデータから計算できる尤度に反映する。可視化するとこんな感じのベータ分布である。

prior_alpha <- 2
prior_beta <- 10

tibble(
  x = seq(0, 1, 0.01),
  y = dbeta(x, prior_alpha, prior_beta)
) %>% 
  ggplot(., aes(x, y)) +
  geom_line() + 
  geom_vline(xintercept = 0.10, linetype = "dashed") + 
  scale_x_continuous(breaks = seq(0, 1, 0.05)) + 
  labs(title = "Beta Distribution with alpha = 2, beta = 10", x = "Conversion Rate", y = "Density") + 
  theme_bw()

rstanでMCMCを行って、パラメタを推定する前に、モデルを定義しておく必要がある。ここでのモデルを説明する。sA(sB)はパラメータrateA(rateB)を成功確率として持つ二項分布に従うと仮定しており、nA(nB)回の試行中でrateA(rateB)の確率で成功する回数がsA(sB)回となる確率を計算している。また、rateAとrateBに対して、Beta(2, 10)を事前分布として設定。これは、成功確率が低いことを反映した分布。

model <- "
data {
  // Number of imp
  int nA;
  int nB;
  // Number of cv
  int sA;
  int sB;
}

parameters {
  real<lower=0, upper=1> rateA;
  real<lower=0, upper=1> rateB;
}

model {
  // 前回の事前分布の設定
  rateA ~ beta(2, 10);
  rateB ~ beta(2, 10);
  sA ~ binomial(nA, rateA);
  sB ~ binomial(nB, rateB); 
}
generated quantities {
  real rate_diff;
  rate_diff = (rateA - rateB) / rateB;
}
"

あとは観測されたデータを渡し、

data <- list(
  nA = df_summary$cnt[1],
  nB = df_summary$cnt[2], 
  sA = df_summary$cv[1], 
  sB = df_summary$cv[2]
  )

data

## $nA
## [1] 15575
## 
## $nB
## [1] 15182
## 
## $sA
## [1] 1737
## 
## $sB
## [1] 1586

コンパイルして、事後分布を推定する。チェインは4回、試行は1万回とする。

stan_samples <- stan(
  model_code = model, 
  data = data,
  chains = 4, 
  iter = 10000,
  seed = 1989
  )

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include   -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '4e512128d358f3d12936fb6a0b4cb873' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 6e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.022791 seconds (Warm-up)
## Chain 1:                0.024821 seconds (Sampling)
## Chain 1:                0.047612 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '4e512128d358f3d12936fb6a0b4cb873' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.023059 seconds (Warm-up)
## Chain 2:                0.026099 seconds (Sampling)
## Chain 2:                0.049158 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '4e512128d358f3d12936fb6a0b4cb873' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.022928 seconds (Warm-up)
## Chain 3:                0.027642 seconds (Sampling)
## Chain 3:                0.05057 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '4e512128d358f3d12936fb6a0b4cb873' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.02178 seconds (Warm-up)
## Chain 4:                0.026929 seconds (Sampling)
## Chain 4:                0.048709 seconds (Total)
## Chain 4:

パラメータの事後分布の収束を確認しておく。サンプリングチェーンが安定して、トレースが描かれているので、収束していると思われる。また、自己相関なども確認できないので問題はなさそう。rate_diffのトレースプロットをみるとわかるが、パターンAとパターンBの成功率の相対比にグループAに偏っていることがわかる。

traceplot(stan_samples)

相対比のところは、トレースプロットでみなくても事後分布のプロットを見たほうが早い。

plot(stan_samples)

最後に、前回もまとめているAとBの相対的な優秀さ度合いを可視化しておく。まずは、サンプリング結果を取り出す。

df_posterior <- as.data.frame(stan_samples)
head(df_posterior)

##       rateA     rateB  rate_diff      lp__
## 1 0.1108897 0.1057191 0.04890879 -10540.51
## 2 0.1085480 0.1047004 0.03674785 -10541.09
## 3 0.1116404 0.1058727 0.05447805 -10540.51
## 4 0.1069296 0.1055698 0.01288031 -10542.18
## 5 0.1123825 0.1037077 0.08364578 -10540.46
## 6 0.1125160 0.1000795 0.12426630 -10542.07

rate_diffはstanのモデルのgenerated quantitiesブロックで定義しておいたので、これが0より大きい回数を計算して平均すれば割合がわかる。結果として97%でパターンAが優れていることがわかる。

ab_result <- mean(sign(df_posterior$rate_diff > 0))
ab_result

## [1] 0.9769

事後分布の信用区間も計算しておく。5%:0.01はAをBと比較した際に、Aが101%以下の効果を出す確率が5%で、95%:0.12は、AをBと比較した際に、Aが112%以上の効果を出す確率が5%。つまり、AをBと比較した際に、Aが90%の確率で101%から112%の効果を出すだろうと解釈できる。

# 90% Credible Interval
quantile(df_posterior$rate_diff, probs = c(0.05, 0.95))

##         5%        95% 
## 0.01131337 0.12629568

最後に可視化しておく。

df_posterior %>% 
  mutate(col = if_else(rate_diff > 0 , 'high', 'low')) %>% 
  ggplot(aes(rate_diff, fill = col)) + 
  geom_histogram(bins = 100, alpha = 1/2) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_text(
    data = data.frame(),
    aes(x = 0.03, y = 50, label = percent(ab_result, accuracy = 0.01)), inherit.aes = FALSE, size = 5
    ) +
  scale_x_continuous(labels = percent_format(accuracy = 0.01)) +
  scale_fill_brewer(palette = "Set1") + 
  labs(subtitle = "Histogramn of (A - B) / B Samples: Probability", x = "(A-B)/B", y = "Density") + 
  theme_bw()

おまけ

(x - y) / yの数式のお気持ちをおまけとしてまとめておく。突然ベイジアンABテストの文脈で出てくるので、どういうことなんだろうかと疑問に思ったが、よく見ると、この式自体は売上高成長率(Sales Growth Rate)と同じである。xとyの差をyで割ることによって、xとyの相対的な大きさや関係性を表現しようとしているだけではあるが、おさらいを兼ねてメモ書きを残しておく。

式の分子(x - y)は、xとyの差を表している。これにより、xがyよりも大きい場合は正の値が、xがy よりも小さい場合は負の値が得られる。つまり、この差が意味しているのは、どれだけの大きさの差があるか、またその方向性は正なのか負なのかを意味する。

次に、分母yは除数で、差をどれだけの大きさで割るかを決める。yの値が大きいほど、差をより小さな比率で表すことになり、yの値が小さいと、差の影響がより強く反映される。

式全体で見ると(x - y) / yは、xとyの相対的な関係を示すことになる。例としてxが100でyが50の場合、絶対的な差は50。しかし、相対的な比率を計算すると、(100 - 50)/50 = 1となる。この場合、xはyの倍の値を持っていることが分かる。(x - y)/y = -0.5は、xがyの半分以下であることを示す。このように、相対的な比率を用いることで、数値間の関係性を相対的に評価できる。

ベイジアンABテストのはなし(rstan編)

はじめに

サンプルデータの準備

ベイジアンフレームワークのABテスト

おまけ

参考文献