Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

ベイズ統計学その20〜ベイズ推定を用いたt検定〜

□Bayesian Estimation Supersedes the t Test(BEST)

本日はベイズで平均値の差の検定を行っていきます。使用パッケージはBESTですが、JAGSもインストールしておく必要がありますので、ご注意ください。また、以下2本の論文と説明書を参考にしております。

John Kruschke(2013) Bayesian Estimation Supersedes the t Test, Journal of Experimental Psychology: Vol. 142, No. 2, 573– 603

Mike Meredith and John Kruschke(2015): Bayesian Estimation Supersedes the t-Test

Packages "BEST"

平均値の比較といえば、Student’s t testですが、正規性とか等分散性という制約がついていて、正規性がないと検定精度が下がってしまうので、正規性があるけど、等分散性がないならWelch’s t test、反対に正規性がないけど、等分散性があるならMann-Whitney U test、もはや両方ないならBrunner-Munzel testを行うなど、検定するまでに検討すべきポイントが多くて毎度毎度、平均値の比較には困らされます(個人的に・・・)。最近は、等分散も正規性も仮定しなくてよく、検定精度も高いBrunner-Munzel testにお世話になっております。

 

が、やはり検定なのでp値うんぬんで有意かどうかを判断しなければいけず、そもそも危険率5%(0.05)という判断基準が作為的なので・・・・みたいなことをグサッとケチをつけてくる人がたまにおりますので、そんな時こそマルコフ連鎖モンテカルロ法を利用してt testをベイズ推定します(きっと静かになってくれるはずです。まー最高事後密度区間も作意的ではありますが、、、)。また、事前の情報を反映させることができるので、不確かさをモデルに取り込むことができる点もベイズの魅力ですね。

Brunner-Munzel testの理論や検定精度については、以下のHPを参考にください。

http://d.hatena.ne.jp/hoxo_m/20150217/p1

https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.php

ベイズ推定におけるt testのモデル

基本的なモデルの説明に参ります。観測されたデータがt分布に従って分布すると考え、そのt分布のパラメタに事前分布を仮定する、階層性を持ったモデルです。詳しく見ていくと以下の図のようになります。

f:id:teruaki-sugiura:20160327105648p:plain

※John Kruschke(2013)を参考に作成

 

f:id:teruaki-sugiura:20160327113455p:plain

t分布を構成するパラメタ(μ、σ、ν)について、μには正規分布、σには一様分布、νにはずらし指数分布(はじめて聞きました・・・)を仮定します。そして、各t分布においてνだけ共通の事前分布となります。これが、Kruschke(2013)の論文で示されているモデルであり、事前分布に無情報事前分布を仮定している状態です。また、これはこの後に実行するBESTmcmcのpriors = NULLの状態です。ではRで動かしていきましょう。これから実践するのはKruschke(2015)で示されているモデルです。つまり、σとνにガンマ分布を事前分布として仮定しているモデルです。

f:id:teruaki-sugiura:20160327105824p:plain

※John Kruschke(2015)を参考に作成

> y1 <- c(5, 5, 4, 4, 3, 4, 7, 4, 6, 8)
> y2 <- c(2, 3, 5, 4, 3, 2, 3, 1, 4, 2)
> y <- c(y1, y2)
# priorsにlist形式で事前分布を渡すことができます。
# priors = NULLだとKruschke(2013)のモデルとなります
> priors <- list(muM = mean(y), muSD = sd(y)*5, sigmaMode = sd(y), sigmaSD = sd(y)*5, nuMean = 30, nuSD = 30)
> bayes.t.test.result <- BESTmcmc(y1, y2, priors=priors, thinSteps = 3, burnInSteps = 2000, parallel=FALSE)

f:id:teruaki-sugiura:20160327105921p:plain

> plot(bayes.t.test.result)
> meanDiff <- (bayes.t.test.result$mu1 - bayes.t.test.result$mu2)
> prob <- mean(meanDiff > 0)
> prob

f:id:teruaki-sugiura:20160327110036p:plain

平均値の差の平均が2.05で、95%HDIは[0.572~3.55]となっております。また、平均値の差が0以上の確率は99.5%となっています。SDはこんな感じです。

> plot(bayes.t.test.result, which="sd")

f:id:teruaki-sugiura:20160327110154p:plain

パラメタの事後分布も確認します。effSzは効果量のことです。

> summary(bayes.t.test.result)

f:id:teruaki-sugiura:20160327110225p:plain

事後分布が収束しているかも確認します。Rhatが1なので収束していると考えても良さそうです。

> print(bayes.t.test.result)

f:id:teruaki-sugiura:20160327110245p:plain

□ROPE(Region Of Practical Equivalence)

次はROPEについてみていきましょう。John Kruschke(2010)のDoing Bayesian Data Analysis: A Tutorial Introduction with RのP85には以下のように記述されています。

The 95% HDI is also one way for declaring which values of the parameter are deemed “credible”.

「95%Highest Density Intervalは、そのパラメタの値が信頼できるかを明言する方法の1つでもある。」

実質的に同じ値とみなす範囲=Region Of Practical Equivalence 

いまいち理解が進まないですね(訳が下手で申し訳ない・・・)。1つ例を考えてみましょう。例えば、θ=0.5である場合、ROPEの範囲を[0.45~0.55]と設定します。ベイズではパラメタを確率変数と考えますので、値が散らばるので、ぴったりθ=0.5にはなりません。しかし、95%HDIの中に、ROPEの範囲[0.45~0.55]が納まるのであれば、θを「実質的に同じ値とみなし、θ=0.5と確信する」という考えがこそが、ROPEの基本的な考えです。 

> plot(bayes.t.test.result, compVal=1, ROPE=c(1, 3))

f:id:teruaki-sugiura:20160327110424p:plain

> plotPostPred(bayes.t.test.result)

f:id:teruaki-sugiura:20160327110451p:plain

これはMCMCでの30個のランダムステップを選んで作られており、実際のデータもプロットされています。モデルがサンプルデータにうまいこと当てはまっているかどうかを視覚的に確認できます。

> plotAll(bayes.t.test.result)

すべてのパラメタの事後分布を確認できますが、ノートパソコンだと画面が小さすぎて(出力されるウインドウがでかい)エラーが出てしまいます。

f:id:teruaki-sugiura:20160327110523p:plain

以上でおしまい。

広告を非表示にする