ベイジアンABテストのためにARPUのモデリングに挑戦してみた

はじめに

この記事は、DMMグループ Advent Calendar 2022 の24日目の記事になります。

はじめまして。マーケティング本部データ戦略部事業アナリシスグループ所属の石村遼汰(@Xc6Da)と申します。21卒エンジニア職としてDMMに就職し、データアナリストとして日々業務に取り組んでいます。

DMMには60を超える多種多様なサービスが展開されています。我々のチームでは、各アナリストに担当サービスが割り当てられ、担当サービスに対してデータを活用したグロース支援を行います。具体的には、データを分析して得られた示唆を元にした機能改善やキャンペーン等の施策提案、ABテストの実施、効果検証が主な業務となります。

この記事では、上で取り上げた業務の中からABテストについて話を掘り下げます。
最近よく耳にするベイジアンABテストについて、我々のチームで導入することが技術的に可能であるかを確かめました。我々がABテストの際にメインKPIとすることが多い ARPU(来訪ユーザーあたり売上) のモデリングについて、その試行錯誤の記録として本記事を書きました。

記事を書いたモチベーション

普段全く技術記事を書かない自分ですが、重い腰を上げて筆を執った理由は以下の通りです。

我々のチームでは連続値をとる変数をKPIにすることが多い
日本語で書かれたベイジアンABテストに関する記事では、CVRやクリック率など割合を推定した記事が多く、連続値をとる変数の実装例が少ない

特に2点目が大きなモチベーションになりました。会社のデータでベイズ推論を試してみようと「ベイジアンAB Python」とGoogle検索すると、1ページ目の検索結果9件中8件が二項分布を仮定したモデリングの実装が書かれており、連続値のモデリング記事がほとんどヒットしませんでした。ベイズ統計も統計モデリングも得意でない自分は諦めてしまいました。

しばらく経ち、改めて英語で「bayesian ab testing」と検索すると、連続値を対象としたモデリングの実装について書かれた英語記事がいくつかヒットしました。よくよく読んでみると意外と自分でもできそうだし、ちょうどアドベントカレンダーの時期だからということで、これらの英語記事を参考に試行錯誤し、記事を執筆することにしました。

数少ない日本語記事としてお役に立てれば嬉しいです。

なぜベイジアンABテスト?

一般的なABテストの効果検証では、頻度主義に基づいた統計的仮説検定(t検定、カイ二乗検定など)を用いることが多く、テストの結果得られた群間の差分が統計的に有意かそうでないかを調べます。これに対してベイジアンABテストでは、ベイズ推論を用いて群間の差分の事後分布(のパラメータ)を計算し、各群が他の群より優れている確率を直接的に算出することができます。

手法としての差分は上述の通りですが、それぞれについての詳細や比較はGoogle Optimize(GoogleのABテストツール)のヘルプページに分かりやすく記載されているので、そちらを参照ください。他にも最近は色々な会社でベイジアンABテストについてブログが書かれており、ソウゾウさん、AI Shfitさん、NIKKEIさんが各社それぞれ非常に分かりやすく説明されているので、そちらを読んでいただけると良いと思います。

ARPUの統計モデルを考える

ここからはモデリングや実装に入っていきます。
今回はPPLとして PyMC3 を利用しています。Pythonでベイズ推論を行うことができるライブラリは他にも NumPyro や TensorFlow Probability などがありますが、自分の手元で環境構築がうまくできたのが PyMC3 だったから、という理由で採用しています。

課題設定

単に指標の話だけでは想像がつきにくいと思いますので、今回はECサイトにおいて以下のようなABテストを実施したと仮定し、その効果検証を行う想定で進めていきます。

売上改善を目的としたUI改善の2群のABテスト
- テスト群に新しいUI、コントロール群に現行のUIが表示される
テスト期間は2週間
- サンプルサイズは数百万 ~ 数千万ほどの規模
メインKPIはARPU(来訪ユーザーあたり売上)

メインKPIがARPUなので、2群のARPUの差分を推定することが目的になります。今回はベイズ推論を用いて 2群のARPUの差分の事後分布 を推定します。

データの分布

テストの結果、ARPUとARPPU(購買ユーザーあたり売上)について、各群で以下のようなデータが得られたとします。

ゼロが過剰に多く、裾が重いという特徴がある分布になります。(こうした分布はECサイト等Webサービスのデータではよく見られる分布なのではないでしょうか。)

データの集計単位

今回はデータの集計単位を以下の3パターンに考え、それぞれモデリングを実施/比較しました。

未集計個別データ(以下、個別データ)
サンプルサイズ = 数百万 * 2
全期間のデータを集計したデータ(以下、全期間データ)
サンプルサイズ = 1 * 2
日ごとに集計したデータ(以下、日ごとデータ)
サンプルサイズ = 14

上述の通り、観測されるデータは数百万ほどの規模のため、全く集計せず全データを利用するのは集計データに比べて時間/金銭的コストがかかってしまいます。もし1 ~ 3で結果に差がないようだったら、可能な限り集計したデータを使いたい気持ちがあるため、この3パターンでモデリングを試しました。

ゼロが多いデータに対するモデリング

上で確認したような、ゼロが多いデータに対するモデリング方法のうち、今回参考になりそうなものは以下の2つが挙げられます。

ハードルモデル(Hurdle model)
データを 0 とそれ以外で別個の分布から生成されたと考えます。まず 0 か 0より大きい値かどうかの二値についてベルヌーイ分布から生成され、0が得られなかった場合は 0より大きい値をとる確率分布Dから生成されたとするモデルです。
つまり観測されたデータのうち、0はベルヌーイ分布から、0より大きい値は確率分布Dから生成されたデータとみなすことになります。
ゼロ過剰モデル(Zero-Inflated model)
0以上の値をとる確率分布Dから生成されたデータについて、ベルヌーイ分布によって一定の確率で0になるというモデルです。
つまり観測されたデータのうち、0はベルヌーイ分布と確率分布Dから、0より大きい値は確率分布Dから生成されたデータとみなすことになります。

今回モデリングの対象となるARPUについて、ARPU = CVR * ARPPU とKPI分解して考察することが多いこと、ARPUの生成メカニズムとしてハードルモデルの考え方の方が自然であることから、ハードルモデルを採用することにしました。

CVRのモデリング

まずはCVRの推定を行います。
は各群のCVRをパラメータとするベルヌーイ分布から得られると仮定します。また、CVRの事前分布として、ベータ分布(今回のケースは一様分布)を採用します。は割当を表す変数で、1がテスト群、0がコントロール群とします。

実装は以下の通りになります。なるべく要約統計量を使用したかったので、二項分布に従うとしてモデリングをしています。

import pymc3 as pm

with pm.Model() as model_cvr:
    # test-group
    cvr_t = pm.Beta('cvr_t', alpha=1.0, beta=1.0)
    cv_t = pm.Binomial('cv_t', n=data_t['visit_uu_cnt'], p=cvr_t, observed=data_t['purchase_uu_cnt'])
    
    # control-group
    cvr_c = pm.Beta('cvr_c', alpha=1.0, beta=1.0)
    cv_c = pm.Binomial('cv_c', n=data_c['visit_uu_cnt'], p=cvr_c, observed=data_c['purchase_uu_cnt'])

    # delta
    cvr_delta = pm.Deterministic('cvr_delta', cvr_t - cvr_c)

    # sampling
    trace_cvr = pm.sample(5000, chains=4, start=pm.find_MAP(), return_inferencedata=False)

結果は以下の通りです。r_hat < 1.1 であること、推定された事後分布に偏りがないことからCVRは問題なく推定されていると言えます。

ARPPUのモデリング

続いて、今回の肝となるARPPUのモデリングに入ります。
裾の重い分布のモデリングでは、指数分布、ガンマ分布、対数正規分布を利用することが多いようです。今回は指数分布、ガンマ分布を仮定したモデリング結果を紹介します。

対数正規分布についても試してみたのですが、推定されるARPUの差分が実際のデータから計算できるARPUの差分と乖離してしまったため、ARPPUの分布として使えないと判断しました。

他にも、指数分布とガンマ分布に対して最小値を制限した打ち切り分布も試してみましたが、手間の割に大きくモデリング結果が変わらなかったため、こちらも不採用としました。

指数分布

ARPPUが指数分布に従うとしてモデリングを行います。
なお、指数分布は 分散 = 期待値 **2 という性質があります。今回のデータでは裾がかなり重たいせいで 分散 >> 期待値 ** 2 という状態でしたので、指数分布を仮定するのは不適切な気がしましたが、とりあえず一発目ということで試してみました。

個別データ

個別データでは以下のようなモデルを仮定します。の事前分布は本当であれば共役なガンマ分布が良いと思いますが、手軽な無情報事前分布ということで広めの一様分布を採用しています。

実装は以下の通りです。個別データではレコード数が非常に多いせいか、デフォルトのサンプリングでは64GBのメモリがあってもjupyter notebookのカーネルがメモリエラーで死んでしまいます。これを回避するために return_inferencedata=False, idata_kwargs={'log_likelihood': False} を設定しています。

with pm.Model() as model_ex:
    # test-group
    lam_t = pm.Uniform('lam_t', 0, 10000)
    arppu_t_obs = pm.Exponential('arppu_t_obs', lam=lam_t, observed=data_t)
    arppu_t = pm.Deterministic('arppu_t', 1 / lam_t)
    
    # control-group
    lam_c = pm.Uniform('lam_c', 0, 10000)
    arppu_c_obs = pm.Exponential('arppu_c_obs', lam=lam_c, observed=data_c)
    arppu_c = pm.Deterministic('arppu_c', 1 / lam_c)

    # sampling
    trace_ex = pm.sample(5000, chains=4, start=pm.find_MAP(), return_inferencedata=False, idata_kwargs={'log_likelihood': False})

集計データ

全期間データ、日ごとデータの集計したデータでは以下のようなモデルを仮定します。
指数分布は再生性を持ち、指数分布の和はガンマ分布に従います。よって下記のように、売上の和、購買ユーザー数といった集計された値のみを用いてモデリングすることができます。

実装は以下の通りです。

with pm.Model() as model_ex2:
    # test-group
    lam_t = pm.Uniform('lam_t', 0, 10000)
    arppu_t_obs = pm.Gamma('arppu_t_obs', alpha=data_t['purchase_uu_cnt'], beta=lam_t, observed=data_t['sum_sales'])
    arppu_t = pm.Deterministic('arppu_t', 1 / lam_t)
    
    # control-group
    lam_c = pm.Uniform('lam_c', 0, 10000)
    arppu_c_obs = pm.Gamma('arppu_c_obs', alpha=data_c['purchase_uu_cnt'], beta=lam_c, observed=data_c['sum_sales'])
    arppu_c = pm.Deterministic('arppu_c', 1 / lam_c)

    # sampling
    trace_ex2 = pm.sample(5000, chains=4, start=pm.find_MAP(), return_inferencedata=False, idata_kwargs={'log_likelihood': False})

ガンマ分布

指数分布のセクションでも述べましたが、指数分布は 分散 = 期待値 ** 2 という性質があります。しかし今回のデータは明らかに過分散なデータのため、ARPPUが指数分布に従うとは言い難いです。
これに対してガンマ分布はパラメータが2つあるため、指数分布よりも柔軟なモデリングが可能と考えました。具体的には、今回のデータのような過分散なデータを再現できるのではという期待がありました。

個別データ

個別データでは以下のようなモデルを仮定します。α,βの事前分布は広めの一様分布を採用しています。

実装は以下の通りです。

with pm.Model() as model_gm:
    # test-group
    alpha_t = pm.Uniform('alpha_t', 0, 10000)
    beta_t = pm.Uniform('beta_t', 0, 10000)
    arppu_t_obs = pm.Gamma('arppu_t_obs', alpha=alpha_t, beta=beta_t, observed=data_t['sales'])
    arppu_t = pm.Deterministic('arppu_t', alpha_t / beta_t)
    
    # control-group
    alpha_c = pm.Uniform('alpha_c', 0, 10000)
    beta_c = pm.Uniform('beta_c', 0, 10000)
    arppu_c_obs = pm.Gamma('arppu_c_obs', alpha=alpha_c, beta=beta_c, observed=data_c['sales'])
    arppu_c = pm.Deterministic('arppu_c', alpha_c / beta_c)
    
    # sampling
    trace_gm = pm.sample(5000, chains=4, start=pm.find_MAP(), return_inferencedata=False, idata_kwargs={'log_likelihood': False})

集計データ

全期間データ、日ごとデータの集計したデータでは以下のようなモデルを仮定します。
ガンマ分布もまた再生性を持ち、ガンマ分布の和はガンマ分布に従います。よって下記のように、売上の和、購買ユーザー数といった集計された値のみを用いてモデリングすることができます。

実装は以下の通りです。

with pm.Model() as model_gm2:
    # test-group
    alpha_t = pm.Uniform('alpha_t', 0, 10000)
    beta_t = pm.Uniform('beta_t', 0, 10000)
    arppu_t_obs = pm.Gamma('arppu_t_obs', alpha=alpha_t * data_t['purchase_uu_cnt'], beta=beta_t, observed=data_t['sum_sales'])
    arppu_t = pm.Deterministic('arppu_t', alpha_t / beta_t)
    
    # control-group
    alpha_c = pm.Uniform('alpha_c', 0, 10000)
    beta_c = pm.Uniform('beta_c', 0, 10000)
    arppu_c_obs = pm.Gamma('arppu_c_obs', alpha=alpha_c * data_c['purchase_uu_cnt'], beta=beta_c, observed=data_c['sum_sales'])
    arppu_c = pm.Deterministic('arppu_c', alpha_c / beta_c)

    # sampling
    trace_gm2 = pm.sample(5000, chains=4, start=pm.find_MAP(), return_inferencedata=False, idata_kwargs={'log_likelihood': False})

結果

得られたサンプリング結果について、全てを表示すると煩雑になってしまうので、以下の観点で表形式にまとめました。

収束性
- r_hat が 1.1 より小さいかどうか
- 推定された事後分布に偏りが見られないか(目視で確認…)
実データとの適合
- 推定されたパラメータから生成したARPPUの分布と、実データのARPPUの分布が似ているか(目視で確認…)
計算速度
- サンプリング完了までの時間

収束性の確認については、 pm.summary や arviz.plot_trace 等の関数を用いました。
なお、個別データのサンプリング結果で上述の関数を使用するとメモリエラーでカーネルが落ちてしまうため、以下のような関数で結果を見られるようにしました。

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd


def summary(trace):
    res = []
    varnames = list(filter(lambda x: '__' not in x, trace.varnames))
    for varname in varnames:
        mean = trace[varname].mean()
        sd = trace[varname].std()
        hdi_low, hdi_high = az.hdi(trace[varname], hdi_prob=0.94)
        r_hat = az.rhat(np.array(trace.get_values(varname, combine=False)))
        res.append([varname, mean, sd, hdi_low, hdi_high, r_hat])
    res = pd.DataFrame(res, columns=['varname', 'mean', 'sd', 'hdi_3%', 'hdi_97%', 'r_hat'])
    res = res.set_index('varname')
    return res


def plot_trace(trace, min_height=4):
    varnames = list(filter(lambda x: '__' not in x, trace.varnames))
    figsize = [min_height * 3 * 2, len(varnames) * min_height]
    fig, axes = plt.subplots(len(varnames), 2, figsize=figsize)
    for varname, ax in zip(varnames, axes):
        chains = trace.get_values(varname, combine=False)
        for chain_samples in chains:
            # calculate and plot kde
            grid, pdf = az.kde(chain_samples)
            ax[0].plot(grid, pdf)
            ax[0].set_title(varname)
            
            # plot sampling result
            ax[1].plot([i for i in range(len(chain_samples))], chain_samples, alpha=0.3)
            ax[1].set_title(varname)
            ax[1].set_xlim(0, len(chain_samples) - 1)
    return axes

指数分布

指数分布の結果は以下のようになりました。どの集計方法でも結果に大差はないことから、指数分布を採用するのであれば計算量の観点から全期間データで集計してモデリングをするのが良いでしょう。
データとの適合については、以下の図が分かりやすいですが、指数分布の形でフィッティングさせるには難しい分布だなと感じさせられます。一例として、全期間データでの結果を掲載します。

	収束	データとの適合	計算速度
個別データ	◎	△	◎
全期間データ	◎	△	◎
日ごとデータ	◎	△	◎

ガンマ分布

ガンマ分布の結果は以下のようになりました。
個別データは時間こそかかったものの、収束性は問題ありませんでした。ただしデータとの適合について、こちらが期待していたような分散の大きい分布にはならず、また指数分布のモデリング結果と似通った結果になってしまいました。αが 1 に近い値に推定されていることからも、指数分布として推定されてしまったことが伺えました。

指数分布とは対象的に、ガンマ分布では集計データでのモデリングは上手くいきませんでした。日ごとデータでは収束こそしたものの、得られたパラメータの分布は分散が非常に大きく、またサンプリングされる値も実データと大きく乖離がありました。全期間データでは全てが☓ということで、望ましくないことが分かりました。

	収束	データとの適合	計算速度
個別データ	◎	△	△
全期間データ	×	×	×
日ごとデータ	◎	×	○

失敗例として、全期間データの結果を掲載します。ARPPUは収束しているように見えるものの、各種パラメータは全く収束していないことが伺えます。

ARPUの差分の事後分布

最後に、本稿の目的であるARPUの差分の事後分布を可視化しました。 arviz.kde でカーネル密度推定を行った結果を描画しました。なお収束しなかったガンマ分布-全期間データのモデリングは除外しています。

ガンマ分布-日ごとデータはかなり分散が広くなっているものの、その他については推定されたパラメータが似通っていたため、ARPUの差分の事後分布もほぼ同様の結果となりました。

よって今回試したモデルの中から1つ選ぶなら、計算量/時間の観点で優れている指数分布-全期間データを採用することになるでしょう。

まとめ

本稿では、ゼロが多く裾の重いARPUの統計モデリングを考えました。具体的には、ARPUをCVRとARPPUに分解し、それぞれ別の確率分布を仮定しモデリングを行うことでARPUの差分の事後分布を推定しました。

初めてベイズ推論をやってみて、いくつか疑問が生じました。

MCMCのパラメータはどれくらいに設定すべきか?
今回、MCMCはサンプリング数は5,000、バーンイン期間はデフォルトの1,000、チェーン数は4としてモデリングを行いました。これらの値に特に理由はなく、なんとなくで設定してしまいましたが、実データでの分析ではどのくらいのパラメータを設定することが多いのか勘所が分からないです。推論したいモデルのパラメータ数や使用するデータのサンプルサイズにも依ると思いますが、ガチャガチャ数値を弄ってみて勘所を身に着けていくものなのでしょうか…?
事前分布は無情報にすべきか?
今回パラメータの事前分布は、なるべく恣意性が入らないよう広めの一様分布を採用しました。しかし結果のセクションにもあるように、実データとの適合はイマイチかなと感じました。
時期にも依りますが、ARPUは値が極端にブレることはないため、過去のデータやAAテスト時のデータからある程度パラメータにあたりをつけることができます。もう少し値を絞りこむような事前分布を設定しても良いのかなと思いました。一方で、テストごとKPIごとに毎回事前分布を調整することになると、ABテスト実施に少なくない工数がかかってしまうことが予想されるため、あまり複雑にならない範囲で情報を盛り込んだ事前分布を設定できたら良いのかなと思いました。
実データとの適合はどれくらい気にすべきか?
上でも触れましたが、個人的には実データとの適合がイマイチな点がどうしても気になってしまいました。もちろん完全一致させることは難しいですし、目的であるARPUの差分については問題なく収束して"それっぽい"事後分布が推定されているので、ABテストの最終的な目標である「正しい意思決定」を考えると実データとの適合はそこまで追求しなくても良いのかもしれません。自分がベイズモデリングの思想が分からないので、有識者の方からコメントを頂けると嬉しいです。

また、上述のようなモデリングの課題をクリアしても、実際にベイジアンABテストを導入/運用する場合は更に以下のことをクリアする必要があると感じました。

勝敗の基準
- テスト群 > コントロール群となる確率は何%以上であれば勝ちとするのか?
AAテストの方法
- 2群の差がないとする基準は?
ベイズモデリングの知識がないメンバーへの展開
- 結果の解釈は容易な一方で、頻度主義の検定(t検定など)よりも結果を算出するのが大変
  - MCMCでのエラーハンドリングやパラメータの事前分布の設定にコツがいる
- 新たなKPIを設定するたびに、事前分布を探す必要がある
分析環境の構築
- 色々試したがライブラリの依存関係で躓くことが多い

試行錯誤した結果、ぶっちゃけ頻度主義のABテストの方が楽だなと感じてしまいましたが、3群以上のABテストが実行しやすいとか、サンプルサイズを気にしなくて済む等のメリットは大きいため、もう少し納得の行くまで調べてから導入可否を考えようと思いました。

さいごに

現在、事業アナリシスグループではデータアナリストを募集しています。
実データを用いた統計分析や機械学習に興味がある方、データを活用した事業貢献に携わりたい方に来ていただけると嬉しいです！ご興味ある方は下記リンクよりご応募ください。皆様の応募、心よりお待ちしております。

ベイジアンABテストのためにARPUのモデリングに挑戦してみた

はじめに

記事を書いたモチベーション

なぜベイジアンABテスト?

ARPUの統計モデルを考える

課題設定

ゼロが多いデータに対するモデリング

CVRのモデリング

ARPPUのモデリング

ARPUの差分の事後分布

まとめ

さいごに

参考文献

シェア

関連する記事

関連する求人