Power分析：階層型モデル

「MCMC 使用」からの続き。

前回 power 算出に使ったスクリプト Jags-Ydich-Xnom1subj-MbernBeta-Power.R の「hypothetical parameter 仮説パラメータ」は定数だった

omega = 0.70
kappa = 2000

しかし現実的には、仮説パラメータは定数化できるほど明確ではない。そこで「ω, k を変数化」する、つまり「階層型モデル」で power 分析を行う。

左は「階層型モデル：Therapeutic Touch を検証」の階層図（図 9.7）。ω はグループ、θ は個人、k はグループをまたぐ個人の concentration 。

ただし、ここでは A_ω, B_ω, S_k, R_k は使用せず、以下のように、グループの平均と標準偏差を設定する

idealGroupMean = 0.65
idealGroupSD = 0.07

Therapeutic Touch を例にすると、「Therapeutic Touch 能力がある可能性」は、被験者全体をグループとした平均は 65% 、標準偏差は 0.07。つまり「58% から 82%」と仮定した。

この「仮説を支援する度合い」を決めるのが次の変数

idealNsubj = 10
idealNtrlPerSubj = 10

これは「仮説を支援するデータ数」となる。つまり、被験者の数（idealNsubj）と被験者毎の試行回数（idealNtrlPerSubj）。

次は図 13.2 で、前回の図 13.1 に "Real or Hypothetical World" の上段が加わったもの。

以下はスクリプト Jags-Ydich-XnomSsubj-MbinomBetaOmegaKappa-Power.R の主要部分で、図 13.2 との対応を色付けした。

idealGroupMean = 0.65
idealGroupSD = 0.07
idealNsubj = 100
idealNtrlPerSubj = 100
betaAB = betaABfromMeanSD( idealGroupMean , idealGroupSD )
theta = rbeta( idealNsubj , betaAB$a , betaAB$b )
theta = ((theta-mean(theta))/sd(theta))*idealGroupSD + idealGroupMean
theta[ theta >= 0.999 ] = 0.999 # must be between 0 and 1
theta[ theta <= 0.001 ] = 0.001 # must be between 0 and 1
z = round( theta*idealNtrlPerSubj ) 
dataMat=matrix(0,ncol=2,nrow=0,dimnames=list(NULL,c("y","s")))
for ( sIdx in 1:idealNsubj ) {
  yVec = c(rep(1,z[sIdx]),rep(0,idealNtrlPerSubj-z[sIdx]))
  dataMat = rbind( dataMat , cbind( yVec , rep(sIdx,idealNtrlPerSubj) ) )
}
idealDatFrm = data.frame(dataMat)
mcmcCoda = genMCMC( data=idealDatFrm , saveName=NULL , 
                    numSavedSteps=2000 , thinSteps=20 )
mcmcMat = as.matrix(mcmcCoda)
# Making the ω and k graphs
openGraph(width=7,height=3)
layout(matrix(1:2,ncol=2))
　　　　　　... 省略 ...                       
Nsubj = 2*7 ; NtrlPerSubj = 47  # 658 flips total
nSimulatedDataSets = min(500,NROW(mcmcMat))
simCount=0
for (simIdx in ceiling(seq(1,NROW(mcmcMat),length=nSimulatedDataSets))) {
  simCount=simCount+1
  genOmega = mcmcMat[simIdx,"omega"]
  genKappa = mcmcMat[simIdx,"kappa"]
  genTheta = rbeta(Nsubj, genOmega*(genKappa-2)+1,(1-genOmega)*(genKappa-2)+1)
  dataMat=matrix(0,ncol=2,nrow=0,dimnames=list(NULL,c("y","s")))
  for ( sIdx in 1:Nsubj ) {
    z = rbinom( 1 , size=NtrlPerSubj , prob=genTheta[sIdx] )
    yVec = c(rep(1,z),rep(0,NtrlPerSubj-z))
    dataMat = rbind(dataMat,cbind(yVec,rep(sIdx,NtrlPerSubj)))
  }
  goalAchieved = goalAchievedForSample( data.frame(dataMat) )
  goalTally = rbind( goalTally , goalAchieved )
}


グループと個人レベル、不確実な仮説

ここでは Jags-Ydich-XnomSsubj-MbinomBetaOmegaKappa-Power.R の実行結果を、本書に沿った 3 パターンで示す。尚、本書の subject-level を「個人レベル」と訳した。

Dec 24, 2015：本書掲載の実行結果と、ここに記載した私が得た実行結果とは若干異なる。「サンプリング誤差」と容認するには大きな誤差だが、本章のテーマに大きく影響しないので今のところは無視する。

実験#1：グループ優位

スクリプト中の変数 Nsubj, NtrlPerSub はそれぞれ、個人の数、個人毎試行数で、図 13.2 下段の "Simulated sample of data" に影響する。

idealGroupMean = 0.65
idealGroupSD = 0.07
idealNsubj = 100
idealNtrlPerSubj = 100

Nsubj = 2*7 ; NtrlPerSubj = 47

[1] "omegaAboveROPE: Est.Power=0.99; Low Bound=0.978; High Bound=0.996"
[1] "omegaNarrowHDI: Est.Power=0.998; Low Bound=0.991; High Bound=1"
[1] "thetasAboveROPE: Est.Power=0.998; Low Bound=0.991; High Bound=1"
[1] "thetasNarrowHDI: Est.Power=0; Low Bound=0; High Bound=0.006"

omega がグループレベル、thetas が 個人レベル。今回の結果は、グループレベルの power が、個人レベルに優っている。


実験#2：個人優位

個人数を半減、個人毎試行回数を 2 倍に変更。つまり合計試行数は変更せず、個人データ数を増やした。

Nsubj = 7 ; NtrlPerSubj = 2*47

[1] "omegaAboveROPE: Est.Power=0.812; Low Bound=0.776; High Bound=0.845"
[1] "omegaNarrowHDI: Est.Power=0.268; Low Bound=0.23; High Bound=0.308"
[1] "thetasAboveROPE: Est.Power=0.998; Low Bound=0.991; High Bound=1"
[1] "thetasNarrowHDI: Est.Power=0.998; Low Bound=0.991; High Bound=1"

power は omega で低下して、thetas で増加。

この結果も示すように「a general trend in hierarchical estimates（階層型見積りの一般的傾向）」は

個人レベルの見積りで高い精度を求めるなら、「個人のデータ数」を増やす必要がある。グループレベルで高い精度を求めるなら、「個人の数」を増やす必要がある（注：必ずしも「個人のデータ数」を多く要しないが、多いほど良い）。

つまり、個人を高い精度で見積る場合は、個人に関するデータ（個人の打撃数、など）が多く必要。グループの場合は、個人の数（メジャーリーグの選手の総数、など）が多く必要。

これは直感的にも納得する。


実験#3：不確実な仮説

「idealized hypothesis（理想的な仮説）」を不確実なものとしたスクリプトを実行する。「不確実さ」を、少ない「個人数」と、少ない「試行数」で定義して、スクリプトでは以下のように 100 から 10 に減少（脚色部分） 。

idealGroupMean = 0.65
idealGroupSD = 0.07
idealNsubj = 10
idealNtrlPerSubj = 10

これ以外の設定は実験#1と同じ。

上記グラフは、上段が実験#1の結果、下段が今回（本書の図 13.3 に相当）。明らかに下段で ω は不確実になり、k も縮小。上段と idealGroupMean, idealGroupSD は同じなので、個人数と試行数がこの結果の原因。


[1] "omegaAboveROPE: Est.Power=0.708; Low Bound=0.667; High Bound=0.747"
[1] "omegaNarrowHDI: Est.Power=0.698; Low Bound=0.657; High Bound=0.737"
[1] "thetasAboveROPE: Est.Power=0.854; Low Bound=0.821; High Bound=0.883"
[1] "thetasNarrowHDI: Est.Power=0; Low Bound=0; High Bound=0.006"

power もグループレベル と 個人レベルで小さくなった。

個人数と試行数が少ないと

the less certain hypothesis has reduced the power

つまり、「仮説の不確実性が power を減少させる」

This influence of the uncertainty of the hypothesis is a central feature of the Baysian approach to power analysis.

仮説に対する不確実性の影響こそが、Baysian による power 分析の大きな特徴。