Power分析：MCMC未使用

「Power とは Goal の達成確率」からの続き。

ここでは power とサンプル数を MCMC を使わない mathematical derivation で求める。

「データ分析はデータあってこそ」は、至極当然のことだが、こうやって数値的に捉えると更に納得する。そもそも偏ったサンプルは論外だが、自説を主張するには、最低限必要なサンプル数（データ数）というものがある。「そのサンプル数で大丈夫？」と疑うことは大切。

これ以降、本書の goal を「目的」と記す。

モード ω と power から最低サンプル数を算出する。ただし、サンプル数は「95%HDI は、0.50 近辺の ROPE を含まない」「95%HDI の幅は最大 0.2」の「目的」を達成するものである。

左は表 13.1 で、目的が「95%HDI 区間が ROPE[0.48, 0.52] を含まない」、例えば「コインに偏りがある」など。

表の読み方は、行 0.9 と 列 0.65 のサンプル数 150 の場合は

「95%HDI区間が ROPE[0.48, 0.52] を含まない」目的を達する確率 90% を望む場合、最低 150 のサンプル数が必要

つまり、サンプル数はモード ω（列）に反比例して、power（行）に比例する。これは直感的に理解できる。ω の増大は偏りの増大で少ないサンプル数で済み、サンプル数が多いほど power は高まる。

左は表 13.2 で、目的は「95%HDI の幅が最大 0.2」。最大幅を 0.1 に狭めると、0.2 より多くのサンプル数を要する。


以下、本書提供のスクリプト minNforHDIpower-example.R の実行結果で、表 13.1, 13.2 と同じ値を得る。

# Table 13.1
> show(sampSizeExclNull)
     Gen.Mode
Power 0.6 0.65 0.7 0.75 0.8 0.85
  0.7 238   83  40   25  16    7
  0.8 309  109  52   30  19   14
  0.9 430  150  74   43  27   16

# Table 13.2
> show(sampSizeMaxHdiWid)
     Gen.Mode
Power 0.6 0.65 0.7 0.75 0.8 0.85
  0.7  91   90  88   86  81   75
  0.8  92   92  91   90  87   82
  0.9  93   93  93   92  91   89

スクリプトでは、表の行と列のループで関数 minNforHDIpower を呼び出す。末尾に minNforHDIpower-example.R の主要部分を記す。

「MCMC 使用」に続く。


R サンプルコード

95%HDI区間が ROPE[0.48, 0.52] を含まない
# Specify desired powers:
desPow = seq(.7,.9,by=0.1)
# Specify generating modes:
genMod = seq(0.60,0.85,by=0.05)

# HDI excludes ROPE around nullVal:
# Declare matrix for storing results:
sampSizeMatrix = matrix( NA ,
        nrow=length(desPow) , ncol=length(genMod) ,
        dimnames=list("Power"=desPow,"Gen.Mode"=genMod) )

for ( desPowIdx in 1:length(desPow) ) {
  for ( genModIdx in 1:length(genMod)  ) {
    sampSize = minNforHDIpower(
       genPriorMode=genMod[genModIdx], genPriorN=2000,
       HDImaxwid=NULL, nullVal=0.5, ROPE=c(0.48,0.52),
       desiredPower=desPow[desPowIdx] , 
       audPriorMode=0.5 , audPriorN=2 ,
       HDImass=0.95 , initSampSize=5 , verbose=FALSE )
    sampSizeMatrix[desPowIdx,genModIdx ] = sampSize
    show(sampSizeMatrix)
  }
}

95%HDIの幅が最大 0.2
for ( desPowIdx in 1:length(desPow) ) {
  for ( genModIdx in 1:length(genMod)  ) {
    sampSize = minNforHDIpower(
       genPriorMode=genMod[genModIdx], genPriorN=10,
       HDImaxwid=0.20, nullVal=NULL, ROPE=NULL,
       desiredPower=desPow[desPowIdx] , 
       audPriorMode=0.5 , audPriorN=2 ,
       HDImass=0.95 , initSampSize=50 , verbose=FALSE )
    sampSizeMatrix[desPowIdx,genModIdx ] = sampSize
    show(sampSizeMatrix)
  }
}