今回は「Tutorial 5: Clinical trial」をやります。データは以下の通り。
outcome | treat |
---|---|
0 | 0 |
0 | 0 |
1 | 0 |
0 | 0 |
0 | 0 |
1 | 1 |
0 | 1 |
1 | 1 |
1 | 1 |
1 | 1 |
2列目が薬の処理の 有/無、1列目はその結果で 1が治った/0が治らなかった となっています。問題は薬の処理は効果があるかどうか、すなわち治る確率を変えるかを判定することです。今回は以下のように判定します。
効果のある場合をmodel1
として「治る確率は薬の処理によって変化し、処理有の時にp.treated
、処理無の時にp.control
である」というモデル、効果の無い場合をmodel2
として「治る確率は薬の処理によらず常にある確率(p.common
)である」というモデルとしとします。このmodel1
とmodel2
のどちらが正しいかを無情報事前分布を与えてモデル選択をさせるという判定の仕方です。
実装の唯一の難しい点はBUGSコードではInfer.netやSAS(MCMCプロシージャ)のように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コードは省略します。結果は以下の通りです。
mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | Rhat | n.eff | |
---|---|---|---|---|---|---|---|---|---|
p.control | 0.34 | 0.22 | 0.04 | 0.17 | 0.29 | 0.45 | 0.91 | 1.00 | 1200 |
p.treated | 0.66 | 0.21 | 0.15 | 0.54 | 0.70 | 0.82 | 0.96 | 1.00 | 1200 |
p.common | 0.50 | 0.26 | 0.03 | 0.30 | 0.50 | 0.70 | 0.97 | 1.01 | 1200 |
model.index | 1.24 | 0.42 | 1 | 1 | 1 | 1 | 2 | 1.00 | 1200 |
model1
の方が優勢ですね。その時の治る確率は 処理有だと0.7ぐらい、処理無だと0.3ぐらいであることが分かりました。