StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

infer.netの例題シリーズ2 Tutorial 5: Clinical trial

今回は「Tutorial 5: Clinical trial」をやります。データは以下の通り。

outcometreat
00
00
10
00
00
11
01
11
11
11

2列目が薬の処理の 有/無、1列目はその結果で 1が治った/0が治らなかった となっています。問題は薬の処理は効果があるかどうか、すなわち治る確率を変えるかを判定することです。今回は以下のように判定します。

効果のある場合をmodel1として「治る確率は薬の処理によって変化し、処理有の時にp.treated、処理無の時にp.controlである」というモデル、効果の無い場合をmodel2として「治る確率は薬の処理によらず常にある確率(p.common)である」というモデルとしとします。このmodel1model2のどちらが正しいかを無情報事前分布を与えてモデル選択をさせるという判定の仕方です。

実装の唯一の難しい点はBUGSコードではInfer.netやSASMCMCプロシージャ)のようにIf文が書けない点です。BUGSではこういう時はstep()関数やequals()関数で頑張ります。 BUGSコードは以下になりました。

model {
   for (i in 1:N) {
      Y[i] ~ dbern(p[i])
      p[i] <- equals(model.index, 1)*q1[Treat[i]] + equals(model.index, 2)*q2
   }
   q1[1] <- p.control
   q1[2] <- p.treated
   q2    <- p.common
   p.control ~ dbeta(1,1)
   p.treated ~ dbeta(1,1)
   p.common  ~ dbeta(1,1)

   model.index ~ dcat(p.model[])
   p.model[1] <- 0.5
   p.model[2] <- 0.5
}

ポイントは4行目です。equals(e1,e2)という関数はe1 == e2の時のみ1、その他は0を取る関数です。model.indexが1か2 の値を取るのでこのようにしてIf文と同じことをさせています。また、確率の事前分布のところ、dbeta(1,1)の代わりにdunif(0,1)でよいのですが、Infer.netの例題にあわせました。

キックするRコードは省略します。結果は以下の通りです。

meansd2.5%25%50%75%97.5%Rhatn.eff
p.control0.34 0.22 0.04 0.17 0.29 0.45 0.91 1.00 1200
p.treated0.66 0.21 0.15 0.54 0.70 0.82 0.96 1.00 1200
p.common0.50 0.26 0.03 0.30 0.50 0.70 0.97 1.01 1200
model.index1.24 0.42 111121.00 1200

model1の方が優勢ですね。その時の治る確率は 処理有だと0.7ぐらい、処理無だと0.3ぐらいであることが分かりました。