Power分析：MCMC使用

「MCMC 未使用」からの続き。

MCMC で power を見積り、MCMC 未使用の結果と近似することを確認する。

以下は以前記した図 13.1 の下段と、図に対応する R スクリプトの実装。スクリプト中には、図の各フェーズに対応して色分けした。

omega = 0.70
kappa = 2000
sampleN = 74
nSimulatedDataSets = 1000
for ( simIdx in 1:nSimulatedDataSets ) {
  genTheta = rbeta(1,omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)
  sampleZ = rbinom(1,size=sampleN,prob=genTheta)
  simulatedData = c(rep(1,sampleZ),rep(0,sampleN-sampleZ))
  goalAchieved = goalAchievedForSample(simulatedData)
  goalTally = rbind(goalTally,goalAchieved)
}

nSimulatedDataSets がシミュレーション回数で、図中の「雲（"Representative...", "Simulated..."）」と「正方形（"Posterior..."）」の「多重度合い」に相当する。

関数 goalAchievedForSample 中で、図の "Prior for skeptical audience" はハードコードされている（厳密には関数 genMCMC の中、dbeta(1, 1) ）。

以下が、上記コードを含む Jags-Ydich-Xnom1subj-MbernBeta-Power.R の実行結果

[1] "ExcludeROPE: Est.Power=0.886; Low Bound=0.865; High Bound=0.905"

[1] "NarrowHDI: Est.Power=0.347; Low Bound=0.318; High Bound=0.377"

この ExcludeROPE の結果 Est.Power = 0.886 は、左の表 13.1 の ω = 0.70, サンプル数 = 74 の Power = 0.9  に対応。つまり 0.886 ≈ 0.9 より MCMC の結果と近似することを確認。

ここで、スクリプト中の変数の変化、power の集計方法を見る。

関数 goalAchievedForSample の一回の呼び出し結果は、例えば

> goalAchieved

$ExcludeROPE

[1] TRUE

$NarrowHDI

[1] TRUE

これは 図 13.1 の右下段 "Posterior distribution from simulated data sample: Are goals achieved?" の TRUE / FALSE 。これらを goalTally に集めて

> NROW(goalTally); colnames(goalTally)
[1] 1000
[1] "ExcludeROPE" "NarrowHDI"  

スクリプト実行結果の Est.Power = 0.886, 0.347 は、TRUE の数とサンプル数 1000 との割合

> sum(unlist(goalTally[,1])); sum(unlist(goalTally[,2])) 
[1] 886
[1] 347


次に kappa, sampleN を以下のように変更して

omega = 0.70
kappa = 10 #2000
sampleN = 91  #74

次のスクリプト実行結果を得た

[1] "ExcludeROPE: Est.Power=0.665; Low Bound=0.635; High Bound=0.694"
[1] "NarrowHDI: Est.Power=0.862; Low Bound=0.84; High Bound=0.882"

この NarrowHDI の結果 Est.Power = 0.862 は、左の表 13.2 の ω = 0.70, サンプル数 = 91 の Power = 0.8  に対応。つまり 0.862 ≈ 0.8 より MCMC の結果と近似することを確認。


「階層型モデル」に続く。