📑

ベイズファクターでの群間比較の評価

2024/01/11に公開

Python

Stan

tech

はじめに

先日、ベイズ流の多重比較法はないかと調べていたとき、岡田（2013）に出会いました。この論文では、ベイズ統計による情報仮説の評価（Bayesian evaluation of informative hypothesis）の枠組みを利用して、群間比較を評価する方法を紹介しています。ベイズ流A/Bテストなどの2群間比較の評価を紹介した日本語の記事はネット上に複数存在しますが、3群以上の群間比較をコード付きで紹介しているものは見当たりませんでしたので、自身の理解の整理も兼ねてまとめました。

この記事では、岡田（2013）とその参考文献であるHoijthink (2011) の内容を紹介し、2値データの3群比較を例として紹介します。また、岡田（2013）のOpenBUGSでのプログラムを参考にPythonとStanを使ったサンプルコードを作成しました。

参考文献

岡田謙介（2013）．「ベイズ統計による情報仮設の評価は分散分析にとって代わるのか?」『基礎心理学研究』32，223-231．

Hoijtink， H. (2011). Informative Hypotheses: Theory and Practice for Behavioral and Social Scientists, Chapman and Hall/CRC.
https://www.amazon.co.jp/dp/1439880514

例：2値データの3群比較

webページのデザインをdesign1、design2、design3の3つ用意して、ページを訪れたユーザーにランダムに3つのデザインのうち1つを表示してクリック率を測定します。その結果から各デザインを評価することを考えます。

design1が表示されたユーザー数を $n_1$ 、実際にクリックしたユーザー数を $x_1$ としたとき、design1でのクリック率を $\hat{\theta}_1=x_1/n_1$ とします。同様にdesign2のクリック率は $\hat{\theta}_2=x_2/n_2$ 、design3のクリック率は $\hat{\theta}_3=x_3/n_3$ と表します。このとき、 $\hat{\theta}_1 \sim \hat{\theta}_3$ に対応する確率変数^[1] $\theta_1 \sim \theta_3$ の大小関係に興味があるとします。具体的には、 $\hat{\theta}_1 < \hat{\theta}_2 < \hat{\theta}_3$ であったとき、次の4つの仮説のうちどれを支持すれば良いのか知りたいとします。

$H_1:$ design1よりもdesign2でのクリック率が高く、さらにそれよりもdesign3が高いという仮説

H_1:\theta_1 < \theta_2 < \theta_3.

H_2:

design2とdesign3のクリック率がdesign1よりも高いという仮説

H_2:\theta_1 < \{\theta_2, \theta_3\}.

H_3:

design1とdesign2のクリック率よりもdesign3が高いという仮説

H_3:\{\theta_1,\theta_2\} < \theta_3.

H_u:

design1、design2、design3のクリック率には不等式制約がない（大小の制約がない）という仮説

H_u:\theta_1,\theta_2, \theta_3.

上記の4つの仮説は、ベイズファクターやそれを変換した事後モデル確率を利用することで評価することができます。

なお、ベイズ統計による情報仮説の評価の枠組みでは、 $H_1 \sim H_3$ のようなパラメータに不等式制約（inequality constraint）をおく仮説を情報仮説（informative hypothesis）、 $H_u$ のようにパラメータに制約をおかないものは無制約仮説（unconstrained hypothesis）と呼びます。

ベイズファクターと事後モデル確率

ベイズファクター

ベイズ統計では仮説やモデルがどれだけデータが支持するのかを評価したい場合、ベイズファクター（bayes factor）を用いて評価します。データを $\bm{d}$ 、 $i(=1,\ldots,I)$ 番目の情報仮説を $H_i$ としたとき、仮説 $H_i$ と $H_u$ を比較するベイズファクターは

\textit{BF}_{iu} = \frac{p(\bm{d} | H_i)}{p(\bm{d} | H_u)}

と表せます。これはベイズの定理によって

\begin{aligned} \textit{BF}_{iu} &= \frac{p(\bm{d} | H_i)}{p(\bm{d} | H_u)}\\ &= \frac{p(H_i | \bm{d})p(\bm{d}) / p(H_i)}{p(H_u | \bm{d})p(\bm{d}) / p(H_u)}\\ &= \frac{p(H_i | \bm{d}) / p(H_u | \bm{d})}{ p(H_i)/p(H_u)} \end{aligned}

と書き直すことができます。分母の $p(H_i)/p(H_u)$ は2つの仮説の事前確率の比であるので事前オッズ（prior odds）、分子の $p(H_i | \bm{d}) / p(H_u | \bm{d})$ は事後確率の比なので事後オッズ（posterior odds）と呼ばれています。つまり、ベイズファクターは、仮説 $H_u$ に対して仮説 $H_i$ を支持する事前から事後のオッズの変化を表しています^[2]。

情報仮説の評価の枠組みでのベイズファクター

情報仮説の評価の枠組みでは、ベイズファクターを用いて情報仮説と無制約仮説を比べ、その良さを評価します。仮説 $H_u$ と $H_i$ を比較するベイズファクターは、以下の式で表すことができます。

\textit{BF}_{iu} = \frac{f_i}{c_i}

事前オッズ $c_i$ は無制約仮説 $H_u$ の事前分布のうち、情報仮説 $H_i$ と整合的な確率密度の割合を表します。これは仮説 $H_i$ の複雑さ（complexity）と呼ばれます。一方、事後オッズ $f_i$ は無制約仮説 $H_u$ の事後分布のうち、情報仮説 $H_i$ と整合的な確率密度の割合を表し、モデルへの当てはまり（fit）と呼ばれます。

複雑さ $c_i$ は単純な計算で算出できます。例えば、 $H_1$ は3つのパラメータに2つの不等式制約をおくものであるので、その場合の数は $3!=6$ 通りになります。 $H_1$ はそのうちの1通りなので確率は $1/6$ となります。無制約仮説の場合は、不等式仮説をおかないので1通りのうちの1通りで、確率は $1$ です。そのため、 $c_1 = p(H_1)/p(H_u) = (1/6)/ 1= 1/6$ となります。一方、 $H_2$ と $H_3$ は3つのパラメータに1つの不等式制約をおくもので、 ${}_3 \mathrm{C}_1=3$ なので3通り、そのうちの1通りなので確率は $1/3$ です。無制約仮説の事後確率は同様に1であるので、 $c_2$ と $c_3$ は $1/3$ となります。

モデルへの当てはまりは、事後分布からサンプリングしたデータのうち仮説に適合している割合を集計することで求めることができます。

事後モデル確率

情報仮説の評価の枠組みでのベイズファクターは事後モデル確率（posterior model probability, PMP）へ変換することが出来ます。 $i$ 番目の情報仮説のPMPは、仮説 $H_u$ に対する仮説 $I$ 個のベイズファクターと無制約仮説のベイズファクター $\textit{BF}_{uu}(=1)$ 　から

\textit{PMP}_i = \frac{\textit{BF}_{iu}}{1+\sum_{i}{\textit{BF}_{iu}}}

と表すことが出来ます。また、無制約仮説のPMPについては次のようになります。

\textit{PMP}_u = \frac{1}{1+\sum_{i}{\textit{BF}_{iu}}}

2値データの3群比較（トイデータへの適用）

Google Colaboratoryを使って、webページのデザインのトイデータを作成し、実際にベイズ統計による情報仮説の評価を行ってみます。

ライブラリのインストール

今回利用するライブラリをインストールします。

!pip install pystan==3.8.0
!pip install nest_asyncio==1.5.8
!pip install arviz==0.16.1

トイデータの作成

乱数を発生させてデータを作成します。

import numpy as np
import pandas as pd

np.random.seed(seed=123)

design_list = ["design1", "design2", "design3"]
n_list = [6000, 5970, 5950]
theta_list = [0.20, 0.22, 0.23]
x_list = [np.random.binomial(n, prob) for n, prob in zip(n_list, theta_list)]

data = pd.DataFrame({"design": design_list, "n": n_list, "x": x_list})
data["theta_hat"] = data["x"] / data["n"]
print(data)
#     design     n     x  theta_hat
# 0  design1  6000  1243   0.207167
# 1  design2  5970  1297   0.217253
# 2  design3  5950  1406   0.236303

事後分布の推定

クリック率は二項分布を仮定し、事前分布はベータ分布として事後分布をStanを使って推定します^[3]。このとき、ベータ分布のパラメータは $\alpha=\beta=1$ として、無情報事前分布を指定しています。Stanのコードは次のようになります。なお、A/Aテストや過去のテストで情報がある場合は、事前分布に反映しても良いです。

stan_code = """
data {
  int<lower=0> treat_n; //デザイン数
  array[treat_n] int<lower=0> x; //クリック数
  array[treat_n] int<lower=0> n; //サンプルサイズ
}
parameters {
  array[treat_n] real<lower=0, upper=1> theta; //クリック率
}
transformed parameters {
}
model {
for (i in 1:treat_n) {
  theta[i] ~ beta(1, 1);
  x[i] ~ binomial(n[i], theta[i]);
  }
}
generated quantities{
}
"""

作成したトイデータをデータフレーム形式からStanで扱える形式に変換します。

stan_data = {'treat_n': len(data['design']), 'x': data['x'].to_list(), 'n': data['n'].to_list()}
print(stan_data)
# {'treat_n': 3, 'x': [1243, 1297, 1406], 'n': [6000, 5970, 5950]}

最後にStanで事後分布を推定します。

import nest_asyncio
nest_asyncio.apply()
import stan
import arviz

posterior = stan.build(stan_code, data=stan_data, random_seed=111)
fit = posterior.sample(num_samples=2000, num_warmup=1000, num_chains=5, num_thin=1)
result = fit.to_frame()
# Building...

# Building: found in cache, done.Sampling:   0%
# Sampling:  20% (3000/15000)
# Sampling:  40% (6000/15000)
# Sampling:  60% (9000/15000)
# Sampling:  80% (12000/15000)
# Sampling: 100% (15000/15000)
# Sampling: 100% (15000/15000), done.
# Messages received during sampling:
#   Gradient evaluation took 1e-05 seconds
#   1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
#   Adjust your expectations accordingly!
#   Gradient evaluation took 9e-06 seconds
#   1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
#   Adjust your expectations accordingly!
#   Gradient evaluation took 6e-06 seconds
#   1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
#   Adjust your expectations accordingly!
#   Gradient evaluation took 7e-06 seconds
#   1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
#   Adjust your expectations accordingly!
#   Gradient evaluation took 7e-06 seconds
#   1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
#   Adjust your expectations accordingly!

パラメータの事後分布とトレースプロットは以下のように表示できます。マルコフ連鎖は定常状態に収束していそうです。

fig = arviz.plot_trace(fit)

図1. パラメータの事後分布（左）とトレースプロット（右）： $\theta_1$ は青、 $\theta_2$ はオレンジ、 $\theta_3$ は緑色で描画されている。また、線の種類はマルコフ連鎖の各chainを表している。

arviz.summary(fit, hdi_prob=0.95)

推定結果は次のようになっています。mean、sd、HDI(highest density interval)は問題なさそうです。Vehtari et al. (2021)^[4]の付録^[5]では、数値実験から収束判断の基準を示しています。この基準によれば、ESS（effective samples size）を表すess_bulkおよび、連鎖の裾が重い場合を考慮したess_tailはchain数×100（今回の場合は500）より大きく、r_hat（split- $\hat{R}$ ^[6]）は1.01より小さければ収束したと判断できる^[7]ので、収束したとしました。

index	mean	sd	hdi_2.5%	hdi_97.5%	ess_bulk	ess_tail	r_hat
theta[0]	0.207	0.005	0.197	0.217	10800.0	7557.0	1.0
theta[1]	0.217	0.005	0.207	0.228	11529.0	7745.0	1.0
theta[2]	0.236	0.005	0.226	0.247	10123.0	7343.0	1.0

ベイズファクターと事後モデル確率の計算

ベイズファクターと事後モデル確率は、事後分布のサンプリングの結果から次のように簡単に計算できます。

import numpy as np
import pandas as pd

# fit
# design1 < design2 < design3
f1 = np.mean(
    (result["theta.1"] < result["theta.2"]) & (result["theta.2"] < result["theta.3"])
)

# design1 < {design2, design3}
f2 = np.mean(
    (result["theta.1"] < result["theta.2"]) & (result["theta.1"] < result["theta.3"])
)

# {design1, design2} < design3
f3 = np.mean(
    (result["theta.1"] < result["theta.3"]) & (result["theta.2"] < result["theta.3"])
)

# complexity
c1 = 1 / 6
c2 = c3 = 1 / 3

# bayes factor
bf1 = f1 / c1
bf2 = f2 / c2
bf3 = f3 / c3
bf = [bf1, bf2, bf3] + [1]

# bayes factor & posterior model probability
labels = ["H_1", "H_2", "H_3", "H_u"]
bf_pmp = pd.DataFrame({"hypothesis": labels, "bf_iu": bf, "pmp": bf / sum(bf)})
bf_pmp

bf_pmpの中身

index	hypothesis	bf_iu	pmp
0	H_1	5.4276	0.4470324674255028
1	H_2	2.7309	0.22492463801538537
2	H_3	2.9829	0.24568006984367532
3	H_u	1.0	0.08236282471543643

ベイズファクターと事後モデル確率の解釈

ベイズファクター

無制約仮説 $H_u$ と比べて、各情報仮説の事後オッズは事前オッズから $H_1$ は5.4倍、 $H_2$ は2.7倍、 $H_3$ は3.0倍に変化しました。

事後モデル確率

各仮説の事後モデル確率は、 $H_1$ は45%、 $H_2$ は22%、 $H_3$ は25%、 $H_u$ は8%となるので、この中から仮説を一つ選択するならば、 $H_1$ が選択されます。

おわりに

今回、岡田（2013）やHoijthink (2011) を読み、情報仮説の評価でのベイズファクターを利用する方法を実際に手を動かして計算することで、今まで曖昧であったベイズファクターの知識が自身の中で定着したように感じました。Hoijthink (2011) には今回のような単純なケースだけでなく、もう少し複雑なケースも紹介しているようなので、そちらも読んでみようかと思います。

脚注

ベイズ統計ではデータ（ $\hat{\theta}_1 \sim \hat{\theta}_3$ ）を定数、母数（パラメータ）を確率変数として扱います。 ↩︎
ベイズ統計における尤度比検定に相当します。 ↩︎
今回、尤度関数は二項分布、事前分布に共役事前分布のベータ分布を指定しているので、事後分布はベータ分布となるため、実はMCMCで計算する必然性はないです。 ↩︎
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P. C. (2021). Rank-Normalization, Folding, and Localization: An Improved for Assessing Convergence of MCMC (with Discussion), Bayesian Analysis, 16(2), 667–718, doi:10.1214/20-BA1221. ↩︎
https://avehtari.github.io/rhat_ess/rhat_ess.html ↩︎
Gelman-Rubin 統計量 ↩︎
GelmanのBayesian Data Analysis Third Edition (2013) の基準より厳しくなってます。 ↩︎

はじめに

参考文献

例：2値データの3群比較

ベイズファクターと事後モデル確率

ベイズファクター

情報仮説の評価の枠組みでのベイズファクター

事後モデル確率

2値データの3群比較（トイデータへの適用）

ライブラリのインストール

トイデータの作成

事後分布の推定

ベイズファクターと事後モデル確率の計算

ベイズファクターと事後モデル確率の解釈

ベイズファクター

事後モデル確率

おわりに

Discussion