ベイズファクターを用いた仮説検定~平均値の差の検定~

◆この記事のテーマ

ベイズファクターを使って平均値の差を検定する方法を整理します。 この手法は、例えば下図のようなシーンに適用できます。下図左側では、交通事故対策Aの実施前と実施後の単位時間あたり事故発生件数のデータを蓄積し、対策前後で事故発生件数の平均値に差があるかを検定しています。また、下図右側では交通事故対策Aが実施された交差点と実施されていない交差点の単位時間あたり事故発生件数の平均値に差があるかを検定しています。

※上記の事例は整数データであるからポアソン分布や負の二項分布等を使用するのが理想ですが、ここでは事故発生件数の平均値が十分大きく正規近似が可能であると仮定します。

※左側は対応のある(one-sample)グループ間の検定、右側は対応のない(two-sampe)グループ間の検定です。設定上、後者の設計では対策実施効果以外の要因(異なる交差点を比較に用いることによるランダム効果)が結果に絡む可能性が高くなるので、対策Aの効果を確かめることが分析の目的であれば前者の対応のある検定の設計がより望ましいといえます。

参考文献は以下の通り。ベイズファクターを使って平均値の差の検定が英語なら文献が多い。

effect size δの導入

ここでは、対応のある場合、対応のない場合に分けて検定の対象となるパラメータδを導入します。

対応のある場合

同じ対象に対し、2つの異なる状態で観測された値をそれぞれ確率変数XiYiからの観測値xiyiとおき、独立同分布の確率変数の組Xn=(X1,X2,,Xn)Yn=(Y1,Y2,,Yn)から観測値xn=(x1,x2,,xn)yn=(y1,y2,,yn)を得たとします。

まずは、簡単に確率変数XiYiの差が平均μの正規分布に従うと記載しちゃいます。

(1)XiYiNormal(μ,σ2),   i=1,,n

(1)式の唐突に現れた分散σ2(XiYi)のスケールを表現するパラメータです。これはそのときどきの検定対象によって異なる値をとります。このとき、平均μσの間に何の関係も与えない場合、μの値の大小がどれくらいの意味をもつものなのか把握が困難になってしまいます。そこで、平均μをeffect sizeと呼ばれる値δとスケールを指定する値に分け、スケールを指定する値をσとしちゃいます。

(2)μ=δσ

これでようやくスケールが単位あたりに調整されたeffect sizeδに着目すれば良いことになり、以降の議論を一般的な形で記述できるようになりました。


Xi及びYiδの関係のモデル化(対応のある場合) (3)XiYiNormal(δσ,σ2),   i=1,,n

対応のない場合

グループ1、グループ2に所属する対象からの観測値をそれぞれ確率変数XiYiからの観測値xiyiとおき、独立同分布の確率変数の組Xn=(X1,X2,,Xn)Ym=(Y1,Y2,,Ym)から観測値xn=(x1,x2,,xn)ym=(y1,y2,,ym)を得たとします。

まずは、XiYjがともに分散σの正規分布に従い、また2つの分布の平均がαだけ離れているとおきます。

(4){XiNormal(μα2,σ2),   i=1,,nYjNormal(μ+α2,σ2),   j=1,,m

この場合もやはり同様に、2つの分布の平均値の差αと分散σの間に何の関係も与えない場合、αの値の大小がどれくらいの意味をもつものなのか把握が困難になってしまいます。そこで、やはり同様に

(5)α=δσ

とおいてやれば、ようやくスケールが単位あたりに調整されたeffect sizeδに着目すれば検定が可能になり、こちらも以降の議論を一般的な形で記述できるようになります。


Xi及びYiδの関係のモデル化(対応のない場合) (6){XiNormal(μδσ2,σ2),   i=1,,nYjNormal(μ+δσ2,σ2),   j=1,,m

仮説の設定

(3)(6)式のδに着目し、従来通りに帰無仮説と対立仮説を設定します。


◆帰無仮説と対立仮説 従来の統計的仮説検定に対応する仮説は、 δの正負が明らかでない両側検定の場合、 Unrestricted  Model:{H0:δ=0H1:δ0 δが0以上の値をとることが明らかな片側検定の場合、 Orderrestricted  Model:{H0:δ=0H2:δ>0

以降では説明の簡略化のため、両側検定の場合の仮説を前提にします。

事前分布の設定

さて、ベイズの枠組みではすべての母数に事前分布を設定してやる必要があります。この事前分布の役割は、分析者が前もって知りうる情報を反映させるものである…と思っていませんか?
この考え方は主観ベイズsubjective Bayes)と呼ばれており、間違ってはいません。しかし今回の目的は検定です。それもベイズファクターを用いた検定であり、事前分布の設定が検定結果に大きく左右されることからも、なるべく事前分布に主観を与えたくありません。このような態度を、客観ベイズobjective Bayes)と呼び、ベイズファクターを用いた仮説検定の場合は客観ベイズの立場で事を考えることが重要になってきます(Rouder et al.(2009))。

古今東西、客観ベイズの分野(それもなぜか心理学研究分野)では、研究者の間で合意のとれた検定のための汎用性のある事前分布の設定方法が盛んに研究されてきました。平均値の差の検定のための汎用的な事前分布もいくつか開発されています。ここでは天下り的に、JZSの事前分布JZS proir)と単位情報事前分布unit information proir)について整理します。

◆JZSの事前分布

JZSの事前分布は検定対象のパラメータδまわりの事前分布を以下のように設定する。

(7){δNormal(0,g)gInverse Gamma(12,r22)π(σ2)1σ2

ここで、rはscale factorとよばれ、effect sizeのスケールを設定するために分析者が指定するパラメータである。またπ(x)はパラメータxの事前分布である。

これは、下記設定に等しい。

(8){δCauchy(0,r)π(σ2)1σ2

以上が対応のある平均値の差の検定における事前分布の設定である。対応のない平均値の差の検定ではさらにμに以下のとおり事前分布を設定する。

(9)π(μ)1


memo

コーシー分布は下記式で定義される。 Cauthy(y|μ,σ)=1πσ11+((yμ)/σ)2 μ=0σ=1としたコーシー分布を標準コーシー分布と呼び、下記式で定義される。 Cauthy(y|μ=1,σ=1)=1π11+y2 逆ガンマ分布は下記式で定義される。 Inverse Gamma(y|α,β)=βαΓ(α)yα1exp(βy)   (y>0)


◆単位情報事前分布

単位情報事前分布は検定対象のパラメータαまわりの事前分布を以下のように設定する。

(10){δNormal(0,r2)σα21σ2

以上が対応のある平均値の差の検定における事前分布の設定である。対応のない平均値の差の検定ではさらにμに以下のとおり事前分布を設定する。

(11)π(μ)1

JZSの事前分布と単位情報事前分布の唯一の違いはeffect size δの事前分布(正規分布)の分散の設定です。JZSの事前分布は分散gが逆カイ二乗分布に従うとした一方、単位情報事前分布は分散固定です。その意味では、単位情報事前分布の方がモデルに与える事前情報が多いことになります。またr=1のとき、2つの設定はδの事前分布を標準正規分布にするか、標準コーシー分布にするかの違いに帰着することが分かります。

標準正規分布と標準コーシー分布を比較した図が以下です。

図からも、標準コーシー分布は標準正規分布と比較し裾が重く、0から外れた値もとりやすい分布になっていることから、r=1の場合、JZSの事前分布は単位情報事前分布より事前情報が少ないことが確認できます。

また、両設定に共通のパラメータσ及び、対応のある検定の時のみ必要なμについては、ともにJeffeysの事前分布と呼ばれる事前分布を用いています。Jeffreysの事前分布はフィッシャー情報量から導かれるもので、事前分布と事後分布のカルバック・ライブラー情報量を最大化する効果がある、という点で無情報事前分布として使えるとされているようです。このあたりは詳しくないですが、Appendixに定義から正規分布の平均及び分散に対するJeffreysの事前分布を導出する方法を載せています。

ちなみに、σのように検定の対象となるパラメータとは別に存在し、かつ両モデルに共通の役割をもつパラメータはnuisance parameterと称され、こうしたパラメータはベイズファクターの値の決定にクリティカルでないことが証明されています。これに関しては以前の記事を参照してくだせえ。

ベイズファクターの解法

前節でモデルの設定を整理したので、本節ではベイズファクターの求め方を整理します。

解析的に求める

先に説明したJZSの事前分布と単位情報事前分布は解析的に計算することが可能です。ここでは導出や計算式までは深堀りしません。あまり面白くないので…。式を確認したい場合は参考文献を参照のこと。この文献の著者らはJZSの事前分布や単位情報事前分布の解析解を求めることができるWebページRパッケージを開発しているので、実践編で解析解を求める際はこのあたりを活用させてもらいます。

MCMC法で近似的に求める

本記事と実践編の本題はこちらです。これにはネストされたモデル同士の比較で有効なSavage-Dickey法を使ってベイズファクターを計算します。Savage-Dickey法の概要を以下に示します。詳細は以前の記事を参照のこと。


◆Savage-Dickey法

モデルM0M1に共通のパラメータδ=(δ1,,δm)があるとき、M0M1

(12)p(D|M0)=p(D|δ=δ0,M1)

とネストされた関係にある場合、ベイズファクターBF01は下記のとおり計算できる。

(13)BF01=p(δ=δ0|D,M1)p(δ=δ0|M1)

ここで、Dは全データセットである。

(11)(12)式において、H0H1にそれぞれ対応するモデルをM0M1と読み替えれば、H1に対応するモデルM1における対象パラメータの事前分布と事後分布の一点比較をするだけでベイズファクターBF01を求められることが分かります。

事前分布p(δ=δ0|M1)の方は、MCMCの結果からではなく、設定したモデルから値を算出します。たとえば、従来通りの仮説検定の枠組みでは、JZSの事前分布、単位情報事前分布の双方ともにδ=0が帰無仮説に対応するパラメータ空間上の点ですから、

(14){JZSp(δ=0|M1)=1πr    BF01=πr×p(δ=0|D,M1) p(δ=0|M1)=12πr2    BF01=2πr2×p(δ=0|D,M1) 

これが求めたいベイズファクターになります。

あとはp(δ=0|D,M1)が分かればベイズファクターを求められるので、以降はこれを求めることを考えます。

δの事後分布の求め方

条件の無い周辺事後分布を使う

こちらは全パラメータの同時分布からδ以外のパラメータψを周辺化消去することでδの事後分布を求めています。JZSの事前分布の場合、ψ=(σ,g,(μ))、単位情報事前分布の場合、ψ=σですね。rは分析者が事前に設定するものなのでこれに含まれません。

(15)p(δ=0|D,M1)=ψp(δ=0,ψ|D,M1)dψ=Eψ[p(δ=0,ψ|D,M1)]

(14)式はMCMCのサンプルを用いて近似的に求めるなら、ψを無視し、δのサンプルの分布をカーネル密度推定、もしくは以前の記事のように対数スプライン推定をするだけでイケます。簡単で便利です。

条件付き周辺事後分布を使う

こちらはψに関する条件付きのδの事後分布を使ってp(δ=0|xn,M1)を求める方法です。参考文献に倣ってConditional MArginal Density Estimator を略し、CMDEと呼びます。

(16)p(δ=0|D,M1)=ψp(δ=0|ψ,D,M1)dψ=Eψ[p(δ=0,|ψ,D,M1)]

(14)式の計算では切り捨てていたψについて、その事後分布もxnについての情報を含んでいるのだから、ψも考慮してp(δ=0|xn,M1)を求めてやるべきだ、というのがCMDEの発想ですが、実際、(16)式による計算は(15)式よりもはるかにp(δ=0|D,M1)の推定精度が良いことが参考文献による実験で解明されています。その理由はラオ・ブラックウェルの定理により明らかとのこと。勉強できたらAppendixに載せたいと思います。

(16)式は、MCMCサンプルを用いて次のように近似的に計算されます。

(17)Eψ[p(δ=0,|ψ,D,M1)]1Tt=1Tp(δ=0|ψ(t),D,M1)

(17)式は、条件付き事後分布p(δ|ψ,D,M1)の関数形が明らかであれば解析的に計算可能です。そこで、JZSの事前分布、単位情報事前分布それぞれの場合の条件付き事後分布p(δ|ψ,D,M1)及びその導出過程を、対応のある場合とない場合に分けて提示します。


◆対応のある平均値の差の検定におけるCMDE

JZSの事前分布を設定したときの条件付き事後分布p(δ|ψ=(σ,g),D=(xn,yn),M1)は以下のとおりである。

(18)p(δ|ψ,D,M1)=Normal(i=1n(xiyi)σ(n+1g),1n+1g)

単位情報事前分布を設定したときの条件付き事後分布p(δ|ψ=σ,D=(xn,yn),M1)は、

(19)p(δ|ψ,D,M1)=Normal(i=1n(xiyi)σ(n+1r2),1n+1r2)


◆導出 (18)式を導出する。 ベイズルールより、 p(δ|D,ψ,M1)p(D|δ,ψ,M1)p(δ|ψ) (3)式、(7)式を用いるため、 {p(D|δ,ψ,M1)=i=1np(xiyi|δ,σ)=i=1nNormal(xiyi|δσ,σ2)p(δ|ψ)=p(δ|g)=Normal(δ|0,g) よって、 p(δ|D,ψ,M1)=i=1nNormal(xiyi|δσ,σ2)×Normal(δ|0,g)=i=1n12πσexp(((xiyi)δσ)22σ2)×12πgexp(δ22g)exp(nσ2δ22i=1n(xiyi)δσ+i=1n(xiyi)22σ2δ22g)exp(12σ2nσ2(δ1nσi=1n(xiyi))2δ22g)exp(12((n+1g)δ221σi=1n(xiyi)δ))exp(121n+1g(δ1σi=1n(xiyi)n+1g)2)Normal(i=1n(xiyi)σ(n+1g),1n+1g) (19)式は(3)(10)式を用いるので、上述の導出でgr2に置き換えることで同様に導出できる。


◆対応のない平均値の差の検定におけるCMDE

JZSの事前分布を設定したときの条件付き事後分布p(δ|ψ=(σ,g,μ),D=(xn,ym),M1)は以下のとおりである。

(20)p(δ|ψ,D,M1)=Normal(2g(j=1myji=1nxi+(nm)μ)(4+(n+m)g)σ,4g4+(n+m)g)

単位情報事前分布を設定したときの条件付き事後分布p(δ|ψ=(σ,μ),D=(xn,ym),M1)は、

(21)p(δ|ψ,D,M1)=Normal(2r2(j=1myji=1nxi+(nm)μ)(4+(n+m)r2)σ,4r24+(n+m)r2)


◆導出 (20)式を導出する。 ベイズルールより、 p(δ|D,ψ,M1)p(D|δ,ψ,M1)p(δ|ψ) (6)式、(7)式を用いるため、 {p(D|δ,ψ,M1)=i=1nNormal(xi|μδσ2,σ2)×j=1mNormal(yi|μ+δσ2,σ2)p(δ|ψ)=p(δ|g)=Normal(δ|0,g) よって、 p(δ|D,ψ,M1)=i=1nNormal(xi|μδσ2,σ2)×j=1mNormal(yi|μ+δσ2,σ2)×Normal(δ|0,g)=i=1n12πσexp{(xi(μδσ2))22σ2}×j=1m12πσexp{(yj(μ+δσ2))22σ2}×12πgexp(δ22g)exp{12σ2(δσ(i=1nxij=1myj)(nm)μσδ+n+m4σ2δ2)12gδ2}exp{12((1g+n+m4)δ2+i=1nxij=1myj(nm)μσδ)}exp{124g4+(n+m)g(δ2g(j=1myji=1nxi+(nm)μ)(4+(n+m)g)σ)2}Normal(2g(j=1myji=1nxi+(nm)μ)(4+(n+m)g)σ,4g4+(n+m)g) (21)式は(6)(10)式を用いるので、上述の導出でgr2に置き換えることで同様に導出できる。

これでベイズファクターを使って平均値の差の検定をする準備が整いました。実践編は実践編で!

Appendix

Jaffrayの事前分布

◆Jeffreysの事前分布

Jeffreysの事前分布pJ(θ)は、下記のとおり定義される。

(22)PJ(θ)=I(θ)

ここで、I(θ)はフィッシャー情報量とよばれ、θは(通常確率変数の)パラメータである。

(23)I(θ)=E[(θln L(θ|x))2|θ]


この定義に沿って正規分布の平均・分散に対するJeffreysの事前分布を導出する。


正規分布の平均μに対するJeffreysの事前分布の導出 正規分布に従う確率変数Xについて考える。 L(μ|x)=12πσexp{(xμ)22σ2} より、 ln L(μ|x)=log(2πσ)(xμ)22σ2 だから、 μln L(μ|x)=(xμ)2σ2 フィッシャー情報量I(μ)は、 I(μ)=E[(xμσ2)2]=1σ2E[(xμσ)2]=αX,2σ2 ここで、αX,2は確率変数Xの2次の標準化積率である。標準化積率は αX,2=E[XμXσX2]=μX,2σX2 と表現されること(ここでμX,kXの中心積率、σXXの標準偏差)、μX,2Xの分散であることから、 I(μ)=αX,2σ2=1σ2 よって、 PJ(μ)=I(μ)=1σ1


正規分布の分散σ2に対するJeffreysの事前分布の導出 正規分布に従う確率変数Xについて考える。 L(σ2|x)=12πσexp{(xμ)22σ2} ここでσ2=gとおいて、 L(g|x)=12πgexp{(xμ)22g} このとき gln L(g|x)=(xμ)2g2g2 より、 フィッシャー情報量I(σ)は、 I(g)=E[((xμ)2g2g2)2]=14g2E[(xμ)4g2]21g2E[(xu)2g]+14g=14g2αX,421g2αX,2+14g=12g2 途中、αX,4=3αX,2=1を用いた。よって PJ(σ2)=PJ(g)=I(g)=12σ1σ

コメントを書く


※ コメントは承認されると表示されます

承認されたコメント一覧